PaperlayerView exact revision on arXiv

Source-linked AI summary

Efficient phase-factor evaluation in quantum signal processing

Yulong Dong, Xiang Meng, K. Birgitta Whaley, Lin Lin

arXiv:2002.11649v2quant-phmath.OCphysics.comp-ph

TL;DR

QSP phase-factor computation is difficult and existing approaches require variable precision arithmetic or become unreliable at high polynomial degree. The paper uses stable optimization over SU(2) matrix products and demonstrates accurate phase factors for high-degree polynomials using double precision.

  • Problem

    Finding the phase factors associated with a given QSP polynomial is extremely difficult, while existing direct methods require variable precision arithmetic and struggle at large degree.

  • Method

    The paper iteratively optimizes a loss function over phase factors, using numerically stable SU(2) matrix multiplications instead of recursively constructing complementing polynomials.

  • Results

    The optimization method accurately determines phase factors for polynomial degrees above 10,000 with function-approximation error below 10^-12 using standard double precision arithmetic.

  • Takeaways & Limitations

    The approach provides an efficient and accurate route to evaluating phase factors needed to build QSP circuits for unitary representations of non-unitary operations.

  • Takeaways & Limitations

    The optimization landscape is highly non-convex, so random initialization can trap the solver in a local minimum and the zero-phase initialization is a nonzero stationary point.

Abstract

from arXiv · show

Quantum signal processing (QSP) is a powerful quantum algorithm to exactly implement matrix polynomials on quantum computers. Asymptotic analysis of quantum algorithms based on QSP has shown that asymptotically optimal results can in principle be obtained for a range of tasks, such as Hamiltonian simulation and the quantum linear system problem. A further benefit of QSP is that it uses a minimal number of ancilla qubits, which facilitates its implementation on near-to-intermediate term quantum architectures. However, there is so far no classically stable algorithm allowing computation of the phase factors that are needed to build QSP circuits. Existing methods require the usage of variable precision arithmetic and can only be applied to polynomials of relatively low degree. We present here an optimization based method that can accurately compute the phase factors using standard double precision arithmetic operations. We demonstrate the performance of this approach with applications to Hamiltonian simulation, eigenvalue filtering, and the quantum linear system problems. Our numerical results show that the optimization algorithm can find phase factors to accurately approximate polynomials of degree larger than $10,000$ with error below $10^{-12}$.

I. INTRODUCTION

QSP offers a general, ancilla-efficient route to representing matrix functions, but computing its phase factors remains numerically difficult. The paper proposes optimization-based phase-factor evaluation that avoids high-precision recursive methods and supports high-degree approximations.

  • Motivation: QSP represents broad classes of matrix functions through eigenvalue or singular-value transformations while using few ancilla qubits.Applications include Hamiltonian simulation, matrix inversion, and eigenvalue filtering.
  • Challenge: Finding phase factors for a desired polynomial is an extremely difficult inverse problem that has limited practical QSP applications.Earlier work established existence, but constructive procedures remained challenging.
  • Limitations of prior methods: Previous methods require high-precision root solving or variable-precision arithmetic, with bit complexity O(d log(d/ϵ)) and observed runtime O(d3).These methods are therefore not classically numerically stable for high-degree polynomials.
  • Proposed approach: The proposed optimization method uses standard double precision and reaches L∞ error ≲10^-12 for polynomial degrees ≳10,000.It minimizes a loss function through numerically stable SU(2) matrix multiplications and uses carefully chosen initial phases.
  • Results: Numerical observations indicate O(d2) computational cost while the number of classical bits remains constant at 64 as degree increases.The method is demonstrated for Hamiltonian simulation, eigenstate filtering, and matrix inversion.
  • Polynomial construction: Remez exchange directly constructs minimax polynomial approximations, avoiding analytic expansions and their additional approximations.The paper reports that this can streamline phase-factor construction and produce significantly lower-degree polynomials.
  • Additional structure: Phase factors can decay away from the center, enabling phase padding that initializes high-degree factors from lower-degree solutions.The decay is linked to Chebyshev-coefficient decay.

