Variational Methods

What is the cheapest family of states whose output we can tune? Take a parameterized circuit \(U(\theta)\): a fixed gate sequence in which some rotation angles \(\theta=(\theta_1,\ldots,\theta_m)\) are adjustable. Applied to a simple input, it produces a family of trial states \(\ket{\psi(\theta)}=U(\theta)\ket{0^n}\). The choice of circuit structure is called an ansatz, a German word for an educated guess about the form of the answer.

Here is the running example. Apply \(R_y(\theta)=e^{-i\theta Y/2}\) to the first of two qubits in \(\ket{00}\), then a CNOT:

\[\ket{\psi(\theta)}=\cos\tfrac\theta2\,\ket{00}+\sin\tfrac\theta2\,\ket{11}.\]

At \(\theta=\pi/2\) this is the Bell state \(\frac{1}{\sqrt2}(\ket{00}+\ket{11})\). But the family contains only states with real, nonnegative-ratio amplitudes on \(\ket{00}\) and \(\ket{11}\). It cannot produce \(\ket{01}\), or \(\frac{1}{\sqrt2}(\ket{00}-\ket{11})\), or any state outside a one-dimensional curve in the space of two-qubit states.

The expressiveness of an ansatz is the set of states it can reach. A more expressive ansatz, with more gates and parameters, is more likely to contain a good answer, but it costs more to run, has more parameters to optimize, and, as we will find, can become harder to train. The ansatz is an assumption about where good answers lie. If that assumption is wrong, no optimizer can compensate. Physically motivated ansätze, built from the structure of the problem’s Hamiltonian, encode this assumption deliberately; hardware-efficient ansätze, built from whatever gates a device does well, make it implicitly.

How do we score a trial state? A variational algorithm minimizes a cost function, usually the expectation value of a Hamiltonian in the trial state:

\[E(\theta)=\bra{\psi(\theta)}H\ket{\psi(\theta)}.\]

The quantum computer prepares the state and measures; the classical computer averages. For \(H=\sum_jh_jP_j\), each Pauli string is measured separately: rotate each qubit so its Pauli becomes \(Z\) (a measurement basis change: a Hadamard gate for \(X\), \(S^\dagger\) then a Hadamard gate for \(Y\)), measure in the computational basis, and multiply the outcomes \(\pm1\). The average of that product estimates the expectation value \(\langle P_j\rangle\).

Take \(H=-Z\otimes Z-X\otimes X\), whose ground state is the Bell state with energy \(-2\). In the trial state, \(Z\otimes Z\) gives \(+1\) on both \(\ket{00}\) and \(\ket{11}\), so \(\langle ZZ\rangle=1\). And \(X\otimes X\) swaps \(\ket{00}\leftrightarrow\ket{11}\), so \(\langle XX\rangle=2\cos\frac\theta2\sin\frac\theta2=\sin\theta\). Therefore

\[E(\theta)=-1-\sin\theta,\]

minimized at \(\theta=\pi/2\), where \(E=-2\): the trial family contains the exact ground state.

Pauli strings that commute qubit by qubit, such as \(Z_1Z_2\) and \(Z_1\), can be measured together from the same shots. Grouping terms this way reduces the number of distinct circuits. For molecules with thousands of terms, grouping and other measurement strategies change the cost by orders of magnitude.

How many measurements does one energy estimate need? Each measurement of a Pauli string gives \(\pm1\), so its variance is \(1-\langle P\rangle^2\le1\). Averaging \(R\) independent shots gives an estimate with variance at most \(1/R\). This statistical uncertainty is shot noise. For a sum \(H=\sum_jh_jP_j\) estimated term by term with \(R_j\) shots for term \(j\),

\[\mathrm{Var}(\hat E)=\sum_j\frac{h_j^2\,\sigma_j^2}{R_j},\qquad\sigma_j^2=1-\langle P_j\rangle^2.\]

How should a fixed total \(R=\sum_jR_j\) be divided? Write \(a_j=|h_j|\sigma_j\). The Cauchy–Schwarz inequality gives

