Quantum Signal Processing Phase Factorization via Logarithmic Derivatives

AI summary

Overview

Research area: Quantum algorithms, specifically quantum signal processing (QSP) and the classical preprocessing step that computes QSP phase factors from a target polynomial.

Technical level: Advanced. The paper assumes familiarity with block encodings, Laurent/Chebyshev polynomial lifts, contour integration and residue calculus, and the Newton–Girard identities.

Scope: A single-technique paper that replaces Prony's method in an existing stable QSP phase-factor algorithm with a logarithmic-derivative construction, giving an O(d log²d) method with numerical errors around 10⁻¹² on four standard application benchmarks.

What This Paper Is About

Quantum signal processing implements a polynomial transformation of a matrix by interleaving fixed rotations with queries to a block encoding, and the circuit is specified by a list of angles called phase factors. Computing those phase factors reliably for high-degree polynomials is the bottleneck. This paper modifies a known stable factorization method (attributed to Ying, reference [19]) so that one of its steps — recovering a polynomial from weighted power sums of its roots via Prony's method — is replaced by a logarithmic derivative whose Fourier coefficients give unweighted power sums, after which the Newton–Girard identities recover the polynomial directly.

Key Contributions

1. Elimination of Prony's method. The paper observes that the logarithmic derivative l(z) = h′(z)/h(z) of the associated polynomial h(z) = z^{2d}(1 − a(z)²) has residue 1 at every simple root, so its Fourier coefficients yield the unweighted power sums of the roots inside the unit disk rather than weighted ones. This removes the need to determine unknown residues.

2. Direct polynomial recovery via Newton–Girard. Because the power sums S₁, …, S_{2d} are obtained directly, the coefficients of the monic polynomial m(z) = ∏_{|ξ_j|<1}(z − ξ_j) follow from the Newton–Girard recursion without ever computing roots.

3. Accelerated coefficient recovery. The O(d²) Newton–Girard recursion is replaced by a fast power-series computation: m̃(x) = exp(−Σ_{k≥1} (S_k/k) x^k) evaluated modulo x^{2d+1} by Newton iteration with FFT-based multiplication, costing O(d log d). Combined with the divide-and-conquer phase stripping of reference [15], the overall algorithm is O(d log²d).

4. Simplification by setting b = 0. Because no Prony matrix is formed, the auxiliary polynomial b(t) needed in [19] to make that matrix rank 2d is unnecessary; the authors set b(t) = 0, reducing the problem to constructing a single complementary polynomial. The paper notes the resulting factorization coincides with that of Berntson and Sünderhauf [3], differing only in how it is computed.

Main Findings

• Accuracy: Across Hamiltonian simulation, eigenstate filtering, matrix inversion, and the Fermi–Dirac operator, phase factors are computed with maximum errors on the order of 10⁻¹². Specifically, errors remain below 2 × 10⁻¹² for Hamiltonian simulation, below 5 × 10⁻¹³ for eigenstate filtering, below 3 × 10⁻¹² for matrix inversion even at degree 51633, and below 2 × 10⁻¹² for the Fermi–Dirac examples.

• Runtime: All tested instances complete within a few seconds on a 2.00 GHz AMD Ryzen AI 7 350 laptop in double precision, with runtime growing nearly linearly in d, consistent with the O(d log²d) prediction. The method remains numerically stable at polynomial degrees exceeding 10⁴.

• Comparison table (Table 1): Maximum phase factor error and maximum runtime for the Prony-based algorithm [19] versus the Newton–Girard algorithm: Hamiltonian simulation 5 × 10⁻¹² and 20 s versus 1.4 × 10⁻¹² and 0.3 s; eigenstate filtering 2 × 10⁻¹² and 30 s versus 4.5 × 10⁻¹³ and 0.3 s; matrix inversion 2 × 10⁻¹¹ and 300 s versus 2.4 × 10⁻¹² and 3.5 s; Fermi–Dirac 4.5 × 10⁻¹² and 70 s versus 1.1 × 10⁻¹² and 0.6 s. The paper states the Prony values are taken from [19] rather than from a reimplementation, so they were obtained on different hardware.

