Skip to main content

View on GitHub

Open this notebook in GitHub to run it yourself
This notebook is based on Ref. [1]. Given an efficient block-encoding for a Hamiltonian, this algorithm preforms an efficient Quantum Phase Estimation. The core quantum function of this model uses almost all Qmod built-in operations: control, power, within_apply, and invert. The algorithm assumes we have the block-encoding of a matrix HH U(s,m)−H=(H/s∗∗∗),U_{(s,m)-H} =\begin{pmatrix} H/s & * \\ * & * \end{pmatrix}, with mm being the size of the block variable and ss some scaling factor. Given this quantum function, we can define the following unitary (usually called the Szegedy quantum walk operator [2]): W≡Π∣0⟩mU(s,m)−H,W\equiv \Pi_{|0\rangle_m} U_{(s,m)-H}, where Π∣0⟩m\Pi_{|0\rangle_m} is a reflection operator about the block state
  1. The spectrum of the walk operator has a nice relation to the spectrum of the block-encoded Hamiltonian [3]:
{eigenvalues:}e{±iarccos⁡(λ/s)},{witheigenvectors:}∣φ{±}{λ}⟩≡{1}{{2}}(∣v{λ}⟩∣0⟩m±i∣⊥{λ}⟩),(1) \text\{eigenvalues: \} e^\{\pm i \arccos(\lambda/s)\}, \text\{ with eigenvectors: \} |\varphi^\{\pm\}_\{\lambda\}\rangle \equiv \frac\{1\}\{\sqrt\{2\}\}\left(|v_\{\lambda\}\rangle |0\rangle_m \pm i|\perp_\{\lambda\}\rangle\right), \quad (1) where ∣vλ⟩|v_\lambda\rangle is an eigenstate of the Hamiltonian HH with an eigenvalue λ\lambda. Namely, the eigenphases of HH are related by some nonlinear function (arccos⁡\arccos) to the eigenvalues of HH. The algorithm works under the assumption that the block-encoding unitary itself is also Hermitian, that is, U(s,m)−HU_{(s,m)-H} is Unitary and Hermitian.

​
Preliminaries

We start with defining some utility functions that are not implemented as part of Classiq, and might be included in the future. These functions are used for the specific block-encoding used in this notebook. 
!pip install -qq "classiq[chemistry]"


import matplotlib.pyplot as plt
import numpy as np

from classiq import *
from classiq.applications.chemistry.op_utils import qubit_op_to_qmod


from sympy import fwht

from classiq.open_library.functions.state_preparation import apply_phase_table


def get_graycode(size, i) -> int:
    if i == 2**size:
        return get_graycode(size, 0)
    return i ^ (i >> 1)


def get_graycode_angles_wh(size, angles):
    transformed_angles = fwht(np.array(angles) / 2**size)
    return [transformed_angles[get_graycode(size, j)] for j in range(2**size)]


def get_graycode_ctrls(size):
    return [
        (get_graycode(size, i) ^ get_graycode(size, i + 1)).bit_length() - 1
        for i in range(2**size)
    ]


@qfunc
def multiplex_ra(a_y: float, a_z: float, angles: list[float], qba: QArray, ind: QBit):
    assert a_y**2 + a_z**2 == 1
    # TODO support general (0,a_y,a_z) rotation
    assert (
        a_z == 1.0 or a_y == 1.0
    ), "currently only strict y or z rotations are supported"
    size = max(1, (len(angles) - 1).bit_length())
    extended_angles = angles + [0] * (2**size - len(angles))
    transformed_angles = get_graycode_angles_wh(size, extended_angles)
    controllers = get_graycode_ctrls(size)

    for k in range(2**size):
        if a_z == 0.0:
            RY(transformed_angles[k], ind)
        else:
            RZ(transformed_angles[k], ind)

        skip_control(lambda: CX(qba[controllers[k]], ind))


