Variational Quantum Eigensolver

The Variational Quantum Eigensolver is the algorithm that made quantum chemistry a near-term ambition rather than a fault-tolerant one. Its target is the oldest hard problem in quantum simulation: find the lowest eigenvalue E0E_0 of a Hermitian operator HH — the electronic Hamiltonian of a molecule, or the Hamiltonian of a lattice magnet. Exact diagonalization is hopeless because the Hilbert-space dimension is 2n2^n in the qubit number nn, and the textbook quantum answer, quantum phase estimation, needs deep coherent circuits that implement controlled time evolution e−iHte^{-iHt} — far past the reach of today's noisy, intermediate-scale (NISQ) hardware. VQE sidesteps both walls by turning eigenvalue estimation into optimization. It leans on a theorem older than quantum computing itself, the Rayleigh–Ritz variational principle: for any normalized trial state ∣ψ⟩|\psi\rangle, the expectation ⟨ψ∣H∣ψ⟩\langle\psi|H|\psi\rangle can never dip below E0E_0. So one parameterizes a family of states ∣ψ(θ)⟩=U(θ)∣0⟩⊗n|\psi(\boldsymbol{\theta})\rangle=U(\boldsymbol{\theta})|0\rangle^{\otimes n} — the ansatz — and slides the parameters downhill until the energy bottoms out. The quantum device does only what shallow circuits do well: prepare an entangled state and measure local observables; every decision about where to step next is offloaded to a classical optimizer. To make the energy measurable, HH is expanded as a weighted sum of Pauli strings H=∑kckPkH=\sum_k c_k P_k, each read out by rotating into its eigenbasis and sampling. VQE thus lives at the crossing of quantum algorithms, classical optimization, and quantum chemistry, and it is the prototype of the whole family of variational quantum algorithms, QAOA included.

Key takeaways
  • VQE estimates the ground-state energy of a Hamiltonian HH by minimizing E(θ)=⟨ψ(θ)∣H∣ψ(θ)⟩E(\boldsymbol{\theta})=\langle\psi(\boldsymbol{\theta})|H|\psi(\boldsymbol{\theta})\rangle over a parameterized circuit U(θ)U(\boldsymbol{\theta}), using the variational bound E(θ)≥E0E(\boldsymbol{\theta})\ge E_0.
  • It is a hybrid quantum–classical loop: a quantum processor prepares ∣ψ(θ)⟩|\psi(\boldsymbol{\theta})\rangle and measures Pauli expectation values; a classical optimizer proposes the next θ\boldsymbol{\theta}.
  • The Hamiltonian is written as a weighted sum of Pauli strings H=∑kckPkH=\sum_k c_k P_k, and the energy reassembles by linearity as E=∑kck⟨Pk⟩E=\sum_k c_k\langle P_k\rangle; the measurement cost scales as (∑k∣ck∣)2/ϵ2(\sum_k|c_k|)^2/\epsilon^2.
  • Fermionic (electronic-structure) problems are mapped to qubits by the Jordan–Wigner or Bravyi–Kitaev encoding; common ansätze are UCCSD (chemistry-tailored) and hardware-efficient (shallow but landscape-prone).
  • Gradients come exactly and cheaply from the parameter-shift rule, but expressive ansätze suffer barren plateaus where the gradient variance decays as 2−n2^{-n}.
  • Key idea: the quantum computer evaluates the energy E(θ)=⟨ψ(θ)∣H∣ψ(θ)⟩E(\boldsymbol{\theta})=\langle\psi(\boldsymbol{\theta})|H|\psi(\boldsymbol{\theta})\rangle term-by-term over Pauli strings while a classical optimizer descends the surface; because E(θ)≥E0E(\boldsymbol{\theta})\ge E_0 always, the minimum is a rigorous upper bound on the ground-state energy.
The core insight

Replace deep coherent eigenvalue extraction with a hybrid loop: the quantum computer evaluates E(θ)=⟨ψ(θ)∣H∣ψ(θ)⟩E(\boldsymbol{\theta})=\langle\psi(\boldsymbol{\theta})|H|\psi(\boldsymbol{\theta})\rangle term-by-term over Pauli strings, and a classical optimizer descends this surface. The variational principle guarantees E(θ)≥E0E(\boldsymbol{\theta})\ge E_0 always, so the minimum is a rigorous upper bound on the ground-state energy.

Where this lives today

VQE is the workhorse of the NISQ era, and its story since 2014 is a candid mix of milestones and hard lessons. The first demonstration solved HeH+\mathrm{HeH}^+ on a two-photon chip; Google and IBM soon traced full dissociation curves of H2\mathrm{H_2}, LiH\mathrm{LiH}, and BeH2\mathrm{BeH_2} on superconducting qubits, and in 2020 Google's Sycamore ran a 12-qubit Hartree–Fock VQE on hydrogen chains and diazene isomerization — the largest chemistry simulation on a quantum computer at the time. VQE is also the canonical proving ground for error mitigation: zero-noise extrapolation, probabilistic error cancellation, and symmetry verification were all sharpened on variational chemistry, because the variational bound gives a clean target to restore. Beyond molecules, VQE attacks Hubbard and Heisenberg lattice models relevant to high-temperature superconductivity, and excited-state variants reach spectral gaps and reaction barriers. The open question is whether any of this beats a classical computer. Barren plateaus (2018) showed that expressive ansätze train exponentially poorly, and a growing body of "de-quantization" results shows many barren-plateau-free circuits are also classically simulable, sharpening the question of where a real advantage could hide; the measurement wall of equation (8) compounds it. The live frontiers — problem-tailored and symmetry-preserving ansätze, adaptive circuit growth (ADAPT-VQE), classical-shadow measurement, and rigorously-guaranteed error mitigation — all aim to keep VQE's shallow-circuit promise while escaping its two exponential traps. Whether it delivers practical quantum advantage or serves as the instructive prototype that taught the field what a variational algorithm must avoid, VQE remains the entry point to the whole variational program, QAOA included.

