# Quantum linear system problems

### Overview

The quantum linear system problem (QLSP) is a quantum variant of the linear system problem $$Ax=b$$. Given the access to an invertible matrix $$A$$ and to a normalized quantum state $$|b\rangle$$, QLSP asks a construction of the solution quantum state

$$
|x\rangle := A^{-1}|b\rangle / \Vert A^{-1} |b\rangle\Vert\_2
$$

with bounded error. 

Suppose the condition number of the coefficient matrix is $$\kappa := \Vert A\Vert\_2 / \Vert A^{-1}\Vert\_2$$ and $$A \in \mathbb{C}^{N \times N}$$. The recent progress reveals some near-optimal quantum algorithms for solving QLSP with bounded error $$\epsilon$$ and complexity $$\tilde{O}\left(\kappa\log\left(1/\epsilon\right)\mathrm{polylog}(N)\right)$$. In this example, we will deploy a much simpler (but sacrificing the total complexity) construction of the solution to QLSP.

### The equivalent problem in function approximation

The immediate idea for solving QLSP is approximately implementing the map $$A \mapsto A^{-1}$$. Hence, in light of QSP, the implementation boils down to a scalar function $$f(x) = x^{-1}$$. Without loss of generality, we assume the matrix is normalized so that $$\Vert A \Vert\_2 \le 1$$. Then, the condition number implies that the eigenvalues of $$A$$ lie in the interval $$D\_\kappa := \[-1, -\kappa] \cup \[\kappa, 1]$$*.* Hence, we have to find an odd polynomial approximation to $$f(x)$$ on the interval $$D\_\kappa$$. Here, the interval is set to $$D\_\kappa$$ rather than $$\[-1,1]$$ to exclude the singularity of $$f(x)$$ at $$x = 0$$.

<figure><img src="https://40428858-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FznDlSgNqsJLx79OHu47P%2Fuploads%2F8gIMxkNdAxUZAFPFRCGl%2Fquantum_linear_system_problems.jpg?alt=media&#x26;token=e0ecf1b4-bfc1-4e04-953b-e72eed52ba62" alt=""><figcaption><p>Approximating the matrix inversion within the interval of condition number.</p></figcaption></figure>

### Setup parameters

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

$$
f(x) = \frac{1}{2\kappa x}, \quad\max\_{x \in D\_\kappa} |f(x)| = \frac{1}{2}.
$$

This improves the numerical stability.

```matlab
kappa = 10;
targ = @(x) (1/(2*kappa))./x;
% approximate f(x) by a polynomial of degree deg
deg = 151; 
parity = mod(deg, 2);
```

### Polynomial approximation using convex-optimization-based method

To numerically find the best polynomial approximating $$f(x)$$ on the interval $$D\_\kappa$$, we use a subroutine which solves the problem using convex optimization. We first set the parameters of the subroutine.

```matlab
% set the parameters of the solver
opts.intervals=[1/kappa,1];
opts.objnorm = Inf;
opts.epsil = 0.2;
opts.npts = 500;
opts.fscale = 1; % disable further rescaling of f(x)ab
```

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.

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

To visualize this polynomial approximation, we plot it with the target function. From the figure, we can find both agree very well on $$D\_\kappa$$ but have distinct regularity otherwise.

<details>

<summary>visualize polynomial approximation</summary>

```matlab
% convert the Chebyshev coefficients to Chebyshev polynomial
func = @(x) ChebyCoef2Func(x, coef, parity, true);
figure()
tiledlayout(1,2)
nexttile
hold on
xlist1 = linspace(0.5/kappa,1,500)';
targ_value1 = targ(xlist1);
plot(xlist1,targ_value1,'b-')
plot(-xlist1,-targ_value1,'b-')
xlist1 = linspace(-1,1,1000)';
func_value1 = func(xlist1);
plot(xlist1,func_value1,'-')
hold off
xlabel('$x$', 'Interpreter', 'latex')
ylabel('$f(x)$', 'Interpreter', 'latex')
legend('target function', '', 'polynomial approximation')

nexttile
plot(xlist,func_value-targ_value)
xlabel('$x$', 'Interpreter', 'latex')
ylabel('$f_\mathrm{poly}(x)-f(x)$', 'Interpreter', 'latex')
print(gcf,'quantum_linear_system_problem_polynomial.png','-dpng','-r500');
```

</details>

<figure><img src="https://40428858-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FznDlSgNqsJLx79OHu47P%2Fuploads%2FKctce4Tio5qEo8MNFB2y%2Fquantum_linear_system_problem_polynomial.png?alt=media&#x26;token=fdcdc8e3-3902-4393-8e7f-7b94087d87c5" alt=""><figcaption><p>Left: The odd polynomial approximation to the target function solved using the convex-optimization-based method. Right: The point-wise approximation error.</p></figcaption></figure>

### Solving the phase factors for QLSP

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

```matlab
% 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;
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 $$1,000$$ equally spaced points, the residual error is $$6.2172 \times 10^{-15}$$ which attains almost machine precision. We also plot the point-wise error.

```matlab
xlist = linspace(0.1,1,1000)';
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,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')
print(gcf,'quantum_linear_system_problem_error.png','-dpng','-r500');
```

<figure><img src="https://40428858-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FznDlSgNqsJLx79OHu47P%2Fuploads%2FnN22VyMelt4nuP22X6qF%2Fquantum_linear_system_problem_error.png?alt=media&#x26;token=2c2b883f-153a-4dc0-afdc-1db8908a25bf" alt="" width="563"><figcaption><p>The point-wise error of the solved phase factors.</p></figcaption></figure>

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

<details>

<summary>Output of the code</summary>

```
norm error = 7.64691e-07
max of solution = 0.800546
iter          err
   1  +1.2889e-01 
   2  +8.8228e-03 
   3  +6.4969e-05 
   4  +3.7459e-09 
   5  +1.4694e-15 
Stop criteria satisfied.
The residual error is
   6.2172e-15
```

</details>