@qfunc
def lcu_paulis_graycode(terms: list[SparsePauliTerm], data: QArray, block: QArray):
    n_qubits = data.len
    n_terms = len(terms)
    table_z = np.zeros([n_qubits, n_terms])
    table_y = np.zeros([n_qubits, n_terms])
    probs = [abs(term.coefficient) for term in terms] + [0.0] * (2**block.len - n_terms)
    hamiltonian_coeffs = np.angle([term.coefficient for term in terms]).tolist() + [
        0.0
    ] * (2**block.len - n_terms)
    accumulated_phase = np.zeros(2**block.len).tolist()

    for k in range(n_terms):
        for pauli in terms[k].paulis:
            if pauli.pauli == Pauli.Z:
                table_z[pauli.index, k] = -np.pi
                accumulated_phase[k] += np.pi / 2
            elif pauli.pauli == Pauli.Y:
                table_y[pauli.index, k] = -np.pi
                accumulated_phase[k] += np.pi / 2
            elif pauli.pauli == Pauli.X:
                table_z[pauli.index, k] = -np.pi
                table_y[pauli.index, k] = np.pi
                accumulated_phase[k] += np.pi / 2

    def select_graycode(block: QArray, data: QArray):
        for i in range(n_qubits):
            multiplex_ra(0, 1, table_z[i, :], block, data[i])
            multiplex_ra(1, 0, table_y[i, :], block, data[i])
        apply_phase_table(
            [p1 - p2 for p1, p2 in zip(hamiltonian_coeffs, accumulated_phase)], block
        )

    within_apply(
        lambda: inplace_prepare_state(probs, 0.0, block),
        lambda: select_graycode(block, data),
    )

​
Defining a specific usecase: a Molecule and its block-encoding Hamiltonian function

molecule_H2_geometry = [("H", (0.0, 0.0, 0)), ("H", (0.0, 0.0, 0.735))]


from openfermion.chem import MolecularData
from openfermionpyscf import run_pyscf

basis = "sto-3g"  # Basis set
multiplicity = 1  # Singlet state S=0
charge = 0  # Neutral molecule
molecule = MolecularData(molecule_H2_geometry, basis, multiplicity, charge)

molecule = run_pyscf(
    molecule,
    run_fci=True,  # relevant for small, classically solvable problems
)


from classiq.applications.chemistry.mapping import FermionToQubitMapper
from classiq.applications.chemistry.problems import FermionHamiltonianProblem

# Define a  Hamiltonian in an active space
problem = FermionHamiltonianProblem.from_molecule(molecule=molecule)
mapper = FermionToQubitMapper()

qubit_hamiltonian = mapper.map(problem.fermion_hamiltonian)
print("Your Hamiltonian is", qubit_hamiltonian, sep="\n")
Output:

Your Hamiltonian is
  (-0.09057898608834769+0j) [] +
  (0.04523279994605784+0j) [X0 X1 X2 X3] +
  (0.04523279994605784+0j) [X0 X1 Y2 Y3] +
  (0.04523279994605784+0j) [Y0 Y1 X2 X3] +
  (0.04523279994605784+0j) [Y0 Y1 Y2 Y3] +
  (0.17218393261915538+0j) [Z0] +
  (0.12091263261776627+0j) [Z0 Z1] +
  (0.16892753870087907+0j) [Z0 Z2] +
  (0.1661454325638241+0j) [Z0 Z3] +
  (-0.2257534922240238+0j) [Z1] +
  (0.1661454325638241+0j) [Z1 Z2] +
  (0.17464343068300453+0j) [Z1 Z3] +
  (0.1721839326191554+0j) [Z2] +
  (0.12091263261776627+0j) [Z2 Z3] +
  (-0.22575349222402386+0j) [Z3]
Finally, we calculate the ground state energy as a reference solution to the quantum solver 
classical_sol = molecule.fci_energy
print(f"Expected energy: {classical_sol} Ha")
Output:

Expected energy: -1.1373060357533995 Ha
  


mol_hamiltonian = qubit_op_to_qmod(qubit_hamiltonian)
num_qubits = mol_hamiltonian.num_qubits
be_scaling = sum(np.abs(term.coefficient) for term in mol_hamiltonian.terms)
normalized_mol_hamiltonian = mol_hamiltonian * (1 / be_scaling)


data_size = normalized_mol_hamiltonian.num_qubits
num_terms = len(normalized_mol_hamiltonian.terms)
block_size = (num_terms - 1).bit_length() if num_terms != 1 else 1


