Halving the Cost of Quantum Algorithms with Randomization
John M. Martyn Patrick Rall



 
Center for Theoretical Physics, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA
IBM Quantum, MIT-IBM Watson AI Lab, Cambridge, Massachusetts 02142, USA

 
IBM Quantum, MIT-IBM Watson AI Lab, Cambridge, Massachusetts 02142, USA

§ ABSTRACT 
    Quantum signal processing (QSP) provides a systematic framework for implementing a polynomial transformation of a linear operator, and unifies nearly all known quantum algorithms. In parallel, recent works have developed randomized compiling, a technique that promotes a unitary gate to a quantum channel and enables a quadratic suppression of error (i.e., ϵ→ O(ϵ^2)) at little to no overhead. Here we integrate randomized compiling into QSP through Stochastic Quantum Signal Processing. Our algorithm implements a probabilistic mixture of polynomials, strategically chosen so that the average evolution converges to that of a target function, with an error quadratically smaller than that of an equivalent individual polynomial. Because nearly all QSP-based algorithms exhibit query complexities scaling as O(log(1/ϵ))—stemming from a result in functional analysis—this error suppression reduces their query complexity by a factor that asymptotically approaches 1/2. By the unifying capabilities of QSP, this reduction extends broadly to quantum algorithms, which we demonstrate on algorithms for real and imaginary time evolution, phase estimation, ground state preparation, and matrix inversion.



MIT-CTP/5756
Halving the Cost of Quantum Algorithms with Randomization
    Patrick Rall
    September 9, 2024
=========================================================





§ INTRODUCTION
 

Classical randomness plays a pivotal role in the design of quantum protocols and algorithms. In the near-term, randomized benchmarking <cit.> is central to calibrating and assessing the quality of quantum gates, and quasi-probability methods like probabilistic error cancellation and noise twirling can help reduce noise <cit.>. Similarly, random circuit sampling is central to quantum supremacy experiments <cit.>, and randomized measurements provide a powerful probe into the properties of complex quantum systems <cit.>. As we progress towards quantum advantage and early fault-tolerant quantum hardware, many lines of research aim to reduce the requirements of traditional quantum algorithms by incorporating classical randomness <cit.>.

Randomized compiling is a key example of leveraging classical randomness to improve quantum computation <cit.>. As its name suggests, this process randomly compiles gates at execution time, or equivalently, promotes a unitary gate to a quantum channel that is a probabilistic mixture of unitaries. Remarkably, randomized compiling can quadratically suppresses gate errors without increasing the cost of circuit synthesis. Yet, applications of this technique to quantum algorithms have so far been restricted to Trotterized Hamiltonian simulation <cit.> and phase estimation <cit.>, leaving a vacuum of applications to other algorithms.  

In an effort to fill this gap of randomized quantum algorithms, we propose using quantum signal processing (QSP) <cit.> as a medium for achieving widespread advantage of randomized compiling. QSP prepares a polynomial transformation of a linear operator, and has been shown to encompass nearly all quantum algorithms, from Hamiltonian simulation and quantum search, to matrix inversion and fast integer factoring <cit.>. 




In this work, we achieve exactly this goal. We merge randomized compiling with QSP by developing Stochastic Quantum Signal Processing (Stochastic QSP). By virtue of randomized compiling, our construction quadratically suppresses the error in a QSP polynomial approximation of a target function. To study how this suppression impacts the cost of QSP-based algorithms, we show that an elementary result in the approximation of smooth functions implies that nearly all QSP-based algorithms achieve a query complexity that scales with the error ϵ as O(log(1/ϵ)), which we also empirically confirm. Hence the quadratic suppression of error afforded by stochastic QSP translates to an asymptotic halving of the cost of a QSP-based algorithm over their deterministic counterparts (asymptotic in the limit of log(1/ϵ) dominating the cost). In realizing this cost reduction, we “combine the strengths of QSP and randomization," as hypothesized in Ref. <cit.>. 



An outline of this work is as follows. After reviewing QSP and other preliminary topics in Sec. <ref>, we present stochastic QSP in Sec. <ref>. We then demonstrate the versatility of our scheme by showing its compatibility with several generalizations and variants of QSP in Sec. <ref>. Finally, in Sec. <ref>, we benchmark the performance of stochastic QSP for various polynomials relevant to quantum algorithms, including those for real and imaginary time evolution, phase estimation, ground state preparation, and matrix inversion. 
In Appendix <ref> we review some useful results on Fourier and Chebyshev expansions of smooth functions, and in Appendix <ref> we prove some extensions of randomized compilation techniques.




§ PRELIMINARIES

We begin by discussing the preliminary topics of this work, including notation and background concepts (Sec. <ref>), quantum signal processing (QSP) (Sec. <ref>), polynomial approximations to smooth functions (Sec. <ref>), and randomized compiling (Sec. <ref>). 




 §.§ Notation and Background Concepts

We will study functions F(x) on the domain x ∈ [-1,1], where it will be convenient to define the function norm

    F _[-1,1] := max_x∈[-1,1] |F(x)| .

Of particular interest will be functions bounded as F _[-1,1]≤ 1. 

A convenient set of functions on this domain are the Chebyshev polynomials. The order n Chebyshev polynomial is defined as T_n(x) = cos(n arccos(x)) for integer n≥ 0. It can be shown that T_n(x) is a polynomial of degree n with parity n  2 (i.e., either even or odd, depending on n) and bounded magnitude T_n_[-1,1] = 1. An important property of the Chebyshev polynomials is that they furnish an orthogonal basis in which an arbitrary function on x∈ [-1,1] can be expanded: 

    F(x) = c_0/2 + ∑_n=1^∞ c_n T_n(x) ,     c_n = 2/π∫_-1^1 F(x) T_n(x)/√(1-x^2) dx ,

where c_n are the Chebyshev coefficients. If F(x) is a degree d polynomial, then this series terminates at order n=d. 


In this work we will also study unitary and non-unitary transformations. We will denote an operator by a Latin character, say A, and the associated channel by the corresponding calligraphic character: 𝒜(ρ) = A ρ A^†. 

In analyzing these operators, we will consider the spectral norm (equivalently, the operator norm), defined as

    A  = sup_|ψ⟩ A |ψ⟩ ,

where the supremum is taken over normalized states ⟨ψ | ψ⟩ = 1. It can be shown that this norm equates to the maximum singular value of A. We will also use the trace norm (equivalently, the Schatten 1-norm)

    A _1 = ( √(A^† A)) ,

which equates to the sum of singular values of A. 


This norm is associated with the trace distance between two density matrices:

    d_tr (ρ,σ) = 1/2ρ - σ_1 .

This bounds the discrepancy in expectation values as | (ρ O) - (σ O) | ≤ O ρ - σ_1  = 2  O  d_tr (ρ,σ). 


Another relevant metric is the diamond norm, defined for a channel ℰ as:

    ℰ_♢ = sup_ρ (ℰ⊗ℐ)(ρ) _1 ,

where this supremum is taken over normalized density matrices ρ that live in a possibly enlarged Hilbert space, and ℐ is the identity channel ℐ(ρ) = ρ. The diamond norm induces the diamond distance between two channels:

    d_♢ (ℰ, ℱ)  :   = 1/2ℰ - ℱ_♢
       = 1/2sup_ρ (ℰ⊗ℐ)(ρ) - (ℱ⊗ℐ)(ρ) _1.

For channels 𝒜(ρ) = A ρ A^† and ℬ(ρ) = B ρ B^† with spectral norms A, B ≤ 1, the diamond distance is upper bounded as (see Lemma 4 of Ref. <cit.> for proof):

    d_♢ (𝒜, ℬ) = 1/2𝒜 - ℬ_♢≤ A - B  .






 §.§ QSP

The quantum signal processing (QSP) algorithm is a systematic method of implementing a polynomial transformations on a quantum subsystem <cit.>. QSP works by interleaving a signal operator U, and a signal processing operator S, both taken to be SU(2) rotations about different axes. Conventionally, the signal operator is an x-rotation through a fixed angle and the signal processing operator a z-rotation through a variable angle ϕ:

    U(x) = [         x i√(1-x^2); i√(1-x^2)         x ],      
        S(ϕ) = e^iϕ Z.

Then, by selecting a set of d+1 QSP phases ϕ⃗ = (ϕ_0, ϕ_1, ... , ϕ_d) ∈ℝ^d+1, one can construct the following QSP sequence as an interleaved product of U and S, whose matrix elements are polynomials in x:

    U_ϕ⃗(x)    = S(ϕ_0) ∏_i=1^d  U(x)  S(ϕ_i) 
       = [            P(x)   iQ(x)√(1-x^2); iQ^*(x)√(1-x^2)          P^*(x) ],

where P(x) and Q(x) are polynomials parameterized by ϕ⃗ that obey:

    1.  deg(P) ≤ d,   deg(Q) ≤ d-1 
       2.  P(x) has parity  d  2,  and  Q(x) has parity 
          (d-1) 2 
       3.  |P(x)|^2 + (1-x^2) |Q(x)|^2 = 1,  ∀ x ∈ [-1,1].


This result implies that one can prepare polynomials in x by projecting into a block of U_ϕ⃗, e.g. ⟨ 0| U_ϕ⃗ | 0⟩ = P(x). While this class of polynomials is limited by the conditions of Eq. (<ref>), one can prove that by projecting into other bases (e.g., the | + ⟩⟨ + | basis), and incorporating linear-combination-of-unitaries circuits <cit.>, QSP can encode an arbitrary degree-d polynomial that need only obey P_[-1,1]≤ 1 <cit.>. For any such polynomial, the corresponding QSP phases ϕ⃗ can be efficiently determined classically <cit.>, thus amounting to merely a pre-computation step. In addition, as per Eq. (<ref>), the cost of realizing such a degree-d polynomial is d queries to U(x).


Remarkably, QSP can be generalized to implement polynomial transformations of linear operators through its extension to the quantum eigenvalue transformation (QET) <cit.> and quantum singular value transformation (QSVT) <cit.>. This is achieved analogous to QSP: provided access to an unitary that block-encodes a linear operator A, we can construct an operation that encodes a polynomial transformation P(A):

    U[A]=[ A ·; · · ]↦ U_ϕ⃗[A] = 
        [ P(A)    ·;    ·    · ] ,

