Polynomial Transformations

A block encoding of a matrix puts it in the top-left block of a unitary. What is the smallest case? A one-qubit unitary whose top-left entry is a real number \(x\in[-1,1]\), the unknown signal. What must the other entries be? The first column must have length one, so the entry below \(x\) has magnitude \(\sqrt{1-x^2}\), and the same holds along the first row. A symmetric choice with imaginary off-diagonal entries completes the matrix:

\[W(x)=\begin{pmatrix}x&i\sqrt{1-x^2}\\i\sqrt{1-x^2}&x\end{pmatrix}.\]

Its columns are orthonormal, so it is unitary. Why this completion?

Writing \(x=\cos\theta\), this is \(W=\cos\theta\,I+i\sin\theta\,X=e^{i\theta X}\), a rotation about the \(X\) axis by an angle that encodes \(x\), and its powers are easy to compute. We call \(W(x)\) the signal operator. Its top-left entry, \(\bra0W\ket0=x\), is a one-qubit block encoding of the number \(x\).

What does applying \(W\) twice encode? Multiply the matrices: the top-left entry of \(W^2\) is \(x\cdot x+(i\sqrt{1-x^2})(i\sqrt{1-x^2})=2x^2-1\). It is not \(x^2\). The second application lets amplitude that left \(\ket0\) on the first call return, and it returns with a minus sign. Since \(W^d=e^{id\theta X}\), the top-left entry of \(W^d\) is \(\cos(d\theta)\):

\[\bra0W(x)^d\ket0=T_d(x)=\cos(d\arccos x).\]

These are the Chebyshev polynomials: \(T_0=1\), \(T_1=x\), \(T_2=2x^2-1\), \(T_3=4x^3-3x\), and in general \(T_{d+1}=2xT_d-T_{d-1}\). Each is a genuine polynomial in \(x\) of degree \(d\), bounded by \(1\) on \([-1,1]\), and oscillating between \(\pm1\) there.

So \(d\) uses of a block encoding of \(x\) give a block encoding of a degree-\(d\) polynomial of \(x\). But Chebyshev polynomials are only one family. Can we reach others with the same number of calls?

What else can we do between two signal calls? Insert \(Z\)-rotations \(e^{i\phi Z}=\mathrm{diag}(e^{i\phi},e^{-i\phi})\). For two calls with one phase in the middle,

\[\bra0W(x)\,e^{i\phi Z}\,W(x)\ket0=x^2e^{i\phi}-(1-x^2)e^{-i\phi}=\cos\phi\,(2x^2-1)+i\sin\phi.\]

The phase reshapes the response. At \(\phi=0\) we recover \(T_2\). At \(\phi=\pi/2\) the two paths combine to the constant \(i\): the signal dependence cancels entirely. In between, the real part is \(T_2\) scaled by \(\cos\phi\). With \(d\) signal calls and \(d+1\) phases \(\Phi=(\phi_0,\ldots,\phi_d)\), the sequence

\[U_\Phi(x)=e^{i\phi_0Z}\prod_{k=1}^{d}W(x)\,e^{i\phi_kZ}\]

is called quantum signal processing, and the angles are its phase factors.

Which polynomials can appear as \(P(x)=\bra0U_\Phi(x)\ket0\)? The answer is a theorem, which we quote (Low and Chuang, 2017; Gilyén, Su, Low, and Wiebe, 2019). A polynomial \(P\) is achievable with \(d\) calls if and only if

(i) \(\deg P\le d\); (ii) \(P\) has the parity of \(d\), even or odd; (iii) there is a polynomial \(Q\) of degree at most \(d-1\) and opposite parity with \(|P(x)|^2+(1-x^2)|Q(x)|^2=1\) on \([-1,1]\).

Condition (iii) says that \(P\) can be completed to a unitary; \(Q\) is the off-diagonal entry. In practice one usually wants a real target, and then the condition simplifies: any real polynomial of degree \(d\) with parity \(d\) and \(|P(x)|\le1\) on \([-1,1]\) is the real part of an achievable response. This is part of the same quoted theorem. The real part can be isolated with one extra qubit, by an LCU that averages the sequence with the sequence whose response is the complex conjugate. Phase factors for a given target are computed classically; numerically stable algorithms exist for degrees in the thousands.

The theorem’s conditions shape what algorithms can do. Where do they come from, and what do they forbid? Two matter constantly.