print(f"The block size is {block_size}, and the scaling factor s is : {be_scaling}")
Output:

The block size is 4, and the scaling factor s is : 1.9850721353060015
  


class BlockEncodedState(QStruct):
    data: QNum[data_size]
    block: QNum[block_size]


@qfunc
def be_hamiltonian(state: BlockEncodedState):
    lcu_paulis_graycode(normalized_mol_hamiltonian.terms, state.data, state.block)

​
Defining a Walk Operator

We use the reflect_around_zero function from Classiq’s open library to define a walk operator function with the declaration below. This function implements I−2∣0⟩⟨0∣I-2|0\rangle\langle 0|, so we must insert a minus phase (This can be done by adding a minus sign using the phase function). 
from classiq.qmod.symbolic import pi


@qfunc
def my_reflect_about_zero(qba: QNum):
    control(qba == 0, lambda: phase(pi))
    phase(pi)


@qfunc
def walk_operator(
    be_qfunc: QCallable[BlockEncodedState], state: BlockEncodedState
) -> None:
    be_qfunc(state)
    my_reflect_about_zero(state.block)
We define a classical function that takes the eigenphases of the Walk operator, and returns the (scaled) eigenvalues of HH, according to equation (1) above. 
def post_process_walk_phases(w_eigphase, be_scaling):
    return np.cos(2 * np.pi * w_eigphase) * be_scaling
We also define a utility function for ploting the results: 
def get_qpe_walk_result(df, be_scaling, to_plot=True):
    filtered_block_res = df[df["block"] == 0].copy()
    block_prob = df.loc[df["block"] == 0, "probability"].sum()
    print(f"probability to measure the block variable at state zero: {block_prob}")
    filtered_block_res["post_processed_phases"] = post_process_walk_phases(
        filtered_block_res["phase_var"], be_scaling
    )

    max_prob_energy = filtered_block_res.loc[
        filtered_block_res["probability"].idxmax(), "post_processed_phases"
    ]
    print(f"\nEnergy with maximal probability: {max_prob_energy} Ha")

    if to_plot:
        plt.plot(
            filtered_block_res["post_processed_phases"],
            filtered_block_res["probability"],
            "o",
        )
        plt.xlabel("Energy (Ha)", fontsize=16)
        plt.ylabel("P(Energy)", fontsize=16)
        plt.tick_params(axis="both", labelsize=16)
        plt.title("Energy Histogram from QPE")

    return filtered_block_res

​
Setting initial state and QPE size

from classiq.applications.chemistry.hartree_fock import get_hf_state

# We take the Hartree Fock state
hf_state = get_hf_state(problem, mapper)

QPE_SIZE = 5

​
A Naive QPE

Before going to the optimized implementation, designing a specific QPE that operates on a walk operator, we start with a naive QPE implementation. 
@qfunc
def main(
    block: Output[QNum[block_size]],
    phase_var: Output[QNum[QPE_SIZE, SIGNED, QPE_SIZE]],
) -> None:

    data = QNum(size=data_size)
    prepare_basis_state(hf_state, data)
    allocate(block)

    allocate(phase_var)
    qpe(lambda: walk_operator(be_hamiltonian, [data, block]), phase_var)
    drop(data)


qprog_qpe_naive = synthesize(main)
show(qprog_qpe_naive)
Output:

Quantum program link: https://platform.classiq.io/circuit/38VyDBs3xcSfRT772EkonhMeqwA
  


results_qpe_naive = execute(qprog_qpe_naive).get_sample_result()


df_qpe_naive = results_qpe_naive.dataframe
print(f"Classical solution:, {classical_sol} Ha")
post_processed_df_qpe_naive = get_qpe_walk_result(
    df_qpe_naive, be_scaling, to_plot=True
)
Output:

Classical solution:, -1.1373060357533995 Ha
  probability to measure the block variable at state zero: 0.50732421875

  Energy with maximal probability: -1.102846988772674 Ha
output

​
Optimized QPE design for the walk operator

