Requires |a_ij| ≤ 1. The authors state the gate complexity is "bounded by O(N^2) gates with a modest prefactor of 2 for real-valued matrices (4 for complex-valued matrices), and a limited polylogarithmic overhead", and that "this gate complexity scales exponentially in the number of qubits for generic dense matrices". So this is a compression technique for structured matrices — Heisenberg, Hubbard, Laplacian — and not a way to load dense unstructured data cheaply. The paper states only that the circuit parameters can be easily generated for problems up to fifteen qubits.

An access model for $A$ — sparse-access oracles, a Pauli or LCU decomposition, a purification, or an explicit arithmetic description — plus a target precision $\varepsilon$.

1. Matrix you can query → Block-encodingsolve one Walsh-Hadamard for all angles

   Every entry is stored as a rotation amplitude by one oracle, $O_{A}|0\rangle|i\rangle|j\rangle = a_{ij}|0\rangle|i\rangle|j\rangle + \sqrt{1-|a_{ij}|^{2}}|1\rangle|i\rangle|j\rangle$, which for real $A$ is a controlled $R_{y}$ through $\theta_{ij} = \arccos(a_{ij})$. assumption: $|a_{ij}| \le 1$, the domain on which that arccosine is real — the paper states the bound as a condition on the matrix. Synthesising $N^{2}$ separate multi-controlled rotations would cost $O(N^{4})$, which the authors call "quadratically worse than the classical representation cost"; instead the $N^{2}$ angles are treated as one uniformly controlled rotation and the per-gate angles $\hat\theta$ come from a single linear system $\hat{H}^{\otimes 2n}P_{G}\hat\theta = \theta$ — a Walsh-Hadamard transform in Gray-code order, solvable classically in $O(N^{2}\log N^{2})$ — leaving $O(N^{2})$ one- and two-qubit gates. Hadamard layers either side and a SWAP project the row superposition back onto $\langle i|$, leaving $a_{ij}/2^{n}$. approximation: compression discards entries of $\hat\theta$ below a cutoff $\delta_{c}$, which sparsifies in the transformed angle domain — "a sparse representation of $\hat\theta$ does not typically mean that $\theta$ and $A$ are sparse in the usual sense".

   approximationassumption

A unitary $U$ on $s+a$ qubits, its subnormalization $\alpha$, and its ancilla/flag count $a$. Because $\lVert U\rVert = 1$, Gilyén, Su, Low and Wiebe's Definition 43 forces $\lVert A\rVert \le \alpha + \varepsilon$.

given  A with |a_ij| <= 1, and no oracle assumption at all

build the block-encoding circuit directly from:
    uniformly controlled Ry rotations    # magnitude
    uniformly controlled Rz rotations    # phase
    Hadamards and a SWAP, the rotations sitting between them

# Camps and Van Beeumen Theorem 1: that circuit is a (1/2^n, n+1)-block-
# encoding of A when |a_ij| <= 1

threshold the rotation angles at a cutoff compression threshold delta_c
cancel the CNOT chains the thresholding leaves behind

return the compressed circuit

# Camps and Van Beeumen Theorem 2: the cutoff delta_c gives a
# (1/2^n, n+1, N^3*delta_c)-block-encoding of A, up to third order in delta_c

# the 1/2^n above is transcribed as the paper states it, and contested on this
# record: it is inconsistent with the paper's own Definition 1, which sets
# A = alpha * A-tilde and therefore ||A||_2 <= alpha, the same orientation as
# Gilyen et al.'s Definition 43 and not its reciprocal, while the Theorem 1
# proof derives a block equal to a_ij/2^n. Under either paper's definition the
# subnormalization is alpha = 2^n = N. This applies to the arXiv text as
# served; a journal version may differ

# gate complexity, in the authors' words: bounded by O(N^2) gates with a
# modest prefactor of 2 for real-valued matrices (4 for complex-valued
# matrices), and a limited polylogarithmic overhead -- and this gate complexity
# scales exponentially in the number of qubits for generic dense matrices

# so this is a compression technique for structured matrices -- Heisenberg,
# Hubbard, Laplacian -- and not a way to load dense unstructured data cheaply

# the paper states only that the circuit parameters can be easily generated
# for problems up to fifteen qubits

