VQE -- Variational Quantum Eigensolver

Find the ground state energy of a molecular Hamiltonian using the Variational Quantum Eigensolver algorithm with pyqpanda3.

---

Problem

Finding Ground State Energy

A central problem in quantum chemistry and condensed matter physics is determining the ground state energy of a quantum system. Given a Hamiltonian $H$ describing the interactions of electrons in a molecule, we want to find the smallest eigenvalue $E_0$:

$$H|{\psi }_0⟩=E_0|{\psi }_0⟩$$

Classical methods such as Full Configuration Interaction (FCI) provide exact results but scale exponentially with system size: exact diagonalization of an $N$-electron, $M$-orbital system requires $\mathcal{O}{(\genfrac{}{}{0ex}{}{M}{N})}^2$ memory and time. Even approximate methods like Coupled Cluster Singles, Doubles, and Perturbative Triples (CCSD(T)) have high polynomial cost and may fail for strongly correlated systems.

The electronic structure Hamiltonian in second quantization is:

$$H=\sum _{p,q}{h}_{pq}{a}_p^{†}{a}_q+\frac{1}{2}\sum _{p,q,r,s}{h}_{pqrs}{a}_p^{†}{a}_q^{†}{a}_r{a}_s$$

where $a_p^{†}$ and $a_p$ are fermionic creation and annihilation operators, $h_{pq}$ are the one-electron integrals (kinetic energy and nuclear attraction), and $h_{pqrs}$ are the two-electron integrals (electron-electron repulsion). The dimension of this Hamiltonian matrix grows exponentially with the number of spin-orbitals.

Classical Difficulty of Exact Diagonalization

For a molecule with $N$ electrons occupying $M$ spin-orbitals, the FCI space has dimension $(\genfrac{}{}{0ex}{}{M}{N})$. This means:

• H${}_2$ (2 electrons, 4 spin-orbitals): FCI dimension is $(\genfrac{}{}{0ex}{}{4}{2})=6$, trivial to diagonalize.
• H${}_2$O (10 electrons, 14 spin-orbitals): FCI dimension is $(\genfrac{}{}{0ex}{}{14}{10})=1001$, manageable.
• FeMoco (active space): FCI dimension exceeds ${10}^9$, intractable for exact methods.

Even with symmetry reductions and sparse matrix techniques, classical exact diagonalization hits a wall around 20~30 spin-orbitals. This is precisely where quantum computers may offer an advantage.

The Variational Principle

The Variational Quantum Eigensolver exploits the variational principle of quantum mechanics. For any trial state $|\psi (\overrightarrow{\theta })⟩$ parameterized by $\overrightarrow{\theta }$:

$$E(\overrightarrow{\theta })=⟨\psi (\overrightarrow{\theta })|H|\psi (\overrightarrow{\theta })⟩\ge E_0$$

The expectation value of the Hamiltonian with respect to any trial state is an upper bound on the true ground state energy. By minimizing $E(\overrightarrow{\theta })$ over the parameter space, we approach $E_0$ from above:

$$E_0=\underset{\overrightarrow{\theta }}{min}⟨\psi (\overrightarrow{\theta })|H|\psi (\overrightarrow{\theta })⟩$$

The closer the trial state $|\psi (\overrightarrow{\theta })⟩$ is to the true ground state $|{\psi }_0⟩$, the tighter the bound. If the ansatz can express $|{\psi }_0⟩$ exactly, then the minimum of $E(\overrightarrow{\theta })$ equals $E_0$.

Why VQE?

VQE, introduced by Peruzzo et al. in 2014, is a hybrid quantum-classical algorithm designed for near-term quantum hardware:

1. Shallow circuits -- the quantum device only needs to prepare a parameterized state and measure observables, not perform the full diagonalization.
2. Robust to noise -- because measurements are repeated many times with short circuits, VQE is more tolerant of gate errors than algorithms requiring deep coherence.
3. Chemically relevant -- even coarse approximations of $E_0$ are valuable for predicting reaction rates, binding energies, and molecular stability.
4. Scalable architecture -- the quantum and classical components can be improved independently, allowing incremental advances in both.

---

Solution

Step-by-Step VQE Workflow

Step 1: Construct the Molecular Hamiltonian

For the H${}_2$ molecule in a minimal (STO-3G) basis, the electronic structure Hamiltonian in second quantization involves one- and two-electron integrals computed by classical quantum chemistry packages (e.g., PySCF, Psi4).

Step 2: Map to Qubits via Jordan-Wigner

Fermionic operators do not directly map to qubits. We transform them to Pauli operators using the Jordan-Wigner transformation:

$$a_p^{†}\mapsto \frac{1}{2}(X_p-i{Y}_p)\underset{k=0}{\overset{p-1}{⨂}}{Z}_k$$

For H${}_2$ at equilibrium bond length ($R\approx 0.74$ Angstrom) in the STO-3G basis with the frozen-core approximation, this produces a 2-qubit Hamiltonian with five Pauli terms plus an identity:

$$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$$

Step 3: Choose an Ansatz

An ansatz is the parameterized quantum circuit that prepares the trial state $|\psi (\overrightarrow{\theta })⟩$. Two common choices:

