How to do Hamiltonian simulation with GQSP in PennyLane

The Hamiltonian and the exact target

We will use a small two-qubit Heisenberg-type Hamiltonian. This is a convenient choice for this demo since it is already a sum of Pauli terms (the form Qubitization block-encodes directly), its terms do not all commute (so \(e^{-iHt}\) is non-trivial), and it is small enough that we can form the full matrix \(e^{-iHt}\) and compare against it directly. In general, though, the methods in this demo work for any Pauli Hamiltonian. \(\lambda=\sum_k|c_k|\) is the normalization used by the qubitization walk.

import warnings

warnings.filterwarnings("ignore", message=".*JAX.*")  # JAX is unused here
import numpy as np
import pennylane as qp
from scipy.linalg import expm
from scipy.special import jv  # Bessel functions of the first kind
import matplotlib.pyplot as plt

coeffs = [0.5, 0.3, 0.4, 0.2]
obs = [qp.X(0) @ qp.X(1), qp.Y(0) @ qp.Y(1), qp.Z(0) @ qp.Z(1), qp.Z(0)]
H = qp.Hamiltonian(coeffs, obs)
lam = sum(abs(c) for c in coeffs)  # l1-norm of the coefficients
t = 0.7  # evolution time

H_matrix = qp.matrix(H, wire_order=[0, 1])
U_exact = expm(-1j * H_matrix * t)

# show the Hamiltonian explicitly
print("H =", H)
print(f"{len(coeffs)} terms on 2 qubits   lambda = sum|c_k| = {lam}   t = {t}")
print("H as a matrix:")
print(np.round(H_matrix, 3))

H = 0.5 * (X(0) @ X(1)) + 0.3 * (Y(0) @ Y(1)) + 0.4 * (Z(0) @ Z(1)) + 0.2 * Z(0)
4 terms on 2 qubits   lambda = sum|c_k| = 1.4000000000000001   t = 0.7
H as a matrix:
[[ 0.6+0.j  0. +0.j  0. +0.j  0.2+0.j]
 [ 0. +0.j -0.2+0.j  0.8+0.j  0. +0.j]
 [ 0. +0.j  0.8+0.j -0.6+0.j  0. +0.j]
 [ 0.2+0.j  0. +0.j  0. +0.j  0.2+0.j]]

Step 1: The qubitization walk

qp.Qubitization(H, control) builds the Low and Chuang [2] qubitization walk operator

which applies a reflection about the control \(|0\rangle\) and then the block encoding of \(H/\lambda\). Its eigenphases are

for each eigenvalue \(E\) of \(H\), so a walk eigenvalue \(z = e^{i\theta}\) satisfies \(\cos\theta = E/\lambda\). The control register needs \(\lceil \log_2 L \rceil\) qubits for \(L\) Hamiltonian terms.

The printout below confirms this, since every \(\arccos(E/\lambda)\) appears among the walk’s eigenphases. The walk also carries a few extra phases (here \(0\) and \(\pi\)) coming from the complementary subspace where the control ancillas are not \(|0\rangle\). These encode no information about \(H\) and are projected out when we read off the top-left block in Step 3.

n_ctrl = int(np.ceil(np.log2(len(coeffs))))
anc = [f"a{i}" for i in range(n_ctrl)]


@qp.qnode(qp.device("default.qubit"))
def walk():
    qp.Qubitization(H, control=anc)
    return qp.state()


W = qp.matrix(walk, wire_order=anc + [0, 1])()
walk_phases = np.sort(np.unique(np.round(np.abs(np.angle(np.linalg.eigvals(W))), 4)))
arccos_E = np.sort(
    np.unique(np.round(np.arccos(np.clip(np.linalg.eigvalsh(H_matrix) / lam, -1, 1)), 4))
)
print("walk eigenphases   :", walk_phases)
print("arccos(E / lambda) :", arccos_E)

walk eigenphases   : [0.     1.0613 1.2626 1.487  2.6357 3.1416]
arccos(E / lambda) : [1.0613 1.2626 1.487  2.6357]

Step 2: The Jacobi-Anger polynomial

We need a polynomial \(P\) with \(P(e^{i\theta}) = e^{-i\lambda t\cos\theta}\). The Jacobi-Anger expansion provides exactly this, as the Laurent series in \(z=e^{i\theta}\):

with \(J_k\) being the Bessel functions of the first kind. The series converges super-exponentially once \(K \gtrsim a\), so we truncate at \(|k|\le K\), where \(K\) is the truncation order (the highest power of \(z\) we keep).

Two practical points that the GQSP machinery requires:

• poly_to_angles() needs a polynomial in non-negative powers, so we shift the Laurent series by \(z^{K}\) (we undo this shift in the circuit later).