Camps and Van Beeumen's Theorem 1 states the circuit is an "$(1/2^n, n+1)$-block-encoding" of $A$ when $|a_ij| \leq 1$, and Theorem 2 that a cutoff compression threshold $\delta_c$ gives an $(1/2^n, n+1, N^3\cdot\delta_c)$-block-encoding of $A$ up to third order in $\delta_c$.

The stated parameter $1/2^n$ is inconsistent with the paper's own Definition 1, which sets A = α·Ã and therefore $||A||_2 \leq \alpha$ — the same orientation as Gilyén et al.'s Definition 43, not its reciprocal — while the Theorem 1 proof derives a block equal to $a_ij/2^n$. Under either paper's definition the subnormalization is $\alpha = 2^n = N$. This applies to the arXiv text as served; a journal version may differ.

• FABLE reference implementation (Qiskit and QCLAB)

  • FABLE: Fast Approximate Quantum Circuits for Block-Encodings

    Daan Camps, Roel Van Beeumen · 2022

  About

  Camps and Van Beeumen built FABLE at Lawrence Berkeley National Laboratory — Camps at the National Energy Research Scientific Computing Center, Van Beeumen in the Applied Mathematics and Computational Research Division. The problem they state is that "Many of these algorithms achieve optimal complexity in terms of black box matrix oracle queries, but so far the problem of computing quantum circuit implementations for block-encodings of matrices has been under-appreciated" — "these algorithms" being the ones derived from the quantum singular value transformation, and "many", not all, of them. The artefact is the tool that closes that gap: a routine that takes a matrix and returns an explicit circuit of one- and two-qubit gates, plus the compression pass that shrinks it.

  Methods

  The release is two routines, and they differ in ways a caller meets on the first call, so each step below is attributed. Both rescale $A$ only when the infinity norm of the vectorised matrix exceeds 1: the Python routine sets $\alpha = \max_{ij}|a_{ij}| + \sqrt{\epsilon_{\mathrm{mach}}}$ and divides by it — the $\sqrt{\epsilon_{\mathrm{mach}}}$ term is what makes the rescaled entries strictly less than 1, and $\alpha$ so defined is the number the function hands back — while the MATLAB routine uses $\alpha = \max_{ij}|a_{ij}|$ exactly; otherwise both set $\alpha = 1$ and leave $A$ alone. On shape they diverge outright: the Python routine zero-pads a non-square or non-power-of-two input out to the next power of two, the MATLAB routine asserts `N == 2^n` and `N == size(A,2)` and fails on anything else. For real data both form the angle vector $2\arccos(a_{ij})$. For complex data both build a magnitude and a phase vector and compose two uniformly controlled rotations, $R_y$ for the magnitude and $R_z$ for the phase, but not from the same quantities: Python takes the magnitude angles as $2\arccos|a_{ij}|$ and the phase angles as $-2\arg(a_{ij})$ throughout, while MATLAB keeps the signed real part as the magnitude and sets the phase to zero for any entry whose imaginary part falls below $10^{-14}$. Each then solves the uniformly-controlled-rotation angle system $\hat{H}^{\otimes 2n} P_G \hat{\theta} = \theta$ with a scaled fast Walsh–Hadamard transform followed by the Gray-code permutation — the paper states this is "efficiently solved by a classical algorithm in $O(N^2 \log N^2)$ using a fast Walsh–Hadamard transform [20] which is implemented in the reference implementation of FABLE". Compression is where the two part again. Both carry the cutoff criterion, and both implement it as a test on magnitude — coefficients with $|\hat{\theta}_i|$ at or below the cutoff are set to zero, where the paper's step (i) is written "Remove all rotation gates for angles $\hat{\theta}_i \leq \delta_c$" — and in Python that step runs only when a truthy `epsilon` is passed, so the default call compresses nothing. MATLAB carries a second criterion the Python routine does not have, `'percentage'`: sort the $N^2$ coefficients by magnitude and zero the smallest $p\%$ of them. The surviving gates are emitted by the compressed-uniform-rotation builder, which applies the second half of the paper's rule: parity-check each run of consecutive CNOTs and keep one CNOT per control qubit that occurs an odd number of times. The assembled circuit is Hadamards on the row register, the oracle, a SWAP pairing the two $n$-qubit registers, and Hadamards again; the Python routine then closes with `circ.reverse_bits()`, the little-endian convention that puts the leading block in the order its own test asserts.

  Data

  The benchmark inputs are stated in the paper. The accuracy sweep of Fig. 4 uses "randomly generated matrices with entries drawn from the standard normal distribution, the matrices are 2 to 7 qubit operators such that the FABLE circuits require 5 to 15 qubits". The Hamiltonian benchmarks are a Heisenberg XXX chain with $J_x = J_y = J_z$ and $h_z = 0$ on 2 to 7 qubits; Fermi–Hubbard Hamiltonians "generated through OpenFermion [22] for the case $t = 1$, $U = 0$, i.e., non-interacting fermions on a 1D and 2D lattices with two spin orbitals per site and non-periodic boundary conditions", then "mapped to qubits using the Bravyi–Kitaev transformation"; and finite-difference discretised Laplacians — 1D on 2 to 7 qubits with periodic and non-periodic Dirichlet boundary conditions, and 2D five-point-stencil Kronecker sums on rectangular grids "requiring at most 7 qubits".

  Code

  Both implementations are released together at https://github.com/QuantumComputingLab/fable; the paper's introduction states that "Implementations of the FABLE algorithm built on top of QCLAB [12, 13] and Qiskit [14] are made publicly available" there, and the README describes them as "two reference implementations of the FABLE algorithm". The Python version is `py/fable`, depends on `qiskit>=0.19.1` and `numpy>=1.20.3`, and is distributed on PyPI as `fable-circuits` (version 1.0.2, uploaded 2025-01-03); its entry point is `fable(a, epsilon=None)`, returning a Qiskit `QuantumCircuit` and the subnormalization factor `alpha`. The MATLAB version is the single file `fable-qclab/fable.m`, built on QCLAB (cited as version 0.1.3) and called as `[circuit, OA, alpha, info] = fable( A, 'cutoff', 1e-4, logging )`; it also accepts `'percentage', 80`, and the README notes that "The `'percentage'` and `logging` options are only available in the MATLAB version of FABLE." The licence is a three-clause BSD notice with LBNL's standard enhancements clause: "FABLE Copyright (c) 2022, The Regents of the University of California, through Lawrence Berkeley National Laboratory". The whole repository tree is nine files — the two synthesis routines and the Python helper module, a two-case Python test file, `pyproject.toml`, `setup.cfg`, the README and the licence — so the drivers that produced Figs. 4–6 and Table I are not in it.

  Results

  Everything reported in the arXiv version read here (v2, dated 29 July 2022) is classical simulation and gate counting; the paper reports no run on quantum hardware. Fig. 4 shows "the result of noise-free QCLAB [12] simulations for compressed FABLE circuits" — the one benchmark the paper attributes to a named implementation, and it attributes it to the MATLAB/QCLAB half — for randomly generated real-valued 2-qubit ($4 \times 4$) to 7-qubit ($128 \times 128$) matrices at three compression cutoffs $\delta_c = 10^{-6}, 10^{-4}, 10^{-2}$: "the bound from Theorem 2 always holds but is overly pessimistic", and "Not shown in Fig. 4 are the majority of random realizations with an error close to the $10^{-16}$ machine precision." The remaining benchmarks name no implementation. For the Heisenberg XXX chain on 2 to 7 qubits with the compression threshold set to machine precision, "FABLE can give an accurate encoding of a Heisenberg XXX model with just 50% of the CNOT gates and 25% of the Ry rotations". Table I gives the Hubbard counts against the maximum $4^n$: 1D, 2 sites, 4 qubits — 130 CNOT (50.8%) and 65 Ry (25.4%); 1D, 6 sites, 12 qubits — 180,234 CNOT (1.1%) and 81,921 Ry (0.5%); 2D, 2×3 sites, 12 qubits — 252,626 CNOT (1.5%) and 90,113 Ry (0.5%). Turning the interaction on removes the gain: "We ran the same experiment for interacting Hubbard Hamiltonians ($U \neq 0$) but for this case, the matrices do not compress well and the maximum number of gates is required to encode the Hamiltonians." For the Laplacians, "the 2D Laplacians can be compressed better than 1D Laplacians and require in some cases fewer than 5% of the maximum number of gates". The acknowledgment states the work "used resources of the National Energy Research Scientific Computing Center (NERSC)". Separately, the shipped Python tests record what the released code recovers: `test_fable_real` and `test_fable_complex` build a random real and a random complex $8 \times 8$ matrix respectively, run the circuit on Qiskit Aer's `AerSimulator(method="unitary")`, and assert that `a/alpha/2**n` equals the leading $2^n \times 2^n$ block of the returned unitary — its real part, in the real-valued test. The README's usage block states the same relation the other way round, `np.linalg.norm(alpha * N * unitary.data[0:N, 0:N] - A)/np.linalg.norm(A)`. So in the released code the target matrix is recovered from that block by multiplying it by $\alpha\,2^n$, with $\alpha$ the rescaling factor the routine returns.

