Variational Quantum Eigensolver in Python — Deep Dive

Implement VQE from scratch in Python — molecular Hamiltonians, ansatz design, gradient strategies, and noise-resilient optimization

Building a VQE Pipeline

A complete VQE implementation requires four components: molecular Hamiltonian construction, qubit mapping, ansatz circuit design, and a classical-quantum optimization loop.

Molecular Hamiltonian with Qiskit Nature

from qiskit_nature.second_q.drivers import PySCFDriver
from qiskit_nature.second_q.mappers import JordanWignerMapper
from qiskit_nature.second_q.problems import ElectronicStructureProblem

# Define the molecule (H2 at 0.74 Angstrom bond length)
driver = PySCFDriver(
    atom="H 0.0 0.0 0.0; H 0.0 0.0 0.74",
    basis="sto-3g",
    charge=0,
    spin=0
)

# Compute the electronic structure problem
problem = driver.run()

# Get the second-quantized Hamiltonian
hamiltonian = problem.hamiltonian.second_q_op()

# Map to qubit operators using Jordan-Wigner transformation
mapper = JordanWignerMapper()
qubit_op = mapper.map(hamiltonian)

print(f"Number of qubits: {qubit_op.num_qubits}")
print(f"Number of Pauli terms: {len(qubit_op)}")
# H2 in STO-3G: 4 qubits, ~15 Pauli terms

Qubit Mapping Strategies

The choice of fermion-to-qubit mapping affects circuit depth and qubit count:

from qiskit_nature.second_q.mappers import (
    JordanWignerMapper,
    ParityMapper,
    BravyiKitaevMapper
)

# Jordan-Wigner: simplest, each orbital → one qubit
jw_op = JordanWignerMapper().map(hamiltonian)

# Parity: enables qubit reduction for symmetries
parity_mapper = ParityMapper(num_particles=problem.num_particles)
parity_op = parity_mapper.map(hamiltonian)

# Bravyi-Kitaev: balances locality of different operators
bk_op = BravyiKitaevMapper().map(hamiltonian)

print(f"JW qubits: {jw_op.num_qubits}")
print(f"Parity qubits: {parity_op.num_qubits}")  # Often 2 fewer
print(f"BK qubits: {bk_op.num_qubits}")

Parity mapping with two-qubit reduction is standard for small molecules — it exploits particle number and spin symmetry to eliminate qubits.

Ansatz Design

UCCSD Ansatz

The chemically-motivated Unitary Coupled Cluster ansatz:

from qiskit_nature.second_q.circuit.library import UCCSD, HartreeFock

# Reference state: Hartree-Fock configuration
num_spatial_orbitals = problem.num_spatial_orbitals
num_particles = problem.num_particles

hf_state = HartreeFock(
    num_spatial_orbitals=num_spatial_orbitals,
    num_particles=num_particles,
    qubit_mapper=mapper
)

# UCCSD ansatz
ansatz = UCCSD(
    num_spatial_orbitals=num_spatial_orbitals,
    num_particles=num_particles,
    qubit_mapper=mapper,
    initial_state=hf_state,
    reps=1
)

print(f"UCCSD parameters: {ansatz.num_parameters}")
print(f"Circuit depth: {ansatz.decompose().depth()}")

Hardware-Efficient Ansatz

Shorter circuits that match hardware connectivity:

from qiskit.circuit.library import EfficientSU2

ansatz = EfficientSU2(
    num_qubits=qubit_op.num_qubits,
    reps=2,
    entanglement='linear',  # Match chip topology
    insert_barriers=True
)

# Prepend Hartree-Fock initialization
full_ansatz = hf_state.compose(ansatz)

ADAPT-VQE: Iterative Ansatz Construction

ADAPT-VQE builds the ansatz one operator at a time, selecting the most impactful operation at each step:

from qiskit_nature.second_q.circuit.library import UCC