where the unspecified entries ensure unitarity. Paralleling Eq. (<ref>), this operation is constructed as an interleaved sequence of U[A] and parameterized rotations. In essence, this applies QSP within each eigenspace (or singular value space) of A, such that that the resulting sequence encodes a degree-d polynomial transformation P(A) acting on the eigenvalues (singular values) of A. The cost of realizing the polynomial P(A) is d queries to the block-encoding, translating to a runtime O(d). Lastly, while Eq. (<ref>) caters to an encoding in the |0⟩⟨ 0| matrix element, one can more generally take A to be accessed by orthogonal projectors Π, Π' as A=Π U[A] Π'. To wit, Eq. (<ref>) corresponds to the conventional choice Π = Π' = |0⟩⟨ 0 | ⊗ I.

QET and QSVT are powerful algorithms, shown to unify and simplify most known quantum algorithms, while maintaining near-optimal query complexities <cit.>. Explicitly, an algorithm can be cast into the language of QET/QSVT by constructing a polynomial approximation to a matrix function that solves the problem of interest. For instance, in Hamiltonian simulation, one can design a polynomial P(H) ≈ e^-iHt to simulate time evolution <cit.>. Algorithms encompassed in this framework include the primoridal algorithms of search, simulation, and phase estimation <cit.>, as well as more modern algorithms, like matrix inversion <cit.>, and ground state preparation <cit.>.

Finally, throughout this work, we use the term “QSP" in place of “QET" and “QSVT" for simplicity, following the conventional parlance. However, this should be understood to be QET/QSVT when acting on a linear operator rather than a scalar.  







 §.§ Polynomial Approximations to Smooth Functions


As emphasized above, the utility of QSP lies in generating matrix functions without the need to unitarily diagonalize the underlying matrix. Specifically, QSP enables the approximation of a matrix function F(A), while remaining agnostic to the eigenvalues of A. This is achieved by selecting a polynomial approximation to the target function P(x) ≈ F(x) and implementing P(A) with QSP, where the accuracy of this approximation can be tuned by increasing the degree of P(x). Because the cost of a QSP-based algorithms scales with the polynomial degree, their complexity rests on results in approximation theory. 

To make this connection concrete, let us consider a few examples encountered in QSP-based algorithms. First, consider the decaying exponential function e^-β (x+1) for a parameter β>0, rescaled so that e^-β (x+1)_[-1,1] = 1. In QSP, this function is employed to prepare thermal states and estimate partition functions at inverse temperature β <cit.>. It is well established that e^-β (x+1) can be approximated to within additive error ϵ over x∈[-1,1] by a polynomial of degree d = O(√(β)log(1/ϵ)) <cit.>. As anticipated, the degree increases with decreasing error ϵ and increasing β. 

Similarly, consider the error function erf(kx) = 2/√(π)∫_0^kx e^-t^2 dt for a parameter k>0, which is naturally bounded as erf(kx)_[-1,1] < 1. This function is used in QSP to estimate the step function by using a large value of k, with notable applications to phase estimation <cit.> and ground state preparation <cit.>. Prior work has proven that erf(kx) can be approximated to within additive error ϵ over x∈[-1,1] by a polynomial of degree d = O(k log(1/ϵ)) <cit.>. As before, the degree grows with decreasing error and increasing k. 



Observe that in both of these examples, the degree of the polynomial approximation scales with the error as O(log(1/ϵ)). This scaling is a generic feature of polynomial approximations to smooth functions. To understand this phenomenon, consider expanding a function on x∈[-1,1] in the basis of Chebyshev polynomials as in Eq. (<ref>): 

    F(x) = c_0/2 + ∑_n=1^∞ c_n T_n(x).
 
As we prove in Appendix <ref>, if F(x) is C^∞ function (i.e., continuous and infinitely differentiable), then its Chebyshev coefficients decay super-polynomially as |c_n| = e^-O(n^r) for some exponent r > 0.

For a large class of smooth functions, it is found that r=1 <cit.>, such that |c_n| = e^-O(n) decays geometrically. In this case, a truncation of the Chebyshev series at order d, P(x) = ∑_n=0^d c_n T_n(x), furnishes a degree d polynomial approximation to F(x) that suffers error 

    max_x∈ [-1,1]| P(x) - F(x) | = max_x ∈ [-1,1]|∑_n=d+1^∞ c_n T_n(x) | 
        ≤∑_n=d+1^∞ |c_n| = ∑_n=d+1^∞ e^-O(n) = e^-O(d).

Hence, to guarantee an error at most ϵ, it suffices to choose a degree d = O(log(1/ϵ)).

In practice many polynomial approximations are constructed via truncated Chebyshev series. This includes polynomial approximations to a wide range of functions relevant to quantum algorithms, including the decaying exponential e^-β x, trigonometric functions sin(tx), cos(tx), the step function Θ(x), the inverse function 1/x,[Although the step function Θ(x) and the inverse function 1/x exhibit singularities at x=0, and thus are not C^∞ functions, they can however be approximated by C^∞ functions by excluding a small region around their singularity. This strategy is used in practice, and renders these function amenable to results on polynomial approximations to smooth functions.] and beyond. Accordingly, these approximations all exhibit degrees that scale with the error as d = O(log(1/ϵ)), which carries over to the complexity of their corresponding QSP-based algorithms. 







 §.§ Randomized Compiling and the Mixing Lemma

In order to incorporate randomization into QSP, we will use the concept of randomized compiling. Formally introduced in Ref. <cit.>, randomized compiling can be understood as follows. In a quantum algorithm, one repeatedly executes a quantum circuit and measures the output to obtain useful information. For instance, in Hamiltonian simulation, a quantum state is repeatedly time-evolved and then measured to extract an expectation value. In randomized compiling, instead of executing the same circuit at each iteration, one executes a circuit sampled from a distribution, where each sample is drawn independently. This process can be viewed as replacing a unitary operation with a a quantum channel that is a probabilistic mixture of unitaries. 

Remarkably, if this mixture is chosen strategically, randomized compiling enables a quadratic suppression of error: if an individual gate approximates a target unitary with error ϵ, the randomly compiled channel can approximate the corresponding target channel with error O(ϵ^2). This error suppression is achieved at little to no increase in overhead, requiring only the ability to classically sample a distribution and implement gates on the fly. From a physical point of view, randomized compiling achieves this quadratic suppression by turning coherent errors into incoherent errors: whereas N coherent errors can add constructively to O(N), incoherent errors essentially perform a random walk and collectively average to O(√(N)) <cit.>. 

The precise error suppression achieved by randomized compiling is quantified by the Hastings-Campbell mixing lemma, independently proven in Refs. <cit.> in the context of gate synthesis. 

    Let V be a target unitary operator, and 𝒱(ρ) = V ρ V^† the corresponding channel. Suppose there exist m unitaries { U_j }_j=1^m and an associated probability distribution p_j that approximate V as
    
    U_j - V ≤ a  for all  j,
       ∑_j=1^m p_j U_j - V ≤ b ,

    for some a,b > 0. Then, the corresponding channel Λ(ρ) = ∑_j=1^m p_j U_j ρ U_j^† approximates 𝒱 as
    
    Λ - 𝒱_♢≤ a^2 + 2b .



To see how this enables a quadratic suppression of error, suppose that one can determine an ensemble of unitaries { U_j } that each achieve spectral error a = ϵ, and a distribution p_j such that b=O(ϵ^2). Then, whereas an individual unitary U_j suffers diamond norm error O(ϵ), the mixing lemma establishes that the channel Λ achieves a quadratically suppressed diamond norm error ≤ a^2 + 2b = O(ϵ^2). Importantly, because Λ is a probabilistic mixture of the unitaries { U_j }, the cost of simulating Λ is no more expensive than the cost of sampling p_j and implementing an individual unitary U_j. 


The mixing lemma has been leveraged to improve a variety of quantum protocols via randomized compiling. Noteworthy examples include reducing the cost of gate synthesis <cit.>, tightening fault-tolerance thresholds for general noise models <cit.>, and enhancing the precision of state preparation <cit.>. On the algorithmic side, the mixing lemma has been merged with Trotterization to significantly reduce the complexity of chemical simulations <cit.>, double the order of Trotter formulae <cit.>, and accelerate imaginary time evolution <cit.>. Here we continue this campaign by extending randomized compiling to quantum signal processing. As QSP unifies nearly all quantum algorithms <cit.>, this opens the floodgates to a new suite of randomized quantum algorithms with reduced query complexities.




§ STOCHASTIC QUANTUM SIGNAL PROCESSING

Our goal is to integrate randomized compiling into QSP, and thereby establish a framework for designing randomized quantum algorithms that achieve reduced query complexities. We begin in Sec. <ref> by introducing an extension of the mixing lemma for operators block-encoded in unitaries. Then, in Sec. <ref>, we use this result to develop stochastic QSP: we replace a single QSP polynomial with a channel that is a probabilistic mixture of QSP polynomials, each strategically crafted to exploit the mixing lemma and quadratically suppress error. As we show, this furnishes randomized QSP-based algorithms with roughly half the cost of their deterministic counterparts.




 §.§ The Mixing Lemma for Block-Encodings

As emphasized in Sec. <ref>, QSP polynomials are constructed as block-encodings. That is, a QSP polynomial P(A) is encoded in a block of a higher dimensional unitary U, and accessed as P(A) = Π U Π' for some orthogonal projectors Π, Π'. Conventionally, the projectors are taken to be Π = Π' = |0⟩⟨ 0 | ⊗ I, such that P(A) is encoded in the |0⟩⟨ 0 | block of the unitary. 


To apply the mixing lemma to QSP, it is therefore necessary to establish a variant of the mixing lemma for operators block-encoded in unitary transformations. For the sake of simplicity, we present this theorem for an operator encoded in the |0⟩⟨ 0| block of a unitary:
 
    Let V be a unitary that block-encodes a (possibly non-unitary) target operator S as S = (⟨ 0 | ⊗ I) V ( |0 ⟩⊗ I). Suppose there exist m unitaries { U_j }_j=1^m that block-encode operators R_j as R_j = (⟨ 0 | ⊗ I)  U_j ( | 0 ⟩⊗ I). Also suppose there exists an associated probability distribution p_j such that
    
    R_j  - S ≤ a  for all  j , 
       ∑_j=1^m p_j R_j - S ≤ b .

    Then, the corresponding unitary channel Λ(ρ) = ∑_j=1^m p_j U_j ρ U_j^† approximates the action of the channel 𝒱 as
    
    Λ̅ - 𝒱̅_♢≤ a^2 + 2b ,

    where Λ̅ and 𝒱̅ are channels that access the block-encodings of Λ and 𝒱 by appending an ancilla qubit and projecting onto |0⟩:
    
    Λ̅(ρ) = (⟨ 0 | ⊗ I) ·Λ( | 0 ⟩⟨ 0 | ⊗ρ) · ( |0 ⟩⊗ I) 
       𝒱̅(ρ) = (⟨ 0 | ⊗ I) ·𝒱( | 0 ⟩⟨ 0 | ⊗ρ) · ( |0 ⟩⊗ I).