The physical picture

The whole algorithm is a division of labour along the line each machine is good at. A quantum processor can cheaply hold and manipulate an entangled state on 2n2^n amplitudes but cannot think — only prepare and measure. A classical computer can optimize but cannot store the state. VQE couples them: the quantum side is a black-box energy oracle that, handed a parameter vector θ\boldsymbol{\theta}, prepares ∣ψ(θ)⟩|\psi(\boldsymbol{\theta})\rangle, measures every Pauli term of HH, and reports a number E(θ)E(\boldsymbol{\theta}); the classical side takes that number (and its gradient, via the parameter-shift rule) and decides where to step next. The first figure is this loop.

Figure
The VQE hybrid loop. The quantum processor prepares the ansatz state ψ(θ) and measures each Pauli term ⟨P_k⟩; the classical optimizer assembles the energy E(θ) = Σ c_k⟨P_k⟩ and updates θ. For H₂ in the STO-3G basis this is 4 qubits and about 15 Pauli strings, with term count growing as O(n⁴).
One loop, two machines: the quantum side measures energy, the classical side steers QUANTUM PROCESSOR 1 · Prepare the ansatz state |ψ(θ)⟩ = U(θ) |0⟩⊗ⁿ 2 · Measure each Pauli term ⟨P_k⟩ , sampled in its basis shallow circuit · repeated N times CLASSICAL OPTIMIZER 3 · Assemble the energy E(θ) = Σ_k c_k ⟨P_k⟩ 4 · Update the parameters θ ← θ − η ∇E(θ) gradient from parameter shift ⟨P_k⟩ new θ loop until E(θ) converges to the variational minimum ≥ E₀ H₂ / STO-3G: n = 4 qubits · ~15 Pauli strings · cost ~ O(n⁴)

Two features deserve a second look. First, the quantum computer never holds coherence for long — it prepares a shallow state, measures, and resets, so decoherence has little time to accumulate; that is precisely why VQE is a NISQ algorithm and phase estimation is not. Second, the energy is not read off in one shot but as a weighted sum of Pauli expectations, each estimated by many repetitions. The variational principle is what makes the scheme worth the trouble: whatever state the circuit prepares, its energy is an upper bound that can only improve as the optimizer works. The second figure shows that floor for the simplest case — a single qubit, where the whole landscape is one sinusoid and the bound is a horizontal line the curve just kisses.

Figure
The variational landscape for a single qubit, H = Z + ½X, with the R_y(θ) ansatz so E(θ) = cos θ + ½ sin θ. The optimizer slides downhill to the minimum E(θ) = −√1.25 = −1.11803 at θ = 3.605 rad, which exactly equals the true ground-state energy E₀ = −√1.25 — the variational floor is touched, not crossed.
The variational floor: E(θ) never drops below E₀ = −1.11803, and the minimum touches it +1.0 0 −1.0 0 π 2π parameter θ (rad) energy E(θ) variational floor E₀ = −√1.25 = −1.11803 min at θ* = 3.605 rad E(θ) = cos θ + ½ sin θ

Two warnings before the equations. The landscape is only this friendly for one qubit; in general E(θ)E(\boldsymbol{\theta}) is a rugged, high-dimensional, non-convex surface whose most dangerous pathology is the barren plateau — for expressive random circuits the gradient shrinks like 2−n2^{-n}, so the surface goes exponentially flat and the optimizer is blind. And measuring HH is really measuring a budget of Pauli strings; because commuting Paulis can be read out together, grouping them into compatible sets is the algorithm's version of caching and can cut the runtime by an order of magnitude.

The core equations

Everything begins with the variational cost function — the Rayleigh quotient of HH in the trial state. The inequality on the right is the entire justification for the method: the measured energy is a rigorous upper bound on E0E_0, so lower is always better and never wrong.

E(θ)=⟨ψ(θ)∣H∣ψ(θ)⟩⟨ψ(θ)∣ψ(θ)⟩≥E0(1)E(\boldsymbol{\theta})=\frac{\langle\psi(\boldsymbol{\theta})|H|\psi(\boldsymbol{\theta})\rangle}{\langle\psi(\boldsymbol{\theta})|\psi(\boldsymbol{\theta})\rangle}\ge E_0 \qquad (1)

