PaperlayerView exact revision on arXiv

Source-linked AI summary

Product Decomposition of Periodic Functions in Quantum Signal Processing

Jeongwan Haah

arXiv:1806.10236v4quant-ph

TL;DR

Computing long sequences of interspersing rotations for quantum signal processing was numerically unstable, while earlier polynomial-time claims did not analyze that instability. The algorithm approximates the target function as a matrix element of an SU(2)-valued polynomial function and decomposes that function into primitive factors using bounded-precision computation. O(N^3 polylog(N/ϵ)) random-access-machine time and O(N log(N/ϵ)) bits of precision suffice for a reliable decomposition.

  • Problem

    Computing long sequences of interspersing rotations for quantum signal processing was numerically unstable, while earlier polynomial-time claims did not analyze that instability.

    Source: §1 · pp. 1–2

  • Method

    The algorithm approximates the target function as a matrix element of an SU(2)-valued polynomial function and decomposes that function into primitive factors using bounded-precision computation.

    Source: §1 · p. 2, §4 · p. 8

  • Results

    O(N^3 polylog(N/ϵ)) random-access-machine time and O(N log(N/ϵ)) bits of precision suffice for a reliable decomposition.

    Source: §1 · p. 2, §4.3 · p. 14

  • Takeaways & Limitations

    For any Laurent polynomial function F:U(1)→SU(2), t↦F(t²) belongs to some Eₙ and has a unique product decomposition, covering prior achievable-function results [2] [4].

    Source: §3 · p. 4

    The algorithm is not numerically stable in the usual sense: it can require poly(N log(1/ϵ)) rather than polylog(N/ϵ) bits, with instability arising during early stages of iterative degree reduction.

    Source: §7 · p. 19

Abstract

from arXiv · show

We consider an algorithm to approximate complex-valued periodic functions $f(e^{i\theta})$ as a matrix element of a product of $SU(2)$-valued functions, which underlies so-called quantum signal processing. We prove that the algorithm runs in time $\mathcal O(N^3 \mathrm{polylog}(N/\epsilon))$ under the random-access memory model of computation where $N$ is the degree of the polynomial that approximates $f$ with accuracy $\epsilon$; previous efficiency claim assumed a strong arithmetic model of computation and lacked numerical stability analysis.

1 Introduction

Quantum signal processing uses classical preprocessing to find interspersing single-qubit rotations, but numerical instability left the efficiency of computing long sequences unclear. The paper refines the decomposition algorithm and bounds its precision and running time.

1 Introduction

  • Quantum signal processing transforms eigenvalues through a matrix element and has supported gate-efficient Hamiltonian simulation and been proposed for linear equations [3, 4, 9].The cited work also describes applications to sparse black-box and explicitly described local Hamiltonian simulation [5] [6] [7].Source: §1 · p. 1
  • Prior preprocessing instability made it unclear whether long sequences of interspersing rotations could be computed efficiently; one report called sequences of length ≳30 prohibitive for Jacobi–Anger expansions [7, App. H.3].The usefulness of quantum signal processing depends on computing long sequences of these rotations.Source: §1 · p. 1
  • Earlier methods expanded large-degree polynomials found from another polynomial’s roots, risking instability in both convolution computations and precision tracking from root distributions.These numerical problems had not previously been addressed.Source: §1 · p. 2
  • O(N^3 polylog(N/ϵ)) time on a sequential random-access machine and O(N log(N/ϵ)) bits of precision suffice for a reliable decomposition.The algorithm avoids too-small numbers by sacrificing initial approximation quality and replaces polynomial expansions with a Fourier transform.Source: §1 · p. 2
  • The refined algorithm establishes that classical preprocessing can be done in polynomial time on a Turing machine, with O(N log(N/ϵ)) bits sufficient for the desired accuracy.N is the degree of the polynomial approximating the transformation function to additive error ϵ.Source: §1 · p. 2
  • A rudimentary Mathematica experiment produced a Jacobi–Anger sequence longer than 2,000 in a few hours; its empirical time complexity appeared quadratic in N.The authors report that reduced precision and a lower-level implementation later expedited the computation by orders of magnitude.Source: §1 · p. 2

2 Quantum Signal Processing

Controlled applications of W and W† transfer an eigenphase to the control qubit, where single-qubit operations shape the resulting transformation. Alternating these operations yields an SU(2)-valued function whose matrix element can implement the target transformation after post-selection.