def adapt_vqe(hamiltonian, mapper, problem, max_iterations=20, threshold=1e-5):
    """
    Simplified ADAPT-VQE: iteratively grow the ansatz.
    """
    num_so = problem.num_spatial_orbitals
    num_particles = problem.num_particles

    # Operator pool: all single and double excitations
    pool = UCC(
        num_spatial_orbitals=num_so,
        num_particles=num_particles,
        qubit_mapper=mapper,
        excitations='sd'
    )

    selected_ops = []
    current_params = []

    for iteration in range(max_iterations):
        # Compute gradient of each pool operator
        gradients = []
        for op_idx in range(len(pool.operators)):
            grad = compute_operator_gradient(
                op_idx, pool, selected_ops, current_params, hamiltonian
            )
            gradients.append(abs(grad))

        # Select operator with largest gradient
        best_idx = max(range(len(gradients)), key=lambda i: gradients[i])

        if gradients[best_idx] < threshold:
            print(f"Converged at iteration {iteration}")
            break

        selected_ops.append(best_idx)
        current_params.append(0.0)  # Initialize new parameter

        # Re-optimize all parameters
        current_params = optimize_parameters(
            selected_ops, current_params, pool, hamiltonian
        )

    return selected_ops, current_params

The Optimization Loop

Standard VQE with Qiskit

from qiskit_algorithms import VQE
from qiskit_algorithms.optimizers import COBYLA, SPSA, L_BFGS_B
from qiskit.primitives import Estimator

# Estimator computes expectation values
estimator = Estimator()

# Choose optimizer
optimizer = COBYLA(maxiter=500, tol=1e-6)

# Run VQE
vqe = VQE(
    estimator=estimator,
    ansatz=ansatz,
    optimizer=optimizer,
    initial_point=np.zeros(ansatz.num_parameters)
)

result = vqe.compute_minimum_eigenvalue(qubit_op)

print(f"VQE energy: {result.eigenvalue.real:.8f} Hartree")
print(f"Iterations: {result.cost_function_evals}")
print(f"Optimal parameters: {result.optimal_parameters}")

Manual VQE Loop (More Control)

import numpy as np
from scipy.optimize import minimize

def vqe_cost_function(params, ansatz, hamiltonian, estimator):
    """Evaluate energy for given parameters."""
    bound_circuit = ansatz.assign_parameters(params)
    job = estimator.run(bound_circuit, hamiltonian)
    energy = job.result().values[0]
    return energy.real

def run_vqe_manual(ansatz, hamiltonian, estimator, method='COBYLA'):
    """VQE with full control over the optimization."""
    energies = []

    def callback(params):
        energy = vqe_cost_function(params, ansatz, hamiltonian, estimator)
        energies.append(energy)
        if len(energies) % 10 == 0:
            print(f"  Step {len(energies)}: E = {energy:.8f}")
        return energy

    # Multiple random starts to avoid local minima
    best_result = None
    for trial in range(5):
        x0 = np.random.uniform(-np.pi, np.pi, ansatz.num_parameters)
        result = minimize(
            callback,
            x0,
            method=method,
            options={'maxiter': 200, 'rhobeg': 0.5}
        )
        if best_result is None or result.fun < best_result.fun:
            best_result = result

    return best_result, energies

Gradient-Based Optimization

For ansätze where gradients are available:

from qiskit_algorithms.gradients import ParamShiftEstimatorGradient

gradient = ParamShiftEstimatorGradient(estimator)

# Use gradient-based optimizer
optimizer = L_BFGS_B(maxiter=100)

vqe = VQE(
    estimator=estimator,
    ansatz=ansatz,
    optimizer=optimizer,
    gradient=gradient,
    initial_point=np.zeros(ansatz.num_parameters)
)

SPSA for Noisy Hardware

Simultaneous Perturbation Stochastic Approximation works well with noisy measurements:

optimizer = SPSA(
    maxiter=300,
    learning_rate=0.05,
    perturbation=0.1,
    # Calibrate learning rate and perturbation
    # based on initial energy landscape sampling
)