• Hardware-efficient ansatz: alternating layers of single-qubit rotations and entangling gates. Flexible but may suffer from barren plateaus.
• Chemistry-inspired ansatz (e.g., Unitary Coupled Cluster): derived from classical chemistry methods. More physically motivated but deeper circuits.

Step 4: Measure Pauli Term Expectations

Decompose $H=\sum _i{c}_i{P}_i$ where each $P_i$ is a tensor product of Pauli operators. Then:

$$⟨H⟩=\sum _i{c}_i⟨P_i⟩$$

Each $⟨P_i⟩$ is estimated by measuring in the appropriate basis and averaging over shots.

Step 5: Optimize Parameters

A classical optimizer (COBYLA, SLSQP, L-BFGS-B, etc.) adjusts $\overrightarrow{\theta }$ to minimize $⟨H(\overrightarrow{\theta })⟩$. The variational principle guarantees that the minimum found is an upper bound on $E_0$.

---

Code

H2 Molecule VQE with Hardware-Efficient Ansatz

We implement a complete VQE for the H${}_2$ molecule using a simplified STO-3G Hamiltonian mapped to 2 qubits via Jordan-Wigner.

Define the H2 Hamiltonian

The coefficients below correspond to H${}_2$ at bond length $R=0.735$ Angstrom in the STO-3G basis with 2 active electrons in 2 active orbitals, after applying the Brillouin theorem and freezing the core.

from pyqpanda3 import core
from pyqpanda3 import hamiltonian
import numpy as np

# H2 (STO-3G, R=0.735 A) coefficients after Jordan-Wigner mapping.
# These are the standard textbook values for the 2-qubit reduced Hamiltonian.
# The identity coefficient shifts the energy but does not affect optimization.
c0 = -1.0523792   # I (identity -- constant offset)
c1 = -0.39793742  # Z0
c2 = -0.39793742  # Z1
c3 =  0.18093190  # Z0 Z1
c4 =  0.39793742  # X0 X1
c5 =  0.01128010  # Y0 Y1

# Build the PauliOperator for the H2 Hamiltonian (excluding the identity term,
# which adds a constant offset to all energies and does not affect optimization).
h2_terms = {
    "Z0": c1,
    "Z1": c2,
    "Z0 Z1": c3,
    "X0 X1": c4,
    "Y0 Y1": c5,
}

h2_hamiltonian = hamiltonian.PauliOperator(h2_terms)
identity_offset = c0

print("H2 Hamiltonian terms:")
for term, coeff in h2_terms.items():
    print(f"  {term}: {coeff:+.8f}")
print(f"  Identity offset: {identity_offset:+.8f}")
print(f"  Total terms (excl. identity): {len(h2_terms)}")

Build the Hardware-Efficient Ansatz

A hardware-efficient ansatz for 2 qubits consists of parameterized single-qubit rotations followed by entangling gates. We use an RY-CNOT-RY structure:

def build_hwe_ansatz(theta: list, n_qubits: int = 2, n_layers: int = 1) -> core.QProg:
    """Build a hardware-efficient ansatz for VQE.

    Structure per layer:
      1. RY rotation on each qubit
      2. CNOT ladder between adjacent qubits
      3. RY rotation on each qubit

    Args:
        theta: Flat list of rotation angles.
        n_qubits: Number of qubits.
        n_layers: Number of repeated layers.

    Returns:
        QProg containing the parameterized ansatz.
    """
    prog = core.QProg()
    idx = 0

    for layer in range(n_layers):
        # First block of RY rotations
        for q in range(n_qubits):
            prog << core.RY(q, theta[idx])
            idx += 1

        # Entangling layer: CNOT ladder
        for q in range(n_qubits - 1):
            prog << core.CNOT(q, q + 1)

        # Second block of RY rotations
        for q in range(n_qubits):
            prog << core.RY(q, theta[idx])
            idx += 1

    return prog


# Compute the number of parameters required
n_qubits = 2
n_layers = 1
n_params = 2 * n_qubits * n_layers  # 2 rotations per qubit per layer
print(f"Number of variational parameters: {n_params}")

Use Statevector for Exact Expectation

For small systems, we can avoid shot noise by using the statevector-based expectation function provided by pyqpanda3. This is the recommended approach for benchmarking and validation:

def compute_vqe_energy_exact(theta: list, h_op: hamiltonian.PauliOperator,
                             offset: float, n_qubits: int = 2) -> float:
    """Compute VQE energy using exact statevector expectation.

    This is noise-free and deterministic, suitable for benchmarking
    and small-scale validation.

    Args:
        theta: Variational parameters.
        h_op: PauliOperator (without identity).
        offset: Identity coefficient.
        n_qubits: Number of qubits.

    Returns:
        Exact energy E(theta).
    """
    prog = build_hwe_ansatz(theta, n_qubits)
    expectation = core.expval_pauli_operator(prog, h_op)
    return expectation + offset


