# Singular value threshold projector

In this example, We introduce the implementation of singular value threshold projector.

### What is singular value threshold projector

Singular value threshold projectors project out singular vectors with singular values below/above a certain threshold. Suppose that we have access to a unitary $$U$$, its inverse $$U^\dagger$$ and the controlled reflection operators $$(2\Pi-I)$$, $$(2\tilde{\Pi}-I)$$.  $$A:=\tilde{\Pi} U\Pi$$ is of our interest and has a singular value decomposition $$A=W\Sigma V$$. For $$S\subset R$$ let $$\Sigma\_S$$ be the matrix obtained from $$\Sigma$$ by replacing all diagonal entries $$\Sigma\_{ii}\in S$$ by 1 and all diagonal entries $$\Sigma\_{ii}\notin S$$ by 0. We define $$\Pi\_S:=\Pi V\Sigma\_S V^{\dagger}\Pi$$ and similarly $$\tilde{\Pi}\_S:=\tilde{\Pi} W\Sigma\_S W^{\dagger}\tilde{\Pi}$$. For example, set $$S=\[0,\delta]$$ and $$\Pi\_S$$ projects out right singular vectors with singular value at most $$\delta$$. These threshold projectors play a major role in quantum algorithms.

### The equivalent problem in function approximation

As an example, we implement $$\Pi\_{\[0,0.5]}$$. We need to implement singular value transformation of $$A$$ for the following rectangle function. We call a subroutine to find the best even polynomial approximating $$f(x)$$ on the interval $$D\_{\delta}=\[0, 0.5-\delta]\cup \[0.5+\delta,1].$$ We solve the problem by convex optimization. Here are the parameters set for the subroutine. 

### Set up parameters

```matlab
opts.intervals=[0,0.5-delta,0.5+delta,1];
opts.objnorm = Inf;
opts.epsil = 0.1;
opts.npts = 500;
opts.isplot= true;
% scaling factor for the infinity norm
opts.fscale = 0.9;
opts.maxiter = 100;
opts.criteria = 1e-12;
opts.useReal = false;
opts.targetPre = true;
opts.method = 'Newton';
```

### Solving polynomial approximation

```matlab
targ = @(x) rectangularPulse(-0.5,0.5,x);
% Compute its Chebyshev coefficients. 
parity = 0;
deg = 250;
coef_full=cvx_poly_coef(targ, deg, opts);
% The solver outputs all Chebyshev coefficients while we have to post-select 
% those of odd order due to the parity constraint.
coef = coef_full(1+parity:2:end);
```

<figure><img src="https://40428858-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FznDlSgNqsJLx79OHu47P%2Fuploads%2FHeskwwvMCPvDMCR9MiMy%2Fsingular_value_threshold_projectors_polynomial.png?alt=media&#x26;token=f7c92609-083b-4e2f-935d-8ab22b003747" alt=""><figcaption><p>Polynomial approximation of the target function</p></figcaption></figure>

<figure><img src="https://40428858-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FznDlSgNqsJLx79OHu47P%2Fuploads%2FhiZvaM0dZ2qkm9lZophc%2Fsingular_value_threshold_projectors_polynomial_error.png?alt=media&#x26;token=5f40cc66-a14e-4926-a5f4-0715cc2cbb15" alt=""><figcaption><p>Polynomial approximation error</p></figcaption></figure>

### Solving phase factors by running the solver

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

### Verifying the solution

```matlab
xlist1 = linspace(0,0.5-delta,500)';
xlist2 = linspace(0.5+delta,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,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,'singular_value_threshold_projectors_error.png','-dpng','-r500');
```

<figure><img src="https://40428858-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FznDlSgNqsJLx79OHu47P%2Fuploads%2F8I7VM7Oh9oKxIeyA5fLZ%2Fsingular_value_threshold_projectors_error.png?alt=media&#x26;token=b699d95c-fa43-4943-bcd4-4028c146b0ed" alt=""><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 = 6.5514e-08
max of solution = 0.9
iter          err
   1  +4.8369e-01 
   2  +1.9791e-01 
   3  +2.7179e-02 
   4  +3.0099e-04 
   5  +3.2263e-08 
Stop criteria satisfied.
The residual error is
   3.3529e-14
```

</details>