Tensor hypercontraction block encoding

What it fills

• Block-encode a matrix

  Wrap an operator $A$ inside a larger unitary $U$ so that $A/α$ sits in $U$'s top-left block, giving every routine above it one uniform way to touch the matrix. The subnormalization $α$ and the ancilla count are outputs of this layer, not free parameters.

When it applies

Lee et al. state the contribution as a circuit and a complexity, not as a solve: "We describe quantum circuits with only Õ(N) Toffoli complexity that block encode the spectra of quantum chemistry Hamiltonians in a basis of $N$ arbitrary (e.g., molecular) orbitals" (abstract). What a consumer does with it is stated as an option rather than as this paper's work — "With O(λ/ε) repetitions of these circuits one can use phase estimation to sample in the molecular eigenbasis" — and §II is titled "Tensor Hypercontraction Representations for Quantum Simulation", a representation rather than an algorithm. The mechanism is the one its siblings on this slot use: "Our approach to encoding the eigenspectra of the THC representation... will use the linear combination of unitaries (LCU) query model", with the coefficient loading done by "a three step procedure where we first prepare an equal superposition over the μ and ν registers, then perform coherent alias sampling, then swap" — which is why this node steps into `state-preparation`.

Requires

These do not move the route along. The method needs each of them alongside its own work, and the cost of getting them is part of what the method costs.

• Prepare an input state

  Map $|0\ldots0⟩$ to a state whose amplitudes are proportional to a specified vector $b$, to within $\varepsilon$. The cost is set by which description of $b$ you hold, not by the algorithm that consumes it.

  THC's PREPARE builds the coefficient state over the factor-index registers $\mu,\nu \le M$ (the THC rank) and the one-body index $\ell \le N/2$, not over the $N$-orbital or $N^4$ two-electron-integral indices — that register choice is where the saving over a generic dense load comes from. The target is $\frac{1}{\sqrt{\lambda}}\Big(|+\rangle|+\rangle\sum_{\ell=1}^{N/2}\sqrt{|t_{\ell}|}\,|\ell\rangle|M+1\rangle+\frac{1}{\sqrt{2}}\sum_{\mu,\nu=1}^{M}\sqrt{|\zeta_{\mu\nu}|}\,|\mu\rangle|\nu\rangle\Big)$, built from $d = N/2 + M(M+1)/2$ coefficients in three steps: Hadamards and inequality tests give an equal superposition on $\mu \le \nu \le M+1$; a QROM outputs alternate values, sign qubits and a keep register for coherent alias sampling; then the $\mu$ and $\nu$ registers are swapped, controlled on an ancilla in $|+\rangle$ and on $\nu \ne M+1$. Each qubitized walk step invokes this preparation once and its inverse once, the two entering the per-step Toffoli count beside select's.

Example

given  quantum chemistry Hamiltonian H on N spin orbitals, two-electron
       integrals V_pqrs (Eq. 1)
       a THC factorization, computed OFFLINE and classically:
           G_pqrs = sum_{mu,nu=1..M} chi_p^(mu) chi_q^(mu) zeta_{mu,nu}
                                      chi_r^(nu) chi_s^(nu)             (Eq. 4)
           chi^(mu) taken normalized for each mu WLOG -- constant
               factors are absorbed into zeta_{mu,nu}
           M is fit per molecule by an L2-norm minimization (ISDF initial
               guess, then L-BFGS-B, then AdaGrad), NOT chosen by a proven
               bound                                             (Sec. II B)
       target energy accuracy epsilon, split as eps_pea + eps_thc <= epsilon
       between the phase-estimation error and the THC truncation error (Eq. 25)

requires  the "obvious" block encoding -- qubitizing G directly in the
    p,q,r,s orbital basis -- is NOT the construction below; that variant is
    relegated to Appendix E because its 1-norm
        lambda_thc = sum_{p,q,r,s} | sum_{mu,nu} chi_p^(mu) chi_q^(mu)
                                      zeta_{mu,nu} chi_r^(nu) chi_s^(nu) |  (Eq. 9)
    puts the sum over mu,nu OUTSIDE the absolute value, so "lambda_thc is
    much larger than even lambda_V, the 1-norm of the original Hamiltonian"
    (Sec. II C) -- avoiding that blow-up is the reason for stage 1 below