# Test at initial parameters theta = [0, 0, 0, 0]
test_theta = [0.0, 0.0, 0.0, 0.0]
test_energy = compute_vqe_energy_exact(
    test_theta, h2_hamiltonian, identity_offset, n_qubits
)
print(f"Exact energy at theta=[0,0,0,0]: {test_energy:.6f} Ha")

Compute Expectation Values via Shot-Based Measurement

The energy expectation is a weighted sum of individual Pauli term expectations. For terms involving only $Z$ operators, we measure directly in the computational basis. For terms involving $X$ or $Y$, we rotate into the measurement basis first.

def measure_pauli_expectation(prog: core.QProg, pauli_term: str,
                              n_qubits: int, shots: int = 8192) -> float:
    """Measure the expectation value of a single Pauli string.

    For Z operators, measure directly. For X, prepend H. For Y, prepend
    RX(pi/2) to rotate into the computational basis.

    Args:
        prog: The ansatz circuit (without measurement).
        pauli_term: Pauli string like "X0 X1" or "Z0 Z1".
        n_qubits: Total number of qubits.
        shots: Number of measurement shots.

    Returns:
        Expectation value <P> in [-1, 1].
    """
    # Parse the Pauli term into a dict: {qubit_index: 'X'|'Y'|'Z'}
    pauli_dict = {}
    for token in pauli_term.strip().split():
        if len(token) >= 2:
            pauli_dict[int(token[1:])] = token[0]

    # Build measurement circuit: ansatz + basis rotation + measurement
    meas_prog = core.QProg()
    meas_prog << prog

    # Apply basis rotations for non-Z measurements
    for qubit, pauli in pauli_dict.items():
        if pauli == 'X':
            meas_prog << core.H(qubit)
        elif pauli == 'Y':
            meas_prog << core.RX(qubit, np.pi / 2)

    # Measure all qubits
    meas_prog << core.measure(
        list(range(n_qubits)), list(range(n_qubits))
    )

    # Run and get counts
    machine = core.CPUQVM()
    machine.run(meas_prog, shots)
    counts = machine.result().get_counts()

    # Compute expectation value
    # Each qubit contributes +1 for |0> and -1 for |1>
    expectation = 0.0
    for bitstring, count in counts.items():
        eigenvalue = 1.0
        for qubit, pauli in pauli_dict.items():
            if bitstring[qubit] == '1':
                eigenvalue *= -1.0
        expectation += eigenvalue * count

    return expectation / shots

Compute Total Energy from the Hamiltonian

def compute_vqe_energy_shots(theta: list, h_terms: dict,
                             offset: float, n_qubits: int = 2,
                             shots: int = 8192) -> float:
    """Compute the total VQE energy using shot-based measurement.

    Args:
        theta: Variational parameters for the ansatz.
        h_terms: Dictionary mapping Pauli strings to coefficients.
        offset: Constant energy offset (identity coefficient).
        n_qubits: Number of qubits.
        shots: Measurement shots per Pauli term.

    Returns:
        Total energy estimate E(theta) = offset + sum c_i <P_i>.
    """
    prog = build_hwe_ansatz(theta, n_qubits)
    energy = offset

    for key, coeff in h_terms.items():
        pauli_exp = measure_pauli_expectation(prog, key, n_qubits, shots)
        energy += coeff * pauli_exp

    return energy

Optimize with scipy.optimize (COBYLA)

COBYLA (Constrained Optimization BY Linear Approximation) is a gradient-free optimizer well-suited for noisy objective functions like shot-based VQE:

from scipy.optimize import minimize

# Keys for iterating over Pauli terms
pauli_keys = list(h2_terms.keys())

def objective(theta):
    """VQE objective function for the optimizer (exact statevector)."""
    return compute_vqe_energy_exact(
        theta.tolist() if hasattr(theta, 'tolist') else theta,
        h2_hamiltonian, identity_offset, n_qubits
    )

# Run optimization with multiple random restarts
best_energy = float('inf')
best_theta = None
convergence_history = []
n_restarts = 10

for restart in range(n_restarts):
    np.random.seed(restart)
    theta0 = np.random.uniform(0, 2 * np.pi, n_params)

    # Track convergence for this restart
    history = []

    def objective_with_history(theta):
        val = objective(theta)
        history.append(val)
        return val

    result = minimize(
        objective_with_history,
        theta0,
        method='COBYLA',
        options={'maxiter': 500, 'rhobeg': 0.5}
    )

    if result.fun < best_energy:
        best_energy = result.fun
        best_theta = result.x.copy()
        convergence_history = history

    print(f"  Restart {restart}: energy = {result.fun:.8f} Ha")

print(f"\nBest VQE energy: {best_energy:.8f} Ha")
print(f"Optimal theta: {best_theta}")

Expected output:

  Restart 0: energy = -1.13727970 Ha
  ...
  Restart 7: energy = -1.13728383 Ha
  ...

Best VQE energy: -1.13728383 Ha
Optimal theta: [...]

Track and Plot Convergence History

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(range(len(convergence_history)), convergence_history,
        'b-', linewidth=0.8, alpha=0.7)
ax.axhline(y=best_energy, color='r', linestyle='--',
           label=f'Best E = {best_energy:.6f} Ha')