For brevity, we defer the proof of this theorem to Appendix <ref>; there, we also showcase an analogous result for arbitrary block-encodings accessed by projectors Π, Π'. Lemma <ref> indicates that by implementing the probabilistic mixture of block-encoding unitaries Λ(σ) = ∑_j=1^m p_j U_j σ U_j^† on an input state σ = |0⟩⟨ 0 | ⊗ρ, and projecting the block-encoding qubit onto |0⟩ (hence the projectors | 0 ⟩⊗ I and ⟨ 0 | ⊗ I), one can reproduce evolution under the target operator S with diamond norm error a^2 + 2b. In short, this result is proven by showing that the channel Λ̅(ρ) is equal to the probabilistic mixture of block-encoded operators ∑_j=1^m p_j R_j ρ R_j^†, to which the mixing lemma applies. Parallel to the usual mixing lemma, this result enables a quadratic suppression of error by selecting operators R_j and an associated probability distribution p_j such that a=ϵ and b = O(ϵ^2).






 §.§ Stochastic QSP

Lemma <ref> very naturally applies to QSP. In this context, the target operation is a matrix function: S = F(A), yet the operators we have access to are QSP polynomials: R_j = P_j(A). A common goal is to simulate evolution under F(A) as ℱ_A(ρ) = F(A)ρ F(A)^†, which encompasses algorithms such as time evolution, linear systems solvers, and ground state preparation, among many others. Traditionally, one achieves this goal with QSP by finding a suitable polynomial approximation to F(x) as |P(x)-F(x)|≤ϵ, such that evolving under this polynomial with QSP as 𝒫_A(ρ) = P(A) ρ P(A)^† suffers error 𝒫_A - ℱ_A _♢≤ O(ϵ). If P(x) is a degree d polynomial, this procedure requires d queries to the block-encoding of A. 


Here we exploit the mixing lemma to approximate evolution under F(A) to the same level of accuracy, but at asymptotically half the number of queries to the block-encoding. We achieve this by designing an ensemble of polynomials that each approximate F(A) as P_j(A) - F(A)≤ O(√(ϵ)), and an associated probability distribution that obeys ∑_j p_j P_j(A) - F(A) ≤ O(ϵ). Then, Lemma <ref> readily implies that the channel Λ_A(ρ) = ∑_j P_j(A) ρ P_j(A)^† suffers error Λ_A - ℱ_A _♢≤ O(ϵ). We also show that implementing Λ_A requires a number of queries to the block-encoding ≈ d/2 + O(1), a cost reduction stemming from the fact that polynomial approximations of smooth functions tend to have degrees that scale as d = O(log(1/ϵ)) (see Sec. <ref>). Intuitively, this scaling implies that a polynomial that achieves error O(√(ϵ)) (e.g., P_j(x)) has a degree that is asymptotically half that of a polynomial that achieves error O(ϵ) (e.g., P(x)). 


Therefore, rather than implement a degree d polynomial, one can instead sample over an ensemble of polynomials of average degree ≈ d/2 + O(1), while retaining the same level of precision. As the corresponding channel is constructed as a probabilistic mixture of QSP sequences, we term this algorithm Stochastic Quantum Signal Processing:

  
Suppose that F(x) is a bounded function F_[-1,1]≤ 1 with a Chebyshev expansion

    F(x) = ∑_n=0^∞ c_n T_n(x)

on the domain x ∈ [-1,1], where for some degree d ≥ 2 the coefficients decay as

    |c_n| ≤  C e^- q n for all n ≥⌈ d/2 ⌉ +1,

for some constants C, q > 0, and this bound is assumed to be tight at n=⌈ d/2 ⌉ +1: |c_⌈ d/2 ⌉ +1| = C e^-q(⌈ d/2 ⌉ +1). Suppose furthermore that a degree-d truncation of F(x), P(x) = ∑_n=0^d c_n T_n(x), achieves an approximation error of ϵ as:

    | F(x) - P(x) | ≤∑_n=d+1^∞ |c_n| =: ϵ,

such that a QSP implementation of the channel 𝒫_A(ρ) = P(A) ρ P(A)^† for some operator A deviates from the target channel ℱ_A(ρ) = F(A) ρ F(A)^† by an error

    𝒫_A - ℱ_A _♢≤ 2ϵ = O(ϵ),

while making d queries to the block-encoding of A. 

Then there exists an ensemble of ⌊ d/2 ⌋ polynomials {P_j(x)}_j=1^⌊ d/2 ⌋ of degree deg(P_j) ≤ d, and an associated probability distribution p_j such that:

    |  P_j(x)  - F (x) | ≤ 2√(Cϵ/1-e^-q)= O(√(ϵ)) 
       | ∑_j=1^⌊ d/2 ⌋ p_j P_j(x)  - F(x) | ≤ϵ = O(ϵ),

while the average degree of these polynomials is

    d_avg   := ∑_j=1^⌊ d/2 ⌋ p_j deg(P_j) 
       ≤⌈d/2⌉(1 - 2 √(ϵ (1-e^-q)C)) + 2 √(ϵ (1-e^-q)C) + 1/1-e^-q
        = d/2(1 - O( √(ϵ)) ) + O(1) 
       ≈ d/2.

Therefore, according to the mixing lemma for block-encodings (Lemma <ref>), the channel Λ_A(ρ) = ∑_j=1^⌊ d/2 ⌋ p_j P_j(A) ρ P_j(A)^† suffers error

    Λ_A - ℱ_A _♢≤(4C/1-e^-q + 2) ϵ = O(ϵ),

while making d_avg≲ d/2 queries to a block-encoding of A in expectation. PW: Why is it OK to bound the error by O(ϵ)? What if the hidden constant is huge? A referee might object that the method is not practical unless the hidden constant is show to be sufficiently small???


    The underlying goal of this theorem is to simulate evolution under F(A) according to the channel ℱ_A(ρ) = F(A) ρ F(A)^†. For notational convenience, let us define the degree t polynomial truncation of the Chebyshev series of F(x) as
    
    P^[t](x) := ∑_n=0^t c_n T_n(x).

    By the assumption that |c_n| ≤ C e^-qn for n ≥⌈ d/2 ⌉ +1, we find that for t ≥⌈ d/2 ⌉ +1, this polynomial truncation suffers error
    
    | P^[t](x) - F(x) | = |∑_n=t+1^∞ c_n T_n(x) | 
       ≤∑_n=t+1^∞ |c_n| ≤∑_n=t+1^∞ C e^-qn = Ce^-qt/1-e^-q .

    Therefore, the error suffered by the degree d truncation P^[d](x) = P(x) is at most
    
    |P^[d](x) - F(x) | ≤Ce^-qd/1-e^-q =: ϵ.

    
    To estimate evolution under ℱ_A, one can use QSP to implement the channel 𝒫_A(ρ) = P^[d](A) ρ P^[d](A)^†. Using the diamond norm bound of Eq. (<ref>), we see that 𝒫_A suffers error 
    
    𝒫_A - ℱ_A _♢   ≤ 2  P^[d](A) - F(A) 
       ≤ 2 max_x∈[-1,1]| P^[d](x) - F(x) | 
       ≤ 2 ϵ.

    As P^[d](A) is a degree d polynomial, the cost of implementing 𝒫_A is d queries to the block-encoding of A. 
    

    Next, let us define our ensemble of polynomials as
    
    P_j(x) = 
            
                P^[ ⌈ d/2 ⌉ ](x) + c_⌈ d/2 ⌉+j/p_j T_⌈ d/2 ⌉ +j(x)     |c_⌈ d/2 ⌉ + j| > 0 
     
                0     c_⌈ d/2 ⌉ +j =0

    for j=1,2,..., ⌊ d/2 ⌋, where p_j is the associated probability distribution defined as
    
    p_j = |c_⌈ d/2 ⌉ +j|/∑_k=1^⌊ d/2 ⌋ |c_⌈ d/2 ⌉ +k|.

    Intuitively, each polynomial P_j(x) consists of the degree ⌈ d/2 ⌉ truncation P^[ ⌈ d/2 ⌉ ](x), and an additional higher order Chebyshev polynomial chosen such that the average polynomial is the degree d truncation: 
    
    ∑_j=1^⌊ d/2 ⌋ p_j P_j(x) = P^[d](x),
 
    where we've used the identity ⌈ d/2 ⌉ + ⌊ d/2 ⌋ = d. The distribution p_j is chosen such that terms with larger Chebyshev coefficients are given more probability mass and preferentially sampled. Each polynomial in this ensemble is guaranteed to suffer error
    
    |P_j(x) - F(x)|    ≤ |P^[ ⌈ d/2 ⌉ ](x) - F(x) | + |c_⌈ d/2 ⌉ +j|/p_j
       ≤∑_n= ⌈ d/2 ⌉ +1^∞ |c_n| + ∑_k=1^⌊ d/2 ⌋ |c_⌈ d/2 ⌉ +k|
       ≤ 2 ∑_n=⌈ d/2 ⌉ +1^∞ |c_n|  ≤  2 ∑_n=⌈ d/2 ⌉ +1^∞ C e^-qn
        = 2Ce^-q ⌈ d/2 ⌉/1-e^-q≤2Ce^-q d/2 /1-e^-q
        = 2√(C ϵ/1-e^-q) = O(√(ϵ)).

    On the other hand, the average polynomial suffers error
    
    | ∑_j=1^⌊ d/2 ⌋ p_j P_j(x) - F(x) | = | P^[d](x) - F(x) | ≤ϵ.

    
    In the language of the mixing lemma for block-encodings (Lemma <ref>), this corresponds to values a = 2√(Cϵ/1-e^-q) and b = ϵ. Accordingly, the channel Λ_A(ρ) = ∑_j=1^⌊ d/2 ⌋ p_j P_j(A) ρ P_j(A)^† suffers error
    
    Λ_A - ℱ_A _♢≤ a^2 + 2b = (4C/1-e^-q + 2 ) ϵ = O(ϵ) .


    Lastly, the expected cost of instantiating Λ_A(ρ) is the average degree d_avg = ∑_j=1^⌊ d/2 ⌋ p_j (P_j), which corresponds to the average number of queries to the block-encoding. Note that the degrees of the polynomials are (P_j) = ⌈d/2⌉ +j ≤ d. To evaluate the average degree, recall that this theorem assumes that 
    
    |c_n| ≤  C e^- q n for all n ≥⌈ d/2 ⌉ +1,

    where this bound is tight at n= ⌈ d/2 ⌉ + 1: |c_⌈ d/2 ⌉ +1| = C e^-q(⌈ d/2 ⌉ +1). This implies that the mean of the distribution p_j = |c_⌈ d/2 ⌉ +j|/∑_k=1^⌊ d/2 ⌋ |c_⌈ d/2 ⌉ +k| is upper bounded by the mean of the geometric distribution p̃_j = e^-q(⌈ d/2 ⌉ +j)/∑_k=1^⌊ d/2 ⌋ e^-q( ⌈ d/2 ⌉ +k) = e^-qj/∑_k=1^⌊ d/2 ⌋ e^-q k. Hence, we may upper bound d_avg as  
    
    d_avg   = ∑_j=1^⌊ d/2 ⌋ p_j (P_j) = ⌈d/2⌉ + ∑_j=1^⌊ d/2 ⌋ j p_j 
       ≤⌈d/2⌉ + ∑_j=1^⌊ d/2 ⌋ j p̃_j 
        = ⌈d/2⌉ - ⌊d/2⌋e^-q ⌊ d/2 ⌋/1-e^-q ⌊ d/2 ⌋ + 1/1-e^-q
       ≤⌈d/2⌉ - ⌊d/2⌋ 2e^-q⌊ d/2 ⌋ + 1/1-e^-q
       ≤⌈d/2⌉ - ( ⌈d/2⌉-1 ) 2e^-q d/2  + 1/1-e^-q
        = ⌈d/2⌉(1 - 2 √(ϵ (1-e^-q)C)) + 2 √(ϵ (1-e^-q)C) + 1/1-e^-q
        = d/2(1 - O( √(ϵ)) ) + O(1) ,

    where line 3 follows from evaluating mean of the geometric distribution (i.e., evaluating a geometric series), the inequality of line 4 holds for e^-qd/2 = √(ϵ)((1-e^-q)/C )^1/2≤ 1/2 (i.e., sufficiently small ϵ), line 5 follows from ⌊ d/2 ⌋ = ⌈ d/2 ⌉ -1 and e^-q⌊ d/2 ⌋≥ e^-qd/2, and line 6 follows from Eq. (<ref>).
    







  