# --- stage 1: rotate into a non-orthogonal auxiliary basis where the ---
# --- Coulomb operator is diagonal in the number operators --------------
c_mu,sigma^dagger = sum_p chi_p^(mu) a_p,sigma^dagger                  (Eq. 10)
    # c_mu is NOT a proper fermionic mode -- non-unitary, non-orthogonal,
    #   a projection into an M-dimensional (per spin) auxiliary space,
    #   roughly 4-10x the number of orbital qubits                 (Sec. II C)
G = (1/2) sum_{alpha,beta in up,down} sum_{mu,nu=1..M}
        zeta_{mu,nu} n_mu,alpha n_nu,beta                              (Eq. 12)
    # "provides a diagonal form of the Coulomb operator" -- diagonal in
    #   n_mu,sigma = c_mu,sigma^dagger c_mu,sigma, in the ENLARGED basis,
    #   not diagonal in the original orbital basis                (Sec. II C)
substituting n = (1 - Z)/2 splits G into an energy shift, a residual
    ONE-body piece T^(2->1), and the two-body Z-Z piece            (Eq. 15, 16)
T' = T + T^(2->1),  T'_pq = T_pq + sum_{r=1..N/2} V_pqrr                (Eq. 18)
    # the exact integrals V_pqrr are used here in place of the fitted
    #   G_pqrr, so the THC approximation enters the TWO-body term only
    #   -- the paper's own choice                                      (Eq. 17)
diagonalize T' -> eigenvalues t_l, l = 1..N/2, rotations U_T,l         (Eq. 19)
H = -(1/2) sum_sigma sum_{l=1..N/2} t_l  U_T,l^dagger Z_1,sigma U_T,l
    + (1/8) sum_{alpha,beta} sum_{mu,nu=1..M}
        zeta_{mu,nu}  U_mu^dagger Z_1,alpha U_mu  U_nu^dagger Z_1,beta U_nu   (Eq. 26)
    # every term is a single Z sandwiched between basis-rotation unitaries
    #   U -- this is the form SELECT below implements
lambda_zeta = (1/2) sum_{mu,nu} |zeta_{mu,nu}|                         (Eq. 13)
    # CAUTION: lambda_zeta is only "a contribution from two-body terms",
    #   NOT the walk's 1-norm; it is what replaces lambda_thc of Eq. 9
lambda = sum_{l=1..N/2} |t_l| + lambda_zeta = O(lambda_zeta)           (Eq. 20)
    # lambda -- one-body plus two-body -- is the quantity PREPARE
    #   normalizes by and phase estimation pays for; the paper states its
    #   relation to lambda_zeta as an ORDER bound, not an equality, and it
    #   would be an equality only if every t_l vanished              (Eq. 20)

# --- stage 2: PREPARE -- coefficient state over mu, nu, l, NOT p,q,r,s --
target |PREPARE> =
    (1/sqrt(lambda)) ( |+>|+> sum_{l=1..N/2} sqrt(|t_l|) |l>|M+1>
        + (1/sqrt(2)) sum_{mu,nu=1..M} sqrt(|zeta_{mu,nu}|) |mu>|nu> )    (Eq. 27)
    # the SAME lambda as Eq. 20 -- both the one-body t_l and the two-body
    #   zeta_{mu,nu} coefficients live in this one state
    # d = N/2 + M(M+1)/2 coefficients total -- register size scales with
    #   the THC rank M, not with N^4                                (Eq. 28)
2a: Hadamards + inequality tests (nu<=M+1, mu<=nu, mu<=N/2 if nu=M+1)
    -> equal superposition over mu<=nu<=M+1
    # cost 10*ceil(log(M+1)) + 2*b_r - 9 Toffolis, b_r ~ 7 bits of
    #   rotation precision                                             (Fig. 3)
2b: QROM outputs alternate values, sign qubits, and an aleph-bit keep
    register -- coherent alias sampling                                (Eq. 30, 31)
    # "is equivalent to discretizing the squared amplitudes to the
    #   nearest 1/(2^aleph d)" -- the loaded coefficients are APPROXIMATED
    #   here; this is one of the paper's four named error sources
2c: swap the mu, nu registers, controlled on an ancilla in |+> and on
    nu != M+1                                                          (Fig. 4)

