Pith. sign in

REVIEW 3 major objections 5 minor 25 references

This paper proposes a Fourier-extension LCU decomposition that approximates any non-unitary operator by 4m unitaries with error decaying exponentially in m, while keeping block-encoding subnormalization O(||A||_2 log log 1/ε).

Reviewed by Pith at T0; open to challenge. T0 means a machine referee read the full paper against a public rubric. the ladder, T0–T4 →

T0 review · deepseek-v4-flash

2026-08-03 08:06 UTC pith:56LQK6JN

load-bearing objection Fourier-LCU is a genuine new decomposition with strong small-m numerics, but the flagship alpha=O(loglog 1/epsilon) scaling is an openly empirical extrapolation from m<=16; treat it as plausible, not proven. the 3 major comments →

arxiv 2601.18024 v2 pith:56LQK6JN submitted 2026-01-25 quant-ph

Fourier extensions for matrix-function block encodings with error-independent subnormalization bounds

classification quant-ph PACS 03.67.Ac
keywords linear combination of unitariesblock encodingFourier extensionsubnormalizationquantum numerical linear algebramatrix functionsopen quantum systemsquantum singular value transformation
verification ladder T0 review T1 audit T2 compute T3 formal T4 reserved

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The paper's central thesis is that the long-standing trade-off between approximation error and subnormalization in general LCU block encodings can be broken by replacing the local Taylor/sine approximation of the map τ with a Fourier extension—a sine series on a larger periodic domain that approximates τ on a smaller interior interval. For any non-unitary A=H1+iH2, this yields A as a linear combination of 4m unitaries whose error decays exponentially in m, and whose subnormalization is α=(2η/π)max(||H1||_2,||H2||_2)Σ|a_k|. Observing numerically that Σ|a_k|=O(log m) for the least-squares coefficients, the paper concludes α=O(||A||_2 log log 1/ε), near-optimal since α must always grow with ||A||_2. It also shows that by regularizing the least-squares problem with an L1 penalty, one can traverse a Pareto front that lowers α further at a fixed error budget. A sympathetic reader should care because this removes the polynomial penalty that has limited practical use of non-unitary block encodings in quantum simulation and QSVT workflows.

Core claim

For A=H1+iH2, the Fourier-LCU decomposition is A ≈ (1/(2τ)) Σ_{k=1}^m a_k [ i e^{-ikτH1} − e^{-ikτH2} − i e^{ikτH1} + e^{ikτH2} ], with τ=π/(η max[||H1||_2,||H2||_2]). The coefficients a_k are the least-squares Fourier-extension coefficients for f(τ)=τ on [−π/η,π/η], and each sine is written as a difference of two unitaries, so the total is 4m unitaries. The approximation error converges exponentially in m, while the subnormalization is α=(2η/π)max[||H1||_2,||H2||_2]Σ|a_k|. On the evidence of m≤16 computations, Σ|a_k|=O(log m), giving α=O(||A||_2 log log 1/ε) when m=O(log 1/ε). The paper further introduces a regularized fitting objective whose L1 term directly penalizes α, demonstrating nume

What carries the argument

The central object is the Fourier-extension sine series f(τ)≈Σ_{k=1}^m a_k sin(kτ), solved by continuous least squares on the interior interval [−π/η,π/η] of a [−π,π] periodic domain. Extension factor η>1 creates a buffer zone that makes the periodic extension smooth, so coefficients decay exponentially; η=2+0.460m^{-0.319} minimizes the cost mα. The sine terms are rewritten via e^{ikτH}−e^{-ikτH}=2i sin(kτH), converting the series into an LCU. The subnormalization is controlled by S_m=Σ|a_k|, and the overcomplete dictionary (η>1) admits a null space that regularized fitting exploits to reduce S_m at fixed error.

Load-bearing premise

The double-logarithmic subnormalization claim rests entirely on the empirical observation that the L1 norm of the least-squares Fourier-extension coefficients satisfies Σ|a_k|=O(log m) for m≤16, and on the assumption that this growth persists for all m needed to reach the target error; the paper says it cannot conclude α becomes independent of m.

