Operator Embedding

Suppose we want to multiply a qubit’s amplitude vector by

\[A=\begin{pmatrix}1&0\\0&\tfrac12\end{pmatrix}.\]

On \(\ket0\) it does nothing; on \(\ket1\) it halves the amplitude. A unitary preserves the length of every vector, but \(A\ket1=\frac12\ket1\) has length \(\frac12\). So no unitary on this qubit equals \(A\), even up to a phase: \(A\) is a non-unitary matrix. Worse, \(A\) shrinks different inputs by different amounts: \(A\ket+=\frac{1}{\sqrt2}(\ket0+\frac12\ket1)\) has squared length \(\frac58\). Normalizing the output afterwards does not make a linear operation.

How much can a matrix stretch a vector? The operator norm, which we used to measure errors, answers this:

\[\|A\|=\max_{\|\psi\|=1}\|A\ket\psi\|.\]

For our diagonal \(A\), \(\|A\|=1\), attained on \(\ket0\). For a diagonal matrix it is the largest absolute entry, and for a Hermitian matrix it is the largest absolute eigenvalue. In general, \(\|A\ket\psi\|^2=\bra\psi A^\dagger A\ket\psi\), and \(A^\dagger A\) is Hermitian with nonnegative eigenvalues. Expanding \(\ket\psi\) in its eigenvectors shows that the largest value of \(\bra\psi A^\dagger A\ket\psi\) is its largest eigenvalue, so \(\|A\|\) is the square root of that eigenvalue. A unitary has \(\|U\|=1\) and also shrinks nothing.

The obstacle, then, is that we want to apply a linear map that changes lengths, using only length-preserving operations. The idea that resolves it is to let the missing length go somewhere else: into a branch of a larger system that we later discard or recycle.

Can a larger unitary carry \(A\) as one of its parts? Add one ancilla qubit. Order the basis with the ancilla first, so a two-qubit matrix splits into four \(2\times2\) blocks according to the ancilla’s input and output values. Take

\[U=\begin{pmatrix}A&B\\B&-A\end{pmatrix},\qquad A=\begin{pmatrix}1&0\\0&\tfrac12\end{pmatrix},\quad B=\begin{pmatrix}0&0\\0&\tfrac{\sqrt3}{2}\end{pmatrix}.\]

Because \(A\) and \(B\) are diagonal they commute, and \(A^2+B^2=I\). Multiplying out, \[U^2=\begin{pmatrix}A^2+B^2&AB-BA\\BA-AB&B^2+A^2\end{pmatrix}=I.\] Since \(A\) and \(B\) are real and diagonal, \(U\) is real and symmetric, so \(U^\dagger=U\). Therefore \(U^\dagger U=U^2=I\): \(U\) is unitary, and its own inverse. On an input with the ancilla in \(\ket0\),

\[U\,\ket0\ket\psi=\ket0\,A\ket\psi+\ket1\,B\ket\psi.\]

The ancilla-\(0\) branch carries exactly \(A\ket\psi\). The length that \(A\) removed has gone into the ancilla-\(1\) branch, where \(B\ket\psi\) supplies it: \(\|A\psi\|^2+\|B\psi\|^2=1\).

In general, a unitary \(U\) on \(a\) ancilla qubits and a system is a block encoding of a matrix \(A\) if

\[(\bra{0^a}\otimes I)\,U\,(\ket{0^a}\otimes I)=\frac{A}{\alpha}.\]

The top-left block of \(U\) is \(A/\alpha\). The number \(\alpha\ge\|A\|\) is the subnormalization: a block of a unitary cannot stretch any vector, so it must have norm at most \(1\), and a matrix with larger norm must be scaled down. If the block equals \(A/\alpha\) only to within \(\epsilon\) in operator norm, \(U\) is an \((\alpha,a,\epsilon)\)-block encoding. Every later algorithm that processes matrices assumes block encodings as input, so \(\alpha\), \(a\), \(\epsilon\), and the gate cost of \(U\) are the parameters a resource estimate must report.

How do we get \(A\ket\psi\) out as a state of its own? A block encoding applies \(A/\alpha\) only on the ancilla-\(0\) branch. To obtain \(A\ket\psi\) as a state, measure the ancilla and keep the run if the outcome is \(0\). This is postselection. The probability of success is