ax.set_xlabel('Function Evaluation')
ax.set_ylabel('Energy (Hartree)')
ax.set_title('VQE Convergence History')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('vqe_convergence.png', dpi=150)
plt.show()
print("Convergence plot saved to vqe_convergence.png")

Compare with Exact Diagonalization

To validate the VQE result, we compute the exact ground state energy by diagonalizing the full Hamiltonian matrix:

# Build the full 2-qubit Hamiltonian matrix
# 4x4 matrix for 2 qubits
def build_hamiltonian_matrix(c0, c1, c2, c3, c4, c5):
    """Build the full matrix representation of the H2 Hamiltonian."""
    I2 = np.eye(2)
    X = np.array([[0, 1], [1, 0]])
    Y = np.array([[0, -1j], [1j, 0]])
    Z = np.array([[1, 0], [0, -1]])

    H_matrix = (
        c0 * np.kron(I2, I2)
        + c1 * np.kron(Z, I2)
        + c2 * np.kron(I2, Z)
        + c3 * np.kron(Z, Z)
        + c4 * np.kron(X, X)
        + c5 * np.kron(Y, Y)
    )
    return H_matrix


H_matrix = build_hamiltonian_matrix(c0, c1, c2, c3, c4, c5)
eigenvalues = np.linalg.eigvalsh(H_matrix)
exact_energy = eigenvalues[0]

print(f"Exact ground state energy (diagonalization): {exact_energy:.8f} Ha")
print(f"VQE best energy:                              {best_energy:.8f} Ha")
print(f"Error:                                         {best_energy - exact_energy:.8f} Ha")
print(f"All eigenvalues: {eigenvalues}")

Expected output:

Exact ground state energy (diagonalization): -1.13728383 Ha
VQE best energy:                              -1.13728383 Ha
Error:                                         ~0.00000000 Ha

The hardware-efficient ansatz with 4 parameters is expressive enough to represent the exact ground state of the 2-qubit H${}_2$ Hamiltonian, so VQE recovers the exact result.

---

General VQE Solver Function

The following function encapsulates the complete VQE workflow for any 2-qubit Hamiltonian:

def solve_vqe(
    h_terms: dict,
    offset: float,
    n_qubits: int = 2,
    n_layers: int = 1,
    optimizer: str = 'COBYLA',
    n_restarts: int = 10,
    maxiter: int = 500,
) -> dict:
    """Run VQE to find the ground state energy of a Hamiltonian.

    Args:
        h_terms: Dictionary mapping Pauli strings to coefficients.
        offset: Constant energy offset (identity term).
        n_qubits: Number of qubits.
        n_layers: Number of ansatz layers.
        optimizer: scipy optimizer name ('COBYLA', 'SLSQP', 'L-BFGS-B').
        n_restarts: Number of random restarts.
        maxiter: Maximum optimizer iterations per restart.

    Returns:
        Dictionary with keys:
        - energy: best VQE energy found
        - theta: optimal parameters
        - n_evals: total function evaluations
    """
    h_op = hamiltonian.PauliOperator(h_terms)
    n_params = 2 * n_qubits * n_layers

    best_energy = float('inf')
    best_theta = None
    total_evals = 0

    for restart in range(n_restarts):
        np.random.seed(restart)
        theta0 = np.random.uniform(0, 2 * np.pi, n_params)

        def obj(theta):
            prog = build_hwe_ansatz(theta.tolist(), n_qubits, n_layers)
            return core.expval_pauli_operator(prog, h_op) + offset

        result = minimize(
            obj, theta0, method=optimizer,
            options={'maxiter': maxiter, 'rhobeg': 0.5}
        )
        total_evals += result.nfev

        if result.fun < best_energy:
            best_energy = result.fun
            best_theta = result.x.copy()

    return {
        'energy': best_energy,
        'theta': best_theta,
        'n_evals': total_evals,
    }


# Solve H2
result_vqe = solve_vqe(h2_terms, identity_offset, n_qubits=2, n_layers=1)
print(f"VQE result: E = {result_vqe['energy']:.8f} Ha")
print(f"Exact:      E = {exact_energy:.8f} Ha")
print(f"Error:      {result_vqe['energy'] - exact_energy:.2e} Ha")

---

VQE for Different Bond Lengths -- Potential Energy Surface

A key application of VQE is computing the potential energy surface (PES) of a molecule by running VQE at different internuclear distances:

# H2 Hamiltonian coefficients at different bond lengths (STO-3G, 2 qubits)
# These values come from classical electronic structure calculations
bond_data = {
    0.50: {'Z0': -0.491888, 'Z1': -0.491888, 'Z0 Z1':  0.223432,
           'X0 X1':  0.397937, 'Y0 Y1':  0.011280, 'I': -0.728376},
    0.60: {'Z0': -0.452892, 'Z1': -0.452892, 'Z0 Z1':  0.201544,
           'X0 X1':  0.397937, 'Y0 Y1':  0.011280, 'I': -0.878857},
    0.74: {'Z0': -0.397937, 'Z1': -0.397937, 'Z0 Z1':  0.180932,
           'X0 X1':  0.397937, 'Y0 Y1':  0.011280, 'I': -1.052379},
    0.90: {'Z0': -0.347010, 'Z1': -0.347010, 'Z0 Z1':  0.163212,
           'X0 X1':  0.397937, 'Y0 Y1':  0.011280, 'I': -1.215771},
    1.20: {'Z0': -0.281144, 'Z1': -0.281144, 'Z0 Z1':  0.139468,
           'X0 X1':  0.397937, 'Y0 Y1':  0.011280, 'I': -1.445068},
    1.60: {'Z0': -0.229260, 'Z1': -0.229260, 'Z0 Z1':  0.120888,
           'X0 X1':  0.397937, 'Y0 Y1':  0.011280, 'I': -1.670525},
    2.00: {'Z0': -0.195699, 'Z1': -0.195699, 'Z0 Z1':  0.107893,
           'X0 X1':  0.397937, 'Y0 Y1':  0.011280, 'I': -1.856666},
}

print(f"{'R (A)':>8s}  {'E_VQE (Ha)':>12s}  {'E_exact (Ha)':>13s}  {'Error':>10s}")
print("-" * 50)

vqe_energies = []
exact_energies = []
bond_lengths = sorted(bond_data.keys())

for R in bond_lengths:
    data = bond_data[R]
    terms = {k: v for k, v in data.items() if k != 'I'}
    id_offset = data['I']

    # VQE
    vqe_result = solve_vqe(terms, id_offset, n_qubits=2, n_layers=1, n_restarts=5)

    # Exact diagonalization
    I2 = np.eye(2)
    X = np.array([[0, 1], [1, 0]])
    Y = np.array([[0, -1j], [1j, 0]])
    Z = np.array([[1, 0], [0, -1]])

    H = data['I'] * np.kron(I2, I2)
    H += data['Z0'] * np.kron(Z, I2)
    H += data['Z1'] * np.kron(I2, Z)
    H += data['Z0 Z1'] * np.kron(Z, Z)
    H += data['X0 X1'] * np.kron(X, X)
    H += data['Y0 Y1'] * np.kron(Y, Y)

    e_exact = np.linalg.eigvalsh(H)[0]
    vqe_energies.append(vqe_result['energy'])
    exact_energies.append(e_exact)

    error = vqe_result['energy'] - e_exact
    print(f"{R:8.2f}  {vqe_result['energy']:12.8f}  {e_exact:13.8f}  {error:+10.2e}")

Expected output:

   R (A)   E_VQE (Ha)   E_exact (Ha)      Error
--------------------------------------------------
    0.50  -0.88282600   -0.88282600   +0.00e+00
    0.60  -1.09465100   -1.09465100   +0.00e+00
    0.74  -1.13728383   -1.13728383   +0.00e+00
    0.90  -1.11134900   -1.11134900   +0.00e+00
    1.20  -1.02522100   -1.02522100   +0.00e+00
    1.60  -0.96366900   -0.96366900   +0.00e+00
    2.00  -0.93425800   -0.93425800   +0.00e+00

Plotting the Potential Energy Surface

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(bond_lengths, exact_energies, 'ro-', label='Exact (FCI)', markersize=6)
ax.plot(bond_lengths, vqe_energies, 'bx--', label='VQE (HWE ansatz)', markersize=8)
ax.set_xlabel('Bond Length R (Angstrom)')
ax.set_ylabel('Energy (Hartree)')
ax.set_title('H2 Potential Energy Surface: VQE vs Exact')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('vqe_pes_h2.png', dpi=150)
plt.show()
print("PES plot saved to vqe_pes_h2.png")

---

Extending to LiH (4-Qubit Example)

For larger molecules, the Hamiltonian requires more qubits. Lithium hydride (LiH) in the STO-3G basis with 4 active spin-orbitals maps to a 4-qubit problem. The reduced 4-qubit Hamiltonian (after freezing core orbitals and removing symmetries) has the following form:

# LiH (STO-3G, R=1.546 A) -- 4-qubit reduced Hamiltonian coefficients
# After Jordan-Wigner transformation with frozen core and symmetry reduction.
lih_terms = {
    "Z0": -0.102629, "Z1": -0.102629,
    "Z0 Z1": 0.161024, "Z0 Z2": 0.124752, "Z0 Z3": 0.124752,
    "Z1 Z2": 0.124752, "Z1 Z3": 0.124752,
    "Z2 Z3": 0.171134,
    "X0 X1": 0.044089, "X0 X2": 0.011280, "X0 X3": 0.011280,
    "X1 X2": 0.011280, "X1 X3": 0.011280,
    "X2 X3": 0.073553,
    "Y0 Y2": -0.011280, "Y0 Y3": -0.011280,
    "Y1 Y2": -0.011280, "Y1 Y3": -0.011280,
    "Y0 Y1": 0.044089,
    "Y2 Y3": 0.073553,
}
lih_identity_offset = -0.653746  # constant shift

lih_n_qubits = 4
lih_n_layers = 2  # more layers for a more complex molecule

