Download Notebook - zixy_phase_estimation.ipynb
Trotterized phase estimation with Zixy¶
This notebook demonstrates trotterized QPE workflows using building blocks
from guppyalgos, with classical comparisons kept local for context and debugging.
Flow:
Canonical one-qubit \(0.5\,Z_0\) toy model.
Small commuting example.
Small noncommuting example.
H2 STO-3G Jordan-Wigner demo.
import numpy as np
import zixy.qubit.pauli as zqp
from scipy.linalg import expm
from typing import no_type_check
from guppylang import guppy
from guppylang.defs import GuppyFunctionDefinition
from guppylang.std.angles import angle
from guppylang.std.builtins import array, comptime, output, nat
from guppylang.std.quantum import (
discard_array,
measure_array,
qubit,
ry,
x,
collect_measurements
)
from guppyalgos.algorithms.phase_estimation import qpe
from guppyalgos.primitives.state_preparation import uniform_state
from guppyalgos.algorithms.time_evolution.trotter import cntrl_trotter_first_order
from guppyalgos.utils import (
fixed_point_to_float,
phase_distance_mod_2,
phase_to_energy_qpe,
qarray,
)
from zixy.fermion.mappings import JordanWignerMapper
from zixy.fermion.operator.general import String as FermionString
from zixy.qubit.pauli import RealTermSum as RealPauliOperator
Toy model¶
For \(H=\tfrac{1}{2}Z\) and \(t=-4\phi\),
\[\begin{split}
e^{-i \frac{\pi}{2} t H}
=
e^{i\pi\phi Z}
=
\begin{pmatrix}
e^{i\pi\phi} & 0 \\
0 & e^{-i\pi\phi}
\end{pmatrix}.
\end{split}\]
# Classical analysis
phi = 0.75 # phase to encode
rz_ham_op = RealPauliOperator.from_str("(0.5, Z0)")
n_rz_qubits = len(rz_ham_op.qubits)
rz_prepared_state = np.array([1.0, 0.0], dtype=complex)
### tools for trotter analysis (no error here)
# Under the current pauli_exp / trotter convention, choosing time_step = -4 * phi
# reproduces the same one-qubit unitary as CRz(-2 * pi * phi) in the canonical example.
rz_time_step = -4 * phi
rz_analysis = analyze_qubit_pauli_operator(
rz_ham_op, rz_time_step, rz_prepared_state, little_endian=True
)
rz_analysis_row = rz_analysis.iloc[0]
rz_target_phase = float(rz_analysis_row["target_phase"])
rz_ham_mat = rz_analysis_row["ham_mat"]
rz_unitary_from_ham = rz_analysis_row["exact_step"]
rz_unitary_direct = np.diag([
np.exp(1j * np.pi * phi),
np.exp(-1j * np.pi * phi),
])
assert np.allclose(rz_unitary_from_ham, rz_unitary_direct)
Phase Kickback (Single Eigenphase)¶
For an eigenstate \(|\psi\rangle\) with phase \(\phi\) (mod 2),
\[
U^{k}|\psi\rangle = e^{i\pi k\phi}|\psi\rangle
\quad\Rightarrow\quad
\mathrm{C}U^{k}\big(|+\rangle\otimes|\psi\rangle\big)
=
\frac{1}{\sqrt{2}}\Big(|0\rangle\otimes|\psi\rangle + e^{i\pi k\phi}|1\rangle\otimes|\psi\rangle\Big).
\]
QPE Wrapper¶
Make an abstract QPE program template once, then reuse it for each operator example below.
from guppylang.std.qsystem.helios import collect_measurements
n_ancilla = 4
def make_trotter_qpe_program[n_state_q: nat](
n_ancilla: int,
n_state_qubits: int,
state_preparation: GuppyFunctionDefinition[[array[qubit, n_state_q]], None],
power_oracle: GuppyFunctionDefinition[[qubit, array[qubit, n_state_q], int], None],
) -> GuppyFunctionDefinition[[], None]:
"""Construct a canonical QPE entrypoint for a single system register."""
ancilla_prep_function = uniform_state(2**n_ancilla)
@guppy
@no_type_check
def canonical_trotter_qpe_program() -> None:
state_reg = qarray(n_state_qubits)
phase_reg = qarray(n_ancilla)
state_preparation(state_reg)
ancilla_prep_function(phase_reg)
qpe(phase_reg, state_reg, power_oracle)
output("qpe_bitstring", collect_measurements(measure_array(phase_reg)))
discard_array(state_reg)
return canonical_trotter_qpe_program
# Guppy methods for simple wrapped Rz ham
controlled_rz_trotter_step = cntrl_trotter_first_order(rz_ham_op, n_rz_qubits)
# Powers of U are implemented as repeated applications of the controlled Trotter steps.
@guppy
def power_oracle(
control: qubit,
state_reg: array[qubit, n_rz_qubits],
power: int,
) -> None:
'QPE-compatible power oracle built from repeated controlled Trotter steps.'
for _ in range(power):
controlled_rz_trotter_step(control, state_reg, rz_time_step)
@guppy
def state_preparation(
system_register: array[qubit, n_rz_qubits],
) -> None:
"""Prepare the |0> eigenstate of the toy Rz Hamiltonian."""
ry(system_register[0], angle(0)) #identity!
# Assemble the simple QPE program
rz_qpe_program = make_trotter_qpe_program(
n_ancilla=n_ancilla,
n_state_qubits=n_rz_qubits,
state_preparation=state_preparation,
power_oracle=power_oracle,
)
# Run and inspect the simple qpe program
n_shots = 500
sim_result = (
rz_qpe_program.emulator(n_qubits=n_ancilla + n_rz_qubits)
.with_seed(5)
.with_shots(n_shots)
.run()
)
rz_counts = sim_result.register_counts()["qpe_bitstring"]
dominant_key, dominant_counts, dominant_phase = dominant_measured_phase(rz_counts)
rz_dominant_trotter_phase, rz_dominant_trotter_weight = dominant_trotter_phase(
rz_analysis
)
rz_analysis_summary = rz_analysis[
[
"target_phase",
"dominant_trotter_phase",
"dominant_trotter_overlap",
"trotter_step_error",
]
]
rz_counts_df = measurement_dataframe(rz_counts)
print(rz_analysis_summary.to_string(index=False, float_format=lambda value: f"{value:.4f}"))
print(rz_counts_df.to_string(index=False, float_format=lambda value: f"{value:.4f}"))
target_phase dominant_trotter_phase dominant_trotter_overlap trotter_step_error
0.75 0.75 1.0 0.0
bitstring counts phase empirical_probability
0110 500 0.75 1.0
First-Order Trotterization¶
We approximate time evolution with a first-order product formula:
\[
U(t) = e^{-i \frac{\pi}{2} t \sum_k c_k P_k}
\;\approx\;
\prod_k e^{-i \frac{\pi}{2} t\, c_k P_k},
\]
where each \(P_k\) is a Pauli string and \(c_k\in\mathbb{R}\).
See the Trotterized Hamiltonian simulation notebook for more details.
In this case we do controlled-U using the conjugation pattern.
Three Commuting Terms¶
We can also use a small Hamiltonian with three commuting Pauli terms,
Z0, Z1, and Z0 Z1. In this case the first-order Trotter step is exact,
so there is no Trotter error to confuse the phase-estimation example.
Commuting Case¶
If all \(P_k\) commute, the product formula is exact:
\[
\prod_k e^{-i \frac{\pi}{2} t\, c_k P_k}
=
e^{-i \frac{\pi}{2} t \sum_k c_k P_k}.
\]
# Classical analysis for the commuting ham
commuting_ham_op = RealPauliOperator.from_str("(0.5, Z0), (0.25, Z1), (0.25, Z0 Z1)")
n_commuting_qubits = len(commuting_ham_op.qubits)
commuting_prepared_state = np.array([0.0, 0.0, 0.0, 1.0], dtype=complex)
commuting_time_step = 1.0
commuting_analysis = analyze_qubit_pauli_operator(
commuting_ham_op,
commuting_time_step,
commuting_prepared_state,
little_endian=True,
)
commuting_analysis_row = commuting_analysis.iloc[0]
commuting_exact_step = commuting_analysis_row["exact_step"]
commuting_trotter_step = commuting_analysis_row["trotter_step"]
commuting_target_phase = float(commuting_analysis_row["target_phase"])
commuting_energies = commuting_analysis_row["energies"]
print(commuting_analysis_row["trotter_step_error"])
1.9238135806952475e-16
# Guppy methods for simple commuting ham
controlled_commuting_trotter_step = cntrl_trotter_first_order(
commuting_ham_op, n_commuting_qubits
)
# Powers of U are implemented as repeated applications of the controlled Trotter steps.
@guppy
def power_oracle_commuting(
control: qubit,
state_reg: array[qubit, n_commuting_qubits],
power: int,
) -> None:
for _ in range(power):
controlled_commuting_trotter_step(
control, state_reg, commuting_time_step
)
@guppy
def commuting_state_preparation(
system_register: array[qubit, n_commuting_qubits],
) -> None:
'Prepare the |11> eigenstate of the commuting Hamiltonian.'
x(system_register[0])
x(system_register[1])
# Assemble the commuting hamiltonian QPE program
commuting_qpe_program = make_trotter_qpe_program(
n_ancilla=n_ancilla,
n_state_qubits=n_commuting_qubits,
state_preparation=commuting_state_preparation,
power_oracle=power_oracle_commuting,
)
# Run and inspect commuting hamiltonian qpe program results
sim_result = (
commuting_qpe_program.emulator(n_qubits=n_ancilla + n_commuting_qubits)
.with_seed(5)
.with_shots(n_shots)
.run()
)
commuting_counts = sim_result.register_counts()["qpe_bitstring"]
dominant_key, dominant_counts, dominant_phase = dominant_measured_phase(
commuting_counts
)
commuting_dominant_trotter_phase, commuting_dominant_trotter_weight = dominant_trotter_phase(
commuting_analysis
)
commuting_summary = commuting_analysis[
[
"target_phase",
"dominant_trotter_phase",
"dominant_trotter_overlap",
"trotter_step_error",
]
]
commuting_counts_df = measurement_dataframe(commuting_counts)
print(commuting_summary.to_string(index=False, float_format=lambda value: f"{value:.4f}"))
print(commuting_counts_df.to_string(index=False, float_format=lambda value: f"{value:.4f}"))
target_phase dominant_trotter_phase dominant_trotter_overlap trotter_step_error
0.25 0.25 1.0 1.923814e-16
bitstring counts phase empirical_probability
0100 500 0.25 1.0
Noncommuting Terms and Trotter Error¶
Now switch to a genuinely noncommuting Hamiltonian. The exact evolution and the first-order Trotter step no longer match, so the trotterized QPE oracle estimates the eigenphases of the approximate step, not the exact Hamiltonian phases.
# Classical analysis of non-commuting ham
noncommuting_ham_op = RealPauliOperator.from_str(
"(0.6708203932499369, X0), (-1.3416407864998738, Z0)"
)
n_noncommuting_qubits = len(noncommuting_ham_op.qubits)
# This keeps the same X:Z ratio as before, but sets the
# trotterized eigenphase near 0.5 (with 4 ancilla bits).
noncommuting_time_step = 2.236
noncommuting_ham_mat = noncommuting_ham_op.to_sparse_matrix(True).toarray()
noncommuting_energies, noncommuting_eigenvectors = np.linalg.eigh(noncommuting_ham_mat)
noncommuting_prepared_state = canonical_real_single_qubit_state(
noncommuting_eigenvectors[:, 0]
)
noncommuting_state_prep_angle = ry_angle_for_real_single_qubit_state(
noncommuting_prepared_state
)
noncommuting_prepared_state_from_ry = np.array(
[
np.cos(noncommuting_state_prep_angle / 2),
np.sin(noncommuting_state_prep_angle / 2),
]
)
assert np.allclose(noncommuting_prepared_state_from_ry, noncommuting_prepared_state)
noncommuting_analysis = analyze_qubit_pauli_operator(
noncommuting_ham_op,
noncommuting_time_step,
noncommuting_prepared_state,
little_endian=True,
)
noncommuting_analysis_row = noncommuting_analysis.iloc[0]
noncommuting_target_phase = float(noncommuting_analysis_row["target_phase"])
noncommuting_exact_step = noncommuting_analysis_row["exact_step"]
noncommuting_trotter_step = noncommuting_analysis_row["trotter_step"]
noncommuting_trotter_phases = noncommuting_analysis_row["trotter_phases"]
noncommuting_trotter_overlaps = noncommuting_analysis_row["trotter_overlaps"]
assert not np.allclose(noncommuting_trotter_step, noncommuting_exact_step)
print(noncommuting_analysis_row["trotter_step_error"])
2.479548026725464
# Guppy methods for the non-commuting example
controlled_noncommuting_trotter_step = cntrl_trotter_first_order(
noncommuting_ham_op, n_noncommuting_qubits
)
@guppy
def noncommuting_state_preparation(
system_register: array[qubit, n_noncommuting_qubits],
) -> None:
'Prepare the exact ground state of the scaled noncommuting Hamiltonian.'
ry(system_register[0], angle(noncommuting_state_prep_angle))
@guppy
def power_oracle_noncommuting(
control: qubit,
state_reg: array[qubit, n_noncommuting_qubits],
power: int,
) -> None:
for _ in range(power):
controlled_noncommuting_trotter_step(
control, state_reg, noncommuting_time_step
)
# Assemble the non-commuting hamiltonian QPE program
noncommuting_qpe_program = make_trotter_qpe_program(
n_ancilla=n_ancilla,
n_state_qubits=n_noncommuting_qubits,
state_preparation=noncommuting_state_preparation,
power_oracle=power_oracle_noncommuting,
)
# Run and inspect non-commuting hamiltonian qpe program results
sim_result = (
noncommuting_qpe_program.emulator(
n_qubits=n_ancilla + n_noncommuting_qubits
)
.with_seed(5)
.with_shots(n_shots)
.run()
)
noncommuting_counts = sim_result.register_counts()["qpe_bitstring"]
dominant_key, dominant_counts, dominant_phase = dominant_measured_phase(
noncommuting_counts
)
noncommuting_dominant_trotter_phase, noncommuting_dominant_trotter_weight = dominant_trotter_phase(
noncommuting_analysis
)
noncommuting_summary = noncommuting_analysis[
[
"target_phase",
"trotter_phases",
"trotter_overlaps",
"dominant_trotter_phase",
"dominant_trotter_overlap",
"trotter_step_error",
]
]
noncommuting_counts_df = measurement_dataframe(noncommuting_counts)
print(noncommuting_summary.to_string(index=False, float_format=lambda value: f"{value:.4f}"))
print(noncommuting_counts_df.to_string(index=False, float_format=lambda value: f"{value:.4f}"))
target_phase trotter_phases trotter_overlaps dominant_trotter_phase dominant_trotter_overlap trotter_step_error
1.677 [0.49996775775044966, 1.5000322422495505] [0.816182459549431, 0.18381754045056925] 0.499968 0.816182 2.479548
bitstring counts phase empirical_probability
0010 257 0.5 0.514
0011 243 1.5 0.486
H2 STO-3G Trotterized QPE¶
Finally, apply the same trotterized QPE pattern to the full H2 STO-3G Jordan–Wigner Hamiltonian defined at the start of the notebook. We use the Hartree–Fock basis state as the input guess and compare it against the exact diagonalization before running QPE.
Phase to Energy¶
With \(U(t)=e^{-i \frac{\pi}{2} t H}\), a measured phase \(\phi\) corresponds to
\[
E = -\frac{2\phi}{t}\quad (\text{mod } 2/t).
\]
# prepare the trotterized hamiltonian
h2_ham_qubits=4
jw = JordanWignerMapper(h2_ham_qubits, mode_ordering=None)
H2_STO3G_HF_ENERGY = -1.1175058842043306
H2_STO3G_FCI_ENERGY = -1.136846575472054
H2_STO3G_FERMION_OPERATOR_ZIXY = [
(0.7430177069924181, ()),
(-1.2702927243904383, ((0, 1), (0, 0))), # F0^ F0
(-1.2702927243904380, ((1, 1), (1, 0))),
(-0.45680735030941033, ((2, 1), (2, 0))),
(-0.45680735030941020, ((3, 1), (3, 0))),
(0.6800618575841273, ((0, 1), (0, 0), (1, 1), (1, 0))), # F0^ F0 F1^ F1 #create/annihilate in mode 0, then create/annihilate in mode 1
(0.48890859745047305, ((0, 1), (0, 0), (2, 1), (2, 0))),
(0.6685772770134887, ((0, 1), (0, 0), (3, 1), (3, 0))),
(0.6685772770134887, ((1, 1), (1, 0), (2, 1), (2, 0))),
(0.48890859745047305, ((1, 1), (1, 0), (3, 1), (3, 0))),
(0.7028135332762804, ((2, 1), (2, 0), (3, 1), (3, 0))),
(-0.35933735912603115, ((0, 1), (1, 1), (2, 0), (3, 0))), # F0^ F1^ F2 F3 #create in modes 0,1, then annihilate in modes 2,3
(-0.35933735912603115, ((0, 1), (1, 0), (3, 1), (2, 0))), # F0^ F1 F3^ F2
]
h2_ham_op = RealPauliOperator.from_str("(0.0, I0)", h2_ham_qubits)
for coeff, ops in H2_STO3G_FERMION_OPERATOR_ZIXY:
if len(ops) == 0:
h2_ham_op += RealPauliOperator.from_str(f"({coeff}, I0)", h2_ham_qubits)
else:
fermion_string = FermionString(
h2_ham_qubits, [(mode, bool(adj)) for mode, adj in ops]
)
h2_ham_op += jw.encode_real(fermion_string, coeff)
assert np.isclose(
np.linalg.eigh(h2_ham_op.to_sparse_matrix(False).toarray())[0][0],
H2_STO3G_FCI_ENERGY,
)
n_h2_qubits = len(h2_ham_op.qubits)
n_h2_ancilla = 6
h2_time_step = 0.27508092748651924 # 0.25
# The Hartree-Fock computational basis state is |0011> in the big-endian matrix basis.
# In Guppy's little-endian qarray convention, that means flipping q[0] and q[1].
h2_hf_state = np.zeros(2**n_h2_qubits, dtype=complex)
h2_hf_state[3] = 1.0
h2_analysis = analyze_qubit_pauli_operator(
h2_ham_op,
h2_time_step,
h2_hf_state,
little_endian=False,
)
h2_analysis_row = h2_analysis.iloc[0]
h2_ground_energy = float(h2_analysis_row["energies"][0])
h2_ground_phase = float(h2_analysis_row["target_phase"])
h2_hf_ground_overlap = float(h2_analysis_row["exact_overlaps"][0])
h2_trotter_step_error = float(h2_analysis_row["trotter_step_error"])
# Guppy methods
controlled_h2_trotter_step = cntrl_trotter_first_order(h2_ham_op, n_h2_qubits)
@guppy
def power_oracle_h2(
control: qubit,
state_reg: array[qubit, n_h2_qubits],
power: int,
) -> None:
for _ in range(power):
controlled_h2_trotter_step(control, state_reg, h2_time_step)
@guppy
def h2_hf_state_preparation(
system_register: array[qubit, n_h2_qubits],
) -> None:
'Prepare the JW Hartree-Fock basis state |0011>.'
x(system_register[0])
x(system_register[1])
# Assemble the H2 trotterized QPE program
h2_qpe_program = make_trotter_qpe_program(
n_ancilla=n_h2_ancilla,
n_state_qubits=n_h2_qubits,
state_preparation=h2_hf_state_preparation,
power_oracle=power_oracle_h2,
)
# Run and inspect
sim_result = (
h2_qpe_program.emulator(n_qubits=n_h2_ancilla + n_h2_qubits)
.with_seed(5)
.with_shots(n_shots)
.run()
)
h2_counts = sim_result.register_counts()["qpe_bitstring"]
dominant_key, dominant_counts, dominant_phase = dominant_measured_phase(h2_counts)
h2_dominant_trotter_phase, h2_dominant_trotter_weight = dominant_trotter_phase(
h2_analysis
)
h2_summary = h2_analysis[
[
"target_phase",
"dominant_trotter_phase",
"dominant_trotter_overlap",
"trotter_step_error",
]
]
h2_counts_df = measurement_dataframe(h2_counts)
h2_counts_df["approximate energy (negative branch)"] = h2_counts_df["phase"].map(
lambda phase: phase_to_energy_qpe(phase, h2_time_step)
)
print(h2_summary.to_string(index=False, float_format=lambda value: f"{value:.4f}"))
print(h2_counts_df.to_string(index=False, float_format=lambda value: f"{value:.4f}"))
target_phase dominant_trotter_phase dominant_trotter_overlap trotter_step_error
0.156362 0.15625 0.988019 0.038557
bitstring counts phase empirical_probability approximate energy (negative branch)
101000 497 0.15625 0.994 -1.136029
011111 2 1.93750 0.004 -14.086764
101111 1 1.90625 0.002 -13.859558