Hamiltonian simulation

Hamiltonian simulation builds a circuit for exp(−iHt), with H = Σ_γ h_γ P_γ a sum of Pauli strings. QSO does this two ways: product formulas, which come with rigorous error bounds, and qDRIFT, a randomised method. Both use Circuit.pauli_rotation, which appends exp(−iθP/2) for any Pauli string P. One-qubit strings and XX, YY, ZZ become a single gate. Longer strings become a basis change, a CNOT ladder and one rz.

Product formulas

trotter_circuit(h, n_qubits, t, steps, order) splits the time into steps steps of length τ = t/steps:

trotter_error_bound returns an upper bound on the spectral-norm distance between the circuit and exp(−iHt). The commutators in the bound are computed exactly in the Pauli algebra:

• order 1: τ²/2 · Σ_γ ‖[Σ_{γ'>γ} H_γ', H_γ]‖ (Childs, Su, Tran, Wiebe and Zhu, Phys. Rev. X 11, 011020 (2021), Proposition 15);
• order 2: the nested-commutator bound of their Proposition 16;
• order 4 and higher: the smaller of a stage-by-stage second-order bound and a Taylor-remainder bound.

Each bound is multiplied by the number of steps.

import numpy as np
import metatron_qso as q

h = q.heisenberg_xxz(4, 1.0, 0.5)
w, v = np.linalg.eigh(h.to_matrix(4))
exact = (v * np.exp(-1j * w)) @ v.conj().T          # exp(-iH), t = 1

def unitary(c):
    cols = []
    for k in range(16):
        prep = q.Circuit(4)
        for b in range(4):
            if k >> b & 1:
                prep.x(b)
        cols.append(q.simulate(prep.append(c)).amplitudes)
    return np.array(cols).T

for order in (1, 2, 4):
    c = q.trotter_circuit(h, 4, 1.0, steps=10, order=order)
    error = np.linalg.norm(unitary(c) - exact, 2)
    bound = q.trotter_error_bound(h, 1.0, steps=10, order=order)
    assert error <= bound
    print(order, len(c), f"{error:.2e} <= {bound:.2e}")

The identity term of H only adds a global phase, so circuits leave it out. In the tests (Heisenberg, Ising and random four-qubit Hamiltonians), the first- and second-order bounds are 3.5 to 13 times the actual error. The bound for order 4 and higher is valid but loose: it shrinks only as r⁻² with the number of steps r until the steps are short, while the actual error shrinks as r⁻⁴. Doubling the steps divides the actual error by 2^order, and the tests check this.

qDRIFT

qdrift_circuit(h, n_qubits, t, samples, seed) draws N = samples terms at random, with probability |h_γ|/λ (λ = Σ|h_γ|). Each drawn term is applied for time λt/N (Campbell, Phys. Rev. Lett. 123, 070503 (2019)). The average over seeds is a channel within diamond distance qdrift_error_bound = 2λ²t²/N · e^{2λt/N} of the exact evolution. The tests compute that channel exactly. Its error falls as 1/N and stays below the bound.

import metatron_qso as q

h = q.PauliSum("0.8 X0 X1 + 0.5 Z1 + 0.3 Y1 Y2 - 0.6 Z0 Z2")
c = q.qdrift_circuit(h, 3, 1.0, samples=200, seed=7)
print(len(c), q.qdrift_error_bound(h, 1.0, 200))

qDRIFT pays off when H has many small terms, because its cost depends on λ and not on the number of terms.

On the command line

$ qso trotter --hamiltonian "Z0 Z1 + 0.5 X0 + 0.5 X1" --time 1 --steps 10 --order 2
order 2, 10 steps, t = 1: 41 gates, depth 31, error bound 5.000e-3, largest column error 1.938e-3 (a lower bound on ‖U − exp(−iHt)‖)

--output FILE writes the circuit as OpenQASM 2.0.