return  0  # PREPARE alone does not return a value -- it hands a prepared
        #   register to SELECT below

# --- stage 3: SELECT -- apply the Z through N/2 Givens rotations per ---
# --- factor, not a full N-qubit basis change ----------------------------
# "only N/2 Givens rotations are needed, instead of O(N^2) if all N basis
#   vectors were being rotated at once" -- because the register controls
#   rotating ONE number operator's basis at a time      (Eq. 51 of [10], Sec. II D)
for register in (mu, nu):                        # performed twice
    swap spin-up/spin-down halves, controlled on the spin qubit    # cost N/2
    QROM out N/2 rotation angles indexed by this register
        # mu-step: M + N/2 - 2 Toffolis (one-body and two-body share it)
        # nu-step: M - 2 Toffolis (two-body term only)
    apply the rotation into a phase-gradient register, cost N*(i-2)
        # i = bits of angle precision
    apply Z_1 on the rotated qubit                     # Clifford, no Toffoli
    uncompute the rotation, cost N*(i-2)
    erase the QROM by a variable-radix scheme:
        mu-register cost ceil(M/k_r1) + ceil(N/(2*k_r1)) + k_r1        (Eq. 34)
        nu-register cost ceil(M/k_r2) + k_r2                            (Eq. 35)
    swap back                                                       # cost N/2
# "the select operation needs to be made self-inverse, and the operation
#   as depicted is not self-inverse" -- fixed by an extra NOT on the
#   mu/nu swap-control qubit, verified algebraically term by term (Eqs. 36-42)

# --- stage 4: one qubitized walk step ------------------------------------
W = reflection(about |PREPARE>|0>) . SELECT
    # a single controlled-Z on the all-zero |l> state realizes the
    #   reflection; SELECT is self-inverse by stage 3                  (Fig. 1)
alpha = lambda                # the block-encoding subnormalization, the
                              #   general LCU lambda = sum_l |omega_l|  (Eq. 22)
                              #   evaluated for this Hamiltonian        (Eq. 20)
per-step Toffoli count ~ C_S + C_P + C_P^dagger + log(L) + O(1)        (Eq. 23)

return  block encoding (W, alpha = lambda, ancilla count a)
# consumed by phase estimation with ceil(pi*lambda/(2*eps_pea))
#   repetitions of W -- lambda, not lambda_zeta                        (Eq. 23)

# end-to-end Toffoli complexity, repetitions folded in:
#   O~(N * lambda_zeta / epsilon), space O~(N) qubits -- this FORM is what
#   the paper proves; it is NOT the same claim as a specific exponent on N
#   (lambda_zeta appears here rather than lambda only because
#    lambda = O(lambda_zeta) by Eq. 20)
# the exponents actually quoted -- N^3.1 toward the continuum limit,
#   N^2.1 toward the thermodynamic limit -- come from linear fits on a log
#   scale to HYDROGEN CHAIN data, not a proof
# likewise the FeMoco totals (M=350: 2142 qubits, 5.3e9 Toffolis on 108
#   spin orbitals; M=450: 2196 qubits, 3.2e10 Toffolis on 152 spin
#   orbitals, both at eps=0.0016 Hartree) are NUMERICAL results for two
#   specific molecules with a numerically FIT M, not an asymptotic bound
#   on M -- the paper states only that M = O(N polylog(1/eps_thc)) is
#   empirically observed, and calls its own resource numbers "upper
#   bounds on the cost of the most efficient possible implementations"
# preparation of the initial state |phi> that phase estimation walks on
#   is explicitly assumed elsewhere in the paper, not part of this method

Cost, as the source states it

The conclusions state Toffoli complexity $\widetilde{\mathcal O}(N\lambda_\zeta/\epsilon)$ and space complexity $\widetilde{\mathcal O}(N)$ qubits, with $N$ the number of spin-orbitals, $\epsilon$ the target energy accuracy, and $\lambda_\zeta = \frac{1}{2}\sum_{\mu,\nu}|\zeta_{\mu\nu}|$ the one-norm of the THC core tensor over its $M$ auxiliary indices; the walk's normalization $\lambda$, one-body plus two-body, is $\mathcal O(\lambda_\zeta)$. That is the phase-estimation total with the repetitions already multiplied in, not the per-circuit encoding price. Nothing here bounds $M$: it is fitted numerically per molecule, and the $M$ between $2N$ and $3N$ the paper mentions is relayed from prior work's classical $\mathcal O(N^4)$ fitting algorithm. Measured, not proved, at $\epsilon = 0.0016$ Hartree: $M=350$ gives 2,142 logical qubits and $5.3\times10^9$ Toffolis on the 108-spin-orbital Reiher FeMoco Hamiltonian, $M=450$ gives 2,196 qubits and $3.2\times10^{10}$ Toffolis on the 152-spin-orbital Li FeMoco Hamiltonian — both end-to-end phase-estimation totals.

