Block Encoding: From LCU to QSVT¶

Download Notebook - block_encoding_demo.ipynb

Encode a two-term Hamiltonian with LCU.
Reuse the block encoding for qubitization and a polynomial transformation with QSVT.
Extend the construction to eight terms and alias-sampling preparation with configurable QROM fanout.
The examples use guppy. Run from a source checkout with development dependencies installed; the checks use guppyalgos.tests.helpers.
The repository default is little endian.

Projecting out the encoded block¶

A block encoding stores a matrix inside a larger unitary. Initialize the preparation qubits at zero and project them back onto zero to extract the system block:

\[ M=(\langle0^a|\otimes I)\,U\,(|0^a\rangle\otimes I). \]

Each example passes an explicit projection dictionary to get_unitary_projected. Its entries follow the circuit’s preparation-register order; QSVT also projects its signal qubit onto zero.
Projection preserves the block’s scale. For normalized \(|\psi\rangle\), the post-selection probability is \(\|M|\psi\rangle\|^2\).
Compare the extracted matrix with the expected block, allowing an overall global phase.
Reference polynomials use Zixy Pauli sums: A * A is the operator product. Convert to a matrix only when comparing with the simulated block.

1. LCU: encode two terms¶

Choose \(H=0.6X+0.4Y\), with normalization \(\lambda=0.6+0.4=1\). One preparation qubit selects between X and Y:

\[ P|0\rangle=\sqrt{0.6}|0\rangle+\sqrt{0.4}|1\rangle, \qquad S=|0\rangle\langle0|\otimes X+|1\rangle\langle1|\otimes Y. \]

LCU applies PREPARE, SELECT, then UNPREPARE. The resulting unitary \(B\) contains the Hamiltonian in its zero-ancilla block:

\[\begin{split} \langle0|B|0\rangle=H/\lambda =\begin{pmatrix}0&0.6-0.4i\\0.6+0.4i&0\end{pmatrix}. \end{split}\]

hamiltonian = zqp.RealTermSum.from_str("(0.6, X0), (0.4, Y0)", 1)
data = LCUData.from_hamiltonian(hamiltonian)
prepare = multiplexor_prep(data.amplitudes)
select = build_single_cntrl_select(data)
A = zqp.ComplexTermSum.from_str(hamiltonian.to_str(), 1) / data.l1_norm


@guppy
def unprepare(prep_qreg: array[qubit, 1]) -> None:
    with dagger:
        prepare(prep_qreg)


@guppy
def encode(prep_qreg: array[qubit, 1], state_qreg: array[qubit, 1]) -> None:
    LCU(prepare, select, unprepare).compose(prep_qreg, state_qreg)


encoded_block = get_unitary_projected(
    encode, 1, {"prep": [False]},
    endianness=Endianness.LITTLE,
    n_extra_qubits=2,
)
assert_allclose_ignorephase(
    encoded_block, A.to_sparse_matrix().toarray(), threshold=1e-7,
)
print("LCU block matches H / normalization.")

LCU block matches H / normalization.

2. Qubitization: reuse the LCU¶

Add a reflection on the preparation qubit to make a quantum walk:

\[ W=RB,\qquad R=I-2|0\rangle\langle0|. \]

Qubitization takes the same LCU and the reflection. Two steps encode

\[ \langle0|W^2|0\rangle=T_2(A)=2A^2-I=0.04I, \qquad A=H/\lambda. \]

Use .compose(...) for one step or .power(..., d) for \(d\) steps. Keep the preparation register coherent between them.

@guppy
def walk_squared(prep_qreg: array[qubit, 1], state_qreg: array[qubit, 1]) -> None:
    Qubitization(
        LCU(prepare, select, unprepare), Reflection[1, 0](cnx)
    ).power(prep_qreg, state_qreg, 2)


walk_block = get_unitary_projected(
    walk_squared, 1, {"prep": [False]},
    endianness=Endianness.LITTLE,
    n_extra_qubits=2,
)

identity = zqp.ComplexTermSum.from_str("(1, I0)", 1)
walk_polynomial = 2 * (A * A) - identity
assert_allclose_ignorephase(
    walk_block, walk_polynomial.to_sparse_matrix().toarray(),
    threshold=1e-7,
)
print("Two walk steps encode 2A² - I.")

Two walk steps encode 2A² - I.

3. QSVT: choose a polynomial¶

QSVT combines the LCU, its adjoint and a phase sequence. Here \(B^\dagger=B\), so both arguments use the same LCU. The three half-turn phases [0.5, 0.5, 0.5] give the cubic

\[ p(x)=\frac{x-2x^3}{\sqrt{2}},\qquad \langle0|_s\langle0|_p V|0\rangle_s|0\rangle_p=p(A)=-\frac{0.04}{\sqrt{2}}A. \]