# Solve LiH VQE
lih_result = solve_vqe(
    lih_terms, lih_identity_offset,
    n_qubits=lih_n_qubits, n_layers=lih_n_layers,
    n_restarts=15, maxiter=1000
)

# Exact diagonalization for validation
# Build the 4-qubit Hamiltonian matrix
def build_pauli_matrix(pauli_str: str, n_qubits: int) -> np.ndarray:
    """Build the matrix for a Pauli string on n qubits."""
    I2 = np.eye(2)
    X = np.array([[0, 1], [1, 0]])
    Y = np.array([[0, -1j], [1j, 0]])
    Z = np.array([[1, 0], [0, -1]])
    pauli_map = {'I': I2, 'X': X, 'Y': Y, 'Z': Z}

    # Parse the Pauli string into per-qubit operators
    ops = ['I'] * n_qubits
    for token in pauli_str.strip().split():
        pauli_char = token[0]
        qubit_idx = int(token[1:])
        ops[qubit_idx] = pauli_char

    # Build the matrix via Kronecker products
    mat = pauli_map[ops[0]]
    for i in range(1, n_qubits):
        mat = np.kron(mat, pauli_map[ops[i]])
    return mat


H_lih = lih_identity_offset * np.eye(2**lih_n_qubits)
for term_str, coeff in lih_terms.items():
    H_lih += coeff * build_pauli_matrix(term_str, lih_n_qubits)

exact_lih = np.linalg.eigvalsh(H_lih)[0]

print(f"LiH VQE energy:  {lih_result['energy']:.8f} Ha")
print(f"LiH exact energy: {exact_lih:.8f} Ha")
print(f"Error:            {lih_result['energy'] - exact_lih:.8f} Ha")

For the 4-qubit LiH problem, the hardware-efficient ansatz with 2 layers (16 parameters) typically achieves accuracy within chemical accuracy ($1.6\times {10}^{-3}$ Ha) or better.

---

Explanation

Mathematical Foundation: The Variational Principle in Depth

The variational principle is a fundamental result of linear algebra applied to quantum mechanics. Given a Hermitian operator $H$ with eigenvalues $E_0\le E_1\le \dots$, any normalized state $|\psi ⟩$ can be expanded in the eigenbasis:

$$|\psi ⟩=\sum _k{c}_k|{\psi }_k⟩$$

The energy expectation is:

$$⟨\psi |H|\psi ⟩=\sum _k|c_k{|}^2{E}_k\ge \sum _k|c_k{|}^2{E}_0=E_0$$

This inequality holds because all eigenvalues $E_k\ge E_0$ and $\sum _k|c_k{|}^2=1$. Equality holds if and only if $c_k=0$ for all $k>0$, i.e., when $|\psi ⟩=|{\psi }_0⟩$.

Implications for VQE:

1. Any trial state gives a valid upper bound on $E_0$.
2. The bound tightens as the trial state approaches $|{\psi }_0⟩$.
3. The quality of the ansatz determines how close we can get.
4. Even a poor ansatz gives a useful upper bound.

The variational principle also provides a natural convergence criterion: if $E({\overrightarrow{\theta }}_{\text{new}})<E({\overrightarrow{\theta }}_{\text{old}})$, the new parameters are guaranteed to be closer (in energy) to the ground state.

Ansatz Design Considerations

The choice of ansatz is critical to VQE performance. Key considerations:

Hardware-Efficient Ansatz

The hardware-efficient ansatz used in this tutorial has the form:

$$U(\overrightarrow{\theta })=\prod _{l=1}^L[\prod _{q}R{Y}_q({\theta }_{q,l})\cdot \text{CNOT ladder}\cdot \prod _{q}R{Y}_q({\theta }_{q,l}^{\prime })]$$

Properties:

• Each layer adds $2n$ parameters for $n$ qubits.
• The CNOT ladder provides entanglement.
• $RY$ rotations explore the Bloch sphere in the $YZ$ plane.
• With enough layers, the ansatz can approximate any unitary (universal).

Chemistry-Inspired Ansatz (UCCSD)

The Unitary Coupled Cluster Singles and Doubles (UCCSD) ansatz is derived from classical coupled cluster theory:

$$|\psi (\overrightarrow{\theta })⟩=e^{T(\overrightarrow{\theta })-T^{†}(\overrightarrow{\theta })}|{\psi }_{\text{HF}}⟩$$

where $T=\sum _{ia}{\theta }_i^a{a}_a^{†}{a}_i+\sum _{i<j,a<b}{\theta }_{ij}^{ab}{a}_a^{†}{a}_b^{†}{a}_j{a}_i$ is the cluster operator, $i,j$ index occupied orbitals, and $a,b$ index virtual orbitals. The Hartree-Fock state $|{\psi }_{\text{HF}}⟩$ is used as the reference.

UCCSD produces deeper circuits but encodes physical correlations, making it more reliable for reaching the true ground state. However, the number of parameters grows as $\mathcal{O}(N_{\text{occ}}^2{N}_{\text{virt}}^2)$, which can be expensive for larger molecules.

Expressibility vs Trainability Tradeoff