• PennyLane qml.FABLE template

  About

  Xanadu's PennyLane shipped FABLE as a built-in template. Its v0.36.0 release note says the algorithm "as outlined in arXiv:2205.00081 is now accessible via the `qml.FABLE` template", and places it against PennyLane's existing exact block-encoder: "The usage of `qml.FABLE` is similar to `qml.BlockEncode` but provides a more efficient circuit construction at the cost of a user-defined approximation level, `tol`."

  Methods

  The template takes a matrix and $2n + 1$ wires, where $2^n \times 2^n$ is the size after padding: like the LBNL Python routine it zero-pads an arbitrary $N \times M$ input, warning each time it does so, and its docstring says as much — "FABLE can be implemented for matrices of arbitrary shape and size. When given a $(N \times M)$ matrix, the matrix is padded with zeroes until it is of $(N \times N)$ dimension, where $N$ is equal to $2^n$, and $n$ is an integer." It reuses the Möttönen state-preparation machinery already in PennyLane — `compute_theta` and `gray_code` from `pennylane/templates/state_preparations/mottonen.py` — to turn $\arccos(a_{ij})$ into the uniformly-controlled-rotation angles, emits an `RY(2*theta)` only when `abs(2 * theta) > tol`, and cancels CNOTs by the same parity rule as the paper: a control wire already pending in the accumulator is deleted rather than emitted twice. The circuit ends with SWAPs pairing the two index registers and a second layer of Hadamards. Where it departs from the LBNL reference implementation is on the input, and only for concrete — non-traced — input: it is real-only, a complex input raising `ValueError("Support for imaginary values has not been implemented.")`, and it does not rescale, raising instead when $\max|a_{ij}| > 1$, with the docstring telling the caller "It is also assumed that the values of the input matrix are within $[-1, 1]$. Apply a subnormalization factor if needed." Neither check runs under JIT tracing: both sit inside `if not qml.math.is_abstract(input_matrix):`, the constructor refusing `tol != 0` in that branch instead, and the `tol` threshold is likewise skipped for abstract angles, which are emitted as `RY` unconditionally.

  Data

  None found yet.

  Code

  `pennylane/templates/subroutines/fable.py` in https://github.com/PennyLaneAI/pennylane — Python, Apache License 2.0, "Copyright 2018-2024 Xanadu Quantum Technologies Inc." The documented signature, at the v0.36.0 tag and on the current stable API page alike, is `FABLE(input_matrix, wires, tol=0, id=None)`. It was added by pull request 5107 ("Implement FABLE as a Template (issue 4848)"), merged on 2024-04-11, and first shipped in PennyLane v0.36.0.

  Results

  The v0.36.0 release note's worked example block-encodes $A = \begin{pmatrix} 0.1 & 0.2 \\ 0.3 & 0.4 \end{pmatrix}$ on three wires of the `default.qubit` state-vector simulator with `tol = 0.001`, and recovers it from the circuit matrix as `2 * mat[0:2, 0:2]`, printing `array([[0.1+0.j, 0.2+0.j], [0.3+0.j, 0.4+0.j]])`. The current API documentation runs the same check on `default.qubit` for $\begin{pmatrix} 0.1 & 0.2 \\ 0.3 & -0.2 \end{pmatrix}$ with `tol=0`, recovering the matrix as `2**s * qp.matrix(example_circuit)().real[0 : 2**s, 0 : 2**s]`.