Implementations

• Qualtran's `PrepareTHC` and `SelectTHC` bloqs

  • Even more efficient quantum computations of chemistry through tensor hypercontraction
  • Expressing and Analyzing Quantum Algorithms with Qualtran

  Qualtran is Google Quantum AI's Python library for expressing and analyzing fault-tolerant quantum algorithms. It carries this record's PREPARE/SELECT pair as two classes in `qualtran/bloqs/chemistry/thc/`, distinct from the library's generic Pauli-LCU and sparse-access block-encoding bloqs — this module is specific to the THC Hamiltonian's own registers ($\mu$, $\nu$, the one-body index) rather than a generic coefficient vector. A separate function, `get_walk_operator_for_thc_ham`, wires the two into a `SelectBlockEncoding` and then a `QubitizationWalkOperator` — the walk-operator layer that the sibling `qubitization-simulation` record's own Qualtran entry already documents one step further downstream, so this entry is confined to what that one does not cover: the PREPARE and SELECT circuits themselves.

  `UniformSuperpositionTHC` builds the equal superposition over $\mu \le \nu \le M+1$ by Hadamards, an amplitude-amplifying Ry rotation and a chain of comparators (`LessThanConstant`, `LessThanEqual`, `EqualsAConstant`, `GreaterThanConstant`), stated by its own docstring to cost "$10 \log(M+1) + 2 b_r - 9$" Toffolis — the same closed form this record's own pseudocode gives for that stage. `PrepareTHC.from_hamiltonian_coeffs` takes the one-body eigenvalues $t_\ell$, THC leaf tensor and central tensor as plain arrays, flattens the upper triangle of $\zeta$ plus $t_\ell$, and calls `preprocess_probabilities_for_reversible_sampling` for the alt/keep values of coherent alias sampling; a comment in the same method states it computes $\lambda$ "using the formula from the reference / OpenFermion: resource_estimates.thc.compute_lambda_thc", and the code beside it visibly follows Eq. 11/12 (normalizing $\zeta$), Eq. 19 for the one-body sum $\lambda_T$ and Eq. 20 for the two-body sum $\lambda_z$; the two are then added into `sum_of_l1_coeffs` on the next line with no equation number attached to that total. `SelectTHC.build_composite_bloq` assembles the system register out of `THCRotations` (applied, then a controlled-$Z$ via `ApplyControlledZs`, then un-applied) sandwiched between controlled spin-swaps, with the docstring noting explicitly that the $\mu/\nu$ swap is "NOT performed as part of SELECT as they're acounted for during Prepare" — a register-accounting choice this record's own pseudocode also makes. `THCRotations` itself, the Givens-rotation network that is stage 3 of this record's theory, has no `build_composite_bloq`: its own docstring calls it "a placeholder waiting for an actual implementation", and its cost is supplied only through `build_call_graph`, an analytic Toffoli count taken "from listings on page 17 of Ref. [1]". `PrepareTHC`, by contrast, does carry a full `build_composite_bloq`.

  No molecule and no physical dataset. The registered `bloq_example`s build small numeric instances: `_thc_uni` takes `num_mu=10, num_spin_orb=4`; `_thc_prep` calls a test helper, `build_random_test_integrals(num_mu=8, num_spat=4, seed=7)`, which draws a symmetrized random one-body matrix and random $\zeta$, $\eta$ matrices from `np.random.RandomState(7)`, then builds `PrepareTHC.from_hamiltonian_coeffs(..., num_bits_state_prep=8, log_block_size=2)`; `_thc_sel` builds a bare `SelectTHC(num_mu=10, num_spin_orb=8, num_bits_theta=12, keep_bitsize=10)` with no Hamiltonian data at all. None of the three reaches molecular scale.

  https://github.com/quantumlib/Qualtran, Python, Apache-2.0. `qualtran/bloqs/chemistry/thc/prepare.py` (classes `UniformSuperpositionTHC`, `PrepareTHC`) and `qualtran/bloqs/chemistry/thc/select_bloq.py` (classes `THCRotations`, `SelectTHC`), glued by `qualtran/bloqs/chemistry/thc/walk_operator.py`'s `get_walk_operator_for_thc_ham`. A notebook, `qualtran/bloqs/chemistry/thc/thc.ipynb`, sits beside the modules. Files read from the repository's default branch `main` on 2026-08-27.

  Classical simulation of the circuit's own bookkeeping, not a molecular run. `test_prepare_alt_keep_vals` in `prepare_test.py` reconstructs the target probability distribution from the alt/keep values `PrepareTHC` computes and checks it against the true normalized $|\zeta|$/$|t_\ell|$ weights to tolerance `eps = 2**-mu / len(flat_data)`, for three parameter triples `(num_mu, num_spat, mu)` of `(10, 4, 10)`, `(40, 10, 17)`, `(72, 31, 27)`, each built from `build_random_test_integrals(..., seed=7)`. `test_prepare_qrom_counts` checks the Toffoli count `PrepareTHC` attributes to data loading against the count its own `QROAMClean` bloq reports for the same data. `bloq_autotester` — run over `_thc_uni`, `_thc_prep` and `_thc_sel` — checks each bloq's decomposition is internally consistent (`assert_valid_bloq_decomposition`) and its declared costs match a call-graph count. None of this executes on hardware or at a molecular scale; the end-to-end FeMoco-scale Toffoli/qubit totals this record's own `cost` field quotes are not reproduced by any test in this module.

