MethodLayer 2

Direct method: root finding, then layer stripping

Compute the complementary polynomial by finding the roots of a high-degree polynomial, then strip off one phase factor at a time from the assembled $SU(2)$ -valued product. Haah's product decomposition is the version of this route that comes with a full arithmetic-model analysis.

Takes

Chebyshev coefficients of a real polynomial $P$ of degree $d$ with definite parity and $|P(x)| \le 1$ on $[-1,1]$ , plus a target accuracy $\varepsilon$ .

Returns

A phase sequence $\Phi \in R^{d+1}$ , often symmetric ( $\phi_j = \phi_{d-j}$ ), together with the classical running time and the arithmetic precision the method requires.

Same contract as the slot it fills.

From Polynomial approximation to QSP phase sequence

A circle is an object you are holding. This method is drawn heavier, opened into its own steps; the other lines between the same two ends are the alternatives recorded for the same slot. Circles are named on hover, and each one is a link.

Nothing drawn here has a recorded way through it that this figure leaves shut. See it on the map

What it fills

QSP phase factors
Given an admissible polynomial, compute the phase sequence $Φ$ that makes the quantum-signal-processing product reproduce it to accuracy $\varepsilon$ in classical finite-precision arithmetic.

When it applies

Root finding is the numerical weak point of the family. Ni and Ying's survey groups this family — the Gilyén-Su-Low-Wiebe construction, Haah, and the halving/capitalization method — as requiring $O(d\,\mathrm{polylog}(d/\varepsilon))$ bits of precision, citing Haah, which means variable-precision arithmetic rather than double precision. They also state that the stability of the layer-stripping process used by most direct methods remains an open question.

Requires

Every step this method names moves its route along, so there is nothing it needs alongside them.

Example

given  a degree-N approximation at accuracy epsilon
       # N is cost's symbol (Haah); the precision claim below is the
       # survey's and is written in d -- the record uses both

# stage 1 -- the complementary polynomial
find the roots of a high-degree polynomial
build the complementary polynomial from them
# root finding is the numerical weak point of the family

# stage 2 -- layer stripping
assemble the SU(2)-valued product
repeat, one phase factor at a time:
    strip a phase factor off the assembled product
# how a factor is stripped, and what ends the stripping, are not stated

return the phase factors

# precision -- this is what separates the route from the double-precision
#   front ends: Ni and Ying's survey groups this family (the
#   Gilyen-Su-Low-Wiebe construction, Haah, and the
#   halving/capitalization method) as requiring O(d polylog(d/epsilon))
#   bits of precision, citing Haah, which means variable-precision
#   arithmetic rather than double precision

# cost -- Haah, whose product decomposition is the version of this route
#   that comes with a full arithmetic-model analysis:
#   O(N^3 polylog(N/epsilon)) time for a degree-N approximation at
#   accuracy epsilon, under the random-access memory model of
#   computation. The model is the point of the result -- earlier
#   efficiency claims had assumed a strong arithmetic model and lacked
#   numerical stability analysis

# open -- the stability of the layer-stripping process used by most
#   direct methods remains an open question

# unreconciled -- the same survey places the halving/capitalization
#   method of Chao, Ding, Gilyen, Huang and Szegedy in this family, while
#   that paper itself reports finding sequences of more than 3000 angles
#   within 5 minutes in standard double precision arithmetic. The two
#   claims have not been reconciled in the literature

Cost, as the source states it

Haah: $O(N³ \mathrm{polylog}(N/\varepsilon))$ time for a degree- $N$ approximation at accuracy $\varepsilon$ , under the random-access memory model of computation. The model is the point of the result — earlier efficiency claims had assumed a strong arithmetic model of computation and lacked numerical stability analysis, and this work replaces that with a realistic one.

Implementations