SPSA estimates gradients using only 2 circuit evaluations regardless of the number of parameters (vs. 2p for parameter-shift).

Measurement Optimization

Grouping Commuting Paulis

Pauli terms that commute can be measured simultaneously:

from qiskit.quantum_info import SparsePauliOp

# Group commuting terms to reduce measurement circuits
grouped_op = qubit_op.group_commuting()

print(f"Total Pauli terms: {len(qubit_op)}")
print(f"Measurement groups: {len(grouped_op)}")
# Grouping typically reduces measurements by 50-70%

Shot Budget Allocation

Not all terms contribute equally. Allocate more shots to high-weight terms:

def adaptive_shots(hamiltonian, total_budget=100000):
    """Allocate shots proportional to coefficient magnitude."""
    coeffs = np.abs(hamiltonian.coeffs)
    weights = coeffs / coeffs.sum()
    shots = (weights * total_budget).astype(int)
    shots = np.maximum(shots, 100)  # Minimum 100 shots per term
    return shots

Potential Energy Surface Scanning

Map how energy changes with molecular geometry:

def potential_energy_surface(distances):
    """Compute H2 energy at various bond lengths."""
    energies_vqe = []
    energies_exact = []

    for d in distances:
        driver = PySCFDriver(
            atom=f"H 0.0 0.0 0.0; H 0.0 0.0 {d}",
            basis="sto-3g"
        )
        problem = driver.run()
        hamiltonian = problem.hamiltonian.second_q_op()
        qubit_op = JordanWignerMapper().map(hamiltonian)

        # VQE
        vqe = VQE(estimator, ansatz, COBYLA(maxiter=200))
        vqe_result = vqe.compute_minimum_eigenvalue(qubit_op)
        energies_vqe.append(vqe_result.eigenvalue.real)

        # Exact (for comparison)
        from qiskit_algorithms import NumPyMinimumEigensolver
        exact = NumPyMinimumEigensolver()
        exact_result = exact.compute_minimum_eigenvalue(qubit_op)
        energies_exact.append(exact_result.eigenvalue.real)

    return energies_vqe, energies_exact

distances = np.arange(0.3, 3.0, 0.1)
vqe_e, exact_e = potential_energy_surface(distances)

Error Mitigation for VQE

Zero-Noise Extrapolation

def zne_vqe(circuit, hamiltonian, noise_factors=[1.0, 1.5, 2.0]):
    """Run VQE at multiple noise levels and extrapolate to zero noise."""
    energies = []
    for factor in noise_factors:
        # Stretch circuit to amplify noise
        stretched = stretch_circuit(circuit, factor)
        energy = run_on_hardware(stretched, hamiltonian)
        energies.append(energy)

    # Richardson extrapolation
    # For linear model: E(λ) = E₀ + aλ
    from numpy.polynomial import polynomial as P
    coeffs = P.polyfit(noise_factors, energies, deg=len(noise_factors)-1)
    zero_noise_energy = P.polyval(0, coeffs)
    return zero_noise_energy

Scaling Considerations

Molecule	Qubits (JW)	UCCSD Parameters	Estimated Gate Count
H₂	4	3	~50
LiH	12	44	~3,000
H₂O	14	52	~5,000
N₂	20	160	~50,000
FeMoCo	~100+	~10,000+	~10⁸+

FeMoCo (the nitrogen fixation catalyst) is the grand challenge of quantum chemistry — solving it would revolutionize fertilizer production. It requires fault-tolerant quantum computers that are likely a decade away.

One thing to remember: VQE is a framework, not a single algorithm — its performance depends entirely on three choices you make: the qubit mapping (how efficiently you encode the molecule), the ansatz (whether it can represent the true ground state), and the optimizer (whether it finds the global minimum in a noisy landscape). Mastering these choices is what separates a textbook VQE from one that produces chemically meaningful results.

pythonquantum-computingvqequantum-chemistry