THC qubitized phase estimation

Download Notebook - thc_qpe.ipynb

This example composes canonical quantum phase estimation (QPE) with the controlled qubitization walk for the BLISS-THC Hamiltonian in Caesura et al., arXiv:2501.06165. It uses 16 spatial orbitals and THC rank 16: a non-trivial instance containing 152 alias-sampled coefficients and 15 neighboring Givens rotations in each orbital transformation.

The implementation is adapted from the PREPARE and Select circuit construction in the original paper. Here, that construction is expressed using the repository’s reusable alias-sampling, SelectTHCCntrl, controlled-LCU, qubitization, and QPE components.

The complete circuit is too large for practical statevector simulation. Its constituent algorithms are tested numerically on smaller systems, while this notebook serves as an integration example: it assembles the full workflow at a realistic size and checks that the resulting Guppy program compiles.

The construction is layered rather than monolithic:

THC coefficients and rotations
        ↓
PREPARE + controlled SELECT + UNPREPARE
        ↓
LCUCntrl
        ↓
controlled qubitization = controlled LCU + controlled reflection
        ↓
QPE powered oracle

Each layer has a typed interface. The classical THC data can come from a chemistry workflow, and any signature-compatible PREPARE or SELECT can replace the generated oracle. All alias and Select registers remain coherent between walk steps. The reflection tests the full alias workspace. The SELECT record is cleared by its inverse QROM lookup before alias UNPREPARE; the phase-gradient resource stays outside the reflection.

from typing import no_type_check

from guppylang import guppy
from guppylang.std.builtins import array, nat, output
from guppylang.std.quantum import (
    collect_measurements,
    discard,
    discard_array,
    h,
    measure_array,
    qubit,
    x,
    cx,
)

from guppyalgos.primitives.gate_decompositions.cnx.cnx import cnx
from guppyalgos.algorithms.block_encoding.lcu import LCUCntrl
from guppyalgos.algorithms.phase_estimation import qpe
from guppyalgos.primitives.gate_decompositions.cnx.cnz import cnz
from guppyalgos.primitives.rotations import GivensCascadePhaseGradient
from guppyalgos.algorithms.state_preparation.alias_sampling import (
    AliasSamplingRegs,
    alias_samp_prep,
)
from guppyalgos.primitives.state_preparation.phase_gradient import Convention, phase_gradient
from guppyalgos.algorithms.block_encoding.thc import (
    SelectTHCCntrl,
    SelectTHCCntrlRegs,
    THCWalkTargetRegs,
    build_thc_lcu_data,
    generate_thc_parameters,
    load_select_registers,
)
from guppyalgos.utils import qarray, transversal

1. Supply classical THC data

THCParameters is the boundary between classical preprocessing and circuit construction. The coefficient arrays specify the one- and two-body weights; each row of the rotation arrays contains the neighboring-Givens angles for one orbital transformation. Angles are expressed as fractions of a full turn.

To place both kinds of coefficient in one alias-sampling table, picture the one-body terms as an extra column appended to the symmetric two-body matrix. The public term table represents entries in this column with nu=None; only the internal QROM encoding replaces None with the paper’s reserved integer nu=M. A one-body term therefore still has only the meaningful index mu.

For this standalone example, generate_thc_parameters supplies deterministic non-zero data. A chemistry application would replace this one call with coefficients and rotations obtained from its classical factorization workflow. With 16 orbitals and rank 16, the flattened LCU contains 136 two-body coefficients and 16 one-body coefficients.

n_orbitals = 16
thc_rank = 16
parameters = generate_thc_parameters(n_orbitals, thc_rank, seed=26)

thc_lcu_data = build_thc_lcu_data(
    parameters,
    rotation_precision_bits=8,
    alias_precision_bits=12,
)

{
    "orbitals": n_orbitals,
    "THC rank": thc_rank,
    "LCU coefficients": len(thc_lcu_data.alias_probabilities),
    "Givens rotations per transform": thc_lcu_data.n_givens,
    "alias index qubits": thc_lcu_data.n_alias_qubits,
}

{'orbitals': 16,
 'THC rank': 16,
 'LCU coefficients': 152,
 'Givens rotations per transform': 15,
 'alias index qubits': 8}

2. Inspect the composable LCU oracles

build_thc_lcu_data is limited to classical preprocessing:

  • It normalizes the THC coefficient magnitudes.

  • It encodes the Select records and Givens angles as QROM tables.

  • It returns the resulting Guppy QROMs and register dimensions.

  • It does not choose a rotation method or construct PREPARE, SELECT, or UNPREPARE.

This notebook makes the quantum choices explicitly: alias-sampling PREPARE and a phase-gradient Givens cascade. Another compatible Select or rotator could use the same preprocessed data. The main returned values are:

  • alias_probabilities contains the normalized coefficient magnitudes used to construct alias-sampling PREPARE below.

  • select_data_loader is the precomputed QROM for the sampled THC indices and flags.

  • qrom_1_and_2_body and qrom_2_body load the combined and two-body Givens angles.

  • The standard AliasSamplingRegs and SelectTHCCntrlRegs types compose into the PREPARE register below.