• Source of the speedup: The authors state that the decrease in runtime may be largely due to the shift to the divide-and-conquer phase stripping algorithm of [15], rather than to the replacement of Prony's method.

• Degree scaling, Hamiltonian simulation: Experiments use τ = 1000, 2000, 3000, 4000, 5000. Degree grows linearly in τ with slope about 1, consistent with the bound d = τ + O(τ^{1/3} log(1/ε)^{2/3}) [14].

• Degree scaling, eigenstate filtering: Experiments use 1/Δ = 12.5, 25, 50, 100, 200, with k = 20/Δ so d = 2k = 40/Δ. Degree growth is exactly linear because the target is itself a polynomial.

• Degree scaling, matrix inversion: Experiments use κ = 16, 64, 256, 1024; κ = 1024 yields over 50000 phase factors, matching the largest instance in [19]. Predicted d ≈ 10κ√(ln(1/ε)) ≈ 53κ, consistent with the observed d/κ ≈ 50.

• Degree scaling, Fermi–Dirac: Experiments use β = 100, 200, 400, 800, 1600. Predicted d ≈ (β/π) log(1/ε), with predicted slope ln(10¹²)/π ≈ 8.8 versus an observed d/β ≈ 7.

• Excluded costs: In all cases the polynomial approximation step takes under 3 seconds and the function approximation error is below 4 × 10⁻¹¹; following [19], both are excluded from reported times and errors. The function approximation error is stated to be below 40ε for threshold ε = 10⁻¹².

• Error measurement caveat: The phase factor error is the relative L∞ error between the target polynomial and the polynomial reconstructed from the computed phase factors, evaluated at 20 random points on the sampling grid plus the maximizer of |a|; with only 21 points, the measured maximum is described as a noisy estimate of the true maximum.

Methodology in Plain English

The target function is first turned into a polynomial approximation by sampling it, applying an FFT, and discarding frequencies above a degree d chosen so that the discarded coefficients are below about 10⁻¹², with d adjusted to match the parity of the function. The polynomial is scaled so its maximum absolute value is 0.3, which keeps the construction well away from singular behavior.

From that polynomial a single auxiliary polynomial h is built whose roots come in reciprocal pairs: exactly 2d inside the unit disk and 2d outside, with none on the circle. Instead of trying to find those roots, the method samples the rational function h′/h on a grid of about 40d points around the unit circle. A key algebraic fact is that this function splits into a sum of simple terms 1/(z − ξ_j), each with coefficient exactly one. Integrating against the unit circle picks out only the roots inside the disk and returns their plain power sums, with no unknown weights.

Those power sums are then converted into polynomial coefficients. The classical Newton–Girard identities do this in O(d²); the paper instead writes the reversed polynomial as an exponential of a power series built from the power sums and computes that exponential efficiently with FFT-based polynomial arithmetic. The recovered polynomial gives the complementary polynomial, from which p and r are assembled, and the phase factors are peeled off one layer at a time. That peeling is done with the divide-and-conquer strategy of [15], which exploits the fact that after k steps only the top 2k+1 coefficients matter.

The implementation is validated against the QSPPACK reference implementation, and the code is released on GitHub.

Why This Matters

Phase-factor computation is the classical preprocessing that gates most practical QSP-based quantum algorithms; if it is slow or numerically fragile, high-degree transformations become impractical regardless of the quantum circuit cost. This work keeps the numerical stability of the factorization approach while removing a complicated step and reducing the asymptotic cost to O(d log²d), and it makes the connection to the FFT-based algorithm of Berntson and Sünderhauf explicit — the two produce the same canonical complementary polynomial but compute it differently.

