Hybrid Quantum-Classical Optimization: A Deep Dive into VQE
A comprehensive guide to the Variational Quantum Eigensolver: ansatz design, optimizer selection, barren plateaus, and a complete H2 ground state calculation using Qiskit's Estimator primitive.
Circuit diagrams
The Core Idea Behind VQE
The Variational Quantum Eigensolver is built on the variational principle from quantum mechanics: for any Hamiltonian H and any normalized quantum state |psi>, the expectation value is always at least as large as the ground state energy E_0:
<psi|H|psi> >= E_0
This is an upper bound. The closer your trial state |psi> is to the true ground state, the tighter the bound. VQE exploits this by parameterizing a family of quantum states (the ansatz) and using a classical optimizer to minimize the expectation value <psi(theta)|H|psi(theta)> over the parameters theta.
The division of labor is deliberate: the quantum computer prepares and measures the state (a task exponentially hard classically), while the classical computer updates the parameters (a task that scales polynomially in the number of parameters). This hybrid loop is the defining feature of near-term quantum algorithms.
The VQE Loop in Detail
- Initialize parameters theta at some starting point (random or informed by chemistry)
- Prepare the parametric circuit |psi(theta)> on the quantum device
- Measure the expectation value E(theta) =
<psi(theta)|H|psi(theta)>by measuring each Pauli term in H - Pass E(theta) to the classical optimizer, which proposes a new theta
- Repeat steps 2-4 until convergence
The bottleneck is step 3: for a Hamiltonian with M Pauli terms, you need M separate circuit executions (or fewer if terms commute and can be co-diagonalized). For large molecular Hamiltonians, M can reach thousands, making each iteration expensive even on a simulator.
Ansatz Design: The Most Critical Choice
The ansatz determines what states are reachable. Too restrictive: VQE cannot find the ground state. Too expressive: the optimization landscape becomes unnavigable (barren plateau problem). Good ansatz design is part art, part chemistry.
Hardware-Efficient Ansatz
The hardware-efficient ansatz (HEA) uses single-qubit rotations and entangling gates arranged in layers. It does not respect any molecular symmetry, which makes it flexible but also means it can explore unphysical regions of Hilbert space.
from qiskit.circuit import QuantumCircuit, ParameterVector
def hardware_efficient_ansatz(n_qubits: int, reps: int) -> QuantumCircuit:
"""
Alternating layers of Ry rotations and CX entangling gates.
Total parameters: n_qubits * reps
"""
n_params = n_qubits * reps
theta = ParameterVector("theta", n_params)
qc = QuantumCircuit(n_qubits)
param_idx = 0
for rep in range(reps):
# Rotation layer
for q in range(n_qubits):
qc.ry(theta[param_idx], q)
param_idx += 1
# Entanglement layer (linear)
for q in range(n_qubits - 1):
qc.cx(q, q + 1)
return qc
hea = hardware_efficient_ansatz(4, reps=2)
print(hea.draw("text"))
print(f"Parameters: {hea.num_parameters}")
UCCSD Ansatz
Unitary Coupled Cluster with Singles and Doubles (UCCSD) is motivated by chemistry. It applies single excitations (one electron hops from an occupied to a virtual orbital) and double excitations (two electrons hop simultaneously). The resulting operator is anti-Hermitian, so its exponentiation gives a unitary:
|psi_UCCSD> = exp(T - T†) |HF>
where T is the cluster operator and |HF> is the Hartree-Fock reference state.
from qiskit_nature.second_q.circuit.library import UCCSD, HartreeFock
from qiskit_nature.second_q.mappers import JordanWignerMapper
# For a 2-orbital, 2-electron system (H2 STO-3G)
n_spatial_orbitals = 2
n_particles = (1, 1) # (alpha electrons, beta electrons)
mapper = JordanWignerMapper()
hf_state = HartreeFock(
num_spatial_orbitals=n_spatial_orbitals,
num_particles=n_particles,
qubit_mapper=mapper
)
uccsd_ansatz = UCCSD(
num_spatial_orbitals=n_spatial_orbitals,
num_particles=n_particles,
qubit_mapper=mapper,
initial_state=hf_state
)
print(f"UCCSD circuit depth: {uccsd_ansatz.decompose().depth()}")
print(f"UCCSD parameters: {uccsd_ansatz.num_parameters}")
UCCSD is chemistry-motivated and gives compact circuits for small molecules. For H2 in STO-3G it needs only 2 parameters. The downside: building the UCCSD circuit requires knowing the molecular structure ahead of time, and it becomes expensive for large active spaces.
Optimizer Comparison
Choosing the right optimizer is as important as choosing the right ansatz. The key question: does your hardware allow gradient computation?
COBYLA (Gradient-Free)
Constrained Optimization By Linear Approximation (COBYLA) constructs a local linear model of the objective and steps toward its minimum. It requires only function evaluations, no gradients. Robust to noise, but converges slowly for many parameters.
from qiskit_algorithms.optimizers import COBYLA
cobyla = COBYLA(maxiter=500, rhobeg=1.0)
# rhobeg controls the initial trust region size
Best for: fewer than ~20 parameters, noisy hardware, or when gradients are unreliable.
SPSA (Stochastic Gradient)
Simultaneous Perturbation Stochastic Approximation (SPSA) estimates the gradient using only 2 circuit evaluations regardless of the number of parameters. It perturbs all parameters simultaneously in a random direction:
gradient_approx = (f(theta + c*delta) - f(theta - c*delta)) / (2*c) * delta^{-1}
where delta is a random +/-1 vector. This is vastly more efficient than finite differences, which would require 2*N evaluations for N parameters.
from qiskit_algorithms.optimizers import SPSA
spsa = SPSA(
maxiter=300,
learning_rate=0.1,
perturbation=0.1,
second_order=False # set True for 2-SPSA (approximates Hessian)
)
Best for: many parameters (50+), hardware with shot noise, situations where gradient estimation is the bottleneck.
L-BFGS-B (Gradient-Based)
Limited-memory Broyden-Fletcher-Goldfarb-Shanno with Bounded variables uses gradient information to build an approximate inverse Hessian. Converges faster than gradient-free methods but requires computing exact gradients (parameter shift rule) and does not handle shot noise gracefully.
from qiskit_algorithms.optimizers import L_BFGS_B
lbfgsb = L_BFGS_B(maxiter=200, epsilon=1e-8)
# Use with statevector simulator where exact gradients are available
Best for: simulators with exact expectation values, fewer than ~50 parameters, when high accuracy is needed.
The Barren Plateau Problem
The most serious obstacle to scaling VQE is the barren plateau: as the number of qubits grows, the gradient of the cost function vanishes exponentially in the number of qubits. In flat regions, the optimizer cannot determine which direction to move.
Mathematically, for a random circuit on n qubits, the variance of the gradient decays as:
Var[d E / d theta_k] ~ O(2^{-n})
This means that with n=50 qubits, you need exponentially many measurements just to estimate the gradient direction reliably. Gradient-free methods like COBYLA suffer equally because the loss landscape is flat everywhere.
Mitigation strategies include:
- Local cost functions: Measure only a subset of qubits; local cost functions avoid barren plateaus but may have different minima
- Layer-by-layer training: Fix all layers except one, optimize that layer, then freeze it and move to the next
- Informed initialization: Start near the Hartree-Fock state rather than randomly; chemistry-motivated initializations avoid the flat regions
- Problem-specific ansatze: ADAPT-VQE grows the ansatz one operator at a time, avoiding large circuits that cause barren plateaus
Full H2 VQE Example
import numpy as np
from qiskit.quantum_info import SparsePauliOp
from qiskit.primitives import StatevectorEstimator
from qiskit_algorithms import VQE
from qiskit_algorithms.optimizers import COBYLA, SPSA
from qiskit.circuit import QuantumCircuit, ParameterVector
# H2 Hamiltonian at equilibrium geometry (STO-3G, JW mapping)
h2_hamiltonian = SparsePauliOp.from_list([
("IIII", -0.8105),
("IIIZ", +0.1722),
("IIZI", -0.2228),
("IIZZ", +0.1209),
("IZII", -0.2228),
("ZIII", +0.1722),
("ZZII", +0.1209),
("IZIZ", +0.1744),
("XXYY", -0.0453),
("YYXX", -0.0453),
("XYYX", +0.0453),
("YXXY", +0.0453),
])
# Chemistry-motivated 2-parameter ansatz
theta = ParameterVector("theta", 2)
ansatz = QuantumCircuit(4)
ansatz.x(0) # Hartree-Fock reference |0011> (qiskit ordering: qubit 0 = rightmost)
ansatz.x(1)
ansatz.ry(theta[0], 0)
ansatz.ry(theta[1], 2)
ansatz.cx(0, 1)
ansatz.cx(2, 3)
ansatz.cx(1, 2)
# Track convergence
energy_history = []
def callback(eval_count, params, value, metadata):
energy_history.append(value)
if eval_count % 50 == 0:
print(f" Eval {eval_count:4d}: E = {value:.6f} Hartree")
estimator = StatevectorEstimator()
optimizer = COBYLA(maxiter=400, rhobeg=0.5)
vqe = VQE(
estimator=estimator,
ansatz=ansatz,
optimizer=optimizer,
initial_point=[0.01, 0.01],
callback=callback,
)
print("Running VQE for H2...")
result = vqe.compute_minimum_eigenvalue(h2_hamiltonian)
print(f"\nFinal VQE energy: {result.eigenvalue:.6f} Hartree")
print(f"Optimal theta: {[f'{p:.4f}' for p in result.optimal_point]}")
Text-Based Convergence Plot
def text_convergence_plot(energies, width=60, height=15):
"""Print a crude ASCII convergence plot to the terminal."""
if not energies:
return
min_e = min(energies)
max_e = max(energies)
e_range = max_e - min_e if max_e > min_e else 1.0
print(f"\nConvergence ({len(energies)} evaluations)")
print(f"Max: {max_e:.4f} Min: {min_e:.4f} Hartree")
print("-" * (width + 10))
for row in range(height, -1, -1):
threshold = min_e + (row / height) * e_range
line = ""
step = max(1, len(energies) // width)
for i in range(0, len(energies), step):
line += "*" if energies[i] >= threshold else " "
label = f"{threshold:8.4f} |"
print(label + line)
print(" " * 9 + "-" * width)
print(f"{'Start':>14}" + " " * (width - 14) + "End")
text_convergence_plot(energy_history)
A well-behaved VQE run shows rapid initial descent followed by slow refinement. If the curve is noisy and never settles, the optimizer is struggling with shot noise or a flat landscape; try increasing shots, switching to SPSA, or reducing the ansatz complexity.
Comparing Optimizers Side by Side
from qiskit_algorithms.optimizers import L_BFGS_B
results = {}
for opt_name, optimizer in [
("COBYLA", COBYLA(maxiter=300)),
("SPSA", SPSA(maxiter=300, learning_rate=0.1, perturbation=0.1)),
("L-BFGS-B", L_BFGS_B(maxiter=200)),
]:
vqe_run = VQE(
estimator=StatevectorEstimator(),
ansatz=ansatz,
optimizer=optimizer,
initial_point=[0.01, 0.01],
)
res = vqe_run.compute_minimum_eigenvalue(h2_hamiltonian)
results[opt_name] = (res.eigenvalue, res.cost_function_evals)
print(f"{opt_name:12s}: E = {res.eigenvalue:.6f} evals = {res.cost_function_evals}")
exact_energy = -1.857275 # FCI/exact diagonalization
print(f"\nExact FCI energy: {exact_energy:.6f} Hartree")
for name, (e, evals) in results.items():
print(f"{name:12s}: error = {abs(e - exact_energy):.2e} Hartree")
Typical results: all three optimizers find the correct energy for H2 because the landscape is simple. On larger problems, COBYLA struggles with many parameters, SPSA handles shot noise best, and L-BFGS-B is fastest on simulators but degrades badly on real hardware.
Practical Recommendations
For research-scale VQE on NISQ hardware:
- Use UCCSD for small molecules where accuracy matters; use HEA for hardware benchmarking
- SPSA is the most hardware-appropriate optimizer: only 2 circuit evaluations per gradient step, robust to shot noise
- Increase shots to 8192 or more for accurate gradient estimates; fewer shots cause the optimizer to chase noise
- Apply zero-noise extrapolation (ZNE) from Mitiq to suppress coherent errors before running VQE on real devices
- For molecules beyond 10 qubits, consider ADAPT-VQE or orbital optimization to keep the circuit shallow
VQE remains the workhorse near-term quantum chemistry algorithm, but reaching chemical accuracy on classically intractable molecules will require fault-tolerant hardware or significantly better error mitigation strategies.
Was this tutorial helpful?