• poly_to_angles() also requires \(|P(e^{i\theta})|\le 1\), so we rescale the coefficients by a constant \(s<1\) (and divide it back out at the end).

def jacobi_anger_poly(a, K):
    "shifted, rescaled Jacobi-Anger coefficients for exp(-i a cos theta), plus the scale s."
    laurent = {k: (-1j) ** k * jv(k, a) for k in range(-K, K + 1)}
    p = [laurent[j - K] for j in range(2 * K + 1)]  # shift to powers 0..2K
    grid = np.exp(1j * np.linspace(0, 2 * np.pi, 400))
    s = 0.99 / max(abs(np.polyval(p[::-1], z)) for z in grid)  # enforce |P| <= 1
    return [c * s for c in p], s


K = 8
poly, s = jacobi_anger_poly(lam * t, K)
print(f"K = {K}: polynomial degree {len(poly) - 1}, scale s = {s:.4f}")

K = 8: polynomial degree 16, scale s = 0.9900

Step 3: Run GQSP and recover \(e^{-iHt}\)

qp.GQSP(U, angles, control) applies the single polynomial \(P(U)\) with angles from qp.poly_to_angles(poly, "GQSP"). Recall how we built that polynomial in Step 2. The truncated Jacobi-Anger series already approximates \(e^{-iHt}\) on the walk eigenvalues, and we made two changes to it: multiplying by \(z^{K}\) to remove the negative powers, and scaling by \(s\) to enforce \(|P|\le 1\). So, the coefficients in poly are those of \(P(z) = s\, z^{K}\, [\text{Jacobi-Anger series}]\), and GQSP applies this single polynomial of \(W\) in one step, giving

The \(s\) and \(W^{K}\) terms are not separate circuits we multiply on afterwards. They are the two corrections already baked into the polynomial, which we now simply undo. We cancel the \(W^{K}\) factor by applying the adjoint walk \(W^\dagger\) a total of \(K\) times after the GQSP block, and we divide out \(s\) at read-out. The evolution then sits in the all-ancilla-zero block (the top-left \(\dim\times\dim\) corner), up to a global phase.

def gqsp_evolution(K):
    # recovered exp(-iH t) on the system, read from the GQSP top-left block.
    poly, s = jacobi_anger_poly(lam * t, K)
    angles = qp.poly_to_angles(poly, "GQSP")

    @qp.qnode(qp.device("default.qubit"))
    def circuit():
        qp.GQSP(qp.Qubitization(H, control=anc), angles, control="g")
        for _ in range(K):  # undo the z^K Laurent shift
            qp.adjoint(qp.Qubitization(H, control=anc))
        return qp.state()

    M = qp.matrix(circuit, wire_order=["g"] + anc + [0, 1])()
    block = M[:4, :4] / s  # all-ancilla-zero (top-left) block
    return block


block = gqsp_evolution(K)
ph = np.exp(-1j * np.angle(block[0, 0] / U_exact[0, 0]))  # match global phase
err = np.linalg.norm(block * ph - U_exact, 2)
print(f"||GQSP block - exp(-iHt)||_2 = {err:.2e}   (K = {K})")
print(
    "recovered block is unitary?  ||B^dag B - I|| =",
    f"{np.linalg.norm(block.conj().T @ block - np.eye(4), 2):.2e}",
)

||GQSP block - exp(-iHt)||_2 = 1.21e-08   (K = 8)
recovered block is unitary?  ||B^dag B - I|| = 7.72e-09

Convergence

The truncation error falls super-exponentially with the order \(K\) once \(K\gtrsim \lambda t\), which is the hallmark of the Jacobi-Anger approach. We plot the 2-norm distance to the exact evolution against \(K\).

def gqsp_error(K, t_val):
    "spectral-norm error of the GQSP evolution at order K and time t_val."
    poly, s = jacobi_anger_poly(lam * t_val, K)
    angles = qp.poly_to_angles(poly, "GQSP")

    @qp.qnode(qp.device("default.qubit"))
    def circuit():
        qp.GQSP(qp.Qubitization(H, control=anc), angles, control="g")
        for _ in range(K):  # undo the z^K Laurent shift
            qp.adjoint(qp.Qubitization(H, control=anc))
        return qp.state()

    M = qp.matrix(circuit, wire_order=["g"] + anc + [0, 1])()
    block = M[:4, :4] / s
    U = expm(-1j * H_matrix * t_val)
    ph = np.exp(-1j * np.angle(block[0, 0] / U[0, 0]))
    return np.linalg.norm(block * ph - U, 2)


Ks = list(range(2, 21))
lambda_t_targets = [1, 3, 5]  # show lambda*t = 1, 3, 5 exactly
t_values = [m / lam for m in lambda_t_targets]  # so lam * t = m