II.2. Quantum signal processing

QSP interleaves phase rotations with signal-processing iterates to realize a target polynomial through block-encoded operators. The construction connects polynomial conditions in SU(2) to a circuit using the phase factors that specify the approximation.

  • Circuit structure: A QSP circuit consists of d iterates and d + 1 rotations arranged in an alternating sequence.The phase vector contains the parameters used by the circuit.
  • Phase conventions: The phase factors ˜φi specify the polynomial f(x) that approximates the desired function F(x).Different phase-factor conventions are related by explicit shifts and equivalence relations.
  • Operator implementation: For a non-unitary Hermitian operator, shifted phase factors ϕi implement the QSP construction, with an optional global phase discarded.The shifts differ for the final phase and for the preceding phases.
  • Block-encoding circuit: The circuit in Fig. 1 block-encodes eiϕUΠ using one signal construction with an (m + 1)-qubit Toffoli control.The three line groups represent 1, m, and n qubits.
  • Implementation cost: The resulting (n+m)-qubit unitary UeΦ can be implemented with one additional ancilla qubit.This connects the phase-factor construction to the circuit implementation.
  • Qubitization connection: Qubitization reduces matrix-function representation to an SU(2) representation on invariant subspaces indexed by eigenvalues.The upper-left element of the restricted unitary becomes a function of the corresponding eigenvalue.
  • Polynomial conditions: The QSP theorem guarantees phase factors for polynomials P and Q satisfying degree, parity, and normalization conditions on [−1, 1].The resulting phase vector has d + 1 entries.
  • Phase-factor computation: The constructive QSP theorem provides an algorithm for computing phase factors from P and Q, subject to numerical-stability considerations.The phase factors are unique up to specified trivial equivalence relations.

II.3. Representing general matrix polynomials

QSP represents general matrix polynomials by combining circuits for parity-separated and real/imaginary components. These constructions use LCU to obtain block-encodings of the desired matrix polynomial.

  • A degree-d polynomial P yields a block-encoding of P(A/α) using d queries to the block-encoding UA.The circuit also uses O((m + 1)d) other primitive quantum gates.
  • Real polynomials without definite parity are decomposed into even and odd parts, then combined with LCU.This construction gives a (2, m + 3, 0)-block-encoding of f(A/α).
  • Complex polynomials are split into real and imaginary parts, whose block-encodings are combined using LCU.The resulting construction gives a (4, m + 4, 0)-block-encoding of f(A/α).
  • The general construction assumes a Hermitian block-encoding, but can be generalized to non-Hermitian block-encodings with one additional ancilla qubit.Quantum singular value transformation offers an alternative that removes this extra ancilla and handles non-Hermitian matrices.

II.4. Direct methods for finding phase factors

The paper contrasts direct phase-factor methods based on polynomial factorization with an optimization strategy that minimizes a distance between the QSP output and the target polynomial. The optimization formulation avoids high-precision root finding and uses numerically stable SU(2) matrix operations.

  • The GSLW method reconstructs complementing polynomials through roots of 1−f(x)^2, but is unstable in standard double precision for large degrees.Its Taylor-expansion-based reduction also becomes numerically unstable as polynomial degree increases.
  • The Haah method maps polynomials to the unit circle, making it equivalent to a Chebyshev expansion that improves numerical stability over the standard basis.Both GSLW and Haah obtain phase factors through direct, single-shot calculations.
  • The optimization method minimizes L(Φ) = dist{Re[⟨0|UΦ(x)|0⟩], f(x)} instead of factoring polynomials.In practice, the distance is represented by mean squared loss over discrete sample points.
  • When the loss reaches zero, its minimizer provides the desired phase factors without constructing complementing polynomials by high-precision root finding.Gradients and Hessians involve SU(2) matrix multiplications, which the paper identifies as numerically stable.
  • Choosing Q(x) real permits symmetric phase factors under inversion, reducing the number of optimization degrees of freedom by a factor of 2.The symmetry follows from Theorem 2 and motivates symmetric phase-factor constructions.

III.2. Choice of objective function