\[\Big(\sum_ja_j\Big)^2=\Big(\sum_j\frac{a_j}{\sqrt{R_j}}\sqrt{R_j}\Big)^2\le\Big(\sum_j\frac{a_j^2}{R_j}\Big)\Big(\sum_jR_j\Big),\]

so the variance is at least \((\sum_ja_j)^2/R\), with equality when \(a_j/\sqrt{R_j}\) is proportional to \(\sqrt{R_j}\), that is, \(R_j\propto|h_j|\sigma_j\). Then

\[\mathrm{Var}(\hat E)=\frac{\big(\sum_j|h_j|\sigma_j\big)^2}{R}\le\frac{\lambda^2}{R},\]

with \(\lambda=\sum_j|h_j|\), the same 1-norm that set the cost of a linear-combination block encoding. To reach precision \(\epsilon\) in the energy takes about \(\lambda^2/\epsilon^2\) shots. This shot allocation rule is the measurement analogue of choosing PREPARE amplitudes proportional to \(\sqrt{|h_j|}\).

The \(1/\epsilon^2\) scaling is the price of measuring instead of using phase estimation, which reaches \(1/\epsilon\) with long coherent circuits. For molecules, energies are measured in the hartree, the natural energy unit of molecular calculations. There \(\lambda\) is typically tens to hundreds of hartree, while the target precision is about \(10^{-3}\) hartree, so \(\lambda/\epsilon\) is about \(10^4\) to \(5\times10^5\), and \(\lambda^2/\epsilon^2\) ranges from about \(10^{8}\) to \(2.5\times10^{11}\) shots per energy evaluation, and an optimization needs many evaluations. Measurement cost is often the bottleneck of variational methods.

Gradient descent updates \(\theta\leftarrow\theta-\eta\,\partial E/\partial\theta\) with a step size \(\eta\). How do we get the derivative from measurements? Finite differences, \((E(\theta+h)-E(\theta-h))/2h\), are biased for finite \(h\) and amplify shot noise as \(h\to0\). There is something better.

Suppose \(\theta\) enters the circuit once, through a gate \(e^{-i\theta P/2}\) with \(P^2=I\) (any Pauli rotation). Then \(E(\theta)\) has the form \(a+b\cos\theta+c\sin\theta\): the gate contributes \(\cos\frac\theta2\) and \(\sin\frac\theta2\) to the amplitudes, and the expectation value is quadratic in those, hence a trigonometric polynomial of frequency one. For such a function, the parameter-shift rule holds exactly:

\[\frac{dE}{d\theta}=\frac{E(\theta+\pi/2)-E(\theta-\pi/2)}{2}.\]

Check on our example: \[\frac12\big(-1-\sin(\theta+\frac\pi2)+1+\sin(\theta-\frac\pi2)\big)=\frac12(-\cos\theta-\cos\theta)=-\cos\theta\], which is the derivative of \(-1-\sin\theta\). The shifts are large, \(\pm\pi/2\), so the two evaluations differ by a quantity of order one and shot noise is not amplified.

Each parameter needs two circuit evaluations per gradient component, so a gradient of \(m\) parameters costs \(2m\) energy estimates, each with its own shot budget. If one parameter controls several gates, apply the rule to each occurrence separately and add the results (the chain rule); applying a single shift to the shared parameter would be wrong.

A gradient estimated from finitely many shots is a random variable. Gradient descent with noisy gradients is stochastic optimization: on average it moves downhill, but each step jitters. Near a minimum, the true gradient becomes small and the noise dominates, so the optimizer wanders. A sensible stopping rule halts when the estimated improvement is smaller than its uncertainty, rather than chasing digits the measurements cannot resolve. Increasing shots as the optimization proceeds is a common compromise.

Hardware noise is different in kind. A device interacts with its surroundings, which can acquire which-path records of the computation, just as garbage did; this unwanted recording is called decoherence. It turns the prepared state into a mixture, and together with gate errors it biases the measured energy toward the energy of that mixture. More shots do not remove a bias. They only estimate the biased value more precisely. Error mitigation techniques, such as extrapolating results from deliberately amplified noise back to zero noise, reduce the bias at a cost in shots that grows rapidly with circuit size. They do not provide the guarantees of error correction.