• OpenFermion's `resource_estimates.thc` cost and factorization module

  • Even more efficient quantum computations of chemistry through tensor hypercontraction

  OpenFermion is Google Quantum AI's electronic-structure-to-quantum-circuits package. Its `resource_estimates` module's own README lists "the [tensor hypercontraction](https://arxiv.org/abs/2011.03494) (THC) method" by name, linking directly to the paper this record cites, alongside the single- and double-factorization methods on the same slot. The README states the dependency shape up front: the module needs PySCF and "is not installed by default" (`pip install openfermion[resources]`), and "For THC factorization, it also requires BTAS and the PyBTAS wrapper, which require their own installation + depends" — so the factorization routine below is opt-in even within an OpenFermion install.

  The README names three subroutines per factorization: `factorize()`, `compute_lambda()`, `compute_cost()`; for THC these are `thc_via_cp3` (aliased `factorize`), `compute_lambda`, and `compute_cost`. `compute_cost(n, lam, dE, chi, beta, M, stps)` reproduces this record's own Toffoli-count derivation term by term: its `cp1 = 2 * (10 * nM + 2 * br - 9)` is twice the record's own "cost 10*ceil(log(M+1)) + 2*b_r - 9 Toffolis" for the equal-superposition stage (the factor of 2 is forward plus inverse), and later terms are individually commented by the equation they implement, e.g. `cs2a`/`cs2b` as the two QROMs of the rotation angles and `costref` as the reflection. `compute_lambda(pyscf_mf, etaPp, MPQ)` is headed "Compute lambdas for THC according to PRX QUANTUM 2, 030305 (2021) Section II. D." — the published form of the same paper — and its own inline comments tag `lambda_z = np.sum(np.abs(MPQ_normalized)) * 0.5` as Eq. 13, `lambda_T = np.sum(np.abs(e))` as Eq. 19, and their sum as Eq. 20, matching this record's own $\lambda_\zeta$ and $\lambda$. `thc_via_cp3` is the classical rank-reduction step, but it is not the ISDF route this record's pseudocode attributes to the paper's Section II B: it SVDs the two-electron-integral matrix, then runs a CP3 tensor decomposition through the external `pybtas` package (`import pybtas`, raising `ImportError` if absent) on the truncated Cholesky factor, and only optionally refines the result with L-BFGS-B (`lbfgsb_opt_thc_l2reg`, under `perform_bfgs_opt`) — a different starting point than the paper's own ISDF-initial-guess procedure, even though the L-BFGS-B refinement stage is shared.

  No physical dataset for `compute_cost`, whose inputs are the scalar parameters $n$, $\lambda$, $\delta E$, $\chi$, $\beta$, $M$ and a step count. `compute_lambda` takes a PySCF mean-field object plus the THC leaf and central tensors as input — molecular integrals the user supplies via PySCF, not data shipped with the module. `factorize` (`thc_via_cp3`) instead takes only a bare four-index electron-repulsion-integral tensor and a target THC rank; the leaf and central tensors it needs are among its own return values, not something supplied to it.

  https://github.com/quantumlib/OpenFermion, Python, Apache-2.0. `src/openfermion/resource_estimates/thc/compute_cost_thc.py` (`compute_cost`), `compute_lambda_thc.py` (`compute_lambda`), and `factorize_thc.py` (`thc_via_cp3`, exposed as `factorize`), all re-exported from `openfermion.resource_estimates.thc`. `compute_lambda` and `factorize` are gated behind `HAVE_DEPS_FOR_RESOURCE_ESTIMATES` in `__init__.py`, consistent with the README's optional-dependency note. Files read from the repository's default branch `main` on 2026-08-27.

  Two exact-value regression tests in `compute_cost_thc_test.py`, both reproducing published per-molecule totals rather than a hardware run. `test_reiher_thc` calls `compute_cost` twice — once to size the step count, once for the final answer — with `N=108, LAM=306.3, DE=0.001, CHI=10, BETA=16, THC_DIM=350` and asserts the return is exactly `(10912, 5250145120, 2142)`, i.e. 10,912 Toffolis per step, 5,250,145,120 Toffolis total, 2,142 logical qubits. `test_li_thc` does the same with `N=152, LAM=1201.5, DE=0.001, CHI=10, BETA=20, THC_DIM=450`, asserting `(16923, 31938980976, 2196)` — 31,938,980,976 Toffolis total, 2,196 qubits. Both magnitudes and both qubit counts land where this record's own `cost` field puts the Reiher and Li FeMoco totals (2,142 and 2,196 qubits exactly; the Toffoli totals round to the same first two significant figures, $5.3\times10^{9}$ and $3.2\times10^{10}$), though the test's `DE=0.001` is not stated here to be the same quantity as the $\epsilon=0.0016$ Hartree this record's own `cost` field names, since the paper splits its target accuracy into a phase-estimation part and a THC-truncation part and it is the former, not the total, that a Toffoli-cost function like this one would consume.