Suppose that F(x) is a bounded function F_[-1,1]≤ 1 with a Chebyshev expansion

    F(x) = ∑_n=0^∞ c_n T_n(x)

on the domain x ∈ [-1,1], where for some degree d ≥ 2 the coefficients decay as

    |c_n| ≤  C e^- q n for all n ≥  d/2,

for some constants C, q > 0. Suppose furthermore that a degree-d truncation of F(x), P(x) = ∑_n=0^d c_n T_n(x), achieves an approximation error of ϵ as:

    | F(x) - P(x) | ≤∑_n=d+1^∞ |c_n| =: ϵ,

such that a QSP implementation of the channel 𝒫_A(ρ) = P(A) ρ P(A)^† for some operator A deviates from the target channel ℱ_A(ρ) = F(A) ρ F(A)^† by an error

    𝒫_A - ℱ_A _♢≤ 2ϵ = O(ϵ),

while making d queries to the block-encoding of A. 

Then there exists an ensemble of polynomials {P_j(x)} of degree deg(P_j) ≤ d, and an associated probability distribution p_j such that:

    |  P_j(x)  - F (x) | ≤ 2√(ϵ)= O(√(ϵ)) 
       | ∑_j=1 p_j P_j(x)  - F(x) | ≤ϵ = O(ϵ),

while the average degree of these polynomials is

    d_avg   := ∑_j=1 p_j deg(P_j) 
       ≤d/2 + log(C)/2q - log(1-e^-q)/2q + 1/2 + 1/1-e^-q
        = d/2 + O(1).

Therefore, according to the mixing lemma for block-encodings (Lemma <ref>), the channel Λ_A(ρ) = ∑_j=1 p_j P_j(A) ρ P_j(A)^† suffers error

    Λ_A - ℱ_A _♢≤ 6 ϵ = O(ϵ),

while making d_avg≤ d/2 + O(1) queries to a block-encoding of A in expectation. The cost reduction realized by this channel is [Here, ≲ neglects the terms independent of d and C in Eq. (<ref>), which are less relevant than log(C)/2q in practice. See examples in Sec. <ref>.] 

    d_avg/d≲1/2(1+ log(C)/qd) ,

which approaches 1/2 in the limit of large d. 



    The underlying goal of this theorem is to simulate evolution under F(A) according to the channel ℱ_A(ρ) = F(A) ρ F(A)^†. For notational convenience, let us define the degree t polynomial truncation of the Chebyshev series of F(x) as
    
    P^[t](x) := ∑_n=0^t c_n T_n(x).

    By the assumption that |c_n| ≤ C e^-qn for n ≥ d/2, we find that for t ≥ d/2, this polynomial truncation suffers error an over x ∈ [-1,1]:
    
    | P^[t](x) - F(x) | = |∑_n=t+1^∞ c_n T_n(x) | 
       ≤∑_n=t+1^∞ |c_n| ≤∑_n=t+1^∞ C e^-qn = Ce^-qt/1-e^-q .

    Therefore, the error suffered by the degree d truncation P^[d](x) = P(x) is at most
    
    |P^[d](x) - F(x) | ≤Ce^-qd/1-e^-q =: ϵ.

    
    To estimate evolution under ℱ_A, one can use QSP to implement the channel 𝒫_A(ρ) = P^[d](A) ρ P^[d](A)^†. Using the diamond norm bound of Eq. (<ref>), we see that 𝒫_A suffers error 
    
    𝒫_A - ℱ_A _♢   ≤ 2  P^[d](A) - F(A) 
       ≤ 2 max_x∈[-1,1]| P^[d](x) - F(x) | 
       ≤ 2 ϵ.

    As P^[d](A) is a degree d polynomial, the cost of implementing 𝒫_A is d queries to the block-encoding of A. 

    Our goal is to reproduce evolution under ℱ_A to diamond norm error O(ϵ) by using a probabilistic mixture of QSP polynomials and invoking the mixing lemma. We construct these polynomials by truncating the original series at a cutoff degree d^*, selected as follows. Because the mixing lemma quadratically suppresses error, and the degree scales as log(1/ϵ), there should exist an ensemble of polynomials of degree around d^* ≈ d/2 with the same error as the original polynomial. We can then readily determine d^* by demanding that it be the smallest integer suffering error at most √(ϵ):
    
    ( Ce^-qd^*/1-e^-q)^2 ≤ϵ = Ce^-qd/1-e^-q  ⇒
        d^* = ⌈d/2 + log(C)/2q - log(1-e^-q)/2q⌉
           = d/2 + O(1).

    Because the error suffered by a single polynomial obtained by truncating at d^* is greater than that at d, we have d^* < d. 

    Next, let us define our ensemble of polynomials as
    
    P_j(x) = 
            
                P^[ d^* ](x) + c_d^*+j/p_j T_d^* +j(x)     |c_d^* + j| > 0 
     
                0     c_d^* +j =0

    for j=1,2,..., d-d^*, where p_j is the associated probability distribution defined as
    
    p_j = |c_ d^* +j|/∑_k=1^d-d^*  |c_d^* +k|.

    Intuitively, each polynomial P_j(x) consists of the degree d^* truncation P^[ d^* ](x), and an additional higher order Chebyshev polynomial chosen such that the average polynomial is the degree d truncation: 
    
    ∑_j=1^ d-d^*  p_j P_j(x) = P^[d](x).
 
    The distribution p_j is chosen such that terms with larger Chebyshev coefficients are given more probability mass and are preferentially sampled. Each polynomial in this ensemble is guaranteed to suffer error
    
    |P_j(x) - F(x)|    ≤ |P^[ d^* ](x) - F(x) | + |c_ d^* +j|/p_j
       ≤∑_n= d^* +1^∞ |c_n| + ∑_k=1^d-d^* |c_ d^* +k|
       ≤ 2 ∑_n=d^* +1^∞ |c_n|  ≤  2 ∑_n=d^* +1^∞ C e^-qn
        = 2Ce^-q d^* /1-e^-q≤ 2√(ϵ) = O(√(ϵ)).

    This bound also implies that  P_j(x) is bounded as P_j _[-1,1]≤ 1+2√(ϵ). On the other hand, the average polynomial suffers error
    
    | ∑_j=1^ d-d^* p_j P_j(x) - F(x) | = | P^[d](x) - F(x) | ≤ϵ.


    In the language of the mixing lemma for block-encodings (Lemma <ref>), this corresponds to values a = 2√(ϵ) and b = ϵ. Accordingly, the channel Λ_A(ρ) = ∑_j=1^d-d^* p_j P_j(A) ρ P_j(A)^† suffers error
    
    Λ_A - ℱ_A _♢≤ a^2 + 2b = 6 ϵ = O(ϵ) .


    Lastly, the expected cost of instantiating Λ_A(ρ) is the average degree d_avg = ∑_j=1^ d-d^*  p_j (P_j), which corresponds to the average number of queries to the block-encoding. Note that the degrees of these polynomials are (P_j) = d^* +j ≤ d. To evaluate the average degree, recall that this theorem assumes that 
    
    |c_n| ≤  C e^- q n for all n ≥ d/2 ,

    This implies that the mean of the distribution p_j := |c_d^* +j|/∑_k=1^d-d^*  |c_d^* +k| is upper bounded by the mean of the geometric distribution p̃_j := e^-q(d^* +j)/∑_k=1^d-d^*  e^-q( d^* +k) = e^-qj/∑_k=1^d-d^* e^-q k. Hence, we may upper bound d_avg as  
    
    d_avg   = ∑_j=1^ d-d^* p_j (P_j) = d^* + ∑_j=1^ d-d^* j p_j 
       ≤ d^* + ∑_j=1^ d-d^* j p̃_j 
        = d^* - (d-d^*) e^-q (d-d^*)/1-e^-q (d-d^*)  + 1/1-e^-q
       ≤ d^* + 1/1-e^-q
        = ⌈d/2 + log(C)/2q - log(1-e^-q)/2q⌉ + 1/1-e^-q
       ≤d/2 + log(C)/2q - log(1-e^-q)/2q + 1/2 + 1/1-e^-q
        = d/2 + O(1)

    where line 3 follows from evaluating mean of the geometric distribution (i.e., evaluating a geometric series), and line 4 from d^* < d. Accordingly, the cost reduction realized by Λ_A is the ratio 
    
    d_avg/d   ≤1/2(1+ log(C)/qd - log(1-e^-q)/qd + 1/d + 2/(1-e^-q) d ) 
       ≈1/2(1+ log(C)/qd)

    where the last line follows from the fact that the log(C) contribution dominates in practice (see Sec. <ref>). In the large d limit, the cost reduction approaches 1/2 inverse-polynomially fast.  
    
    Lastly, as this construction makes no reference to the eigenvalues or singular values of A, stochastic QSP applies equally as well to QET and QSVT, where P(A) acts on the eigenvaleus or singular vectors of A, respectively. In addition, while the presentation of stochastic QSP here is tailored toward functions expressed in the basis of Chebyshev polynomials, we extend this result to more general functions and arbitrary arbitrary bases in Sec. <ref>.
    