Boundedness. The response is an entry of a unitary, so its magnitude is at most \(1\). It must stay within \([-1,1]\) on the whole interval, not just where we care about it. To approximate an unbounded function like \(1/x\), we must restrict attention to a region where it is bounded after scaling, and accept an arbitrary bounded response elsewhere.

Parity. Why must the response be even or odd? Compare the signal operator at \(x\) and \(-x\). Conjugating by \(Z\) flips the sign of the off-diagonal entries, so \(W(-x)=-Z\,W(x)\,Z\). The \(Z\)-rotations commute with \(Z\), and \(Z^2=I\), so in \(U_\Phi(-x)\) the inner pairs of \(Z\) cancel and

\[U_\Phi(-x)=(-1)^d\,Z\,U_\Phi(x)\,Z.\]

Taking the top-left entry, where \(Z\) contributes \(+1\) on both sides, gives \(P(-x)=(-1)^dP(x)\). Each call flips the parity. This parity constraint means that a function with no symmetry, like \(e^{-ixt}\), must be split into even and odd parts, \(\cos(xt)\) and \(-i\sin(xt)\), each implemented separately and recombined with a linear combination of unitaries, at a constant-factor cost in success probability.

Within these constraints, the design problem is classical polynomial approximation: find a low-degree bounded polynomial close to the desired function on the relevant region. The degree is the number of signal calls, so it is the algorithm’s query cost. Chebyshev polynomials are the natural tool, because they are bounded on \([-1,1]\) and form a basis; truncating a function’s Chebyshev expansion gives near-optimal approximations. A smooth function on \([-1,1]\) needs low degree; a function with a sharp jump, or one that varies rapidly, needs high degree.

What does it mean to apply a function to a matrix? Let \(H\) be Hermitian with spectral decomposition \(H=\sum_j\lambda_j\ket{v_j}\bra{v_j}\). A function of \(H\) is defined eigenvalue by eigenvalue:

\[f(H)=\sum_jf(\lambda_j)\ket{v_j}\bra{v_j}.\]

For a polynomial this agrees with the obvious definition, since \(H^k\ket{v_j}=\lambda_j^k\ket{v_j}\). Applying \(f(H)\) to a state \(\ket\psi=\sum_jc_j\ket{v_j}\) multiplies each component by \(f(\lambda_j)\), without measuring which component is present. An eigenvalue transformation is a circuit that block-encodes \(f(H/\alpha)\) given a block encoding of \(H/\alpha\).

Three choices of \(f\) cover most applications. A step function that is \(1\) on some eigenvalues and \(0\) on others projects onto an eigenspace: a filter. The function \(1/x\) inverts \(H\), solving linear systems. The function \(e^{-ixt}\) gives time evolution. Phase estimation, by contrast, measures eigenvalues one sample at a time. Polynomial transformations process all eigenvalues coherently and keep the superposition.

Suppose \(U\) is a block encoding of \(H/\alpha\) that is also Hermitian, \(U^2=I\), as LCU encodings with Hermitian SELECT are. Define

\[W=R\,U,\qquad R=(2\ket{0^a}\bra{0^a}-I)\otimes I,\]

the block encoding followed by a reflection about the ancilla being zero. What does \(W\) do near one eigenvector? Take \(\ket{v_j}\) with eigenvalue \(\lambda_j\), write \(x_j=\lambda_j/\alpha\), and let \(\ket{0_j}=\ket{0^a}\ket{v_j}\). The ancilla-zero block of \(U\) is \(H/\alpha\), so the part of \(U\ket{0_j}\) with ancilla zero is \(x_j\ket{0_j}\). The rest is orthogonal to every ancilla-zero state:

\[U\ket{0_j}=x_j\ket{0_j}+\sqrt{1-x_j^2}\,\ket{\perp_j},\]

where \(\ket{\perp_j}\) is normalized with its ancilla orthogonal to \(\ket{0^a}\) (when \(|x_j|<1\)). Where does \(U\) send \(\ket{\perp_j}\)? Apply \(U\) to both sides and use \(U^2=I\):

\[\ket{0_j}=x_jU\ket{0_j}+\sqrt{1-x_j^2}\,U\ket{\perp_j},\qquad\text{so}\qquad U\ket{\perp_j}=\sqrt{1-x_j^2}\,\ket{0_j}-x_j\ket{\perp_j}.\]

So the plane spanned by \(\ket{0_j}\) and \(\ket{\perp_j}\) is mapped to itself. In the basis \((\ket{0_j},\ket{\perp_j})\), \(R\) is \(\mathrm{diag}(1,-1)\), and