where: E(θ)E(\boldsymbol{\theta}) — variational energy, the objective minimized over the parameters; ∣ψ(θ)⟩|\psi(\boldsymbol{\theta})\rangle — trial (ansatz) state produced by the circuit; HH — Hermitian Hamiltonian whose ground energy is sought; E0E_0 — true ground-state energy, the smallest eigenvalue of HH; θ\boldsymbol{\theta} — vector of variational parameters (gate angles). Any state you can build sits at or above the true floor, with equality only at a ground eigenstate. A quantum computer cannot measure a generic HH directly — its detectors read qubits in the computational (ZZ) basis. So HH is expanded in the one operator basis the hardware can handle: tensor products of single-qubit Paulis.

H=∑kck Pk,Pk∈{I,X,Y,Z}⊗n(2)H=\sum_{k} c_k\, P_k,\qquad P_k\in\{I,X,Y,Z\}^{\otimes n} \qquad (2)

where: ckc_k — real coefficient (Hermiticity forces ckc_k real for Hermitian PkP_k); PkP_k — an nn-qubit Pauli string, a tensor product of single-qubit Paulis; nn — number of qubits. There are at most 4n4^n such strings, but physical Hamiltonians are sparse in this basis — molecular Hamiltonians carry only O(n4)\mathcal{O}(n^4) terms. Linearity of the expectation value then collapses the whole energy into a weighted sum of Pauli expectations, each a number between −1-1 and +1+1 that the hardware estimates by rotating into the string's eigenbasis and sampling.

E(θ)=∑kck ⟨ψ(θ)∣Pk∣ψ(θ)⟩(3)E(\boldsymbol{\theta})=\sum_k c_k\,\langle\psi(\boldsymbol{\theta})|P_k|\psi(\boldsymbol{\theta})\rangle \qquad (3)

where: ⟨Pk⟩\langle P_k\rangle — expectation of the Pauli string PkP_k, with eigenvalues ±1\pm1, estimated by sampling in its eigenbasis. The energy is reassembled classically from these pieces — no single circuit ever measures HH as a whole. The trial state is produced by a layered parameterized circuit acting on a fixed reference; the choice of layers — the ansatz — is where all the art (and all the danger) lives.

∣ψ(θ)⟩=UL(θL)⋯U1(θ1) ∣0⟩⊗n(4)|\psi(\boldsymbol{\theta})\rangle=U_L(\theta_L)\cdots U_1(\theta_1)\,|0\rangle^{\otimes n} \qquad (4)

where: Uℓ(θℓ)U_\ell(\theta_\ell) — parameterized layer (single-qubit rotations interleaved with fixed entanglers); LL — circuit depth; ∣0⟩⊗n|0\rangle^{\otimes n} — reference state, in chemistry the Hartree–Fock determinant. A hardware-efficient ansatz picks UℓU_\ell for shallow depth; a chemistry ansatz picks it to respect physics. The flagship chemistry ansatz is unitary coupled cluster with single and double excitations (UCCSD): it dresses the Hartree–Fock reference with a unitary built from fermionic excitation operators.

∣ψUCC⟩=e T−T† ∣ψHF⟩,T=∑iatia aa†ai+∑ijabtijab aa†ab†ajai(5)|\psi_{\mathrm{UCC}}\rangle=e^{\,T-T^{\dagger}}\,|\psi_{\mathrm{HF}}\rangle,\quad T=\sum_{ia}t_i^a\,a_a^{\dagger}a_i+\sum_{ijab}t_{ij}^{ab}\,a_a^{\dagger}a_b^{\dagger}a_j a_i \qquad (5)

where: TT — cluster operator carrying single (tiat_i^a) and double (tijabt_{ij}^{ab}) excitation amplitudes; a†,aa^{\dagger},a — fermionic creation/annihilation operators; tia,tijabt_i^a,t_{ij}^{ab} — the variational amplitudes (components of θ\boldsymbol{\theta}); ∣ψHF⟩|\psi_{\mathrm{HF}}\rangle — Hartree–Fock reference. The anti-Hermitian generator T−T†T-T^\dagger makes eT−T†e^{T-T^\dagger} unitary and hence circuit-implementable. Fermionic operators do not act on qubits until transformed. The Jordan–Wigner map is the standard bridge: it turns each fermionic mode into a qubit and carries fermionic antisymmetry into a Pauli-ZZ parity string on the preceding qubits.

aj†  ↦  12(Xj−iYj)⨂k<jZk(6)a_j^{\dagger}\;\mapsto\;\tfrac{1}{2}\bigl(X_j-iY_j\bigr)\bigotimes_{k<j}Z_k \qquad (6)

where: aj†a_j^\dagger — fermionic creation operator on mode jj; ⨂k<jZk\bigotimes_{k<j}Z_k — parity string of the preceding modes, enforcing the anticommutation relations; Xj,YjX_j,Y_j — Pauli operators on qubit jj. The Bravyi–Kitaev encoding trades this O(n)\mathcal{O}(n)-length string for O(log⁡n)\mathcal{O}(\log n), at the price of a less transparent mapping. To train the circuit one needs gradients, and here VQE has a small miracle: for gates generated by an involution, the derivative is exact from two energy evaluations at shifted angles — no finite-difference step size to tune, no discretization bias.