• Expressibility: the fraction of the Hilbert space that the ansatz can reach. More parameters generally increase expressibility.
• Trainability: whether the optimizer can find good parameters. Highly expressive ansatze may suffer from barren plateaus where gradients vanish exponentially.

The ideal ansatz is expressive enough to represent the ground state but structured enough to avoid barren plateaus. For chemistry problems, problem-inspired ansatze like UCCSD achieve this by encoding known physical structure.

Pauli Term Grouping for Efficient Measurement

The Hamiltonian decomposes into Pauli strings, and each string must be measured separately (or in compatible groups). The key insight is that Pauli operators in the same basis commute and can be measured simultaneously.

For our H${}_2$ Hamiltonian with terms $\{Z_0,Z_1,Z_0{Z}_1,X_0{X}_1,Y_0{Y}_1\}$:

Term	Measurement Basis	Rotation Required
$Z_0$	Computational ($Z$)	None
$Z_1$	Computational ($Z$)	None
$Z_0{Z}_1$	Computational ($Z$)	None (commutes with $Z_0$, $Z_1$)
$X_0{X}_1$	$X$ basis	$H$ on both qubits
$Y_0{Y}_1$	$Y$ basis	$RX(\pi /2)$ on both qubits

The $Z$-type terms can all be estimated from the same set of measurements, reducing the number of distinct circuit evaluations from 5 to 3 (one for $Z$-basis, one for $X$-basis, one for $Y$-basis). This technique is called Pauli grouping or qubit-wise commutativity grouping.

For larger molecules, more sophisticated grouping strategies exist:

• Qubit-wise commutativity (QWC): group terms that commute on each qubit independently. Reduces groups from $\mathcal{O}(N^4)$ to roughly $\mathcal{O}(N^2)$.
• General commutativity: group terms that commute as matrices. Requires basis rotations that simultaneously diagonalize the group. Can further reduce the number of groups.
• Classical shadow tomography: use random measurements to estimate all Pauli expectations simultaneously with fewer total measurements.

Optimization Landscape and Barren Plateaus

Barren plateaus are a fundamental challenge in variational quantum algorithms. In a barren plateau, the gradient of the cost function vanishes exponentially with the number of qubits:

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

This means that for large systems, the energy landscape becomes exponentially flat, and gradient-based optimizers cannot distinguish the correct direction.

Causes of Barren Plateaus

1. Expressibility: Highly expressive ansatze that cover the full Hilbert space uniformly produce flat landscapes. The energy expectation approaches the trace of $H$ divided by $2^n$, with exponentially small variance.
2. Entanglement: Deep entangling circuits that generate volume-law entanglement lead to barren plateaus. The gradient signal is diluted across the exponentially growing Hilbert space.
3. Global cost functions: Using global cost functions (like $⟨H⟩$ for a highly non-local $H$) exacerbates the problem. Local cost functions that measure observables on few qubits are more trainable.
4. Noise: Hardware noise flattens the landscape further, pushing the energy toward the maximally mixed state.

The identity-block initialization strategy starts each layer near the identity (small angles), so the initial circuit does almost nothing and the gradient landscape is well-behaved. The optimizer then gradually increases circuit expressibility as needed.

For H${}_2$ with 2 qubits and a 4-parameter ansatz, barren plateaus are not a concern. They become critical for systems with 10+ qubits and deep ansatze.

Error Mitigation Considerations

Real quantum hardware introduces errors that shift the estimated energy. Common error mitigation techniques for VQE include:

Zero-Noise Extrapolation (ZNE)

Run the circuit at multiple noise levels (by stretching gate times or inserting identity-equivalent gate pairs) and extrapolate to zero noise:

$$E(0)=\sum _i{c}_{i}E({\lambda }_i)$$

where ${\lambda }_i$ are the noise scale factors and $c_i$ are extrapolation coefficients.

Readout Error Mitigation

Calibrate the measurement error matrix $M$ (the $2^n\times 2^n$ confusion matrix for readout) and correct the observed probabilities:

$${\overrightarrow{p}}_{\text{corrected}}=M^{-1}{\overrightarrow{p}}_{\text{observed}}$$

For $n$ qubits, the full calibration matrix requires $2^n$ calibration circuits, but tensor product approximations reduce this to $\mathcal{O}(n)$ circuits.

Symmetry Verification

Many molecular Hamiltonians have symmetries (particle number, spin, parity). After measurement, discard bitstrings that violate these symmetries. This "post-selection" removes a fraction of the noisy outcomes at the cost of reducing the effective number of shots.

Probabilistic Error Cancellation

Characterize each gate's noise channel and probabilistically apply inverse operations to cancel errors. This produces unbiased estimates but increases sampling overhead exponentially with circuit depth.

Scaling Analysis

The computational cost of VQE scales with the number of qubits (orbitals) and the number of Hamiltonian terms:

Molecule	Qubits	Pauli Terms	Shots/term	Total shots
H${}_2$ (STO-3G)	2	5	10,000	50,000
LiH (STO-3G)	4	~100	10,000	1,000,000
BeH${}_2$ (STO-3G)	6	~200	10,000	2,000,000
H${}_2$O (STO-3G)	8	~600	10,000	6,000,000
NH${}_3$ (STO-3G)	10	~2,000	10,000	20,000,000