plt.style.use("pennylane.drawer.plot")
fig, ax = plt.subplots(figsize=(5.4, 3.4))
for t_val in t_values:
    errs = [gqsp_error(K, t_val) for K in Ks]
    a = lam * t_val
    ax.semilogy(Ks, errs, "o-", label=rf"$\lambda t = {round(a)}$")
    print(f"lambda t = {round(a)}:", [f"{e:.1e}" for e in errs])

ax.set_xticks(range(2, 21, 2))
ax.set_xlabel("Jacobi-Anger order $K$")
ax.set_ylabel(r"$\|U_{\mathrm{GQSP}} - e^{-iHt}\|_2$")
ax.set_title("GQSP error convergence")
ax.legend()
fig.tight_layout()
plt.show()

lambda t = 1: ['3.3e-02', '4.8e-03', '2.9e-04', '3.9e-05', '3.2e-06', '1.5e-07', '1.4e-08', '5.6e-10', '2.9e-11', '1.6e-11', '5.0e-11', '1.7e-11', '3.9e-11', '6.4e-11', '5.6e-11', '1.8e-11', '3.9e-11', '6.4e-11', '5.6e-11']
lambda t = 3: ['6.8e-01', '3.0e-01', '9.4e-02', '2.8e-02', '4.7e-03', '9.1e-04', '1.7e-04', '3.6e-05', '3.7e-06', '7.9e-07', '4.9e-08', '5.3e-09', '7.0e-10', '2.7e-11', '5.3e-12', '1.4e-11', '6.7e-12', '4.7e-11', '3.1e-11']
lambda t = 5: ['1.4e+00', '1.1e+00', '6.2e-01', '4.8e-01', '9.0e-02', '4.0e-02', '8.7e-03', '4.1e-03', '5.8e-04', '2.8e-04', '4.1e-05', '7.8e-06', '1.6e-06', '6.0e-08', '3.4e-08', '5.5e-09', '4.1e-10', '9.2e-11', '4.9e-12']

Figure: GQSP error convergence. Spectral-norm distance between the GQSP evolution and the exact \(e^{-iHt}\), versus the Jacobi-Anger order \(K\), for three evolution times \(\lambda t = 1, 3, 5\). Each curve falls super-exponentially once \(K\) exceeds its \(\lambda t\), then floors at machine precision; larger \(\lambda t\) needs a larger \(K\) to reach the same accuracy.

Discussion and conclusion

Running GQSP on the qubitization walk reproduced \(e^{-iHt}\) to within \(\sim 10^{-8}\) at truncation order \(K=8\) (a polynomial of degree \(2K=16\)), matching scipy.linalg.expm to machine precision. Because of the fast Bessel decay seen in the convergence plot, useful accuracy needs only a truncation order \(K = \mathcal{O}(\lambda t + \log(1/\varepsilon))\), where \(\varepsilon\) is the target error. In other words, the number of walk applications grows only mildly as the accuracy is tightened.

There are two practical things to keep in mind:

• One control qubit, one polynomial. GQSP needs a single ancilla and a single complex polynomial, applied directly rather than as separate real and imaginary parts, which is what makes its angle synthesis simpler and more stable than ordinary QSP for this task.

• Where the cost lives. The accuracy knob is the truncation order \(K\); the polynomial GQSP applies has degree \(2K\), so the circuit uses \(2K\) applications of the walk \(W\) (plus \(K\) adjoint-walk applications to undo the \(z^{K}\) shift). The depth therefore grows quickly with \(K\). The resource-estimation demo is a good reference to explore this cost.

In practice, this makes GQSP a natural choice when you need an accurate \(e^{-iHt}\) over a longer evolution time \(t\). The super-exponential convergence means the order \(K\) (and hence the number of walk applications) grows only mildly as you tighten the error, so the cost stays close to the \(\lambda t\) set by the evolution itself. This points to an accuracy-versus-resource trade-off rather than one method dominating: GQSP (and block-encoding methods generally) reach high accuracy with a query count that grows only mildly as the error shrinks, while product-formula (Trotter) approaches are often cheaper in qubit and gate counts for a coarse target but need rapidly more steps to reach the same accuracy. Within the block-encoding / quantum-signal-processing family, GQSP’s particular appeal is economy of ancillas, since a single control qubit carries one complex polynomial applied directly, instead of building the polynomial’s real and imaginary parts separately as ordinary QSP must.

Now it is your turn: swap in your own Pauli Hamiltonian, pick an evolution time, and run the same recipe to obtain \(e^{-iHt}\). To dig deeper into the building blocks, see the introduction to qubitization, QSVT in practice, and the resource-estimation demo.

References