pyqsp's laurent method: completion by root finding, then divide-and-conquer decomposition
- Finding Angles for Quantum Signal Processing with Machine Precision
pyqsp (github.com/ichuang/pyqsp) is a Python package for quantum signal processing whose QuantumSignalProcessingPhases function offers three phase-finding routes chosen with a method argument: laurent (the default), tf, and sym_qsp. The README states that the laurent method "employs techniques originated in" Chao, Ding, Gilyen, Huang and Szegedy's paper and "extends code from its attached repository" at github.com/alibaba-edu/angle-sequence, and that it "exactly computes phases by studying the properties of the desired polynomials using a divide-and-conquer approach." The README separately names this same decomposition code "the halving algorithm," Chao et al.'s own term for what their abstract calls "a novel component we call halving"; that paper's own experiments distinguish halving by name from Jeongwan Haah's decomposition algorithm, which it calls "carving," a different, sequential one-layer-at-a-time procedure the paper says needs far higher-precision arithmetic than halving's standard 64-bit machine precision throughout. This entry therefore lists only Chao et al.'s paper as a source for pyqsp's laurent route: pyqsp's own CITATION file bundles Haah's paper into one undifferentiated, package-wide reference list alongside papers backing pyqsp's other routes, rather than attributing it to this decomposition code specifically. This is the root-finding route this record's family covers, distinct from pyqsp's own Newton-based sym_qsp route and its TensorFlow-based tf optimization route, each of which is a different method in this capability.
Calling QuantumSignalProcessingPhases with method="laurent" (or calling the lower-level angle_sequence function directly) runs completion by root finding in pyqsp/completion.py's completion_from_root_finding, which for the default Wx/Wz convention dispatches to _fg_completion: it forms the Laurent polynomial $1 - F(w)F(w^{-1})$ , finds its roots with numpy's Polynomial.roots, and keeps one root from each reciprocal pair (those with $|root| < 1$ and non-negative imaginary part) per an inline comment in the same function, not its formal docstring, which states the input is "a real, self-inverse Laurent polynomial with no root on the unit circle" whose roots "come in reciprocal, conjugate pairs." The kept roots are assembled into the completing polynomial's coefficients via an FFT-based product rather than direct multiplication of the linear factors, for numerical stability. Before completion, angle_sequence.py applies a "Capitalization" step that puts a slice of the error budget onto the top Laurent coefficient. The completed Low-algebra element is then turned into phase factors by pyqsp/decomposition.py's angseq, which is not a single-layer-at-a-time peel: it is a recursive divide-and-conquer routine whose decompose helper splits the element's degree in half at each level and solves a linear least-squares system (numpy.linalg.lstsq) for the left factor before recursing on both halves.
Repository github.com/ichuang/pyqsp (Python; default branch master). Entry point pyqsp.angle_sequence.QuantumSignalProcessingPhases(poly, ..., method="laurent"), with "laurent" the default value of method; also the lower-level pyqsp.angle_sequence.angle_sequence(p, eps, suc), which runs the same completion-then-decomposition pipeline directly on Laurent-polynomial coefficients. Completion: pyqsp.completion.completion_from_root_finding and its helpers _fg_completion and _pq_completion. Decomposition: pyqsp.decomposition.angseq and decompose. The Laurent-polynomial and Low-algebra objects these operate on are pyqsp.LPoly.LPoly and pyqsp.LPoly.LAlg.
No degree, error, or timing benchmark is reported for the laurent route specifically, unlike the package's sym_qsp route. The README instead characterizes it qualitatively: by default it is "typically quite fast, but can become unstable for high-degree polynomials due to roundoff errors, requiring some randomization." A 2024-10-24 changelog entry in the same README reports that moving the codebase to work internally in the Chebyshev rather than monomial basis "improves the numerical stability of the laurent method substantially." Elsewhere the README attributes an earlier failure approximating $1/x$ partly to "the instability of direct polynomial completion methods in the laurent approach."

Where the claim is contested

There is a live tension over the precision requirement. The survey above places the halving/capitalization method of Chao, Ding, Gilyén, Huang and Szegedy in the root-finding family said to need $O(d \mathrm{polylog}(d/\varepsilon))$ bits, while that paper itself reports finding sequences of more than 3000 angles within 5 minutes in standard double precision arithmetic. The two claims have not been reconciled in the literature.

What it needs

Nobody has taken this apart yet. That is a gap in this graph, not a claim that the method has no parts.

Other ways to fill the same slot

Different approaches

Least-squares optimization of the phase factors
Instead of constructing the complementary polynomial, minimise the mean squared difference between the QSP response $Re[\langle0|U_Φ(x_j)|0\rangle]$ and the target $f$ , evaluated at the positive roots of the Chebyshev polynomial $T_{2d̃}$ , over symmetric phase sequences. Gradients come from $SU(2)$ matrix products, so root finding is avoided entirely and the method runs in standard double precision — but it degrades as $‖f‖_\infty$ approaches 1.
Newton's method for symmetric QSP
Treat phase-factor finding as a nonlinear system rather than a minimization, and solve it with a Newton iteration built for symmetric QSP. The matrix-product-state structure of symmetric QSP makes computing the Jacobian cost about the same as a single function evaluation.
Stable factorization via Prony's method
Build the complementary polynomial directly, with Prony's method as the key step, then obtain the phase factors by factorization. This avoids root finding of high-degree polynomials, which is the step that forces variable-precision arithmetic elsewhere in the direct family.

In the Atlas

No record in the Atlas covers this yet. The catalogue is circuits and primitives; this part of the literature is not in it.

Sources

Product Decomposition of Periodic Functions in Quantum Signal Processing ↗Jeongwan Haah · 2018theory yes · simulation yes · hardware no
Fast Phase Factor Finding for Quantum Signal Processing ↗Hongkang Ni, Lexing Ying · 2024theory yes · simulation yes · hardware no
Finding Angles for Quantum Signal Processing with Machine Precision ↗Rui Chao, Dawei Ding, Andras Gilyen, Cupjin Huang, Mario Szegedy · 2020theory yes · simulation yes · hardware no

This one, drawn