Getting-started walkthrough
An ordered learning path through the library. Each example builds on the one before it, from the single core idea – block-encode a Hamiltonian and walk it – up to bringing your own encoding. Every example verifies its own claims against an independent dense reference, so each script runs to completion or fails loudly.
See also
The user guide pages cover each primitive in depth; Getting Started sets up the installation these scripts assume.
Quickstart: block-encode a Hamiltonian and walk it
The five-minute tour of the core idea. A Hamiltonian is not a circuit, so you
block-encode it: PauliLCU hides H / alpha inside a larger unitary on a few
ancilla qubits, and the qubitization Walk turns powers of that unitary into
Chebyshev moments checked against a dense matrix.
"""Example 1 — Quickstart: block-encode a Hamiltonian and walk it.
The five-minute tour of the core idea. A Hamiltonian H is not a circuit
(it is not unitary), so you *block-encode* it: PauliLCU hides H / alpha
inside a bigger unitary on a few ancilla qubits. The qubitization Walk
then turns powers of that unitary into Chebyshev polynomials of H --
measured as expectation values, the raw data spectral algorithms consume.
Everything downstream (examples 2-7) is built on these three objects:
PauliLCU (the block encoding), Walk (qubitization), and the observable
`moment` path. Here we build them and check the moments against a dense
matrix, in the house style: every claim verified against an independent
reference.
Run: python3 01_quickstart_block_encoding.py
"""
from __future__ import annotations
import os
import cudaq
import numpy as np
from cudaq_algorithms import PauliLCU, Walk
# The convention: read docs/sphinx/conventions.rst. The walk returns +<T_k(H/alpha)>.
_PAULIS = {
"I": np.eye(2, dtype=complex),
"X": np.array([[0, 1], [1, 0]], dtype=complex),
"Y": np.array([[0, -1j], [1j, 0]]),
"Z": np.diag([1.0, -1.0]).astype(complex),
}
def dense_hamiltonian(terms: dict) -> np.ndarray:
"""Dense matrix of a Pauli sum (little-endian: qubit 0 is the low bit,
so word[0] is the rightmost tensor factor -- matching CUDA-Q)."""
matrix = 0
for word, coefficient in terms.items():
factor = np.array([[1]], dtype=complex)
for label in word:
factor = np.kron(_PAULIS[label], factor)
matrix = matrix + coefficient * factor
return matrix
def chebyshev(x: np.ndarray, k: int) -> np.ndarray:
"""T_k(x) for a matrix argument, by the Chebyshev recurrence."""
if k == 0:
return np.eye(x.shape[0], dtype=complex)
previous, current = np.eye(x.shape[0], dtype=complex), x
for _ in range(2, k + 1):
previous, current = current, 2 * x @ current - previous
return current
def main() -> int:
cudaq.set_target(os.environ.get("CUDAQ_DEFAULT_SIMULATOR", "qpp-cpu"))
# 1. Block-encode a Hamiltonian given as {pauli_word: coefficient}.
terms = {"ZZ": 0.5, "XI": 0.3, "IX": 0.3}
encoding = PauliLCU(terms)
print(f"system qubits : {encoding.num_system}")
print(f"ancillas : {encoding.num_ancilla}")
print(f"alpha (1-norm): {encoding.alpha:.6f} # W encodes H / alpha")
# 2. Build the qubitization walk over that encoding.
walk = Walk(encoding)
# 3. Measure Chebyshev moments <T_k(H/alpha)> from a chosen state.
rng = np.random.default_rng(0)
state = rng.normal(size=1 << encoding.num_system).astype(np.complex128)
state /= np.linalg.norm(state)
moments = walk.moments(state, 4)
# 4. Verify against the dense matrix -- the numbers must agree.
x = dense_hamiltonian(terms) / encoding.alpha
print("\n k walk.moment dense <T_k(H/alpha)>")
for k, measured in enumerate(moments):
reference = float(np.real(state.conj() @ chebyshev(x, k) @ state))
print(f" {k} {measured:+.8f} {reference:+.8f}")
assert abs(measured - reference) < 1e-10
print("\nOK — the walk reproduces Chebyshev polynomials of H to machine "
"precision.")
return 0
if __name__ == "__main__":
raise SystemExit(main())
Hamiltonian simulation, two independent ways
Compute the time evolution exp(-iHt)|psi> with two different primitives and
check both against dense linear algebra: QSVT with Jacobi-Anger phase factors,
and Trotter product formulas with an error-vs-steps trade-off. Because the two
constructions are completely different, their agreement is a strong test.
"""Example 2 — Hamiltonian simulation, two independent ways.
Compute exp(-i H t)|psi> with two different primitives and check both
against dense linear algebra:
* QSVT -- block-encode H, apply a polynomial approximation of exp(-ixt)
built from Jacobi-Anger phase factors (generated by QSPPACK).
* Trotter -- a product formula, no block encoding or ancillas, with the
error-vs-steps trade-off and a resource estimate.
Because the two constructions are completely different, agreeing with the
same dense reference means they cross-validate each other. This is the
"primitives compose into applications" story made concrete.
Prerequisite: qsppack (installed with the wheel). Run: python3 02_hamiltonian_simulation.py
"""
from __future__ import annotations
import os
import cudaq
import numpy as np
from scipy.linalg import expm
from cudaq_algorithms import (PauliLCU, PhaseSequence, QSVT, Trotter,
recover_real_time_evolution)
from cudaq_algorithms import sim_utils as sim
_PAULIS = {
"I": np.eye(2, dtype=complex),
"X": np.array([[0, 1], [1, 0]], dtype=complex),
"Z": np.diag([1.0, -1.0]).astype(complex),
}
def dense(terms, num_qubits):
# Little-endian, matching CUDA-Q: word[i] acts on qubit i, and qubit 0
# is the least-significant bit, so build P[word[0]] as the rightmost
# tensor factor by kron-ing each successive Pauli on the left.
matrix = np.zeros((1 << num_qubits, 1 << num_qubits), dtype=complex)
for word, coefficient in terms.items():
factor = np.array([[1]], dtype=complex)
for label in word:
factor = np.kron(_PAULIS[label], factor)
matrix = matrix + coefficient * factor
return matrix
def jacobi_anger_phases(tau, degree):
"""cos/sin QSP phases for exp(-i tau x) via Jacobi-Anger + QSPPACK."""
import contextlib
import io
import qsppack
from scipy import special
cos_coeffs = np.array([0.5 * special.jv(0, tau)] +
[((-1)**k) * special.jv(2 * k, tau)
for k in range(1, degree // 2 + 1)])
sin_coeffs = 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
}
with contextlib.redirect_stdout(io.StringIO()): # hush the QSP solver
cos_phases, _ = qsppack.solve(cos_coeffs, 0, {
**options, "targetPre": True
})
sin_phases, _ = qsppack.solve(sin_coeffs, 1, {
**options, "targetPre": False
})
return [float(p) for p in cos_phases], [float(p) for p in sin_phases]
def main() -> int:
cudaq.set_target(os.environ.get("CUDAQ_DEFAULT_SIMULATOR", "qpp-cpu"))
terms = {"ZZ": 0.7, "XI": 0.4, "IX": 0.4, "ZI": 0.31}
num_qubits = 2
time = 0.8
matrix = dense(terms, num_qubits)
rng = np.random.default_rng(1)
psi = rng.normal(size=1 << num_qubits).astype(np.complex128)
psi /= np.linalg.norm(psi)
exact = expm(-1j * matrix * time) @ psi
# -- QSVT --------------------------------------------------------------
encoding = PauliLCU(terms)
transformer = QSVT(encoding)
tau = encoding.alpha * time
cos_phases, sin_phases = jacobi_anger_phases(tau, degree=16)
cos_state = sim.transform(transformer, psi,
PhaseSequence(cos_phases, convention="qsp"))
sin_state = sim.transform(transformer, psi,
PhaseSequence(sin_phases, convention="qsp"))
qsvt_state = recover_real_time_evolution(cos_state, sin_state, cos_phases,
sin_phases)
qsvt_error = float(np.max(np.abs(qsvt_state - exact)))
print(f"QSVT (degree 16) : max amplitude error {qsvt_error:.2e}")
# -- Trotter -----------------------------------------------------------
evolution = Trotter(terms)
print("Trotter (order 2) : error shrinks as steps^-2")
for steps in (1, 2, 4, 8, 16):
trotter_state = sim.evolve(evolution, psi, time, steps=steps, order=2)
error = float(np.max(np.abs(trotter_state - exact)))
print(f" steps={steps:2d} max amplitude error {error:.2e}")
resources = evolution.resources(steps=16, order=2)
print(f" resource estimate at 16 steps: "
f"{resources.pauli_rotations} Pauli rotations")
assert qsvt_error < 1e-6
print("\nOK — two different circuits, one physics: both reproduce "
"exp(-iHt)|psi>.")
return 0
if __name__ == "__main__":
raise SystemExit(main())
From a molecule to its ground-state energy
The full classical to quantum to classical loop, end to end: PySCF mean field,
Jordan-Wigner qubit Hamiltonian, PauliLCU block encoding, Walk moments, and
a classical Krylov solve for the ground-state energy checked against FCI.
"""Example 3 — From a molecule to its ground-state energy.
The full classical -> quantum -> classical loop, end to end:
PySCF mean field
-> chemistry.from_pyscf (chemist (pq|rs) MO integrals)
-> chemistry.qubit_hamiltonian (Jordan-Wigner -> SpinOperator)
-> PauliLCU (block-encode H / alpha)
-> Walk.moments (Chebyshev moments <T_k(H/alpha)>)
-> classical Krylov solve (quantum exact Lanczos)
-> ground-state energy, checked against FCI.
This is the example to read first if you care about chemistry: it shows
how integrals become a qubit Hamiltonian, how the qubitization walk turns
into measured moments, and how those moments feed a real ground-state
algorithm. Every number is checked against an independent reference (FCI).
Reference: Kirby, Motta, Mezzacapo, "Exact and efficient Lanczos method
on a quantum computer", Quantum 7, 1018 (2023), arXiv:2208.00567.
Prerequisite: PySCF (`pip install pyscf`). Run: python3 03_chemistry_to_ground_state.py
"""
from __future__ import annotations
import os
import cudaq
import numpy as np
from cudaq_algorithms import PauliLCU, Walk, chemistry
def krylov_matrices(moments: np.ndarray, dimension: int):
"""Overlap S and scaled-Hamiltonian H~ from Chebyshev moments.
The Krylov basis is |phi_i> = T_i(H/alpha)|ref>. With
mu_k = <ref|T_k(H/alpha)|ref> and the Chebyshev product identities
T_i T_j = (T_{i+j} + T_{|i-j|}) / 2 and x T_j = (T_{j+1} + T_{|j-1|}) / 2,
both matrices are just combinations of the measured moments.
"""
mu = moments
overlap = np.empty((dimension, dimension))
scaled = np.empty((dimension, dimension))
for i in range(dimension):
for j in range(dimension):
overlap[i, j] = 0.5 * (mu[i + j] + mu[abs(i - j)])
scaled[i,
j] = 0.25 * (mu[i + j + 1] + mu[abs(i - j - 1)] +
mu[i + abs(j - 1)] + mu[abs(i - abs(j - 1))])
return overlap, scaled
def ground_eigenvalue(scaled, overlap, cutoff=1e-8):
"""Solve the projected generalized eigenproblem, dropping null directions."""
values, vectors = np.linalg.eigh(overlap)
keep = values > cutoff
transform = vectors[:, keep] / np.sqrt(values[keep])
projected = transform.T @ scaled @ transform
return float(np.linalg.eigvalsh(projected).min())
def main() -> int:
cudaq.set_target(os.environ.get("CUDAQ_DEFAULT_SIMULATOR", "qpp-cpu"))
try:
from pyscf import fci, gto, scf
except ImportError:
print("This example needs PySCF: pip install pyscf")
return 0
# 1. Converged restricted mean field for H2 / STO-3G.
geometry = [("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7474))]
mol = gto.M(atom=geometry, basis="sto-3g", symmetry=False)
mean_field = scf.RHF(mol).run(verbose=0)
# 2. Extract chemist-notation MO integrals + nuclear repulsion.
one_body, eri, nuclear_repulsion = chemistry.from_pyscf(mean_field)
print(f"molecule : H2 / STO-3G ({one_body.shape[0]} spatial "
f"orbitals, {2 * one_body.shape[0]} qubits)")
# 3. Jordan-Wigner qubit Hamiltonian (nuclear repulsion kept classical).
hamiltonian = chemistry.qubit_hamiltonian(one_body, eri, scalar_offset=0.0)
# 4. Block-encode it and build the qubitization walk.
encoding = PauliLCU(hamiltonian)
walk = Walk(encoding)
print(f"block encoding : {hamiltonian.term_count} Pauli terms, "
f"alpha = {encoding.alpha:.6f}, {encoding.num_ancilla} ancillas")
# 5. Measure Chebyshev moments from the Hartree-Fock reference.
# (HF determinant: the lowest `num_electrons` qubits occupied.)
reference = np.zeros(1 << encoding.num_system, dtype=np.complex128)
reference[(1 << mol.nelectron) - 1] = 1.0
krylov_dimension = 4
moments = np.asarray(walk.moments(reference, 2 * krylov_dimension))
print(f"Chebyshev moments : {np.array2string(moments, precision=4)}")
# 6. Classical Krylov solve -> QEL energy (add nuclear repulsion back).
overlap, scaled = krylov_matrices(moments, krylov_dimension)
qel_energy = ground_eigenvalue(scaled, overlap) * encoding.alpha \
+ nuclear_repulsion
# 7. Check against FCI.
fci_energy = float(fci.FCI(mean_field).kernel()[0])
print(f"\nQEL ground energy : {qel_energy:.10f} Ha")
print(f"FCI reference : {fci_energy:.10f} Ha")
print(f"|error| : {abs(qel_energy - fci_energy):.2e} Ha")
assert abs(qel_energy - fci_energy) < 1e-8
print("\nOK — the qubitization moments reproduce the exact ground state.")
return 0
if __name__ == "__main__":
raise SystemExit(main())
Double factorization and the BlockEncoding protocol
Double factorization rewrites the two-electron tensor as a sum of low-rank
leaves, giving a compression dial you can watch trade accuracy for a cheaper
block encoding. The same pipeline scales to larger molecules and plugs into the
structural BlockEncoding protocol shared by every consumer.
"""Example 4 — Double factorization and the BlockEncoding protocol.
Three ideas, built on real molecules from PySCF:
A. The compression dial (H2). Double factorization rewrites the
two-electron tensor as a sum of low-rank "leaves". Keep fewer leaves
and you trade accuracy for a cheaper block encoding -- a knob you can
watch, checked against the exact FCI energy at every setting.
B. It scales, and DF is compact (N2, 20 qubits). The same pipeline runs
at 20 qubits. There, a flat Pauli LCU block-encodes the Hamiltonian
as a SELECT over ~3000 Pauli strings; double factorization encodes
the *same* operator with a few dozen structured leaves. Far fewer
building blocks, comparable normalization.
C. Polymorphism through the protocol. `PauliLCU` and
`DoubleFactorizedEncoding` are unrelated classes, but both satisfy
`BlockEncoding` -- so `Walk` and `QSVT` consume either with identical
code.
Prerequisite: PySCF (`pip install pyscf`). Takes a few seconds (the N2
FCI reference is the slow part). Run: python3 04_double_factorization_and_the_protocol.py
"""
from __future__ import annotations
import os
import pathlib
import sys
import cudaq
from cudaq_algorithms import (BlockEncoding, PauliLCU, QSVT, Walk, chemistry)
from cudaq_algorithms import double_factorization as df
# The double-factorized encoding is the worked BlockEncoding example that
# lives beside this script (CUDA-Q kernels need real .py files, so a plain
# path-based import works).
sys.path.insert(0, str(pathlib.Path(__file__).resolve().parent))
from df_encoding import DoubleFactorizedEncoding
def mean_field(atom: str, basis: str = "sto-3g"):
from pyscf import gto, scf
mol = gto.M(atom=atom, basis=basis, symmetry=False)
return mol, scf.RHF(mol).run(verbose=0)
def fci_energy(one_body, eri, num_electrons, nuclear_repulsion) -> float:
"""Exact (FCI) total energy of the given MO integrals, via PySCF."""
from pyscf import fci
electronic, _ = fci.direct_spin1.kernel(one_body, eri, one_body.shape[0],
num_electrons)
return float(electronic) + nuclear_repulsion
def main() -> int:
cudaq.set_target(os.environ.get("CUDAQ_DEFAULT_SIMULATOR", "qpp-cpu"))
try:
import pyscf # noqa: F401
except ImportError:
print("This example needs PySCF: pip install pyscf")
return 0
# -- A. The compression dial (H2) ------------------------------------
mol, mf = mean_field("H 0 0 0; H 0 0 0.7414")
one_body, eri, e_nuc = chemistry.from_pyscf(mf)
exact = fci_energy(one_body, eri, mol.nelectron, e_nuc)
full = df.explicit_double_factorization(eri, threshold=0.0)
print(f"A. H2/STO-3G compression dial (exact FCI = {exact:.6f} Ha)")
print(f" {'leaves':>6} {'tensor error':>13} {'alpha':>8} "
f"{'E error (Ha)':>13}")
for num_leaves in range(1, full.num_leaves + 1):
truncated = df.explicit_double_factorization(eri,
max_num_leaves=num_leaves)
reconstructed_eri = df.reconstruct_eri(truncated)
alpha = PauliLCU(
chemistry.qubit_hamiltonian(one_body,
reconstructed_eri,
scalar_offset=e_nuc)).alpha
energy = fci_energy(one_body, reconstructed_eri, mol.nelectron, e_nuc)
print(
f" {num_leaves:>6} {df.factorization_error(eri, truncated):>13.2e}"
f" {alpha:>8.4f} {energy - exact:>+13.2e}")
# -- B. It scales, and DF is compact (N2, 20 qubits) -----------------
mol, mf = mean_field("N 0 0 0; N 0 0 1.09")
one_body, eri, e_nuc = chemistry.from_pyscf(mf)
hamiltonian = chemistry.qubit_hamiltonian(one_body,
eri,
scalar_offset=e_nuc)
flat = PauliLCU(hamiltonian)
factorization = df.explicit_double_factorization(eri, threshold=0.0)
dfe = DoubleFactorizedEncoding(one_body,
factorization,
scalar_offset=e_nuc)
print(f"\nB. N2/STO-3G on {hamiltonian.qubit_count} qubits: "
f"same Hamiltonian, two block encodings")
print(f" flat PauliLCU : alpha = {flat.alpha:6.2f} over "
f"{hamiltonian.term_count} Pauli terms")
print(f" double factn : alpha = {dfe.alpha:6.2f} over "
f"{factorization.num_leaves} leaves")
print(" -> DF encodes the same operator with ~50x fewer SELECT building "
"blocks,\n at a comparable one-norm (in a minimal basis the DF "
"rank is near\n full, so alpha barely drops; larger bases "
"compress far more).")
# -- C. One protocol, two encodings ----------------------------------
print(
"\nC. both are BlockEncodings, so Walk/QSVT consume either unchanged")
for name, encoding in [("PauliLCU", flat),
("DoubleFactorizedEncoding", dfe)]:
assert isinstance(encoding, BlockEncoding)
Walk(encoding)
QSVT(encoding) # the exact same construction works on both
print(f" {name:26s} system={encoding.num_system} "
f"ancilla={encoding.num_ancilla} -> Walk + QSVT built")
print("\nOK — pyscf integrals in, block encodings out; DF is the compact, "
"scalable one, and both plug into the same primitives (example 6 "
"shows how to add your own).")
return 0
if __name__ == "__main__":
raise SystemExit(main())
State preparation and injection
Every primitive factory (Walk, QSVT, Trotter) accepts a state_prep
kernel argument. Pass one and the factory returns a zero-argument, hardware-shaped
circuit with no statevector crossing the API boundary – the seam that makes the primitives
hardware-ready and the seam a tensor-network state-prep compiler would plug into.
"""Example 5 — State preparation and injection.
Every primitive factory (Walk, QSVT, Trotter) takes a `state_prep`
argument: a `(qubits: qview)` CUDA-Q kernel that prepares the register.
Pass one and the factory returns a *zero-argument* circuit -- no
`cudaq.State` argument, no statevector crossing the API boundary -- that
you can sample directly
and hand to synthesis. This is the seam that makes the primitives
hardware-ready, and the seam an MPS/tensor-network state-prep compiler
would plug into.
Two things this shows:
A. The injection contract: an injected walk is a zero-argument kernel,
and its measured moment matches the data-path (cudaq.State) form.
B. Any `(qubits: qview)` kernel qualifies -- including a Givens
Slater-determinant prep produced from an orbital-coefficient matrix.
Run: python3 05_state_prep_and_injection.py
"""
from __future__ import annotations
import os
import cudaq
import numpy as np
from cudaq_algorithms import PauliLCU, Walk, stateprep
# A state-prep kernel is just a `(qubits: qview)` device kernel. This one
# prepares |...01> -- one electron in the lowest orbital (a 2-qubit
# Hartree-Fock reference).
@cudaq.kernel
def hartree_fock_prep(qubits: cudaq.qview):
x(qubits[0])
def main() -> int:
cudaq.set_target(os.environ.get("CUDAQ_DEFAULT_SIMULATOR", "qpp-cpu"))
encoding = PauliLCU({"ZZ": 0.5, "XI": 0.3, "IX": 0.3, "ZI": 0.2})
walk = Walk(encoding)
# -- A. Injection contract -------------------------------------------
# With state_prep, walk.kernel returns a ZERO-ARGUMENT circuit: prepare
# the reference, then apply W^3. Sample it like any kernel.
injected = walk.kernel(power=3, state_prep=hartree_fock_prep)
counts = cudaq.sample(injected, shots_count=2000)
print("A. injected walk is a zero-argument, directly sampleable circuit")
print(f" sampled {len(counts)} bitstrings, e.g. "
f"{dict(list(counts.items())[:3])}")
# The observable path also accepts state_prep. It must agree with the
# data path, where we hand in the same reference as a statevector.
reference = np.zeros(1 << encoding.num_system, dtype=np.complex128)
reference[1] = 1.0 # |01>, matching x(qubits[0])
moment_injected = walk.moment(None, 3, state_prep=hartree_fock_prep)
moment_data = walk.moment(reference, 3)
print(f" moment via injection : {moment_injected:+.10f}")
print(f" moment via data path : {moment_data:+.10f}")
assert abs(moment_injected - moment_data) < 1e-10
print(" -> identical: injection changes how the state enters, not the "
"physics")
# -- B. Any (qubits: qview) kernel qualifies -------------------------
# A Givens Slater determinant of an orbital-coefficient matrix Q is
# also just a (qubits: qview) kernel -- built from Q, injectable the
# same way. Its computational-basis amplitudes are the minors of Q.
orbital_coefficients = np.array([[0.6, 0.0], [0.8, 0.0], [0.0, 0.6],
[0.0, 0.8]])
schedule = stateprep.make_givens_rotation_schedule(orbital_coefficients)
slater_prep = stateprep.slater_determinant_kernel(schedule)
@cudaq.kernel
def prepare_slater():
qubits = cudaq.qvector(4)
slater_prep(qubits)
slater_counts = cudaq.sample(prepare_slater, shots_count=4000)
print("\nB. a Givens Slater-determinant prep is the same kind of kernel")
print(f" 4 spin-orbitals, 2 electrons -> occupied-pair bitstrings:")
for bits, n in sorted(slater_counts.items(), key=lambda kv: -kv[1])[:4]:
print(f" |{bits}> {n/4000:.3f}")
print(" (these frequencies are |det Q[S,:]|^2 -- an entangled "
"reference, injectable exactly like the HF one)")
return 0
if __name__ == "__main__":
raise SystemExit(main())
Bring your own block encoding
Walk and QSVT are generic over the BlockEncoding protocol: any object
exposing the right members works, with no inheritance. This example implements
the protocol from scratch for a two-term single-qubit LCU and watches the walk
measure correct Chebyshev moments on it.
"""Example 6 — Bring your own block encoding.
`Walk` and `QSVT` are generic over the `BlockEncoding` protocol: any
object exposing the right members works, with no inheritance and no
changes to the primitives. This example implements the protocol from
scratch for the simplest nontrivial case -- a two-term LCU of a single
qubit, H = a*X + b*Z -- and shows the qubitization `Walk` measuring
correct Chebyshev moments on it.
The protocol (see docs and `block_encoding.py`): three sizes
(`num_system`, `num_ancilla`, `alpha`) plus kernel factories. The walk
step convention, mirrored below, is: SELECT, then reflect about the
PREPARE state. Get that right and everything else is inherited.
Run: python3 06_bring_your_own_encoding.py
"""
from __future__ import annotations
import math
import os
import cudaq
import numpy as np
from cudaq_algorithms import BlockEncoding, Walk
from cudaq_algorithms.common_kernels import reflect_about_zero
class TwoTermLCU:
"""Block encoding of H = a*X + b*Z (a, b > 0) on one system qubit.
PREPARE puts the single ancilla into (sqrt(a), sqrt(b)) / sqrt(alpha);
SELECT applies X when the ancilla is |0> and Z when it is |1>. Then
<0|B^dag SELECT B|0> = (a*X + b*Z) / alpha = H / alpha.
"""
def __init__(self, a: float, b: float):
self.a, self.b = float(a), float(b)
self.num_system = 1
self.num_ancilla = 1
self.alpha = self.a + self.b
# PREPARE angle: ry(theta)|0> = (cos, sin) = (sqrt(a), sqrt(b))/sqrt(alpha)
self._theta = 2.0 * math.atan2(math.sqrt(self.b), math.sqrt(self.a))
# -- sizes -----------------------------------------------------------
# num_system, num_ancilla, alpha are plain attributes (set above).
# -- kernel factories (data captured at factory time) ----------------
def prepare_kernel(self):
theta = self._theta
@cudaq.kernel
def prepare(ancilla: cudaq.qview):
ry(theta, ancilla[0])
return prepare
def unprepare_kernel(self):
theta = self._theta
@cudaq.kernel
def unprepare(ancilla: cudaq.qview):
ry(-theta, ancilla[0])
return unprepare
def _select_kernel(self):
@cudaq.kernel
def select(ancilla: cudaq.qview, system: cudaq.qview):
# |0><0| (x) X : X on system, controlled on ancilla == 0
x(ancilla[0])
x.ctrl(ancilla[0], system[0])
x(ancilla[0])
# |1><1| (x) Z : Z on system, controlled on ancilla == 1
z.ctrl(ancilla[0], system[0])
return select
def walk_step_kernel(self):
theta = self._theta
select = self._select_kernel()
@cudaq.kernel
def walk_step(ancilla: cudaq.qview, system: cudaq.qview):
select(ancilla, system) # SELECT
ry(-theta, ancilla[0]) # reflect about PREPARE state:
reflect_about_zero(ancilla) # B (I - 2|0><0|) B^dag
ry(theta, ancilla[0])
return walk_step
def adjoint_walk_step_kernel(self):
theta = self._theta
select = self._select_kernel()
@cudaq.kernel
def adjoint_walk_step(ancilla: cudaq.qview, system: cudaq.qview):
ry(-theta, ancilla[0]) # reflection first (self-adjoint),
reflect_about_zero(ancilla)
ry(theta, ancilla[0])
select(ancilla, system) # then SELECT (self-adjoint)
return adjoint_walk_step
def apply_kernel(self):
theta = self._theta
select = self._select_kernel()
@cudaq.kernel
def apply(ancilla: cudaq.qview, system: cudaq.qview):
ry(theta, ancilla[0])
select(ancilla, system)
ry(-theta, ancilla[0])
return apply
# -- hooks left to the reader ---------------------------------------
# The controlled variants (for QPE) follow the same pattern with the
# control on qubit 0 of a combined register; select_observable enables
# odd moments. Omitted here -- even moments below need none of them.
def controlled_apply_kernel(self):
raise NotImplementedError("exercise: controlled U_A for QPE")
def controlled_walk_step_kernel(self):
raise NotImplementedError("exercise: controlled W for QPE")
def controlled_adjoint_walk_step_kernel(self):
raise NotImplementedError("exercise: controlled W-dagger")
def select_observable(self):
raise NotImplementedError("exercise: enables odd Chebyshev moments")
def dense_moment(a, b, k, state):
"""Reference <T_k(H/alpha)> for H = a*X + b*Z."""
X = np.array([[0, 1], [1, 0]], dtype=complex)
Z = np.diag([1.0, -1.0]).astype(complex)
x = (a * X + b * Z) / (a + b)
previous, current = np.eye(2, dtype=complex), x
if k == 0:
current = np.eye(2, dtype=complex)
for _ in range(2, k + 1):
previous, current = current, 2 * x @ current - previous
return float(np.real(state.conj() @ current @ state))
def main() -> int:
cudaq.set_target(os.environ.get("CUDAQ_DEFAULT_SIMULATOR", "qpp-cpu"))
encoding = TwoTermLCU(a=0.6, b=0.9)
# It satisfies the protocol -- structurally, no base class involved.
print(f"isinstance(encoding, BlockEncoding) : "
f"{isinstance(encoding, BlockEncoding)}")
print(f"num_system / num_ancilla / alpha : "
f"{encoding.num_system} / {encoding.num_ancilla} / {encoding.alpha}")
# Walk consumes it with no changes. Measure even Chebyshev moments
# (even orders use the geometry-derived reflection observable, so they
# need no encoding-specific observable hook).
walk = Walk(encoding)
state = np.array([0.8, 0.6], dtype=np.complex128) # a chosen |psi>
print("\n k Walk.moment dense <T_k(H/alpha)>")
ok = True
for k in (0, 2, 4, 6):
measured = walk.moment(state, k)
reference = dense_moment(encoding.a, encoding.b, k, state)
flag = "ok" if abs(measured - reference) < 1e-10 else "MISMATCH"
ok = ok and flag == "ok"
print(f" {k} {measured:+.8f} {reference:+.8f} {flag}")
assert ok
print("\nOK — a from-scratch encoding drops into Walk unchanged. That is "
"the extension point: bring an encoding, inherit the primitives.")
return 0
if __name__ == "__main__":
raise SystemExit(main())
Solve a linear system with QSVT
The payoff of QSVT’s “pick a polynomial” design: example 2 approximated
exp(-ixt) for time evolution; swap the phase sequence for one
approximating 1/x and the same PauliLCU + QSVT machinery becomes a
quantum linear solver – the QSVT reading of HHL. A 5x5 symmetric
positive-definite system is padded to three qubits, block-encoded, and
inverted with a Childs-Kothari-Somma polynomial (degree set by the condition
number); the good-subspace state is proportional to A^-1 b, with the
scale recovered classically. Verified against numpy.linalg.solve.
Requires qsppack for the phase factors, like example 2.
"""Example 7 — Invert a 5x5 matrix with QSVT (the quantum linear solver).
Time evolution (example 2) applied a QSVT polynomial approximating
exp(-ixt). Swap the polynomial and the same machinery solves A x = b:
apply a polynomial approximating 1/x to the block-encoded matrix and the
good-subspace state is proportional to A^-1 b. This is the QSVT reading of
HHL, and it reuses PauliLCU + QSVT unchanged -- only the phase sequence is
new.
Three ideas carry the example:
* Embedding. A is 5x5; a register holds 2^k amplitudes. Pad A to the
next power of two (8x8, three qubits) with an identity block. The pad
is block-diagonal, so with b supported on the first five coordinates
the solution never leaks into the padding.
* The polynomial. 1/x is singular at 0, so it can only be approximated
away from zero, on the spectrum's domain [x_min, x_max]. We use the
Childs-Kothari-Somma closed form: (1 - (1-x^2)^b)/x is an odd
polynomial (degree 2b-1) with an exact Chebyshev expansion, and it
tracks 1/x wherever |x| is bounded away from 0. Degree grows like
kappa^2 log(1/eps), so a well-conditioned A keeps it small.
* Normalization. A QSVT polynomial must satisfy |p| <= 1 on [-1, 1],
but 1/x is unbounded. So we implement a *scaled* inverse p(x) ~ c/x;
the circuit returns c*alpha*A^-1 b, and the classical constant is
recovered from A and b afterward (standard HHL post-processing).
Verified against numpy.linalg.solve: direction fidelity ~1 and a small
residual ||A x_hat - b||. Run: python3 07_matrix_inversion_qsvt.py
"""
from __future__ import annotations
import contextlib
import io
import os
import cudaq
import numpy as np
from numpy.polynomial import chebyshev
from scipy.special import comb
from cudaq_algorithms import PauliLCU, PhaseSequence, QSVT
from cudaq_algorithms import sim_utils as sim
_PAULIS = {
"I": np.eye(2, dtype=complex),
"X": np.array([[0, 1], [1, 0]], dtype=complex),
"Y": np.array([[0, -1j], [1j, 0]]),
"Z": np.diag([1.0, -1.0]).astype(complex),
}
def pauli_terms(matrix: np.ndarray, num_qubits: int) -> dict:
"""Decompose a Hermitian matrix into a {pauli_word: coefficient} sum.
Little-endian, matching CUDA-Q and example 1: word[0] is the rightmost
tensor factor. Coefficients are c_P = tr(P^dagger M) / 2^n; for a
Hermitian M they come out real.
"""
terms: dict = {}
for index in range(4**num_qubits):
labels, digits = [], index
for _ in range(num_qubits):
labels.append("IXYZ"[digits % 4])
digits //= 4
word = "".join(labels)
factor = np.array([[1]], dtype=complex)
for label in word:
factor = np.kron(_PAULIS[label], factor)
coefficient = np.trace(factor.conj().T @ matrix) / (2**num_qubits)
if abs(coefficient) > 1e-10:
terms[word] = coefficient.real
return terms
def inverse_chebyshev_coeffs(b: int) -> np.ndarray:
"""Chebyshev coefficients of (1 - (1-x^2)^b)/x, the CKS 1/x proxy.
Returns the coefficients of T_1, T_3, ..., T_{2b-1} (odd parity) --
exactly the form qsppack.solve consumes. The identity is
(1-(1-x^2)^b)/x = 4 sum_j (-1)^j [2^-2b sum_{i=j+1}^b C(2b, b+i)] T_{2j+1}.
"""
coeffs = np.zeros(b)
for j in range(b):
tail = sum(comb(2 * b, b + i, exact=True) for i in range(j + 1, b + 1))
coeffs[j] = 4 * ((-1)**j) * tail / (2.0**(2 * b))
return coeffs
def evaluate_odd(coeffs: np.ndarray, x: np.ndarray) -> np.ndarray:
"""Evaluate sum_j coeffs[j] T_{2j+1}(x)."""
value = np.zeros_like(x, dtype=float)
for j, coefficient in enumerate(coeffs):
basis = np.zeros(2 * j + 2)
basis[2 * j + 1] = 1.0
value += coefficient * chebyshev.chebval(x, basis)
return value
def inverse_phase_sequence(x_min: float, epsilon: float):
"""Phase factors for a scaled 1/x on [x_min, 1], and the scale c.
Picks the CKS degree from the error bound |(1-x^2)^b| <= exp(-b x_min^2),
scales the polynomial to sit safely inside |p| <= 1 (so the QSP solve
stays well-conditioned), and returns qsp-convention phases.
"""
import qsppack
b = max(int(np.ceil(np.log(1 / epsilon) / x_min**2)), 3)
coeffs = inverse_chebyshev_coeffs(b)
# 1/x-proxy peaks at ~sqrt(b); rescale so |p| tops out near 0.9. Leaving
# headroom below 1 is what keeps the phase-factor solve well-conditioned.
grid = np.linspace(-1, 1, 4000)
scale = 0.9 / np.max(np.abs(evaluate_odd(coeffs, grid)))
coeffs = coeffs * scale
options = {
"criteria": 1e-8,
"method": "Newton",
"typePhi": "full",
"useReal": True,
"targetPre": True,
}
with contextlib.redirect_stdout(io.StringIO()): # hush the QSP solver
phases, _ = qsppack.solve(coeffs.copy(), 1, options)
return [float(p) for p in phases], scale, b
def main() -> int:
cudaq.set_target(os.environ.get("CUDAQ_DEFAULT_SIMULATOR", "qpp-cpu"))
# 1. A 5x5 symmetric positive-definite matrix with a controlled spectrum
# in [1, 2] (condition number 2). The smallest eigenvalue is 1, so
# padding to 8x8 with an identity block adds no new small eigenvalue.
rng = np.random.default_rng(3)
basis, _ = np.linalg.qr(rng.normal(size=(5, 5)))
spectrum = np.array([1.0, 1.25, 1.5, 1.75, 2.0])
a5 = (basis * spectrum) @ basis.T
a5 = 0.5 * (a5 + a5.T)
a8 = np.eye(8, dtype=complex)
a8[:5, :5] = a5
# 2. Block-encode A (Hermitian) as a Pauli LCU. QSVT applies polynomials
# to x = eigenvalue / alpha, so the working domain is [x_min, x_max].
terms = pauli_terms(a8, num_qubits=3)
encoding = PauliLCU(terms)
transformer = QSVT(encoding)
eigenvalues = np.linalg.eigvalsh(a8)
x = eigenvalues / encoding.alpha
x_min, x_max = float(x.min()), float(x.max())
print(f"matrix : 5x5 SPD, padded to 8x8 (3 system qubits)")
print(f"Pauli LCU terms : {len(terms)} "
f"(system {encoding.num_system} + ancilla {encoding.num_ancilla})")
print(f"alpha (1-norm) : {encoding.alpha:.6f}")
print(
f"condition number : {eigenvalues.max() / eigenvalues.min():.2f}")
print(f"spectral domain x : [{x_min:.4f}, {x_max:.4f}]")
# 3. Build phases for a scaled 1/x over the spectral domain, and check
# the polynomial tracks scale/x there before trusting the circuit.
phases, scale, degree_half = inverse_phase_sequence(x_min, epsilon=1e-2)
domain = np.linspace(x_min, x_max, 200)
polynomial = scale * evaluate_odd(inverse_chebyshev_coeffs(degree_half),
domain)
poly_error = float(np.max(np.abs(polynomial - scale / domain)))
print(f"CKS polynomial degree : {2 * degree_half - 1} "
f"(scaled by c = {scale:.5f})")
print(f"|p(x) - c/x| on domain: {poly_error:.2e}")
print(f"phase factors : {len(phases)}")
# 4. Solve A x = b. Prepare b on the system register, run the QSVT
# sequence, keep the all-zero-ancilla block. The qsp convention makes
# the good-subspace amplitude complex (p(x) + i * complementary(x));
# with a real b and real eigenvectors, the real part -- after removing
# the qsp global phase exp(i * sum phi) -- is p(A/alpha) b.
b = np.zeros(8, dtype=complex)
b[:5] = rng.normal(size=5)
b /= np.linalg.norm(b)
raw = sim.transform(transformer, b, PhaseSequence(phases,
convention="qsp"))
good = (raw * np.exp(-1j * np.sum(phases))).real
# 5. The state is proportional to A^-1 b. Recover the real constant from
# A and b (both known), then check the residual -- HHL post-processing.
# c = <b, A good> / <b, b> = scale * alpha up to the convention sign.
constant = float(np.real(np.vdot(b, a8 @ good)) / np.real(np.vdot(b, b)))
x_hat = good / constant
exact = np.linalg.solve(a8, b).real
fidelity = float(
abs(np.vdot(good / np.linalg.norm(good),
exact / np.linalg.norm(exact))))
residual = float(np.linalg.norm(a8 @ x_hat - b))
solution_error = float(
np.linalg.norm(x_hat - exact) / np.linalg.norm(exact))
print(f"\nrecovered constant c : {constant:+.5f} "
f"(|c| vs scale*alpha = {scale * encoding.alpha:.5f})")
print(f"direction fidelity : {fidelity:.6f}")
print(f"residual ||A x - b|| : {residual:.3e}")
print(f"solution rel. error : {solution_error:.3e}")
assert fidelity > 0.999
assert residual < 1e-2
print("\nOK — QSVT applied 1/x to the block-encoded matrix; the "
"good-subspace state solves A x = b.")
return 0
if __name__ == "__main__":
raise SystemExit(main())