Quantum Algorithms · Chapter 09
Polynomial transformations gave one way to build \(e^{-iHt}\): encode \(H\) in a block, qubitize it, and apply a polynomial approximation of the exponential, with degree close to \(\lambda t\). Simulating quantum systems was the first reason anyone wanted a quantum computer, and it remains the application with the clearest case for an advantage. A classical computer tracking \(n\) qubits needs \(2^n\) numbers; a quantum computer needs \(n\) qubits and a circuit for \(e^{-iHt}\).
Block encodings cost ancillas and elaborate PREPARE and SELECT circuits. A Hamiltonian is usually a sum of simple terms, and each term on its own is easy to exponentiate. Can we build the whole evolution from short evolutions under the terms, and how large is the error when the terms do not commute?
We find that alternating short evolutions is exact for commuting terms and errs by an amount set by the commutator otherwise, which tells us how many steps a given accuracy needs. Making the alternation symmetric cancels the leading error, and recursion gives formulas of any even order. Each short evolution then becomes a small circuit, which lets us compare product formulas fairly with polynomial methods and list the errors that lie outside the circuit altogether. Finally we change the Hamiltonian itself, slowly, to carry an easy ground state into a useful one, and find that success depends on an energy gap that is rarely known.
When the gap is unknown, the remaining option is to give up guarantees and tune a short circuit by measurement.
How can simple available evolutions reproduce a desired dynamics?
Evolution can be assembled from short steps whose error is set by commutators, or steered by slow change whose success is set by a gap.
Engineer time evolution: build e^{−iHt} from short evolutions of simple terms with controlled error, compare simulation strategies fairly, and use slowly changing Hamiltonians to prepare ground states.
Use matrix exponentials and spectral decompositions (Quantum Walks), operator-norm error composition (Phase Encoding), Pauli-sum block encodings (Operator Embedding), and polynomial simulation (Polynomial Transformations).
H = X + Z: first-order Trotter with r steps has error about 0.7/r against the bound 1/r, and the adiabatic path from −X to −Z succeeds once the total time is a few units.
Error bounds are worst-case operator-norm bounds for bounded generators. Hardware noise is not modelled. Adiabatic runtimes depend on gaps that are rarely known for large instances.
9.1
What exactly must a simulation deliver, and why should it be hard? Once the task is stated, the obstacle is that the terms of a Hamiltonian do not commute. Comparing short evolutions to second order shows that the commutator sets the error. Bounding that error over a whole run tells us how many steps are needed, and arranging the steps symmetrically reduces the count.
Hamiltonian simulationno fast-forwarding
What should a simulation take in and hand back? The input to Hamiltonian simulation is a Hamiltonian on \(n\) qubits, a time \(t\), and an error tolerance \(\epsilon\). The output is a circuit \(V\) with \(\|V-e^{-iHt}\|\le\epsilon\). The Hamiltonian is usually given as a sum of terms that are individually simple,
\[H=\sum_{j=1}^{L}H_j,\]
for example Pauli strings with coefficients, or a lattice model where each term acts on two neighbouring sites. Each \(e^{-iH_j\tau}\) is easy to implement exactly, as we will see when we turn Pauli exponentials into circuits. The whole \(e^{-iHt}\) is not, because the terms do not commute.
How must the cost grow with \(t\)? Generically at least linearly. If an arbitrary Hamiltonian could be simulated for time \(t\) with far fewer than \(t\) steps, one could compute things known to require many queries, such as parities of many bits, by encoding them in a Hamiltonian and evolving. This no fast-forwarding theorem (Berry, Ahokas, Cleve, and Sanders, 2007) says that long-time dynamics cannot in general be shortcut. Some special Hamiltonians, those that a classical change of basis diagonalizes, can be fast-forwarded; that is the exception, and its absence is what makes simulation hard classically and valuable quantumly.
The comparison with classical computation is the reason for interest. A classical computer storing the state of \(n\) qubits needs \(2^n\) amplitudes. For many physical systems, no known classical method avoids this exponential cost at the accuracy required, while a quantum simulation needs \(n\) qubits and a gate count polynomial in \(n\), \(t\), and \(1/\epsilon\). This is a comparison with known classical methods, which exploit structure such as low entanglement or approximate symmetry, and it must be made for each target system.
A lattice of \(100\) spins has a Hamiltonian with nearest-neighbour terms on a line. How many terms \(L\), and how many amplitudes would a classical state vector store?
\(L=99\) bond terms (more if each bond has several Pauli terms, for example \(XX+YY+ZZ\) gives \(297\)). The state vector has \(2^{100}\approx1.3\times10^{30}\) amplitudes, far beyond any classical memory.
Give a Hamiltonian that can be fast-forwarded, and say why.
\(H=\sum_jZ_j\). It is diagonal in the computational basis, so \(e^{-iHt}=\bigotimes_je^{-iZ_jt}\) for any \(t\), with \(n\) single-qubit rotations regardless of \(t\). Any Hamiltonian made of commuting terms diagonalized by a known efficient circuit can be fast-forwarded.
Why is the error measured in operator norm rather than as an error in one output state?
The simulation circuit may be applied to any input state, and often to states prepared by other parts of a larger algorithm. Operator norm bounds the error for every input at once, and errors in composed operations add, so a total budget can be split among components.
Why does non-commutation matter? Start with two terms, \(H=A+B\). If \(A\) and \(B\) commute, \(AB=BA\), then the exponential factors exactly: \(e^{-i(A+B)t}=e^{-iAt}e^{-iBt}\), just as for numbers. The proof expands both sides as power series; commutation lets the terms be rearranged freely.
If they do not commute, the factorization fails. Take \(A=X\), \(B=Z\). Because \(XZ=-ZX\), the cross terms in \((X+Z)^2=X^2+XZ+ZX+Z^2\) cancel, and \((X+Z)^2=2I\). So \(X+Z\) behaves like \(\sqrt2\) times a matrix squaring to \(I\), and the power-series argument we used for the two-vertex walk gives
\[e^{-it(X+Z)}=\cos(\sqrt2t)\,I-i\,\frac{\sin(\sqrt2t)}{\sqrt2}\,(X+Z).\]
This is not \(e^{-itX}e^{-itZ}\), which is a rotation about a different axis by a different angle.
How big is the difference for a short time \(\delta\)? Expand both to second order:
\[e^{-i\delta A}e^{-i\delta B}=I-i\delta(A+B)-\tfrac{\delta^2}{2}(A^2+2AB+B^2)+O(\delta^3),\]
\[e^{-i\delta(A+B)}=I-i\delta(A+B)-\tfrac{\delta^2}{2}(A^2+AB+BA+B^2)+O(\delta^3).\]
They agree to first order and differ at second order by \(-\frac{\delta^2}{2}(AB-BA)\). The commutator \([A,B]=AB-BA\) measures exactly how far the terms are from commuting, and it controls the error. Splitting \(e^{-iHt}\) into \(r\) short steps, each approximated by a product of term exponentials, is a product formula, also called Trotterization.
Compute \([X,Z]\) and its operator norm.
\(XZ=\begin{pmatrix}0&-1\\1&0\end{pmatrix}\) and \(ZX=\begin{pmatrix}0&1\\-1&0\end{pmatrix}\), so \([X,Z]=\begin{pmatrix}0&-2\\2&0\end{pmatrix}=-2iY\). Its norm is \(2\).
Do \(Z_1Z_2\) and \(X_1X_2\) commute? Use the rule that two Pauli strings commute when they anticommute on an even number of positions.
They differ (anticommute) on both positions 1 and 2, an even number, so they commute: \[Z_1Z_2\,X_1X_2=(Z_1X_1)(Z_2X_2)=(-X_1Z_1)(-X_2Z_2)=X_1X_2\,Z_1Z_2\].
Verify \(e^{-it(X+Z)}\) at \(t=\pi/(2\sqrt2)\). What gate is it, up to phase?
\(\cos(\pi/2)=0\), \(\sin(\pi/2)=1\), so \(e^{-it(X+Z)}=-i(X+Z)/\sqrt2=-iH\), the Hadamard gate up to the global phase \(-i\).
The first-order product formula uses \(r\) steps of length \(\delta=t/r\):
\[S_1(t)=\big(e^{-i\delta A}e^{-i\delta B}\big)^r.\]
How large is the error of one step, exactly rather than to second order? Let \(U(\tau)=e^{-i\tau(A+B)}\) and compare it with the product through
\[G(\tau)=U(\tau)^\dagger\,e^{-i\tau A}e^{-i\tau B},\qquad G(0)=I.\]
Since \(U\) is unitary, the error of one step is \(\|G(\delta)-I\|\). We use two facts from calculus: differentiating the power series term by term gives \(\frac{d}{d\tau}e^{i\tau M}=iMe^{i\tau M}\), and an operator changes over an interval by at most the interval’s length times its largest rate of change. Differentiating \(G\), the terms in \(A\) cancel and
\[G'(\tau)=i\,U(\tau)^\dagger\big(B-e^{-i\tau A}Be^{i\tau A}\big)e^{-i\tau A}e^{-i\tau B}.\]
The bracket is how much \(B\) changes when conjugated by \(e^{-i\tau A}\). Its rate of change in \(\tau\) is \(e^{-i\tau A}(-i[A,B])e^{i\tau A}\), of norm \(\|[A,B]\|\), so the bracket has norm at most \(\tau\|[A,B]\|\). The unitaries around it do not change norms, so \(\|G'(\tau)\|\le\tau\|[A,B]\|\), and integrating from \(0\) to \(\delta\),
\[\big\|e^{-i\delta A}e^{-i\delta B}-e^{-i\delta(A+B)}\big\|\le\frac{\delta^2}{2}\|[A,B]\|.\]
Errors of \(r\) steps add, as for any product of unitaries. So the Trotter error satisfies
\[\big\|S_1(t)-e^{-i(A+B)t}\big\|\le r\cdot\frac{\delta^2}{2}\|[A,B]\|=\frac{t^2\,\|[A,B]\|}{2r}.\]
To reach error \(\epsilon\), choose \(r\ge t^2\|[A,B]\|/(2\epsilon)\) Trotter steps. For \(A=X\), \(B=Z\), \(\|[A,B]\|=2\) and the bound is \(t^2/r\).
How tight is the bound? For \(t=1\) it gives \(1/r\), while the actual error is close to \(0.70/r\): the right shape, about 40% loose. That is typical. Worst-case commutator bounds are often loose by constant or even polynomial factors for physical systems, and tighter bounds that use the structure of all the nested commutators (Childs, Su, Tran, Wiebe, and Zhu, 2021) or empirical error estimates can reduce the step count substantially. For a general sum of \(L\) terms, the first-order error bound is \(\frac{t^2}{2r}\sum_{j<k}\|[H_j,H_k]\|\): only noncommuting pairs contribute.
For \(H=X+Z\), \(t=1\), how many first-order steps does the bound require for \(\epsilon=0.01\)? About how many suffice in fact?
The bound gives \(r\ge1/0.01=100\). The actual error is about \(0.7/r\), so about \(70\) steps suffice. Each step uses two rotations.
Doubling the simulation time \(t\) at fixed \(\epsilon\) multiplies the number of first-order steps by what factor? The total gate count?
\(r\propto t^2\), so four times the steps, and four times the gates. First-order formulas scale quadratically in time, worse than the linear lower bound; higher-order formulas approach linear scaling.
If \(A\) and \(B\) commute, what does the bound say, and what does that mean for a Hamiltonian whose terms all commute?
The bound is zero: one step is exact for any \(t\). A Hamiltonian of commuting terms can be simulated by applying each term’s exponential once for the full time, a fast-forwardable case.
symmetric Trotter formulaSuzuki formulas
Can we remove the second-order error? It came from an asymmetry: \(A\) always acts before \(B\). So symmetrize. The symmetric Trotter formula (also called Strang splitting) uses
\[S_2(\delta)=e^{-i\delta A/2}\,e^{-i\delta B}\,e^{-i\delta A/2}.\]
Expanding each factor to second order and collecting the \(\delta^2\) terms, the three factors give \(A^2/8+A^2/8+A^2/4=A^2/2\), \(B^2/2\), \(AB/2\) and \(BA/2\), so
\[S_2(\delta)=I-i\delta(A+B)-\tfrac{\delta^2}{2}\big(A^2+AB+BA+B^2\big)+O(\delta^3),\]
which matches \(e^{-i\delta(A+B)}\) through second order. The commutator terms have cancelled, and the first difference from \(e^{-i\delta(A+B)}\) appears at third order, with nested commutators \([A,[A,B]]\) and \([B,[B,A]]\). Over \(r\) steps the error is \(O(t^3/r^2)\). In the X+Z example, \(r=10\) symmetric steps give error about \(0.003\), against \(0.07\) for first order.
The cost per step is almost the same, since adjacent half-steps \(e^{-i\delta A/2}e^{-i\delta A/2}\) merge into one. Suzuki showed (Suzuki, 1990), as a result we quote, how to build formulas of any even order \(2k\) recursively from \(S_2\) (the Suzuki formulas), with error \(O(t^{2k+1}/r^{2k})\). The number of steps for error \(\epsilon\) is then about \(t^{1+1/(2k)}/\epsilon^{1/(2k)}\), approaching linear in \(t\) as \(k\) grows. The price is that each step of order \(2k\) contains \(5^{k-1}\) symmetric steps, so very high orders are wasteful; in practice orders two to six are used.
For the symmetric formula, estimate how many steps are needed for \(\epsilon=10^{-4}\) at \(t=1\), using the figure’s values.
The error at \(r=50\) is about \(1.1\times10^{-4}\) and falls as \(1/r^2\), so \(r\approx55\) steps. First order would need about \(7{,}000\).
Why does the symmetric formula cost little more than the first-order one over many steps?
In \(S_2(\delta)^r\), the last half-step \(e^{-i\delta A/2}\) of one step and the first of the next combine into \(e^{-i\delta A}\). So \(r\) symmetric steps use \(r+1\) exponentials of \(A\) and \(r\) of \(B\), about the same as first order.
A fourth-order Suzuki formula has error \(O(t^5/r^4)\) and uses five symmetric steps per step. For \(t=100\), compare the scaling of required steps with second order, ignoring constants.
Second order: \(r\propto t^{3/2}/\epsilon^{1/2}=1000/\sqrt\epsilon\). Fourth order: \(r\propto t^{5/4}/\epsilon^{1/4}\approx316/\epsilon^{1/4}\), times five for the cost per step. For \(\epsilon=10^{-6}\): second order about \(10^6\), fourth order about \(5\times316\times32\approx5\times10^4\). Higher order wins for long times and high precision, until constants take over.
9.2
A product formula is a recipe in terms of term exponentials. How many gates does it actually cost, and is it better than polynomial simulation? Each term exponential turns out to be a short parity circuit, which gives a gate count. With gate counts on both sides, a comparison is fair only under the same error and access assumptions. Even then the circuit’s error is one item in a longer budget.
Pauli exponentialparity circuit
How do we implement one term exponential? Every term exponential in a Pauli-sum Hamiltonian has the form \(e^{-i\theta P}\) for a Pauli string \(P\). Start with \(P=Z\otimes Z\). On a basis state \(\ket{b_1b_2}\), \(Z\otimes Z\) has eigenvalue \((-1)^{b_1\oplus b_2}\), so \(e^{-i\theta Z\otimes Z}\) multiplies by \(e^{-i\theta(-1)^{b_1\oplus b_2}}\): a phase depending only on the parity. The circuit computes the parity into the second qubit with a CNOT, applies \(R_z(2\theta)=e^{-i\theta Z}\) there, and uncomputes the parity with a second CNOT. This is a parity circuit: compute, use, uncompute once more.
For a longer string such as \(Z_1Z_2Z_3Z_4\), a ladder of CNOTs computes the parity of all four bits into one qubit, and its mirror image uncomputes it. For strings containing \(X\) or \(Y\), first rotate those qubits so the Pauli becomes \(Z\): the Hadamard gate maps \(X\) to \(Z\), and \(S^\dagger\) followed by the Hadamard gate maps \(Y\) to \(Z\). Then apply the \(Z\)-string circuit and undo the basis changes.
A Pauli exponential on a string of weight \(w\) (the number of non-identity positions) therefore costs \(2(w-1)\) CNOTs, one rotation, and up to \(2w\) single-qubit basis changes. A first-order step for a Hamiltonian with \(L\) terms of weight at most \(w\) costs \(O(Lw)\) gates. In fault-tolerant settings, each arbitrary rotation \(R_z\) must itself be synthesized from \(T\) gates and Clifford gates, the gates generated by the Hadamard gate, \(S\) and CNOT. A quoted synthesis result (Ross and Selinger, 2016) needs about \(3\log_2(1/\epsilon_{\mathrm{rot}})\) \(T\) gates for rotation error \(\epsilon_{\mathrm{rot}}\), which adds to the total error budget.
Verify that CNOT, \(R_z(2\theta)\) on the target, CNOT implements \(e^{-i\theta Z\otimes Z}\) on \(\ket{01}\).
CNOT leaves \(\ket{01}\) unchanged (control is 0). \(R_z(2\theta)=\mathrm{diag}(e^{-i\theta},e^{i\theta})\) on the target, which is \(1\), gives \(e^{i\theta}\). CNOT again leaves it. The phase \(e^{i\theta}=e^{-i\theta\cdot(-1)}\) matches \(Z\otimes Z\ket{01}=-\ket{01}\).
Give the circuit for \(e^{-i\theta X\otimes Y}\).
Apply a Hadamard gate to qubit 1, and \(S^\dagger\) then a Hadamard gate to qubit 2 (mapping \(X\otimes Y\) to \(Z\otimes Z\)), then CNOT, \(R_z(2\theta)\), CNOT, then undo: a Hadamard gate on qubit 1, and a Hadamard gate then \(S\) on qubit 2.
A Hamiltonian has \(L=1000\) Pauli terms of weight at most \(4\). Estimate the CNOT count of \(r=100\) first-order Trotter steps.
Each term needs at most \(2(4-1)=6\) CNOTs, so a step needs at most \(6000\) and \(100\) steps need about \(6\times10^5\) CNOTs, plus \(10^5\) rotations. Ordering terms so that adjacent CNOT ladders cancel can reduce this substantially.
Which simulation method should we use? We now have three families, each with a cost formula. Product formulas cost about \(L\cdot r\) term exponentials, with \(r\) depending on commutators, time, and order. Linear combinations of unitaries with truncated Taylor series cost about \(\lambda t\log(\lambda t/\epsilon)\) uses of a block encoding. Qubitization with signal processing costs about \(\lambda t+\log(1/\epsilon)\) uses of a block encoding, where \(\lambda\) is the 1-norm and each block encoding uses a PREPARE and SELECT over all \(L\) terms.
Asymptotically, qubitization is optimal. In practice the comparison depends on constants and structure. Product formulas need no ancillas, have simple circuits, and benefit from commutator structure: for lattice models with local interactions, their errors can be far smaller than worst-case bounds suggest. Block-encoding methods pay for PREPARE and SELECT circuits, but their error dependence on precision is exponentially better. An empirical comparison for spin systems (Childs, Maslov, Nam, Ross, and Su, 2018) found that high-order product formulas with tight error estimates used the fewest gates for the instances studied, while qubitization was preferable for very high precision. For quantum chemistry, where \(L\) is large and the terms are nonlocal, block-encoding methods applied to rewritten Hamiltonians with smaller 1-norm currently give the lowest resource estimates.
A fair comparison fixes the same Hamiltonian, the same time, the same error in the same norm, and the same gate set, and counts everything: ancillas, rotation synthesis, and classical precomputation. Comparing the degree of a polynomial with a number of Trotter steps is not a comparison of costs.
For \(\lambda t=1000\) and \(\epsilon=10^{-3}\), compare the number of block-encoding calls for Taylor-series LCU and qubitization, ignoring constants.
Taylor: about \(1000\times\log(10^6)\approx1000\times14=1.4\times10^4\). Qubitization: about \(1000+7\approx10^3\). Both must be multiplied by the cost of one block encoding.
Why might a product formula beat qubitization for a one-dimensional lattice model with \(n=100\) sites, despite worse asymptotic scaling?
The terms are local, so most pairs commute and commutator bounds are small; tight error estimates allow few steps. Each step is shallow and needs no ancillas. Qubitization needs PREPARE and SELECT over all terms in every call, with ancilla registers and controlled Paulis, a large constant overhead per call.
A paper reports “simulation with \(O(\log(1/\epsilon))\) dependence on precision.” What must be checked before comparing it with a product formula?
The access model (block encoding or sparse oracles) and the gate cost of one access; the subnormalization \(\lambda\), which multiplies the time; whether rotation-synthesis and phase-factor precision are included; ancilla count; and whether both methods are compared at the same error for the same Hamiltonian, not just in asymptotic exponents.
model errorobservable estimation
A simulation circuit with operator-norm error \(\epsilon\) is one source of error among several. Where do the others enter? Follow one estimate, the ground energy of a molecule, from the physical question to the reported number.
First we must write down a Hamiltonian. It is a model of the molecule, with a truncated set of orbitals or a discretized space, and its distance from the true system is model error. No circuit can remove it. Next we approximate the evolution under that model by a product formula or polynomial, which is the algorithmic error bounded above. Each rotation in the circuit is then synthesized from a finite gate set, which adds synthesis error. The hardware running the circuit is noisy, and that error is mitigated or corrected at further cost. Finally the output is an expectation value or phase estimated from repeated measurements, with statistical error: order \(1/\delta^2\) runs for precision \(\delta\) by sampling, or order \(1/\delta\) with amplitude or phase estimation. A resource estimate must say how the total is split among these five.
What does an error \(\epsilon\) in the evolution mean for the number we report? Observable estimation answers this. If the goal is an expectation value \(\bra\psi O\ket\psi\) with \(\|O\|\le1\), an operator-norm error \(\epsilon\) in the evolution changes it by at most \(2\epsilon\). Energies are estimated by phase estimation on \(e^{-iHt}\) or on the qubitized walk, which needs an initial state with large overlap on the ground state; the evolution error and the phase-estimation resolution must both be budgeted. Rescaling \(H\) to shorten \(t\) changes none of this, since only \(\lambda t\) and \(\epsilon\) enter the costs.
A target energy accuracy is \(10^{-3}\) in units where \(\|H\|\approx10\). The model error is estimated at \(5\times10^{-4}\). How much remains for the algorithm and for statistics?
About \(5\times10^{-4}\) in total if errors add, to be split between algorithmic, synthesis, and statistical error, for example \(2\times10^{-4}\) each with some margin. Spending the whole budget on circuit precision while ignoring the model error would be wasted effort.
Show that \(|\bra\psi V^\dagger OV\ket\psi-\bra\psi U^\dagger OU\ket\psi|\le2\|O\|\,\|V-U\|\) for unitaries \(U,V\).
Write \(V=U+E\). The difference is \(\bra\psi(U^\dagger OE+E^\dagger OU+E^\dagger OE)\ket\psi\). Bounding with \(\|U\|=1\) gives \(2\|O\|\|E\|+\|O\|\|E\|^2\). A slightly tighter argument, telescoping \(V^\dagger OV-U^\dagger OU=V^\dagger O(V-U)+(V-U)^\dagger OU\), gives exactly \(2\|O\|\|V-U\|\).
Why does a quantum simulation of a molecule not automatically give chemically accurate energies?
Chemical accuracy (about \(1.6\times10^{-3}\) hartree) requires the model Hamiltonian itself, including the basis set and active space, to be accurate to that level, and an initial state with good overlap with the ground state, in addition to a precise algorithm. The quantum algorithm solves the model; it does not fix the model.
9.3
So far the Hamiltonian has been fixed and we have simulated its dynamics. A different strategy changes the Hamiltonian itself, slowly, from one whose ground state is easy to prepare to one whose ground state encodes an answer. If the change is slow enough, the state follows. How slow is slow enough? A one-qubit example shows the state following and identifies the gap as what protects it. Simulating that example at several speeds shows how the required time grows as the gap shrinks. Applied to optimization, the same condition explains both the promise and the limits of the method.
adiabatic evolutionground state
Can a slowly changing Hamiltonian carry one ground state into another? Consider the one-qubit family
\[H(s)=-(1-s)\,X-s\,Z,\qquad0\le s\le1.\]
At \(s=0\), \(H=-X\), whose ground state, the eigenvector of lowest energy, is \(\ket+\). At \(s=1\), \(H=-Z\), with ground state \(\ket0\). In between, \(H(s)\) has eigenvalues \(\pm\sqrt{(1-s)^2+s^2}\), and its ground state rotates smoothly from \(\ket+\) to \(\ket0\). Adiabatic evolution prepares \(\ket+\), then evolves under \(H(s)\) with \(s\) increasing from \(0\) to \(1\) over a total time \(T\). If \(T\) is large, the state stays close to the instantaneous ground state and ends near \(\ket0\).
The separation between the two eigenvalues is the gap, \(g(s)=2\sqrt{(1-s)^2+s^2}\), smallest at \(s=\frac12\) with value \(\sqrt2\). The gap is what protects the state: transitions to the excited state are suppressed when the Hamiltonian changes slowly compared with the energy scale set by the gap. A gap is an energy, not a probability of error; the error depends on how the gap compares with the rate of change.
Verify that \(H(s)^2=\big((1-s)^2+s^2\big)I\), and hence the eigenvalues.
\[H^2=(1-s)^2X^2+s^2Z^2+s(1-s)(XZ+ZX)=\big((1-s)^2+s^2\big)I\] because \(XZ=-ZX\). So the eigenvalues are \(\pm\sqrt{(1-s)^2+s^2}\).
If the evolution is instantaneous (\(T=0\)), what is the probability of ending in \(\ket0\)?
The state stays \(\ket+\), and \(|\braket0+|^2=\frac12\). Too fast, and the state does not follow the Hamiltonian at all.
Why does the adiabatic method need the initial state to be the ground state of \(H(0)\), rather than any easy state?
The adiabatic theorem says the state follows the instantaneous eigenstate it starts in. Starting in a superposition of eigenstates, each component follows its own eigenstate, and the final state is a superposition with the same weights, not the desired final ground state.
How long must the evolution take? Simulate the one-qubit example for several total times \(T\), with \(s=t/T\), and record the probability of ending in \(\ket0\).
For \(T=0.5\) the result is barely better than the sudden \(0.5\); by \(T=4\) it is \(0.986\); by \(T=8\), \(0.9997\). The transition from failure to success happens when \(T\) becomes large compared with the inverse square of the gap.
The adiabatic theorem makes this quantitative. A commonly quoted sufficient condition (see Jansen, Ruskai, and Seiler, 2007, for rigorous versions) is
\[T\gg\frac{\max_s\|dH/ds\|}{g_{\min}^2},\]
where \(g_{\min}\) is the minimum gap along the path. Rigorous versions add terms involving the second derivative and integrals of \(1/g^3\), and require the gap to stay open. The scaling \(1/g_{\min}^2\) is the important part: halving the gap quadruples the time. The schedule matters too. Moving slowly only near the minimum gap and quickly elsewhere, a local adiabatic schedule, can reduce the time to about \(1/g_{\min}\) in favourable cases.
Adiabatic evolution can itself be simulated on a gate-based computer by product formulas, with \(H(s)\) changing from step to step, so it is also a method for preparing ground states within the circuit model. Conversely, adiabatic computation with general Hamiltonians is polynomially equivalent to the circuit model (Aharonov et al., 2004): it is a different presentation of the same computational power, not more.
For the example, \(\|dH/ds\|=\|X-Z\|=\sqrt2\) and \(g_{\min}=\sqrt2\). What does the heuristic condition suggest?
\(T\gg\sqrt2/2\approx0.7\). The simulation shows good fidelity for \(T\gtrsim4\), consistent with “much greater than” meaning a factor of five or so.
A family of problems has \(g_{\min}=2^{-n/2}\). What runtime does the adiabatic condition suggest, and what does this mean?
\(T\sim2^n\) with a linear schedule, or \(2^{n/2}\) with a local schedule. Exponentially small gaps mean exponential runtime: the adiabatic algorithm is then no better than brute force or Grover search.
Why can the adiabatic runtime not be estimated by computing \(g_{\min}\) for a large instance on a classical computer before running?
Computing the gap of a large Hamiltonian is generally as hard as the problem the adiabatic algorithm is meant to solve; for many-body systems it requires exponentially large classical computations. In practice one estimates gaps on small instances, by theory for special families, or not at all, which makes runtime claims for large instances uncertain.
Can slow evolution solve optimization problems? An optimization problem with cost function \(C(x)\) on bit strings can be encoded in a diagonal Hamiltonian \(H_C=\sum_xC(x)\ket x\bra x\), whose ground state is the optimal \(x\). Adiabatic optimization evolves from \(H(0)=-\sum_jX_j\), with ground state \(\ket+^{\otimes n}\), to \(H(1)=H_C\) (Farhi, Goldstone, Gutmann, and Sipser, 2000). Measuring at the end gives a low-cost string if the evolution was slow enough.
What can be proved? For unstructured search, \(H_C\) marks one string among \(N=2^n\), and the minimum gap is of order \(1/\sqrt N\) (exactly \(1/\sqrt N\) for the standard choice of initial and final Hamiltonians as projectors). A linear schedule then takes time of order \(N\), no better than classical search; the local schedule of Roland and Cerf (2002) achieves order \(\sqrt N\), matching Grover. For hard optimization problems, there are instances where the gap is exponentially small, sometimes because of avoided crossings late in the evolution: points where the two lowest eigenvalues come very close, and the ground state switches abruptly between distant configurations. No proof shows that adiabatic optimization solves NP-hard problems efficiently, and the evidence suggests it does not in general.
The honest summary is that slow evolution is a design principle with a clear condition for success, a large enough gap, that is hard to verify for large instances. When an algorithm’s guarantee depends on an unknown gap, the remaining option is to give up guarantees altogether and search for good parameters empirically.
Write \(H_C\) for MaxCut on a single edge between two qubits, with cost \(-1\) for a cut edge and \(0\) otherwise.
The edge is cut when the bits differ, \(Z_1Z_2=-1\). Cost \(-\frac12(1-Z_1Z_2)\), which is \(-1\) on \(\ket{01},\ket{10}\) and \(0\) on \(\ket{00},\ket{11}\). The ground states are the two cut assignments.
For search with \(N=2^{20}\), compare the runtimes of linear and local adiabatic schedules with the gap \(1/\sqrt N\).
\(g_{\min}=2^{-10}\). Linear: about \(1/g_{\min}^2=2^{20}\approx10^6\), of order \(N\). Local: about \(1/g_{\min}=2^{10}=1024\), of order \(\sqrt N\), matching Grover up to constants.
An experiment reports that an adiabatic device found optimal solutions faster than a classical heuristic on a set of instances. What would establish a genuine algorithmic advantage?
Scaling with instance size, not a fixed size; comparison with the best classical algorithms for the same instances, including tuned heuristics; total time-to-solution including repetitions and calibration; instance families not chosen to favour the device; and ideally a theoretical reason, such as a gap bound, for the observed scaling.
Attempt each question before opening a hint. Explain your reasoning, then compare it with the solution.
For H = Z₁Z₂ + X₁, compute the commutator [Z₁Z₂, X₁] and its norm, and give the first-order step count for t = 2, ε = 0.01.
Z₁ and X₁ anticommute; Z₂ is untouched.
r ≥ t²‖[A,B]‖/(2ε).
Z₁Z₂X₁ − X₁Z₁Z₂ = (Z₁X₁ − X₁Z₁)Z₂ = 2iY₁Z₂, with norm 2. So r ≥ 4·2/(2·0.01) = 400 first-order steps by the bound.
Using the rule that Pauli strings commute when they anticommute on an even number of positions, decide which pairs commute among \(X_1X_2\), \(Y_1Y_2\), \(Z_1Z_2\), and \(X_1Z_2\). What is the first-order Trotter error for the Heisenberg pair Hamiltonian \(H=X_1X_2+Y_1Y_2+Z_1Z_2\)?
Compare the strings position by position: distinct non-identity Paulis anticommute.
The first-order bound sums \(\|[H_j,H_k]\|\) over pairs.
\(X_1X_2\), \(Y_1Y_2\), and \(Z_1Z_2\) differ from each other at both positions, so every pair anticommutes twice and commutes. \(X_1Z_2\) anticommutes once with \(X_1X_2\) (position 2) and once with \(Z_1Z_2\) (position 1), so those pairs anticommute; with \(Y_1Y_2\) it differs at both positions and commutes. For the Heisenberg pair all commutators vanish, so \(e^{-iHt}=e^{-iX_1X_2t}e^{-iY_1Y_2t}e^{-iZ_1Z_2t}\) exactly, with no Trotter error for any \(t\). A single step suffices, which is fast-forwarding. On a chain of many sites, neighbouring bonds share a qubit and do not commute, and the error returns.
For \(H=X+Z\) and target error \(\epsilon=10^{-3}\), how many first-order steps does the commutator bound require at \(t=1\) and at \(t=10\)? How does the second-order formula’s step count, of order \(t^{3/2}/\epsilon^{1/2}\), scale between the same two times?
\(r\ge t^2\|[X,Z]\|/(2\epsilon)=t^2/\epsilon\).
Multiplying \(t\) by \(10\) multiplies \(t^{3/2}\) by about \(31.6\).
The first-order bound needs \(r\ge t^2/\epsilon\): \(1000\) steps at \(t=1\) and \(100\,000\) at \(t=10\), a factor of \(100\) for ten times longer evolution. The second-order count grows by \(10^{3/2}\approx32\). Each first-order step has fixed cost, so the first-order gate count grows as \(t^2\). The no-fast-forwarding theorem requires only linear growth, so the first-order formula pays an extra factor of \(t\). Higher-order formulas, or polynomial methods, close most of that gap.
Write a circuit for e^{−iθ Y₁ Y₂} and count its CNOTs and single-qubit gates.
Map Y to Z with S† followed by H on each qubit.
Then use CNOT, R_z(2θ), CNOT.
On each qubit apply S† then H (mapping Y to Z); CNOT from qubit 1 to qubit 2; R_z(2θ) on qubit 2; CNOT from qubit 1 to 2; then H then S on each qubit. Two CNOTs, one rotation, and eight Clifford single-qubit gates (some of which may cancel with neighbouring terms in a full Trotter step).
For \(H=X+Z\), a single first-order step of length \(\delta=0.1\) has actual error about \(0.00998\), while ten such steps (total time \(1\)) have actual error about \(0.070\). (a) Compare each with the bounds \(\frac{\delta^2}{2}\|[X,Z]\|\) and \(\frac{t^2}{2r}\|[X,Z]\|\). (b) Why is the one-step bound nearly tight while the ten-step bound is loose? (c) Using the actual \(1/r\) scaling, estimate the steps needed for error \(0.01\).
\(\|[X,Z]\|=2\).
The triangle inequality adds error norms; it ignores directions.
(a) The one-step bound is \(0.01\cdot2/2=0.01\), matching \(0.00998\) almost exactly. The ten-step bound is \(1\cdot2/20=0.1\), against an actual \(0.070\). (b) For a single short step the leading error term \(-\frac{\delta^2}{2}[X,Z]\) dominates and the bound captures it. Over many steps the bound adds error norms as if every step erred in the worst direction, but the step errors are rotated by the evolution in between and partly cancel, so the total grows more slowly than the sum. (c) The actual error is about \(0.70/r\), so \(r\approx70\) steps reach \(0.01\), against \(100\) from the bound. That difference translates directly into gate count.
Adiabatic search over \(N\) items uses \(H(s)=I-(1-s)\ket u\bra u-s\ket m\bra m\), where \(\ket u\) is the uniform superposition and \(\ket m\) the marked item. For \(N=4\), compute the gap at \(s=0\), \(s=\frac14\), and \(s=\frac12\) using \(g(s)=\sqrt{1-4(1-1/N)s(1-s)}\). Compare the runtimes suggested by the heuristic \(1/g_{\min}^2\) for \(N=4\) and \(N=2^{20}\).
At \(s=\frac12\), \(g=1/\sqrt N\).
For \(N=4\), \(1-1/N=\frac34\).
For \(N=4\): \(g(0)=1\), \(g(\tfrac14)=\sqrt{1-3\cdot\tfrac{3}{16}}=\sqrt{7/16}\approx0.661\), and \(g(\tfrac12)=\sqrt{1-\tfrac34}=\tfrac12=1/\sqrt N\). The heuristic \(1/g_{\min}^2=N\) suggests time of order \(4\) here and \(2^{20}\approx10^6\) for \(N=2^{20}\), no better than classical search. The gap is small only in a window near \(s=\frac12\), of width about \(1/\sqrt N\). A schedule that slows down only there, the local adiabatic schedule, needs time of order \(\sqrt N\approx10^3\), matching Grover search. The improvement comes from knowing where the gap is small.
A first-order simulation has \(L=50\) Pauli terms and \(r=1000\) steps, so \(5\times10^4\) arbitrary rotations. The total error budget is \(10^{-3}\), half allotted to rotation synthesis. Using about \(3\log_2(1/\epsilon_{\mathrm{rot}})\) \(T\) gates per rotation, find \(\epsilon_{\mathrm{rot}}\), the \(T\) count per rotation, and the total \(T\) count. How would halving the number of rotations change the answer?
Synthesis errors add: \(5\times10^4\,\epsilon_{\mathrm{rot}}\le5\times10^{-4}\).
\(\log_2(10^8)\approx26.6\).
The rotation errors add, so \(\epsilon_{\mathrm{rot}}\le5\times10^{-4}/(5\times10^4)=10^{-8}\). Each rotation then costs about \(3\times26.6\approx80\) \(T\) gates, and the total is about \(5\times10^4\times80=4\times10^6\) \(T\) gates. Halving the number of rotations halves the count and also relaxes each rotation’s budget to \(2\times10^{-8}\), saving about \(3\) \(T\) gates per rotation: roughly \(2\times10^6\) in total. The logarithmic dependence means precision is cheap per gate, but the count of rotations, set by \(L\) and \(r\), drives the cost.
A student rescales a Hamiltonian, \(H'=H/10\), and argues that \(e^{-iH't'}\) with \(t'=10t\) equals \(e^{-iHt}\) while every commutator becomes \(100\) times smaller, so Trotterization gets much cheaper. Check the first-order step count and the qubitization cost before and after rescaling.
\([A/10,B/10]=[A,B]/100\), but the time is \(10t\).
Qubitization costs about \(\lambda t\).
For first order, \(r\ge t'^2\|[A',B']\|/(2\epsilon)=(10t)^2\cdot\|[A,B]\|/100/(2\epsilon)=t^2\|[A,B]\|/(2\epsilon)\): unchanged. For qubitization, \(\lambda'=\lambda/10\) and \(\lambda't'=\lambda t\): unchanged. The evolution operator depends only on the product \(Ht\), and so does every honest cost formula. Rescaling changes units, not difficulty. A claimed speedup from “simulating for less time” must be checked against \(\|H\|t\) or \(\lambda t\), the dimensionless quantities the costs depend on.
An adiabatic algorithm for a family of instances has minimum gap g(n) = 1/n² for n qubits. Estimate the runtime with a linear schedule and with a local schedule, and explain whether this would give a polynomial-time algorithm. What would you need to know to trust the gap estimate?
Linear schedule: T ~ 1/g².
Local schedule: T ~ 1/g in favourable cases.
Linear: T ~ n⁴; local: T ~ n². Both are polynomial, so if the gap really scales as 1/n² for the whole family, the adiabatic algorithm runs in polynomial time (and can be simulated by a polynomial-size circuit). To trust it, one needs a proof of the gap bound for all instances, or at least strong evidence across sizes, because numerical gaps on small n can miss exponentially small avoided crossings that appear only at larger n.
Derive the gap formula used in adiabatic search. (a) Show that \(H(s)=I-(1-s)\ket u\bra u-s\ket m\bra m\) acts as \(I\) on states orthogonal to both \(\ket u\) and \(\ket m\), so the gap is set by the two-dimensional block they span. (b) In that block, find the trace and determinant of \(M=(1-s)\ket u\bra u+s\ket m\bra m\) using \(|\braket um|^2=1/N\). (c) Deduce \(g(s)=\sqrt{1-4(1-1/N)s(1-s)}\) and its minimum.
For a \(2\times2\) Hermitian matrix the eigenvalue difference is \(\sqrt{(\mathrm{tr})^2-4\det}\).
\(\det M=s(1-s)(1-|\braket um|^2)\).
(a) A state orthogonal to \(\ket u\) and \(\ket m\) is annihilated by both projectors, so \(H(s)\) gives eigenvalue \(1\). The other eigenvalues live in the span of \(\ket u,\ket m\), where \(H=I-M\). (b) \(\mathrm{tr}M=(1-s)+s=1\). For the determinant, write \(M\) in an orthonormal basis \(\{\ket m,\ket{m^\perp}\}\) with \(\ket u=\frac{1}{\sqrt N}\ket m+\sqrt{1-\frac1N}\ket{m^\perp}\); a direct calculation gives \(\det M=s(1-s)(1-\frac1N)\). (c) The eigenvalues of \(M\) differ by \(\sqrt{1-4s(1-s)(1-\frac1N)}\), and so do those of \(H=I-M\). The smaller eigenvalue of \(H\), \(\frac12(1-\sqrt{\cdots})\), is the ground energy. The larger, \(\frac12(1+\sqrt{\cdots})\le1\), is the first excited energy, below the eigenvalue \(1\) of the orthogonal states. So the gap is \(g(s)=\sqrt{1-4(1-1/N)s(1-s)}\), minimized at \(s=\frac12\) with \(g_{\min}=1/\sqrt N\).
Where this leads. the Variational Methods chapter drops guaranteed constructions: it optimizes the parameters of short circuits against measured objectives, and asks what such heuristics can and cannot establish.