∂E∂θj=12[E ⁣(θj+π2)−E ⁣(θj−π2)](7)\frac{\partial E}{\partial\theta_j}=\tfrac{1}{2}\Bigl[E\!\left(\theta_j+\tfrac{\pi}{2}\right)-E\!\left(\theta_j-\tfrac{\pi}{2}\right)\Bigr] \qquad (7)

where: ∂E/∂θj\partial E/\partial\theta_j — exact gradient component with respect to parameter θj\theta_j; θj±π/2\theta_j\pm\pi/2 — the same circuit run with that one angle shifted by a quarter turn. It looks like a finite difference but it is not: for a Pauli-generated gate the energy is an exact sinusoid in θj\theta_j, so this combination returns the true derivative. Finally, precision costs measurements. Because each Pauli expectation is a sample mean of ±1\pm1 outcomes, resolving the energy to a standard error ϵ\epsilon demands a shot budget that grows as the square of the coefficient one-norm — the origin of VQE's notorious measurement overhead.

N  ≳  1ϵ2(∑k∣ck∣ σk)2  ≤  1ϵ2(∑k∣ck∣)2,N=∑kNk(8)N\;\gtrsim\;\frac{1}{\epsilon^2}\Bigl(\sum_k|c_k|\,\sigma_k\Bigr)^{2}\;\le\;\frac{1}{\epsilon^2}\Bigl(\sum_k|c_k|\Bigr)^{2},\qquad N=\sum_k N_k \qquad (8)

where: NkN_k — shots allocated to Pauli term PkP_k; σk=1−⟨Pk⟩2≤1\sigma_k=\sqrt{1-\langle P_k\rangle^2}\le1 — its per-shot standard deviation; ϵ\epsilon — target standard error on the total energy; NN — total measurement budget. The 1/ϵ21/\epsilon^2 scaling is the shot-noise (standard quantum limit) tax that phase estimation, at 1/ϵ1/\epsilon, avoids.

Derivation AWhy the Rayleigh quotient can never fall below the ground energy

The single fact that makes VQE trustworthy is that its cost function is bounded below by the very quantity it estimates. The goal here is to prove that inequality cleanly — that for every trial state ∣ψ⟩|\psi\rangle the Rayleigh quotient ⟨ψ∣H∣ψ⟩/⟨ψ∣ψ⟩\langle\psi|H|\psi\rangle/\langle\psi|\psi\rangle is at least E0E_0, with equality only when ∣ψ⟩|\psi\rangle is a ground eigenstate. The proof is nothing but the spectral theorem plus the definition of the smallest eigenvalue, and it assumes nothing about the ansatz, optimizer, or hardware.

Step 1. Because HH is Hermitian, the spectral theorem gives an orthonormal eigenbasis {∣En⟩}\{|E_n\rangle\} with real eigenvalues EnE_n, ordered so that E0≤E1≤⋯E_0\le E_1\le\cdots. Expand the trial state in this basis, ∣ψ⟩=∑nan∣En⟩|\psi\rangle=\sum_n a_n|E_n\rangle with H∣En⟩=En∣En⟩H|E_n\rangle=E_n|E_n\rangle. Why: the eigenbasis is the one frame in which HH acts by simple multiplication, converting operator algebra into ordinary sums.

Step 2. Compute the numerator. Using orthonormality ⟨Em∣En⟩=δmn\langle E_m|E_n\rangle=\delta_{mn}, every cross term drops and ⟨ψ∣H∣ψ⟩=∑nEn∣an∣2\langle\psi|H|\psi\rangle=\sum_n E_n|a_n|^2. The denominator is ⟨ψ∣ψ⟩=∑n∣an∣2\langle\psi|\psi\rangle=\sum_n|a_n|^2. Why: orthonormality is what makes the double sum collapse to a single weighted sum of eigenvalues — the interference between distinct eigenstates cannot contribute to the energy.

Step 3. Bound each eigenvalue from below by the smallest one: since En≥E0E_n\ge E_0 for all nn, replacing every EnE_n by E0E_0 can only decrease the sum, giving ∑nEn∣an∣2≥E0∑n∣an∣2\sum_n E_n|a_n|^2\ge E_0\sum_n|a_n|^2. Why: this is the one inequality in the whole argument, and it is exactly where "there is a lowest eigenvalue" enters — the weights ∣an∣2|a_n|^2 are non-negative probabilities, so a term-by-term lower bound survives summation.

Step 4. Divide by ⟨ψ∣ψ⟩=∑n∣an∣2\langle\psi|\psi\rangle=\sum_n|a_n|^2 to reach ⟨ψ∣H∣ψ⟩/⟨ψ∣ψ⟩≥E0\langle\psi|H|\psi\rangle/\langle\psi|\psi\rangle\ge E_0. Equality holds iff all the weight sits on the ground level, i.e. an=0a_n=0 whenever En>E0E_n>E_0, which means ∣ψ⟩∝∣E0⟩|\psi\rangle\propto|E_0\rangle. Why: normalizing turns the weighted average of eigenvalues into a convex combination, and a convex combination of numbers all ≥E0\ge E_0 equals E0E_0 only when it puts no weight above the floor. *Verified: E(θ) ≥ E₀ for every θ; for the single-qubit H = Z + ½X the minimum of cos θ + ½ sin θ equals −√1.25 = −1.11803, exactly the smaller eigenvalue of H, confirming the bound is tight when the ansatz reaches the ground state.*