What it needs

Every step this method names is listed under Requires above. It walks its own span in one hop and calls out to the rest — that is a fact about the recorded route, not a claim that the span is simple.

Other ways to fill the same slot

Different approaches

• Sparse-access oracle construction

  Given row and column index oracles $O_r$, $O_c$ and an entry oracle $O_A$, prepare uniform superpositions over the sparsity pattern, rotate an ancilla by arcsin of each entry, and swap registers to leave $A$ in the flagged block. This is the standard construction behind the sparse-Hamiltonian line, formalized as a block-encoding by Gilyén, Su, Low and Wiebe.

• Block-encoding from a Pauli decomposition

  Write $A = Σ_j c_j P_j$ over Pauli strings; PREPARE loads amplitudes proportional to $sqrt(|c_j|)$ into an ancilla register, SELECT applies the controlled Pauli strings, and PREPARE$†$ unprepares, leaving $A/||c||_1$ in the block flagged by the all-zeros ancilla. This is the input model chemistry and lattice Hamiltonians supply for free.

• FABLE approximate circuit construction

  Build the block-encoding directly from uniformly controlled Ry (magnitude) and Rz (phase) rotations between Hadamards and a SWAP, with no oracle assumption at all, then threshold the rotation angles and cancel the resulting CNOT chains to compress the circuit.

In the Atlas

• Tensor hypercontraction block encoding for quantum chemistry

  Represent the spectrum of a quantum chemistry Hamiltonian, given in an arbitrary (for example molecular) orbital basis, as a block-encoded quantum circuit cheap enough to support phase estimation of a molecular eigenvalue.

Sources

• Even more efficient quantum computations of chemistry through tensor hypercontraction↗Joonho Lee, Dominic W. Berry, Craig Gidney, William J. Huggins, Jarrod R. McClean, Nathan Wiebe, Ryan Babbush · 2020theory yes · simulation ? · hardware no