The objective function is designed for polynomial targets and exploits inversion symmetry to reduce optimization variables. Chebyshev-node sampling connects the discrete loss to polynomial approximation accuracy, while initialization remains important because the landscape is complex.

  • Non-polynomial target functions must first be approximated by a polynomial before entering the optimization objective.Directly using a non-polynomial target generally prevents the equation L(Φ) = 0 from having a solution.
  • Imposing inversion symmetry leaves ⌈(d + 1)/2⌉ phase degrees of freedom and motivates target polynomials with parity d mod 2.The symmetric parameterization matches the number of adjustable polynomial coefficients.
  • The objective uses mean squared loss evaluated at ⌈(d + 1)/2⌉ sample points to characterize the distance between the QSP output and the target polynomial.The optimization is interpreted as polynomial interpolation in the QSP parameterization.
  • Chebyshev nodes allow the polynomial approximation error to be measured directly through the objective function.The nodes are chosen as positive roots of a Chebyshev polynomial.
  • A good initial guess is necessary because the optimization landscape is complex.The optimization is viewed as a variational quantum circuit with some gates fixed.

III.3. Generating approximation polynomials

The paper constructs approximation polynomials using Fourier–Chebyshev expansion or the Remez algorithm, then optimizes QSP phase factors with a quasi-Newton method. A carefully chosen initial point helps avoid stationary points and supports efficient phase-factor computation.

  • Polynomial approximation: Fourier–Chebyshev expansion uses FFT-based quadrature to efficiently evaluate coefficients for polynomial approximations.
  • Polynomial approximation: The Remez algorithm solves the L∞ optimal-approximation problem by alternating coefficient solving and adjustment of d+2 sample points.It can also incorporate parity constraints.
  • Initial-point selection: The optimization objective is highly non-convex, so random initialization can lead solvers to local minima.The all-zero phase initialization is a stationary point with nonzero loss.
  • Initial-point selection: A modified initial point makes Re[P(x)] = 0, avoiding the bias of the natural all-zero initialization toward Td(x).The associated initialization enabled gradient-based optimization to reach a global minimum in the reported cases.
  • Phase-factor optimization: L-BFGS optimizes symmetry-reduced phase factors without Hessian evaluation, avoiding the O(d3) cost of Newton-type optimization.The complete procedure applies this optimization to polynomial components and returns the corresponding phase factors.

IV. NUMERICAL RESULTS

Numerical tests compare the optimization method with direct phase-factor algorithms for Hamiltonian simulation and eigenstate filtering. The optimization method remains stable with double precision, while direct methods become inaccurate or inefficient at larger polynomial degrees.

  • Comparison setup: The tests compare the optimization method with GSLW and Haah using double-precision optimization versus variable-precision direct methods.The optimization method minimizes a loss function through iterative SU(2) matrix multiplication, whereas the comparison methods use direct procedures.
  • Hamiltonian simulation: For Hamiltonian simulation, the target polynomial uses degree d = 1.4|τ| + log(1/ϵ0) with ϵ0 = 10^-14.The polynomial approximates e^-iτx after truncation of the Jacobi–Anger expansion.
  • Hamiltonian simulation: The optimization method scales as τ^2, while Haah scales as τ^3 and GSLW fails to reach ϵ = 10^-12 when d exceeds 369.The optimization method uses a number of arithmetic bits independent of τ, whereas GSLW and Haah require bit counts linear in τ.
  • Hamiltonian simulation: At τ = 5000, the optimization method handles degree d = 7033 while retaining asymptotic CPU scaling of τ^2.The observed scaling agrees with the O(d^2) cost per optimization iteration.
  • Eigenstate filtering function: For eigenstate filtering, the optimization method remains stable through degree d = 10,000 with CPU time scaling as k^2.The experiments vary ∆ and choose k so that k∆ takes values from 3 to 25.
  • Eigenstate filtering function: Direct methods struggle for k∆ > 3, and Haah fails to achieve ϵ = 10^-12 when ∆ < 0.01.The comparison attributes differing direct-method performance partly to high-degree root-finding difficulty and the larger polynomial degree used by Haah.