Let us take a minute to interpret this result. According to Theorem <ref>, stochastic QSP replaces a deterministic polynomial P(x) with an ensemble of polynomials { P_j(x) } and probability distribution p_j, whose average evolution achieves the same precision up to a constant factor, but at asymptotically half the cost. Importantly, this result is agnostic to the specific polynomial, requiring only that its coefficients decay exponentially according to Eq. (<ref>). As we showed in Sec. <ref>, this condition is generally satisfied by polynomial approximations to smooth functions, rendering stochastic QSP applicable to a wide range of algorithms. In practice the values of C and q in Eq. (<ref>) can be chosen to minimize the ratio d_avg/d, which effectively means minimizing log(C)/q. For visual intuition on this exponential decay, we provide an illustration of our stochastic QSP construction in Fig. <ref>. 







The channel implemented by stochastic QSP is the probabilistic mixture of polynomials Λ_A(ρ) = ∑_j p_i P_j(A) ρ P_j(A). As we discussed in Sec. <ref>, this channel may be realized by implementing an identical probabilistic mixture of unitaries { U_j } that block-encode the polynomials { P_j(A) }, and post-selecting on successfully accessing these block-encodings. In practice, this can be achieved by independently sampling j ∼ p_j, and implementing the QSP sequence that block-encodes the polynomial P_j(A). Because P_j(x) is only bounded as P_j _[-1,1]≤ 1+2√(ϵ) according to Eq. (<ref>), this implementation may require rescaling the polynomials by 1+2√(ϵ), which incurs a measurement overhead ∼ (1+2√(ϵ))^2 that asymptotically approaches 1. This procedure of course requires knowledge of the QSP phases for each polynomial P_j(A). These phases can be determined classically using an efficient phase finding algorithm, such as those of Refs. <cit.>, thus amounting to a classical pre-computation step. 



Moreover, observe from Eq. (<ref>) that each polynomial in the ensemble is constructed as the degree d^* ≈ d/2 truncation plus a higher order term in the Chebyshev expansion, up to degree d. Inclusion of the degree d^* truncation guarantees that the first condition of the mixing lemma is satisfied with error O(√(ϵ)), and sampling the higher order terms allows the second condition of the mixing lemma to be satisfied with error O(ϵ). The specific sampling distribution is chosen to be proportional to the coefficients of the Chebyshev expansion. Because these are assumed to decay exponentially, the corresponding probability mass is concentrated around ≈ d/2, and so too is the average degree. Through this view, the ensemble of stochastic QSP can be seen as a fixed low order term plus higher order correction terms, which are randomly sampled according to the QDrift protocol <cit.> to leverage randomized compiling. An ensemble with a similar structure, albeit in the context of Trotterization, was recently used in Ref. <cit.> to double the order Trotter formulae through randomized compiling. 


Lastly, it is important to note that while stochastic QSP reduces the expected cost to d_avg≈ d/2 +O(1), it does not reduce the maximum degree of the polynomials implemented: some polynomials in the ensemble will have degree greater than d_avg, with one even having degree d. This however is unavoidable. In fact it is necessary for the ensemble to contain polynomials of degree > d/2 in order to attain a level of error equivalent to a degree-d polynomial. To see this, observe that if all the polynomials in the ensemble had degree at most k < d, then their average ∑_j p_j P_j(x) would also be a degree k<d polynomial. This average polynomial would not be able to achieve a precision equivalent to that of a degree d polynomial, thus failing to meet the second condition of the mixing lemma. Note however that the distribution p_j concentrates around small values of j, meaning that degrees much larger than ≈ d/2 are rare.

Through this interpretation, stochastic QSP is similar to the “semi-quantum matrix processing" algorithm of Ref. <cit.> for estimating expectation values and matrix elements of a matrix-valued function. The authors achieve this by decomposing a degree-d polynomial approximation of this function into the Chebyshev basis, and sampling its constituent Chebsyhev polynomials, and measuring an estimator of the sampled polynomial that converges to the desired expectation value/matrix element. Like stochastic QSP, this procedure can reduce the expected cost of certain algorithms to a value < d if the Chebyshev coefficients decay quickly, but it does not reduce the maximal degree because the order d Chebyshev polynomial can still be sampled. However, this method differ from stochastic QSP in two key ways. First, Ref. <cit.> employs a quasi-probability technique where improvements in average degree are traded for additional variance in the measurements, whereas stochastic QSP always implements a normalized probability distribution. Second, the constructions in Ref. <cit.> explicitly consider a target observable O and modify their circuit to a accommodate a parity measurement Z ⊗ O with an ancilla register. Hence, their construction does not constitute an approximate realization of a desired matrix function, whereas stochastic QSP indeed approximates evolution under said matrix function, thus rendering our approach more modular.






§ GENERALIZATIONS

In the previous section, we introduced stochastic QSP tailored specifically to polynomial approximations obtained from truncated Chebyshev expansions. However, it turns out that this specialization is not necessary, and that stochastic QSP is readily generalized to a much broader range of use cases. In this section, we show how our method applies to both definite- and indefinite-parity polynomials, Taylor series, trigonometric polynomials, generalized QSP <cit.>, and implementations of positive operator-valued measures (POVMs) using QSP.






 §.§ Definite and Indefinite Parity


The statement of Theorem <ref> applies to functions F(x) of indefinite parity[Recall that a function is said to have definite parity if it is either even or odd. Otherwise it has indefinite parity.] and produces an ensemble of polynomials { P_j } that also have indefinite parity. However, if F(x) has definite parity, it turns out that this construction preserves the parity.

The parity of the implemented polynomials P_j is important because definite-parity polynomials admit simpler implementations via QSP. Indeed, by construction QSP can only produce polynomials of definite parity (see Eqs. (<ref>) and (<ref>)). In general, indefinite parity polynomials require using a linear combination of unitaries circuit, which demands an extra ancilla qubit and rescales the resulting block-encoding by 1/2. This rescaling can be undesirable in algorithm construction because it may necessitate amplitude amplification to be corrected, and can be avoided when the target function F(x) has definite parity in the first place. We would like to retain this optimized performance in the context of stochastic QSP, which we can show is indeed the case.

 
Consider the setting of Theorem <ref>, but also suppose that F(x) has definite party. Then there exists an ensemble of polynomials {P_j} with the same parity that satisfy the conditions of the theorem.

 If F(x) = ∑_n=0^∞ c_n T_n(x) has even (odd) parity then all c_n for odd (even) n must vanish. In other words, the function is only supported on even (odd) Chebyshev polynomials. Observe that all P_j in the construction are supported only on subsets of the support of F(x). Hence they must also be even (odd).


Note that this corollary also extends to complex-valued functions, which are usually approximated by QSP by decomposition into their real and imaginary components. Therefore, for an arbitrary target function, realizing stochastic QSP requires a circuit no more complicated than a QSP circuit that approximates the target function. 





 §.§ Taylor Series


As we discussed in Sec. <ref>, Chebyshev expansions of a smooth C^∞ functions admit exponentially-decaying coefficients, and thus yield polynomial approximations that meet the requisite conditions for stochastic QSP. Functions of interest in the quantum algorithms literature commonly exhibit this smoothness property (like cos(x) or exp(-β x)), or are well-approximated by functions that do (like how the step function is approximated by erf(kx)). 

As such, we often desire a closed-form expression for the coefficients in the Chebyshev expansion, allowing us to give concrete guarantees for the values of C and q required by Theorem <ref>. However, for certain functions like √(x) and -xln(x), which are only smooth in certain domains, obtaining a closed form expression for the Chebyshev coefficients is cumbersome, and the literature generally works with Taylor series instead (see for example Theorem 68 of <cit.> and its applications). Fortunately, stochastic QSP directly generalizes to Taylor series, and expansions into bases of bounded polynomials more generally. 

  Suppose { B_n (x) } are a collection of basis functions of degree n respectively, which are all bounded as B_n_[-1,1]≤ 1. Then the statement of Theorem <ref> holds with B_n(x) in place of T_n(x).

 This follows from the fact that the only property of the Chebyshev polynomials T_n(x) leveraged in the proof of Theorem <ref> is that they were bounded as T_n_[-1,1] = 1.



Taylor series methods derive their accuracy from the analysis in Corollary 66 of Ref. <cit.>. The basic idea is the following. Suppose F(x) is analytic and bounded on the interval [-1,1], and our goal is to approximate it by a polynomial on the interval [-1+δ, 1-δ] for small δ < 1. If F(x) = ∑_n=0^∞ c_n x^n is the Taylor series of F(x) with coefficients |c_n| ≤ 1, then the error from truncating at degree d is:

    sup_|x| ≤ 1-δ|∑_n=d+1^∞ c_n x^n | ≤∑_n=d+1^∞  (1-δ)^n.

We immediately obtain an exponential decay in truncation error. 

This is a slightly weaker condition than the one required for Corollary <ref>; we require an exponential bound on the coefficients themselves rather than on the truncation error. But if we are willing to approximate the stretched function F( (1-δ) x ) on the interval [-1,1] instead, then substituting the Taylor expansion yields

    F((1-δ)x) = ∑_n=0^∞ c_n (1-δ)^n x^n