Go deeper · Derivation B · 2 minThe parameter-shift rule gives exact gradients from two circuit runs
Derivation BThe parameter-shift rule gives exact gradients from two circuit runs

Gradient-based optimizers need ∂E/∂θj\partial E/\partial\theta_j, but a quantum computer cannot differentiate a circuit symbolically — only run it and measure. The goal here is to show that when θj\theta_j enters through a gate e−iθjG/2e^{-i\theta_j G/2} whose generator is an involution (G2=IG^2=I, true for any Pauli generator), the energy is an exact sinusoid in θj\theta_j, so its derivative is captured without approximation by one difference of two shifted evaluations — no step size, no truncation error, no bias-versus-noise tradeoff.

Step 1. Expand the gate. Because G2=IG^2=I, the Taylor series of the exponential splits into even and odd parts that resum to e−iθG/2=cos⁡(θ/2) I−isin⁡(θ/2) Ge^{-i\theta G/2}=\cos(\theta/2)\,I-i\sin(\theta/2)\,G. Why: an involution behaves under exponentiation exactly like ii under eiθe^{i\theta}, so the gate is a rotation in the two-dimensional algebra spanned by II and GG — the source of every sinusoid that follows.

Step 2. Sandwich HH (dressed by all the other, fixed gates) between eiθG/2e^{i\theta G/2} and e−iθG/2e^{-i\theta G/2} and expand. The four resulting terms carry factors cos⁡2(θ/2)\cos^2(\theta/2), sin⁡2(θ/2)\sin^2(\theta/2), and sin⁡(θ/2)cos⁡(θ/2)\sin(\theta/2)\cos(\theta/2), which combine through the double-angle identities into a single harmonic: E(θ)=Acos⁡θ+Bsin⁡θ+CE(\theta)=A\cos\theta+B\sin\theta+C. Why: only the first harmonic can appear because GG has exactly two distinct eigenvalues (±1\pm1), so the parameter dependence lives entirely at frequency one — no higher harmonics are possible.

Step 3. Differentiate the exact form: E′(θ)=−Asin⁡θ+Bcos⁡θE'(\theta)=-A\sin\theta+B\cos\theta. Now evaluate the energy at the two shifted points a quarter turn apart. Using cos⁡(θ±π2)=∓sin⁡θ\cos(\theta\pm\tfrac{\pi}{2})=\mp\sin\theta and sin⁡(θ±π2)=±cos⁡θ\sin(\theta\pm\tfrac{\pi}{2})=\pm\cos\theta gives E(θ+π2)=−Asin⁡θ+Bcos⁡θ+CE(\theta+\tfrac{\pi}{2})=-A\sin\theta+B\cos\theta+C and E(θ−π2)=Asin⁡θ−Bcos⁡θ+CE(\theta-\tfrac{\pi}{2})=A\sin\theta-B\cos\theta+C. Why: the ±π/2\pm\pi/2 shift is chosen precisely so the cosine and sine trade places, isolating the derivative's ingredients.

Step 4. Subtract and halve: 12[E(θ+π2)−E(θ−π2)]=−Asin⁡θ+Bcos⁡θ=E′(θ)\tfrac{1}{2}[E(\theta+\tfrac{\pi}{2})-E(\theta-\tfrac{\pi}{2})]=-A\sin\theta+B\cos\theta=E'(\theta). The constant CC cancels and the result is the exact derivative. Why: because the energy truly is a pure first harmonic, two samples a quarter period apart determine its slope everywhere — a fact finite differences can only approximate. *Verified: for H = Z + ½X and E(θ) = cos θ + ½ sin θ evaluated at θ = 1.0, the parameter-shift estimate ½[E(1+π/2) − E(1−π/2)] = −0.57132 matches the analytic derivative −sin 1 + ½ cos 1 = −0.57132 to all quoted digits.*

Limiting cases

VQE is best understood at its boundary regimes, where the general machinery reduces to a single clean statement.

CaseResultMeaning
Exact ansatz, global optimizermin⁡θE(θ)=E0\min_{\boldsymbol{\theta}}E(\boldsymbol{\theta})=E_0If the ansatz family contains the true ground state and the optimizer finds the global minimum, VQE recovers E0E_0 exactly in the infinite-shot limit.
Infinite measurement shotsE^(θ)→⟨ψ(θ)∣H∣ψ(θ)⟩≥E0\widehat{E}(\boldsymbol{\theta})\to\langle\psi(\boldsymbol{\theta})|H|\psi(\boldsymbol{\theta})\rangle\ge E_0The estimator becomes deterministic and the rigorous variational bound holds; finite shots make it a random variable around this value.
Barren plateau (expressive random ansatz)Var[∂jE]∼2−n\mathrm{Var}[\partial_j E]\sim 2^{-n}For 2-design-like circuits the gradient variance vanishes exponentially in qubit number, so gradient-based training needs ∼2n\sim 2^n shots per step.
Shot-noise-limited precisionϵ∼1N∑k∣ck∣\epsilon\sim\dfrac{1}{\sqrt{N}}\sum_k|c_k|Sampling gives the standard-quantum-limit 1/N1/\sqrt{N} scaling, worse than phase estimation's Heisenberg-limited 1/N1/N.
Single trivial (L=0L=0) layerE=⟨0∣⊗nH∣0⟩⊗nE=\langle 0|^{\otimes n}H|0\rangle^{\otimes n}With no parameters the circuit returns the reference (Hartree–Fock) energy — the classical starting point VQE aims to beat.