What would settle it

Evaluate the least-squares Fourier-extension coefficients in high-precision arithmetic for m beyond 16 and test whether Σ|a_k| continues to grow like log m; if the sum instead grows polynomially in m, the α=O(log log 1/ε) claim fails even if the exponential error convergence holds. A complementary check is to solve for the minimal achievable α at a fixed ε and compare it against the claimed O(log log 1/ε) bound.

Watch this falsifier — get emailed when new claim-graph text bears on it.

If this is right

  • Replacing finite-difference decompositions with Fourier LCU reduces the postselection penalty from α=poly(1/ε) to α=O(log log 1/ε) for a given accuracy, so amplitude amplification becomes far less demanding.
  • Algorithms for open quantum systems, linear systems, and differential equations that currently use four-unitary decompositions can inherit this improved scaling without structural assumptions on the operator.
  • At a fixed error budget, L1-regularized coefficient selection yields smaller α (about 2 versus 5 in the paper's numerical example), meaning higher signal strength and fewer required postselection attempts.
  • For s-sparse Hamiltonians with oracular access, the gate complexity of preparing one block encoding is O(Q[s+1]||A||_2 log^2(1/ε) log^2 log(1/ε)), where Q encodes input normalization.
  • The same Fourier-extended sine series can encode Hermitian eigenvalue transforms and odd singular-value transforms, making it a direct access point for QSVT and quantum linear systems algorithms.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • Beyond the paper: the α=O(log log 1/ε) guarantee is only as solid as the empirical Σ|a_k|=O(log m) observation; proving this growth from Fourier-extension theory would remove the paper's explicit caveat that it cannot rule out m-dependence.
  • Beyond the paper: because Fourier extensions are known to grow outside the fitted interval, using these block encodings for operators whose spectra approach the interval boundary may require shrinking τ by a constant factor; this is a testable design choice not analyzed in the paper.
  • Beyond the paper: if the regularized-coefficient Pareto fronts continue to flatline at larger m, one could expect near-constant subnormalization in practice even for machine-precision errors, shifting the practical bottleneck to the depth of the controlled Hamiltonian simulations.
  • Beyond the paper: the same coefficient machinery might be reused for other functions of Hermitian operators (for example f(τ)=1/τ or sign functions) by Fourier-extending the appropriate odd function, although the paper only treats the identity map.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit.

Referee Report

3 major / 5 minor

Summary. The paper proposes a linear-combination-of-unitaries (LCU) decomposition for non-unitary operators based on a Fourier extension of the identity map f(τ)=τ on a subinterval of a larger periodic domain. The sine-series coefficients are obtained either by least squares (Section 4.1, Table A) or by a regularized optimization that trades error against subnormalization (Section 4.2, Table B). The authors give a quantum circuit implementation (Section 3) and claim that the resulting block encoding has exponentially decaying error in m and subnormalization α=O(∥A∥_2 loglog 1/ε), a substantial improvement over the polynomial dependence of finite-difference decompositions such as Eq. (5). Numerical demonstrations include coefficient tables and a statevector simulation of a dephasing Lindblad equation (Section 6).

Significance. If the advertised asymptotic scaling were rigorously established, Fourier-LCU would be a valuable primitive for quantum numerical linear algebra, turning generic non-unitary block encodings into a near-optimal resource. The fixed-m decomposition and circuit construction are mathematically sound and clearly presented, and the accompanying coefficient tables and statevector simulations are useful. However, the central complexity claim—the double-logarithmic subnormalization—is not proven; the manuscript itself states that the relevant coefficient norm growth is only an empirical observation for m≤16. Since this scaling is the main advertised advantage over existing methods, the significance of the paper as it stands is substantially reduced.

major comments (3)
  1. [Section 5, final paragraph (after Eq. (14))] The central claim α=O(∥A∥_2 log m)=O(∥A∥_2 loglog 1/ε) is not established. The text explicitly says 'we are unable to conclude that α becomes independent of m' and then relies on the 'empirical observation ∑|a_k|=O(log m) for the range of m≤16 relevant to this study.' Because m=O(log 1/ε) can be arbitrarily large as ε→0, evidence for m≤16 cannot determine the asymptotic scaling. The cited Fourier-extension behavior (Ref. [25]) includes cases with non-decaying or growing coefficients, so the observed exponential decay in Table A is not guaranteed to persist. This is the load-bearing point behind the abstract and conclusion; without a proof or high-precision evidence for substantially larger m, the advertised advantage over Eq. (5) is unsupported.
  2. [Section 4.1, Eq. (18) and Fig. 2d] The prescription η=2+0.460m^{-0.319} is a fitted curve over m≤16, and the associated inferences that η*→2 and S'_∞(2)=0 are not proven. The coefficient tables and the error/subnormalization estimates all depend on this choice of η. Since η is used in Eq. (14) for α and in Eq. (16) for the required m, the empirical nature of this fit propagates into the main resource analysis. A rigorous or at least more extensive numerical characterization of η*(m) is needed before the asymptotic complexity statements can be accepted.
  3. [Section 4.2, Corollary 1.2 and Fig. 3] The claim that regularized coefficients reduce α as m grows rests on numerical Pareto-front computations for m up to 64. Corollary 1.2 only establishes that α*_m(ε) converges for fixed ε; it does not quantify the limit or show that it remains small for the target-accuracy range. The λ-sweep procedure is an optimization heuristic, and the statement that the benefit is 'robust' is stronger than what the numerical evidence supports. This is secondary to the main complexity claim, but it should be tempered.
minor comments (5)
  1. [Table A] The table layout is confusing: rows for different m values are merged into a two-column block, with repeated 'k' headers. Separating the cases (e.g., with clear m labels in each column) would improve readability.
  2. [Eq. (16)] The statement m=O(1/δ(η) log 1/ε) is only an order-of-magnitude estimate. The constant and the range of validity should be specified, especially because the decay rate depends on the fitted η.
  3. [Abstract and Section 5] The notation 'loglog 1/ε' is nonstandard; it should be written as log log(1/ε) or log log(1/ε) with parentheses, to avoid confusion with a logarithm of a product.
  4. [Section 5, second paragraph] The gate count O([sm+log{1/ε}]m log m) is derived under the assumption ∥H_max∥=O(∥A∥_2), which is stated but not justified in detail. A short explanation of this bound would be helpful.
  5. [Section 6, Fig. 4] The comparison between least-squares and regularized coefficients is insightful, but the conclusion that the regularized coefficients are useful for 'practical quantum computation' is based on a single example; a brief discussion of parameter regimes would strengthen the claim.

Circularity Check

1 steps flagged

Central α=O(loglog 1/ε) bound reduces to an empirical Sm=O(log m) fit from the authors' own coefficients; no independent derivation.

specific steps
  1. fitted input called prediction [Section 5, final paragraph before 'Combining the above contributions' and the following paragraph]
    "Therefore, we are unable to conclude that α becomes independent of m. Instead, we consider the empirical observation ∑m k=1|ak|=O(logm) for the range of m≤16 relevant to this study, as supported by Fig. 2. ... Combining the above contributions, α=O(∥A∥2logm)=O(∥A∥2loglog(1/ε))."

    Equation (14) defines α = (2η/π) max[∥H1∥2,∥H2∥2] Sm, so the paper's flagship scaling α=O(loglog 1/ε) is obtained by substituting the claimed Sm=O(log m) into that identity. But Sm=O(log m) is not proved; it is read off from the same least-squares coefficients (Table A, Fig. 2c) used to construct the LCU. The paper explicitly disclaims an actual asymptotic conclusion ('we are unable to conclude...'), then elevates the observed trend for m≤16 to the scaling law that drives the complexity claims. This makes the central asymptotic advantage an extrapolation of the authors' own fitted data, not an independently derived result; if the trend changes at larger m, the claimed α scaling collapses.

full rationale

The LCU decomposition itself is not circular: Eq. (9) is a functional-calculus application of the Fourier sine series to the Hermitian/anti-Hermitian parts of A, and the error and subnormalization identities in Eqs. (13)–(14) are exact algebraic consequences. There is no load-bearing self-citation chain, no imported uniqueness theorem, and no renamed known result. The only circularity-like issue is the asymptotic subnormalization bound: the paper's headline double-logarithmic scaling rests on an empirical observation of Sm for m≤16 rather than on a theorem or an independently validated extrapolation. Because the paper is transparent about this empirical basis, the issue is a partial circularity (fitted input called a prediction) rather than a complete reduction by definition. Score 4 reflects that the main mathematical construction has independent content, but its central performance claim is not independently established.

Axiom & Free-Parameter Ledger

2 free parameters · 4 axioms · 0 invented entities

No new physical entities are introduced. The Fourier-extension dictionary and regularized coefficient selection are algorithmic constructs, not entities requiring independent falsifiable evidence. The load-bearing assumptions are the exponential error decay with buffer thickness, the empirical L1-norm growth, and the Hamiltonian-simulation oracle model.

free parameters (2)
  • Extension factor η(m) = 2 + 0.460 m^{-0.319}
    Fitted to numerical minima of mα over m (Fig. 2d); then used to select coefficients and derive the m=O(log 1/ε) and α scaling.
  • Regularization hyperparameter λ = swept toward 0
    Controls the trade-off between approximation error and subnormalization in Eq. (20); no principled selection rule is provided.
axioms (4)
  • domain assumption Fourier-extension error decays exponentially with buffer thickness: ε=O(e^{-δ(η)m})
    Invoked in Section 4.1 to justify the η choice and m=O(1/δ log 1/ε); the authors note this 'is not always the case' (ref 19).
  • ad hoc to paper Empirical coefficient-norm growth ∑|a_k|=O(logm) holds for all m relevant to target ε
    Section 5 explicitly relies on this for α=O(loglog 1/ε); no proof is given, and the authors state they cannot conclude the behavior for larger m.
  • domain assumption Hermitian parts satisfy max[||H1||2,||H2||2]=O(||A||2) and ∥A∥2=O(1)
    Used in Section 5 to simplify α and the final gate-complexity expression.
  • standard math Hamiltonian simulation of e^{±ikτHj} costs O(s||H||max kτ + log 1/ε) gates with oracular access
    Imported from Low and Chuang (2019); controls the reported circuit-depth estimates.

pith-pipeline@v1.3.0-alltime-deepseek · 14020 in / 13902 out tokens · 131871 ms · 2026-08-03T08:06:47.130642+00:00 · methodology

0 comments
read the original abstract

Block encodings of non-unitary matrix functions are central to quantum numerical linear algebra. Hamiltonian simulation is a natural input model for Hermitian matrices, but accurate block encodings often incur large subnormalization. We decouple the accuracy from the subnormalization by formulating the matrix-function block encoding as a Fourier-extension approximation problem, yielding a linear combination of unitaries for Hermitian matrix inputs. Fourier extensions approximate non-periodic functions by a Fourier series on a larger periodic domain, creating redundant coefficients that can be optimized for their absolute sum, and hence subnormalization. The coefficients may be chosen for optimal subnormalization with algebraic convergence, by tuning the subnormalization bound for increasing rates of exponential convergence, or by Sobolev-regularized fitting to accommodate more general spectral sets. Fourier-extension block encodings apply to eigenvalue transforms of Hermitian matrices, or to odd singular-value transforms of general matrices, including as a quantum linear systems algorithm.

Figures

Figures reproduced from arXiv: 2601.18024 by Ben Adcock, Peter Brearley, Thomas L. Howarth.

Figure 1
Figure 1. Figure 1: Quantum circuit for Fourier LCU, implementing the decomposition in Eq. (9). The unitaries V and W are defined in Eqs. (11) and (12) respectively. Then, the phases are encoded into the first column of a second unitary, W, defined as Wi,0 = κ ∗ i |κi | s |κi | ∑ 4m j=1 |κj | , (12) such that Vi,0W∗ i,0 = κi/α, with α = ∑ 4m j=1 |κj | and the asterisk denoting the complex conjugate. When m is a power of 2, th… view at source ↗
Figure 2
Figure 2. Figure 2: (a) Fourier extension sine series to approximate f(τ) = τ by Eq. (6) on the subinterval τ ∈ [−π/η,π/η], corresponding to η = 2+0.460m −0.319, using the coefficients given in Table A. (b) Error L 2 norm on [−π/η,π/η]. (c) Block-encoding subnormalisation given by Eq. (14). (d) Optimal extension factor η to minimise mα given by Eq. (18). which we evaluate using Gauss-Legendre quadrature. The coefficients for … view at source ↗
Figure 3
Figure 3. Figure 3: Left and centre: The relationship between λ, ε, and α. As λ → 0, the solver trades subnormalisation for higher precision. Right: The Pareto front. The horizontal lines show the cost of standard least-squares fits a LS. The vertical distance to the regularised curves represents the subnormalisation improvement, with the flatlining behaviour demonstrating convergence to the stable limit α ⋆ ∞. While the λ-sw… view at source ↗
Figure 4
Figure 4. Figure 4: Comparison of the coefficient selection strategies in Section 4 using statevector simulations. (a) Error L 2 norm ∥|ρ(t)⟩ −⃗ρ(t)∥2 between the computed statevector |ρ(t)⟩ and the analytical solution ⃗ρ(t) rescaled to norm 1, against the number of terms in the Fourier series m, where the number of unitaries in the LCU is 4m. (b) The subnormalisation α of the resulting block encoding, consistent with the def… view at source ↗

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

25 extracted references · 1 linked inside Pith

  1. [1]

    Childs, A. M. & Wiebe, N. Hamiltonian simulation using linear combinations of unitary operations.Quantum Inf. & Comput.12, 901–924 (2012)

  2. [2]

    & Murua, A

    Blanes, S., Casas, F. & Murua, A. Splitting methods for differential equations.Acta Numer.33, 1–161 (2024)

  3. [3]

    M., Kothari, R

    Childs, A. M., Kothari, R. & Somma, R. D. Quantum algorithm for systems of linear equations with exponentially improved dependence on precision.SIAM J. on Comput.46, 1920–1950 (2017)

  4. [4]

    W., Childs, A

    Berry, D. W., Childs, A. M., Ostrander, A. & Wang, G. Quantum algorithm for linear differential equations with exponentially improved dependence on precision.Commun. Math. Phys.356, 1057–1081 (2017)

  5. [5]

    W., Head-Marsden, K., Sager, L

    Schlimgen, A. W., Head-Marsden, K., Sager, L. M., Narang, P. & Mazziotti, D. A. Quantum simulation of open quantum systems using a unitary decomposition of operators.Phys. Rev. Lett.127, 270503 (2021)

  6. [6]

    Chowdhury, A. N. & Somma, R. D. Quantum algorithms for Gibbs sampling and hitting-time estimation. Quantum Inf. & Comput.17, 41–64 (2017). 11/14

  7. [7]

    & Succi, S

    Sanavio, C. & Succi, S. Lattice boltzmann–carleman quantum algorithm and circuit for fluid flows at moderate reynolds number.AVS Quantum Sci.6(2024). 8.Chakraborty, S., Morolia, A. & Peduri, A. Quantum regularized least squares.Quantum7, 988 (2023)

  8. [9]

    Gilyén, A., Su, Y ., Low, G. H. & Wiebe, N. Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics. InProceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing, 193–204 (2019)

  9. [10]

    M., Rossi, Z

    Martyn, J. M., Rossi, Z. M., Tan, A. K. & Chuang, I. L. Grand unification of quantum algorithms.PRX Quantum 2, 040203 (2021)

  10. [11]

    & Gazda, A

    Koska, O., Baboulin, M. & Gazda, A. A tree-approach Pauli decomposition algorithm with application to quantum computing. InISC High Performance 2024 Research Paper Proceedings (39th International Conference), 1–11 (Prometeus GmbH, 2024)

  11. [12]

    Wan, L.-C.et al.Block-encoding-based quantum algorithm for linear systems with displacement structures. Phys. Rev. A104, 062414 (2021)

  12. [13]

    & Rung, T

    Over, P., Bengoechea, S., Brearley, P., Laizet, S. & Rung, T. Quantum algorithm for the advection-diffusion equation by direct block encoding of the time-marching operator.Phys. Rev. A112, L010401 (2025). 14.Wu, P. Additive combinations of special operators.Banach Cent. Publ.30, 337–361 (1994)

  13. [15]

    Li, X.et al.Toward quantum simulation of non-markovian open quantum dynamics: A universal and compact theory.Phys. Rev. A110, 032620 (2024)

  14. [16]

    Bharadwaj, S. S. & Sreenivasan, K. R. Compact quantum algorithms for time-dependent differential equations. Phys. Rev. Res.7, 023262 (2025)

  15. [17]

    & Tapp, A

    Brassard, G., Hoyer, P., Mosca, M. & Tapp, A. Quantum amplitude amplification and estimation.arXiv preprint quant-ph/0005055(2000)

  16. [18]

    & Tong, Y

    Fang, D., Lin, L. & Tong, Y . Time-marching based quantum solvers for time-dependent linear differential equations.Quantum7, 955 (2023)

  17. [19]

    & Huybrechs, D

    Webb, M., Coppé, V . & Huybrechs, D. Pointwise and uniform convergence of fourier extensions.Constr. approximation52, 139–175 (2020)

  18. [20]

    Grover, L. K. A fast quantum mechanical algorithm for database search. InProceedings of the Twenty-Eighth Annual ACM Symposium on Theory of Computing, 212–219 (1996)

  19. [21]

    S., Donoho, D

    Chen, S. S., Donoho, D. L. & Saunders, M. A. Atomic decomposition by basis pursuit.SIAM review43, 129–159 (2001)

  20. [22]

    Brent, R. P. An algorithm with guaranteed convergence for finding a zero of a function.The computer journal 14, 422–425 (1971). 23.Barenco, A.et al.Elementary gates for quantum computation.Phys. review A52, 3457 (1995). 24.Low, G. H. & Chuang, I. L. Hamiltonian simulation by qubitization.Quantum3, 163 (2019)

  21. [25]

    On the fourier extension of nonperiodic functions.SIAM J

    Huybrechs, D. On the fourier extension of nonperiodic functions.SIAM J. on Numer. Analysis47, 4326–4355 (2010). 12/14

  22. [26]

    W., Childs, A

    Berry, D. W., Childs, A. M., Cleve, R., Kothari, R. & Somma, R. D. Exponential improvement in precision for simulating sparse hamiltonians. InProceedings of the forty-sixth annual ACM symposium on Theory of computing, 283–292 (2014)

  23. [27]

    A., Sanavio, C., Perotto, S

    Zecchi, A. A., Sanavio, C., Perotto, S. & Succi, S. Improved amplitude amplification strategies for the quantum simulation of classical transport problems.Quantum Sci. Technol.10, 035039 (2025)

  24. [28]

    Low, G. H. & Chuang, I. L. Optimal hamiltonian simulation by quantum signal processing.Phys. review letters 118, 010501 (2017)

  25. [29]

    Table A.Fourier coefficients for exponential convergence in Eq

    Hughes, A.et al.Trapped-ion two-qubit gates with> 99.99% fidelity without ground-state cooling.arXiv preprint arXiv:2510.17286(2025). Table A.Fourier coefficients for exponential convergence in Eq. (6) form=1, 2, 4, 8, and 16, using η=2+0.460m −0.319. k a k 1 1.1749265633763890×10 0 1 1.4649300593981140×10 0 2−2.4671045529932240×10 −1 1 1.6867069657827318...