where the new coefficients c_n (1-δ)^n exhibit the desired exponential decay, thus rendering this Taylor series expansion amenable to parallel QSP.



 §.§ Trigonometric Polynomials

Recent works <cit.> have considered a model of quantum computation in which a constant-size control register is strongly coupled to many qubits with an otherwise local connectivity graph. In such an architecture, controlled time evolution can be implemented through the Trotter approximation, but Hamiltonian simulation via QSP techniques remains out of reach due to the small size of the control register. Nonetheless QSP-like transformations of a Hamiltonian H can be implemented through applying QSP to a controlled time evolution operator.

A controlled time evolution operator is a block-encoding of e^i H t. Applying QSP to this block-encoding generates trigonometric polynomials ∑_n c_n e^i n H t. If we select t = π/H and let our variable of interest be θ=Ht = π H/H, then we can approximate functions F : [-π, π] → [-1,1] using degree-d trigonometric polynomials:

    F(θ)    ≈∑_n c_n e^i n θ.

Hence, by selecting B_n(θ) =e^inθ for the basis functions, we see how Corollary <ref> applies in this setting as well. Any method for constructing trigonometric polynomials can be used as long as the coefficients c_n decay exponentially. Conveniently, we show in Appendix <ref> that for C^∞ functions F(θ), their Fourier series have exponentially decaying coefficients. We see that, due to the relationship between Fourier expansions and Chebyshev expansions, stochastic QSP applies equally well to trigonometric polynomials as to regular polynomials.




 §.§ Generalized QSP


Recently a technique was proposed <cit.> for optimizing QSP implementation, specifically when the block-encoded operator U is unitary and is encoded via the controlled-U operation. In this situation the usual constraints on parity can be lifted when synthesizing complex polynomials, enabling polynomials of indefinite parity to be generated directly through QSP and avoiding rescaling from using LCU circuit. For real polynomials it is possible to remove parity constraints using Theorem 56 of <cit.>, but this introduces an undesirable factor of 1/2 as discussed earlier. Using the methods of <cit.> this factor can be avoided.

By this reasoning, stochastic QSP is compatible with the construction of Ref. <cit.>, which halves the asymptotic cost of QSP-based Hamiltonian simulation by using generalized QSP. This is achieved by designing a block encoding of both the quantum walk operator and its inverse, to which generalized QSP may be applied to approximate e^-iHt at roughly half the cost of ordinary QSP. Stochastic QSP could be applied on top of this algorithm to further reduce the cost, by a asymptotic factor of 1/4. 




stochastic QSP is compatible with this optimization. From Theorem 3 of <cit.>, we see the only additional requirement is that all sampled polynomials P_j satisfy |P_j(x)|^2 ≤ 1 for x on the complex unit circle. If this is the case for F(x), then, since the P_j are all √(ϵ)-close (see Eq. (<ref>)), we can achieve |P_j(x)|^2 ≤ 1 through rescaling the polynomials by 1+2√(ϵ). This rescaling is already presented earlier in footnote 3; can probably forgo this here






 §.§ POVMs


QSP methods typically project out a single polynomial from the QSP sequence, conventionally taken to be the polynomial P(A) encoded in the |0⟩⟨ 0| block (see Eqs. (<ref>) and (<ref>)). Indeed stochastic QSP is concerned with approximating evolution under a channel by projecting out this component. However, in some situations <cit.>, we are also interested in the other elements of the QSP sequence– particularly those that influence the measurements of the ancilla qubit(s) used to construct the block-encoding. In the conventional encoding, if the ancilla qubit is initialized as |0⟩, then by invoking unitarity and drawing analogy to Eq. (<ref>), this measurement implements a POVM with operators |P(A)|^2 and (1-A^2)|Q(A)|^2 = I-|P(A)|^2. In practice this is desired to approximate a target POVM with operators |F(A)|^2 and I- |F(A)|^2. 


We find that stochastic QSP also applies to the approximation of such a target POVM. We can formalize this by considering an ideal quantum map of the following form. Suppose A has a spectral value decomposition ∑_iλ_i |λ_i⟩⟨λ_i|, and that there is some function G(x) satisfying |F(x)|^2 + (1-x^2)|G(x)|^2 = 1. Then our goal is to approximate the following map:

    |0⟩⊗|λ_i⟩→|0⟩⊗ F(λ_i)|λ_i⟩ + |1⟩⊗ G(λ_i) √(1-λ_i^2)|λ_i⟩.


  
In the setting of Theorem <ref>, consider the ensemble of quantum circuits that implement stochastic QSP (i.e. the QSP circuits that implement {P_j(A)}), upon leaving the QSP ancilla qubit(s) unmeasured. Denoting this ensemble by {U_j}, the quantum channel ∑_j p_j U_j ρ U_j^† approximates the map in Eq. (<ref>) to error O(ϵ) in diamond norm.


In the sake of brevity, we defer this proof to Appendix <ref>. We also note that this result naturally extends to QSVT, in which case the desired mapping is analogous to Eq. (<ref>), but acting on the singular values and singular vectors. 


    




§ APPLICATIONS




Here we apply stochastic QSP to several common polynomials used in the quantum algorithms literature to assess its performance in practice. We study four polynomials:


    
  * The Jacobi-Anger expansion of cosine <cit.>:
    
    cos(tx) = J_0(t) + 2∑_n=1^∞ (-1)^n J_2n(t) T_2n(x),

    where the J_2n(t) are the Bessel functions of the first kind. To achieve additive error at most ϵ, this series may be truncated at degree d = O(|t| + log(1/ϵ)/log(e + log(1/ϵ)/|t|)). This polynomial, in conjunction with the analogous expansion for sin(tx), furnishes an algorithm for Hamiltonian simulation with near-optimal query complexity <cit.>. 
    
    
  * The Jacobi-Anger expansion of an exponential decay <cit.>:
    
    e^-β(x+1) = e^-β[ I_0(β) + 2∑_n=1^∞  I_n(β) T_b(-x)],

    where the I_n(β) is the modified Bessel functions of the first kind. This expansion may be truncated at degree d = O(√((β + log(1/ϵ))log(1/ϵ))) = O(√(β)log(1/ϵ)) to achieve error at most ϵ. Naturally, the resulting polynomial is commonly used for imaginary time evolution <cit.>, thermal state preparation <cit.>, and partition function estimation <cit.>. 
    
    
  * A smooth approximation of 1/x in a domain away from the origin <cit.>:
    
    1/x   ≈1- (1-x^2)^b/x
       = 4∑_n=0^b-1 (-1)^n 2^-2b∑_m=n+1^b 2bb+m T_2n+1(x)

    for an even integer b ≫ 1. While this series is a degree O(b) polynomial, its coefficients decay rapidly, such that it can be further truncated to degree d= O( √(b log(b/ϵ))) while guaranteeing error at most ϵ. The resulting polynomial is particularly useful for inverting matrices and thus solving linear systems of equations. If we take the (non-zero) eigenvalues of the matrix to be lower-bounded as |λ| ≥λ_min, then in order to ensure that the polynomial approximation behaves as ≈ 1/x over the range of eigenvalues, it suffices to choose chooses b = O ( 1/λ_min^2log(1/λ_minϵ) ). This corresponds to a truncation degree d = O(κlog(κ/ϵ)), where κ := 1/λ_min is the condition number. For completeness, we note that algorithms with improved performance have very recently been discovered <cit.>.
    
    
    
    
  * An approximation of (kx) obtained from integrating the Jacobi-Anger expansion of a Gaussian <cit.>:
    
    (kx) = 2ke^-k^2/2/√(π)[I_0(k^2/2)x 
           + ∑_n=1^∞ I_0(k^2/2) (-1)^n (  T_2n+1(x)/2n+1 - T_2n-1/2n-1) ].

    To achieve error ϵ, it suffices to truncate this series at degree d = O(√((k^2 + log(1/ϵ))log(1/ϵ))) = O(k log(1/ϵ)). In practice, this polynomial is used to approximate the step function by selecting k ≫ 1. Notable applications of this approximation include the search problem  <cit.>, phase estimation <cit.> and ground state preparation <cit.>.




All four of these functions feature a cost parameter, namely t,β,b, and k respectively, whose value determines the truncation degree necessary to achieve an accurate approximation.

We apply stochastic QSP to these polynomials, and illustrate the cost reduction ratio d_avg/d as a function of d in Figure <ref>. We rely on the following procedure to compute d_avg. First, we select an integer n_1, and then numerically determine values of C,q such that c_n ≤ Ce^-qn holds for n ≥ n_1. This is achieved by selecting another integer n_2 > n_1 and computing C,q such that Ce^-qn goes through the points (n_1, c_n_1) and (n_2, c_n_2). In doing so, we choose n_1,n_2 to guarantee c_n ≤ Ce^-qn indeed holds for all n ≥ n_1, and also to minimize log(C)/q so as to heuristically reduce the dominant term in the bound on d_avg in Eq. (<ref>). We also select n_1,n_2 independent of the degree d; we find that n_1 naturally converges to the degree at which |c_n| starts to decay exponentially. For cos(tx) and (kx) we find that this regime sets in later as the respective cost parameter increases.
Lastly, for each d, the cutoff degree d^* is computed from Eq. (<ref>). Then d_avg is obtained by explicitly computing the probabilities p_j from Eq. (<ref>) and calculating the corresponding weighted average of degrees. 

In Figure <ref> we observe the desired phenomenon: the cost reduction ratio d_avg/d approaches 1/2 as d increases, with a discrepancy scaling as O(1/d). Additionally, as the cost parameter increases, the magnitude of the error terms ∼log(C)/qd can also increase, resulting in a later approach to 1/2. This is expected because a larger cost parameter requires a larger degree to maintain the same level of approximation. These results indicate that indeed, stochastic QSP reduces the query complexity of QSP by approximately 1/2 in regimes of practical interest. 

A more surprising phenomenon is that for some functions and values of the cost parameter, d_avg/d approaches 1/2 from below, resulting in improved performance for small d. This arises because we determine d_avg by choosing values C,q so as to minimize the deviation log(C)/q. In some cases, this can result in a value C<1, or equivalently log(C)<0, which causes the ratio d_avg/d to deviate from 1/2 by a negative value as per Eq. (<ref>). Moreover, the sawtooth behavior observed for each function in Figure <ref> is explained by the presence of the ceiling function in the definition of d^* in Eq. (<ref>). When d is small then the effect of this rounding is more pronounced.