Worked examples

1 · A single-qubit VQE, worked end to end

Let H=Z+12XH=Z+\tfrac{1}{2}X act on one qubit. Use the one-parameter ansatz ∣ψ(θ)⟩=Ry(θ)∣0⟩=cos⁡θ2∣0⟩+sin⁡θ2∣1⟩|\psi(\theta)\rangle=R_y(\theta)|0\rangle=\cos\tfrac{\theta}{2}|0\rangle+\sin\tfrac{\theta}{2}|1\rangle. Find the variational energy, its minimum, the optimal angle, and check the parameter-shift gradient.

Step 1. Build the cost function from Pauli expectations. For the RyR_y ansatz, direct calculation gives ⟨Z⟩=cos⁡θ\langle Z\rangle=\cos\theta and ⟨X⟩=sin⁡θ\langle X\rangle=\sin\theta. By linearity (equation 3), E(θ)=1⋅⟨Z⟩+12⋅⟨X⟩=cos⁡θ+12sin⁡θE(\theta)=1\cdot\langle Z\rangle+\tfrac{1}{2}\cdot\langle X\rangle=\cos\theta+\tfrac{1}{2}\sin\theta. This is the exact single-harmonic form the parameter-shift derivation predicted.

Step 2. Minimize analytically. Write E(θ)=Rcos⁡(θ−ϕ)E(\theta)=R\cos(\theta-\phi) with amplitude R=12+(12)2=1.25=1.11803R=\sqrt{1^2+(\tfrac12)^2}=\sqrt{1.25}=1.11803. The minimum is −R-R, so Emin⁡=−1.25=−1.11803\boxed{E_{\min}=-\sqrt{1.25}=-1.11803}, reached when cos⁡(θ−ϕ)=−1\cos(\theta-\phi)=-1. Setting E′(θ)=−sin⁡θ+12cos⁡θ=0E'(\theta)=-\sin\theta+\tfrac12\cos\theta=0 gives tan⁡θ=12\tan\theta=\tfrac12; the minimizing root in [0,2π)[0,2\pi) is θ∗=π+arctan⁡12=3.605 rad\theta^{*}=\pi+\arctan\tfrac12=\mathbf{3.605\ rad}.

Step 3. Confirm against exact diagonalization. As a 2×22\times2 matrix H=(11212−1)H=\begin{pmatrix}1&\tfrac12\\[2pt]\tfrac12&-1\end{pmatrix} has eigenvalues ±1+14=±1.11803\pm\sqrt{1+\tfrac14}=\pm1.11803. The variational minimum equals the smaller eigenvalue exactly — the ansatz is expressive enough to reach the ground state, so the bound is tight. VQE returns E0=−1.11803E_0=-1.11803, the true ground energy.

Step 4. Test the parameter-shift rule at θ=1.0\theta=1.0. The shifted energies are E(1+π2)E(1+\tfrac{\pi}{2}) and E(1−π2)E(1-\tfrac{\pi}{2}); their half-difference is 12[E(2.5708)−E(−0.5708)]=−0.57132\tfrac12[E(2.5708)-E(-0.5708)]=\mathbf{-0.57132}. The analytic slope is E′(1)=−sin⁡1+12cos⁡1=−0.84147+0.27015=−0.57132E'(1)=-\sin 1+\tfrac12\cos 1=-0.84147+0.27015=\mathbf{-0.57132}. The two agree to every quoted digit — the shift rule is exact, not a finite-difference estimate.

2 · The measurement budget for chemical accuracy

A VQE energy is only useful in chemistry if it reaches "chemical accuracy," ϵ=1.6 mHartree≈1.6×10−3\epsilon=1.6\ \mathrm{mHartree}\approx1.6\times10^{-3}, the error at which reaction rates come out right. Estimate the shot budget for the H=Z+12XH=Z+\tfrac12X problem above, and use it to explain why measurement, not gate depth, is VQE's bottleneck.

Step 1. Read off the coefficient one-norm. Here the Pauli coefficients are cZ=1c_Z=1 and cX=12c_X=\tfrac12, so ∑k∣ck∣=1.5\sum_k|c_k|=1.5. Taking the worst-case σk=1\sigma_k=1 (no knowledge of the expectation values yet) turns the shot-budget bound of equation (8) into N≈(∑k∣ck∣)2/ϵ2N\approx(\sum_k|c_k|)^2/\epsilon^2. That budget is the closed form of a Lagrange-multiplier optimization: to reach total error ϵ\epsilon with the fewest shots, allocate Nk∝∣ck∣σkN_k\propto|c_k|\sigma_k to each term.