Below, alias_samp_prep constructs PREPARE from the normalized magnitudes. The phase-gradient preparation and cascade are also selected here, outside the data builder.

alias_prepare = alias_samp_prep(
    thc_lcu_data.alias_probabilities,
    thc_lcu_data.alias_precision,
)
select_data_qrom = thc_lcu_data.select_data_loader
qrom_1_and_2_body = thc_lcu_data.qrom_1_and_2_body
qrom_2_body = thc_lcu_data.qrom_2_body

n_alias_qubits = thc_lcu_data.n_alias_qubits
n_index_qubits = thc_lcu_data.n_index_qubits
n_keep_qubits = thc_lcu_data.n_keep_qubits
n_modes = thc_lcu_data.n_modes
n_givens = thc_lcu_data.n_givens
n_phase_qubits = thc_lcu_data.rotation_precision_bits

prepare_phase_gradient = phase_gradient(
    n_phase_qubits,
    convention=Convention.Standard,
)

@guppy.struct
class THCPrepareRegs:
    """Persistent preparation registers and phase-gradient resource."""

    alias_sampling: AliasSamplingRegs[n_alias_qubits, n_keep_qubits]
    select: SelectTHCCntrlRegs[n_index_qubits]
    phase_gradient: array[qubit, n_phase_qubits]

@guppy
@no_type_check
def prepare(regs: THCPrepareRegs) -> None:
    """Prepare the alias distribution and load its THC Select record."""
    alias_prepare(
        regs.alias_sampling.index,
        regs.alias_sampling.alternative,
        regs.alias_sampling.keep,
        regs.alias_sampling.comparison,
        regs.alias_sampling.comparison_result,
        False,
    )
    load_select_registers(select_data_qrom, regs.alias_sampling.index, regs.select)

@guppy
@no_type_check
def cntrl_select(
    control: qubit,
    prep_regs: THCPrepareRegs,
    target_regs: THCWalkTargetRegs[n_modes],
) -> None:
    """Apply the THC Select using the prepared phase-gradient state."""
    cascade = GivensCascadePhaseGradient(prep_regs.phase_gradient)
    select = SelectTHCCntrl(qrom_1_and_2_body, qrom_2_body, cascade, cnx)
    select.compose(control, prep_regs.select, target_regs)
    prep_regs.phase_gradient = select.cascade.phase_gradient

@guppy
@no_type_check
def unprepare(regs: THCPrepareRegs) -> None:
    """Reverse the Select-record lookup and alias preparation."""
    load_select_registers(select_data_qrom, regs.alias_sampling.index, regs.select)
    alias_prepare(
        regs.alias_sampling.index,
        regs.alias_sampling.alternative,
        regs.alias_sampling.keep,
        regs.alias_sampling.comparison,
        regs.alias_sampling.comparison_result,
        True,
    )

@guppy.struct
class THCQPERegs:
    """Persistent quantum registers used by THC QPE."""

    prep: THCPrepareRegs
    target: THCWalkTargetRegs[n_modes]

3. Build the controlled block encoding

THCQPERegs retains the complete preparation bundle, spin registers and phase-gradient state. LCUCntrl applies PREPARE, controlled SELECT and UNPREPARE to that bundle.

UNPREPARE reverses the preparation circuit, but SELECT can leave its workspace entangled with the system. It does not generally return the alias workspace to zero. We therefore allocate the registers once and retain them throughout QPE; nothing is discarded between walk steps.

@guppy
@no_type_check
def cntrl_block_encoding(control: qubit, regs: THCQPERegs) -> None:
    LCUCntrl(prepare, cntrl_select, unprepare).compose(
        control, regs.prep, regs.target
    )

4. Reflect the alias preparation state

  • The reflection flips the phase when the full alias workspace is zero and the external control is one:

\[ R=I-2\Pi,\qquad \Pi=|0\rangle\langle0|_{\mathrm{alias}}\otimes I_{\mathrm{SELECT,system,gradient}}. \]

Why SELECT registers are excluded

  • PREPARE loads the SELECT record from the alias index using load_select_registers.

  • SelectTHCCntrl preserves that record: it reverses its temporary index swaps and restores its control flags.

  • UNPREPARE repeats the QROM lookup with the unchanged alias index, clearing the SELECT record before reversing alias preparation:

\[ |j\rangle|d(j)\rangle|\psi_j\rangle \longmapsto |j\rangle|0\rangle|\psi_j\rangle. \]

  • The SELECT registers therefore return to zero at each reflection boundary and need no zero-test flags in this construction.

  • Alias PREPARE entangles the index with its workspace (the “junk” registers). SELECT acts differently on each selected term, so applying alias UNPREPARE afterward does not generally clear that workspace: it can remain entangled with the system.

  • The walk must therefore reflect about the joint all-zero alias state, including the index and every junk register. Reflecting only on the index would also flip components whose index is zero but whose workspace is nonzero, giving a different walk operator.

  • The five flags implement this joint test while keeping the registers separate. They are computed coherently and then uncomputed; measuring them would collapse the superposition. Flags are an implementation choice—the requirement is to test the full alias workspace.

  • The phase-gradient resource is excluded because it is prepared in a phase-gradient state, not the all-zero state.