§ CONCLUSION

By merging quantum signal processing and randomized compiling, we have developed stochastic QSP. As we showed, stochastic QSP enables us to lift a QSP algorithm to its randomized counterpart, and simultaneously shrink the circuit complexity by a factor of 2. We empirically verified this cost reduction for various algorithms of interest, including real/imaginary time evolution, matrix inversion, phase estimation, and ground state preparation. This reduction can also be interpreted as enabling a cost parameter in the underlying function to increase (e.g. a longer time t in Hamiltonian simulation) without changing the query complexity. In aggregate, this result demonstrates that classical randomness is a useful resource for quantum computing, and can help bring quantum algorithms closer to the near-term. 



Moreover, in this work we did not consider noisy gates, but rather assumed the ability to perform QSP exactly. As such, we leveraged randomized compiling to suppress error in QSP polynomial approximations to functions. Nonetheless, as randomized compiling can also suppress noise in erroneous gates <cit.>, this suggests that a practical implementation of stochastic QSP could benefit from also randomizing over the gate set, as a sort of doubly-stochastic channel. Along these lines, it would be interesting to study the requirements and conditions for implementing stochastic QSP on near-term quantum hardware. 

The performance improvement realized through randomized compiling suggests further uses of this technique in quantum information. While here we have applied randomized compiling to quantum algorithms via QSP, it is likely that randomized compiling can confer a similar advantage to the traditional constructions of quantum algorithms (e.g., Grover search <cit.>, or conventional phase estimation via the quantum Fourier transform <cit.>). Further, it would be interesting to search for problems for which randomized compiling (or a variant thereof) confers a super-quadratic suppression of error, translating to a cost reduction by a factor smaller than 1/2 in our context. With an eye toward applications, the mixing lemma could also be used to better understand and generate random unitaries and unitary designs <cit.>, and perhaps even be integrated with randomized measurements <cit.> to improve near-term protocols for studying quantum systems. Likewise, there is an absence of randomized compiling in classical simulation algorithms, which could admit similar improvements for problems aimed to emulate quantum mechanics. Given the ubiquity of random processes in the physical world, it is only natural that we can harness randomness to gain deeper insights into quantum systems.


Acknowledgements: The authors thank Pawel Wocjan, Anirban Chowdhury, Isaac Chuang, and Zane Rossi for useful discussion and feedback. JMM acknowledges the gracious support of IBM Quantum through the IBM internship program. 



apsrev4-2







§ POLYNOMIAL APPROXIMATION TO SMOOTH FUNCTIONS

Here we recall theorems on polynomial approximations to smooth functions, specifically the decay of Fourier coefficients and Chebyshev coefficients of smooth functions. For more details on these results, see Refs. <cit.>.



 §.§ Decay of Fourier Coefficients

Fourier analysis can be used to determine the rate of decay of the Fourier coefficients of smooth functions. The following is a classic result <cit.>, for which we provide a brief proof:

    Let G(θ) be a function with period 2π and Fourier series 
    
    G(θ) = a_0/2 + ∑_n=-∞^∞ a_n e^inθ.

    If G(θ) is a C^k function on θ∈ [0,2π) (i.e., continuous and differentiable through order k), then the Fourier coefficients decay polynomially as 
    
    a_n = o(1/n^k).

    Similarly, if G(θ) is a C^∞ function on θ∈ [0,2π), then the Fourier coefficients decay super-polynomially as 
    
    a_n = O(e^-qn^r)

    for some q,r > 0. 


    The Fourier coefficients are
    
    a_n = 1/2π∫_0^2π G(θ) e^-i nθ dθ .

    Using integration by parts repeatedly, we can write this as 
    
    a_n    = 1/2π∫_0^2π G(θ) e^-i nθ dθ
       = i/n1/2π∫_0^2π G'(θ) e^-inθ dθ
        =   ... 
       = i^k/n^k1/2π∫_0^2π G^(k)(θ) e^-inθ dθ .

    This implies that the Fourier coefficients of the kth derivative G^(k)(θ) are (-in)^k a_n. 

    According to the Riemann-Lebesgue lemma <cit.>, the Fourier coefficients of a smooth function go to 0 as n→∞. Therefore, if G(θ) is a C^k function, then G^(k)(θ) is continuous, and thus its Fourier coefficients decay as lim_n→∞ (-in)^k a_n = 0. This limit implies that a_n = o(1/n^k). 
    
    On the other hand, if G(θ) is a C^∞ function, then by an analogous argument, lim_n→∞ n^k a_n = 0 for all integers k ≥ 1. This implies that a_n decays super-polynomially, i.e. as a_n = O(e^-qn^r) for some q,r>0. 





 §.§ Decay of Chebyshev Coefficients

An analogous result can be derived for the decay of Chebyshev coefficients. Recall that the nth Chebyshev polynomial T_n(x) is a degree n polynomial defined on x∈ [-1,1] as

    T_n(x) = cos(n arccos(x)).

It can be shown that T_n(x) is a degree n polynomial of definite parity (either even or odd, depending on n) and bounded magnitude |T_n(x)|_[-1,1] = 1 <cit.>. The Chebyshev polynomials provide a convenient basis for expanding functions on x∈ [-1,1]. A function F(x) can be expanded as

    F(x) = c_0/2 + ∑_n=1^∞ c_n T_n(x) ,

where

    c_n = 2/π∫_-1^1 F(x) T_n(x)/√(1-x^2) dx

are the Chebyshev coefficients for all n≥ 0. 

By the relation between Chebyshev series and Fourier series, it can be shown that the Chebyshev coefficients decay in a manner analogous to the Fourier coefficients: 

    Let F(x) be a function on x∈ [-1,1] with the Chebyshev series 
    
    F(x) = c_0/2 + ∑_n=1^∞ c_n T_n(x) ,

    If F(x) is a C^k function on x ∈ [-1,1], then the Chebyshev coefficients decay polynomially as 
    
    c_n = o(1/n^k).

    Similarly, if F(x) is a C^∞ function on x ∈ [-1,1], then the Chebyshev coefficients decay super-polynomially as 
    
    c_n = O(e^-qn^r)

    for some q,r > 0. 


    Let θ = arccos(x), or equivalently x = cosθ. Using this change of variables, we can re-express the Chebyshev coefficients as
    
    c_n    = 2/π∫_-1^1 F(x) T_n(x)/√(1-x^2) dx = 2/π∫_0^π F(cos(θ)) cos(n θ) dθ
       = 1/2π∫_0^2π F(cos(θ)) e^-i n θ dθ

    where the last equality stems from F(cos(θ)) being an even function of θ. Comparing with Eq. (<ref>), we see that the Chebyshev coefficients c_n are equal to the Fourier coefficients of the function G(θ) := F(cos(θ)). Therefore, if we can show that F(x) being a C^k function of x implies that G(θ) = F(cos(θ)) is a C^k function of θ, then we can inherit the result of Theorem <ref> to prove the purported decay pattern of the coefficients c_n. 

    This is easy to show. Observe that the first derivative of G(θ) is
    
    G'(θ) = F'(cosθ) (-sinθ) = - F'(x) sinθ .

    If F(x) is a C^1 function of x=cosθ, such that F'(x) is a continuous function of x, then G'(θ) = -F'(cosθ) sin (θ) is a composition of continuous functions of θ. Therefore, G'(θ) is is also a continuous function of θ, which implies that G(θ) is a C^1 function of θ. A similar trend persists at higher orders: by using the chain rule and product rule, the kth derivative G^(k)(θ) can be expressed as a linear combination of products of derivatives F^(j)(cosθ) for j≤ k and trigonometric polynomials sin^a(θ) cos^b(θ) for some integers a,b≥ 0. Thus, if F(x) is a C^k function of x, then this expression for g^(k)(θ) is also continuous, which implies that G(θ) is a C^k function. 

    Therefore, by appealing to Theorem <ref>, we see that if F(x) is a C^k function, then the Chebyshev coefficients c_n decays polynomially as c_n = o(1/n^k). Similarly, if F(x) is a C^∞ function, then c_n decays super-polynomially as c_n = O(e^-qn^r) for some q,r > 0.  







§ EXTENSIONS OF THE MIXING LEMMA

Here we present generalizations of the mixing lemma to non-unitary evolution and block-encodings. The proofs of these variants parallel the proof of the mixing lemma presented in Ref. <cit.>. 



 §.§ The Generalized Mixing Lemma


Consider a mixing lemma for arbitrary operators, including non-unitary evolution. In this scenario, we we wish to implement a channel 𝒮(ρ) = S ρ S^†, where S is not necessarily unitary, yet we only have access to operators R_j that approximate S. Then we can prove the following:

    Let S be a target operator, possibly non-unitary, and 𝒮(ρ) = S ρ S^† the corresponding channel. Suppose there exist operators R_j and a probability distribution p_j that approximate S as
    
    R_j - S ≤ a  for all  j,
       ∑_j p_j R_j - S ≤ b ,

    for some a,b > 0. Then, the corresponding channel Λ(ρ) = ∑_j p_j R_j ρ R_j^† approximates the channel 𝒮 as
    
    Λ - 𝒮_♢≤ a^2 + 2b S .



Paralleling the proof of the mixing lemma in Ref. <cit.>, let δ_j = R_j - R. This obeys δ_j ≤ a and ∑_j p_j δ_j ≤ b by the assumed conditions of the theorem. The action of the channel 𝒮 can then be expanded as

    Λ(ρ)    = ∑_j p_j R_j ρ R_j^†
       = ∑_j p_j (S+δ_j) ρ (S+δ_j)^†
       = Sρ S^† + (∑_j p_j δ_j )ρ S^†
          + S ρ(∑_j p_j δ_j^†) + ∑_j δ_j ρδ_j^†.

We can then bound the 1-norm Λ(ρ) - 𝒮(ρ) _1 via the triangle inequality as 

    Λ(ρ) - 𝒮(ρ) _1 ≤    (∑_j p_j δ_j )ρ S^†_1 
        +  S ρ(∑_j p_j δ_j^†) _1 + ∑_j δ_j ρδ_j^†_1 .


The first term on the right hand side of Eq. (<ref>) can be bounded by appealing to Holder's inequality (specifically, AB ≤AB_1 and AB ≤A_1 B in this context):

    ( ∑_j p_j δ_j ) ρ S^†_1    ≤∑_j p_j δ_j ρ_1  S^†
       ≤∑_j p_j δ_j ρ_1  S 
       ≤ S  b,