IV.3. Matrix inversion

For matrix inversion, the optimization method computes QSP phase factors for polynomial approximations of 1/x on Dκ, including Remez constructions with odd and even parity. It remains robust through degree 4035, with computational cost scaling quadratically in κ and benefits from suitable initialization and scaling.

  • Matrix inversion: The method approximates 1/x on Dκ=[−1,−1/κ]∪[1/κ,1] using either an extended smooth function or direct Remez constructions with odd or even parity.The Remez approach directly targets the restricted interval Dκ, while the alternative first extends 1/x to [−1,1].
  • Matrix inversion: A factor of 2∼3 reduction in polynomial degree is obtained with the Remez method, and even approximations are slightly less expensive than odd ones.The lower cost of even approximations is attributed to their smaller gradient near the origin.
  • Matrix inversion: Quadratic cost in κ is observed for the optimization method, which successfully determines phase factors for the largest tested polynomial degree, d=4035.The study tests κ=10,20,…,50 and compares the optimization method with GSLW and Haah methods for lower-degree polynomials.
  • Impact of the initial point: Random initialization often becomes trapped in local minima, whereas the structured initial point (π/4,0,…,0,π/4) converges within dozens of iterations to an objective value near zero.The difficulty of finding a good solution increases with polynomial degree.
  • Sensitivity analysis: When the target polynomial satisfies ||f||∞=1−η, the Hessian condition number grows as η^−γ with γ>1 as η→0+, indicating increased optimizer sensitivity near unit norm.Scaling the target polynomial so that ||f||∞≤1−η is preferred for numerical stability and can support reuse of phase factors for slightly perturbed targets.
  • Phase-factor decay and padding: For smooth targets, phase-factor decay follows Chebyshev-coefficient decay, motivating phase padding that reuses lower-degree solutions for higher-degree approximations.Symmetric padding preserves the real part of the upper-left QSP unitary component, and the procedure incrementally increases polynomial degree.

VI. DISCUSSION

The discussion concludes that optimization-based phase evaluation efficiently and accurately supplies the phases needed to build QSP circuits for unitary representations of non-unitary operations.

  • Discussion: Optimization-based phase evaluation efficiently and accurately determines the phase factors needed to build QSP circuits for unitary representations of non-unitary operations.

Appendix A: Uniqueness of phase factors

The appendices describe equivalent QSP phase representations, circuit constructions for real and non-Hermitian block-encodings, and the Chebyshev framework used for polynomial expansions. They also show how circuit design can avoid an additional ancilla qubit in the real-polynomial construction.

  • Equivalent phase representations: GSLW’s and Haah’s phase-factor representations are equivalent through a linear transformation, with Haah’s factors unique in [−π,π)^{d+1}.The 2π-periodicity of Haah phase factors produces corresponding ±π shifts in GSLW phase factors.
  • Real-polynomial circuits: The real-polynomial QSP circuit can avoid an additional ancilla qubit by combining it with the first block-encoding ancilla.The resulting circuit provides a (1,m+1,0)-block-encoding of f(A/α).
  • Non-Hermitian block-encodings: For a non-Hermitian block-encoding of a Hermitian matrix, an additional ancilla qubit is introduced to construct a Hermitian block-encoding before applying QSP.The construction uses both UA and U†A and leads to the circuit in Fig. 17.
  • Chebyshev framework: Chebyshev polynomials of the first and second kinds form the polynomial framework, with the first kind defined by Td(x)=cos(d arccos(x)).Their roots are Chebyshev nodes, and the Chebyshev expansion is related to a Fourier series after substituting x=cos θ.

Appendix E: Remez Method

