Gibbs state

In this example, We introduce the preparation of Gibbs state.

What is Gibbs state

Gibbs state is a $n$-qubit mixed quantum state defined as a density matrix

$$\rho := \frac{1}{Z} e^{-\beta H},\ Z = \mathrm{tr}(e^{-\beta H}).$$

Gibbs state is important in the application of quantum algorithms in quantum simulation, quantum optimization, and quantum machine learning.

The equivalent problem in function approximation

The preparation of the Gibbs state can be viewed as implementing the mapping $H \mapsto \exp(- \beta H)$. Thus, the target function seemingly should be set to $f(x) = \exp(-\beta x)$. However, the parity condition is violated by this function and hence its implementation must be separated into the even and odd components respectively. The even or odd component of this function is hyperbolic cosine or sine whose value blows up exponentially. Hence, naively implementing this function would cause an exponentially vanishing unnormalized quantum state. To remedy this, the Hamiltonian can be first shifted and scaled so that it is strictly positive definite, i.e., there is a positive number $\delta > 0$ so that $\delta I \prec H \prec I$. Then, we only needs a target function that agrees with $f$ on the interval $D = [\delta, 1]$

$$\max_{x \in D} |h(x) - f(x)| < \epsilon.$$

Since the function is partially specified in a subinterval $D$, the best polynomial approximation $h(x)$ can be solved using convex-optimization based method.

Setup parameters

beta = 2;
targ = @(x) exp(-beta * x);
delta = 0.2;
deg = 151;

opts.intervals=[delta,1];
opts.objnorm = 2;
% opts.epsil is chosen to ensure that the target function is bounded by 
% 1-opts.epsil on D 
opts.epsil = 0.2;
opts.npts = 500;
opts.fscale = 1;
opts.isplot = true;
opts.maxiter = 100;
opts.criteria = 1e-12;
opts.useReal = false;
opts.targetPre = true;
opts.method = 'Newton';

Approximating the target function by polynomials

coef_full=cvx_poly_coef(targ, deg, opts);
parity = mod(deg, 2);
% only keep coefficients with consistent parity
coef = coef_full(1+parity:2:end);

Polynomial approximation of the target function

Polynomial approximation error

Solving phase factors by running the solver

[phi_proc,out] = QSP_solver(coef,parity,opts);

Verifying the solution

xlist = linspace(delta,1,500)';
func = @(x) ChebyCoef2Func(x, coef, parity, true);
func_value = func(xlist);
QSP_value = QSPGetEntry(xlist, phi_proc, out);
err= norm(QSP_value-func_value,Inf);
disp('The residual error is');

disp(err);
figure()
plot(xlist,QSP_value-func_value)
xlabel('$x$', 'Interpreter', 'latex')
ylabel('$g(x,\Phi^*)-f_\mathrm{poly}(x)$', 'Interpreter', 'latex')

The point-wise error of the solved phase factors.

Reference

1. Gilyén, A., Su, Y., Low, G. H., & Wiebe, N. (2019, June). Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics. In Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing (pp. 193-204).

2. Dong, Y., Meng, X., Whaley, K. B., & Lin, L. (2021). Efficient phase-factor evaluation in quantum signal processing. Physical Review A, 103(4), 042419.

Output of the code

norm error = 1.33551e-07
max of solution = 0.718311
iter          err
   1  +1.6068e-01 
   2  +8.3312e-03 
   3  +2.5267e-05 
   4  +2.2611e-10 
Stop criteria satisfied.
The residual error is
   1.0214e-14