Project the signal and preparation qubits onto zero to extract the transformed block.
Since \(A^2=0.52I\), \(p(A)=-0.04A/\sqrt{2}\): the result is a small, sign-flipped copy of \(A\), not the identity.

phases = [0.5, 0.5, 0.5]


@guppy
def transform(
    prep_qreg: array[qubit, 1], signal_qreg: array[qubit, 1],
    state_qreg: array[qubit, 1],
) -> None:
    QSVT(
        LCU(prepare, select, unprepare),
        LCU(prepare, select, unprepare),
        comptime(phases),
    ).compose(signal_qreg[0], prep_qreg, state_qreg)


qsvt_block = get_unitary_projected(
    transform, 1, {"prep": [False], "signal": [False]},
    endianness=Endianness.LITTLE,
    n_extra_qubits=2,
)
qsvt_polynomial = (A - 2 * (A * A * A)) / sqrt(2)
assert_allclose_ignorephase(
    qsvt_block, qsvt_polynomial.to_sparse_matrix().toarray(),
    threshold=1e-7,
)
print("QSVT block matches (A - 2A³) / √2.")

QSVT block matches (A - 2A³) / √2.

4. Replace the Hamiltonian¶

For an arbitrary Hermitian Pauli sum \(H=\sum_j a_jP_j\), use LCUData to obtain \(\lambda=\sum_j|a_j|\) and the register sizes. build_unary_iteration_select handles the term selection. The QSVT composition stays the same.

For example, take

\[\begin{split} \begin{aligned} H={}&0.25Z_0-0.125X_1+0.125Y_0Y_1+0.125Z_0X_1\\ &+0.125X_0+0.125Z_1+0.0625X_0Z_1-0.0625Z_0Z_1. \end{aligned} \end{split}\]

Eight terms use three preparation qubits and two system qubits.
The LCU block is checked numerically against \(H/\lambda\), a nontrivial four-by-four matrix.
The same phase sequence defines \(p(H/\lambda)=[H/\lambda-2(H/\lambda)^3]/\sqrt{2}\).
The complete QSVT composition is then compiled.

larger_hamiltonian = zqp.RealTermSum.from_str(
    "(0.25, Z0), (-0.125, X1), (0.125, Y0 Y1), (0.125, Z0 X1), "
    "(0.125, X0), (0.125, Z1), (0.0625, X0 Z1), (-0.0625, Z0 Z1)", 2
)
larger_data = LCUData.from_hamiltonian(larger_hamiltonian)
n_prep = larger_data.n_prep_qubits
n_state = larger_data.n_state_qubits
larger_prepare = multiplexor_prep(larger_data.amplitudes)
larger_select = build_unary_iteration_select(larger_data)


@guppy
def larger_unprepare(prep_qreg: array[qubit, n_prep]) -> None:
    with dagger:
        larger_prepare(prep_qreg)


@guppy
def larger_encode(
    prep_qreg: array[qubit, n_prep], state_qreg: array[qubit, n_state],
) -> None:
    LCU(larger_prepare, larger_select, larger_unprepare).compose(
        prep_qreg, state_qreg,
    )


@guppy
def larger_transform(
    prep_qreg: array[qubit, n_prep], signal_qreg: array[qubit, 1],
    state_qreg: array[qubit, n_state],
) -> None:
    QSVT(
        LCU(larger_prepare, larger_select, larger_unprepare),
        LCU(larger_prepare, larger_select, larger_unprepare),
        comptime(phases),
    ).compose(signal_qreg[0], prep_qreg, state_qreg)


larger_A = zqp.ComplexTermSum.from_str(larger_hamiltonian.to_str(), n_state) / larger_data.l1_norm
larger_block = get_unitary_projected(
    larger_encode, n_state, {"prep": [False] * n_prep},
    endianness=Endianness.LITTLE,
    n_extra_qubits=n_prep - 1,
)
assert_allclose_ignorephase(
    larger_block, larger_A.to_sparse_matrix().toarray(),
    threshold=1e-7,
)
larger_transform.compile_function()
print("Eight-term LCU block verified and its QSVT composition compiled.")

Eight-term LCU block verified and its QSVT composition compiled.

5. LCU with alias-sampling PREPARE¶

Alias sampling prepares the weights of a larger table of terms. Unlike rotation-based preparation, it retains workspace \(|g_j\rangle\) alongside the index \(|j\rangle\):

\[ P|0\rangle=\sum_j\sqrt{\widetilde p_j}|j\rangle|g_j\rangle, \qquad p_j=|a_j|/\lambda. \]

SELECT acts on the index and system registers.
UNPREPARE reverses alias sampling on all preparation registers.
Together they give

\[ \langle0|P^\dagger SP|0\rangle =\sum_j\widetilde p_j\operatorname{sgn}(a_j)P_j. \]

