Skip to main content

View on GitHub

Open this notebook in GitHub to run it yourself

Even more efficient quantum computations of chemistry through tensor hypercontraction

The notebook implements the core block of the algorithm described in https://arxiv.org/pdf/2011.03494
import itertools
from typing import Iterator, Sequence

import numpy as np

from classiq import *
from classiq.interface.generator.model.preferences.preferences import (
    OptimizationLevel,
    TranspilationOption,
)
from classiq.qmod.qmod_variable import QVar
from classiq.qmod.symbolic import pi

Variables Definition

This section defines the variables used in the preprocessing code.
All theoretical motivation, physical meaning, and algorithmic justification are given in the reference article; here, we only document the role of each variable in the code.

Resolution and spin parameters

  • N_EPSILON
    Number of qubits used to represent the quantum number resolution of the result.
    Controls numerical precision in state preparation and phase-dependent quantities.
  • N_UD
    Number of spin-orbitals per spin sector (up/down).
counter = itertools.count()

N_EPSILON = 4  # This should be 10
N_UD = 3  # This should be 54
n_mn = (3 * 2 * N_UD - 1).bit_length()
m = max(2**n_mn - 1, 2)
n_d = (N_UD + m * (m // 2) - 1).bit_length()
n_aleph = 14 + N_EPSILON
n_beth = N_EPSILON
n_t = N_EPSILON
t_max = max(2**n_t - 1, 2)
n_rot = n_beth * N_UD
# n_ell = 2 * n_mn + 4 + n_d + 1 + 2 * n_mn + 2 * n_aleph + 1
# n_psi = 2 * N_UD + 2

Randomized preprocessing data

The following variables are randomly generated placeholders with the correct shapes and bit-widths required by the algorithm:
  • R_BR
    Random rotation angle used in state preparation.
  • MU_ALT, NU_ALT
    Boolean tables encoding alternate auxiliary indices.
  • KEEP
    Boolean table encoding keep / discard decisions for coherent alias sampling.
  • THETA
    Boolean array encoding sign or phase information.
  • ROT_MU, ROT_NU
    Boolean tables encoding rotation data for auxiliary indices μ and ν.

Precomputed amplitude data

  • SINE_AMPLITUDES
    Discrete sine amplitude table used for state preparation.
    Length: 2**n_t.

Note

  • All numerical values are synthetic and used only to validate register sizing, indexing, and data flow.
R_BR = np.random.rand() * 2 * np.pi
MU_ALT = np.random.choice((True, False), (m, n_mn)).tolist()
NU_ALT = np.random.choice((True, False), (m, n_mn)).tolist()
KEEP = np.random.choice((True, False), (m, n_aleph)).tolist()
THETA = np.random.choice((True, False), m).tolist()
ROT_MU = np.random.choice((True, False), (m, n_rot)).tolist()
ROT_NU = np.random.choice((True, False), (m, n_rot)).tolist()
SINE_AMPLITUDES = [
    np.sqrt(2 / (2**n_t + 1)) * np.sin(np.pi * (i + 1) / (2**n_t + 1))
    for i in range(2**n_t)
]

High-Level Role of the Quantum Registers

This section explains the conceptual role of each QStruct in the context of the algorithm described in the article.

ELL — Control and State-Preparation Register

ELL groups all registers used to index Hamiltonian terms, prepare amplitudes, and control conditional operations in the linear-combination-of-unitaries (LCU) / qubitization framework. At a high level, this register:
  • Encodes auxiliary indices (μ, ν) labeling terms in the Hamiltonian expansion
  • Stores temporary control and success flags used during preparation and uncomputation
  • Holds auxiliary data needed for phase, sign, and rotation selection
Conceptually, ELL corresponds to the LCU registers used in:
  • the prepare oracle.
  • the select oracle.
It does not represent the physical system itself, but rather the algorithmic machinery required to construct the qubitized quantum walk.

PSI — System (Fermionic State) Register

PSI represents the quantum state of the simulated fermionic system. At a high level, this register:
  • Encodes the occupation of spin-orbitals for both spin sectors (up and down).
  • Stores the many-body electronic state on which the Hamiltonian acts
  • Serves as the target of controlled operations generated by the select oracle
This register corresponds directly to the fermionic system register on which operators such as ZZ, number operators, and basis rotations act. The additional auxiliary qubits (plus1, plus2) are used to control superpositions over spin components.

Summary

  • ELL:
    Algorithmic control, indexing, and state-preparation workspace for qubitization.
  • PSI:
    Physical system register representing the fermionic many-body state.
Together, these registers are the variables of the block-encoded Hamiltonian and quantum walk construction described in the article.
class ELL(QStruct):
    mu: QNum[n_mn]
    nu: QNum[n_mn]
    m1: QBit
    succ: QBit
    r_br: QBit
    theta: QBit
    s: QNum[n_d]
    theta_alt: QBit
    mu_alt: QNum[n_mn]
    nu_alt: QNum[n_mn]
    keep: QNum[n_aleph]
    sigma: QNum[n_aleph]
    aux: QBit


class PSI(QStruct):
    down: QArray[QBit, N_UD]
    up: QArray[QBit, N_UD]
    plus1: QBit
    plus2: QBit

Helper functions

def new_qubit(name: str = "auto") -> QBit:
    return QBit(f"{name}{next(counter)}")


@qfunc
def array_cast(qvar: QArray, action: QCallable[QArray]) -> None:
    action(qvar)


def in_iteration(
    ctrl: QBit, x: QArray[QBit], bit: int, length: int, table: Iterator[QCallable]
) -> None:
    if length == 1:
        next(table)(ctrl)
        return
    half = 2 ** ((length - 1).bit_length() - 1)
    aux = new_qubit("aux")
    allocate(aux)
    # Q: Why don't we have control strings?
    within_apply(
        within=lambda: X(x[bit]),
        apply=lambda: CCX([ctrl, x[bit]], aux),
    )
    in_iteration(aux, x, bit + 1, half, table)
    CX(ctrl, aux)
    in_iteration(aux, x, bit + 1, length - half, table)
    CCX([ctrl, x[bit]], aux)
    free(aux)

Prepare

  • Part 1
Prepare 1
@qperm
def assign_aux(
    mu: QNum[n_mn],
    nu: QNum[n_mn],
    aux1: Output[QBit],
    aux2: Output[QBit],
    aux3: Output[QBit],
    aux4: Output[QBit],
) -> None:
    assign(nu <= m, aux1)
    assign(mu <= nu, aux2)
    assign(mu > n_mn / 2, aux3)
    assign(1, aux4)


@qfunc
def create_equal_superposition(ell: ELL) -> None:
    aux1 = QBit()
    aux2 = QBit()
    aux3 = QBit()
    aux4 = QBit()
    apply_to_all(H, ell.mu)
    apply_to_all(H, ell.nu)
    within_apply(
        within=lambda: (
            RY(R_BR, ell.r_br),
            inplace_xor(ell.nu == m + 1, ell.m1),
            assign_aux(ell.mu, ell.nu, aux1, aux2, aux3, aux4),
            CCX([ell.m1, aux3], aux4),
        ),
        apply=lambda: control([aux1, aux2, aux4], lambda: Z(ell.r_br)),
    )
    within_apply(
        within=lambda: (
            apply_to_all(H, ell.mu),
            apply_to_all(H, ell.nu),
        ),
        apply=lambda: control([ell.mu, ell.nu], lambda: Z(ell.r_br)),
    ),
    inplace_xor(ell.nu == m + 1, ell.m1)
    within_apply(
        within=lambda: (
            assign_aux(ell.mu, ell.nu, aux1, aux2, aux3, aux4),
            CCX([ell.m1, aux3], aux4),
        ),
        apply=lambda: control([aux1, aux2, aux4], lambda: X(ell.succ)),
    )


@qfunc
def main(ell: Output[ELL]):
    allocate(ell)
    create_equal_superposition(ell)


qprog_prepare_part1 = synthesize(main)

def print_qprog_stats(qprog: QuantumProgram) -> None:
    print(f"Total number of qubits is {qprog.data.width}")
    print(f"Total depth is {qprog.transpiled_circuit.depth}")

print_qprog_stats(qprog_prepare_part1)
show(qprog_prepare_part1)
Output:

Total number of qubits is 82
  Total depth is 3842
  Quantum program link: https://platform.classiq.io/circuit/36pzAyiSqMFc6yacBPxwpcme00w
  

Prepare

  • Part 2
prepare_2
def in_uncontrolled(x: QArray[QBit], length: int, table: Iterator[QCallable]) -> None:
    if length == 1:
        raise NotImplementedError
    half = 2 ** ((length - 1).bit_length() - 1)
    within_apply(
        within=lambda: X(x[0]),
        apply=lambda: in_iteration(x[0], x, 1, half, table),
    )
    in_iteration(x[0], x, 1, length - half, table)


@qperm
def shift_right(x: Const[QNum], y: Output[QNum]) -> None:
    lsb = QBit()
    msbs = QNum("msbs", x.size - 1, False, 0)
    within_apply(
        within=lambda: bind(x, [msbs, lsb]),
        apply=lambda: assign(msbs, y),
    )


@qperm
def s_arithmetic(nu: Const[QNum], mu: Const[QNum], s: QNum) -> None:
    half_nu = QNum()
    within_apply(
        within=lambda: shift_right(nu, half_nu),
        apply=lambda: inplace_xor(nu * (half_nu + 1) + mu, s),
    )


@qperm
def single_qrom_access(
    ctrl: Const[QBit], target: QArray[QBit], data: CArray[CBool]
) -> None:
    repeat(data.len, lambda i: if_(data[i], lambda: CX(ctrl, target[i])))


def _yield_s(ell: ELL) -> Iterator[QCallable]:
    for i in range(m):
        yield lambda ctrl: (
            single_qrom_access(ctrl, ell.mu_alt, MU_ALT[i]),
            single_qrom_access(ctrl, ell.nu_alt, NU_ALT[i]),
            single_qrom_access(ctrl, ell.keep, KEEP[i]),
            single_qrom_access(ctrl, ell.theta, [THETA[i]]),
        )


@qperm
def _swap(state1: QArray[QBit], state2: QArray[QBit]) -> None:
    repeat(state1.len, lambda i: SWAP(state1[i], state2[i]))


@qperm
def in_s(ell: ELL) -> None:
    array_cast(ell.s, lambda s: in_uncontrolled(s, m, _yield_s(ell)))


@qfunc
def prepare_after_equal_superposition(ell: ELL) -> None:
    aux = QBit()
    s_arithmetic(ell.nu, ell.mu, ell.s)
    apply_to_all(H, ell.sigma)
    in_s(ell)
    within_apply(
        within=lambda: assign(ell.keep < ell.sigma, aux),
        apply=lambda: (
            CZ(ell.theta_alt, aux),
            within_apply(
                within=lambda: X(aux),
                apply=lambda: CZ(ell.theta, aux),
            ),
            control(aux, lambda: _swap(ell.mu, ell.mu_alt)),
            control(aux, lambda: _swap(ell.nu, ell.nu_alt)),
        ),
    ),
    H(ell.aux)
    # Q: Why don't we have control strings?
    within_apply(
        within=lambda: X(ell.m1),
        apply=lambda: control([ell.aux, ell.m1], lambda: _swap(ell.mu, ell.nu)),
    )


@qfunc
def main(ell: Output[ELL]):
    allocate(ell)
    prepare_after_equal_superposition(ell)


qprog_prepare_part2 = synthesize(main)

print_qprog_stats(qprog_prepare_part2)
show(qprog_prepare_part2)
Output:

Total number of qubits is 92
  Total depth is 5023
  Quantum program link: https://platform.classiq.io/circuit/36pzJPfsNDdGCbBap1mXDZ9yhwt
  

Prepare summary

@qfunc
def prepare(ell: ELL) -> None:
    create_equal_superposition(ell)
    prepare_after_equal_superposition(ell)

Select

select
def _yield_qrom(target: QVar, data: list[bool]) -> Iterator[QCallable]:
    for i in range(len(data)):
        yield lambda ctrl: single_qrom_access(ctrl, target, data[i])


@qfunc
def pauli_x_basis(q: QBit) -> None:
    H(q)


@qfunc
def pauli_y_basis(q: QBit) -> None:
    SDG(q)
    H(q)


@qperm
def _crzz(theta: CReal, ctrl: Const[QBit], target: Const[QArray[QBit, 2]]) -> None:
    within_apply(
        within=lambda: CX(target[0], target[1]),
        apply=lambda: CRZ(theta, ctrl, target[1]),
    )


@qperm
def digital_rotation(rot: Const[QArray], qs: Const[QArray[QBit, 2]]) -> None:
    repeat(rot.len, lambda i: _crzz(2 ** (-i) * pi / 2**rot.len, rot[i], qs))


@qfunc
def controlled_basis_change(rot: Const[QArray[QNum]], state: QArray[QBit]) -> None:
    repeat(
        rot.len - 1,
        lambda i: (
            within_apply(
                within=lambda: (
                    pauli_x_basis(state[i]),
                    pauli_y_basis(state[i + 1]),
                ),
                apply=lambda: digital_rotation(rot[i], [state[i], state[i + 1]]),
            ),
            within_apply(
                within=lambda: (
                    pauli_x_basis(state[i + 1]),
                    pauli_y_basis(state[i]),
                ),
                apply=lambda: invert(
                    lambda: digital_rotation(rot[i], [state[i], state[i + 1]])
                ),
            ),
        ),
    )


def in_controlled(
    ctrl: QBit, x: QArray[QBit], length: int, table: Iterator[QCallable]
) -> None:
    in_iteration(ctrl, x, 0, length, table)


@qperm
def in_mu(mu: QArray[QBit], m1: QBit, rot: QArray[QNum]) -> None:
    in_controlled(m1, mu, m, _yield_qrom(rot, ROT_MU))


@qperm
def in_nu(nu: QArray[QBit], rot: QArray[QNum]) -> None:
    in_uncontrolled(nu, m, _yield_qrom(rot, ROT_NU))


@qperm
def _ciz(ctrl: Const[QBit], target: Const[QBit]) -> None:
    S(ctrl)
    CZ(ctrl, target)


@qperm
def _cciz(ctrl1: Const[QBit], ctrl2: Const[QBit], target: Const[QBit]) -> None:
    control(ctrl1, lambda: S(ctrl2))
    control([ctrl1, ctrl2], lambda: Z(target))


@qfunc
def select(ell: ELL, psi: PSI) -> None:
    rot = QArray("rot", QNum[n_beth], N_UD)
    within_apply(
        within=lambda: (
            allocate(rot),
            in_mu(ell.mu, ell.m1, rot),
            control(psi.plus1, lambda: _swap(psi.down, psi.up)),
            invert(lambda: controlled_basis_change(rot, psi.down)),
        ),
        apply=lambda: _ciz(ell.succ, psi.down[0]),
    )
    within_apply(
        within=lambda: (
            allocate(rot),
            in_nu(ell.nu, rot),
            control(psi.plus2, lambda: _swap(psi.down, psi.up)),
            invert(lambda: controlled_basis_change(rot, psi.down)),
        ),
        # Q: Why don't we have control strings?
        apply=lambda: within_apply(
            within=lambda: X(ell.m1),
            apply=lambda: _cciz(ell.succ, ell.m1, psi.down[0]),
        ),
    )


@qfunc
def main(ell: Output[ELL], psi: Output[PSI]):
    allocate(ell)
    allocate(psi)
    select(ell, psi)


qprog_select = synthesize(main)

print_qprog_stats(qprog_select)
show(qprog_select)
Output:

Total number of qubits is 107
  Total depth is 3151
  Quantum program link: https://platform.classiq.io/circuit/36pzUmQOu8zcaiL1MoiDQkatOHv
  

Walk Operator

walk
@qfunc
def controlled_reflect(ctrl: QBit, ell: ELL) -> None:
    within_apply(
        within=lambda: invert(lambda: prepare(ell)),
        # Q: Why don't we have good control logic? This is a basic reflection.
        apply=lambda: within_apply(
            within=lambda: (
                X(ctrl),
                apply_to_all(X, ell),
            ),
            apply=lambda: array_cast(ell, lambda _ell: control(_ell, lambda: Z(ctrl))),
        ),
    )


def yield_w(ell: ELL, psi: PSI) -> Iterator[QCallable]:
    yield lambda aux: controlled_reflect(aux, ell)
    while True:
        yield lambda aux: (select(ell, psi), controlled_reflect(aux, ell))

Repeated Walk with Recursion

@qfunc
def init_psi(psi: PSI) -> None:
    H(psi.plus1)
    H(psi.plus2)


@qfunc
def prepare_xi(x: QArray[QBit]) -> None:
    inplace_prepare_amplitudes(SINE_AMPLITUDES, 0, x)


@qfunc
def in_t(t: QArray[QBit], ell: ELL, psi: PSI) -> None:
    in_uncontrolled(t, t_max, yield_w(ell, psi))


@qfunc
def main(t: Output[QArray[QBit, n_t]], ell: Output[ELL], psi: Output[PSI]) -> None:
    allocate(ell)
    allocate(t)
    allocate(psi)
    prepare_xi(t)
    init_psi(psi)
    in_t(t, ell, psi)
    invert(lambda: qft(t))


print(f"{N_EPSILON=}")
print(f"{N_UD=}")
Output:

N_EPSILON=4
  N_UD=3
  

qmod = create_model(main)
qmod = set_preferences(
    qmod,
    Preferences(
        transpilation_option=TranspilationOption.NONE,
        optimization_level=OptimizationLevel.NONE,
        debug_mode=False,
        timeout_seconds=1800,
    ),
)

qprog = synthesize(qmod)

show(qprog)
Output:

Quantum program link: https://platform.classiq.io/circuit/36pztXgMc7Ucpgw4g1QZ9MbHLZ0
  

print(f"Total number of qubits is {qprog.data.width=}")
Output:

Total number of qubits is qprog.data.width=114