Real-world applications demonstrated or targeted in the paper:

• Hamiltonian simulation: approximating cos(τx) and sin(τx) for simulating quantum systems with ||H||₂ ≤ 1; the paper tests τ up to 5000.
• Eigenstate filtering: resolving states near a known energy using Chebyshev filtering with gap Δ; tested down to 1/Δ = 200.
• Matrix inversion: approximating 1/x on a spectrum separated from zero by condition number κ; tested up to κ = 1024 and over 50000 phase factors.
• Fermi–Dirac operator: preparing thermal/fermionic functions with inverse temperature parameter β; tested up to β = 1600.

Industry relevance: libraries and compilers for fault-tolerant quantum algorithms need robust, fast phase-factor generation for high-degree polynomials; a method that handles degrees above 5 × 10⁴ in a few seconds on a laptop, in ordinary double precision, lowers the engineering barrier for building such toolchains.

Future Directions

• Rigorous error analysis of the Newton–Girard step. Recovering coefficients from power sums is described as notoriously ill-conditioned, since errors can be amplified when root magnitudes vary widely. The experiments show errors below 3 × 10⁻¹² above degree 5 × 10⁴, suggesting the ||a||∞ ≤ 0.3 scaling keeps roots far from the unit circle, but a forward error bound is called for to explain this and to identify when it fails.

• Relaxing the ||a||∞ ≤ 0.3 normalization. A larger norm improves the success probability of the resulting circuit and reduces amplitude amplification rounds, but pushes roots toward the unit circle, making the logarithmic derivative nearly singular and its Fourier coefficients slow to decay. Proposed remedies are adaptive sample counts instead of a fixed 40d, integration over deformed contours |z| = r between the inner and outer root sets, and extended or mixed precision only in the most sensitive stages. Characterizing achievable accuracy as a function of 1 − ||a||∞ is named as an open goal.

• Generalizations of QSP. Because the method sets b = 0 and needs only power sums of the inner roots, it applies to any setting requiring a complementary polynomial satisfying |P|² + |Q|² = 1 on the unit circle. Natural next steps are complex-valued targets and generalized QSP with SU(2) rotations, which would lift the parity constraint and could remove the LCU step used to combine even and odd parts.

• Broader use of the cepstral/logarithmic-derivative idea. The author expects the Kolmogorov (cepstral) approach to be useful wherever spectral factorization of a positive trigonometric polynomial is required.

Target Audience

Researchers and practitioners in quantum algorithms who build or use QSP-based methods and need to generate phase factors at high degree — particularly those working on Hamiltonian simulation, eigenstate filtering, matrix inversion, and thermal-state or Fermi–Dirac constructions. It is also relevant to numerical analysts interested in spectral factorization of positive trigonometric polynomials, and to quantum software engineers implementing phase-factor toolchains who want a simple, stable alternative to Prony-based or root-finding-based approaches. Readers without background in block encodings, Laurent polynomial lifts, and residue calculus will find the paper demanding.

Authors’ abstract

Quantum signal processing (QSP) requires the efficient computation of phase factors associated with a polynomial transformation. One existing stable factorization method avoids direct root finding but relies on Prony's method to recover a polynomial from weighted power sums of its roots (arXiv:2202.02671). We present a modification of this approach that eliminates Prony's method by using the logarithmic derivative of an associated polynomial. Because the logarithmic derivative has the same residue at every root after accounting for multiplicity, its Fourier coefficients directly yield the unweighted power sums of the roots inside the unit disk. The Newton--Girard identities then recover the coefficients of the corresponding polynomial without computing its roots. The resulting algorithm has $O(d\log^2 d)$ computational complexity and simplifies the factorization procedure of the original method. In numerical experiments on Hamiltonian simulation, eigenstate filtering, matrix inversion, and the Fermi--Dirac operator, our implementation computes phase factors in seconds with maximum errors on the order of $10^{-12}$.