2 Quantum Signal Processing

  • On a fixed eigenstate of W, control-W kicks its eigenvalue eiθ back to the control qubit, inducing a phase operation there.Conjugating control-W by a single-qubit unitary lets the control basis be chosen arbitrarily.Source: §2 · p. 3
  • Alternating an even number n of control-W and control-W† operations with single-qubit changes of basis yields a product of primitive matrices, with t = eiθ/2 and rank-1 projector P = I − Q.The inverse of control-W is assumed implementable in this construction.Source: §2 · p. 3
  • The product defines an SU(2)-valued function F(t); if ⟨+|F(t)|+⟩ approximates f(eiθ), post-selection on |+⟩ enacts the target transformation.The post-selection success probability depends on |f(eiθ)|, and the choice of |+⟩ is a convention; the square-root branch must be used consistently.Source: §2 · p. 3

3 Polynomial functions U(1) →SU(2)

Section 3 characterizes polynomial SU(2)-valued functions on the unit circle through unique factorization, then describes parity constraints and conditions for completing scalar polynomials into matrix elements.

3 Polynomial functions U(1) →SU(2)

  • Every F(t) in E_n of degree n has a unique decomposition into primitive matrices and a unitary, characterizing polynomial U(1) → SU(2) functions and covering prior achievable-function results [2] [4].The characterization is slightly more general than prior results, which implicitly assume Tr(PjZ) = 0.Source: §3 · p. 4

3.1 Parity constraints

  • Parity forces the factored projectors onto the Bloch-sphere equator, but residual E0 can retain a Z component and omitting it can incur quadratic accuracy loss.Under the stated parity constraint, F(t) is specified by n + 1 angles, including the residual phase.Source: §3.1 · pp. 5–6

3.2 Complementing polynomials

  • If a and b are real-on-circle Laurent polynomials of degree at most n, with a² + b² < 1 on U(1) and reciprocal, pure complements c and d exist with a² + b² + c² + d² = 1.The complements are not generally unique; different eligible root lists can yield them.Source: §3.2 · pp. 7–8
  • The reciprocity conditions are proof-driven: another guarantee that complements exist could enlarge the algorithm’s input scope, while flexible quantum signal processing can combine even and odd functions.Any periodic function is the sum of an even and an odd function.Source: §3.2 · p. 7

4 Efficient implementation with bounded precision

The section presents a bounded-precision algorithm that constructs an SU(2)-valued approximation and decomposes it into primitive single-qubit operations. It establishes an O(N^3 polylog(N/ε)) running time in the random-access memory model, with precision and error controlled through the construction.

4 Efficient implementation with bounded precision

  • The algorithm first constructs an SU(2)-valued function whose matrix element approximates the input, then decomposes it into primitive matrices.The construction is tailored to reduce and track numerical error.Source: §4 · p. 8
  • The input is a Fourier-coefficient list for A(e^iϕ)+iB(e^iϕ), with definite parity for each component and A²+B² ≤ 1.The input function must also be sufficiently close to the ultimate target, whose approximation is analyzed separately.Source: §4 · p. 8

4.1 Algorithm

  • Rational coefficient approximation scales the input by 1−10ε, rounds coefficients to additive error ε/N, and replaces smaller coefficients with zero.This can reduce the Laurent-polynomial degree.Source: §4.1 · p. 9
  • The algorithm finds roots of 1−a(z)²−b(z)², evaluates complementary polynomials on a unit-circle grid, then uses an FFT and sequential coefficient updates to recover primitive factors.The grid uses a power-of-two size D greater than 2n+1.Source: §4.1 · p. 9

4.2 Further details and analysis

  • The rational coefficients of a(z) and b(z) are stored exactly, while the algorithm uses R = Θ(N log(N/ε)) bits of arithmetic precision.The precision choice satisfies the assumptions used in the error estimates.Source: §4.2 · pp. 10, 14, §4.3 · p. 14
  • For a real-on-circle Laurent polynomial bounded between m and M, every zero is at least m/(4Md²) from the unit circle.This separation supports bounded-precision root handling in the algorithm.Source: §4.2 · p. 11

4.3 Computational complexity

  • The root-finding stage takes O(N³ log(N/ε)) time, and exponentially increasing precision can avoid choosing a tailored bit count in advance.The repeated-precision schedule has total time within a constant multiple of its last round.Source: §4.3 · p. 14
  • O(N^3 polylog(N/ε)) is the overall running time under the random-access memory model.The bound includes high-precision arithmetic costs.Source: §4 · p. 8, §4.3 · p. 14

5 Application to Hamiltonian simulation

The section applies quantum signal processing to Hamiltonian simulation by approximating the desired eigenvalue transformation with a truncated Jacobi-Anger expansion. It gives a truncation condition for ε-closeness and reports that measured running time scales cubically with τ.