This example reuses the eight-term Hamiltonian above.
Its weights are exactly representable with four probability bits, so \(\widetilde p_j=p_j\) and the block is \(H\) because \(\lambda=1\).
Other weights are rounded at the chosen precision.
See the alias-sampling notebook for the preparation routine in isolation.

from guppyalgos.algorithms.state_preparation.alias_sampling import (
    AliasSamplingRegs, alias_samp_prep,
)
from guppyalgos.primitives.subroutines.fanout import (
    fanout_basic, fanout_measurement_parity,
)

alias_hamiltonian = larger_hamiltonian
alias_data = LCUData.from_hamiltonian(alias_hamiltonian)
alias_probabilities = np.abs(alias_data.coeffs) / alias_data.l1_norm
unary_it_select = build_unary_iteration_select(alias_data)
n_precision_qubits = 4
alias_precision = 1 / 2**n_precision_qubits
n_alias_index = alias_data.n_prep_qubits
n_alias_state = alias_data.n_state_qubits


def build_alias_lcu(fanout_op=fanout_basic):
    alias_prepare = alias_samp_prep(
        alias_probabilities, precision=alias_precision, fanout_op=fanout_op,
    )

    @guppy
    def prepare_alias(
        prep_qregs: AliasSamplingRegs[n_alias_index, n_precision_qubits],
    ) -> None:
        alias_prepare(
            prep_qregs.index, prep_qregs.alternative, prep_qregs.keep,
            prep_qregs.comparison, prep_qregs.comparison_result, False,
        )

    @guppy
    def select_alias(
        prep_qregs: AliasSamplingRegs[n_alias_index, n_precision_qubits],
        state_qreg: array[qubit, n_alias_state],
    ) -> None:
        unary_it_select(prep_qregs.index, state_qreg)

    @guppy
    def unprepare_alias(
        prep_qregs: AliasSamplingRegs[n_alias_index, n_precision_qubits],
    ) -> None:
        alias_prepare(
            prep_qregs.index, prep_qregs.alternative, prep_qregs.keep,
            prep_qregs.comparison, prep_qregs.comparison_result, True,
        )

    @guppy
    def encode_alias(
        prep_qregs: AliasSamplingRegs[n_alias_index, n_precision_qubits],
        state_qreg: array[qubit, n_alias_state],
    ) -> None:
        LCU(prepare_alias, select_alias, unprepare_alias).compose(
            prep_qregs, state_qreg,
        )

    return encode_alias


alias_encode = build_alias_lcu()

alias_encode.compile_function()
print("Alias LCU with persistent preparation workspace compiled.")

Alias LCU with persistent preparation workspace compiled.

AliasSamplingRegs[3, 4] contains:
- a three-qubit index,
- a three-qubit alternative index,
- two four-qubit probability registers, and
- one comparison flag.
These make 15 preparation qubits in total.
Initialize all preparation qubits to zero and retain them through SELECT and UNPREPARE.
The encoded block projects all 15 preparation qubits onto zero.
UNPREPARE does not necessarily clear the workspace after SELECT. For this Hamiltonian acting on \(|00\rangle\), the comparison register remains nonzero with probability \(3/16\).
Keep the workspace through subsequent coherent operations. Discarding it would remove coherence.

6. Change the alias QROM fanout¶

Alias PREPARE uses QROM to load an alternative index and a keep threshold for each address. Its lookup XORs these words into two workspace registers:

\[ |j\rangle|u\rangle|v\rangle \longmapsto |j\rangle|u\oplus\mathrm{alt}_j\rangle|v\oplus\mathrm{keep}_j\rangle. \]

alias_samp_prep(..., fanout_op=...) chooses how the active address flag controls the stored one-bits. Its build_alias_fanout adapter gathers selected bits from both registers into one fanout call.

fanout_basic: sequential CNOTs.
fanout_measurement_parity: measurement-assisted fanout with feed-forward; four or more selected bits use extra ancillas. See the parity LAQCC notebook for the underlying construction.

Both implement the same lookup on arbitrary workspace states. This lets us change the QROM implementation without changing the probabilities, SELECT or LCU interface. UNPREPARE uses the alias routine’s explicit inverse mode, including the same QROM choice.

# Choose the fanout when constructing alias PREPARE.
measurement_prepare = alias_samp_prep(
    alias_probabilities,
    precision=alias_precision,
    fanout_op=fanout_measurement_parity,
)

# The LCU factory forwards the same option to PREPARE and UNPREPARE.
measurement_alias_encode = build_alias_lcu(
    fanout_op=fanout_measurement_parity,
)
measurement_alias_encode.compile_function()
print("Alias LCU with measurement-assisted QROM compiled.")

Alias LCU with measurement-assisted QROM compiled.

fanout_log is another option when every alias-table row selects at least one bit. The current implementation does not support empty fanouts, so it is not a drop-in choice for arbitrary alias tables. Use the basic or parity implementation when rows may be all zero.