We construct the model design according to Ref. [1], shown in the figure below, with RLR_{\mathcal{L}} being the reflection around zero operation, and χm\chi_m is taken as the usual Hadamard transform for the phase initialization: Screenshot 2025-12-12 at 16.20.32.png 
@qfunc
def qpe_on_walk(
    block_encoding: QCallable[BlockEncodedState],
    state: BlockEncodedState,
    phase_var: QArray,
) -> None:

    hadamard_transform(phase_var)
    control(
        phase_var[0],
        lambda: walk_operator(block_encoding, state),
    )
    repeat(
        count=phase_var.len - 1,
        iteration=lambda i: within_apply(
            lambda: control(
                phase_var[i + 1] == 0,
                lambda: control(state.block == 0, lambda: phase(pi)),
            ),
            lambda: power(2**i, lambda: walk_operator(block_encoding, state)),
        ),
    )
    invert(lambda: qft(phase_var))
We now construct the model, synthesize it, and retrieve the ground state of the molecule 
@qfunc
def main(
    block: Output[QNum[block_size]],
    phase_var: Output[QNum[QPE_SIZE, SIGNED, QPE_SIZE]],
) -> None:

    data = QNum(size=data_size)
    prepare_basis_state(hf_state, data)
    allocate(block)
    allocate(phase_var)

    qpe_on_walk(block_encoding=be_hamiltonian, state=[data, block], phase_var=phase_var)
    drop(data)


qprog_qpe_walk = synthesize(main)
show(qprog_qpe_walk)
Output:

Quantum program link: https://platform.classiq.io/circuit/38VyLEYiNJ1PXIRn5XDbpa2V6MN
  


results_qpe_walk = execute(qprog_qpe_walk).get_sample_result()


df_qpe_optimized = results_qpe_walk.dataframe
print(f"Classical solution:, {classical_sol} Ha")
post_processed_df_qpe_optimized = get_qpe_walk_result(
    df_qpe_optimized, be_scaling, to_plot=True
)
Output:

Classical solution:, -1.1373060357533995 Ha
  probability to measure the block variable at state zero: 0.47119140625

  Energy with maximal probability: -1.102846988772674 Ha
output 
print("For the naive QPE on the walk operator:")
print(f"depth: {qprog_qpe_naive.transpiled_circuit.depth}")
print("=" * 40)
print("For the optimized QPE on the walk operator:")
print(f"depth: {qprog_qpe_walk.transpiled_circuit.depth}")
print("=" * 40)
Output:

For the naive QPE on the walk operator:
  depth: 19686
  ========================================
  For the optimized QPE on the walk operator:
  depth: 4724
  ========================================
  


gs_naive = post_processed_df_qpe_naive.loc[
    post_processed_df_qpe_naive["probability"].idxmax(), "post_processed_phases"
]
gs_optimized = post_processed_df_qpe_optimized.loc[
    post_processed_df_qpe_optimized["probability"].idxmax(), "post_processed_phases"
]
If ΔλW=1/2QPE−SIZE\Delta \lambda_W = 1/2^{\rm QPE-SIZE}, then Δλ=Δ(cos⁡(2πλW))s=2πssin⁡(2πλW)ΔλW=2πssin⁡(2π2−QPE−SIZE)2−QPE−SIZE\Delta \lambda = \Delta(\cos(2\pi\lambda_W))s = 2\pi s \sin(2\pi\lambda_W) \Delta\lambda_W = 2\pi s \sin(2\pi 2^{-\rm QPE-SIZE}) 2^{-\rm QPE-SIZE} 
qpe_err = (
    2 * np.pi * be_scaling * np.sin(2 * np.pi * (1 / 2**QPE_SIZE)) * (1 / 2**QPE_SIZE)
)


print(f"QPE error: {qpe_err}")
Output:

QPE error: 0.07603996508423008
  


assert np.abs(gs_naive - classical_sol) < qpe_err
assert np.abs(gs_optimized - classical_sol) < qpe_err

​
References

[1] R. Babbush et. al., Encoding Electronic Spectra in Quantum Circuits with Linear T Complexity. https://arxiv.org/abs/1805.03662 (2018) [2] Szegedy, M. , Quantum speed-up of Markov chain based algorithms. In 5th Annual IEEE Symposium on Foundations of Computer Science (2004) [3] Lin, L., Lecture notes on quantum algorithms for scientific computation. arXiv:2201.08309 quant-ph (2022)