5 Application to Hamiltonian simulation

  • Quantum signal processing implements e^(−iτH) by transforming eigenvalues of a unitary associated with H; the product requires at least 2τ factors.The post-selection success probability is close to 1 because the target function has unit magnitude on the eigenvalues.Source: §5 · p. 15
  • The Jacobi-Anger expansion for exp(iτ sin ϕ) converges absolutely at a superexponential rate.The analysis uses steepest descent to bound its truncation.Source: §5 · p. 15
  • N = Ω(|τ| + log(1/ε)) makes the truncated expansion ε-close to the full expansion, and measured running time scales cubically with τ.The running-time measurements used two fixed values of ε.Source: §5 · p. 16

6 Application to Matrix inversion

For matrix inversion, the paper uses a polynomial eigenvalue transformation that approximates 1/sin(θλ) when |sin(θλ)| is bounded away from zero by 1/κ. The construction yields a Laurent polynomial and an inversion algorithm whose gate complexity accounts for amplitude amplification.

6 Application to Matrix inversion

  • The plotted running-time fits have exponents 3.11 for ϵ = 10^-9 and 3.17 for ϵ = 10^-4 as functions of τ.The top-right data point uses 2172 primitive matrices in the decomposition.Source: Figure 1 · p. 17
  • The inversion target is the eigenvalue map ∓e^±iθλ ↦ 1/sin(θλ), for a norm-1 Hermitian H block-encoded in W with eigenvalues bounded away from zero by 1/κ.Strictly zero eigenvalues are also allowed when seeking a pseudo-inverse.Source: §6 · p. 17
  • The polynomial filter is close to 1 when |sin ϕ| ≥ 1/κ and vanishes at sin ϕ = 0, allowing the near-zero region to be excluded from the inversion approximation.The text says this function can be replaced by a lower-degree polynomial for large b.Source: §6 · p. 18
  • For |sin ϕ| ≥ 1/κ, the construction yields a genuine Laurent polynomial; its degree improves on [24] by a factor of log κ through a different normalization.The paper approximates 1/(κx) rather than 1/x and describes the analysis difference as using exponential rather than trigonometric functions.Source: §6 · p. 18
  • Choosing b′ = ⌈κ ln(8/ϵ)⌉ gives success probability as small as Ω(1/(κ log(1/ϵ))^2), requiring amplitude amplification with gate-complexity factor O(κ log(1/ϵ)).Overall gate complexity is proportional to the polynomial degree times the number of amplitude-amplification iterations.Source: §6 · p. 19

7 Discussion

The decomposition algorithm is not numerically stable in the usual sense, and its iterative degree-reduction step can require more precision than polylogarithmic in N/ϵ. The discussion identifies small leading coefficients as the source of error while noting that the analysis may overstate instability in some cases.

7 Discussion

  • The decomposition algorithm is not numerically stable in the usual sense: it requires poly(N log(1/ϵ)) rather than polylog(N/ϵ) bits of precision.The discussion says the instability appears unavoidable for methods that reduce polynomial degree one at a time, at least early in the reduction.Source: §7 · p. 19
  • Small leading coefficients cause numerical error, but may grow rapidly during Step 5, making the instability analysis potentially loose.The paper gives a schematic projector identity whose decomposition is less unstable than the worst case and says similar situations may occur for some inputs.Source: §7 · p. 19

A Jordan’s Lemma and block encoding of Hamiltonians

Jordan’s lemma decomposes a space acted on by two self-adjoint projectors into common invariant subspaces of dimension one or two. Applied to reflections, it gives a block structure relevant to Hamiltonian simulation.

A Jordan’s Lemma and block encoding of Hamiltonians

  • For arbitrary self-adjoint projectors P and Q, the space decomposes into orthogonal subspaces invariant under both, each of dimension one or two.The proof constructs a nonzero common invariant subspace of dimension at most two and proceeds by induction.Source: A Jordan’s Lemma and block encoding of… · pp. 19–20
  • For two Hermitian unitaries, a basis makes U1 diagonal and U1U2 block-diagonal with blocks of dimension at most two.Each block belongs to U(1) or U(2).Source: A Jordan’s Lemma and block encoding of… · p. 20
  • In the stated block encoding, the product of reflections has eigenvalues ±e^∓iθ with sin θ = λ, connecting the block structure to Hamiltonian eigenvalues.The example takes H = (⟨G| ⊗ I) U2 (|G⟩ ⊗ I) and U1 = (2|G⟩⟨G| − I) ⊗ I.Source: A Jordan’s Lemma and block encoding of… · p. 20
Loading 1806.10236v4…