Short circuits are the reason variational methods were proposed for noisy devices: fewer gates means less accumulated error. But short circuits limit expressiveness, and the measurement overhead of energy estimation, multiplied by the number of optimizer iterations and gradient components, can be enormous. These pressures pull in opposite directions, and their balance is problem dependent.

Are the gradients large enough to follow? For many ansätze, the landscape \(E(\theta)\) becomes extremely flat as the number of qubits grows. McClean, Boixo, Smelyanskiy, Babbush, and Neven (2018) showed that for sufficiently deep random circuits, the variance of any gradient component over random parameter choices decays exponentially in \(n\). At a random starting point, the gradient is then exponentially small with high probability, and estimating it to useful accuracy requires exponentially many shots. Such a landscape is a barren plateau.

The intuition is concentration of measure. In a space of dimension \(2^n\), a random state has nearly the same expectation value for any fixed observable, close to its average over all states. Circuits deep enough to scramble their input behave like random states, so the energy barely depends on the parameters.

Barren plateaus also arise from global cost functions (measuring all qubits at once), from too much entanglement with unmeasured qubits, and from hardware noise, which flattens landscapes even for shallow circuits. Remedies include problem-inspired ansätze, local cost functions, careful initialization near a known good state, and layer-by-layer training. Trainability, whether the optimum can be found with reasonable resources, is a separate property from expressiveness, whether the optimum is in the family at all. A recent line of work (Cerezo et al., 2023) also suggests that ansätze provably free of barren plateaus often have landscapes that can be estimated efficiently classically, which raises the question of whether they offer any quantum advantage.

What does minimizing the energy guarantee? The variational principle answers this. For any normalized state \(\ket\psi=\sum_kc_k\ket{E_k}\) expanded in eigenstates of \(H\),

\[\bra\psi H\ket\psi=\sum_k|c_k|^2E_k\ge E_0,\]

with equality exactly when \(\ket\psi\) is a ground state. So minimizing the energy over a family of trial states gives an upper bound on the ground energy, and the minimizer is the family’s best approximation in energy.

The variational quantum eigensolver (Peruzzo et al., 2014) combines a parameterized circuit, term-by-term energy measurement, and a classical optimizer. In our example the family contains the exact ground state and the optimizer finds it. In general three gaps separate the answer from the truth: the expressiveness gap (the best state in the family has energy above \(E_0\)), the optimization gap (the optimizer finds a local minimum or stops early), and the estimation gap (shot noise and hardware bias in the reported energy).

The variational principle bounds the energy, not the state. How close must the state be? Let \(\Delta\) be the energy gap. Every excited energy is at least \(E_0+\Delta\), so

\[\bra\psi H\ket\psi=\sum_k|c_k|^2E_k\ge|c_0|^2E_0+(1-|c_0|^2)(E_0+\Delta)=E_0+(1-|c_0|^2)\Delta.\]

An energy error \(\delta\) therefore implies a squared overlap with the ground state of at least \(|c_0|^2\ge1-\delta/\Delta\). A trial state with energy close to \(E_0\) can still be far from the ground state if the gap is small. Properties other than energy computed from a variational state inherit this uncertainty.

Can a short alternating circuit approximate slow evolution for optimization? For a combinatorial problem with cost \(C(x)\) on bit strings, the quantum approximate optimization algorithm (QAOA, Farhi, Goldstone, and Gutmann, 2014) prepares

\[\ket{\gamma,\beta}=e^{-i\beta_pB}e^{-i\gamma_pC}\cdots e^{-i\beta_1B}e^{-i\gamma_1C}\ket+^{\otimes n},\qquad B=\sum_jX_j,\]

