Quantum Hamiltonian simulation

Overview

Quantum Hamiltonian simulation is a class of problems targeting at implementing the map$H \mapsto \exp(-i \tau H)$for some Hermitian matrix $H$. Using QSP, the quantum hamiltonian simulation boils down to setting the target scalar function as $f(x) = \exp(-i\tau x)$. Due to the parity constraint, the implementation of the target scalar function is separated into implementing components $f_\mathrm{Re}(x) = \cos(\tau x)$ and $f_\mathrm{Im}(x) = \sin(\tau x)$ respectively. Then, these components are combined to realize $\exp(-i\tau H)$ using linear combination of unitaries (LCU). To illustrate the workflow, we set the evolution time to $\tau = 100$.

tau = 100
% set the parameters of the solver
opts.maxiter = 100;
opts.criteria = 1e-12;
% use the real representation to speed up the computation
opts.useReal = true;
% set the solver to fixed-point iteration
opts.method = 'CM';

Solving for the real component

To increase the numerical stability, the (real component of the) target function is scaled down by a factor $1/2$. Then, the target function is uniformly upper bounded by $1/2$. The Chebyshev coefficients of the target function is obtained by truncating the series to some finite degree. According to the analysis, it suffices to set $d = 1.4 |\tau| + \log\left(\frac{1}{\epsilon_0}\right)$so that the truncation error is lower than $\epsilon_0$.

targ = @(x) 0.5*cos(tau.*x);
parity = 0; % indicating the even parity
% compute the Chebyshev coefficients
d = ceil(1.4*tau+log(1e14));
f = chebfun(targ,d);
coef = chebcoeffs(f);
% discard coefficients of odd orders due to the even parity
coef = coef(parity+1:2:end);

Start the solver to find phase factors.

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

Verifying the solution

We verify the solved phase factors by computing the residual error in terms of the normalized $\ell^1$ norm

$$\mathrm{residual\_error} = \frac{1}{K} \sum_{k=1}^K |g(x_k,\Phi^*) - f(x_k)|.$$

Using $1,000$ equally spaced points, the residual error is $9.4988 \times 10^{-14}$ which attains almost machine precision. We also plot the pointwise error.

xlist = linspace(0, 1, 1000)';
targ_value = targ(xlist);
QSP_value = QSPGetEntry(xlist, phi_proc, out);
err= norm(QSP_value-targ_value,1)/length(xlist);
disp('The residual error is');
disp(err);
% plot the pointwise error
plot(xlist,QSP_value-targ_value)
xlabel('$$x$$', 'Interpreter', 'latex')
ylabel('$$g(x,\Phi^*)-f(x)$$', 'Interpreter', 'latex')

The point-wise error of the solved phase factors.

Reference

1. Low, G. H., & Chuang, I. L. (2017). Optimal Hamiltonian simulation by quantum signal processing. Physical review letters, 118(1), 010501.

2. 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).

3. 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

iter          err
   1  +1.7869e-01 
   2  +2.6363e-02 
   3  +3.4033e-03 
   4  +3.8478e-04 
   5  +4.0011e-05 
   6  +4.7448e-06 
   7  +6.9274e-07 
   8  +1.0545e-07 
   9  +1.5579e-08 
  10  +2.2009e-09 
iter          err
  11  +3.1451e-10 
  12  +4.3592e-11 
  13  +6.0816e-12 
  14  +8.6952e-13 
Stop criteria satisfied.
The residual error is
   9.4988e-14