where we have used that ρ_1 = 1. Analogously, the second term in Eq. (<ref>) can be upper bounded by S  b. Again invoking Holder's inequality, the last term in Eq. (<ref>) can be bounded as 

    ∑_j δ_j ρδ_j^†_1    ≤∑_j p_j δ_j ρδ_j^†_1 
       ≤∑_j p_j δ_j ρ_1 δ_j^†
       ≤∑_j p_j δ_j ρ_1 δ_j^†
       ≤∑_j p_j δ_j δ_j 
       ≤ a^2.

Therefore we have

    Λ(ρ) - 𝒮(ρ) _1 ≤ a^2 + 2b S ,

for all density matrices ρ. This translates to a diamond norm bound:

    Λ - 𝒮_♢   = sup_σ_1 ≤ 1 (Λ⊗ℐ )(σ) - (𝒮⊗ℐ)(σ) _1 
       ≤ a^2 + 2b S  ,

because appending the identity channel ℐ does not alter the proof of 1-norm bound in Eq. (<ref>), which holds true for all density matrices.


Notably, this result is analogous to the original mixing lemma, but modified by the spectral norm of S. For unitary evolution where S  = 1, the generalized mixing lemma reduces to the usual mixing lemma. 










 §.§ The Mixing Lemma for Block-Encodings

We can consider an analog of the mixing lemma for block-encodings. In this scenario, we wish to evolve under an operator S (which may be non-unitary) encoded in a unitary V, yet we only have access to operators R_j block-encoded in unitaries U_j. We first study the case in which the operators are block-encoded in the |0⟩⟨ 0| block of their respective unitaries, which is the simplest form of a block-encoding and most often considered in the literature. Afterwards, we proceed to the general case, in which the operators are encoded through arbitrary projectors Π and Π'. 


First consider encodings in the |0⟩⟨ 0| block:
[Mixing Lemma for Block-Encodings: |0⟩⟨ 0| Block]
    Let V be a unitary that block-encodes a (possibly non-unitary) target operator S as S = (⟨ 0 | ⊗ I) V ( |0 ⟩⊗ I). Suppose there exist unitaries U_j that block-encode operators R_j as R_j = (⟨ 0 | ⊗ I)  U_j ( | 0 ⟩⊗ I). Also suppose there exists a probability distribution p_j such that
    
    R_j  - S ≤ a  for all  j , 
       ∑_j p_j R_j - S ≤ b ,

    for some a,b ≥ 0. Then, the corresponding unitary channel Λ(ρ) = ∑_j p_j U_j ρ U_j^† approximates the action of the channel 𝒱 as
    
    Λ̅ - 𝒱̅_♢≤ a^2 + 2b ,

    where Λ̅ and 𝒱̅ are channels that access the block-encodings of Λ and 𝒱 by appending an ancilla qubit and post-selecting:
    
    Λ̅(ρ) = (⟨ 0 | ⊗ I) ·Λ( | 0 ⟩⟨ 0 | ⊗ρ) · ( |0 ⟩⊗ I) 
       𝒱̅(ρ) = (⟨ 0 | ⊗ I) ·𝒱( | 0 ⟩⟨ 0 | ⊗ρ) · ( |0 ⟩⊗ I).



   
    First, observe that by linearity 
    
    Λ̅(ρ)    =  (⟨ 0 | ⊗ I ) Λ( | 0 ⟩⟨ 0 | ⊗ρ) ( |0 ⟩⊗ I) 
       = ∑_j p_j R_j ρ R_j^†,

    and 
    
    𝒱̅(ρ)    = (⟨ 0 | ⊗ I ) 𝒱( | 0 ⟩⟨ 0 | ⊗ρ) ( |0 ⟩⊗ I) 
        = S ρ S^†.

    Therefore, Λ̅ and 𝒱̅ are analogous to the channels Λ and 𝒮 considered in Theorem <ref>. Likewise, S and R_j obey the requisite conditions of Theorem <ref> as per the assumptions of this theorem. Accordingly, Theorem <ref> implies that
    
    Λ̅ - 𝒱̅_♢≤ a^2 + 2b  S ≤ a^2 + 2b,

    where we have used that S≤ 1 because S is block-encoded in a unitary.




Next, consider the general case in which the operators are block-encoded by arbitrary projectors Π and Π'. Specifically, in this case, an operator is block-encoded as S = Π V Π' for projection operators Π, Π'. In practice, this means that S is accessed by first applying Π' to project the quantum state into the block of interest, then applying V, and finally applying Π' to extract the desired block. With this intuition, we can prove the following:

    Let V be a unitary, that block-encodes a (possibly non-unitary) target operator S as S = Π V Π', where Π, Π' are orthogonal projectors. Suppose there exist unitaries U_j that block-encode operators R_j as R_j = Π U_j Π '. Also suppose there exists a probability distribution p_j such that
    
    R_j  - S ≤ a  for all  j , 
       ∑_j p_j R_j - S ≤ b .

    Then, the corresponding unitary channel Λ(ρ) = ∑_j p_j U_j ρ U_j^† approximates the action of the channel 𝒱 as
    
    Λ̅ - 𝒱̅_♢≤ a^2 + 2b ,

    where Λ̅ and 𝒱̅ are channels that access the block-encodings of Λ and 𝒱 by applying the projectors Π and Π':
    
    Λ̅(ρ) = Π·Λ( Π' ρΠ' ) ·Π
       𝒱̅(ρ) = Π·𝒱( Π' ρΠ' ) ·Π.




    Analogous to the previous proof, observe that by linearity
    
    Λ̅(ρ) = Π·Λ( Π' ρΠ' ) ·Π
        = ∑_j p_j Π U_j Π' ·ρ·Π' U_j^†Π
        = ∑_j p_j R_j ρ U_j^†,

    and 
    
    𝒱̅(ρ)    = Π·𝒱( Π' ρΠ' ) ·Π
        = Π V Π' ·ρ·Π' V^†Π
        = S ρ S^†,

    where we have use the fact the the orthogonal projectors are Hermitian: Π^† = Π,  Π'^† = Π'. Therefore, we again see that Λ̅ and 𝒱̅ are analogous to the channels Λ and 𝒮 considered in Theorem <ref>, which implies that
    
    Λ̅ - 𝒱̅_♢≤ a^2 + 2b  S ≤ a^2 + 2b,

    where we have again used that S≤ 1 because S is block-encoded in a unitary.





§ PROOF OF COROLLARY <REF>

Here we prove Corollary <ref>. Recall that our map of interest is 

    |0⟩⊗|λ_i⟩→|0⟩⊗ F(λ_i)|λ_i⟩ + |1⟩⊗ G(λ_i) √(1-λ_i^2)|λ_i⟩.



In the setting of Theorem <ref>, consider the ensemble of quantum circuits that implement stochastic QSP (i.e. the QSP circuits that implement {P_j(A)}), upon leaving the QSP ancilla qubit(s) unmeasured. Denoting this ensemble by {U_j}, the quantum channel ∑_j p_j U_j ρ U_j^† approximates the map in Eq. (<ref>) to error O(ϵ) in diamond norm.



We consider a function F(x) with definite parity for ease of presentation, but our analysis holds in the general case. When we apply QSP to an operator A with spectral decomposition ∑_iλ_i |λ_i⟩⟨λ_i| to a definite-parity real polynomial P_j, we implement the map

    |0⟩⊗|λ_i⟩→
             |0⟩⊗ P_j(λ_i)|λ_i⟩ + |1⟩⊗ Q_j(λ_i) √(1-λ_i^2)|λ_i⟩,

where Q_j is some polynomial satisfying

    |P_j(x)|^2 + (1-x^2)|Q_j(x)|^2 = 1

for all x ∈ [-1,1]. To show that the ensemble of unitaries implementing these polynomials approximates the desired map, we again rely on Lemma <ref>. The error of any individual map j is bounded by:

    sup_x |   |0⟩⊗( P_j(x) - F(x))  
          + |1⟩⊗( Q_j(x) - G(x))√(1-x^2)| 
    ≤   √( (2√(ϵ))^2  + sup_x |Q_j(x) - G(x) |^2(1-x^2) ),

since we have already established |P_j - F|  ≤ 2√(ϵ) in the proof of Theorem <ref>. Similarly the error of the probabilistic mixture over the maps is similarly bounded by:

    √(ϵ^2  + sup_x |∑_j p_j Q_j(x) - G(x) |^2(1-x^2)).


Our goal is to bound the above two quantities by deriving properties of the Q_j from |P_j(x)|^2 + (1-x^2)|Q_j(x)|^2 = 1. We define for some fixed j:

    α(x) := P_j(x) - F(x)  and β(x) := Q_j(x)  - G(x).

Now we calculate:

    1    =  |F|^2 + (1-x^2) |G|^2
       =  |P_j + α|^2 + (1-x^2) |Q_j + β|^2
       = (|P_j|^2 + α P_j^* + α^* P_j + |α|^2)
       + (1+x^2)(|Q_j|^2 + β Q_j^* + β^* Q_j + |β|^2)
       = 1 + (α P_j^* + α^* P_j + |α|^2)
       + (1+x^2)(β Q_j^* + β^* Q_j + |β|^2),

where ^* denotes the complex conjugate here. 
Thus:

    (1+x^2)|β Q_j^* + β^* Q_j + |β|^2| = |α P_j^* + α^* P_j + |α|^2|.

We have:

    (1+x^2)|β Q_j^* + β^* Q_j + |β|^2| ≥ |β|^2
    
    |α P_j^* + α^* P_j + |α|^2| ≤ 2|α| + |α|^2 ≤ 3|α|

So |β|^2 ≤ 3|α| ≤ 6√(ϵ).  Hence the quantity in Eq. (<ref>) is bounded by √(40ϵ). An identical calculation considering α̅(x) := ∑_j p_j P_j(x) - F(x) and β̅(x) := ∑_j p_j Q_j(x)  - G(x) yields the bound |β̅|^2 ≤ 3|α̅|, hence bounding the quantity in Eq. (<ref>) by √(10)ϵ≤ 4ϵ.

Now we apply Lemma <ref> with a = √(40ϵ) and b = 4ϵ we obtain an error in diamond norm of a^2 + 2b = 48ϵ = O(ϵ). Naturally the constant 48 can be tightened significantly using a more careful analysis.