Step 2. Evaluate. N≈(1.5)2/(1.6×10−3)2=2.25/2.56×10−6=8.79×105N\approx(1.5)^2/(1.6\times10^{-3})^2=2.25/2.56\times10^{-6}=\mathbf{8.79\times10^{5}} shots per energy evaluation. For a toy one-qubit problem that is already nearly a million measurements just to place one point on the landscape to chemical accuracy.

Step 3. Scale it up. A real molecule such as BeH2\mathrm{BeH_2} has hundreds of Pauli terms and a one-norm of tens of Hartree, so (∑k∣ck∣)2(\sum_k|c_k|)^2 balloons by four to six orders of magnitude; times the hundreds of energy evaluations an optimizer needs, one dissociation-curve point can demand 10910^{9}–101210^{12} shots. This 1/ϵ21/\epsilon^2 measurement wall — not circuit depth — is why practical VQE research obsesses over commuting-term grouping, classical shadows, and low-variance estimators.

Physics of the portrait

VQE is a very modern algorithm resting on a very old idea. The engine under its hood is the variational principle — the observation, more than a century old, that you can bound a system's ground energy from above with any trial wavefunction and then improve the bound by tuning it. Everything VQE does is a quantum-hardware realization of that single move.

Origins · 1 minvariational principle and how this was discovered
VP
variational principle
Rayleigh 1870 · Ritz 1909 · to VQE 2014 · the Rayleigh quotient behind every eigensolver, classical and quantum

The variational principle is one of the most productive ideas in physics, and it long predates the machines that now run it. Lord Rayleigh, in the 1870s, showed that a vibrating system's natural frequency could be bounded by an energy ratio evaluated on any admissible trial shape — the Rayleigh quotient — and that refining the shape tightened the bound. Walther Ritz, in 1909, made it systematic: expand the trial function in basis functions with free coefficients and minimize the quotient over them, converting an eigenvalue problem into ordinary calculus. The Rayleigh–Ritz method became the backbone of quantum chemistry — Egil Hylleraas used it in 1929 to compute the ground-state energy of helium to spectroscopic accuracy with a handful of cleverly chosen parameters, an early triumph that showed the new quantum mechanics was quantitatively exact. For most of a century the trial function was a formula and the minimization ran on paper or on classical computers whose cost exploded with the number of correlated electrons. The Variational Quantum Eigensolver, introduced by Peruzzo, McClean, and collaborators in 2014, is the same principle with a quantum twist: the trial function is now a quantum state prepared by a circuit that a classical computer could not even write down, and the Rayleigh quotient is measured rather than integrated. The free coefficients are gate angles; the minimization is a hybrid loop. Ritz's calculus of variations, running on a chip.

History
1870
— Rayleigh bounds vibration frequencies with a trial-function energy ratio
1909
— Walther Ritz formalizes the variational method for eigenvalue problems
1929
— Hylleraas computes the helium ground state variationally to spectroscopic accuracy
2014
— Peruzzo, McClean and colleagues run the first VQE on a photonic chip
2018
— McClean and colleagues identify barren plateaus in the training landscape

Going further

Texts and reviews. Nielsen & Chuang, Quantum Computation and Quantum Information, §2.2 (spectral theorem, expectation values) and §5 (phase estimation, the deep-circuit alternative) · Preskill, Quantum Computing in the NISQ era and beyond, Quantum 2, 79 (2018) · Cerezo et al., Variational Quantum Algorithms, Nat. Rev. Phys. 3, 625 (2021) · McClean, Romero, Babbush, Aspuru-Guzik, The theory of variational hybrid quantum-classical algorithms, New J. Phys. 18, 023023 (2016).

Primary sources. Peruzzo et al., A variational eigenvalue solver on a photonic quantum processor, Nat. Commun. 5, 4213 (2014) · O'Malley et al., Scalable quantum simulation of molecular energies, Phys. Rev. X 6, 031007 (2016) · McClean et al., Barren plateaus in quantum neural network training landscapes, Nat. Commun. 9, 4812 (2018) · Kandala et al., Hardware-efficient variational quantum eigensolver for small molecules, Nature 549, 242 (2017).

Provenance. The variational bound is Rayleigh–Ritz, recalled from Hermitian operators and the spectral theorem (Nielsen & Chuang §2.2). The Pauli decomposition, energy linearity, parameter-shift rule, and shot-allocation scaling follow McClean et al. NJP 18 (2016) and Cerezo et al. Nat. Rev. Phys. 3 (2021); the barren-plateau 2−n2^{-n} variance is McClean et al. Nat. Commun. 9, 4812 (2018); the Jordan–Wigner mapping is standard (Bravyi–Kitaev the low-weight alternative). A verifier should recompute: the Rayleigh-quotient lower bound; the single-harmonic energy underlying the parameter-shift rule and the numeric check 12[E(1+π2)−E(1−π2)]=−0.57132\tfrac12[E(1{+}\tfrac{\pi}{2})-E(1{-}\tfrac{\pi}{2})]=-0.57132; the single-qubit ground energy −1.25=−1.11803-\sqrt{1.25}=-1.11803; the shot budget (1.5)2/(1.6×10−3)2=8.79×105(1.5)^2/(1.6\times10^{-3})^2=8.79\times10^5; and the chemical-accuracy figure 1.6 mHartree1.6\ \mathrm{mHartree}.