\[U=\begin{pmatrix}x_j&\sqrt{1-x_j^2}\\\sqrt{1-x_j^2}&-x_j\end{pmatrix},\qquad W=RU=\begin{pmatrix}x_j&\sqrt{1-x_j^2}\\-\sqrt{1-x_j^2}&x_j\end{pmatrix}.\]

This is the same calculation we did for the Grover iterate: \(W\) is a rotation with eigenvalues \(e^{\pm i\theta_j}\), where

\[\cos\theta_j=\frac{\lambda_j}{\alpha}.\]

Changing the phase of \(\ket{\perp_j}\) by \(i\) turns this matrix into the signal operator \(W(x_j)\), and that change commutes with \(Z\)-rotations, which in each plane are what phase rotations about the ancilla-zero state become.

Each eigenvalue of \(H\) has its own two-dimensional “qubit”, which is the origin of the name qubitization (Low and Chuang, 2019). Within each plane, \(W\) behaves exactly like the signal operator \(W(x)\) with \(x=\lambda_j/\alpha\). So interleaving \(W\) with phase rotations on the ancilla implements \(P(\lambda_j/\alpha)\) on every eigenvector at once: a quantum signal processing sequence acting on all planes in parallel.

The relation \(\cos\theta=\lambda/\alpha\) is the same one we met for Szegedy’s walk, and \(W\) is in fact a quantum walk built from the block encoding. Phase estimation on \(W\) estimates \(\theta_j\), hence \(\lambda_j=\alpha\cos\theta_j\), without ever implementing \(e^{-iHt}\). This is the basis of efficient energy estimation for chemistry: the cost per step is one block encoding, and precision \(\epsilon\) in \(\lambda\) needs about \(\alpha/\epsilon\) steps.

What replaces eigenvalues for a matrix that is not Hermitian, or not even square? Every matrix has a singular value decomposition, a linear-algebra result we assume: orthonormal input directions \(\ket{w_j}\), orthonormal output directions \(\ket{u_j}\), and nonnegative stretch factors \(\sigma_j\), with

\[A\ket{w_j}=\sigma_j\ket{u_j}.\]

The \(\sigma_j\) are the singular values. Since \(A^\dagger A\ket{w_j}=\sigma_j^2\ket{w_j}\), their squares are the eigenvalues of \(A^\dagger A\), so the largest is \(\|A\|\), as we found for the operator norm. For a Hermitian matrix they are the absolute eigenvalues. Geometrically, \(A\) rotates the input directions to output directions and stretches each by \(\sigma_j\).

The quantum singular value transformation (QSVT), a quoted theorem of Gilyén, Su, Low, and Wiebe (2019) applies a polynomial to the singular values. Given a block encoding \(U\) of \(A/\alpha\), alternate \(U\) and \(U^\dagger\) with phase rotations controlled by the projectors that define the block. For an odd polynomial \(P\) of degree \(d\), the result block-encodes

\[P^{(SV)}(A/\alpha)=\sum_jP(\sigma_j/\alpha)\ket{u_j}\bra{w_j},\]

which maps input directions to output directions with stretch factors transformed by \(P\). For even \(P\), it maps input directions back to input directions, \(\sum_jP(\sigma_j/\alpha)\ket{w_j}\bra{w_j}\). The parity decides which space the result lives in.

QSVT unifies the earlier techniques. Amplitude amplification is QSVT applied to a \(1\times1\) block, the single number \(\sin\theta\), with the Chebyshev polynomial that maps \(\sin\theta\) toward \(1\). Fixed-point amplification uses a polynomial that approximates a step. Phase estimation, Hamiltonian simulation, and linear-system solving are further polynomial choices. One framework, with its conditions and costs, gives a common language for all of them.

Suppose the eigenvalues of \(H/\alpha\) are known to avoid an interval \((\mu-\Delta,\mu+\Delta)\), and we want to project a state onto the eigenvectors with eigenvalues below \(\mu\). The ideal function is a step: \(1\) below \(\mu\), \(0\) above. A spectral filter is a bounded polynomial that approximates this step within \(\epsilon\) outside the gap. How high must its degree be? A polynomial cannot jump, so we approximate a smooth S-shaped curve whose rise is confined to the gap. A quoted approximation result (Gilyén, Su, Low, and Wiebe, 2019) gives degree

\[d=O\!\left(\frac1\Delta\log\frac1\epsilon\right).\]