with \(2p\) parameters, and measures to sample bit strings. The cost layer \(e^{-i\gamma C}\) is diagonal: it only changes phases, so by itself it cannot change any measurement probability. The mixer \(e^{-i\beta B}\) does not commute with \(C\), and it converts phase differences into probability differences, the same partnership of phase and interference that drove every algorithm since Deutsch’s circuit. As \(p\to\infty\), suitable parameters approximate adiabatic evolution by a product formula, so QAOA can do at least as well as the adiabatic algorithm in that limit.

Take MaxCut on one edge: \(C=\frac12(I-Z_1Z_2)\), which is \(1\) when the bits differ. What does one layer achieve? The cost layer gives phase \(1\) to \(\ket{00},\ket{11}\) and \(w=e^{-i\gamma}\) to \(\ket{01},\ket{10}\):

\[e^{-i\gamma C}\ket{++}=\tfrac12\big(\ket{00}+\ket{11}\big)+\tfrac w2\big(\ket{01}+\ket{10}\big).\]

The expected cut is \(F=\tfrac12-\tfrac12\langle Z_1Z_2\rangle\) after the mixer. Rather than evolve the state, move the mixer onto the observable. Since \(ZX=-XZ\), \(e^{i\beta X}Ze^{-i\beta X}=Ze^{-2i\beta X}=\cos2\beta\,Z+\sin2\beta\,Y\), using \(ZX=iY\). So after the mixer \(Z_1Z_2\) becomes \((\cos2\beta\,Z_1+\sin2\beta\,Y_1)(\cos2\beta\,Z_2+\sin2\beta\,Y_2)\), evaluated in the state above. There \(\langle Z_1Z_2\rangle=0\), since half the weight has each sign, and \(\langle Y_1Y_2\rangle=0\) likewise. The cross terms do not vanish: applying \(Z_1Y_2\) and taking the overlap gives \(\langle Z_1Y_2\rangle=\tfrac i2(w^*-w)=-\sin\gamma\), and \(\langle Y_1Z_2\rangle=-\sin\gamma\) by symmetry. Hence \(\langle Z_1Z_2\rangle=-2\sin2\beta\cos2\beta\sin\gamma=-\sin4\beta\sin\gamma\), and the expected cut is

\[F(\gamma,\beta)=\tfrac12+\tfrac12\sin(4\beta)\sin\gamma,\]

which equals \(1\) at \(\gamma=\pi/2\), \(\beta=\pi/8\): one layer finds the cut with certainty.

For larger graphs, what is proved is modest. On 3-regular graphs, depth-1 QAOA guarantees a cut at least \(0.6924\) times the maximum (Farhi, Goldstone, and Gutmann, 2014), below the classical guarantee of about \(0.878\) of Goemans and Williamson (1995). At higher depth, performance improves, but proven guarantees are few, and the parameters must be optimized for each instance or transferred from similar ones. Optimizing the expected cost and producing a single high-quality string are related but different goals: the distribution of samples matters, and the best of many samples is what a user keeps.

If there is no proof, what would convince us? Variational algorithms are heuristics: procedures that may work well in practice without a proof that they do. Heuristics are valuable. Most of classical optimization runs on them. But claims about them must rest on evidence, and the evidence must be of the right kind.

Good benchmarking compares against the best classical methods for the same instances, including tuned heuristics and specialized solvers, not a convenient baseline. It reports how performance scales with problem size, since a fixed-size advantage can be overtaken. It counts total resources: shots, optimizer iterations, calibration, and classical preprocessing, and reports time-to-solution or cost at matched solution quality. It uses instance families not chosen to favour the method. And it states what was simulated, what was run on hardware, and what was extrapolated.

Two further cautions. Shallow circuits on restricted geometries can sometimes be simulated classically, for example with tensor networks, which store weakly entangled states compactly, so small instances showing good quantum results do not show a quantum advantage. And the absence of proven guarantees means that negative results matter as much as positive ones. Variational methods have taught much about devices and problems; whether they offer an asymptotic advantage for any practical task is an open question.

Throughout, measurement has served as the last step, or as feedback to a classical optimizer between circuits. Could an outcome obtained in the middle of a circuit steer the rest of it?