Negative power function

What is negative power function

As a generalization of the quantum linear system problem (QLSP), we consider the following problem: given the access to an invertible matrix $A$ and a normalized quantum state $|b\rangle$, we want to construct a quantum state

$$|x\rangle:=\frac{A^{-c}|b\rangle}{\|A^{-c}|b\rangle\|_2},$$

where $c$ is a integer larger than 1.

The equivalent problem in function approximation

Similarly, the implement boils down to a scalar function $f(x)=x^{-c}$. Without loss of generality, we assume the matrix is normalized. We have to find a polynomial with parity approximating $f(x)$ over the interval $D_{\kappa}:=[-1,-\kappa]\cup[-\kappa,1].$

Setup parameters

For numerical demonstration, we set $\kappa = 10$, $c=2$ and scale down the target function by a factor of $1/(2\kappa^c)$ so that

$$f(x) = \frac{1}{2(\kappa x)^c}, \quad\max_{x \in D_\kappa} |f(x)| = \frac{1}{2}.$$

This improves the numerical stability.

kappa = 10;
targ = @(x) 1./x.^2;
% approximate f(x) by a polynomial of degree deg
deg = 150;
parity = 0;
opts.fscale = 1/(2*kappa^2);

Polynomial approximation using convex-optimization-based method

We want to find the best approximation polynomial with degree up to $d$ in terms of $L^{1}$ norm. To achieve it, we call cvx_poly_coef, which solves the optimization problem and outputs the Chebyshev coefficients of the best approximation polynomial.

opts.intervals=[1/kappa,1];
opts.objnorm = Inf;
opts.epsil = 0.2;
opts.npts = 400;
opts.isplot = true;

Then, simply calling the subroutine yields the coefficients of the approximation polynomial in the Chebyshev basis. As a remark, the solver outputs all coefficients while we have to post-select those of odd order due to the parity constraint.

coef_full=cvx_poly_coef(targ, deg, opts);
coef = coef_full(1+parity:2:end);

Polynomial approximation of the target function

Polynomial approximation error

Solving the phase factors for QLSP

We use Newton's method for solving phase factors. The parameters of the solver is initiated as follows.

% set the parameters of the solver
opts.maxiter = 100;
opts.criteria = 1e-14;
opts.useReal = false;
opts.targetPre = true;
opts.method = 'Newton';
% solve 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 $\ell^\infty$ norm

$$\mathrm{residual\_error} = \max_{k = 1, \cdots, K} |g(x_k,\Phi^*) - f_\mathrm{poly}(x_k)|.$$

Using $500$ equally spaced points, the residual error is $6.2172 \times 10^{-15}$ which attains almost machine precision. We also plot the point-wise error.

xlist1 = linspace(-1,-1/kappa,500)';
xlist2 = linspace(1/kappa,1,500)';
xlist = cat(1, xlist1,xlist2);
func = @(x) ChebyCoef2Func(x, coef, parity, true);
targ_value = targ(xlist);
func_value = func(xlist);
QSP_value = QSPGetEntry(xlist, phi_proc, out);
err= norm(QSP_value-func_value,2);
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')
print(gcf,'negative_power_function_error.png','-dpng','-r500');

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 = 3.55695e-05
max of solution = 0.801961
iter          err
   1  +1.4303e-01 
   2  +1.0329e-02 
   3  +8.1428e-05 
   4  +5.3832e-09 
Stop criteria satisfied.
The residual error is
   3.3761e-14