Hamiltonian-simulation examples
Realizing time evolution exp(-iHt) with two contrasting primitives: the
ancilla-based QSVT path and the ancilla-free Trotter product formula. Both
are checked against exact diagonalization.
See also
Trotter (product-formula simulation) and Qubitization and QSVT — the two primitives these examples contrast.
Hamiltonian simulation with PauliLCU and QSVT
Evolves a 4-qubit Pauli Hamiltonian by block-encoding it with PauliLCU and
applying a QSVT polynomial built from QSPPACK-generated phases. The
time-evolved state is checked against exact diagonalization.
"""Hamiltonian simulation with the PauliLCU block encoding and QSVT.
Evolves a 4-qubit Pauli Hamiltonian with QSPPACK-generated phases and checks
the result against exact diagonalization.
Requires qsppack and scipy. Run with:
python3 hamiltonian_simulation_qsvt.py
"""
import os
import numpy as np
import cudaq
from cudaq_algorithms import sim_utils as sim
from cudaq_algorithms import (PauliLCU, PhaseSequence, QSVT,
recover_real_time_evolution)
HAMILTONIAN = {
"ZIII": 0.70,
"IZII": -0.43,
"IIZI": 0.31,
"IIIZ": -0.22,
"XXII": 0.19,
"IYYI": -0.17,
"IZZX": 0.13,
"XYYX": 0.11,
}
TIME = 0.8
DEGREE = 16
def qsppack_phases(tau, degree):
"""cos/sin QSP phases for exp(-i tau x) via Jacobi-Anger + QSPPACK."""
try:
import qsppack
from scipy import special
except ImportError as exc:
raise SystemExit(
"This example needs qsppack and scipy: pip install qsppack scipy"
) from exc
cos_coefficients = np.array([0.5 * special.jv(0, tau)] +
[((-1)**k) * special.jv(2 * k, tau)
for k in range(1, degree // 2 + 1)])
sin_coefficients = np.array([((-1)**k) * special.jv(2 * k + 1, tau)
for k in range(degree // 2)])
options = {
"criteria": 1e-12,
"method": "Newton",
"typePhi": "full",
"useReal": True,
}
cos_phases, _ = qsppack.solve(cos_coefficients, 0, {
**options, "targetPre": True
})
sin_phases, _ = qsppack.solve(sin_coefficients, 1, {
**options, "targetPre": False
})
return [float(p) for p in cos_phases], [float(p) for p in sin_phases]
# Kept inline so the example runs standalone; the canonical copy of
# this dense reference lives in tests/python/dense_references.py.
def dense_matrix(terms, num_qubits):
dimension = 1 << num_qubits
matrix = np.zeros((dimension, dimension), dtype=np.complex128)
for word, coeff in terms.items():
for column in range(dimension):
row, phase = column, complex(coeff)
for qubit, label in enumerate(word):
bit = (column >> qubit) & 1
if label == "X":
row ^= 1 << qubit
elif label == "Y":
row ^= 1 << qubit
phase *= 1.0j if bit == 0 else -1.0j
elif label == "Z":
phase *= 1.0 if bit == 0 else -1.0
matrix[row, column] += phase
return matrix
def main():
cudaq.set_target(os.environ.get("CUDAQ_DEFAULT_SIMULATOR", "qpp-cpu"))
# -- quantum workflow -------------------------------------------------
encoding = PauliLCU(HAMILTONIAN)
transformer = QSVT(encoding)
tau = encoding.alpha * TIME
cos_phases, sin_phases = qsppack_phases(tau, DEGREE)
rng = np.random.default_rng(13)
psi = rng.normal(size=1 << encoding.num_system).astype(np.complex128)
psi /= np.linalg.norm(psi)
cos_state = sim.transform(transformer, psi,
PhaseSequence(cos_phases, convention="qsp"))
sin_state = sim.transform(transformer, psi,
PhaseSequence(sin_phases, convention="qsp"))
evolved = recover_real_time_evolution(cos_state, sin_state, cos_phases,
sin_phases)
# ----------------------------------------------------------------------
matrix = dense_matrix(HAMILTONIAN, encoding.num_system)
eigenvalues, eigenvectors = np.linalg.eigh(matrix)
exact = eigenvectors @ (np.exp(-1.0j * TIME * eigenvalues) *
(eigenvectors.conj().T @ psi))
l2_error = float(np.linalg.norm(evolved - exact))
fidelity = float(abs(np.vdot(exact, evolved))**2)
print("QSVT Hamiltonian simulation (pure-Python prototype)")
print("=" * 56)
print(f"encoding: {encoding}")
print(f"evolution time: {TIME}")
print(f"tau = alpha*t: {tau:.6f}")
print(f"QSPPACK degree: {DEGREE}")
print(f"L2 state error: {l2_error:.3e}")
print(f"fidelity: {fidelity:.12f}")
if l2_error > 1e-10:
raise SystemExit("evolution error exceeded 1e-10")
print("PASS")
return 0
if __name__ == "__main__":
raise SystemExit(main())
Suzuki-Trotter simulation of a chemistry-style Hamiltonian
Simulates a chemistry-style Hamiltonian, hard-coded as Pauli terms, with the
Trotter product-formula primitive. It shows two ways to run the same evolution
and the error-vs-steps trade-off, checked against dense linear algebra.
"""Suzuki-Trotter simulation of a chemistry-style Hamiltonian.
The Hamiltonian is hard-coded as Pauli terms to keep the example focused
on the algorithm primitives. Two ways to run the same evolution are shown:
1. ``sim_utils.evolve(evolution, ket, time, ...)`` — one call, identity
phase included (simulation helper).
2. A user kernel composing ``trotter.apply_trotter`` with state
preparation — the hardware-shaped composition path.
Run with: python3 trotter_chemistry.py
"""
import os
import numpy as np
import cudaq
from cudaq_algorithms import sim_utils, trotter
# A four-qubit molecular-style Pauli Hamiltonian. In a production chemistry
# workflow these terms would come from a fermion-to-qubit mapping.
HAMILTONIAN = {
"IIII": -0.81054798,
"ZIII": 0.17218393,
"IZII": -0.22575349,
"IIZI": 0.17218393,
"IIIZ": -0.22575349,
"ZZII": 0.12091263,
"ZIZI": 0.16892754,
"ZIIZ": 0.16614543,
"YYYY": 0.04523280,
"XXYY": 0.04523280,
"YYXX": 0.04523280,
"XXXX": 0.04523280,
"IZZI": 0.16614543,
"IZIZ": 0.17464343,
"IIZZ": 0.12091263,
}
TIME = 0.6
STEPS = 4
ORDER = 2
@cudaq.kernel
def prepare_state(q: cudaq.qview):
"""Small product-state superposition so non-commuting terms have
visible effect in the output amplitudes."""
ry(0.31, q[0])
rx(-0.27, q[1])
ry(0.19, q[2])
rx(0.23, q[3])
def pauli_matrix(word):
dim = 2**len(word)
matrix = np.zeros((dim, dim), dtype=np.complex128)
for basis in range(dim):
target_col = np.zeros(dim, dtype=np.complex128)
target_col[basis] = 1.0
result = np.zeros(dim, dtype=np.complex128)
for b, amplitude in enumerate(target_col):
if amplitude == 0.0:
continue
target, phase = b, 1.0 + 0.0j
for qubit, op in enumerate(word):
bit = (b >> qubit) & 1
if op == "X":
target ^= 1 << qubit
elif op == "Y":
target ^= 1 << qubit
phase *= -1.0j if bit else 1.0j
elif op == "Z":
phase *= -1.0 if bit else 1.0
result[target] += phase * amplitude
matrix[:, basis] = result
return matrix
def exact_evolve(evolution, time, ket):
matrix = evolution.identity_coefficient * np.eye(ket.size,
dtype=np.complex128)
for coefficient, word in zip(evolution.coefficients, evolution.words):
matrix += coefficient * pauli_matrix(str(word))
eigenvalues, eigenvectors = np.linalg.eigh(matrix)
return eigenvectors @ (np.exp(-1.0j * time * eigenvalues) *
(eigenvectors.conj().T @ ket))
def main():
cudaq.set_target(os.environ.get("CUDAQ_DEFAULT_SIMULATOR", "qpp-cpu"))
evolution = trotter.Trotter(
HAMILTONIAN,
ordering=trotter.TrotterOrdering.COEFFICIENT_MAGNITUDE_DESCENDING)
resources = evolution.resources(steps=STEPS, order=ORDER)
@cudaq.kernel
def prepare_only():
q = cudaq.qvector(4)
prepare_state(q)
ket0 = np.asarray(cudaq.get_state(prepare_only), dtype=np.complex128)
# Path 1: the one-call simulation helper (identity phase included).
evolved = sim_utils.evolve(evolution,
ket0,
time=TIME,
steps=STEPS,
order=ORDER)
exact = exact_evolve(evolution, TIME, ket0)
direct_error = float(np.linalg.norm(evolved - exact))
# Path 2: the escape hatch — compose apply_trotter in a user kernel.
coefficients = evolution.coefficients
words = [str(w) for w in evolution.words]
@cudaq.kernel
def evolve_kernel(coeffs: list[float], paulis: list[cudaq.pauli_word],
t: float, n_steps: int, formula_order: int):
q = cudaq.qvector(4)
prepare_state(q)
trotter.apply_trotter(coeffs, paulis, t, n_steps, formula_order, q)
kernel_state = np.asarray(cudaq.get_state(evolve_kernel, coefficients,
words, TIME, STEPS, ORDER),
dtype=np.complex128)
# The kernel path omits the identity phase; reintroduce it for comparison.
kernel_state = kernel_state * np.exp(
-1.0j * evolution.identity_coefficient * TIME)
paths_agree = float(np.linalg.norm(kernel_state - evolved))
print("Suzuki-Trotter chemistry-style example")
print("=" * 62)
print(f"num_qubits: {evolution.num_qubits}")
print(f"num_terms: {resources.num_terms}")
print(f"identity_coefficient: {evolution.identity_coefficient:+.8f}")
print(f"order: {resources.order}")
print(f"steps: {resources.steps}")
print(f"pauli_rotations: {resources.pauli_rotations}")
print(f"estimated_cx_count: {resources.estimated_cx_count}")
print(f"l2_error_vs_exact: {direct_error:.6e} (identity phase "
f"included -> no phase alignment needed)")
print(f"kernel_vs_evolve: {paths_agree:.6e}")
print("first four amplitudes:")
for idx, amplitude in enumerate(evolved[:4]):
print(f" |{idx:04b}> {amplitude.real:+.8f}{amplitude.imag:+.8f}j")
if direct_error > 5e-3 or paths_agree > 1e-12:
raise SystemExit("Trotter example exceeded tolerance")
print("PASS")
return 0
if __name__ == "__main__":
raise SystemExit(main())