Chapter 19: Variational Quantum Eigensolver (VQE) — Finding Ground State Energies for Chemistry and Materials Science

"The ground state energy of a molecule is the most important number in chemistry. VQE puts it within reach of a quantum computer — even a noisy one."

Learning Objectives

After completing this chapter, you will be able to:

• State the variational principle and explain why it provides an upper bound on the ground state energy.
• Derive the variational principle from first principles and prove the quadratic error bound.
• Design ansatz circuits for VQE, including hardware-efficient and chemically inspired (UCCSD) ansatzes, and compare their trade-offs.
• Map a molecular Hamiltonian to a sum of Pauli strings and group commuting observables for efficient measurement.
• Implement a complete VQE workflow in Qiskit Nature for the H₂ molecule, including dissociation curves and active space reductions.
• Compare classical optimizers (COBYLA, SPSA, gradient-based) and select appropriate ones for VQE.
• Apply error mitigation techniques — readout error mitigation and zero-noise extrapolation — to improve VQE results on real hardware.
• Extend VQE to larger molecules (LiH, H₂O) and understand the scaling challenges.
• Diagnose and mitigate barren plateaus in variational landscapes.
• Understand the historical development of VQE and its relationship to classical quantum chemistry methods.

19.1 The Central Problem of Quantum Chemistry

19.1.1 Why Ground State Energies Matter

The ground state energy of a molecule is not merely an abstract number — it is the key that unlocks virtually every chemical and materials property. Reaction rates are determined by energy differences between reactants and products (activation barriers). Bond strengths are determined by dissociation energies. Spectroscopic transitions are determined by energy gaps between electronic states. Drug binding affinities are determined by interaction energies between molecules. If we could compute ground state energies accurately and efficiently, we could predict chemical behavior from first principles, without recourse to experiment.

This is the promise of quantum chemistry: to compute molecular properties by solving the electronic Schrödinger equation. The challenge is that the computational cost of exact solution grows exponentially with system size.

19.1.2 The Electronic Schrödinger Equation

The behavior of electrons in a molecule is governed by the time-independent Schrödinger equation:

$$H |\Psi\rangle = E |\Psi\rangle$$

where $H$ is the molecular Hamiltonian operator and $|\Psi\rangle$ is the electronic wavefunction. The ground state $|\Psi_0\rangle$ is the eigenstate with the lowest eigenvalue $E_0$:

$$E_0 = \min_{|\Psi\rangle} \frac{\langle\Psi|H|\Psi\rangle}{\langle\Psi|\Psi\rangle}$$

For a molecule with $N$ electrons and $M$ nuclei, the electronic Hamiltonian in atomic units is:

$$H = -\sum_{i=1}^{N} \frac{1}{2}\nabla_i^2 - \sum_{i=1}^{N}\sum_{A=1}^{M} \frac{Z_A}{|\mathbf{r}_i - \mathbf{R}_A|} + \sum_{i

The first term is the kinetic energy of the electrons, the second is the electron-nuclear attraction, the third is the electron-electron repulsion, and the fourth is the nuclear-nuclear repulsion (a constant for fixed nuclear positions within the Born-Oppenheimer approximation).

19.1.3 The Exponential Wall

For a molecule with $N$ electrons, the wavefunction lives in a Hilbert space of dimension exponential in $N$. In a basis set with $K$ spatial orbitals (each holding up to 2 electrons), the number of possible electron configurations is $\binom{2K}{N}$, which grows combinatorially. For example, the FeMoco active site of nitrogenase — the "holy grail" of quantum chemistry — has approximately 54 electrons in 54 orbitals, giving $\binom{108}{54} \approx 10^{32}$ configurations. Exact diagonalization (Full Configuration Interaction, FCI) of the Hamiltonian requires constructing and diagonalizing matrices of this dimension, which is utterly intractable on any classical computer.

Common Misconception: "Quantum computers can solve the Schrödinger equation instantly." No — they still require clever algorithms. The advantage is that a quantum computer represents the wavefunction natively in its exponentially large Hilbert space. Simulating a molecule with $N$ spin-orbitals requires only $O(N)$ qubits, not $O(2^N)$ classical memory. But extracting the ground state energy from this representation still requires algorithmic ingenuity — this is where VQE comes in.

19.1.4 Classical Methods and Their Limitations

Before VQE, classical quantum chemistry developed a hierarchy of approximate methods:

The gap between "tractable" and "accurate" methods is enormous for strongly correlated systems: transition metal complexes, bond-breaking processes, superconducting materials, and catalytic reactions. These are the problems VQE targets.

19.1.5 Why a Quantum Approach?

This is the killer app for quantum computing. A quantum computer represents the wavefunction natively in its exponentially large Hilbert space. Simulating a molecule with $N$ spin-orbitals requires only $O(N)$ qubits, not $O(2^N)$ classical memory. The challenge is extracting the ground state energy from a noisy, intermediate-scale quantum processor. VQE is the leading approach.

Recurring Theme: Quantum is Linear Algebra, Not Magic. VQE works because it exploits a fundamental property of linear algebra: the Rayleigh quotient $\langle\psi|H|\psi\rangle/\langle\psi|\psi\rangle$ provides an upper bound on the smallest eigenvalue. This is a mathematical theorem, not a quantum mystery. The quantum advantage comes from the fact that we can compute $\langle\psi|H|\psi\rangle$ on a quantum computer using resources that scale polynomially with the system size, even though $H$ is exponentially large as a matrix.

19.2 The Variational Principle

The variational principle is the foundation of VQE. It states that for any trial wavefunction $|\psi(\boldsymbol{\theta})\rangle$ parameterized by real parameters $\boldsymbol{\theta} = (\theta_1, \ldots, \theta_p)$, the expectation value of the Hamiltonian provides an upper bound on the true ground state energy:

$$E_0 \leq \frac{\langle\psi(\boldsymbol{\theta})|H|\psi(\boldsymbol{\theta})\rangle}{\langle\psi(\boldsymbol{\theta})|\psi(\boldsymbol{\theta})\rangle}$$

If $|\psi(\boldsymbol{\theta})\rangle$ is normalized (as all quantum states are), the denominator is 1 and we have:

$$E_0 \leq \langle\psi(\boldsymbol{\theta})|H|\psi(\boldsymbol{\theta})\rangle \equiv E(\boldsymbol{\theta})$$

19.2.1 Proof of the Variational Principle

Proof. Expand $|\psi(\boldsymbol{\theta})\rangle$ in the eigenbasis of $H$: $|\psi(\boldsymbol{\theta})\rangle = \sum_j c_j |\Psi_j\rangle$ where $H|\Psi_j\rangle = E_j |\Psi_j\rangle$ and $E_0 \leq E_1 \leq \cdots$. Then:

$$\langle\psi|H|\psi\rangle = \sum_{j,k} c_j^* c_k \langle\Psi_j|H|\Psi_k\rangle = \sum_j |c_j|^2 E_j$$

Since $E_j \geq E_0$ for all $j$ and $\sum_j |c_j|^2 = 1$ (normalization), we have:

$$\langle\psi|H|\psi\rangle = \sum_j |c_j|^2 E_j \geq E_0 \sum_j |c_j|^2 = E_0$$

with equality if and only if $|\psi\rangle = |\Psi_0\rangle$ (up to a phase). $\square$

19.2.2 The Quadratic Error Bound: Why VQE Can Achieve Chemical Accuracy

A crucial property of the variational principle is that the energy error is second order in the wavefunction error. If the trial state has a small overlap with the first excited state:

$$|\psi\rangle = \sqrt{1 - \epsilon^2}|\Psi_0\rangle + \epsilon|\Psi_1\rangle$$

then the energy expectation value is:

$$\langle\psi|H|\psi\rangle = (1 - \epsilon^2)E_0 + \epsilon^2 E_1 = E_0 + \epsilon^2(E_1 - E_0)$$

The energy error is $\Delta E = \epsilon^2(E_1 - E_0)$, which is $O(\epsilon^2)$, not $O(\epsilon)$. This means that even a somewhat imperfect wavefunction can yield a very accurate energy — a property that is fundamental to why VQE works as well as it does.

Try It Yourself: Verify this quadratic bound numerically. For a two-level system with $H = \text{diag}(E_0, E_1)$, compute $\langle\psi|H|\psi\rangle$ for $|\psi\rangle = \cos\theta|0\rangle + \sin\theta|1\rangle$. Plot the energy error $\Delta E = \langle\psi|H|\psi\rangle - E_0$ against the wavefunction error $\epsilon = \sin\theta$ and confirm the parabolic relationship.

19.2.3 Extended Derivation: The Variational Principle for Excited States

The variational principle can be extended to excited states. The hylleraas-undheim theorem states that if we have a set of orthogonal trial states $\{|\psi_1\rangle, |\psi_2\rangle, \ldots, |\psi_k\rangle\}$, then the $k$-th lowest eigenvalue of the variational energies satisfies:

$$E(\psi_k) \geq E_k$$

where $E_k$ is the $k$-th exact eigenvalue. This follows from the fact that the variational energies are the eigenvalues of the Hamiltonian projected onto the span of the trial states, and the Cauchy interlacing theorem guarantees that these projected eigenvalues interlace with the exact eigenvalues.

Step-by-step derivation:

Let $S = \text{span}\{|\psi_1\rangle, \ldots, |\psi_k\rangle\}$ be the $k$-dimensional subspace spanned by the trial states. The projected Hamiltonian is $P_S H P_S$, where $P_S$ is the projector onto $S$. The eigenvalues of $P_S H P_S$ are the variational energies $E^{(1)} \leq E^{(2)} \leq \cdots \leq E^{(k)}$.

By the min-max theorem (Courant-Fischer), for the $j$-th eigenvalue:

$$E_j = \min_{\dim(V)=j} \max_{|\psi\rangle \in V} \langle\psi|H|\psi\rangle$$

Since $S$ is a $k$-dimensional subspace, the variational eigenvalue $E^{(j)}$ satisfies $E^{(j)} \geq E_j$ for all $j \leq k$. In particular, $E^{(k)} \geq E_k$.

This result has a profound implication for VQE: we can target excited states by optimizing in a subspace orthogonal to previously found ground states. The subspace-search VQE algorithm exploits this by:

1. Finding the ground state $|\psi_0\rangle$ with energy $E_0$ via standard VQE.
2. Constraining the next search to the subspace orthogonal to $|\psi_0\rangle$, yielding $|\psi_1\rangle$ with $E_1 \geq E_1^{\text{exact}}$.
3. Continuing for higher excited states.

Common Misconception: "The variational principle only applies to ground states." While the basic variational principle gives an upper bound only on the ground state energy, the Hylleraas-Undheim theorem extends this to excited states. Moreover, the variance $\langle H^2\rangle - \langle H\rangle^2$ provides a quality metric for any state — not just the ground state — and can be used as an additional objective function in excited-state VQE.

Worked Example 19.1: Variational Principle for a Two-Level System

Consider the Hamiltonian $H = \begin{pmatrix} -1 & 0.5 \\ 0.5 & 2 \end{pmatrix}$ with eigenvalues $E_0 \approx -1.118$ and $E_1 \approx 2.118$.

Using the trial state $|\psi(\theta)\rangle = \cos\theta|0\rangle + \sin\theta|1\rangle$:

$$E(\theta) = \cos^2\theta \cdot (-1) + \sin^2\theta \cdot 2 + 2\cos\theta\sin\theta \cdot 0.5$$

$$= -\cos^2\theta + 2\sin^2\theta + \sin(2\theta) \cdot 0.5$$

Minimizing: $dE/d\theta = 0$ gives $\theta^* \approx 0.224$ radians, yielding $E(\theta^*) \approx -1.118$, which equals the exact ground state energy (to three decimal places). The variational bound $E(\theta) \geq E_0$ is satisfied for all $\theta$.

The variational principle transforms the eigenvalue problem into an optimization problem:

$$\boldsymbol{\theta}^* = \arg\min_{\boldsymbol{\theta}} E(\boldsymbol{\theta})$$

The quantum computer evaluates $E(\boldsymbol{\theta})$ for given parameters; a classical optimizer updates $\boldsymbol{\theta}$ to minimize $E$. This hybrid loop is the essence of VQE.

┌──────────────────────────────────────────────────────────────┐
│                     CLASSICAL OPTIMIZER                       │
│   θ ← θ - η∇E(θ)    or    θ ← COBYLA/SPSA step               │
│                          │                                    │
│          θ               │           E(θ) = ⟨ψ(θ)|H|ψ(θ)⟩     │
│  ┌───────────────────────┴──────────────────────────────┐     │
│  │                 QUANTUM PROCESSOR                     │     │
│  │                                                       │     │
│  │  |0⟩^n ──[ U(θ) ]── Measure in Pauli bases ──→ ⟨H⟩   │     │
│  │                                                       │     │
│  │  Repeat for each Pauli term in H                      │     │
│  └──────────────────────────────────────────────────────┘     │
└──────────────────────────────────────────────────────────────┘

19.2.4 Historical Context: From Ritz to Peruzzo

The variational principle has deep roots in physics and mathematics. Walter Ritz introduced the variational method in 1908 for solving eigenvalue problems in classical mechanics. It became a cornerstone of quantum mechanics after the Schrödinger equation was formulated in 1926. The Ritz variational method — choosing a parameterized trial wavefunction and minimizing $\langle\psi|H|\psi\rangle$ — has been the workhorse of computational chemistry for nearly a century.

VQE, introduced by Peruzzo et al. in 2014, is the quantum realization of this classical idea. The key insight is that while $\langle\psi|H|\psi\rangle$ is exponentially expensive to compute classically, a quantum computer can estimate it efficiently by preparing $|\psi(\boldsymbol{\theta})\rangle = U(\boldsymbol{\theta})|0\rangle^{\otimes n}$ and measuring in the appropriate basis. The classical optimization loop is unchanged — only the energy evaluation is outsourced to the quantum processor.

19.3 The Ansatz: Designing the Trial Wavefunction

The ansatz $U(\boldsymbol{\theta})$ is a parameterized quantum circuit that prepares the trial state:

$$|\psi(\boldsymbol{\theta})\rangle = U(\boldsymbol{\theta}) |0\rangle^{\otimes n}$$

The choice of ansatz is the most critical design decision in VQE. It must be expressive enough to contain (or closely approximate) the true ground state, yet shallow enough to run on NISQ hardware. This is the fundamental expressibility-noise trade-off: more expressive ansatzes are deeper and accumulate more noise.

19.3.1 Hardware-Efficient Ansatz

The hardware-efficient ansatz is designed to respect the native gate set and connectivity of the target device. It consists of alternating layers of single-qubit rotations and entangling gates:

Layer structure (repeated d times):

q_0: ──[R_y(θ_0)]──[R_z(θ_1)]──●────────────────────────
                                  │
q_1: ──[R_y(θ_2)]──[R_z(θ_3)]──┼──●─────────────────────
                                  │  │
q_2: ──[R_y(θ_4)]──[R_z(θ_5)]──┼──┼──●──────────────────
                                  │  │  │
q_3: ──[R_y(θ_6)]──[R_z(θ_7)]──┼──┼──┼──────────────────
                                  │  │  │
...  (entangling pattern repeats with different parameters)

Each qubit receives a parameterized $R_y(\theta) R_z(\theta')$ rotation, followed by a ladder of CNOT gates (or the device's native two-qubit gate). The circuit is repeated $d$ times (the "depth" or number of "repetitions"), with independent parameters in each layer.

Advantages: Respects hardware connectivity, shallow, few parameters per layer. Disadvantages: No chemical intuition; may require many layers to reach chemical accuracy; prone to barren plateaus at scale.

Worked Example 19.2: Counting Parameters in a Hardware-Efficient Ansatz

For an $n$-qubit hardware-efficient ansatz with depth $d$: - Each layer has $2n$ single-qubit rotation parameters (one $R_y$ and one $R_z$ per qubit). - Total parameters: $p = 2nd$.

For a 4-qubit, 3-layer ansatz: $p = 2 \times 4 \times 3 = 24$ parameters. This is already a non-trivial optimization problem — a 24-dimensional landscape with noise.

Common Misconception: "More parameters always means a better ansatz." In principle, a more expressive ansatz can represent a wider range of states. But in practice, more parameters means a harder optimization problem. A compact, chemically motivated ansatz with 10 parameters may outperform a hardware-efficient ansatz with 100 parameters because the former's landscape has a clear global minimum while the latter has a barren plateau.

19.3.2 Symmetry-Informed Ansatzes

Before discussing UCCSD, it is worth noting that molecular symmetries can drastically reduce the size of the ansatz. The electronic Hamiltonian commutes with several symmetry operators:

• Particle number: $N_e = \sum_p a_p^\dagger a_p$ (total electron count is conserved).
• Spin parity: $S_z = \sum_p (-1)^p a_p^\dagger a_p$ and $S^2$ (total spin is conserved).
• Spatial symmetries: Point group symmetries of the molecule (e.g., $D_{2h}$ for linear molecules).

By restricting the ansatz to preserve these symmetries, we eliminate parameters that correspond to symmetry-breaking excitations. For H₂ in the STO-3G basis, spin symmetry reduces UCCSD from 4 parameters to just 1 — the double excitation from occupied to virtual orbitals with the same spin configuration.

This symmetry reduction is not merely a computational convenience. Symmetry-breaking in the ansatz can lead to convergence to unphysical states, spurious local minima, and increased measurement noise. Preserving symmetries ensures that the variational search remains within the physically relevant subspace.

Try It Yourself: For the H₂ molecule, enumerate all possible single and double excitations from the Hartree-Fock reference. Show that spin and spatial symmetries reduce the count from 4 single excitations and 1 double excitation to just 1 parameter (the spin-adapted double excitation). Verify that this single parameter captures the full correlation energy at equilibrium geometry.

19.3.3 Unitary Coupled Cluster (UCCSD)

The gold standard of classical quantum chemistry is Coupled Cluster with Singles and Doubles (CCSD). Its unitary variant, UCCSD, is the most widely used chemically inspired ansatz for VQE.

The UCCSD ansatz is:

$$|\psi(\boldsymbol{\theta})\rangle = e^{\hat{T} - \hat{T}^\dagger} |\Phi_{\text{HF}}\rangle$$

where $|\Phi_{\text{HF}}\rangle$ is the Hartree-Fock reference state and $\hat{T} = \hat{T}_1 + \hat{T}_2$ is the cluster operator:

$$\hat{T}_1 = \sum_{i \in \text{occ}} \sum_{a \in \text{virt}} \theta_i^a \, a_a^\dagger a_i$$

$$\hat{T}_2 = \sum_{i>j}^{occ} \sum_{a>b}^{virt} \theta_{ij}^{ab} \, a_a^\dagger a_b^\dagger a_j a_i$$

Here $a^\dagger$ and $a$ are fermionic creation and annihilation operators, "occ" denotes occupied orbitals, and "virt" denotes virtual (unoccupied) orbitals. The parameters $\theta_i^a, \theta_{ij}^{ab}$ are the coupled cluster amplitudes.

Why unitary? Standard coupled cluster uses the non-unitary exponential $e^{\hat{T}}$, which is not implementable on a quantum computer (unitary operations only). The unitary form $e^{\hat{T} - \hat{T}^\dagger}$ ensures the ansatz is a valid quantum operation.

Derivation: The Trotterization of UCCSD

The exponential $e^{\hat{T} - \hat{T}^\dagger}$ cannot be directly implemented on a quantum computer. It must be approximated via Trotterization (Suzuki-Trotter decomposition):

$$e^{\hat{T} - \hat{T}^\dagger} \approx \prod_k e^{\theta_k (\hat{\tau}_k - \hat{\tau}_k^\dagger)}$$

where each $\hat{\tau}_k$ is a single or double excitation operator. Each factor $e^{\theta_k (\hat{\tau}_k - \hat{\tau}_k^\dagger)}$ is then mapped to a quantum circuit via the Jordan-Wigner or Bravyi-Kitaev transformation.

For a first-order Trotter decomposition with $N_t$ excitation operators, the error is $O(\theta^2/N_t)$, where $\theta = \max_k |\theta_k|$. Higher-order Trotter formulas reduce this error at the cost of deeper circuits.

Advantages: Chemically motivated, systematically improvable, compact parameterization. Disadvantages: Deep circuits (each excitation operator requires $O(n)$ gates), not hardware-aware.

Worked Example 19.3: UCCSD Circuit for H₂

For H₂ in the STO-3G basis (2 electrons in 4 spin-orbitals), the occupied orbitals are indices 0 and 1, and the virtual orbitals are indices 2 and 3. The number of single excitations is $2 \times 2 = 4$, and the number of double excitations is $1$ (only one pair of occupied orbitals and one pair of virtual orbitals). However, spin and spatial symmetries reduce the relevant excitations. The dominant excitation is:

$$\hat{T}_2^{02,13} |\Phi_{\text{HF}}\rangle = a_2^\dagger a_3^\dagger a_1 a_0 |\Phi_{\text{HF}}\rangle$$

After Jordan-Wigner transformation, this single double excitation maps to a relatively compact circuit. The full UCCSD ansatz for H₂ has only 1 parameter (the double excitation amplitude), making it extremely efficient.

Worked Example 19.3b: UCCSD Circuit Depth Analysis

Let us analyze the circuit depth of the UCCSD ansatz for a general molecule with $N$ occupied and $M$ virtual spin-orbitals.

The number of amplitudes (parameters) in UCCSD is: - Singles: $N_{\text{occ}} \times N_{\text{virt}}$ (from each occupied to each virtual spin-orbit) - Doubles: $\binom{N_{\text{occ}}}{2} \times \binom{N_{\text{virt}}}{2}$ (pairs of occupied to pairs of virtual)

For H₂O with 10 electrons in 14 spin-orbitals (7 spatial orbitals), after freezing the oxygen 1s orbital (2 electrons in 1 spatial orbital), the active space has 8 electrons in 12 spin-orbitals. This gives: - Singles: $8 \times 4 = 32$ amplitudes - Doubles: $\binom{8}{2} \times \binom{4}{2} = 28 \times 6 = 168$ amplitudes - Total: 200 amplitudes

Each excitation operator, after Jordan-Wigner transformation, requires $O(n)$ CNOT gates due to the parity strings. The total circuit depth is therefore $O(n \times N_{\text{params}}) = O(n \times N^4)$ for UCCSD, where $n$ is the number of qubits.

For H₂O with $n = 12$ active qubits and 200 amplitudes, this gives roughly $12 \times 200 = 2400$ CNOT gates — far beyond current NISQ coherence times. This is why hardware-efficient ansatzes or ADAPT-VQE (which selects only the most important excitations) are preferred for larger molecules.

19.3.4 Other Ansatz Designs

19.3.5 ADAPT-VQE: Growing the Ansatz Iteratively

ADAPT-VQE (Adaptive Derivative-Assembled Pseudo-Trotter) addresses the ansatz selection problem by building the circuit one operator at a time:

1. Start with the Hartree-Fock state (or an empty circuit).
2. Compute the gradient $\frac{\partial E}{\partial \theta_k}$ for every candidate operator $\hat{\tau}_k$ in the operator pool.
3. Append the operator with the largest gradient to the ansatz, with its parameter initialized to zero.
4. Re-optimize all parameters.
5. Repeat until the largest gradient falls below a threshold.

This procedure ensures that only operators that significantly affect the energy are included, producing compact ansatzes. The cost is the gradient computation at each iteration, but for small molecules this is feasible.

Try It Yourself: Implement a simplified ADAPT-VQE for H₂. Use a pool of single and double excitation operators. At each iteration, compute the gradient of each operator using the parameter-shift rule, add the largest-gradient operator to the ansatz, and re-optimize. How many operators does ADAPT-VQE select before convergence?

Worked Example 19.3c: ADAPT-VQE Gradient Computation

The gradient of the energy with respect to a new operator $\hat{\tau}_k - \hat{\tau}_k^\dagger$ appended to the current ansatz $U(\boldsymbol{\theta}^*)|\Phi_{\text{HF}}\rangle$ is:

$$\frac{\partial E}{\partial \theta_k}\bigg|_{\theta_k=0} = \langle\psi(\boldsymbol{\theta}^*)|[\hat{\tau}_k - \hat{\tau}_k^\dagger, H]|\psi(\boldsymbol{\theta}^*)\rangle$$

By the parameter-shift rule, this gradient can be computed exactly from two circuit evaluations:

$$\frac{\partial E}{\partial \theta_k}\bigg|_{\theta_k=0} = \frac{1}{2}\Big[E(\theta_k = +\pi/2) - E(\theta_k = -\pi/2)\Big]$$

where $E(\theta_k = \pm\pi/2)$ is the energy of the circuit with the operator $e^{\pm i\frac{\pi}{2}(\hat{\tau}_k - \hat{\tau}_k^\dagger)}$ appended.

For H₂ with the single double excitation operator $\hat{\tau} = a_2^\dagger a_3^\dagger a_1 a_0$, the gradient at the Hartree-Fock state is:

$$\frac{\partial E}{\partial \theta}\bigg|_{\theta=0} = \langle\Phi_{\text{HF}}|[\hat{\tau}^\dagger - \hat{\tau}, H]|\Phi_{\text{HF}}\rangle = 2\,\text{Im}\langle\Phi_{\text{HF}}|H\hat{\tau}|\Phi_{\text{HF}}\rangle$$

This is nonzero whenever the HF state is not an eigenstate of $H$, which is the case for any correlated molecule. ADAPT-VQE exploits the fact that the gradient magnitude indicates how important each operator is for lowering the energy.

The algorithm terminates when the largest gradient falls below a threshold $\epsilon_{\text{grad}}$:

$$\max_k \left|\frac{\partial E}{\partial \theta_k}\right| < \epsilon_{\text{grad}}$$

For chemical accuracy (1 mHartree), typical thresholds are $\epsilon_{\text{grad}} \sim 10^{-3}$ Hartree/bohr. ADAPT-VQE for H₂ typically selects only 1–3 operators before convergence, producing an extremely compact ansatz.

19.4 The Molecular Hamiltonian as Pauli Strings

To measure $\langle H \rangle$ on a quantum computer, we must express the molecular Hamiltonian in terms of Pauli operators. This involves three mappings.

19.4.1 Second Quantization

In second quantization, the electronic Hamiltonian is:

$$H = \sum_{pq} h_{pq} a_p^\dagger a_q + \frac{1}{2} \sum_{pqrs} h_{pqrs} a_p^\dagger a_q^\dagger a_r a_s$$

where $h_{pq}$ are one-electron integrals (kinetic energy + nuclear attraction) and $h_{pqrs}$ are two-electron integrals (electron-electron repulsion). The indices $p, q, r, s$ run over spin-orbitals.

Derivation: The One- and Two-Electron Integrals

The one-electron integrals are:

$$h_{pq} = \int \phi_p^*(\mathbf{r}) \left(-\frac{1}{2}\nabla^2 - \sum_A \frac{Z_A}{|\mathbf{r} - \mathbf{R}_A|}\right) \phi_q(\mathbf{r}) \, d\mathbf{r}$$

The two-electron integrals are:

$$h_{pqrs} = \int \int \frac{\phi_p^*(\mathbf{r}_1) \phi_q^*(\mathbf{r}_2) \phi_r(\mathbf{r}_1) \phi_s(\mathbf{r}_2)}{|\mathbf{r}_1 - \mathbf{r}_2|} \, d\mathbf{r}_1 \, d\mathbf{r}_2$$

For a molecule with $K$ spatial orbitals, there are $O(K^2)$ one-electron integrals and $O(K^4)$ two-electron integrals. The quartic scaling of the two-electron integrals is the fundamental bottleneck.

19.4.2 Fermion-to-Qubit Mapping

Fermionic operators obey anticommutation relations: $\{a_p, a_q^\dagger\} = \delta_{pq}$. Qubit operators (Pauli matrices) obey commutation relations. We need a mapping that preserves the algebraic structure.

Jordan-Wigner Transformation:

$$a_p^\dagger = \frac{1}{2}(X_p - iY_p) \otimes Z_{p-1} \otimes \cdots \otimes Z_0$$

$$a_p = \frac{1}{2}(X_p + iY_p) \otimes Z_{p-1} \otimes \cdots \otimes Z_0$$

The string of $Z$ operators enforces the correct anticommutation relations. The cost: a single fermionic operator maps to $O(n)$ Pauli operators, making the Jordan-Wigner transformation inefficient for operators acting on high-index qubits.

Worked Example 19.4: Jordan-Wigner for H₂

For H₂ in the STO-3G basis with 4 spin-orbitals, the relevant operators and their Jordan-Wigner mappings are:

$$a_0^\dagger = \frac{1}{2}(X_0 - iY_0) \otimes I_1 \otimes I_2 \otimes I_3$$

$$a_1^\dagger = \frac{1}{2}(X_1 - iY_1) \otimes Z_0 \otimes I_2 \otimes I_3$$

$$a_2^\dagger = \frac{1}{2}(X_2 - iY_2) \otimes Z_1 \otimes Z_0 \otimes I_3$$

$$a_3^\dagger = \frac{1}{2}(X_3 - iY_3) \otimes Z_2 \otimes Z_1 \otimes Z_0$$

Note how the parity strings grow with qubit index. The Hamiltonian after Jordan-Wigner mapping contains 15 non-trivial Pauli terms (after removing terms that vanish by symmetry).

Bravyi-Kitaev Transformation: A more efficient mapping that reduces the weight of Pauli strings from $O(n)$ to $O(\log n)$ by using a binary-tree encoding of occupation numbers. This is preferred for larger systems.

Parity Mapping: Encodes the parity of occupation numbers, with similar $O(\log n)$ scaling. The parity mapping also enables two-qubit reductions by exploiting $Z_2$ symmetries (particle number and spin parity).

19.4.3 The Final Form

After the mapping, the Hamiltonian becomes a weighted sum of Pauli strings:

$$H = \sum_{j=1}^{M} c_j P_j, \qquad P_j \in \{I, X, Y, Z\}^{\otimes n}$$

where $c_j \in \mathbb{R}$ are coefficients and each $P_j$ is a tensor product of $n$ Pauli matrices. For a molecule with $N$ spin-orbitals, $M$ scales as $O(N^4)$ — the number of terms grows quartically with system size.

Example: H₂ in the STO-3G basis (4 spin-orbitals, 2 electrons). After symmetry reductions (freezing core orbitals and removing qubits corresponding to fixed symmetries), the Hamiltonian in the Bravyi-Kitaev mapping with two-qubit reduction becomes:

$$H = c_0 I + c_1 Z_0 + c_2 Z_1 + c_3 Z_0 Z_1 + c_4 X_0 X_1 + c_5 Y_0 Y_1$$

with only 6 non-trivial terms. This is extremely compact — VQE for H₂ is essentially a 1-parameter optimization (since only one double excitation is relevant after symmetry reduction).

The expectation value is:

$$\langle H \rangle = \sum_{j=1}^{M} c_j \langle P_j \rangle$$

Each $\langle P_j \rangle$ is measured by rotating into the appropriate Pauli basis and sampling.

19.5 Measuring Expectation Values

19.5.1 Single Pauli Measurement

To measure $\langle Z \rangle$, we simply measure in the computational basis:

$$\langle Z \rangle = P(0) - P(1)$$

where $P(0)$ and $P(1)$ are the probabilities of measuring $|0\rangle$ and $|1\rangle$.

To measure $\langle X \rangle$, we apply a Hadamard gate before measurement (since $HZH = X$):

$$\langle X \rangle = \langle \psi | X | \psi \rangle = \langle \psi | H Z H | \psi \rangle$$

To measure $\langle Y \rangle$, we apply $S^\dagger$ then $H$ (since $HS^\dagger ZSH = Y$):

$$\langle Y \rangle = \langle \psi | Y | \psi \rangle = \langle \psi | HS^\dagger ZSH | \psi \rangle$$

For a multi-qubit Pauli string $P = P_1 \otimes P_2 \otimes \cdots \otimes P_n$, we apply the appropriate basis rotation to each qubit and measure all qubits simultaneously.

Worked Example 19.5: Measuring $\langle X \otimes Y \otimes Z \rangle$

For the Pauli string $X_0 Y_1 Z_2$ on three qubits:

q_0: ──[H]──── Measure → bit 0
q_1: ──[S†]──[H]── Measure → bit 1
q_2: ──────────── Measure → bit 2

⟨X₀Y₁Z₂⟩ = E[(-1)^(b₀⊕1) × (-1)^(b₁⊕1) × (-1)^b₂ × i^(-something)]
          = E[s₀ × s₁ × s₂]

where $s_i = (-1)^{b_i}$ for computational basis measurement, with appropriate sign changes from the basis rotations. After $S$ shots, the estimate is $\hat{O} = \frac{1}{S}\sum_{k=1}^S \text{sign}(P, b^{(k)})$.

19.5.2 Naive Measurement Cost

If we measure each of the $M$ Pauli terms independently, the total number of circuit executions is $M \times S$, where $S$ is the number of shots per term. For H₂O in a minimal basis, $M \sim 10^4$; for larger molecules, $M \sim 10^6$–$10^8$. This is prohibitive.

Derivation: Measurement Precision

To achieve precision $\epsilon$ on $\langle H \rangle$, the number of shots per term scales as:

$$S_j \geq \frac{c_j^2}{\epsilon^2}$$

for term $j$ with coefficient $c_j$. The total number of shots is:

$$N_{\text{shots}} = \sum_{j=1}^{M} S_j = \frac{1}{\epsilon^2} \sum_{j=1}^{M} c_j^2$$

For the Hamiltonian $H = \sum_j c_j P_j$, the quantity $\sum_j c_j^2 = \text{Tr}(H^2)/2^n$ (up to normalization). This sets the fundamental measurement cost of VQE.

19.5.3 Grouping Commuting Observables

The key optimization: Pauli terms that commute can be measured simultaneously. Two Pauli strings $P$ and $Q$ commute if and only if they commute qubit-wise: for each qubit $i$, the Pauli operators $P_i$ and $Q_i$ commute.

Qubit-wise commuting (QWC) groups: Terms that share a common eigenbasis. For example, $IZ$, $ZI$, and $ZZ$ all commute qubit-wise and can be measured in a single computational basis measurement.

General commuting (GC) groups: Terms that commute as operators but may not commute qubit-wise. For example, $XX$ and $YY$ commute ($[XX, YY] = 0$) but do not commute qubit-wise. They can be measured simultaneously using a shared Clifford transformation.

The grouping problem — partition $M$ Pauli terms into the minimum number of commuting groups — is NP-hard in general, but greedy heuristics work well in practice.

from qiskit.quantum_info import SparsePauliOp
import numpy as np

def greedy_commuting_groups(hamiltonian: SparsePauliOp):
    """
    Greedy grouping of Pauli terms into qubit-wise commuting sets.
    Each group can be measured with a single circuit execution.
    """
    paulis = hamiltonian.paulis
    coeffs = hamiltonian.coeffs
    remaining = list(range(len(paulis)))
    groups = []

    while remaining:
        group = [remaining[0]]
        for i in remaining[1:]:
            if all(paulis[group[0]].commutes(paulis[i]) for _ in [0]):
                # Check qubit-wise commutativity
                if paulis[group[0]].commutes(paulis[i]):
                    group.append(i)
        groups.append(group)
        remaining = [i for i in remaining if i not in group]

    return groups

Measurement cost reduction. If $M$ terms are partitioned into $G$ groups, the number of distinct circuits is reduced from $M$ to $G$. For typical molecules, $G \ll M$, often by a factor of 10–100.

Advanced Grouping: Diagonal Unitary Approach

A more sophisticated approach uses Clifford unitaries to rotate commuting groups into simultaneously measurable bases. Given a group of commuting Pauli operators $\{P_1, P_2, \ldots, P_k\}$, find a Clifford unitary $C$ such that $C P_1 C^\dagger, C P_2 C^\dagger, \ldots$ are all diagonal (i.e., tensor products of $I$ and $Z$). Then measure in the $C^\dagger$ basis.

Step-by-step procedure: 1. For each commuting group $\{P_1, \ldots, P_k\}$, find a Clifford unitary $C$ such that $C P_j C^\dagger \in \{I, Z\}^{\otimes n}$ for all $j$. 2. Measure the state $C|\psi\rangle$ in the computational basis. 3. Compute $\langle P_j \rangle$ by classical post-processing of the measurement results.

The number of distinct measurement settings equals the number of groups $G$, which can be as small as $O(\log M)$ for certain Hamiltonian structures.

Common Misconception: "Commuting groups always give $O(1)$ measurement settings." While commuting groups reduce the number of settings from $M$ to $G$, the group-mining problem is NP-hard in general. Practical heuristics (greedy, largest-first, sorted-insertion) achieve $G \sim O(M/\sqrt{n})$ for molecular Hamiltonians, which is a significant improvement but not constant.

Worked Example 19.6: Grouping the H₂ Hamiltonian

The reduced H₂ Hamiltonian has terms: $Z_0$, $Z_1$, $Z_0 Z_1$, $X_0 X_1$, $Y_0 Y_1$.

• Group 1 (Z-basis): $Z_0$, $Z_1$, $Z_0 Z_1$ — all measurable in the computational basis.
• Group 2 (X-basis): $X_0 X_1$ — requires Hadamard on both qubits.
• Group 3 (Y-basis): $Y_0 Y_1$ — requires $S^\dagger H$ on both qubits.

Three measurement settings instead of five. For larger molecules, the savings are much more dramatic.

19.6 The Classical Optimizer

The classical optimizer navigates the parameter landscape $E(\boldsymbol{\theta})$ to find the minimum. The choice of optimizer significantly affects VQE performance.

19.6.1 Gradient-Free Optimizers

COBYLA (Constrained Optimization By Linear Approximation): A simplex-based method that builds a linear model of the objective function. Robust to noise, no gradient required. Default in many VQE implementations.

Nelder-Mead: Another simplex method. Simple but can stagnate in high dimensions.

Powell's method: Performs line searches along conjugate directions. Effective for smooth, low-noise landscapes.

19.6.2 Gradient-Based Optimizers

When gradients are available (via the parameter-shift rule), gradient-based methods often converge faster.

The Parameter-Shift Rule: For a gate of the form $G(\theta) = e^{-i\theta P/2}$ where $P$ is a Pauli operator, the gradient is:

$$\frac{\partial E}{\partial \theta} = \frac{1}{2}\Big[E(\theta + \pi/2) - E(\theta - \pi/2)\Big]$$

Derivation of the Parameter-Shift Rule:

Consider the expectation value $E(\theta) = \langle\psi|U^\dagger(\theta) O U(\theta)|\psi\rangle$ where $U(\theta) = e^{-i\theta P/2}$.

Using the identity $e^{-i(\theta \pm \pi/2) P/2} = e^{-i\theta P/2} e^{\mp i\pi P/4}$ and the fact that $P^2 = I$ (for Pauli operators), we can expand:

$$E(\theta + \pi/2) = \langle\psi|U^\dagger(\theta) e^{i\pi P/4} O e^{-i\pi P/4} U(\theta)|\psi\rangle$$

$$E(\theta - \pi/2) = \langle\psi|U^\dagger(\theta) e^{-i\pi P/4} O e^{i\pi P/4} U(\theta)|\psi\rangle$$

Subtracting and using $e^{\pm i\pi P/4} = \frac{1}{\sqrt{2}}(I \pm iP)$:

$$E(\theta + \pi/2) - E(\theta - \pi/2) = 2\langle\psi|U^\dagger(\theta) \frac{iP}{\sqrt{2}} O \frac{1}{\sqrt{2}} U(\theta)|\psi\rangle + \text{c.c.}$$

This simplifies to $\frac{\partial E}{\partial \theta}$ as claimed. $\square$

This is exact — no finite-difference approximation. It requires two circuit evaluations per parameter. For $p$ parameters, $2p$ circuit executions per gradient evaluation.

SPSA (Simultaneous Perturbation Stochastic Approximation): Approximates the gradient using only two circuit evaluations regardless of $p$:

$$\frac{\partial E}{\partial \theta_k} \approx \frac{E(\boldsymbol{\theta} + c\boldsymbol{\Delta}) - E(\boldsymbol{\theta} - c\boldsymbol{\Delta})}{2c \Delta_k}$$

where $\boldsymbol{\Delta}$ is a random perturbation vector (typically $\Delta_k \in \{+1, -1\}$). SPSA is highly noise-robust and is the optimizer of choice for noisy VQE.

Recurring Theme: Noise Is the Enemy. On NISQ hardware, the energy evaluations are noisy. Gradient-based optimizers like L-BFGS-B are sensitive to this noise, while SPSA and COBYLA are more robust. The choice of optimizer is not merely a matter of convergence speed — it is intimately tied to the noise characteristics of the hardware.

19.6.3 Optimizer Comparison

19.6.4 Optimizer Selection Strategy

A practical strategy for VQE optimization:

1. Start with COBYLA for an initial rough optimization (robust to noise, requires no gradient).
2. Switch to SPSA if COBYLA stagnates (better scaling with parameter count).
3. Use parameter-shift gradients with Adam or L-BFGS-B for fine-tuning on a noiseless simulator (where gradient quality is high).
4. Always use multiple random restarts — VQE landscapes have many local minima.

19.6.5 Natural Gradient and Quantum Fisher Information

A powerful alternative to standard gradient descent is quantum natural gradient descent, which uses the quantum Fisher information (QFI) matrix to precondition the gradient. The update rule is:

$$\boldsymbol{\theta}_{k+1} = \boldsymbol{\theta}_k - \eta \, g^{-1}(\boldsymbol{\theta}_k) \nabla E(\boldsymbol{\theta}_k)$$

where $g(\boldsymbol{\theta})$ is the Fubini-Study metric tensor (quantum Fisher information matrix):

$$g_{ij}(\boldsymbol{\theta}) = \text{Re}\left[\frac{\partial \langle\psi|}{\partial \theta_i}\frac{\partial |\psi\rangle}{\partial \theta_j} - \frac{\partial \langle\psi|}{\partial \theta_i}|\psi\rangle\langle\psi|\frac{\partial |\psi\rangle}{\partial \theta_j}\right]$$

Step-by-step derivation of the metric tensor:

Consider the fidelity between two nearby states: $F(\boldsymbol{\theta}, \boldsymbol{\theta} + d\boldsymbol{\theta}) = |\langle\psi(\boldsymbol{\theta})|\psi(\boldsymbol{\theta} + d\boldsymbol{\theta})\rangle|^2$. For small $d\boldsymbol{\theta}$:

$$F = 1 - d\boldsymbol{\theta}^T g(\boldsymbol{\theta}) \, d\boldsymbol{\theta} + O(|d\boldsymbol{\theta}|^3)$$

The natural gradient $g^{-1} \nabla E$ accounts for the geometry of the parameter space, ensuring that equal steps in parameter space correspond to equal changes in the quantum state. This avoids the problem of "slow directions" in the optimization landscape.

For a parameterized circuit where $\theta_i$ parameterizes a gate $e^{-i\theta_i H_i/2}$ with $H_i^2 = I$ (a Pauli rotation), the metric tensor elements can be computed using the parameter-shift rule:

$$g_{ij}(\boldsymbol{\theta}) = \frac{1}{4}\Big[\langle\psi(\boldsymbol{\theta}^{i+})|\psi(\boldsymbol{\theta}^{j+})\rangle + \langle\psi(\boldsymbol{\theta}^{i-})|\psi(\boldsymbol{\theta}^{j-})\rangle - \langle\psi(\boldsymbol{\theta}^{i+})|\psi(\boldsymbol{\theta}^{j-})\rangle - \langle\psi(\boldsymbol{\theta}^{i-})|\psi(\boldsymbol{\theta}^{j+})\rangle\Big]$$

where $\boldsymbol{\theta}^{i\pm}$ denotes $\boldsymbol{\theta}$ with $\theta_i \to \theta_i \pm \pi/2$.

Advantages of natural gradient: - Invariant to reparameterization of the ansatz. - Avoids pathological convergence along "valleys" in the energy landscape. - For VQE with UCCSD, natural gradient descent can converge in $O(1)$ iterations for quadratic landscapes, compared to $O(\kappa)$ iterations for vanilla gradient descent (where $\kappa$ is the condition number of $g$).

Disadvantages: - Computing the full metric tensor requires $O(p^2)$ circuit evaluations per iteration, which is expensive for large parameter counts. - The metric tensor can be singular or near-singular, requiring regularization. - On noisy hardware, the metric tensor estimates are unreliable due to shot noise.

Try It Yourself: Implement quantum natural gradient descent for the H₂ VQE. Compare the convergence curve (energy vs. iteration) with standard gradient descent and COBYLA. How many fewer iterations does natural gradient require to reach chemical accuracy?

19.7 Complete VQE Implementation for H₂

We now implement a complete VQE calculation for the hydrogen molecule using Qiskit Nature. This is a real quantum chemistry calculation — the same workflow used in research papers.

import numpy as np
from qiskit import QuantumCircuit
from qiskit.circuit import Parameter
from qiskit.quantum_info import SparsePauliOp
from qiskit_aer import AerSimulator
from qiskit_aer.primitives import EstimatorV2 as Estimator
from qiskit_nature.second_q.drivers import PySCFDriver
from qiskit_nature.second_q.mappers import JordanWignerMapper, ParityMapper
from qiskit_nature.second_q.formats.molecule_info import MoleculeInfo
from qiskit_nature.second_q.transformers import ActiveSpaceTransformer
from qiskit_nature.second_q.circuit.library import HartreeFock, UCCSD
from qiskit_nature.second_q.algorithms import GroundStateEigensolver
from qiskit_nature.second_q.problems import ElectronicStructureProblem
from scipy.optimize import minimize

# ── Step 1: Define the molecule ──
# H₂ at bond distance 0.735 Å
molecule = MoleculeInfo(
    symbols=["H", "H"],
    coords=[(0.0, 0.0, 0.0), (0.0, 0.0, 0.735)],
    multiplicity=1,  # singlet
    charge=0,
)

# ── Step 2: Run electronic structure (PySCF) to get integrals ──
driver = PySCFDriver.from_molecule_info(molecule, basis="sto3g")
problem = driver.run()

# Extract the electronic Hamiltonian
hamiltonian = problem.hamiltonian
print(f"Hamiltonian has {len(hamiltonian)} terms")
print("First 5 terms:")
for label, coeff in list(hamiltonian.items())[:5]:
    print(f"  {label}: {coeff:.6f}")

# ── Step 3: Map fermions to qubits ──
mapper = JordanWignerMapper()
qubit_hamiltonian = mapper.map(hamiltonian)
print(f"\nQubit Hamiltonian: {len(qubit_hamiltonian)} Pauli terms")
print(f"Number of qubits: {qubit_hamiltonian.num_qubits}")

# ── Step 4: Build the ansatz ──
n_qubits = qubit_hamiltonian.num_qubits
n_electrons = problem.num_particles

# Hartree-Fock initial state
hf_state = HartreeFock(
    num_spatial_orbitals=n_qubits // 2,
    num_particles=n_electrons,
    qubit_mapper=mapper,
)

# UCCSD ansatz
uccsd = UCCSD(
    num_spatial_orbitals=n_qubits // 2,
    num_particles=n_electrons,
    qubit_mapper=mapper,
    initial_state=hf_state,
)

# Combine: HF initial state + UCCSD excitations
ansatz = QuantumCircuit(n_qubits)
ansatz.compose(hf_state, inplace=True)
ansatz.compose(uccsd, inplace=True)

print(f"\nAnsatz parameters: {ansatz.num_parameters}")
print(f"Ansatz depth: {ansatz.depth()}")
print(f"Ansatz 2-qubit gates: {ansatz.num_nonlocal_gates()}")

# ── Step 5: Define the cost function ──
def cost_function(params, ansatz, hamiltonian, estimator):
    """Evaluate E(θ) = ⟨ψ(θ)|H|ψ(θ)⟩."""
    bound_ansatz = ansatz.assign_parameters(params)
    job = estimator.run([(bound_ansatz, hamiltonian)])
    result = job.result()
    return result[0].data.evs  # expectation value

# ── Step 6: Run VQE ──
estimator = Estimator()
backend = AerSimulator()

# Initial parameters (small random values around zero)
np.random.seed(42)
initial_params = np.random.uniform(-0.1, 0.1, ansatz.num_parameters)

# COBYLA optimization
result = minimize(
    cost_function,
    initial_params,
    args=(ansatz, qubit_hamiltonian, estimator),
    method='COBYLA',
    options={'maxiter': 200, 'disp': True},
)

optimized_params = result.x
optimized_energy = result.fun

print(f"\n{'='*60}")
print(f"VQE Results for H₂ (STO-3G)")
print(f"{'='*60}")
print(f"Optimized energy: {optimized_energy:.8f} Hartree")
print(f"Exact energy (FCI): -1.13728383 Hartree")
print(f"Error: {abs(optimized_energy - (-1.13728383)):.8f} Hartree")
print(f"Error: {abs(optimized_energy - (-1.13728383)) * 627.509:.4f} kcal/mol")
print(f"Chemical accuracy (1 kcal/mol): "
      f"{'ACHIEVED' if abs(optimized_energy - (-1.13728383)) * 627.509 < 1 else 'NOT ACHIEVED'}")

19.7.1 Dissociation Curve

A key test of VQE is computing the dissociation curve — the ground state energy as a function of interatomic distance. This reveals bond lengths, bond energies, and the quality of the ansatz at strong correlation (large bond distances).

def compute_dissociation_curve(distances, molecule_base):
    """Compute VQE energies for a range of bond distances."""
    energies = []
    for d in distances:
        mol = MoleculeInfo(
            symbols=molecule_base.symbols,
            coords=[(0.0, 0.0, 0.0), (0.0, 0.0, d)],
            multiplicity=molecule_base.multiplicity,
            charge=molecule_base.charge,
        )
        driver = PySCFDriver.from_molecule_info(mol, basis="sto3g")
        problem = driver.run()
        hamiltonian = problem.hamiltonian
        qubit_ham = mapper.map(hamiltonian)

        # Rebuild ansatz for this geometry
        hf = HartreeFock(
            num_spatial_orbitals=qubit_ham.num_qubits // 2,
            num_particles=problem.num_particles,
            qubit_mapper=mapper,
        )
        ucc = UCCSD(
            num_spatial_orbitals=qubit_ham.num_qubits // 2,
            num_particles=problem.num_particles,
            qubit_mapper=mapper,
            initial_state=hf,
        )
        ansatz = QuantumCircuit(qubit_ham.num_qubits)
        ansatz.compose(hf, inplace=True)
        ansatz.compose(ucc, inplace=True)

        init = np.random.uniform(-0.1, 0.1, ansatz.num_parameters)
        res = minimize(
            cost_function, init,
            args=(ansatz, qubit_ham, estimator),
            method='COBYLA', options={'maxiter': 200}
        )
        energies.append(res.fun)
        print(f"  d = {d:.3f} Å, E = {res.fun:.6f} Hartree")

    return energies

# Compute dissociation curve
distances = np.linspace(0.4, 3.0, 14)
# energies = compute_dissociation_curve(distances, molecule)
# (Commented out — this takes several minutes to run)

Why the dissociation curve matters: At equilibrium, Hartree-Fock provides a good reference, and UCCSD captures the remaining correlation. At large distances (bond-breaking), the Hartree-Fock reference breaks down — the true wavefunction becomes multi-reference, with significant contributions from multiple Slater determinants. This is the regime where VQE can potentially outperform classical single-reference methods.

Dissociation curve (schematic):

E (Hartree)
  |
  |           *  *  *
  |         *        *
  |       *            *        HF fails here
  |      *               *     (wrong dissociation limit)
  |     *                  *
  |    *                     *
  |   *   VQE follows FCI   *  ←  HF dissociation limit
  |  *     (correct)          *   is too high
  | *                          *
  +------------------------------→ R (Å)
  0.4  0.7  1.0  1.5  2.0  2.5  3.0
       ↑
    equilibrium

19.8 Error Mitigation for VQE

On real hardware, VQE results are degraded by noise. Two error mitigation techniques are particularly important.

Recurring Theme: Noise Is the Enemy. Error mitigation is distinct from error correction. Mitigation reduces the effect of noise on expectation values without encoding logical qubits into physical ones. It is a NISQ-era strategy — effective for shallow circuits, but not a substitute for full quantum error correction.

19.8.1 Readout Error Mitigation

Measurement errors — reading $|0\rangle$ as $|1\rangle$ or vice versa — are among the largest error sources. The measurement error can be characterized by a confusion matrix $M$:

$$M_{ij} = P(\text{measure } i \mid \text{prepare } j)$$

For $n$ qubits, $M$ is a $2^n \times 2^n$ matrix. In practice, we assume tensor-product structure (uncorrelated readout errors):

$$M = M_1 \otimes M_2 \otimes \cdots \otimes M_n$$

where each $M_k$ is a $2 \times 2$ matrix for qubit $k$. The calibration procedure:

1. Prepare each computational basis state $|j\rangle$.
2. Measure and record the distribution of outcomes.
3. Build $M$ from the observed frequencies.
4. Given noisy measurement probabilities $\mathbf{p}_{\text{noisy}}$, estimate the true probabilities:

$$\mathbf{p}_{\text{true}} = M^{-1} \mathbf{p}_{\text{noisy}}$$

from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator

def calibrate_readout_error(backend, n_qubits, shots=10000):
    """Build the readout confusion matrix for n qubits."""
    M = np.zeros((2**n_qubits, 2**n_qubits))
    for j in range(2**n_qubits):
        qc = QuantumCircuit(n_qubits, n_qubits)
        # Prepare state |j⟩
        for i in range(n_qubits):
            if (j >> i) & 1:
                qc.x(i)
        qc.measure(range(n_qubits), range(n_qubits))

        sim = AerSimulator()
        counts = sim.run(qc, shots=shots).result().get_counts()
        for bitstring, count in counts.items():
            i = int(bitstring, 2)
            M[i, j] = count / shots
    return M

def mitigate_readout(counts, M):
    """Apply readout error mitigation to measurement counts."""
    n = int(np.log2(M.shape[0]))
    p_noisy = np.zeros(2**n)
    for bitstring, count in counts.items():
        p_noisy[int(bitstring, 2)] = count / sum(counts.values())

    p_true = np.linalg.solve(M, p_noisy)
    # Clip to [0, 1] and renormalize
    p_true = np.clip(p_true, 0, 1)
    p_true /= p_true.sum()

    mitigated_counts = {}
    for i in range(2**n):
        if p_true[i] > 1e-10:
            mitigated_counts[format(i, f'0{n}b')] = int(p_true[i] * sum(counts.values()))
    return mitigated_counts

19.8.2 Zero-Noise Extrapolation (ZNE)

ZNE exploits the fact that noise scales with circuit depth (or gate duration). By running the circuit at multiple noise levels and extrapolating to zero noise, we can estimate the ideal expectation value.

Noise scaling methods: - Gate stretching: Increase the duration of gates (only possible with pulse-level control). - Unitary folding: Insert identity-equivalent gate sequences ($G G^\dagger$) to effectively increase depth. - Local folding: Fold individual gates: $G \to G (G^\dagger G)^n$.

Extrapolation models: - Linear: $E(\lambda) = E_0 + a\lambda$ - Polynomial: $E(\lambda) = E_0 + a\lambda + b\lambda^2$ - Exponential: $E(\lambda) = E_0 e^{-b\lambda} + E_\infty(1 - e^{-b\lambda})$

Derivation: Linear ZNE with Unitary Folding

Consider a circuit $U = U_m \cdots U_2 U_1$ with noise channel $\mathcal{E}$ after each gate. The noisy expectation value at noise scale factor $\lambda = 1$ is:

$$E(1) = \langle O \rangle_{\text{noisy}}$$

By folding the circuit $U \to U(U^\dagger U)^k$, we increase the effective noise by a factor of $\lambda = 2k + 1$. The expectation values at $\lambda = 1, 3, 5$ give:

$$E(1) = E_0 + a, \quad E(3) = E_0 + 3a, \quad E(5) = E_0 + 5a$$

Linear extrapolation: $E_0 = 3E(1) - \frac{1}{2}[E(3) + E(5)]$... or more precisely:

$$E_0 = \frac{\lambda_2 E(\lambda_1) - \lambda_1 E(\lambda_2)}{\lambda_2 - \lambda_1}$$

for two noise scale factors $\lambda_1$ and $\lambda_2$.

def zero_noise_extrapolation(ansatz, hamiltonian, estimator,
                              noise_levels, backend):
    """
    Perform ZNE by unitary folding at multiple noise levels.
    noise_levels: list of scale factors (1.0 = original circuit)
    """
    energies = []
    for scale in noise_levels:
        # Fold the circuit: each gate G → G (G† G)^(scale-1)
        folded_circuit = fold_circuit(ansatz, scale)
        job = estimator.run([(folded_circuit, hamiltonian)])
        result = job.result()
        energies.append(result[0].data.evs)

    # Extrapolate to zero noise
    noise_levels = np.array(noise_levels)
    energies = np.array(energies)

    # Linear extrapolation
    coeffs = np.polyfit(noise_levels, energies, 1)
    e_zero_noise = coeffs[1]  # intercept

    return e_zero_noise, noise_levels, energies

def fold_circuit(circuit, scale):
    """Naive global folding: U → U (U† U)^(scale-1)."""
    folded = circuit.copy()
    for _ in range(int(scale) - 1):
        folded.compose(circuit.inverse(), inplace=True)
        folded.compose(circuit, inplace=True)
    return folded

Common Misconception: "Error mitigation solves the noise problem." Error mitigation can significantly reduce systematic errors, but it does not eliminate all noise. ZNE assumes a particular noise model (typically depolarizing), and extrapolation errors can be large if the model is wrong. Readout mitigation requires calibrating the confusion matrix, which drifts over time. Error mitigation is essential for NISQ, but it is not a substitute for error correction.

19.8.3 Probabilistic Error Cancellation (PEC)

Beyond ZNE and readout mitigation, probabilistic error cancellation (PEC) provides a systematic method to remove bias from noisy expectation values. PEC works by decomposing the noisy channel into a linear combination of inverse noise channels.

Mathematical formulation:

Let $\mathcal{E}$ be the actual (noisy) channel acting on the quantum state. If we can express the ideal (noiseless) channel $\mathcal{I}$ as:

$$\mathcal{I} = \sum_i \alpha_i \mathcal{G}_i$$

where $\{\mathcal{G}_i\}$ are quasi-probability operations (implementable on the quantum computer) and $\{\alpha_i\}$ are real coefficients (possibly negative), then we can estimate $\langle O \rangle_{\text{ideal}}$ by:

1. Sampling operations $\mathcal{G}_i$ with probability proportional to $|\alpha_i|$.
2. Weighting each measurement outcome by the sign-corrected factor $\text{sign}(\alpha_i) \gamma^{-1}$, where $\gamma = \sum_i |\alpha_i|$ is the sampling overhead.

The estimate is:

$$\langle O \rangle_{\text{ideal}} = \gamma^{-1} \sum_i \text{sign}(\alpha_i) \langle O \rangle_{\mathcal{G}_i}$$

Step-by-step PEC example for a single depolarizing gate:

Consider a single-qubit gate followed by a depolarizing channel with parameter $p$. The ideal channel $\mathcal{I}$ can be expressed as:

$$\mathcal{I} = \frac{1}{1-p}\mathcal{E}_{\text{dep}} - \frac{p}{3(1-p)}(X\mathcal{E}_{\text{dep}}X + Y\mathcal{E}_{\text{dep}}Y + Z\mathcal{E}_{\text{dep}}Z)$$

The coefficients are $\alpha_0 = \frac{1}{1-p}$, $\alpha_X = \alpha_Y = \alpha_Z = -\frac{p}{3(1-p)}$, and $\gamma = \frac{1}{1-p} + 3 \cdot \frac{p}{3(1-p)} = \frac{1+p}{1-p}$.

To estimate the ideal expectation value, we: 1. Run the circuit as-is (with probability $\propto |\alpha_0| = 1/(1-p)$). 2. Run the circuit with an additional $X$, $Y$, or $Z$ gate before the noisy gate (each with probability $\propto p/(3(1-p))$). 3. Combine the results with appropriate sign weights.

The sampling overhead scales as $\gamma^d$ where $d$ is the circuit depth, which grows exponentially for deep circuits. PEC is therefore practical only for shallow circuits, but it can provide exact (unbiased) error cancellation.

Common Misconception: "PEC eliminates all noise perfectly." While PEC provides an unbiased estimator, it comes at an exponential sampling cost. The overhead $\gamma^d$ means that for circuits with more than ~20 noisy gates, PEC requires more samples than are practically feasible. PEC and ZNE are complementary: PEC is exact but expensive, while ZNE is approximate but scalable.

19.9 Scaling to Larger Molecules: LiH and H₂O

19.9.1 Lithium Hydride (LiH)

LiH has 4 electrons and requires 12 spin-orbitals in the STO-3G basis (12 qubits). The Hamiltonian contains approximately 630 Pauli terms. Key challenges:

• Qubit count: 12 qubits is within current hardware capabilities but pushes the limits of classical simulation for verification.
• Ansatz depth: UCCSD for LiH requires hundreds of CNOT gates — challenging for NISQ fidelities.
• Optimization landscape: More parameters (dozens for UCCSD) mean a higher-dimensional optimization problem.

Symmetry reduction is essential: the number of qubits can be reduced by exploiting $Z_2$ symmetries (particle number, spin parity). With the ParityMapper and two-qubit reduction, LiH can be mapped to 8 qubits.

Worked Example 19.7: LiH Ground State by VQE

The LiH molecule (4 electrons, 6 spin-orbitals in minimal basis) has the following Hamiltonian structure after symmetry reduction:

Using the ParityMapper with two-qubit reduction, LiH maps to 8 qubits (from the original 12). The number of Pauli terms is approximately 630, which can be grouped into approximately 45 commuting groups for simultaneous measurement.

The UCCSD ansatz for LiH has approximately 20 parameters (after removing symmetry-forbidden excitations). The circuit depth for the UCCSD ansatz is approximately 1500 CNOT gates — challenging but feasible on near-term hardware with error mitigation.

Key steps for a LiH VQE calculation:

1. Define the molecule at various bond distances: $R \in [1.0, 4.0]$ Bohr.
2. At each $R$, compute the molecular integrals using PySCF.
3. Apply the ParityMapper with two-qubit reduction.
4. Group the Pauli terms into commuting sets.
5. Run VQE with COBYLA optimizer, measuring each commuting group separately.
6. Apply readout error mitigation and ZNE to the results.

At equilibrium ($R \approx 3.0$ Bohr), the Hartree-Fock energy is approximately $-7.863$ Hartree and the FCI energy is approximately $-7.882$ Hartree, giving a correlation energy of about 19 mHartree. VQE with UCCSD should recover 95%+ of this correlation energy.

The bond-breaking challenge: At large bond distances ($R > 4.0$ Bohr), the LiH molecule dissociates into Li + H. The ground state becomes strongly multi-reference: the wavefunction has significant contributions from multiple Slater determinants. Single-reference methods like CCSD break down (the CCSD energy diverges), while UCCSD with a Hartree-Fock reference also struggles. This is a regime where ADAPT-VQE or multi-reference methods are essential.

19.9.2 Water (H₂O)

H₂O has 10 electrons and requires 14 spin-orbitals in a minimal basis. After symmetry reduction, approximately 10 qubits. The Hamiltonian contains ~10,000 Pauli terms. This is at the edge of what NISQ VQE can handle.

Active space approximation: Freeze core orbitals (oxygen 1s) and discard high-energy virtual orbitals, reducing the active space to, say, 6 electrons in 8 orbitals (8 qubits). This is the standard approach for making larger molecules tractable.

from qiskit_nature.second_q.transformers import ActiveSpaceTransformer

# Define active space: 6 electrons in 8 spin-orbitals
transformer = ActiveSpaceTransformer(
    num_electrons=6,
    num_spatial_orbitals=4,  # 4 spatial = 8 spin-orbitals
)

# Apply to the problem
reduced_problem = transformer.transform(problem)
print(f"Reduced qubits: {reduced_problem.num_spin_orbitals}")

19.9.3 The Scaling Challenge

Molecule    Spin-Orbitals  Qubits (JW)  Qubits (reduced)  Pauli Terms  UCCSD Parameters
H₂          4               4            2                  15           1
LiH         12              12            8                  ~630         ~20
H₂O         14              14           10                  ~10,000      ~50
NH₃         16              16           12                  ~30,000      ~100
FeMoco*     108            108           ~54                 ~10⁶         ~10⁴

* FeMoco is the nitrogenase active site, a long-term target for quantum chemistry

Recurring Theme: We're at the Beginning. VQE has been demonstrated for small molecules (H₂, LiH, BeH₂) on current hardware. Scaling to chemically relevant molecules like FeMoco or complex catalysts requires hardware with hundreds to thousands of high-fidelity qubits — likely beyond the NISQ era. The current state is analogous to the ENIAC era of classical computing: the principle is proven, but practical impact awaits hardware maturation.

19.10 Challenges and Limitations

19.10.1 Barren Plateaus

For many ansatz designs, the gradient variance vanishes exponentially with the number of qubits:

$$\text{Var}\left[\frac{\partial E}{\partial \theta_k}\right] \sim O\left(\frac{1}{2^n}\right)$$

This means that for $n \gtrsim 20$, the optimization landscape is essentially flat everywhere — the optimizer cannot find a direction to descend, and training requires exponentially many measurements to resolve the gradient from statistical noise.

Detailed Analysis of Barren Plateaus

The barren plateau phenomenon was first rigorously analyzed by McClean et al. (2018). They showed that for any ansatz that forms a 2-design (a set of unitaries whose first and second moments match those of the Haar measure), the gradient variance is:

$$\text{Var}[\partial_k E] = \frac{\text{Tr}(H^2) - \frac{(\text{Tr}H)^2}{2^n}}{2^n(2^n + 1)}$$

This decays as $O(1/2^n)$, making gradient-based optimization infeasible for large systems.

Mitigations: - Use problem-inspired ansatzes (UCCSD) rather than hardware-efficient ones — they are not 2-designs. - Employ layerwise training strategies — train one layer at a time, freezing previous layers. - Use correlated parameter initialization — start near the Hartree-Fock state. - Employ local cost functions — measuring only a few qubits rather than all. - Initialize with classically computed parameters (e.g., from CCSD amplitudes).

Common Misconception: "Barren plateaus make VQE impossible for large systems." Barren plateaus are a serious concern, but they are not an absolute barrier. Problem-inspired ansatzes, clever initialization, and local cost functions can mitigate them. However, the scaling challenge remains real — VQE for large molecules is an active research frontier.

19.10.2 Measurement Overhead

Even with commuting group optimization, the number of measurements required for chemical accuracy ($\sim$1 kcal/mol $\approx$ 1.6 mHartree) grows rapidly with system size. For a molecule requiring $G$ groups and $S$ shots per group, the total shot count is $G \times S$. Achieving 1 mHartree precision may require $10^8$–$10^{10}$ shots for medium-sized molecules.

Derivation: Scaling of Measurement Overhead

The variance of the energy estimator is:

$$\text{Var}[\hat{E}] = \sum_{j=1}^{M} \frac{c_j^2 \text{Var}[\langle P_j \rangle]}{S_j}$$

where $S_j$ is the number of shots allocated to measuring Pauli term $j$. Minimizing the total number of shots subject to a precision constraint $\text{Var}[\hat{E}] \leq \epsilon^2$ gives the optimal allocation:

$$S_j = \frac{|c_j|}{\epsilon^2} \left(\sum_k |c_k|\right)$$

The total number of shots is $N = \frac{1}{\epsilon^2}\left(\sum_k |c_k|\right)^2$, which scales as $O(M^2/\epsilon^2)$ for typical molecular Hamiltonians.

19.10.3 Classical Competitors

Classical quantum chemistry methods — particularly coupled cluster with perturbative triples (CCSD(T)) and density matrix renormalization group (DMRG) — achieve high accuracy for many molecules. VQE must demonstrate advantage on problems where these methods fail, such as strongly correlated systems and transition metal complexes.

Recurring Theme: Quantum Advantage Is Problem-Specific. VQE does not need to outperform classical methods on all molecules. It needs to excel on specific, hard problems — strongly correlated systems, bond-breaking processes, and transition metal chemistry — where classical methods struggle. The path to advantage is through carefully chosen problem instances, not general superiority.

19.11 Second Qiskit Example: VQE with Noise Modeling and Error Mitigation

import numpy as np
from qiskit import QuantumCircuit
from qiskit.circuit import Parameter
from qiskit.quantum_info import SparsePauliOp
from qiskit_aer import AerSimulator
from qiskit_aer.noise import NoiseModel, depolarizing_error, thermal_relaxation_error
from qiskit_aer.primitives import EstimatorV2 as Estimator
from scipy.optimize import minimize

# ── Define a simple 2-qubit Hamiltonian (reduced H₂) ──
H2_hamiltonian = SparsePauliOp.from_list([
    ("II", -1.052),
    ("IZ",  0.397),
    ("ZI", -0.397),
    ("ZZ", -0.011),
    ("XX",  0.180),
])

# ── Build a hardware-efficient ansatz ──
theta = ParameterVector('θ', 4)
ansatz = QuantumCircuit(2)
ansatz.ry(theta[0], 0)
ansatz.ry(theta[1], 1)
ansatz.cx(0, 1)
ansatz.ry(theta[2], 0)
ansatz.ry(theta[3], 1)
ansatz.cx(0, 1)

# ── Define noise model ──
noise_model = NoiseModel()
p1q = 0.001  # single-qubit gate error
p2q = 0.01   # two-qubit gate error
noise_model.add_all_qubit_quantum_error(
    depolarizing_error(p1q, 1), ['ry', 'rz', 'h']
)
noise_model.add_all_qubit_quantum_error(
    depolarizing_error(p2q, 2), ['cx']
)

# ── Run VQE with and without noise ──
estimator_ideal = Estimator()
estimator_noisy = Estimator(
    backend_options=dict(noise_model=noise_model)
)

def cost_function(params, ansatz, hamiltonian, estimator):
    bound = ansatz.assign_parameters(dict(zip(theta, params)))
    job = estimator.run([(bound, hamiltonian)])
    return job.result()[0].data.evs

# Ideal optimization
np.random.seed(42)
x0 = np.random.uniform(-0.1, 0.1, 4)

result_ideal = minimize(
    cost_function, x0,
    args=(ansatz, H2_hamiltonian, estimator_ideal),
    method='COBYLA', options={'maxiter': 200}
)

result_noisy = minimize(
    cost_function, x0,
    args=(ansatz, H2_hamiltonian, estimator_noisy),
    method='COBYLA', options={'maxiter': 200}
)

print(f"Ideal VQE energy:  {result_ideal.fun:.8f}")
print(f"Noisy VQE energy:  {result_noisy.fun:.8f}")
print(f"Exact energy:      -1.13728383 Hartree")
print(f"Ideal error:        {abs(result_ideal.fun - (-1.13728)):.6f}")
print(f"Noisy error:        {abs(result_noisy.fun - (-1.13728)):.6f}")