\[p_{\text{succ}}=\left\|\frac{A\ket\psi}{\alpha}\right\|^2,\]

and on success the system holds the normalized state \(A\ket\psi/\|A\ket\psi\|\).

In the example with \(\ket\psi=\ket+\), success has probability \(\frac58\) and leaves \(\frac{1}{\sqrt5}(2\ket0+\ket1)\). With input \(\ket1\), success has probability \(\frac14\) and leaves \(\ket1\): the output state is the same as the input, but it was obtained only a quarter of the time. That input dependence is how a non-unitary map shows up in a unitary world: it is paid for in probability.

Without further help, the expected number of attempts is \(1/p_{\text{succ}}=\alpha^2/\|A\psi\|^2\), and each failure destroys the input state, which must be prepared again. If the preparation is a known circuit, amplitude amplification reduces the cost to order \(\alpha/\|A\psi\|\) uses of \(U\), the preparation, and their inverses. Either way, two quantities control the cost: the subnormalization \(\alpha\), and how much \(A\) shrinks the particular input. A block encoding with an unnecessarily large \(\alpha\) wastes success probability quadratically.

Let \(A=\sum_{j=0}^{L-1}a_jU_j\) with unitaries \(U_j\) and coefficients \(a_j>0\) (signs and phases can be absorbed into the \(U_j\)). Put \(\alpha=\sum_ja_j\). How can one circuit apply all the terms with the right weights? Let an ancilla register of \(\lceil\log_2L\rceil\) qubits hold a term label \(j\). One circuit, called PREPARE, loads the weights into its amplitudes:

\[\mathrm{PREP}\ket0=\sum_j\sqrt{a_j/\alpha}\,\ket j.\]

A second circuit, called SELECT, applies the \(j\)th unitary on the branch where the label is \(j\):

\[\mathrm{SEL}=\sum_j\ket j\bra j\otimes U_j.\]

The circuit \(W=(\mathrm{PREP}^\dagger\otimes I)\,\mathrm{SEL}\,(\mathrm{PREP}\otimes I)\) is a block encoding of \(A/\alpha\). To see this, compute the ancilla-\(0\) block:

\[(\bra0\otimes I)W(\ket0\otimes I)=\sum_{j,k}\overline{\braket{k}{\mathrm{PREP}\,0}}\,\braket{j}{\mathrm{PREP}\,0}\,\braket kj\,U_j=\sum_j\frac{a_j}{\alpha}U_j=\frac A\alpha.\]

The ancilla is prepared in a superposition of term labels; SELECT applies every term on its own branch; un-preparing projects back onto the starting state, and each term’s amplitude \(\sqrt{a_j/\alpha}\) appears twice, once going in and once coming out, giving weight \(a_j/\alpha\).

This is the linear combination of unitaries (LCU) method of Childs and Wiebe (2012). Its subnormalization is \(\alpha=\sum_j|a_j|\), the 1-norm of the coefficients, which is at least \(\|A\|\) and can be much larger when terms partly cancel. Its cost is one PREPARE, one PREPARE\(^\dagger\), and one SELECT. SELECT on \(L\) terms is a multiplexed controlled operation; built naively, it costs about \(L\) controlled unitaries plus the overhead of decoding \(j\).

Does the construction really add the terms as amplitudes? Take \(A=\frac12(I+X)\) on one system qubit, so \(a_0=a_1=\frac12\), \(\alpha=1\). PREPARE is a Hadamard gate on one ancilla. SELECT is a controlled-\(X\) from ancilla to system. The block encoding is Hadamard, CNOT, Hadamard on the ancilla.

What to observe

Probability \(\frac12\), leaving the system in \(\ket+\). Indeed \(A\ket0=\frac12(\ket0+\ket1)=\frac{1}{\sqrt2}\ket+\), with squared length \(\frac12\).

Track the state. The first Hadamard on the ancilla gives \(\frac{1}{\sqrt2}(\ket0+\ket1)\ket0\). CNOT gives \(\frac{1}{\sqrt2}(\ket0\ket0+\ket1\ket1)\). The second Hadamard gives \(\frac12\big(\ket0(\ket0+\ket1)+\ket1(\ket0-\ket1)\big)\). The ancilla-\(0\) branch is \(\frac12(\ket0+\ket1)=A\ket0\).