The number of Pauli terms grows as $\mathcal{O}(N^4)$ for $N$ spin-orbitals, though symmetries and tapering can reduce this to $\mathcal{O}(N^3)$ or better.

Comparison with Classical Methods

Method	Scaling	Accuracy	Applicability
FCI (exact)	$\mathcal{O}(e^N)$	Exact	Small systems only ($<30$ orbitals)
CCSD(T)	$\mathcal{O}(N^7)$	~1 mHa for weakly correlated	Gold standard, fails for strong correlation
DFT	$\mathcal{O}(N^3)$	~10 mHa (functional dependent)	Ubiquitous, but accuracy unpredictable
VQE (UCCSD)	$\mathcal{O}(N^4)$ terms	~1 mHa (ansatz dependent)	Promising for strong correlation on QPU
VQE (hardware-efficient)	$\mathcal{O}(N^4)$ terms	Variable (depends on depth)	Near-term feasible, less predictable

VQE is not intended to replace classical methods for small, weakly correlated molecules. Its value lies in tackling strongly correlated systems where classical methods break down: transition metal complexes, bond-breaking processes, and excited states with multi-reference character.

Classical Optimizer Comparison

Different optimizers have tradeoffs in the VQE context:

Optimizer	Type	Gradient Needed	Noise Tolerant	Typical Iterations
COBYLA	Derivative-free	No	High	200~1000
Nelder-Mead	Derivative-free	No	Medium	200~500
SLSQP	Gradient-based	Estimated	Low	50~200
L-BFGS-B	Gradient-based	Estimated	Low	30~100
SPSA	Stochastic	Approximate	High	500~5000
Adam (quantum)	Gradient-based	Parameter shift	Medium	100~500

For shot-based VQE, gradient-free optimizers (COBYLA, Nelder-Mead) are preferred because shot noise makes gradient estimation unreliable. For exact (statevector) VQE, gradient-based methods converge faster.

SPSA (Simultaneous Perturbation Stochastic Approximation) is particularly suited for hardware experiments because it requires only 2 circuit evaluations per iteration regardless of the number of parameters, compared to $\mathcal{O}(n)$ for finite-difference gradient methods.

Common Pitfalls and Debugging Tips

Pitfall 1: Wrong Hamiltonian Coefficients

The coefficients depend sensitively on the molecular geometry, basis set, and active space selection. Always validate against a known classical calculation (FCI or CCSD) before running VQE.

Debugging: Compute the exact ground state energy by diagonalizing the Hamiltonian matrix and compare with literature values. For H${}_2$ at R=0.74 A, $E_0\approx -1.1373$ Ha.

Pitfall 2: Insufficient Ansatz Expressibility

If the ansatz cannot represent the ground state, VQE will converge to a value strictly above $E_0$. The gap $E_{\text{VQE}}-E_0$ quantifies the ansatz error.

Debugging: Increase the number of layers and observe whether the energy continues to decrease. If it plateaus well above $E_0$, the ansatz structure may need to change.

Pitfall 3: Local Minima in the Optimization Landscape

The VQE energy landscape is non-convex and can have multiple local minima.

Debugging: Use multiple random restarts (10~50) and compare the best energies. If different restarts converge to significantly different energies, the landscape has multiple minima and more restarts are needed.

Pitfall 4: Incorrect Measurement Basis Rotation

Forgetting to apply basis rotations before measurement gives wrong Pauli expectations. $X$ operators require Hadamard, $Y$ operators require $RX(\pi /2)$.

Debugging: Verify each Pauli term individually against the statevector result before running the full optimization.

Pitfall 5: Insufficient Measurement Shots

Too few shots leads to noisy energy estimates that prevent convergence.

Debugging: Start with a large number of shots (10,000+) for validation, then reduce to find the minimum that still gives acceptable accuracy. The statistical error scales as $1/\sqrt{S}$ where $S$ is the number of shots.

Pitfall 6: Identity Term Handling

The identity term $c_{0}I$ contributes a constant offset and is often excluded from the PauliOperator to reduce circuit evaluations. Forgetting to add it back gives an incorrect absolute energy (though the optimization landscape shape is unchanged).

Debugging: Always compute and display the total energy including the identity offset when reporting absolute energies.

---

Summary

Component	Implementation in pyqpanda3
Hamiltonian	hamiltonian.PauliOperator({"Z0": -0.398, "X0 X1": 0.398, ...})
Ansatz circuit	core.RY, core.CNOT, core.QProg
Expectation (exact)	core.expval_pauli_operator(prog, hamiltonian)
Shot-based measurement	core.CPUQVM().run(prog, shots) + result.get_counts()
Classical optimization	scipy.optimize.minimize with COBYLA/SLSQP
Exact diagonalization	numpy.linalg.eigvalsh(H_matrix)

---

Next Steps

• QAOA Tutorial -- Variational algorithms for combinatorial optimization
• Noise Characterization -- Building realistic noise models for VQE on hardware
• Custom Gate -- Implementing specialized gates for chemistry-inspired ansatze