VQE and QAOA

Variational algorithms on exact simulations, compared with exact answers: the full-CI energy of H₂, Lanczos ground energies of spin chains, and brute-force maximum cuts.

import matplotlib.pyplot as plt
import numpy as np

import metatron_qso as q

H₂ in the STO-3G basis

Two qubits after the parity mapping with two-qubit reduction. VQE with a one-layer RealAmplitudes ansatz reaches the full-CI energy.

h2 = q.h2_sto3g()
result = q.run_vqe(h2["hamiltonian"], q.real_amplitudes(2, 1), restarts=4, seed=0)
energy = result["energy"] + h2["nuclear_repulsion"]
exact = h2["hamiltonian"].ground_state_energy() + h2["nuclear_repulsion"]
print(f"VQE   {energy:.10f} Ha")
print(f"exact {exact:.10f} Ha (full CI of the basis: -1.137306)")
print("iterations:", result["iterations"], "energies of the restarts:", np.round(result["restart_energies"], 8))

VQE   -1.1373060358 Ha
exact -1.1373060358 Ha (full CI of the basis: -1.137306)
iterations: 12 energies of the restarts: [-1.85727503 -1.85727503 -1.85727503 -1.85727503]

plt.figure(figsize=(5, 3))
plt.semilogy(np.array(result["history"]) - h2["hamiltonian"].ground_state_energy() + 1e-16)
plt.xlabel("iteration")
plt.ylabel("energy error (Ha)")
plt.tight_layout()
plt.show()

Depth matters: a transverse-field Ising chain

Six spins at the critical point J = h = 1. The error of VQE falls with the number of ansatz layers.

h = q.transverse_field_ising(6, 1.0, 1.0)
exact = h.ground_state_energy()
layers = [1, 2, 3, 4, 5, 6]
errors = []
for L in layers:
    r = q.run_vqe(h, q.efficient_su2(6, L), restarts=4, seed=1)
    errors.append(r["energy"] - exact)
    print(f"{L} layers: E = {r['energy']:.8f}, error {errors[-1]:.2e}")

plt.figure(figsize=(5, 3))
plt.semilogy(layers, np.maximum(errors, 1e-12), "o-")
plt.xlabel("layers")
plt.ylabel("E_VQE − E_0")
plt.tight_layout()
plt.show()

1 layers: E = -7.23420165, error 6.20e-02
2 layers: E = -7.29130728, error 4.92e-03
3 layers: E = -7.29213077, error 4.10e-03
4 layers: E = -7.29593020, error 3.00e-04


5 layers: E = -7.29620836, error 2.14e-05
6 layers: E = -7.29622473, error 5.08e-06

QAOA for MaxCut

On a ring, depth-p QAOA reaches exactly (2p+1)/(2p+2) of the maximum cut. Layer-wise optimisation (INTERP) carries each optimum to the next depth.

graphs = {
    "ring C10": q.Graph.cycle(10),
    "Petersen": q.Graph.petersen(),
    "Metatron's Cube": q.Graph.metatron_cube(),
}
ratios = {}
for name, g in graphs.items():
    ratios[name] = [r["approximation_ratio"] for r in q.optimize_qaoa_layerwise(g, 4, seed=0)]
    print(name, np.round(ratios[name], 5))
print("ring, theory:", [round((2 * p + 1) / (2 * p + 2), 5) for p in range(1, 5)])

plt.figure(figsize=(5, 3))
for name, r in ratios.items():
    plt.plot(range(1, 5), r, "o-", label=name)
plt.xlabel("depth p")
plt.ylabel("approximation ratio")
plt.legend()
plt.tight_layout()
plt.show()

ring C10 [0.75    0.83333 0.875   0.9    ]
Petersen [0.86556 0.9178  0.97074 0.98906]


Metatron's Cube [0.99282 0.9991  0.97005 0.99419]
ring, theory: [0.75, 0.83333, 0.875, 0.9]

The most likely bit string of the optimised state is an optimal cut:

qaoa = q.MaxCutQaoa(q.Graph.petersen(), 4)
best = qaoa.optimize(restarts=8, seed=0)
print(best["most_likely_bitstring"], "cut", best["most_likely_cut"], "of", qaoa.max_cut)
print("probability of an optimal cut:", round(best["optimal_probability"], 4))

1001101101 cut 12.0 of 12.0
probability of an optimal cut: 0.9767