Contrast this with choosing randomly: with probability \(\frac12\) apply \(I\), otherwise apply \(X\). On \(\ket0\) that produces \(\ket0\) or \(\ket1\) at random, a classical mixture, which a measurement in the \(\pm\) basis gives \(+\) only half the time. The LCU output \(\ket+\) gives \(+\) every time. The difference is that LCU adds amplitudes coherently, \(\frac12\ket0+\frac12\ket1\), while random selection adds probabilities. A coherent sum can interfere; a random choice cannot.

Which Hamiltonians can LCU encode? The most common case is a Hamiltonian written as a weighted sum of Pauli strings, tensor products of \(I,X,Y,Z\):

\[H=\sum_jh_jP_j.\]

Can every Hermitian matrix be written this way? Start with one qubit. The products of distinct matrices among \(I,X,Y,Z\) have zero trace, while each squares to \(I\), so \(\operatorname{Tr}(P_aP_b)=2\) if \(a=b\) and \(0\) otherwise. Four mutually orthogonal matrices in the four-dimensional space of \(2\times2\) matrices form a basis, so any \(M\) is \(\sum_ah_aP_a\) with \(h_a=\operatorname{Tr}(P_aM)/2\). If \(M\) is Hermitian, \(h_a^*=\operatorname{Tr}((P_aM)^\dagger)/2=\operatorname{Tr}(MP_a)/2=h_a\), so the coefficients are real. On \(n\) qubits the \(4^n\) Pauli strings satisfy \(\operatorname{Tr}(P_aP_b)=2^n\) or \(0\) in the same way, so the same argument applies. Physical Hamiltonians usually have only polynomially many nonzero terms. Each \(P_j\) is unitary and self-inverse. Write \(h_jP_j=|h_j|\cdot(\mathrm{sign}(h_j)P_j)\), absorbing the sign into the unitary. LCU gives a block encoding of \(H/\lambda\) with

\[\lambda=\sum_j|h_j|,\]

the 1-norm of the coefficients. For example, \(H=0.5\,Z\otimes Z-0.3\,X\otimes I\) has \(\lambda=0.8\). PREPARE maps \(\ket0\) to \(\sqrt{5/8}\ket0+\sqrt{3/8}\ket1\); SELECT applies \(Z\otimes Z\) on ancilla \(0\) and \(-X\otimes I\) on ancilla \(1\).

The 1-norm \(\lambda\) is the quantity that every later algorithm pays for. Simulating \(e^{-iHt}\) with block encodings costs about \(\lambda t\) uses of the encoding, and phase estimation of energies to precision \(\epsilon\) costs about \(\lambda/\epsilon\). For molecules, \(\lambda\) can be far larger than \(\|H\|\), and much research on quantum chemistry algorithms reduces \(\lambda\) by better representations of the same Hamiltonian, for example by rewriting the terms that describe pairs of interacting electrons.

Suppose \(U_A\) is a block encoding of \(A/\alpha\) with \(a\) ancillas and \(U_B\) one of \(B/\beta\) with \(b\) ancillas. Apply \(U_B\), then \(U_A\), each with its own ancilla register. On the branch where both registers are \(0\) at the end,

\[(\bra{0^a}\bra{0^b}\otimes I)\,U_AU_B\,(\ket{0^a}\ket{0^b}\otimes I)=\frac{A}{\alpha}\cdot\frac{B}{\beta}.\]

The reason is that \(U_A\) does not touch the \(b\) register, so the projection onto \(\ket{0^b}\) can be moved past it and applied right after \(U_B\), and similarly for \(\ket{0^a}\). This is a product of block encodings, with subnormalization \(\alpha\beta\) and \(a+b\) ancillas.

Reusing a single ancilla register for both is tempting, since it saves qubits, but it is generally wrong. After \(U_B\), the ancilla is not in \(\ket0\) on the failure branch, and \(U_A\) can map part of that failure amplitude back into the ancilla-\(0\) block, adding unwanted terms to the product. (A later construction, qubitization, does reuse the ancillas, but only with a reflection in between that makes the reuse correct.)

Subnormalizations multiply, so long products lose success probability quickly: a product of \(k\) encodings each with \(\alpha=2\) has subnormalization \(2^k\). Algorithms that apply polynomials of a matrix avoid this with that more economical construction.