Builds on: Hermitian operators and the spectral theorem · Pauli operators and the Pauli group · Expectation values and projective measurement · Parameterized quantum circuits · Second quantization and fermionic operators. Leads to: Quantum Approximate Optimization Algorithm (QAOA) · Magic State Distillation · Quantum Phase Estimation.

Ask the tutor: ask me to walk you through the single-qubit VQE by hand — building E(θ)E(\theta) from ⟨Z⟩\langle Z\rangle and ⟨X⟩\langle X\rangle, minimizing it, and checking the parameter-shift gradient — or to explain why barren plateaus and the measurement wall are the two exponentials VQE has to outrun.

Is VQE guaranteed to find the true ground-state energy?

No — this is the caveat to internalize. VQE returns only an upper bound E(θ)≥E0E(\boldsymbol{\theta})\ge E_0, and two things can hold it above the true floor. First, expressibility: if the ansatz simply cannot produce the ground state — too few layers, the wrong symmetry sector, missing excitations — no amount of optimization reaches E0E_0. Second, optimization: even an expressive ansatz can trap the optimizer in a local minimum or a barren plateau. VQE tells you the best energy achievable within your ansatz family as explored by your optimizer — a genuine bound, not a certificate of exactness. Comparing against a richer ansatz or a classical benchmark is how practitioners gain confidence.

Why not just use quantum phase estimation, which is exact?

Quantum phase estimation (QPE) extracts eigenvalues to Heisenberg-limited 1/ϵ1/\epsilon precision, but pays with circuit depth: it implements controlled time evolution e−iHte^{-iHt} over many ancilla-controlled steps, building coherent circuits far deeper than today's qubits survive. VQE deliberately trades depth for shallowness — decoherence has little time to act, and the cost migrates to classical optimization and measurement repetition. The price is precision: VQE's shot-noise-limited 1/N1/\sqrt{N} is quadratically worse than QPE's 1/N1/N. Rule of thumb: QPE is the fault-tolerant answer for when error-corrected qubits exist; VQE is the near-term answer for the noisy hardware we actually have.

What exactly is the ansatz, and how do I choose one?

The ansatz is the parameterized circuit U(θ)U(\boldsymbol{\theta}) that turns a fixed reference into the trial state you optimize, and it is where VQE's central tension lives. Chemistry-inspired ansätze like UCCSD are built from the molecule's excitation operators, so they respect physics and their parameters carry meaning — but they compile to deep circuits. Hardware-efficient ansätze are shallow stacks of native rotations and entanglers fitted to a chip — but they ignore problem structure and are the most prone to barren plateaus. The frontier is in between: symmetry-preserving ansätze that conserve particle number and spin, and adaptive methods (ADAPT-VQE) that grow the circuit operator by operator. Choosing an ansatz is choosing a tradeoff between expressibility, trainability, and depth.

Where does the barren plateau come from, physically?

From too much expressiveness. If the ansatz is deep and unstructured enough that its outputs scatter like Haar-random states (a "2-design"), averaging the gradient over the circuit's randomness gives zero mean and a variance that shrinks as 2−n2^{-n}. Intuitively, a random state has essentially the same near-maximal-entropy energy as any nearby random state, so the landscape is exponentially flat almost everywhere with no local signal pointing downhill. Resolving a gradient against shot noise then costs ∼2n\sim 2^n measurements per step, defeating scaling. Mitigations reduce expressiveness or localize the cost: shallow problem-tailored ansätze, local rather than global cost functions, layerwise training, and initializations near the identity.

Why does measurement dominate the runtime instead of gate operations?

Because the energy is not a single measurable quantity — it is a weighted sum of many Pauli expectations, each estimated by sampling ±1\pm1 outcomes with shot noise falling only as 1/N1/\sqrt{N}. Equation (8) shows the total scales as (∑k∣ck∣)2/ϵ2(\sum_k|c_k|)^2/\epsilon^2. For chemical accuracy ϵ=1.6\epsilon=1.6 mHartree the 1/ϵ21/\epsilon^2 alone is ∼4×105\sim4\times10^5; multiply by a molecular one-norm of tens of Hartree squared and by the hundreds of energy evaluations an optimizer needs, and one result can require 10910^9–101210^{12} shots. Each shot is a full state preparation plus readout, so measurement — not the shallow gate sequence — sets the wall-clock time. Hence the intense interest in commuting-term grouping, classical shadows, and variance-reduced estimators.

How does a fermionic molecule become a qubit problem at all?

Through a fermion-to-qubit encoding. The molecular Hamiltonian is first written in second quantization from creation and annihilation operators a†,aa^\dagger,a, which obey fermionic anticommutation relations qubits do not natively share. The Jordan–Wigner transform (equation 6) fixes this by mapping each mode to one qubit and attaching a string of Pauli-ZZ operators on all preceding qubits; that parity string reproduces the minus signs of fermionic exchange. The Bravyi–Kitaev encoding does the same with only O(log⁡n)\mathcal{O}(\log n)-length strings, lowering the weight of the resulting terms. After the mapping, H=∑kckPkH=\sum_k c_k P_k is ready for the VQE loop and the reference is the Hartree–Fock determinant as a computational-basis string.

Where this leads