The Remez appendix formulates best uniform approximation over parity-restricted function bases satisfying the Haar condition. Its exchange method generates polynomials converging uniformly to the optimal approximation.

  • Remez Method: The approximation problem minimizes the L∞ norm over a parity-restricted basis rather than the full monomial basis.The basis dimension is chosen as N=⌈(d+1)/2⌉ for a degree-d polynomial with definite parity.
  • Remez Method: Bases such as odd or even Chebyshev polynomials satisfy the Haar condition on suitable domains, ensuring linear independence at every set of N points.Examples include gj(x)=T2j−1(x) and gj(x)=T2j−2(x) on X⊂(0,1].
  • Remez Method: The Remez exchange method produces approximation polynomials that converge uniformly to the optimal polynomial, generally linearly and sometimes quadratically.The improved quadratic rate applies for a large range of functions.

Appendix F: L-BFGS Algorithm

This appendix introduces BFGS and L-BFGS, then outlines a Remez method for computing a best-approximation polynomial through iterative residual maximization.

  • L-BFGS background: BFGS is a quasi-Newton method for unconstrained optimization that stores and updates a dense inverse-Hessian approximation.Its approximation is updated by a rank-two update using gradient information.
  • L-BFGS background: L-BFGS approximates BFGS while using a limited amount of computer memory.
  • Remez method: The Remez method takes an interval, target function, Haar-condition basis, and N + 1 initial points as input.
  • Remez method: Each iteration solves for polynomial coefficients and a residual, locates residual roots, and maximizes signed residuals between consecutive roots.
  • Remez method: The method replaces an initial point when needed, repeats the exchange process, and outputs an approximation to the best-approximation polynomial.

Appendix G: Implementation details of the direct methods for finding phase factors

This appendix describes direct GSLW and Haah implementations for constructing QSP phase factors, including numerical preprocessing, complementary-polynomial construction, and matrix-factor decomposition. The procedures produce phase-factor representations with approximation error bounded by ϵ.

  • Implementation choices: The implementations use Julia and PolynomialRoots for polynomial roots, while GSLW additionally rounds near-real roots and computes coefficients in the Chebyshev basis with FFT.The root tolerances use machine precision ϵ = 10^-16 in the implementation.
  • Phase-factor extraction: For direct phase-factor extraction, complementary polynomials B and C are chosen so that P and Q satisfy |P(x)|^2 + (1 − x^2)|Q(x)|^2 = 1 on [−1, 1].
  • Phase-factor extraction: The phase factors are obtained recursively from leading coefficients, followed by a final phase choice and fixed alternating phases.The leading coefficients determine φt through e^(2iφt) = pt/qt−1.
  • Haah method: The Haah method uses real-on-circle Laurent polynomials with bounded magnitude and definite parity, and can set B(z) = 0 when the second polynomial is not needed.
  • GSLW method: GSLW preprocessing rationally approximates scaled Laurent-polynomial coefficients, removes sufficiently small coefficients, preserves parity, and finds all roots of 1 − a(z)^2 − b(z)^2.The threshold is ϵ/(30d), and the resulting polynomials are stored with rational coefficients.
  • GSLW method: The GSLW construction forms complementary polynomials from factorized roots, computes matrices by FFT, and iteratively generates E0, E1(z), ..., E2n(z).
  • GSLW method: The resulting matrices satisfy an approximation bound of ϵ for A(z^2) + iB(z^2) over the unit circle.

Appendix H: Proof of Theorem 5

This appendix proves Theorem 5 by expanding symmetrically partitioned QSP phase factors and analyzing the resulting unitary products. The proof treats even and odd degree cases separately.

  • Even degree case: The proof first handles d′ = 2d + 2 by reducing the claim to an equivalent statement and deriving bounds on the phase-factor expansion.
  • Symmetric decomposition: The QSP phase factors are divided into symmetric left and right groups, allowing the unitary to be expressed as a product of two QSP matrices.
  • Expansion: Each QSP unitary is rewritten using tangent variables tj := tan(φj), which are treated as O(φj) in the expansion.
  • Expansion: Only odd orders in the phase factors contribute to the relevant imaginary component, using the transpose relation between the symmetric unitary factors.
  • Parity cases: The argument establishes Eq. (39) for even d′ and then begins the corresponding proof for d′ = 2d + 1.
  • Odd degree case: For odd d′, the proof again defines symmetric left and right phase-factor groups before applying the analogous expansion.
Loading 2002.11649v2…