How do errors in the blocks affect a product? Real block encodings are approximate: rotation angles are synthesized to finite precision, and PREPARE states are loaded approximately. If the blocks are \(A/\alpha+E\) and \(B/\beta+F\) with \(\|E\|\le\epsilon\) and \(\|F\|\le\delta\), the product block is

\[\Big(\frac A\alpha+E\Big)\Big(\frac B\beta+F\Big)=\frac{AB}{\alpha\beta}+\frac A\alpha F+E\frac B\beta+EF.\]

Since \(\|A/\alpha\|,\|B/\beta\|\le1\), the approximation error is at most \(\delta+\epsilon+\epsilon\delta\). For small errors, errors add, just as they did for products of unitaries.

Two error sources are easy to miss. First, a PREPARE state with amplitude errors changes the encoded coefficients: if each \(\sqrt{a_j/\alpha}\) is off by \(\eta\), each weight \(a_j/\alpha\) is off by about \(2\eta\sqrt{a_j/\alpha}\), and the block error is bounded by summing these over terms. Second, errors relative to \(A/\alpha\) are errors relative to \(A\) multiplied by \(\alpha\): to apply \(A\) itself with error \(\epsilon_A\), the block must be accurate to \(\epsilon_A/\alpha\). Large subnormalizations demand high precision as well as many repetitions.

Postselection fails sometimes, and amplitude amplification fixes that only if we can reflect about the input state, which requires a preparation circuit for it. What if the input is an unknown state handed to us by an earlier part of the algorithm?

When the encoded matrix is proportional to a unitary, there is a way out. Suppose the ancilla-\(0\) block of \(W\) is \(\sin\theta\,V\) for a unitary \(V\) and a fixed angle \(\theta\). Then for every input,

\[W\ket{\Psi}=\sin\theta\,\ket{\Phi}+\cos\theta\,\ket{\Phi^\perp},\qquad\ket\Psi=\ket0\ket\psi,\quad\ket\Phi=\ket0V\ket\psi,\]

where \(\ket{\Phi^\perp}\) has its ancilla orthogonal to \(\ket0\). Amplification would need a reflection about \(\ket\Psi\), which depends on the unknown \(\psi\). Can the reflection \(R=(2\ket0\bra0-I)\otimes I\) about the ancilla being \(0\), which does not depend on \(\psi\), take its place?

To find out, we need the input-side partner of \(\ket{\Phi^\perp}\). Define \(\ket{\Psi^\perp}=W^\dagger\big(\cos\theta\ket\Phi-\sin\theta\ket{\Phi^\perp}\big)\), so that \(W\) maps the pair \(\ket\Psi,\ket{\Psi^\perp}\) to the pair \(\ket\Phi,\ket{\Phi^\perp}\) by a fixed two-by-two matrix. Substituting \(\ket{\Phi^\perp}=(W\ket\Psi-\sin\theta\ket\Phi)/\cos\theta\) gives

\[\ket{\Psi^\perp}=\frac{1}{\cos\theta}\,W^\dagger\ket\Phi-\frac{\sin\theta}{\cos\theta}\,\ket\Psi.\]

What is its ancilla-\(0\) part? The block of \(W^\dagger\) is the adjoint of the block of \(W\), namely \(\sin\theta\,V^\dagger\). So the ancilla-\(0\) part of \(W^\dagger\ket\Phi\) is \(\sin\theta\,V^\dagger V\ket\psi=\sin\theta\ket\psi\), and the ancilla-\(0\) part of \(\ket{\Psi^\perp}\) is \(\frac{\sin\theta}{\cos\theta}\ket\psi-\frac{\sin\theta}{\cos\theta}\ket\psi=0\). So \(R\) fixes \(\ket\Psi\) and \(\ket\Phi\), and negates \(\ket{\Psi^\perp}\) and \(\ket{\Phi^\perp}\).

That is exactly the situation of amplitude amplification, with \(W\) as the preparation and \(\ket\Phi\) as the good state, and with the same angle \(\theta\) for every \(\psi\). The iteration \(-WRW^\dagger R\) therefore rotates by \(2\theta\), and \(k\) iterations after \(W\) give amplitude \(\sin((2k+1)\theta)\) on \(\ket0V\ket\psi\), whatever the input. This is oblivious amplitude amplification (Berry, Childs, Cleve, Kothari, and Somma, 2014).