Five temporary zero-test flags

Flag

Alias register tested

0

Index.

1

Alternative index.

2

Keep threshold.

3

Comparison register.

4

Comparison-result qubit.

  • Each flag starts at zero and becomes one exactly when its associated register is zero:

\[ |r_j\rangle|0\rangle_{f_j}\longmapsto |r_j\rangle|[r_j=0]\rangle_{f_j}. \]

  • zero_flag uses X gates to turn an all-zero condition into an all-one condition for cnx, then restores the input bits. A single qubit uses X on its flag followed by CX.

  • Compute the five flags without measuring the preparation registers.

  • Apply cnz(flags, control, cnx), contributing the phase \((-1)^{c f_0f_1f_2f_3f_4}\).

  • Uncompute by repeating the flag routine, restoring the flags to zero before discarding them. The phase change remains.

@guppy
def zero_flag_array[n: nat](qreg: array[qubit, n], flag: qubit) -> None:
    transversal(x, qreg)
    cnx(qreg, flag)
    transversal(x, qreg)


@guppy
def zero_flag_qubit(q: qubit, flag: qubit) -> None:
    x(flag)
    cx(q, flag)


@guppy.overload(zero_flag_array, zero_flag_qubit)
def zero_flag():
    """Toggle a flag when a qubit or register is zero."""


@guppy
def preparation_zero_flags(regs: THCPrepareRegs, flags: array[qubit, 5]) -> None:
    zero_flag(regs.alias_sampling.index, flags[0])
    zero_flag(regs.alias_sampling.alternative, flags[1])
    zero_flag(regs.alias_sampling.keep, flags[2])
    zero_flag(regs.alias_sampling.comparison, flags[3])
    zero_flag(regs.alias_sampling.comparison_result, flags[4])



@guppy
def reflect_preparation(control: qubit, regs: THCPrepareRegs) -> None:
    flags = qarray(5)
    preparation_zero_flags(regs, flags)
    cnz(flags, control, cnx)
    preparation_zero_flags(regs, flags)
    discard_array(flags)


@guppy
def cntrl_walk_power(control: qubit, unitary_regs: THCQPERegs, power: int) -> None:
    for _ in range(power):
        cntrl_block_encoding(control, unitary_regs)
        reflect_preparation(control, unitary_regs.prep)

5. Compose with canonical QPE

A three-qubit phase register asks the powered oracle for walk powers 1, 2, and 4 before applying the inverse QFT. A useful energy estimate requires preparing an eigenstate, or a state with significant overlap with an eigenstate, in the spin registers. This example leaves them in the computational zero state because its purpose is to demonstrate and type-check the complete THC-QPE composition.

The resulting circuit is intentionally too large for a useful notebook statevector test. check() validates types and qubit ownership. It does not prove that a register is zero or that the compiled circuit implements the intended Hamiltonian. Registers are discarded only after QPE and the phase measurement have finished.

n_qpe_qubits = 3


@guppy
@no_type_check
def main() -> None:
    phase_qreg = qarray(n_qpe_qubits)
    alias_regs = AliasSamplingRegs(
        qarray(n_alias_qubits), qarray(n_alias_qubits),
        qarray(n_keep_qubits), qarray(n_keep_qubits), qubit(),
    )
    select_regs = SelectTHCCntrlRegs(
        qubit(), qubit(), qarray(n_index_qubits), qarray(n_index_qubits),
    )
    gradient_qreg = qarray(n_phase_qubits)
    prepare_phase_gradient(gradient_qreg)
    unitary_regs = THCQPERegs(
        THCPrepareRegs(alias_regs, select_regs, gradient_qreg),
        THCWalkTargetRegs(qarray(n_modes), qarray(n_modes)),
    )
    transversal(h, phase_qreg)
    qpe(phase_qreg, unitary_regs, cntrl_walk_power)
    output("phase", collect_measurements(measure_array(phase_qreg)))

    # QPE is finished: these registers are no longer needed.
    discard_array(unitary_regs.prep.alias_sampling.index)
    discard_array(unitary_regs.prep.alias_sampling.alternative)
    discard_array(unitary_regs.prep.alias_sampling.keep)
    discard_array(unitary_regs.prep.alias_sampling.comparison)
    discard(unitary_regs.prep.alias_sampling.comparison_result)
    discard(unitary_regs.prep.select.one_body_flag)
    discard(unitary_regs.prep.select.coefficient_sign)
    discard_array(unitary_regs.prep.select.first_index_qreg)
    discard_array(unitary_regs.prep.select.second_index_qreg)
    discard_array(unitary_regs.prep.phase_gradient)
    discard_array(unitary_regs.target.spin_up)
    discard_array(unitary_regs.target.spin_down)


main.check()
"THC QPE with persistent preparation workspace type checked successfully."

'THC QPE with persistent preparation workspace type checked successfully.'