Mathematical Formulation of Variational Quantum Eigensolver (VQE) for Molecular Hamiltonian Mapping
The Variational Quantum Eigensolver is a hybrid quantum-classical algorithm used to approximate the ground-state energy of a quantum system. In molecular simulation, the method begins with an electronic structure Hamiltonian, maps fermionic operators to qubit operators, prepares a parameterized trial state, and minimizes the measured expectation value of the Hamiltonian. This document defines the mathematical formulation and provides a Python reference implementation for a minimal VQE execution loop.
The electronic Hamiltonian in second quantization can be written as:
\[ H = \sum_{p,q} h_{pq} a_p^\dagger a_q + \frac{1}{2} \sum_{p,q,r,s} h_{pqrs} a_p^\dagger a_q^\dagger a_r a_s \]
Here, \(a_p^\dagger\) and \(a_p\) are fermionic creation and annihilation operators. The one-electron integrals \(h_{pq}\) and two-electron integrals \(h_{pqrs}\) are computed from molecular orbitals. A fermion-to-qubit transformation such as Jordan-Wigner or Bravyi-Kitaev maps the fermionic Hamiltonian to a weighted sum of Pauli strings:
\[ H_q = \sum_{j=1}^{M} c_j P_j \]
where \(c_j \in \mathbb{R}\) and:
\[ P_j \in \{I,X,Y,Z\}^{\otimes n} \]
VQE prepares a parameterized state \( |\psi(\theta)\rangle \) using an ansatz circuit. The objective function is the Hamiltonian expectation value:
\[ E(\theta)= \langle\psi(\theta)|H_q|\psi(\theta)\rangle = \sum_{j=1}^{M} c_j \langle\psi(\theta)|P_j|\psi(\theta)\rangle \]
The variational principle guarantees:
\[ E_0 \leq E(\theta) \]
where \(E_0\) is the exact ground-state energy. The optimization problem is:
\[ \theta^\star=\arg\min_{\theta}E(\theta) \]
Reference Implementation
from __future__ import annotations
from dataclasses import dataclass
from typing import List, Sequence, Tuple
import numpy as np
from scipy.optimize import minimize
from qiskit import QuantumCircuit, transpile
from qiskit.quantum_info import SparsePauliOp, Statevector
@dataclass(frozen=True)
class VQEResult:
energy: float
parameters: np.ndarray
evaluations: int
def build_ansatz(num_qubits: int, parameters: Sequence[float]) -> QuantumCircuit:
circuit = QuantumCircuit(num_qubits)
cursor = 0
for qubit in range(num_qubits):
circuit.ry(parameters[cursor], qubit)
cursor += 1
for qubit in range(num_qubits - 1):
circuit.cx(qubit, qubit + 1)
for qubit in range(num_qubits):
circuit.rz(parameters[cursor], qubit)
cursor += 1
return circuit
def expectation_value(
hamiltonian: SparsePauliOp,
circuit: QuantumCircuit
) -> float:
state = Statevector.from_instruction(circuit)
value = state.expectation_value(hamiltonian)
return float(np.real(value))
def make_h2_like_hamiltonian() -> SparsePauliOp:
labels = [
"II",
"IZ",
"ZI",
"ZZ",
"XX",
]
coefficients = [
-1.052373245772859,
0.39793742484318045,
-0.39793742484318045,
-0.01128010425623538,
0.18093119978423156,
]
return SparsePauliOp.from_list(list(zip(labels, coefficients)))
def run_vqe(
hamiltonian: SparsePauliOp,
num_qubits: int,
initial_parameters: np.ndarray
) -> VQEResult:
evaluations = 0
def objective(theta: np.ndarray) -> float:
nonlocal evaluations
evaluations += 1
circuit = build_ansatz(num_qubits, theta)
compiled = transpile(circuit, basis_gates=["u3", "cx"], optimization_level=1)
return expectation_value(hamiltonian, compiled)
result = minimize(
objective,
initial_parameters,
method="COBYLA",
options={
"maxiter": 150,
"rhobeg": 0.25,
"tol": 1e-6,
},
)
return VQEResult(
energy=float(result.fun),
parameters=np.array(result.x),
evaluations=evaluations,
)
if __name__ == "__main__":
hamiltonian = make_h2_like_hamiltonian()
num_qubits = 2
initial = np.array([0.1, -0.2, 0.05, 0.3], dtype=float)
result = run_vqe(
hamiltonian=hamiltonian,
num_qubits=num_qubits,
initial_parameters=initial
)
print("estimated ground-state energy:", result.energy)
print("optimized parameters:", result.parameters)
print("objective evaluations:", result.evaluations)The reference Hamiltonian is a compact two-qubit operator used to demonstrate the VQE control flow. A full molecular workflow would obtain the coefficients from an electronic structure package, select an active space, choose a fermion-to-qubit mapping, and optionally apply qubit tapering. The Pauli expansion remains the central execution format because each term can be measured through basis rotation and computational-basis sampling.
For a Pauli string containing \(X\), \(Y\), and \(Z\) operators, measurement requires rotating the qubit basis before sampling. The expectation value of a Pauli term is estimated as:
\[ \langle P_j\rangle \approx \frac{1}{N} \sum_{k=1}^{N} m_k \]
where \(m_k\in\{-1,+1\}\) is the measured eigenvalue for shot \(k\), and \(N\) is the number of shots. The full energy estimator is:
\[ \hat{E}(\theta)= \sum_{j=1}^{M} c_j\hat{\mu}_j \]
Its sampling variance is:
\[ \operatorname{Var}[\hat{E}] = \sum_{j=1}^{M} c_j^2 \operatorname{Var}[\hat{\mu}_j] + 2\sum_{i<j}c_ic_j\operatorname{Cov}(\hat{\mu}_i,\hat{\mu}_j) \]
In practical implementations, commuting Pauli terms may be grouped to reduce measurement cost. Ansatz selection is equally important. A shallow hardware-efficient ansatz may execute on near-term devices but can suffer from barren plateaus or insufficient expressivity. A chemically motivated ansatz may encode domain structure more accurately but can produce deeper circuits.
The optimizer should be selected according to the noise model. Gradient-free optimizers are often used for shot-based estimates, while analytic or parameter-shift gradients may be used when measurement noise is controlled. The parameter-shift rule for a gate generated by a Pauli operator can be written as:
\[ \frac{\partial E(\theta)}{\partial \theta_k} = \frac{1}{2} \left[ E\left(\theta_k+\frac{\pi}{2}\right) – E\left(\theta_k-\frac{\pi}{2}\right) \right] \]
A complete VQE record should store the molecular geometry, basis set, mapping method, Hamiltonian coefficients, ansatz definition, optimizer settings, initial parameters, measurement grouping strategy, shot count, and final convergence trace. Without these fields, the reported energy value is not reproducible.