A clean special case: if \(\sin\theta=\frac12\), one iteration gives exactly \(\ket0V\ket\psi\), since \(3\theta=\pi/2\). LCU makes this easy to arrange. If \(V\approx\sum_ja_jU_j\) with \(\sum_ja_j=2\), the LCU block encoding has \(\sin\theta=\frac12\), and one round of oblivious amplification applies \(V\) deterministically, with three uses of the LCU circuit. How can this simulate \(e^{-iHt}\)? Split the time into short segments of length \(\delta\), and truncate the exponential’s power series after \(K\) terms:

\[e^{-iH\delta}\approx\sum_{k=0}^K\frac{(-iH\delta)^k}{k!}.\]

With \(H=\sum_jh_jP_j\), each power \(H^k\) expands into products of \(k\) Pauli strings, which are unitaries. So the truncated series is a linear combination of unitaries, with coefficient 1-norm \(\sum_{k\le K}(\lambda\delta)^k/k!\le e^{\lambda\delta}\). Choosing \(\delta\) so that this sum is \(2\), about \(\lambda\delta=\ln2\), makes \(\sin\theta=\tfrac12\), and one round of oblivious amplification applies the segment almost deterministically. The truncated sum is only nearly unitary, so a robust version of the argument, which we quote, is needed. This truncated Taylor series method (Berry, Childs, Cleve, Kothari, and Somma, 2015) needs about \(\lambda t/\ln2\) segments, each with \(K\) of order \(\log(\lambda t/\epsilon)\) terms, for a total of order \(\lambda t\log(\lambda t/\epsilon)\) uses of the Pauli terms, up to smaller logarithmic factors.

Suppose a matrix has few nonzero entries per row, and we can compute them. Can we encode it without storing it? A matrix is \(s\)-sparse if each row and column has at most \(s\) nonzero entries. Sparse access provides two oracles: one that, given a row \(i\) and an index \(k\le s\), returns the column of the \(k\)th nonzero entry in that row, and one that, given \(i,j\), returns the entry \(H_{ij}\) to some number of bits. For a matrix defined by a rule, such as a lattice Hamiltonian or a graph’s adjacency matrix, both oracles are efficient circuits.

From these oracles a block encoding of \(H/(s\,\|H\|_{\max})\) can be built with a constant number of oracle calls, a quoted construction (Gilyén, Su, Low, and Wiebe, 2019), where \(\|H\|_{\max}\) is the largest entry magnitude. The construction prepares, for each row, a uniform superposition over its \(s\) nonzero positions, and uses the entry oracle to rotate an ancilla by an angle encoding \(H_{ij}\). The subnormalization \(s\|H\|_{\max}\) is at least \(\|H\|\) and can exceed it by a factor up to \(s\), for example when each row has one large entry and many tiny ones.

Sparse access is the setting of many simulation and linear-systems results. It is natural when the matrix has a short description. It is not free when the matrix is stored as data: then the oracles must read memory coherently.

What if the matrix is an arbitrary \(N\times N\) table of numbers stored classically? A block encoding needs coherent access to it: a circuit that, on a superposition of addresses, returns the stored values in superposition. Such a device is often called QRAM. Every stored value must be able to influence the output, so its circuit size grows at least linearly with the amount of data; known designs reach logarithmic depth (Giovannetti, Lloyd, and Maccone, 2008). Either way, loading the data once costs order \(N^2\) for a dense matrix, which already matches or exceeds many classical algorithms on it.

Claims of exponential speedup for problems on classical data, such as recommendation systems, principal component analysis, or regression, usually assume that the data is already available through a QRAM-like structure in which it can also be sampled efficiently. Tang (2018) showed that when a classical algorithm is given the comparable sampling access to the same data, it can match the quantum recommendation algorithm up to polynomial factors. This dequantization program has since removed the exponential gap for several proposed quantum machine-learning algorithms with low-rank inputs, while leaving large polynomial gaps in some cases.

The lesson is the one we drew when we first specified access: a speedup is a comparison between methods with compatible access. A block encoding is a precise and powerful interface, and it makes clear what an algorithm consumes. Whether supplying that interface is cheap depends on the matrix.