A narrower gap requires a steeper polynomial and higher degree; higher precision costs only logarithmically.

Filters prepare ground states. If a trial state has overlap \(\gamma\) with the ground state and the ground energy is separated from the rest by an energy gap, filtering below the midpoint of the gap leaves the ground state with success amplitude about \(\gamma\). With amplitude amplification the total cost is about \(\frac1\gamma\cdot\frac{\alpha}{\Delta_E}\log\frac1\epsilon\) block encodings, where \(\Delta_E\) is the energy gap (Lin and Tong, 2020). The same two quantities appear as in phase estimation on a superposition of eigenstates: the overlap of the initial state, and the spectral resolution required. Filtering does not remove the need for a good trial state; it uses it more efficiently than repeated phase estimation.

The quantum linear systems problem: given a block encoding of an invertible matrix \(A\) and a preparation circuit for \(\ket b\), produce a state proportional to \(A^{-1}\ket b\). Assume \(A\) has singular values in \([1/\kappa,1]\) after scaling; \(\kappa\) is the condition number, the ratio of largest to smallest singular value.

The function \(1/x\) is unbounded, so approximate the scaled function \(\frac{1}{\kappa x}\) on \(1/\kappa\le|x|\le1\), where it lies in \([1/\kappa,1]\), and let the polynomial do anything bounded on \(|x|<1/\kappa\). A quoted approximation result (Childs, Kothari, and Somma, 2017) gives an odd polynomial of degree \(O(\kappa\log(\kappa/\epsilon))\). Applying it gives a block encoding of \(A^{-1}/\kappa\) (roughly), and postselection succeeds with probability \(\|A^{-1}b\|^2/\kappa^2\), at least \(1/\kappa^2\). Amplitude amplification brings the number of repetitions to order \(\kappa\), for a total of order \(\kappa^2\log(\kappa/\epsilon)\) block-encoding calls; more refined methods reach the optimal order \(\kappa\log(1/\epsilon)\).

The first algorithm of this kind was Harrow, Hassidim, and Lloyd (2009), using phase estimation; polynomial methods improved its dependence on precision. The output is a quantum state, and the caveats about state outputs apply in full. Reading all \(N\) entries of the solution would cost order \(N\) repetitions. The algorithm is useful when the next step needs only a few numbers, such as an expectation value \(\bra xM\ket x\), and when \(A\) has an efficient block encoding with small \(\kappa\) and \(\ket b\) is cheap to prepare. Aaronson’s “read the fine print” (Nature Physics, 2015) lists these conditions; each one can erase the advantage.

To simulate \(e^{-iHt}\), apply \(f(x)=e^{-i\alpha tx}\) to the block encoding of \(H/\alpha\). Split it into \(\cos(\alpha tx)\) and \(-i\sin(\alpha tx)\), each bounded and of definite parity. Their Chebyshev expansions are known in closed form, the Jacobi–Anger expansion:

\[\cos(\tau x)=J_0(\tau)+2\sum_{k\ge1}(-1)^kJ_{2k}(\tau)T_{2k}(x),\qquad\sin(\tau x)=2\sum_{k\ge0}(-1)^kJ_{2k+1}(\tau)T_{2k+1}(x),\]

where the coefficients \(J_k(\tau)\) are standard named functions, the Bessel functions. We need only one property, which we quote: \(J_k(\tau)\) becomes tiny once \(k\) exceeds \(\tau\), and then decays faster than exponentially. For \(\tau=10\), \(|J_{10}|\approx0.21\), \(|J_{14}|\approx0.012\), \(|J_{18}|\approx1.5\times10^{-4}\) and \(|J_{25}|\approx7\times10^{-9}\). Truncating at degree \(d\approx e\tau/2+\log(1/\epsilon)\) gives error \(\epsilon\), again as a quoted bound. With \(\tau=\alpha t\),

\[d=O\!\left(\alpha t+\log\frac1\epsilon\right).\]

This additive dependence, time plus log precision rather than a product, is optimal: quoted lower bounds show that simulation for time \(t\) needs order \(\alpha t\) queries, and order \(\log(1/\epsilon)/\log\log(1/\epsilon)\) for precision. Qubitization with QSP achieves it with only a few ancillas (Low and Chuang). The cost is stated in block-encoding calls; each call’s gate count, the subnormalization \(\alpha\), and the phase precision must be added for a full count. This optimal simulation must still be compared with simpler methods that use no block encodings at all.