# An Improved Classical Singular Value Transformation
for Quantum Machine Learning

Ainesh Bakshi
<EMAIL_ADDRESS>
MIT  Ewin Tang
<EMAIL_ADDRESS>
University of Washington

###### Abstract

Quantum machine learning (QML) has shown great potential to produce large
quantum speedups for computationally intensive linear algebra tasks. The
quantum singular value transformation (QSVT), introduced by Gilyén, Su, Low
and Wiebe [GSLW19], is a unifying framework to obtain QML algorithms. We
provide a classical algorithm that matches the performance of QSVT on low-rank
inputs, up to a small polynomial overhead. Under efficient quantum-accessible
memory assumptions, given a bounded matrix $A\in\mathbb{C}^{m\times n}$, a
vector $b\in\mathbb{C}^{n}$, and a bounded degree-$d$ polynomial $p$, QSVT can
output a measurement from the state $|p(A)b\rangle$ in $\mathcal{O}(d\lVert
A\rVert_{F})$ time after linear-time pre-processing. We show that, in the same
setting, for any $\varepsilon>0$, we can output a vector $v$ such that
$\|v-p(A)b\|\leqslant\varepsilon\|b\|$ in
$\mathcal{O}(d^{9}\MT@delim@RIGHT\Norm@star{A}_{F}^{4}/\varepsilon^{2})$ time
after linear-time pre-processing. This improves upon the best known classical
algorithm [CGL+20a], which requires
$\mathcal{O}(d^{22}\MT@delim@RIGHT\Norm@star{A}_{F}^{6}/\varepsilon^{6})$
time.

Instantiating the aforementioned algorithm with different polynomials, we
obtain fast _quantum-inspired_ algorithms for regression, recommendation
systems, and Hamiltonian simulation. We improve in numerous parameter settings
on prior work, including those that use problem-specialized approaches, for
quantum-inspired regression [CGL+20b, CGL+20a, SM21, GST22, CCH+22] and
quantum-inspired recommendation systems [Tan19, CGL+20a, CCH+22].

Our key insight is to combine the _Clenshaw recurrence_ , an iterative method
for computing matrix polynomials, with sketching techniques to simulate QSVT
classically. The tools we introduce in this work include (a) a _non-oblivious_
matrix sketch for approximately preserving bi-linear forms, (b) a _non-
oblivious_ asymmetric approximate matrix product sketch based on
$\ell_{2}^{2}$ sampling, (c) a new stability analysis for the Clenshaw
recurrence, and (d) a new technique to bound arithmetic progressions of the
coefficients appearing in the Chebyshev series expansion of bounded functions,
each of which may be of independent interest.

###### Contents

  1. 1 Introduction
  2. 2 Our Results
    1. 2.1 Applications: Dequantizing QML
  3. 3 Technical overview
  4. 4 Related work
  5. 5 Preliminaries
    1. 5.1 Linear algebra
    2. 5.2 Polynomials and the Chebyshev Basis
    3. 5.3 Sampling and query access
  6. 6 Extending the Sketching Toolkit
    1. 6.1 The Bi-Linear Entry-wise Sampling Transform
    2. 6.2 Approximate Matrix Product via $\ell^{2}_{2}$ Sampling
  7. 7 Sums of Chebyshev coefficients
  8. 8 Properties of the Clenshaw recursion
    1. 8.1 Deriving the Clenshaw recursions
    2. 8.2 Evaluating even and odd polynomials
  9. 9 Stability of the scalar Clenshaw recursion
    1. 9.1 Analyzing error propagation
    2. 9.2 Bounding the iterates of the Clenshaw recurrence
  10. 10 Computing matrix polynomials
    1. 10.1 Computing odd matrix polynomials
    2. 10.2 Generalizing to even polynomials
    3. 10.3 Bounding iterates of singular value transformation
  11. 11 Dequantizing algorithms
    1. 11.1 Recommendation systems
    2. 11.2 Linear regression
    3. 11.3 Hamiltonian simulation
  12. Acknowledgments

##  1 Introduction

Quantum machine learning (QML) has rapidly developed into a highly active
field of study, with numerous proposals for speeding up machine learning tasks
on quantum computers ([HHL09, RML14, KP17, WKS14, WKS16, BKL+19], and see
[DW20, CHI+18] for a survey). These proposals include quantum algorithms for
several basic tasks in machine learning, including regression, principal
component analysis, support vector machines, recommendation systems,
Hamiltonian simulation and semi-definite programming. The central goal in QML
is to demonstrate a problem on which quantum computers obtain a substantial
practical speedup over classical computers. A successful resolution of this
goal would provide compelling motivation to further invest resources into
developing scalable quantum computers (i.e. be a _killer app_ for quantum
computing).

The quantum singular value transformation (QSVT) framework [LC17, CGJ19,
GSLW19] uses ideas from signal processing to present a unified approach to
design algorithms for quantum linear algebra and, by extension, QML. This
framework is known to capture essentially all linear algebraic QML techniques
[MRTC21] (and, more generally, the most prominent examples of quantum
advantage), so it will be the focus of our investigation in this work.
Informally, the QSVT framework defines a central primitive known as the
_block-encoding_ , and shows that, given a matrix $A\in\mathbb{C}^{m\times n}$
with bounded operator norm as a block-encoding, one can generate a block-
encoding for a degree-$d$ polynomial applied to that matrix, $p(A)$, defined
in the appropriate sense, with only $O(d\log(mn))$ overhead in gate
complexity.111We will generally not concern ourselves with $\log(mn)$ factors,
since quantum algorithms typically count bit complexity where classical
algorithms count word RAM complexity, which muddles any comparison of such
$\log(mn)$ and $\log\frac{1}{\varepsilon}$ factors. Efficient block-encodings
do not exist in general, but they do exist for two broad classes of matrices,
assuming appropriately strong forms of coherent access: matrices with low
sparsity (a typical setting for quantum simulation) and matrices with low
stable rank (a typical setting for quantum machine learning). We treat the
latter case; specifically, we assume that the classical input data is in a
quantum-accessible data structure, which allows for efficient (as in low-
depth) preparation of a block-encoding of $A/\lVert A\rVert_{F}$, where
$\lVert A\rVert_{F}=\sum_{i,j}\MT@delim@RIGHT\abs@star{A_{i,j}}^{2}$ denotes
the Frobenius norm of $A$.

This type of block-encoding is the one commonly used for quantum linear
algebra algorithms on classical data [KLLP19, CdW21], since it works for
arbitrary matrices and vectors, paying only a $\lVert
A\rVert_{F}/\MT@delim@RIGHT\norm@star{A}$ in sub-normalization. The main
consequence of the QSVT framework for QML is that, given a matrix
$A\in\mathbb{C}^{m\times n}$ with bounded operator norm, a vector
$b\in\mathbb{C}^{n}$ of unit norm, and a degree-$d$ polynomial $p$ with
$\MT@delim@RIGHT\abs@star{p(x)}\leqslant 1$ for $x\in[-1,1]$, we can output a
sample222There is a cost of post-selection which we ignore in this exposition.
More formally, we consider $|p(\frac{A}{\lVert A\rVert_{F}})\rangle$ to be an
un-normalized state, and define a sample from such a state $|\psi\rangle$ to
be a measurement of $|\psi\rangle$ in the computational basis with probability
$\MT@delim@RIGHT\norm@star{\psi}^{2}$, and “failure” otherwise. from the state
$|p(\frac{A}{\lVert A\rVert_{F}})b\rangle$ in time
$O\MT@delim@RIGHT\Paren@star{d\log(mn)}$. Since in applications we care about
applying polynomials to the singular values of $A$, which in $\frac{A}{\lVert
A\rVert_{F}}$ lie in the range $[0,\frac{\MT@delim@RIGHT\norm@star{A}}{\lVert
A\rVert_{F}}]$, we need to pay an additional overhead of $\frac{\lVert
A\rVert_{F}}{\MT@delim@RIGHT\norm@star{A}}$ to amplify these singular values
[GSLW19, Theorem 30] to the full region $[0,1]$, making the gate complexity to
produce a state $\varepsilon$-close to $|p(A)b\rangle$,
$\mathcal{O}(d\frac{\lVert
A\rVert_{F}}{\MT@delim@RIGHT\norm@star{A}}\log(mn)\log(1/\varepsilon))$. Note
that, though the quantum algorithm has only logarithmic dependence on
$\varepsilon$, if we wish to determine a property of the output state (say,
whether its norm is large or whether it has high overlap with another state),
this incurs an $\mathcal{O}(1/\varepsilon)$, since distinguishing a state from
one $\varepsilon$-far in trace distance requires $\Omega(1/\varepsilon)$
additional overhead, even when given an oracle efficiently preparing that
state. Therefore, the effective running time for QSVT is
$\Omega\MT@delim@RIGHT\Paren@star{\frac{d\lVert
A\rVert_{F}}{\MT@delim@RIGHT\norm@star{A}\varepsilon}\log(mn)}$.

At first glance, it appears that QSVT obtains an exponential speed-up over
classical algorithms, since it scales logarithmically with input size.
However, as first argued in the work of [Tan19], to obtain a fair comparison,
the classical algorithms should be equipped with a $\ell_{2}^{2}$ sampling
data structure to access the input (see Section 5.3 for the precise model).
Given access to such a data structure, and an accuracy parameter
$\varepsilon>0$, Chia, Gilyén, Li, Lin, Tang and Wang [CGL+20a] showed that
there exists a classical algorithm that outputs a vector $v$ such that
$\MT@delim@RIGHT\norm@star{p(A)b-v}\leqslant\varepsilon$, in
$\mathcal{O}(d^{22}\MT@delim@RIGHT\norm@star{A}_{F}^{6}/\varepsilon^{6})$
time. As a result, they obtain _quantum-inspired_ algorithms for quantum-
inspired versions of several fundamental problems in machine learning, albeit
with a running time that is a large polynomial in the degree of the polynomial
$p$, the Frobenius norm of $A$, and the accuracy parameter $\varepsilon$.

Multiple papers [KLLP19, KP22] have conjectured that the large polynomial
running time (significantly larger than quartic) for classical algorithms may
be inherent, and thus can demonstrate a practical speedup for several problems
in QML. This is borne out in prior work on quantum-inspired algorithms, which
is dominated by the cost of computing a singular value decomposition of a
matrix with
$\Omega((\frac{\MT@delim@RIGHT\norm@star{A}_{F}}{\MT@delim@RIGHT\norm@star{A}\varepsilon})^{2})$
rows and columns, immediately incurring a power-six dependence in
$\varepsilon$ and the Frobenius norm of $A$. Therefore, the central question
we address in this paper is as follows:

> _Is the large polynomial running time for dequantizing QSVT inherent?_

##  2 Our Results

We answer the central question above in the negative for all parameters except
polynomial degree, and show that there are indeed better classical algorithms
that simulate QSVT. Our main result is as follows:

######  Theorem 2.1 (Classical Singular Value Transformation, informal Theorem
10.1).

Suppose we are given $A\in\mathbb{C}^{m\times n}$ with
$\MT@delim@RIGHT\norm@star{A}\leqslant 1$ and $b\in\mathbb{C}^{n}$, and an
accuracy parameter $\varepsilon$. Then, after
$O(\operatorname{nnz}(A)+\operatorname{nnz}(b))$ preprocessing time to create
a data structure333If we are already given $A$ and $b$ in the QRAM data
structures needed to prepare a block-encoding of $A/\lVert A\rVert_{F}$ and a
quantum state of $b/\MT@delim@RIGHT\norm@star{b}$, this preprocessing can be
done in $O(d^{8}\log^{8}(d)\log^{2}(n)\lVert A\rVert_{F}^{4}/\varepsilon^{4})$
time., for an even or odd degree-$d$ polynomial $p$ such that
$\MT@delim@RIGHT\abs@star{p(x)}\leqslant 1$ for $x\in[-1,1]$, we can output a
description of a vector $y\in\mathbbm{C}^{n}$ such that with probability at
least $0.9$, $\|y-p(A)b\|\leqslant\varepsilon\MT@delim@RIGHT\norm@star{b}$ in
$O\MT@delim@RIGHT\Paren@star{d^{9}\log^{4}(d)\|A\|_{F}^{4}\log(n)/\varepsilon^{2}}$
time. This description gives $y$ as either $Ax$ or $A^{\dagger}x$ for a sparse
vector $x$, depending on the parity of $p$, and allows us to compute entries
of $y$ in
$\tilde{O}\MT@delim@RIGHT\Paren@star{d^{4}\MT@delim@RIGHT\norm@star{A}_{F}^{2}/\varepsilon^{2}}$
time or obtain an $\ell_{2}^{2}$ sample from $y$ in
$\tilde{O}\MT@delim@RIGHT\Paren@star{d^{6}\MT@delim@RIGHT\norm@star{A}_{F}^{4}/(\varepsilon^{2}\MT@delim@RIGHT\norm@star{y}^{2})}$
time.

######  Remark 2.2 (No large speedup for low-degree QSVT circuits).

Recall the setting of QSVT for QML, where we have $A$ with bounded operator
norm in a quantum-accessible data structure, can apply a QSVT circuit, and
wish to learn a simple linear algebraic property of the output state (say, the
inner product $\langle v|p(A)b\rangle$ for a given vector $v$). As discussed,
the quantum gate complexity of this QSVT algorithm444This is not a “lower
bound” statement, merely a statement about the gate complexity of QSVT. Note
that it’s possible to embed quantumly hard problems into even this sort of
QSVT, such that, for very particular linear algebraic properties, there is a
large quantum speedup. For example, one could ask to compute the Forrelation
[AA18] of the output vector with a given vector; this would require
$\Omega(n)$ queries to the input classically but runs in time polynomial in
$\log(n)$. is $\Omega(d\lVert A\rVert_{F}/\varepsilon)$.

Classically, we can use Theorem 2.1 to compute a description for $p(A)b$ in
$O(d^{9}\log^{4}(d)\lVert A\rVert_{F}^{4}\log(n)/\varepsilon^{2})$ time, which
we can then use to compute an entry or estimate an overlap $\langle
v|p(A/\lVert A\rVert_{F})b\rangle$. The runtime is dominated by the cost of
the initial algorithm, and so the gap between quantum and classical is 1-to-9
for the degree $d$, 1-to-4 for $\lVert A\rVert_{F}$ (which we think of as
square root of stable rank $\frac{\lVert
A\rVert_{F}}{\MT@delim@RIGHT\norm@star{A}}$), and 1-to-2 for $\varepsilon$.
So, considering constant degree, the gap is merely quartic, which is the type
of speedup that prior work suggests might just suffice to achieve a quantum
advantage for intermediate-term quantum computers [BMN+21].

######  Remark 2.3 (Comparison to [CGL+20a, JLGS20, GLG22]).

There are three papers that get similar results about “dequantizing the
quantum singular value transformation”. The work of Chia, Gilyén, Li, Lin,
Tang, and Wang [CGL+20a] gives a runtime of $\widetilde{O}(d^{22}\lVert
A\rVert_{F}^{6}/\varepsilon^{6})$ (after renormalizing so that
$\MT@delim@RIGHT\norm@star{A}=1$ instead of $\lVert A\rVert_{F}=1$). We
improve in all parameters over this work: degree of the polynomial $p$,
Frobenius norm of $A$, and accuracy parameter $\varepsilon$.

The work of Jethwani, Le Gall, and Singh [JLGS20] provides two algorithms for
applying $f(A)b$ for Lipschitz-continuous $f$ and $A$ with condition number
$\kappa$. (Standard results in function approximation state that such
functions can be approximated by polynomials of degree $O(L)$, where $L$ is
the Lipschitz constant of $f$ and the tail is either polynomial or exponential
depending on how smooth $f$ is [Tre19].) In particular, they achieve a runtime
of $O(\lVert A\rVert_{F}^{6}\kappa^{20}(d^{2}+\kappa)^{6}/\varepsilon^{6})$ to
apply a degree-$d$ polynomial. Again, we improve in all parameters over this
work: degree of the polynomial $p$, Frobenius norm of $A$, and accuracy
parameter $\varepsilon$. We also do not incur condition number dependence.

Finally, the work of Gharibian and Le Gall [GLG22] considers QSVT when the
input matrices are sparse, which is the relevant regime for quantum chemistry
and other problems in many-body systems. Here, the matrices of interest are
local Hamiltonians, which are row- and column-sparse. Their main contribution
is distinguishing between constant-degree QSVT circuits, which can be
simulated in polynomial time classically, and polynomial-size QSVT circuits
(in the number of qubits), which can solve BQP-complete problems. We deal with
a different setting, where all circuits can be simulated efficiently in
polynomial time.

Next, we describe the implications of Theorem 2.1 to specific problems in QML
and show that we obtain faster _quantum-inspired_ algorithms for several such
problems.

###  2.1 Applications: Dequantizing QML

We begin by stating the de-quantization result we obtain for regression.

######  Corollary 2.4 (De-quantizing Regression, informal Corollary 11.5).

Given an $\ell_{2}^{2}$-sampling oracle for $A\in\mathbb{C}^{m\times n}$ and
$b\in\mathbb{C}^{n}$ such that
$\MT@delim@RIGHT\Norm@star{A},\MT@delim@RIGHT\norm@star{b}\leqslant 1$, a
parameter $0<\sigma<1$, and an accuracy parameter $\varepsilon>0$, we can
output the representation of a vector $y\in\mathbb{C}^{n}$ such that with
probability at least $0.9$,
$\MT@delim@RIGHT\norm@star{y-A^{+}_{\geqslant\sigma}b}\leqslant\varepsilon$,
where $A_{\geqslant\sigma}^{+}$ denotes the function on $A$ that is the
inverse for singular values that are $\geqslant\sigma$ and smoothly thresholds
away all singular values below $\sigma$. This algorithm runs in
$\tilde{O}\MT@delim@RIGHT\Paren@star{\MT@delim@RIGHT\norm@star{A}_{F}^{4}/(\varepsilon^{2}\sigma^{11})}$
time, where $\sigma$ is the chosen threshold for $A$. We can also output a
sample from $y$ in the same running time.

######  Remark 2.5 (Comparison with [CGL+20b, CCH+22, SM21]).

We note that Chia et. al. [CGL+20b] get a running time of
$O\MT@delim@RIGHT\Paren@star{\MT@delim@RIGHT\norm@star{A}_{F}^{6}/(\varepsilon^{6}\sigma^{28})}$.
We improve over this result in all parameters.

Chepurko et. al. [CCH+22] get a running time of $\tilde{O}(\lVert
A\rVert_{F}^{4}\log(d)/(\sigma^{8}\varepsilon^{4}))$, where $\sigma$ is the
minimum singular value of $A$, assuming that we perform regression with some
sizable regularization $\lambda>0$. We improve over their work in
$\varepsilon$ dependence, and match it in $\lVert A\rVert_{F}$ dependence. As
for the $\sigma$ dependence, this achieves error to
$\varepsilon\MT@delim@RIGHT\norm@star{A}\MT@delim@RIGHT\norm@star{b}$ error,
which is worse than the $\varepsilon\MT@delim@RIGHT\norm@star{x^{*}}$ bound
that we actually achieve; making the conversion costs an additional
$1/\sigma^{4}$ in overhead. So, for this low-error setting, we achieve better
$\sigma$ dependence, but if only worse error is required, their algorithm
performs better.

Shao and Montanaro [SM21] get a running time of
$O\MT@delim@RIGHT\Paren@star{\lVert
A\rVert_{F}^{6}/\sigma^{8}\varepsilon^{2}}$ which matches our
$\varepsilon$-dependence, obtains a better $\sigma$-dependence and a worse
Frobenius norm dependence. Additionally they requires that $b$ is in the image
of $A$.

In the context of Recommendation systems, the goal is to output a sample from
the rows of a low-rank approximation to $A$, denoted by $A_{\geqslant\sigma}$,
where we zero out all the singular values smaller than $\sigma$ (see Section
11.1 for a formal problem description). We then obtain the following
corollary:

######  Corollary 2.6 (De-quantizing Recommendation Systems, informal
Corollary 11.3).

Given a matrix $A\in\mathbb{C}^{m\times n}$ such that
$\MT@delim@RIGHT\norm@star{A}\leqslant 1$, an accuracy parameter
$\varepsilon$, and an $i\in[n]$, we can produce a data structure in
$O(\operatorname{nnz}(A))$ time such that, we can compute a vector $x_{i}$
such that with probability at least $0.9$,
$\MT@delim@RIGHT\Norm@star{A^{\dagger}x-[A_{\geqslant\sigma}]_{i,*}}\leqslant\varepsilon\MT@delim@RIGHT\norm@star{A_{i,*}}$
in $\tilde{O}\MT@delim@RIGHT\Paren@star{\lVert
A\rVert_{F}^{4}/\MT@delim@RIGHT\Paren@star{\sigma^{9}\varepsilon^{2}}}$ time.
Further, we can $\ell_{2}^{2}$-sample from $x_{i}$ in
$\tilde{O}\MT@delim@RIGHT\Paren@star{\lVert
A\rVert_{F}^{4}/\MT@delim@RIGHT\Paren@star{\sigma^{6}\varepsilon^{2}\MT@delim@RIGHT\norm@star{A^{\dagger}x}^{2}}}$
time.

######  Remark 2.7 (Comparison to [CGL+20a, CCH+22]).

Chia et. al. [CGL+20a] achieve a runtime of $\tilde{O}(\frac{\lVert
A\rVert_{F}^{6}}{\sigma^{16}\varepsilon^{6}})$. We improve upon it in every
parameter, including error $\varepsilon$, the threshold $\sigma$, and the
Frobenius norm $\lVert A\rVert_{F}$.

Chepurko et. al [CCH+22, Theorem 26] achieves a runtime that is at least
$\Omega(k^{3}/\varepsilon^{6})$ to get a low-rank approximation of an input
vector $M$ with the guarantee that $\lVert
A-M\rVert_{F}^{2}\leqslant(1+\varepsilon)\lVert A-A_{k}\rVert_{F}^{2}$. The
authors use that $\ell_{2}^{2}$ importance sampling sketches oversample ridge
leverage score sketch in certain parameter regimes, so this requires certain
addition assumptions on the size of $\lVert A_{k}\rVert_{F}$ and the residual
$\lVert A-A_{k}\rVert_{F}$. Work on quantum recommendation systems [KP17,
GSLW19] require a singular value threshold $\sigma$ instead of a rank
threshold $k$, and the standard way to convert between this “sketching”-style
error bound and the “QSVT”-style error bound is to bound $k\leqslant\lVert
A\rVert_{F}^{2}/\sigma^{2}$. Upon doing this, we see that the runtime is
$\tilde{O}(\frac{\lVert A\rVert_{F}^{6}}{\sigma^{6}\varepsilon^{6}})$. Our
work improves upon this in the $\lVert A\rVert_{F}$ and $\varepsilon$
parameters, but we lose a factor of $\sigma^{3}$.

Next, we state the de-quantization result we obtain for Hamiltonian
simulation:

######  Corollary 2.8 (Hamiltonian Simulation, informal Corollary 11.8).

Given a symmetric Hamiltonian $H\in\mathbb{C}^{n\times n}$ with
$\MT@delim@RIGHT\norm@star{H}\leqslant 1$ and a vector $b\in\mathbb{C}^{n}$,
we can output a description of a vector $v$ such that, with probability
$\geqslant 0.9$,
$\MT@delim@RIGHT\norm@star{v-e^{iHt}b}\leqslant\varepsilon\MT@delim@RIGHT\norm@star{b}$
in $\tilde{O}(t^{9}\lVert H\rVert_{F}^{4}/\varepsilon^{2})$.

We note that the only prior work [CGL+20a] we are aware of in the low-rank
regime obtains a running time
$O\MT@delim@RIGHT\Paren@star{t^{16}\MT@delim@RIGHT\norm@star{H}_{F}^{6}/\varepsilon^{6}}$,
and we improve upon it in every parameter.

##  3 Technical overview

In this section, we describe our classical framework for simulating QSVT and
provide an overview of our key new contributions. In brief: our main
conceptual contribution, using an iterative method (the Clenshaw recurrence)
instead of a pure sketching algorithm, is enough to achieve
$O(d^{13}\MT@delim@RIGHT\norm@star{A}_{F}^{4}/\varepsilon^{4})$, corresponding
to $d$ iterations of matrix-vector products of size $\lVert
A\rVert_{F}^{2}/(\varepsilon/d^{3})^{2}$ by $\lVert
A\rVert_{F}^{2}/(\varepsilon/d^{3})^{2}$. With insights into the stability of
the Clenshaw recurrence and sums of Chebyshev coefficients, we improve the
$\varepsilon/d^{3}$ to an $\varepsilon/(d^{2}\log^{2}(d))$. In the worst case,
we may need to rescale $\varepsilon$ to $\varepsilon/d^{2}$; in other words,
our stability analysis is tight up to $\log(d)$ factors. With insights into
matrix sparsification, we improve the $\varepsilon^{4}$ to an
$\varepsilon^{2}$, which is clearly tight. Together,555It is natural to wonder
here why the complexity is not something like
$d\MT@delim@RIGHT\norm@star{A}_{F}^{4}/(\varepsilon/(d^{2}\log^{2}(d)))^{2}=d^{5}\log^{2}(d)\MT@delim@RIGHT\norm@star{A}_{F}^{4}/\varepsilon^{2}$.
Such a runtime is conceivable, but our analysis essentially replaces two
factors of $1/\varepsilon$ with factors of $1/d^{2}$, so our sketching ideas
do not save any factors of $d$. this gives the final runtime of
$O(d^{9}\log^{4}(d)\MT@delim@RIGHT\norm@star{A}_{F}^{4}/\varepsilon^{2})$. We
begin by providing an informal description of the input access model.

#### Input model: Oversampling and query access.

Our goal is to give a classical version of the quantum singular value
transformation. In the quantum setting, we have a matrix
$A\in\mathbb{C}^{m\times n}$ with $\MT@delim@RIGHT\norm@star{A}\leqslant 1$
(given as a _block-encoding_) and a vector $b$ with
$\MT@delim@RIGHT\norm@star{b}=1$ (given encoded into the amplitudes of a
quantum state), and wish to compute $f(A)b$, where $f:[-1,1]\to\mathbb{R}$ is
an appropriately smooth function. The block-encoding access model is
inherently quantum, but there are standard ways to construct a block-encoding
for a matrix from classical data, including by assuming the matrix is sparse
with efficiently computable entries [GSLW19, Lemma 48] and given as a quantum
circuit [GSLW19, Definition 44]. Both of these support universal quantum
computation, and therefore cannot be simulated classically unless BPP=BQP.

However, block-encodings can also be achieved given an arbitrary matrix in a
data structure placed in quantum random access memory [GSLW19, Lemma 50], a
proposal for a quantum hardware architecture theorized to give the ability for
quantum computers to efficiently access stored memory in superposition
[GLM08]. As noted in prior work, this model incurs a square root of stable
rank, i.e. $O\MT@delim@RIGHT\Paren@star{\lVert
A\rVert_{F}/\MT@delim@RIGHT\norm@star{A}}$ overhead, allowing it to be
simulated classically with only polynomial overhead with sketching algorithms
[CGL+20a].

Chia et. al. [CGL+20a] proceed by introducing the access model of
_oversampling and query access_ , which can be interpreted as the classical
analogue of the block-encoding and amplitude-encoded quantum state in the
quantum setting. In particular, given a vector $v\in\mathbb{C}^{n}$,
_oversampling and query access_ corresponds to: (1) given an index $i\in[n]$,
outputting $v_{i}$, (2) sampling an index $j$ with probability
$\MT@delim@RIGHT\abs@star{v(j)}^{2}/\MT@delim@RIGHT\norm@star{v}^{2}$, and (3)
outputting $\MT@delim@RIGHT\norm@star{v}^{2}$. Similarly for a matrix
$A\in\mathbb{C}^{m\times n}$, _oversampling and query access_ corresponds to
having _oversampling and query access_ to all rows of A, as well as the vector
of the row norms. We point the reader to [CGL+20a] for an explanation of why
this model is the right classical analogue to benchmark quantum algorithms
against. In short, this model has closure properties very similar to that of
the block-encoding, and can be achieved whenever quantum states and generic
block-encodings can be prepared efficiently.

#### Computing polynomials of matrices.

The main primitive of QSVT is to take a block-encoding of a matrix $A$ and
give a block-encoding of an even or odd polynomial of that matrix, $p(A)$,
with an overhead of $\deg(p)$. When $A$ is asymmetric, we can interpret QSVT
as applying the matrix function that applies $p$ to each singular value of $A$
(Definition 5.1). In this way, quantum linear algebra algorithms can apply
generic functions to the singular values of a matrix, provided that they are
smooth enough to be approximated well by a low-degree polynomial.

In order to simulate QSVT classically, given a matrix
$A\in\mathbbm{C}^{m\times n}$, a vector $b\in\mathbbm{C}^{n}$, and a
polynomial $p:[-1,1]\to\mathbb{R}$, our goal is to compute some description
$p(A)b$. Specifically, we aim for our algorithm to run in
$\operatorname{poly}(\lVert A\rVert_{F},\frac{1}{\varepsilon},d,\log(mn))$
time after $O(\operatorname{nnz}(A)+\operatorname{nnz}(b))$ preprocessing, and
to output sampling and query access to $x$, where
$\MT@delim@RIGHT\norm@star{x-p(A)b}\leqslant\varepsilon\lVert
p\rVert_{\mathrm{sup}}\MT@delim@RIGHT\norm@star{b}$, where $\lVert
p\rVert_{\mathrm{sup}}=\max_{x\in[-1,1]}\MT@delim@RIGHT\abs@star{p(x)}$. We
note that as a byproduct, we also obtain a
$O\MT@delim@RIGHT\Paren@star{\operatorname{nnz}(A)+\operatorname{nnz}(b)+\operatorname{poly}(\lVert
A\rVert_{F},1/\varepsilon,d,\log(mn))}$ algorithm for the task of outputting
$x$.

We require the running time of the algorithm to be independent of input
dimension (after the preprocessing) and therefore are compelled to create
sketches of $A$ and $b$ and work with these sketches. We note that prior work
[Tan19, Tan21, CGL+20a] stops here, and directly computes a SVD of the
sketches, and applies the relevant polynomial $p$ to the singular values of
the sketch. As noted in the previous section, this approach loses large
polynomial factors in the relevant parameters.

#### Combining sketches with iterative algorithms.

Our main conceptual insight is to run iterative algorithms on the resulting
sketches of $A$ and $b$ in order to approximate matrix polynomials. In the
canonical numerical linear algebra regime (working with matrices directly,
instead of sketches), there are two standard methods to achieve this goal: (1)
compute the polynomial explicitly through something like a Clenshaw recurrence
[Cle55], or (2) use the Lanczos method [Lan50] to create a roughly-orthogonal
basis for the Krylov subspace $\\{b,Ab,A^{2}b,\ldots\\}$ and then apply the
function exactly to the matrix in this basis, in which $A$ is tridiagonal,
implicitly using a polynomial approximation in the analysis. We note that in a
world where we are allowed to pay $O(\operatorname{nnz}(A))$ (or even
$O(m+n)$) time per-iteration, we can simply run either of these algorithms and
call it a day. The main challenge in the setting we consider is that each
iterative step must run in time that is dimension-independent.

Given that we sketch each iterate down to a size that is dimension-
independent, we introduce additive error at each step. So, in essence, we must
perform a stability analysis of an iterative method, where the error
introduced in each iteration is from truncating the iterates. Finite-
precision/stability analysis of Clenshaw and Lanczos iteration are well-
understood [Cle55, Pai76, MMS18], and therefore one might expect to use these
analysis in a black-box manner. However, unlike typical finite-precision
analyses, which is concerned with error either at the granularity of “number
of bits to maintain”, and so is fine with polynomial loss, our final runtime
depends polynomially on the quality of our error analysis. This constraint
requires us to form a more intimate understanding of these iterative methods.

Folklore intuition suggests that Lanczos is a stabler algorithm for applying
matrix functions, but state-of-the-art analyses of it rely on the stability of
the Clenshaw recurrence as a subroutine (see [MMS18]), and therefore gives
strictly worse error-accumulation bounds than the Clenshaw recurrence. In
particular, if we wish to compute a generic $p(A)b$ to $\varepsilon$ error in
the regime where every matrix-vector product $Ax$ incurs an error of
$\varepsilon\MT@delim@RIGHT\norm@star{A}\MT@delim@RIGHT\norm@star{x}$ using
Lanczos, the stability analysis of Musco, Musco, and Sidford suggests that the
error of the output is $O(d^{5.5}\varepsilon)$, which would introduce a
$d^{11}$ in our setting [MMS18].666 This computation arises from taking Lemma
9 of [MMS18] to evaluate a degree-$d$ polynomial, say, bounded by 1 in
$[-1,1]$. A generic example of such polynomial is only by a constant in
$[-1-\eta,1+\eta]$ when $\eta=O(1/d^{2})$ (Lemma 5.4), and has Chebyshev
coefficients bounded only by a constant, without decaying. Thus, the bound
from Lemma 9 becomes $O(d^{5}\MT@delim@RIGHT\norm@star{E})$. Continuing the
analysis into [Pai76], $E$ is the matrix whose $i$th column is the error
incurred in the $i$th iteration; each column has norm
$\varepsilon\MT@delim@RIGHT\norm@star{A}\MT@delim@RIGHT\norm@nostar{v_{j+1}}=\varepsilon$
in our setting where we only incur error in matrix-vector products, since
$\MT@delim@RIGHT\norm@star{A}=1$ by normalizing and
$\MT@delim@RIGHT\norm@nostar{v_{j+1}}=1$ because the algorithm normalizes it
to unit norm, and we assume that scalar addition and multiplication can be
performed exactly. We have no control over the direction of error, so we can
at best bound $\MT@delim@RIGHT\norm@star{E}$ with the column-wise bound,
giving $\MT@delim@RIGHT\norm@star{E}\leqslant\lVert
E\rVert_{F}\leqslant\sqrt{k}\varepsilon$. So, our version of [MMS18, Equation
16] gives $\MT@delim@RIGHT\norm@star{E}\leqslant\sqrt{d}\varepsilon$, which
gives us the asserted $O(d^{5.5}\varepsilon)$ bound.  To incur less error, we
do not use Lanczos and analyze the Clenshaw recurrence directly.

#### Stability of the scalar Clenshaw recurrence.

The Chebyshev polynomials $\\{T_{\ell}(x)\\}_{\ell}$ form a basis and
therefore any degree-$d$ polynomial $p(x)$ can be written as a linear
combination of Chebyshev polynomials, i.e.
$p(x)=\sum_{\ell=0}^{d}a_{\ell}T_{\ell}(x)$. The Clenshaw recurrence computes
$p(x)$ through the iteration computing $q_{d}$ through to $q_{0}$:

$\displaystyle q_{d+1},q_{d+2}$ $\displaystyle\coloneqq 0;$ $\displaystyle
q_{k}$ $\displaystyle\coloneqq 2xq_{k+1}-q_{k+2}+a_{k};$ $\displaystyle p(x)$
$\displaystyle=\tfrac{1}{2}(a_{0}+q_{0}-q_{2}).$

For example, in the randomized numerical linear algebra (RNLA) literature,
this is often applied in the case where $a_{d}=1$ and $a_{k}=0$ otherwise, to
evaluate a degree-$d$ Chebyshev polynomial $T_{d}(x)$. Note that by Markov’s
inequality, a bounded polynomial has derivative bounded by $d^{2}$ [Sch41].
This is achieved for a Chebyshev polynomial, $p(x)=T_{d}(x)$. So, if our error
was only in changing $x$ to some value in $(x-\varepsilon,x+\varepsilon)$,
this error must cascade to a $O(d^{2}\varepsilon)$ worst-case error in the
output. This rough argument suggests that, in some sense, a $O(d^{2})$
overhead is the best we could hope for.

Our first technical contribution is an analysis showing that the Clenshaw
algorithm gives this optimal overhead, up to a logarithmic overhead. This
proceeds by showing that the overhead can be upper bounded by the size of the
largest Clenshaw iterate $\MT@delim@RIGHT\abs@star{q_{k}}$ (see Proposition
9.1), and then bounding the size of iterate $\MT@delim@RIGHT\abs@star{q_{k}}$
by $O(d\log(d)\lVert p\rVert_{\mathrm{sup}})$ (see Theorem 9.4). The main
lemma we prove states that for a bounded polynomial
$p(x)=\sum_{\ell=0}^{d}a_{\ell}T_{\ell}(x)$, sums of the form
$a_{\ell}+a_{\ell+2}+\cdots$ are all bounded by $O(\log(\ell)\lVert
p\rVert_{\mathrm{sup}})$ (7.3). This statement follows from facts in Fourier
analysis (in particular, this $\log(\ell)$ is the same $\log$ as the one
occurs when bounding the $L^{1}$ norm of the Dirichlet kernel). Finally, we
note that just using individual bounds on coefficients does not suffice, since
this would only give a bound of $O(d)$ on these sums, which can be
exponentially worse than the bound we obtain, and would give $O(d^{3})$
overhead (and a $d^{13}$ in the final runtime instead of a $d^{9}$).

We note that we are not aware of any prior work where such a sharp analysis
appears in the literature. The standard literature either considers an
additive error (where, unlike usual models like floating-point arithmetic,
each multiplication incurs identical error regardless of magnitude) [Cle55,
FP68, MH02] or eventually boils down to bounding
$\MT@delim@RIGHT\abs@star{a_{i}}$ (since their main concern is dependence on
$x$) [Oli77, Oli79], which is insufficient to get our $d^{2}\log(d)$ stability
bound. The modern work we are aware of shows a $O(d^{2})$ bound only for
Chebyshev polynomials [BKM22], sometimes used to give a $O(d^{3})$ bound for
computing generic bounded polynomials [MMS18], since a degree-$d$ polynomial
can be written as a linear combination of $T_{k}(x)$ with bounded
coefficients.

#### Extending Clenshaw to compute matrix polynomials.

So far, we have only considered the Clenshaw recurrence for a scalar $x$.
Generalizing this to matrices is fairly straightforward: for a matrix $A$ (not
necessarily square), we wish to apply an even or odd polynomial to it, and
then apply it to a vector $b$. Then, we can use a corresponding variant of the
Clenshaw recurrence to compute it. For example, the matrix Clenshaw recurrence
for computing $p(x)=\sum_{\ell=0}^{d}a_{2\ell+1}T_{2\ell+1}(x)$ is

$\displaystyle u_{d+1},u_{d+2}$ $\displaystyle\coloneqq 0;$ $\displaystyle
u_{k}$ $\displaystyle\coloneqq 2(2AA^{\dagger}-I)u_{k+1}-u_{k+2}+2a_{2k+1}Ab;$
$\displaystyle p(A)b=u$ $\displaystyle\coloneqq\tfrac{1}{2}(u_{0}-u_{1}).$

The question becomes how to perform this iteration efficiently and stably, in
a dimension-independent regime. We begin by sketching down our matrix and
vector: we show that it suffices to maintain a sparse description of $u_{k}$
of the form $u_{k}=Av_{k}$ where $v_{k}$ is sparse. In particular, we produce
sketches $S\in\mathbb{C}^{n\times s}$ and $T\in\mathbb{C}^{t\times m}$ such
that

  1. 1.

$\MT@delim@RIGHT\norm@nostar{AS\MT@delim@RIGHT\Paren@star{AS}^{\dagger}-AA^{\dagger}}\leqslant\varepsilon$,

  2. 2.

$\MT@delim@RIGHT\norm@nostar{ASS^{\dagger}b-Ab}\leqslant\varepsilon$,

  3. 3.

$\MT@delim@RIGHT\norm@nostar{TAS\MT@delim@RIGHT\Paren@star{TAS}^{\dagger}-AS\MT@delim@RIGHT\Paren@star{AS}^{\dagger}}\leqslant\varepsilon$

Sketches that satisfy the above type of guarantees are called _approximate
matrix product_ (AMP) sketches, and are standard in the quantum-inspired
algorithms literature [CGL+20a]. We observe that if we have sampling and query
access to $A$ and $b$, then we can produce these sketches of size
$s,t=O(\frac{\lVert
A\rVert_{F}^{2}}{\varepsilon^{2}}\log(n)\log\frac{1}{\delta})$, and then
compute $TAS$ in $O(st)$ time. We also note that the first two guarantees
follow from observing that AMP sketches oversample the symmetric approximate
matrix product sketch by a factor of $2$, and thus both guarantees hold
simultaneously. The third guarantee is straight-forward and does not require
asymmetry. Using these guarantees we can sketch the iterates as follows:

$\displaystyle u_{k}$
$\displaystyle=2(2AA^{\dagger}-I)u_{k+1}-u_{k+2}+2a_{2k+1}Ab$
$\displaystyle=4AA^{\dagger}Av_{k+1}-2Av_{k+1}-Av_{k+2}+2a_{2k+1}Ab$
$\displaystyle\approx
AS[4(TAS)^{\dagger}(TAS)v_{k+1}-2v_{k+1}-v_{k+2}+2a_{2k+1}S^{\dagger}b].$

Therefore, we can interpret Clenshaw iteration as the recursion on the
dimension-independent term $v_{k}\approx
4(TAS)^{\dagger}(TAS)v_{k+1}-2v_{k+1}-v_{k+2}+2a_{2k+1}S^{\dagger}b$, and then
applying $AS$ on the left to lift it back to $m$ dimensional space. We can
then analyze the sketched recursion to obtain an algorithm that runs in
$O(\lVert A\rVert_{F}^{4}/\varepsilon^{4})$ time per-iteration, not including
the loss from rescaling $\varepsilon$, whereas we wish to achieve a
$O(1/\varepsilon^{2})$ dependence per-iteration.

Though so far we have only used a very limited subset of the sketching
toolkit—namely, $\ell_{2}^{2}$ importance sampling—we remark that it’s not
clear how, for example, oblivious sketches [NN13] or the connection between
importance sampling and leverage score sampling [CCH+22] help us, since our
choices of sketches are optimal up to log factors for the guarantees we
desire. To get the additional improvement, we need a new sketching technique.

#### Improving the $\varepsilon$-dependence.

A natural next step to improve per-iteration runtime is to sparsify the matrix
$TAS$, in order to make matrix-vector products more efficient. If we can
sparsify $TAS$ to $O(1/\varepsilon^{2})$ non-zero entries, then we get the
desired quadratic savings in per-iteration cost.

There is significant literature on sparsifying the entries of a matrix [AM07,
KD14, BKKS21]. However, it does not suffice for us to use these as a black
box. For example, consider the sketch given by Drineas and Zouzias [DZ11]: for
a matrix $M\in\mathbb{R}^{n\times n}$, zero out every entry smaller than
$\frac{\varepsilon}{2n}$, then sample entries proportional to their
$\ell_{2}^{2}$ magnitude, and consider the corresponding unbiased estimator of
$M$, denoted $\tilde{M}$. The guarantee is that the operator norms are close,
$\MT@delim@RIGHT\norm@star{M-\tilde{M}}\leqslant\varepsilon$, and the sparsity
is $O(n\log(n)\frac{\lVert A\rVert_{F}^{2}}{\varepsilon^{2}})$. This is
precisely the guarantee we need and the sketch can be performed efficiently;
however, this does not sparsify the matrix for us, since in our setting, our
matrices $TAS$ have dimension $\lVert A\rVert_{F}^{2}/\varepsilon^{2}$, so the
sparsity guarantee is only $O(\lVert
A\rVert_{F}^{4}/\varepsilon^{4}\ln(n))=O(st\ln(n))$. In other words, this
sketch gives us no improvement on sparsity!

#### Bi-linear entry-wise sampling transform.

Our second main technical contribution is to bypass this barrier by noticing
that we don’t need a guarantee as strong as a spectral norm bound. Instead, we
only need to achieve approximations of the form

$\MT@delim@RIGHT\norm@star{ASMv-AS\tilde{M}x_{k}}\leqslant\varepsilon,$ (1)

for various different choices of $x_{k}$. So, we only need our sketch
$\tilde{M}$ to approximate $M$ in some directions with good probability,
instead of approximating it in all directions. We define a simple sketch,
which we call the Bilinear Entry Sampling Transform (best777That is, we
pronounce ${\textsc{best}}(A)$ as “best of $A$”. We make no normative claim
about our sketch vis-a-vis other sketches.), that is an unbiased estimator for
$A$ linearly (rather than the product $A^{\dagger}A$) and achieves the
aforementioned guarantee (Equation (1)). This sketch is sampled from the same
distribution as the one of Drineas and Zouzias from above, but without the
zero-ing out small entries. In particular, we define

$M^{(k)}=\frac{1}{p_{i,j}}A_{i,j}e_{i}e_{j}^{\dagger}\hskip 14.45377pt\textrm{
with probability
}p_{i,j}=\frac{\MT@delim@RIGHT\abs@star{A_{i,j}}}{\MT@delim<EMAIL_ADDRESS>
Then, to get the sketch with sparsity $r$, we take the average of $r$
independent copies of the above random matrix:

${\textsc{best}}(A)=\frac{1}{r}\sum_{k\in[r]}M^{(k)}.$

This definition is not new: for example, one of the earliest papers on matrix
sparsification for linear algebra algorithms briefly considers this sketch
[AM07], and this sketch has been used implicitly in prior work on dequantizing
QML algorithms [Tan21]. However, as far as we know, our analysis of this
sketch for preserving bi-linear forms and saving a factor of
$1/\varepsilon^{2}$ is novel.

We show that ${\textsc{best}}(A)$ satisfies the following guarantees: taking
$r=\Omega(\lVert M\rVert_{F}^{2}/\varepsilon^{2})$, this sketch preserves the
bilinear form $u^{\dagger}Mv$ to
$\varepsilon\MT@delim@RIGHT\norm@star{u}\MT@delim@RIGHT\norm@star{v}$ error
with probability $\geqslant 0.9$ (Lemma 6.2). Second, taking $r=\Omega(\lVert
M\rVert_{F}^{2}n)$, this sketch preserves the norm
$\MT@delim@RIGHT\norm@star{Mv}$ to $0.1\MT@delim@RIGHT\norm@star{v}$ error
with probability $\geqslant 0.9$. So, by taking $r=\Omega(\lVert
M\rVert_{F}^{2}(n+\frac{1}{\varepsilon^{2}}))$, we can get both properties.
However, with this choice of $r$, ${\textsc{best}}(M)$ does not preserve the
spectral norm $\MT@delim@RIGHT\norm@star{M}$, even to constant error.888 At
this point, one might wonder whether one can simply zero out small entries to
get these guarantees with the additional spectral norm bound. This is not the
case; if we only zero out entries smaller than $\varepsilon/(2n)$, this
threshold is too small to improve the matrix Bernstein tail bound, whereas if
we zero out entries smaller than, say, $1/(100n)$, the matrix Bernstein tail
bound goes through, but the bilinear form is not $\varepsilon$-preserved.
This sketch can be interpreted as a relaxation of the approximate matrix
product and the sparsification sketches above, since we use it in regimes
where $\MT@delim@RIGHT\norm@star{{\textsc{best}}(A)}$ is unbounded and
$\MT@delim@RIGHT\norm@star{{\textsc{best}}(A)b}$ is not $\varepsilon$-close to
$\MT@delim@RIGHT\norm@star{Ab}$. However, if we use it in the “interior” of a
matrix product expression, as a proxy for approximate matrix product, it can
be used successfully to improve the $n/\varepsilon^{2}$ dependence of spectral
norm entry sparsification to something like $n+1/\varepsilon^{2}$. In our
setting, this corresponds to an improvement of $1/\varepsilon^{4}$ to
$1/\varepsilon^{2}$.

#### Applying best to Clenshaw Iteration.

Returning to our odd Clenshaw iteration, we had our approximate iterate

$\displaystyle u_{k}$ $\displaystyle\approx
AS[4(TAS)^{\dagger}(TAS)v_{k+1}-2v_{k+1}-v_{k+2}+2a_{2k+1}S^{\dagger}b].$ We
can approximate this by taking $B={\textsc{best}}(TAS)$ and
$B_{\dagger}={\textsc{best}}((TAS)^{\dagger})$ with sparsity $r=\Theta(\lVert
A\rVert_{F}^{4}/\varepsilon^{2})$ to get that the bilinear forms like
$[AS]_{i,*}(TAS)^{\dagger}((TAS)v_{k+1})$ are preserved. This allows us to
successfully approximate with sparsified matrices, $\displaystyle\approx
AS[4B_{\dagger}Bv_{k+1}-2v_{k+1}-v_{k+2}+2a_{2k+1}S^{\dagger}b].$
$\displaystyle\text{so }v_{k}$ $\displaystyle\approx
4B_{\dagger}Bv_{k+1}-2v_{k+1}-v_{k+2}+2a_{2k+1}S^{\dagger}b.$

This recurrence in $v_{k}$ will be our algorithm (Algorithm 10.4): compute $S$
and $T$ in the pre-processing phase, then compute the iteration of $v_{k}$’s,
pulling fresh copies of $B$ and $B_{\dagger}$ each time, and output it as the
description of the output $u\approx p(A)b$. What remains is the error
analysis, which similar to the scalar case, requires bounding the size of the
iterates $\MT@delim@RIGHT\norm@star{v_{k}}$. We used the
$\varepsilon$-approximation of the bilinear form to show that the
sparsifications $B$ and $B_{\dagger}$ successfully approximate $u_{k}$; we use
the constant-error approximation of the norm to show that the
$\MT@delim@RIGHT\norm@star{v_{k}}$’s are bounded.

#### Challenges in extending finite-precision Clenshaw iteration.

We note here that the above error does not directly follow from the error of
the sketches along with the finite-precision scalar Clenshaw iteration. The
error we incur in iteration $k$ is not
$\varepsilon\MT@delim@RIGHT\norm@star{u_{k+1}}$, but
$\varepsilon\MT@delim@RIGHT\norm@star{v_{k+1}}$, so our error analysis
requires understanding both the original Clenshaw iteration along with the
“dual” Clenshaw iteration $v_{k}$, which requires a separate error analysis.

#### Sums of Chebyshev coefficients.

One final technical wrinkle remains, which is to bound the error accumulation
of the matrix Clenshaw recurrences. In the same way that the original scalar
Clenshaw algorithm required bounding arithmetic progressions of Chebyshev
coefficients with step size two, $a_{\ell}+a_{\ell+2}+\cdots$ by the norm of
the corresponding polynomial $p=\sum_{\ell=0}^{d}a_{\ell}T_{\ell}(x)$, to
prove stability for the even and odd versions, we need to bound arithmetic
progressions with step size four. Surprisingly, this is significantly more
challenging.

We give a thorough explanation for this in Remark 7.8, but in brief, some
arithmetic progressions of Chebyshev coefficients arise naturally as linear
combinations of polynomial evaluations. For example, $\sum_{k\geqslant
0}a_{2k}=\tfrac{1}{2}(p(-1)+p(1))$, so we can conclude that
$\MT@delim@RIGHT\abs@star{\sum_{k\geqslant 0}a_{2k}}\leqslant\lVert
p\rVert_{\mathrm{sup}}$. Through Fourier analysis, this can be generalized to
progressions with different step sizes, but breaks for progressions with
certain offsets. The sum $\sum_{k\geqslant 0}a_{4k+1}$ is one such example,
and this is a quantity we need to bound to give a $\tilde{O}(d^{2})$ bound on
the iterates of odd matrix Clenshaw. A naïve bound of $O(d)$ follows from
bounding each coefficient separately, but this results in a significantly
worse running time down the line.

In Lemma 7.4, we show that we can bound the above sum $\sum_{k\geqslant
0}a_{4k+1}$ by $O(\log^{2}(d)\lVert p\rVert_{\mathrm{sup}})$. This shows that
this sum has near-total cancellation: despite our guarantee on coefficient
being merely that it is bounded by a constant,999 In fact, for any degree $d$,
there are polynomials of that degree which are bounded and yet
$\sum_{\ell=0}^{d}\MT@delim@RIGHT\abs@star{a_{\ell}}=\Omega(d)$ [Tre19,
Theorems 8.2 and 8.3].  the sum of $O(d)$ of these coefficients is only poly-
logarithmic in magnitude. The proof proceeds by considering many different
polynomial evaluations $p(x_{0}),p(x_{1}),\ldots,p(x_{d})$, and trying to
write the arithmetic progression as a linear combination of these evaluations
$\sum c_{k}p(x_{k})$. This can be thought of as a linear system in the
$c_{k}$’s, which we can then prove has a solution, and then bound the
$\ell_{1}$ norm of this solution. To do this, we argue that we can bound its
solution $A^{-1}b$ by the solution of a different system, $C^{-1}b$, where $C$
is an matrix that is entrywise larger than $A$. Since this is not true
generally, we use strong properties of the particular matrix $A$ in the linear
system at hand to prove this.

##  4 Related work

#### Quantum machine learning.

Starting with the breakthrough work of Harrow, Hassidim and Lloy [HHL09] for
solving sparse linear systems, QML has seen a growing number of proposals.
Several works assume access to _quantum random access memory_ and provide
quantum algorithms for Principal Component Analysis [LMR14], support vector
machines [RML14], k-means clustering [LMR13], quantum recommendation systems
[KP17], and neural networks [WKS16, WKS14]. We refer the reader to [BWP+17]
for a comprehensive survey. In particular, we note that the current proposals
that resist de-quantization and potentially obtain a super-polynomial quantum
speedup include Zhao, Fitzsimons, and Fitzsimons on Gaussian process
regression [ZFF19] and Lloyd, Garnerone, and Zanardi on topological data
analysis [LGZ16].

#### Quantum-inspired algorithms.

Since the breakthrough work of Tang [Tan19] for recommendation systems, there
has been a flurry of work obtaining _quantum-inspired_ classical algorithms
for various problems in machine learning and numerical linear algebra. The
works most closely related to ours are the quantum-inspired framework for QSVT
by Chia, Gilyén, Li, Lin, Tang and Wang [CGL+20a], and a de-quantization of
QSVT when the input matrices are sparse by Gharibian and Le Gall [GLG22]. A
series of papers also focus on dequantizing specific problems, such as
regression [CGL+20b, CGL+20a, SM21, GST22, CCH+22], recommendation system
[Tan19, CGL+20a, CCH+22] and principal component analysis [Tan21]. Classical
algorithms for quantum simulation include Van den Nest’s work on simulating
restricted classes of quantum circuits using probabilistic methods [VdN11] and
Rudi, Wossnig, Ciliberto, Rocchetto, Pontil, and Severini’s work on using the
Nyström method to simulate a sparse Hamiltonian $H$ on a sparse input state
[RWC+20]. .

#### Randomized numerical linear algebra.

Our work is draws upon several ideas from the randomized numerical linear
algebra literature and we refer the reader to the following surveys: [Mah11,
Woo14] for the relevant background. The asymmetric approximate matrix product
sketch we introduce is closely related to the approximate matrix product
sketches considered in [MI11, MZ11, CEM+15] (non-oblivious) and [CNW15]
(oblivious). We note that in our setting, we cannot afford to use oblivious
sketches, since we lose the ability to maintain sampling and query access to
the resulting vectors. Cohen, Nelson, and Woodruff obtain an AMP guarantee in
the symmetric case using ridge leverage score sampling, whereas we work with
asymmetric matrices and would need to pay condition number dependent
oversampling factors to sample from the ridge leverage score distribution.
Magdon-Ismail [MI11] and Magen and Zouzias [MZ11] construct asymmetric
approximate matrix products, but require sampling from
$\Set{\MT@delim@RIGHT\norm@star{A_{i}}^{2}+\MT@delim@RIGHT\norm@star{B_{i}}^{2}/(\sum_{i\in[n]}\MT@delim@RIGHT\norm@star{A_{i}}^{2}+\MT@delim@RIGHT\norm@star{B_{i}}^{2})}$
and
$\Set{\MT@delim@RIGHT\norm@star{A_{i}}\MT@delim@RIGHT\norm@star{B_{i}}/\sum_{i\in[n]}\MT@delim@RIGHT\norm@star{A_{i}}\MT@delim@RIGHT\norm@star{B_{i}}}$
respectively. We note that we do not have efficient access to these
distributions in the _sampling and query access_ model. Drineas, Kannan and
Mahoney [DKM06] use $\ell_{2}^{2}$ sampling to get the assymetric approximate
matrix product guarantee, but under the Frobenius norm instead of spectral
norm, which yet again does not suffice for our purposes. Regarding the bi-
linear entry-wise sampling transform, creating sketches by sampling entries
proportional to their squared-magnitude are well-known in the literature but
are typically used to bound the operator norm of the matrix rather than any
particular bi-linear form. Bounding the operator norm directly happens to be
too restrictive, as we discussed in the technical overview. Regardless, entry-
wise sampling sketches were introduced by Achlioptas and McSherry [AM07] to
speed up low-rank approximation. Arora, Hazan and Kale used them in the
context of faster algorithms for SDP solving [AHK06]. Since then, a series of
works have provided sharper guarantees for such sketches [GT09, DZ11, KD14,
KDMI17, BKKS21]. Finally, there has been a recent flurry of work on sublinear
time algorithms for low-rank approximation under various structural
assumptions on the input [MW17, BW18, IVWW19, BCW20], however all such
algorithms still incur at least a linear dependence on the number of rows and
columns of the input matrix.

In addition to the usual sketching and sampling toolkit, we also use ideas
from iterative algorithms in the numerical linear algebra literature.
Iterative methods, such as power iteration, Krylov subspace methods, Golub-
Kahan Bidiagonalization, Arnoldi iteration, and the Lanczos iteration are
ubiquitous in scientific computing and are used for matrix inversion, solving
linear systems, linear programming, low-rank approximation, and numerous other
fundamental linear algebra primitives [Saa81, TBI97]. Our work is closely
related to iterative algorithms that use the sparsification sketches and
exploit singular values gaps, such as [AZL16, GS18, MM15, BCW22], however
these results incur dimension dependent factors, which are crucial to avoid in
our setting.

##  5 Preliminaries

We use the notation $f\lesssim g$ to denote the ordering
$f=\mathcal{O}\left(g\right)$ (and respectively for $\gtrsim$ and $\eqsim$).

###  5.1 Linear algebra

For vectors $v\in\mathbb{C}^{n}$, $\|v\|$ denotes standard Euclidean norm (so
$\|v\|\coloneqq(\sum_{i=1}^{n}\MT@delim@RIGHT\abs@star{v_{i}}^{2})^{1/2}$).
For a matrix $A\in\mathbb{C}^{m\times n}$, the Frobenius norm of $A$ is
$\lVert
A\rVert_{F}\coloneqq(\sum_{i=1}^{m}\sum_{j=1}^{n}\MT@delim@RIGHT\abs@star{A_{i,j}}^{2})^{1/2}$
and the spectral norm of $A$ is
$\MT@delim@RIGHT\norm@star{A}\coloneqq\sup_{x\in\mathbbm{C}^{n},\|x\|=1}\|Ax\|$.

A _singular value decomposition_ (SVD) of $A$ is a representation
$A=UDV^{{\dagger}}$, where for $N\coloneqq\min(m,n)$,
$U\in\mathbbm{C}^{m\times N}$ and $V\in\mathbbm{C}^{n\times N}$ are isometries
and $D\in\mathbbm{R}^{N\times N}$ is diagonal with $\sigma_{i}\coloneqq
D_{i,j}$ and
$\sigma_{1}\geqslant\sigma_{2}\geqslant\cdots\geqslant\sigma_{N}\geqslant 0$.
We can also write this decomposition as
$A=\sum_{i=1}^{N}\sigma_{i}u_{i}v_{i}^{\dagger}$, where $u_{i}\coloneqq
U_{*,i}$ and $v_{i}\coloneqq V_{*,i}$.

We denote the set of singular values of $A$ by
$\operatorname{Spec}(A)\coloneqq\\{\sigma_{k}\\}_{k\in[N]}$. For a Hermitian
matrix $A\in\mathbbm{C}^{n\times n}$ and a function
$f:\mathbbm{R}\to\mathbbm{C}$, $f(A)$ refers to applying $f$ to the
eigenvalues of $A$.

###  5.2 Polynomials and the Chebyshev Basis

We consider polynomials with complex coefficients, $p\in\mathbbm{C}[x]$. For a
Hermitian matrix $A$, $p(A)$ refers to evaluating the polynomial with $x$
replacing $A$; this is equivalent to applying $p$ to the eigenvalues of $A$.
The right definition for applying $p$ to a general non-square matrix is
subtle; as done in QSVT, we restrict to settings where the matrix formed by
evaluating $p$ on the singular values of $A$ coincides with the evaluation of
a corresponding polynomial in $A$.

######  Definition 5.1 (Definition 6.1 of [CGL+20a]).

For a matrix $A\in\mathbbm{C}^{m\times n}$ and degree-$d$ polynomial
$p(x)\in\mathbbm{C}[x]$ of parity-$d$ (i.e., even if $d$ is even and odd if
$d$ is odd), we define the notation $p(A)$ in the following way:

  1. 1.

If $p$ is _even_ , meaning that we can express $p(x)=q(x^{2})$ for some
polynomial $q(x)$, then

$p(A)\coloneqq q(A^{\dagger}A)=p(\sqrt{A^{\dagger}A}).$

  2. 2.

If $p$ is _odd_ , meaning that we can express $p(x)=x\cdot q(x^{2})$ for some
polynomial $q(x)$, then

$p(A)\coloneqq A\cdot q(\sqrt{A^{\dagger}A}).$

For example, if $p(x)=x^{2}+1$, then $p^{(\mathrm{QV})}(A)=A^{\dagger}A+I$,
and if $p(x)=x^{3}+x$, then $p^{(\mathrm{QV})}(A)=AA^{\dagger}A+A$. Looking at
a singular value decomposition $A=\sum\sigma_{i}u_{i}v_{i}^{\dagger}$,
$p^{(\mathrm{QV})}(A)=\sum p(\sigma_{i})u_{i}v_{i}^{\dagger}$ when $p$ is odd
and $p^{(\mathrm{QV})}(A)=\sum p(\sigma_{i})v_{i}v_{i}^{\dagger}$ when $p$ is
even, thus making this definition coincide with the singular value
transformation as given in [GSLW19, Definition 16].

We work in the Chebyshev basis of polynomials throughout. Let $T_{\ell}(x)$
and $U_{\ell}(x)$ denote Chebyshev polynomials of the first and second kind,
respectively. They can be defined on $[-1,1]$ via

$\displaystyle T_{k}(\cos(\theta))$ $\displaystyle=\cos(n\theta)$ (2)
$\displaystyle U_{k}(\cos(\theta))$ $\displaystyle=\sin((n+1)x)/\sin(x),$ (3)

but we will give attention to their recursive definitions, since we will use
them for computation.

$\displaystyle T_{0}(x)$ $\displaystyle=1$ $\displaystyle U_{0}(x)$
$\displaystyle=1$ $\displaystyle T_{1}(x)$ $\displaystyle=x$ $\displaystyle
U_{1}(x)$ $\displaystyle=2x$ (4) $\displaystyle T_{k}(x)$
$\displaystyle=2x\cdot T_{k-1}(x)-T_{k-2}(x)$ $\displaystyle U_{k}(x)$
$\displaystyle=2x\cdot U_{k-1}(x)-U_{k-2}(x)$

For a function $f:[-1,1]\to\mathbb{R}$, we denote $\lVert
f\rVert_{\mathrm{sup}}\coloneqq\sup_{x\in[-1,1]}\MT@delim@RIGHT\abs@star{f(x)}$.
In this norm, the Chebyshev polynomials have $\lVert
T_{k}(x)\rVert_{\mathrm{sup}}=1$ and $\lVert
U_{k}(x)\rVert_{\mathrm{sup}}=n+1$.

Any Lipschitz continuous function101010 We call a function
$f:[-1,1]\to\mathbb{R}$ Lipschitz continuous if there exists a constant $C$
such that $\MT@delim@RIGHT\abs@star{f(x)-f(y)}\leqslant
C\MT@delim@RIGHT\abs@star{x-y}$ for $x,y\in[-1,1]$.  $f:[-1,1]\to\mathbb{R}$
can be written as a (unique) linear combination of Chebyshev polynomials,
$f(x)=\sum_{\ell}a_{\ell}T_{\ell}(x)$ (where we interpret $T_{\ell}(x)\equiv
0$ for negative $\ell$). When $f$ is a degree-$d$ polynomial, then
$a_{\ell}=0$ for all $\ell>d$. A common way to approximate a function is by
truncating the polynomial expansion; we denote this operation by
$f_{k}(x)=\sum_{\ell=0}^{k}a_{\ell}T_{\ell}(x)$, and we denote the remainder
to be
$\bar{f}_{k}(x)=f(x)-f_{k}(x)=\sum_{\ell=k+1}^{\infty}a_{\ell}T_{\ell}(x)$.

We use the following well-known properties of Chebyshev polynomials from Mason
and Handscomb [MH02].

$\displaystyle T_{i}(x)$ $\displaystyle=\tfrac{1}{2}(U_{i}(x)-U_{i-2}(x))$ (5)
$\displaystyle U_{i}(x)$ $\displaystyle T_{jk}(x)$
$\displaystyle=T_{j}(T_{k}(x))$ (6) $\displaystyle U_{2k+1}(x)$
$\displaystyle=U_{k}(T_{2}(x))U_{1}(x)=U_{k}(T_{2}(x))2x$ (7)
$\displaystyle\tfrac{d}{dx}T_{k}(x)$ $\displaystyle=kU_{k-1}(x)$ (8)

We recommend the book by Trefethen, and rely on its results to prove our
bounds. We use a couple of standard results about Chebyshev polynomials.

######  Lemma 5.2 (Coefficient bound, consequence of [Tre19, Eq. (3.12)]).

Let $f:[-1,1]\to\mathbb{R}$ be a Lipschitz continuous function. Then all its
Chebyshev coefficients $a_{k}$ are bounded:
$\MT@delim@RIGHT\abs@star{a_{k}}\leqslant 2\lVert f\rVert_{\mathrm{sup}}$.

######  Lemma 5.3 ([Tre19, Theorems 8.1 and 8.2]).

Let a function $f$ analytic in $[-1,1]$ be analytically continuable to the
open Bernstein ellipse
$E_{\rho}=\\{\frac{z+z^{-1}}{2}\mid\MT@delim@RIGHT\abs@star{z}<\rho\\}$, where
it satisfies $\MT@delim@RIGHT\abs@star{f(x)}\leqslant M$ for some $M$. Then
its Chebyshev coefficients satisfy $\MT@delim@RIGHT\abs@star{a_{0}}\leqslant
M$ and

$\displaystyle\MT@delim@RIGHT\abs@star{a_{k}}\leqslant 2M\rho^{-k},\quad
k\geqslant 1.$

Consequently,

$\displaystyle\lVert\bar{f}_{n}\rVert_{\mathrm{sup}}\leqslant\frac{2M\rho^{-n}}{\rho-1}.$

######  Lemma 5.4.

For a degree-$d$ polynomial $p$, and $\delta=\frac{1}{4d^{2}}$,

$\displaystyle\sup_{x\in[-1-\delta,1+\delta]}\MT@delim@RIGHT\abs@star{p(x)}\leqslant
e\lVert p\rVert_{\mathrm{sup}}.$

###### Proof.

Without loss of generality, take $\lVert p\rVert_{\mathrm{sup}}=1$. By
Proposition 2.4 in Sachadeva and Vishnoi [SV+14], and basic properties of
Chebyshev polynomials,

$\displaystyle\sup_{x\in[-1-\delta,1+\delta]}\MT@delim@RIGHT\abs@star{p(x)}\leqslant\sup_{x\in[-1-\delta,1+\delta]}\MT@delim@RIGHT\abs@star{T_{d}(x)}=T_{d}(1+\delta).$

Further, by Proposition 2.5 in [SV+14], we can evaluate $T_{d}(1+\delta)$ via
the formula

$\displaystyle T_{d}(x)$
$\displaystyle=\frac{1}{2}\Big{(}x+\sqrt{x^{2}-1}\Big{)}^{d}+\frac{1}{2}\Big{(}x-\sqrt{x^{2}-1}\Big{)}^{d}$
$\displaystyle T_{d}(1+\delta)$
$\displaystyle=\frac{1}{2}\Big{(}1+\delta+\sqrt{2\delta+\delta^{2}}\Big{)}^{d}+\frac{1}{2}\Big{(}1+\delta-\sqrt{2\delta+\delta^{2}}\Big{)}^{d}$
$\displaystyle\leqslant\exp\Big{(}d(\delta+\sqrt{2\delta+\delta^{2}})\Big{)}$
$\displaystyle\leqslant\exp\Big{(}\tfrac{1}{4d}+\sqrt{\tfrac{1}{2}+\tfrac{1}{16d^{2}}}\Big{)}\leqslant
e$

∎

###  5.3 Sampling and query access

We now introduce the “quantum-inspired” access model, following prior
exposition [CGL+20a]. We refer the reader there for a more thorough
investigation of this access model. From a sketching perspective, this model
encompasses “the set of algorithms that can be performed in time independent
of input dimension, using only $\ell_{2}^{2}$ sampling”, and is a decent
classical analogue for the input given to a quantum machine learning
algorithms operating on classical data.

######  Definition 5.5 (Sampling and query access to a vector, [CGL+20a,
Definition 3.2]).

For a vector $v\in\mathbbm{C}^{n}$, we have $\operatorname{SQ}(v)$, _sampling
and query access_ to $v$, if we can:

  1. 1.

query for entries of $v$;

  2. 2.

obtain independent samples $i\in[n]$ where we see $i$ with probability
$|v(i)|^{2}/\MT@delim@RIGHT\norm@star{v}^{2}$;

  3. 3.

query for $\|v\|$.

These samples are called something like “$\ell_{2}^{2}$ importance samples” in
the randomized numerical linear algebra literature. Quantumly, these can be
considered to be simulated measurements of the quantum state
$|v\rangle\coloneqq\frac{1}{\|v\|}\sum v_{i}|i\rangle$ in the computational
basis. Sampling and query access is closed under taking linear combinations,
once we introduce slack in the form of oversampling.

######  Definition 5.6 (Oversampling and query access, [CGL+20a, Definition
3.4]).

For $v\in\mathbbm{C}^{n}$ and $\phi\geqslant 1$, we have
$\operatorname{SQ}_{\phi}(v)$, _$\phi$ -oversampling and query access_ to $v$,
if we have the ability to query $v$ and $\operatorname{SQ}(\tilde{v})$ for
$\tilde{v}\in\mathbbm{C}^{n}$ a vector satisfying
$\|\tilde{v}\|^{2}=\phi\|v\|^{2}$ and
$\MT@delim@RIGHT\abs@star{\tilde{v}(i)}^{2}\geqslant\MT@delim@RIGHT\abs@star{v(i)}^{2}$
for all $i\in[n]$.

The $\ell_{2}^{2}$ distribution over $\tilde{v}$ $\phi$-oversamples the
distribution over $v$:

$\displaystyle\frac{\MT@delim@RIGHT\abs@star{\tilde{v}_{i}}^{2}}{\|\tilde{v}\|^{2}}=\frac{\MT@delim@RIGHT\abs@star{\tilde{v}_{i}}^{2}}{\phi\|v\|^{2}}\geqslant\frac{1}{\phi}\frac{\MT@delim<EMAIL_ADDRESS>
Intuitively speaking, estimators that use $\mathcal{D}_{v}$ can also use
$\mathcal{D}_{\tilde{v}}$ via rejection sampling at the expense of a factor
$\phi$ increase in the number of utilized samples. From this observation we
can prove that oversampling access implies an approximate version of the usual
sampling access:

######  Lemma 5.7 (Oversampling to sampling, [CGL+20a, Lemma 3.5]).

Suppose we are given $\operatorname{SQ}_{\phi}(v)$ and some $\delta\in(0,1]$.
We can sample from $\mathcal{D}_{v}$ with probability $\geqslant 1-\delta$ in
$\mathcal{O}\left(\phi\log\frac{1}{\delta}\right)$ queries to
$\operatorname{SQ}_{\phi}(v)$. We can also estimate $\|v\|$ to $\nu$
multiplicative error for $\nu\in(0,1]$ with probability $\geqslant 1-\delta$
in $\mathcal{O}\left(\frac{\phi}{\nu^{2}}\log\frac{1}{\delta}\right)$ queries.
Both of these algorithms take linear time in the number of queries.

Generally, compared to a quantum algorithm that can output (and measure) a
desired vector $|v\rangle$, our algorithms will output
$\operatorname{SQ}_{\phi}(u)$ such that $\|u-v\|$ is small. So, to dequantize
a quantum algorithm, we will need to bound $\phi$ to show that we can output
samples from $|v\rangle$. As for error, bounds on $\|u-v\|$ imply that
measurements from $u$ and $v$ follow distributions that are close in total
variation distance [Tan19, Lemma 4.1]. Now, we show that oversampling and
query access of vectors is closed under taking small linear combinations.

######  Lemma 5.8 (Linear combinations, [CGL+20a, Lemma 3.6]).

Given $\operatorname{SQ}_{\varphi_{t}}(v_{t})\in\mathbbm{C}^{n}$ and
$\lambda_{t}\in\mathbbm{C}$ for all $t\in[\tau]$, we have
$\operatorname{SQ}_{\phi}(\sum_{t=1}^{\tau}\lambda_{t}v_{t})$ for
$\phi=\tau\frac{\sum\varphi_{t}\|\lambda_{t}v_{t}\|^{2}}{\|\sum\lambda_{t}v_{t}\|^{2}}$.
After paying the pre-processing cost of querying for each of the norms of the
$\tilde{v}_{t}$’s, the cost of any query is equal to the cost of sampling from
any of the $v_{t}$’s plus the cost of querying an entry from all of the
$v_{t}$’s.

######  Definition 5.9 (Oversampling and query access to a matrix, [CGL+20a,
Definition 3.7]).

For a matrix $A\in\mathbbm{C}^{m\times n}$, we have $\operatorname{SQ}(A)$ if
we have $\operatorname{SQ}(A(i,\cdot))$ for all $i\in[m]$ and
$\operatorname{SQ}(a)$ for $a\in\mathbbm{R}^{m}$ the vector of row norms
($a(i)\\!\coloneqq\\!\|A(i,\cdot)\|$).

We have $\operatorname{SQ}_{\phi}(A)$ if we have $\operatorname{Q}(A)$ and
$\operatorname{SQ}(\tilde{A})$ for $\tilde{A}\in\mathbbm{C}^{m\times n}$
satisfying $\lVert\tilde{A}\rVert_{F}^{2}=\phi\lVert A\rVert_{F}^{2}$ and
$\MT@delim@RIGHT\abs@star{\tilde{A}(i,j)}^{2}\geqslant\MT@delim@RIGHT\abs@star{A(i,j)}^{2}$
for all $(i,j)\in[m]\times[n]$.

######  Remark 5.10.

QML algorithms achieve their largest speedups over classical algorithms when
given state preparation in $\operatorname{polylog}(mn)$ time. So, we concern
ourselves with this setting.

We can get $\operatorname{SQ}$ access to input matrices and vectors in input-
sparsity time. Given $v\in\mathbbm{C}^{n}$ in the standard RAM model, the
alias method [Vos91] takes $\Theta(\operatorname{nnz}(v))$ pre-processing time
to output a data structure that uses $\Theta(\operatorname{nnz}(v))$ space and
can sample from $v$ in $\Theta(1)$ time. In other words, we can get
$\operatorname{SQ}(v)$ with constant-time queries in
$\mathcal{O}\left(\operatorname{nnz}(v)\right)$ time, and by extension, for a
matrix $A\in\mathbbm{C}^{m\times n}$, $\operatorname{SQ}(A)$ with constant-
time queries in $\mathcal{O}\left(\operatorname{nnz}(v)\right)$
time.111111This holds in the word-RAM model, with an additional $\log$
overhead when considering bit complexity..

Therefore, the quantum-inspired setting can be directly translated to a basic
randomized numerical linear algebra algorithm. More precisely, with this data
structure, a fast quantum-inspired algorithm (say, one running in time
$\mathcal{O}\left(T\operatorname{\mathbf{sq}}(A)\right)$ for $T$ independent
of input size) implies an algorithm in the standard computational model
(running in $\mathcal{O}\left(\operatorname{nnz}(A)+T\right)$ time).

######  Corollary 5.11.

Suppose we are given sampling and query access to a matrix
$A\in\mathbbm{C}^{m\times n}$ and a vector $b\in\mathbbm{C}^{n}$, where we can
respond to queries in $O(1)$ time. Further suppose we have a vector $u\in
C^{n}$ implicitly represented by $v\in\mathbbm{C}^{m}$ and $\eta$, with
$u=A^{\dagger}v+\eta b$. Then by Lemma 5.8, we have
$\operatorname{SQ}_{\phi}(u)$ for

$\displaystyle\phi=(\lVert
v\rVert_{0}+1)\frac{\sum_{k}\|v_{k}A_{k}\|^{2}+\|\eta b\|^{2}}{\|u\|^{2}}$

and a query cost of $O(\lVert v\rVert_{0})$. In particular, by Lemma 5.7 we
can draw one sample from $u$ with probability $\geqslant 1-\delta$ in
$O(\lVert v\rVert_{0}\phi\log\frac{1}{\delta})$ time.

This factor is only large when the linear combination has significantly
smaller norm than the components $v_{t}$ in the sum suggest. Usually, in our
applications, we can intuitively think about this overhead being small when
the desired output vector mostly lies in a subspace spanned by singular
vectors with large singular values in our low-rank input. Quantum algorithms
also have the same kind of overhead. Namely, the QSVT framework encodes this
in the subnormalization constant $\alpha$ of block-encodings, and the overhead
from the subnormalization appears during post-selection [GSLW19]. When this
cancellation is not too large, the resulting overhead typically does not
affect too badly the runtime of our applications.

##  6 Extending the Sketching Toolkit

In this section, we show how to extend the modern sketching toolkit (see e.g.
[Woo14]) in two ways: (a) we provide a sub-sampling sketch that preserves bi-
linear forms with only a inverse quadratic dependence on $\varepsilon$ and (b)
a non-oblivious, $\ell_{2}^{2}$ sampling based asymmetric approximate matrix
product sketch.

###  6.1 The Bi-Linear Entry-wise Sampling Transform

######  Definition 6.1 (Bi-linear Entry-wise Sparsifying Transform).

For a matrix $A\in\mathbbm{C}^{m\times n}$, the best of $A$ with parameter $T$
is a matrix sampled as follows: for all $k\in[T]$,

$M^{(k)}=\frac{1}{p_{i,j}}A_{i,j}e_{i}e_{j}^{\dagger}\hskip 14.45377pt\textrm{
with probability
}p_{i,j}=\frac{\MT@delim@RIGHT\abs@star{A_{i,j}}}{\MT@delim@RIGHT\Norm@star{A}_{F}^{2}}$

Then,

${\textsc{best}}{}(A)=\frac{1}{T}\sum_{k\in[T]}M^{(k)}.$

######  Lemma 6.2 (Basic Properties of the Bi-Linear Entry-wise Sparsifying
Transform).

For a matrix $A\in\mathbbm{C}^{m\times n}$, let $M={\textsc{best}}{}(A)$ with
parameter $T$. Then, for $X\in\mathbbm{C}^{m\times m}$, $u\in\mathbbm{C}^{m}$
and $v\in\mathbbm{C}^{n}$, we have

$\displaystyle\operatorname{nnz}(M)$ $\displaystyle\leqslant T$ (9)
$\displaystyle\operatorname*{\mathbf{E}}\left[M\right]$ $\displaystyle=A$ (10)
$\displaystyle\operatorname*{\mathbf{E}}\left[M^{\dagger}XM-A^{\dagger}XA\right]$
$\displaystyle=\frac{1}{T}\MT@delim@RIGHT\Paren@star{\operatorname{Tr}(X)\|A\|_{F}^{2}I-A^{\dagger}XA}$
(11)
$\displaystyle\operatorname*{\mathbf{E}}\left[\MT@delim@RIGHT\Paren@star{u^{\dagger}Mv-\operatorname*{\mathbf{E}}\left[u^{\dagger}Mv\right]}^{2}\right]$
$\displaystyle=\frac{1}{T}\Big{(}\|A\|_{F}^{2}\|u\|^{2}\|v\|^{2}-(u^{\dagger}Av)^{2}\Big{)}$
(12)

###### Proof.

Observe, since $M$ is an average of $T$ sub-samples, each of which are
$1$-sparse, $M$ has at most $T$ non-zero entries. Next,

$\operatorname*{\mathbf{E}}\left[M\right]=\frac{1}{T}\sum_{k\in
T}\operatorname*{\mathbf{E}}\left[M^{(k)}\right]=\sum_{i\in[m]}\sum_{j\in[n]}p_{i,j}\frac{A_{i,j}}{p_{i,j}}e_{i}e_{j}^{\dagger}=A$
(13)

Similarly,

$\begin{split}\operatorname*{\mathbf{E}}\left[M^{\dagger}XM\right]&=\frac{1}{T^{2}}\operatorname*{\mathbf{E}}\left[\MT@delim@RIGHT\Paren@star{\sum_{k\in[T]}M^{(k)}}^{\dagger}X\MT@delim@RIGHT\Paren@star{\sum_{k\in[T]}M^{(k)}}\right]\\\
&=\frac{1}{T^{2}}\operatorname*{\mathbf{E}}\left[\MT@delim@RIGHT\Paren@star{\sum_{k,k^{\prime}\in[T]}\MT@delim@RIGHT\Paren@star{M^{(k)}}^{\dagger}XM^{(k^{\prime})}}\right]\\\
&=\frac{1}{T^{2}}\MT@delim@RIGHT\Paren@star{\MT@delim@RIGHT\Paren@star{\sum_{k\neq
k^{\prime}\in[T]}\operatorname*{\mathbf{E}}\left[M^{(k)}\right]^{\dagger}\cdot
X\cdot\operatorname*{\mathbf{E}}\left[M^{(k^{\prime})}\right]}+\MT@delim@RIGHT\Paren@star{\sum_{k\in[T]}\operatorname*{\mathbf{E}}\left[\MT@delim@RIGHT\Paren@star{M^{(k)}}^{\dagger}XM^{(k)}\right]}}\\\
&=\MT@delim@RIGHT\Paren@star{1-\frac{1}{T}}A^{\dagger}XA+\frac{1}{T}\sum_{i\in[m],j\in[n]}p_{i,j}\frac{A_{i,j}^{2}}{p_{i,j}^{2}}e_{j}e_{i}^{\dagger}Xe_{i}e_{j}^{\dagger}\\\
&=\MT@delim@RIGHT\Paren@star{1-\frac{1}{T}}A^{\dagger}XA+\frac{\MT@delim@RIGHT\Norm@star{A}_{F}^{2}}{T}\sum_{i\in[m],j\in[n]}X_{i,i}e_{j}e_{j}^{\dagger}\\\
&=\MT@delim@RIGHT\Paren@star{1-\frac{1}{T}}A^{\dagger}XA+\frac{\MT@delim@RIGHT\Norm@star{A}_{F}^{2}\operatorname{Tr}(X)}{T}I.\end{split}$
(14)

Finally, observe, for $i,i^{\prime}\in[m]$ and $j,j^{\prime}\in[n]$,

$M^{(k)}_{i,j}M^{(k)}_{i^{\prime},j^{\prime}}=\begin{cases}M_{i,j}^{2}\hskip
14.45377pt\textrm{if }(i,j)=(i^{\prime},j^{\prime})\\\ 0\hskip
28.90755pt\textrm{otherwise }\end{cases},$

and therefore,

$\begin{split}\operatorname*{\mathbf{E}}\left[\MT@delim@RIGHT\Paren@star{u^{\dagger}Mv-\operatorname*{\mathbf{E}}\left[u^{\dagger}Mv\right]}^{2}\right]&=\operatorname*{\mathbf{E}}\left[\MT@delim@RIGHT\Paren@star{u^{\dagger}(M-A)v}^{2}\right]\\\
&=\operatorname*{\mathbf{E}}\left[\sum_{i,i^{\prime}\in[m],j,j^{\prime}\in[n]}\MT@delim@RIGHT\Paren@star{M-A}_{i,j}\MT@delim@RIGHT\Paren@star{M-A}_{i^{\prime},j^{\prime}}u_{i}v_{j}u_{i^{\prime}}v_{j^{\prime}}\right]\\\
&=\operatorname*{\mathbf{E}}\left[\sum_{i\in[m],j\in[n]}\MT@delim@RIGHT\Paren@star{M}^{2}_{i,j}u_{i}^{2}v_{j}^{2}-\MT@delim@RIGHT\Paren@star{u^{\dagger}Av}^{2}\right]\\\
&=\frac{1}{T}\MT@delim@RIGHT\Paren@star{\MT@delim@RIGHT\norm@star{A}_{F}^{2}\MT@delim@RIGHT\norm@star{u}^{2}\MT@delim@RIGHT\norm@star{v}^{2}-\MT@delim@RIGHT\Paren@star{u^{\dagger}Av}^{2}}\end{split}$
(15)

∎

We list a simple consequence of these bounds that we use later.

######  Corollary 6.3.

For a matrix $A\in\mathbbm{C}^{m\times n}$, let $M={\textsc{best}}(A)$ with
parameter $T$. Then, for matrices $X\in\mathbbm{C}^{\ell\times m}$ and
$Y\in\mathbbm{C}^{n\times d}$,

$\displaystyle\operatorname{\mathbf{Pr}}\Big{[}\lVert
XMY\rVert_{F}\geqslant\lVert XAY\rVert_{F}+\frac{\lVert X\rVert_{F}\lVert
A\rVert_{F}\lVert Y\rVert_{F}}{\sqrt{\delta T}}\Big{]}\leqslant\delta$

###  6.2 Approximate Matrix Product via $\ell^{2}_{2}$ Sampling

In this subsection, we extend the well-known approximate matrix product (see
for instance [Woo14]) to the setting where we have an $\ell_{2}^{2}$-sampling
oracle. Typically, the approximate matrix product guarantee is achieved in the
oblivious sketching model, however, we cannot extend the oblivious sketching
to the quantum setting. Quantum-inspired algorithms follows row subsampling
literature, which can be performed in this setting. Formally, we show the
following:

######  Definition 6.4.

Given two matrices $A\in\mathbb{C}^{m\times n}$ and $B\in\mathbb{C}^{n\times
d}$, along with a probability distribution $p\in\mathbb{R}_{\geqslant 0}^{n}$,
we define the _Asymmetric Approximate Matrix Product_ of sketch size $s$,
denoted ${\textsc{AMP}}_{s}(A,B,p)$, to be the $n\times s$ matrix whose
columns are i.i.d. sampled according to the law

$\displaystyle[{\textsc{AMP}}_{s}(A,B,p)]_{*,j}=\frac{e_{k}}{\sqrt{s\cdot
p_{k}}}\text{ with probability }p_{k}$

For an $S={\textsc{AMP}}_{s}(A,B,p)$, we will typically consider the
expression $ASS^{\dagger}B$, which can be written as the sum of independent
rank-one matrices $\frac{1}{s\cdot p_{k}}A_{*,k}B_{k,*}$.

######  Theorem 6.5 (Asymmetric Approximate Matrix Multiplication).

Given matrices $A\in\mathbb{R}^{m\times n}$ and $B\in\mathbb{R}^{n\times d}$
such that $\MT@delim@RIGHT\Norm@star{A}=1,\MT@delim@RIGHT\Norm@star{B}=1$, let
${\textsc{AMP}}(A,B)\in\mathbb{R}^{m\times d}$ denote a sketch, where
${\textsc{AMP}}(A,B)=\frac{1}{k}\sum_{i\in[k]}x_{i}\otimes y_{i}$ such that
$x_{i}\otimes y_{i}=\frac{1}{p_{i}}a_{i}\otimes b_{i}$ with probability
$p_{i}\geqslant
c\MT@delim@RIGHT\Norm@star{a_{i}}^{2}/(2\MT@delim@RIGHT\Norm@star{A}^{2}_{F})+\MT@delim@RIGHT\Norm@star{b_{i}}^{2}/(2\MT@delim@RIGHT\Norm@star{B}^{2}_{F})$,
for a fixed constant $c$. Then, with probability at least $1-\delta$

$\|\tilde{AB}-AB\|\leqslant\sqrt{\frac{c^{\prime}\log(n)\log(1/\delta)\MT@delim@RIGHT\Paren@star{\|A\|_{F}^{2}+\|B\|_{F}^{2}}}{k}},$

for a fixed universal constant $c^{\prime}$.

To obtain this result, we prove the following key lemma:

######  Lemma 6.6 (Concentration of Asymmetric Random Outer Products).

Let $(x,y)$ be a tuple of random vectors, where $x\in\mathbb{R}^{m}$ and
$y\in\mathbb{R}^{d}$ such that

$\MT@delim@RIGHT\Paren@star{\max_{i\in[n]}\MT@delim@RIGHT\norm@star{x_{i}}^{2}\MT@delim@RIGHT\norm@star{\operatorname*{\mathbf{E}}\left[yy^{\top}\right]}+\max_{i\in[n]}\MT@delim@RIGHT\norm@star{y_{i}}^{2}\MT@delim@RIGHT\norm@star{\operatorname*{\mathbf{E}}\left[xx^{\top}\right]}}\leqslant
M^{2}.$

Let $\\{(x_{i},y_{i})\\}_{i\in[n]}$ be independent copies of the random
variable $(x,y)$. Then, for any $t\in(0,1)$,

$\operatorname{\mathbf{Pr}}\left[\MT@delim@RIGHT\Norm@star{\frac{1}{n}\sum_{i\in[n]}x_{i}\otimes
y_{i}-\operatorname*{\mathbf{E}}\left[x\otimes y\right]}\geqslant
t\right]\leqslant(m+d)\exp\MT@delim@RIGHT\Paren@star{-\frac{ct^{2}n}{M^{2}+\max_{i\in[n]}\MT@delim@RIGHT\norm@star{x_{i}\otimes
y_{i}}t}}$

###### Proof.

For $i\in[n]$, let $Z_{i}=\frac{1}{n}\MT@delim@RIGHT\Paren@star{x_{i}\otimes
y_{i}-\operatorname*{\mathbf{E}}\left[x\otimes y\right]}$. Further,
$\MT@delim@RIGHT\norm@star{Z_{i}}\leqslant\frac{2}{n}\MT@delim@RIGHT\norm@star{x_{i}\otimes
y_{i}}\leqslant\frac{2M}{n}$. Next, we bound the variance :

$\begin{split}\sigma^{2}&=\max\MT@delim@RIGHT\Paren@star{\underbrace{\MT@delim@RIGHT\norm@star{\sum_{i\in[n]}\operatorname*{\mathbf{E}}\left[Z_{i}Z_{i}^{\top}\right]}}_{(i)},\underbrace{\MT@delim@RIGHT\norm@star{\sum_{i\in[n]}\operatorname*{\mathbf{E}}\left[Z_{i}^{\top}Z_{i}\right]}}_{(ii)}}\end{split}$

We bound term $(i)$ as follows:

$\begin{split}\MT@delim@RIGHT\norm@star{\sum_{i\in[n]}\operatorname*{\mathbf{E}}\left[Z_{i}Z_{i}^{\top}\right]}&=\MT@delim@RIGHT\norm@star{\frac{1}{n}\operatorname*{\mathbf{E}}\left[\MT@delim@RIGHT\Paren@star{x_{i}y_{i}^{\top}-\operatorname*{\mathbf{E}}\left[x\otimes
y\right]}\MT@delim@RIGHT\Paren@star{x_{i}y_{i}^{\top}-\operatorname*{\mathbf{E}}\left[x\otimes
y\right]}^{\top}\right]}\\\
&\leqslant\frac{2}{n}\MT@delim@RIGHT\norm@star{\operatorname*{\mathbf{E}}\left[\MT@delim@RIGHT\norm@star{y_{i}}^{2}x_{i}x_{i}^{\top}\right]}\\\
&\leqslant\frac{\max_{i\in[n]}\MT@delim@RIGHT\norm@star{y_{i}}^{2}}{n}\operatorname*{\mathbf{E}}\left[xx^{\top}\right]\end{split}$
(16)

We bound term $(ii)$ as follows:

$\begin{split}\MT@delim@RIGHT\norm@star{\sum_{i\in[n]}\operatorname*{\mathbf{E}}\left[Z_{i}^{\top}Z_{i}\right]}&=\MT@delim@RIGHT\norm@star{\frac{1}{n}\operatorname*{\mathbf{E}}\left[\MT@delim@RIGHT\Paren@star{x_{i}y_{i}^{\top}-\operatorname*{\mathbf{E}}\left[x\otimes
y\right]}^{\top}\MT@delim@RIGHT\Paren@star{x_{i}y_{i}^{\top}-\operatorname*{\mathbf{E}}\left[x\otimes
y\right]}\right]}\\\
&\leqslant\frac{2}{n}\MT@delim@RIGHT\norm@star{\operatorname*{\mathbf{E}}\left[\MT@delim@RIGHT\norm@star{x_{i}}^{2}y_{i}y_{i}^{\top}\right]}\\\
&\leqslant\frac{\max_{i\in[n]}\MT@delim@RIGHT\norm@star{x_{i}}^{2}}{n}\operatorname*{\mathbf{E}}\left[yy^{\top}\right]\end{split}$
(17)

Let
$M^{2}\geqslant\MT@delim@RIGHT\Paren@star{\max_{i\in[n]}\MT@delim@RIGHT\norm@star{x_{i}}^{2}\MT@delim@RIGHT\norm@star{\operatorname*{\mathbf{E}}\left[yy^{\top}\right]}+\max_{i\in[n]}\MT@delim@RIGHT\norm@star{y_{i}}^{2}\MT@delim@RIGHT\norm@star{\operatorname*{\mathbf{E}}\left[xx^{\top}\right]}}$.
Applying Matrix Bernstein (see Fact 6.7),

$\begin{split}\operatorname{\mathbf{Pr}}\left[\MT@delim@RIGHT\Norm@star{\sum_{i\in[n]}Z_{i}}\geqslant
t\right]&=\operatorname{\mathbf{Pr}}\left[\MT@delim@RIGHT\Norm@star{\frac{1}{n}\sum_{i\in[n]}x_{i}\otimes
y_{i}-\operatorname*{\mathbf{E}}\left[x\otimes y\right]}\geqslant t\right]\\\
&\leqslant(m+d)\exp\MT@delim@RIGHT\Paren@star{-\frac{ct^{2}n}{M^{2}+\max_{i\in[n]}\MT@delim@RIGHT\norm@star{x_{i}\otimes
y_{i}}t}},\end{split}$

for a fixed constant $c\geqslant 1$. ∎

######  Fact 6.7 (Matrix Bernstein, Theorem 1.6 [Tro12]).

Given a sequence of independent $d_{1}\times d_{2}$ random matrices
$\Set{Z_{i}}_{i\in[k]}$ such that for all $i\in[k]$,
$\operatorname*{\mathbf{E}}\left[Z_{i}\right]=0$ and
$\MT@delim@RIGHT\norm@star{Z_{i}}\leqslant L$ almost surely, let

$\sigma^{2}=\max\MT@delim@RIGHT\Paren@star{\MT@delim@RIGHT\norm@star{\sum_{i\in[k]}\operatorname*{\mathbf{E}}\left[Z_{i}Z_{i}^{\top}\right]},\MT@delim<EMAIL_ADDRESS>
Then, for any $t\geqslant 0$,

$\operatorname{\mathbf{Pr}}\left[\MT@delim@RIGHT\norm@star{\sum_{i\in[k]}Z_{i}}\geqslant
t\right]\leqslant\MT@delim@RIGHT\Paren@star{d_{1}+d_{2}}\exp\MT@delim@RIGHT\Paren@star{-\frac{t^{2}/2}{\sigma^{2}+Lt/3}}$

It is now straight-forward to prove Theorem 6.5 using the aforementioned
lemma:

###### Proof of Theorem 6.5.

Setting $x_{i}=\frac{1}{\sqrt{p_{i}}}\cdot a_{i}$ and
$y_{i}=\frac{1}{\sqrt{p_{i}}}\cdot b_{i}$, where
$p_{i}=\frac{\MT@delim@RIGHT\norm@star{a_{i}}^{2}}{2\MT@delim@RIGHT\norm@star{A}_{F}^{2}}+\frac{\MT@delim@RIGHT\norm@star{b_{i}}^{2}}{2\MT@delim@RIGHT\norm@star{B}_{F}^{2}}$
, we have

$\begin{split}\MT@delim@RIGHT\norm@star{x_{i}}&\leqslant
2\MT@delim@RIGHT\norm@star{A}_{F}\\\
\MT@delim@RIGHT\norm@star{y_{i}}&\leqslant
2\MT@delim@RIGHT\norm@star{B}_{F}\\\ \MT@delim@RIGHT\norm@star{x_{i}\otimes
y_{i}}^{1/2}&\leqslant\MT@delim@RIGHT\norm@star{\frac{\MT@delim@RIGHT\norm@star{A}_{F}\MT@delim@RIGHT\norm@star{B}_{F}}{\MT@delim@RIGHT\norm@star{a_{i}}\MT@delim@RIGHT\norm@star{b_{i}}}a_{i}\otimes
b_{i}}^{1/2}\leqslant\sqrt{\MT@delim@RIGHT\Norm@star{A}_{F}\cdot\MT@delim<EMAIL_ADDRESS>
Further,

$\begin{split}\operatorname*{\mathbf{E}}\left[x\otimes
x\right]&=\sum_{i\in[n]}\frac{1}{2}\MT@delim@RIGHT\Paren@star{\frac{\MT@delim@RIGHT\norm@star{a_{i}}^{2}}{\MT@delim@RIGHT\norm@star{A}_{F}^{2}}+\frac{\MT@delim@RIGHT\norm@star{b_{i}}^{2}}{\MT@delim@RIGHT\norm@star{B}_{F}^{2}}}\cdot\frac{a_{i}\otimes
a_{i}}{p_{i}}\\\ &=AA^{\top},\end{split}$ (18)

and similarly

$\begin{split}\operatorname*{\mathbf{E}}\left[y\otimes
y\right]&=\sum_{i\in[n]}\frac{1}{2}\MT@delim@RIGHT\Paren@star{\frac{\MT@delim@RIGHT\norm@star{a_{i}}^{2}}{\MT@delim@RIGHT\norm@star{A}_{F}^{2}}+\frac{\MT@delim@RIGHT\norm@star{b_{i}}^{2}}{\MT@delim@RIGHT\norm@star{B}_{F}^{2}}}\cdot\frac{b_{i}\otimes
b_{i}}{p_{i}}b_{i}\otimes b_{i}\\\ &=BB^{\top}.\end{split}$ (19)

Observe, setting
$M^{2}=2\MT@delim@RIGHT\Paren@star{\MT@delim@RIGHT\norm@star{A}_{F}^{2}+\MT@delim@RIGHT\norm@star{B}_{F}^{2}}$
suffices, and applying Matrix Bernstein, we have

$\operatorname{\mathbf{Pr}}\left[\MT@delim@RIGHT\norm@star{AB-\tilde{AB}}\geqslant
t\right]\leqslant(m+d)\exp\MT@delim@RIGHT\Paren@star{-\frac{ct^{2}k}{\MT@delim@RIGHT\norm@star{A}_{F}^{2}+\MT@delim<EMAIL_ADDRESS>(20)

Setting
$t=\sqrt{\frac{c\log(1/\delta)\log(m+d)\MT@delim@RIGHT\Paren@star{\MT@delim@RIGHT\norm@star{A}_{F}^{2}+\MT@delim@RIGHT\norm@star{B}_{F}^{2}}}{k}}$,
we know that with probability at least $1-\delta$

$\MT@delim@RIGHT\norm@star{AB-\tilde{AB}}\leqslant\sqrt{\frac{c\log(1/\delta)\log(m+d)\MT@delim@RIGHT\Paren@star{\MT@delim@RIGHT\norm@star{A}_{F}^{2}+\MT@delim@RIGHT\norm@star{B}_{F}^{2}}}{k}},$

as desired. ∎

##  7 Sums of Chebyshev coefficients

To give improved stability bounds for the Clenshaw recurrence, we need to
bound various sums of Chebyshev coefficients. Since we aim to give bounds that
hold for all degree-$d$ polynomials, we use no property of the function beyond
that it has a unique Chebyshev expansion; of course, for any particular choice
of function $f$, the bounds in this section can be improved by explicitly
computing its Chebyshev coefficients, or in some cases, by using smoothness
properties of the function [Tre19, Theorems 7.2 and 8.2].

Let $f:[-1,1]\to\mathbb{R}$ be a Lipschitz continuous function. Then it can be
expressed uniquely as a linear combination of Chebyshev polynomials
$f(x)=\sum_{i=0}^{\infty}a_{i}T_{i}(x)$. A broad topic of interest in
approximation theory is bounds for linear combinations of these coefficients,
$\sum a_{i}c_{i}$, in terms of $\lVert f\rVert_{\mathrm{sup}}$; this was one
motivation of Vladimir Markov in proving the Markov brothers’ inequality
[Sch41, p575]. Our goal for this section will be to investigate this question
in the case where these sums are arithmetic progressions of step four. This
will be necessary for later stability analyses, and is one of the first non-
trivial progressions to bound. We begin with some straightforward assertions
(see [Tre19] for background).

######  Fact 7.1.

Let $f:[-1,1]\to\mathbb{R}$ be a Lipschitz continuous function. Then its
Chebyshev coefficients $\\{a_{\ell}\\}_{\ell}$ satisfy

$\begin{split}\MT@delim@RIGHT\abs@nostar[\Big{]}{\sum_{\ell}a_{\ell}}&=\MT@delim@RIGHT\abs@star{f(1)}\leqslant\lVert
f\rVert_{\mathrm{sup}}\\\
\MT@delim@RIGHT\abs@nostar[\Big{]}{\sum_{\ell}(-1)^{\ell}a_{\ell}}&=\MT@delim@RIGHT\abs@star{f(-1)}\leqslant\lVert
f\rVert_{\mathrm{sup}}\\\
\MT@delim@RIGHT\abs@nostar[\Big{]}{\sum_{\ell}a_{\ell}\MT@delim@RIGHT\iver@star{\ell\text{
is
even}}}&=\MT@delim@RIGHT\abs@nostar[\Big{]}{\sum_{\ell}a_{\ell}\frac{1}{2}(1+(-1)^{\ell})}\leqslant\lVert
f\rVert_{\mathrm{sup}}\\\
\MT@delim@RIGHT\abs@nostar[\Big{]}{\sum_{\ell}a_{\ell}\MT@delim@RIGHT\iver@star{\ell\text{
is
odd}}}&=\MT@delim@RIGHT\abs@nostar[\Big{]}{\sum_{\ell}a_{\ell}\frac{1}{2}(1-(-1)^{\ell})}\leqslant\lVert
f\rVert_{\mathrm{sup}}\end{split}$

We use the following result on Lebesgue constants to bound truncations of the
Chebyshev coefficient sums.

######  Lemma 7.2 ([Tre19, Theorem 15.3]).

Let $f:[-1,1]\to\mathbb{R}$ be a Lipschitz continuous function, let
$f_{k}(x)=\sum_{\ell=0}^{k}a_{\ell}T_{\ell}(x)$, and let optimal degree-$k$
approximating polynomial to $f$ be denoted $f_{k}^{*}$. Then

$\displaystyle\lVert f-f_{k}\rVert_{\mathrm{sup}}$
$\displaystyle\leqslant\Big{(}4+\frac{4}{\pi^{2}}\log(k+1)\Big{)}\lVert
f-f_{k}^{*}\rVert_{\mathrm{sup}}$
$\displaystyle\leqslant\Big{(}4+\frac{4}{\pi^{2}}\log(k+1)\Big{)}\lVert
f\rVert_{\mathrm{sup}}.$

Similarly,

$\displaystyle\lVert f_{k}\rVert_{\mathrm{sup}}$ $\displaystyle\leqslant\lVert
f-f_{k}\rVert_{\mathrm{sup}}+\lVert
f\rVert_{\mathrm{sup}}\leqslant\Big{(}5+\frac{4}{\pi^{2}}\log(k+1)\Big{)}\lVert
f\rVert_{\mathrm{sup}}.$

This implies bounds on sums of coefficients.

######  Fact 7.3.

Consider a function $f(x)=\sum_{\ell}a_{\ell}T_{\ell}(x)$. Then

$\displaystyle\MT@delim@RIGHT\abs@nostar[\Big{]}{\sum_{\ell=k}^{\infty}a_{\ell}\MT@delim@RIGHT\iver@star{\ell-k\text{
is even}}}\leqslant\lVert
f-f_{k-2}\rVert_{\mathrm{sup}}\leqslant\Big{(}4+\frac{4}{\pi^{2}}\log(k-1)\Big{)}\lVert
f\rVert_{\mathrm{sup}},$

where the inequalities follow from 7.1 and Lemma 7.2. When $k=0,1$, then the
sum is bounded by $\lVert f\rVert_{\mathrm{sup}}$, as shown in 7.1.

Now, we prove similar bounds in the case that $f(x)$ is an odd function. In
particular, we want to obtain a bound on alternating signed sums of the
Chebshyev coefficients and we incur a blowup that scales logarithmically in
the degree.

######  Lemma 7.4.

Let $f:[-1,1]\to\mathbb{R}$ be an _odd_ Lipschitz continuous function with
Chebyshev coefficients $\\{a_{\ell}\\}_{\ell}$, so that $a_{k}=0$ for all even
$k$. Then the Chebyshev coefficient sum is bounded as

$\displaystyle\MT@delim@RIGHT\abs@nostar[\Big{]}{\sum_{\ell=0}^{d}(-1)^{\ell}a_{2\ell+1}}$
$\displaystyle\leqslant(\ln(d)+2)\max_{0\leqslant k\leqslant 2d+1}\lVert
f_{k}\rVert_{\mathrm{sup}}$
$\displaystyle\leqslant(\ln(d)+2)\Big{(}5+\frac{4}{\pi^{2}}\ln(2d+2)\Big{)}\lVert
f\rVert_{\mathrm{sup}}$
$\displaystyle\leqslant\Big{(}16+4\ln^{2}(d+1)\Big{)}\lVert
f\rVert_{\mathrm{sup}}.$

We first state the following relatively straight-forward corollary:

######  Corollary 7.5.

Lemma 7.4 gives bounds on arithmetic progressions with step size four. Let
$f:[-1,1]\to\mathbb{R}$ be a Lipschitz continuous function, and consider
nonnegative integers $c\leqslant d$. Then

$\displaystyle\MT@delim@RIGHT\abs@nostar[\Big{]}{\sum_{\ell=c}^{d}a_{\ell}\MT@delim@RIGHT\iver@star{\ell-c\equiv
0\ (\mathrm{mod}\ 4)}}$ $\displaystyle\leqslant(32+8\ln^{2}(d+1))\lVert
f\rVert_{\mathrm{sup}}$

###### Proof.

Define $f^{\mathrm{odd}}\coloneqq\tfrac{1}{2}(f(x)-f(-x))$ and
$f^{\mathrm{even}}\coloneqq\tfrac{1}{2}(f(x)+f(-x))$ to be the odd and even
parts of $f$ respectively. Triangle inequality implies that $\lVert
f^{\mathrm{odd}}\rVert_{\mathrm{sup}},\lVert
f^{\mathrm{even}}\rVert_{\mathrm{sup}}\leqslant\lVert f\rVert_{\mathrm{sup}}$.
Suppose $c,d$ are odd. Then

$\displaystyle\MT@delim@RIGHT\abs@nostar[\Big{]}{\sum_{\ell=c}^{d}a_{\ell}\MT@delim@RIGHT\iver@star{\ell-c\equiv
0\ (\mathrm{mod}\ 4)}}$
$\displaystyle=\frac{1}{2}\MT@delim@RIGHT\abs@nostar[\Big{]}{\sum_{\ell=0}^{\lfloor(d-c)/2\rfloor}a_{c+2\ell}(1\pm(-1)^{\ell})}$
$\displaystyle\leqslant\frac{1}{2}\Big{(}\MT@delim@RIGHT\abs@nostar[\Big{]}{\sum_{\ell=0}^{\lfloor(d-c)/2\rfloor}a_{c+2\ell}}+\MT@delim@RIGHT\abs@nostar[\Big{]}{\sum_{\ell=0}^{\lfloor(d-c)/2\rfloor}(-1)^{\ell}a_{c+2\ell}}\Big{)}$
$\displaystyle\leqslant\frac{1}{2}\Big{(}\lVert
f_{c}^{\mathrm{odd}}\rVert_{\mathrm{sup}}+\lVert
f_{d}^{\mathrm{odd}}\rVert_{\mathrm{sup}}+2(\ln(d)+2)\max_{0\leqslant
k\leqslant d}\lVert f_{k}^{\mathrm{odd}}\rVert_{\mathrm{sup}}\Big{)}$
$\displaystyle\leqslant(32+8\ln^{2}(d+1))\lVert
f^{\mathrm{odd}}\rVert_{\mathrm{sup}}$
$\displaystyle\leqslant(32+8\ln^{2}(d+1))\lVert f\rVert_{\mathrm{sup}}$

The case when $c$ is even is easy, and it follows from 7.3: for , we know that

$\displaystyle\Big{\lVert}\sum_{\ell}a_{2\ell}T_{\ell}(x)\Big{\rVert}_{\mathrm{sup}}=\Big{\lVert}\sum_{\ell}a_{2\ell}T_{\ell}(T_{2}(x))\Big{\rVert}_{\mathrm{sup}}=\Big{\lVert}\sum_{\ell}a_{2\ell}T_{2\ell}(x)\Big{\rVert}_{\mathrm{sup}}=\Big{\lVert}f^{\mathrm{even}}(x)\Big{\rVert}_{\mathrm{sup}}\leqslant\lVert
f\rVert_{\mathrm{sup}},$

so

$\begin{split}\MT@delim@RIGHT\abs@nostar[\Big{]}{\sum_{\ell\geqslant
c}a_{\ell}\MT@delim@RIGHT\iver@star{\ell-c\equiv 0\ (\mathrm{mod}\
4)}}&=\MT@delim@RIGHT\abs@nostar[\Big{]}{\sum_{\ell\geqslant
c/2}a_{2\ell}\MT@delim@RIGHT\iver@star{\ell-c/2\text{ is even}}}\\\
&\leqslant\Big{(}4+\frac{4}{\pi^{2}}\log(c/2-1)\Big{)}\Big{\lVert}\sum_{\ell}a_{2\ell}T_{\ell}(x)\Big{\rVert}_{\mathrm{sup}}\\\
&\leqslant\Big{(}4+\frac{4}{\pi^{2}}\log(c/2-1)\Big{)}\lVert
f\rVert_{\mathrm{sup}},\end{split}$

and combining these two bounds for $c$ and $d$ gives the desired statement. ∎

We note that Lemma 7.4 will be significantly harder to prove. See Remark 7.8
for an intuitive explanation why. We begin with two structural lemmas on how
the solution to a unitriangular linear system behaves, which might be on
independent interest.

######  Lemma 7.6 (An entry-wise positive solution).

Suppose that $A\in\mathbb{R}^{d\times d}$ is an upper unitriangular matrix
such that, for all $i\leqslant j$, $a_{ij}>0$, $a_{ij}>a_{i-1,j}$. Then
$A^{-1}\vec{1}$ is a vector with positive entries.

The same result holds when $A$ is a lower unitriangular matrix such that, for
all $i\geqslant j$, $a_{ij}>0$, $a_{ij}>a_{i+1,j}$.

###### Proof.

Let $x=A^{-1}\vec{1}$. Then $x_{d}=1\geqslant 0$. The result follows by
induction:

$\displaystyle x_{i}$ $\displaystyle=1-\sum_{j=i+1}^{d}A_{ij}x_{j}$
$\displaystyle=\sum_{j=i+1}^{d}(A_{i+1,j}-A_{ij})x_{j}+1-\sum_{j=i+1}^{d}A_{i+1,j}x_{j}$
$\displaystyle=\sum_{j=i+1}^{d}(A_{i+1,j}-A_{ij})x_{j}+1-[Ax]_{i+1}$
$\displaystyle=\sum_{j=i+1}^{d}(A_{i+1,j}-A_{ij})x_{j}$ $\displaystyle>0$

For lower unitriangular matrices, the same argument follows. The inverse
satisfies $x_{1}=1$ and

$\displaystyle x_{i}$ $\displaystyle=1-\sum_{j=1}^{i-1}A_{ij}x_{j}$
$\displaystyle=\sum_{j=i+1}^{d}(A_{i-1,j}-A_{ij})x_{j}+1-\sum_{j=i+1}^{d}A_{i-1,j}x_{j}>0$

∎

Next, we characterize how the solution to a unitriangular linear system
behaves when we consider a partial ordering on the matrices.

######  Lemma 7.7.

Let $A$ be a nonnegative upper unitriangular matrix such that
$A_{ij}>A_{i-1,j}$ and $A_{ij}>A_{i,j+1}$ for all $i\leqslant j$. Let $B$ be a
matrix with the same properties, such that $A\geqslant B$ entrywise. By Lemma
7.6, $x^{(A)}=A^{-1}\vec{1}$ and $x^{(B)}=B^{-1}\vec{1}$ are nonnegative. It
further holds that
$\sum_{i=1}^{d}[A^{-1}\vec{1}]_{i}\leqslant\sum_{i=1}^{d}[B^{-1}\vec{1}]_{i}$.

###### Proof.

We consider the line between $A$ and $B$, $A(t)=A(1-t)+Bt$ for $t\in[0,1]$.
Let $x(t)=A^{-1}\vec{1}$; we will prove that $\vec{1}^{T}x(t)$ is
monotonically increasing in $t$. The gradient of $x(t)$ has a simple form
[Tao13]:

$\displaystyle A(t)x(t)$ $\displaystyle=\vec{1}$
$\displaystyle\partial[A(t)x(t)]$ $\displaystyle=\partial_{t}[\vec{1}]$
$\displaystyle(B-A)x(t)+A(t)\partial_{t}x(t)$ $\displaystyle=0$
$\displaystyle\partial_{t}x(t)$ $\displaystyle=A^{-1}(t)(A-B)x(t).$

So,

$\displaystyle\vec{1}^{T}\partial_{t}x(t)$
$\displaystyle=\vec{1}^{T}A^{-1}(t)(A-B)A^{-1}(t)\vec{1}$
$\displaystyle=[A^{-T}(t)\vec{1}]^{T}(A-B)[A^{-1}(t)\vec{1}].$

Since $A$ and $B$ satisfy the entry constraints, so do every matrix along the
line. Consequently, the column constraints in Lemma 7.6 are satisfied for both
$A$ and $A^{T}$, so both $A^{-T}(t)\vec{1}$ and $A^{-1}(t)\vec{1}$ are
positive vectors. Since $A\geqslant B$ entrywise, this means that
$\vec{1}^{T}\partial_{t}x(t)$ is positive, as desired. ∎

###### Proof of Lemma 7.4.

We first observe that the following sorts of sums are bounded. Let
$x_{k}:=\cos(\frac{\pi}{2}(1-\frac{1}{2k+1}))$. Then, using that
$T_{\ell}(\cos(x))=\cos(\ell x)$,

$\displaystyle f_{2k+1}(x_{k})$
$\displaystyle=\sum_{\ell=0}^{2k+1}a_{\ell}T_{\ell}(x_{k})$
$\displaystyle=\sum_{\ell=0}^{k}a_{2\ell+1}T_{2\ell+1}(x_{k})$
$\displaystyle=\sum_{\ell=0}^{k}a_{2\ell+1}\cos\Big{(}\frac{\pi}{2}\Big{(}2\ell+1-\frac{2\ell+1}{2k+1}\Big{)}\Big{)}$
$\displaystyle=\sum_{\ell=0}^{k}(-1)^{\ell}a_{2\ell+1}\sin\Big{(}\frac{\pi}{2}\frac{2\ell+1}{2k+1}\Big{)}.$

We have just shown that

$\displaystyle\MT@delim@RIGHT\abs@nostar[\Big{]}{\sum_{\ell=0}^{k}a_{2\ell+1}(-1)^{\ell}\sin\Big{(}\frac{\pi}{2}\frac{2\ell+1}{2k+1}\Big{)}}\leqslant\lVert
f_{2k+1}\rVert_{\mathrm{sup}}.$ (21)

We now claim that there exist non-negative $c_{k}$ for
$k\in\\{0,1,\ldots,d\\}$ such that

$\displaystyle\sum_{\ell=0}^{d}(-1)^{\ell}a_{2\ell+1}$
$\displaystyle=\sum_{k=0}^{d}c_{k}f_{2k+1}(x_{k}).$ (22)

The $f_{2k+1}(x_{k})$’s can be bounded using Lemma 7.2. The rest of the proof
will consist of showing that the $c_{k}$’s exist, and then bounding them.

To do this, we consider the coefficient of each $a_{2\ell+1}$ separately; let
$A^{(k)}\in[0,1]^{d+1}$ (index starting at zero) be the vector of coefficients
associated with $p_{2k+1}(x_{k})$:

$\displaystyle
A^{(k)}_{\ell}=\sin\Big{(}\frac{\pi}{2}\frac{2\ell+1}{2k+1}\Big{)}\text{ for
}0\leqslant\ell\leqslant k,\,0\text{ otherwise}$ (23)

Note that the $A^{(k)}_{\ell}$ is always non-negative and increasing with
$\ell$ up to $A_{k}^{(k)}=1$. Then Eq. 22 holds if and only if

$\displaystyle c_{0}A^{(0)}+\cdots+c_{d}A^{(d)}=\vec{1},$

or in other words, the equation $Ac=\vec{1}$ is satisfied, where $A$ is the
matrix with columns $A^{(k)}$ and $c$ is the vector of $c_{\ell}$’s. Since $A$
is upper triangular (in fact, with unit diagonal), this can be solved via
backwards substitution: $c_{d}=1$, then $c_{d-1}$ can be deduced from $c_{d}$,
and so on. More formally, the $s$th row gives the following constraint that
can be rewritten as a recurrence.

$\displaystyle\sum_{t=s}^{d}\sin\Big{(}\frac{\pi}{2}\frac{2s+1}{2t+1}\Big{)}c_{t}=1$
(24) $\displaystyle
c_{s}=1-\sum_{t=s+1}^{d}\sin\Big{(}\frac{\pi}{2}\frac{2s+1}{2t+1}\Big{)}c_{t}$
(25)

Because the entries of $A$ increase in $\ell$, the $c_{\ell}$’s are all
positive.

Invoking Lemma 7.6 with the matrix A establishes that such $c_{s}$ exist; our
goal now is to bound them. Doing so is not as straightforward as it might
appear: since the recurrence Eq. 25 _subtracts_ by $c_{t}$’s, an upper bound
on $c_{t}$ for $t\in[s+1,d]$ does not give an upper bound on $c_{s}$; it gives
a lower bound. So, an induction argument to show bounds for $c_{s}$’s fails.
Further, we were unable to find any closed form for this recurrence. However,
since all we need to know is the sum of the $c_{s}$’s, we show that we _can_
bound this via a generic upper bound on the recurrence.

Here, we apply Lemma 7.7 to $A$ as previously defined, and the bounding matrix
is (for $i\leqslant j)$

$\displaystyle B_{ij}$
$\displaystyle=\frac{i}{j}\leqslant\frac{2i+1}{2j+1}\leqslant\sin\Big{(}\frac{\pi}{2}\frac{2i+1}{2j+1}\Big{)}=A_{ij},$

using that $\sin(\frac{\pi}{2}x)\geqslant x$ for $x\in[0,1]$. Let
$\hat{c}=B^{-1}\vec{1}$. Then $\hat{c}_{i}=\frac{1}{i+1}$ for $i\neq d$ and
$\hat{c}_{d}=1$.

$\displaystyle[B\hat{c}]_{i}=\sum_{j=i}^{d}B_{ij}\hat{c}_{j}=\sum_{j=i}^{d-1}\frac{i}{j}\frac{1}{j+1}+\frac{i}{d}=i\sum_{j=i}^{d-1}\Big{(}\frac{1}{j}-\frac{1}{j+1}\Big{)}+\frac{i}{d}=i\Big{(}\frac{1}{i}-\frac{1}{d})+\frac{i}{d}=1$

By Lemma 7.7, $\sum_{i}c_{i}\leqslant\sum_{i}\hat{c}_{i}\leqslant\ln(d)+2$.
So, altogether, we have

$\begin{split}\MT@delim@RIGHT\abs@nostar[\Big{]}{\sum_{\ell=0}^{d}(-1)^{\ell}a_{2\ell+1}}&=\MT@delim@RIGHT\abs@nostar[\Big{]}{\sum_{k=0}^{d}c_{k}f_{2k+1}(x_{k})}\\\
&\leqslant\sum_{k=0}^{d}c_{k}\lVert f_{2k+1}\rVert_{\mathrm{sup}}\\\
&\leqslant\Big{(}\sum_{k=0}^{d}c_{k}\Big{)}\max_{0\leqslant k\leqslant
d}\lVert f_{2k+1}\rVert_{\mathrm{sup}}\\\
&=\Big{(}\sum_{k=0}^{d}c_{k}\Big{)}\max_{0\leqslant k\leqslant 2d+1}\lVert
f_{k}\rVert_{\mathrm{sup}}\\\ &\leqslant(\ln(d)+2)\max_{0\leqslant k\leqslant
2d+1}\lVert f_{k}\rVert_{\mathrm{sup}}\end{split}$

∎

######  Remark 7.8.

A curious reader will (rightly) wonder whether this proof requires this level
of difficulty. Intuition from the similar Fourier analysis setting suggests
that arithmetic progressions of any step size at any offset are easily
bounded. We can lift to the Fourier setting by considering, for an
$f:[-1,1]\to\mathbb{R}$, a corresponding $2\pi$-periodic
$g:[0,2\pi]\to\mathbb{R}$ such that

$\displaystyle g(\theta)\coloneqq
f(\cos(\theta))=\sum_{k=0}^{\infty}a_{k}T_{k}(\cos(\theta))=\sum_{k=0}^{\infty}a_{k}\cos(k\theta)=\sum_{k=0}^{\infty}a_{k}\frac{e^{ik\theta}+e^{-ik\theta}}{2}$

This function has the property that
$\MT@delim@RIGHT\abs@star{g(\theta)}\leqslant\lVert f\rVert_{\mathrm{sup}}$
and $\widehat{g}(k)=a_{\MT@delim@RIGHT\abs@star{k}}/2$ (except
$\widehat{g}(0)=a_{0}$). Consequently,

$\frac{1}{t}\sum_{j=0}^{t-1}f\Big{(}\cos(\frac{2\pi
j}{t})\Big{)}=\frac{1}{t}\sum_{j=0}^{t-1}g\Big{(}\frac{2\pi
j}{t}\Big{)}=\frac{1}{t}\sum_{j=0}^{t-1}\sum_{k=-\infty}^{\infty}\widehat{g}(k)e^{2\pi
ijk/t}=\sum_{k=-\infty}^{\infty}\widehat{g}(k)\sum_{j=0}^{t-1}\frac{1}{t}e^{2\pi
ijk/t}\\\
=\sum_{k=-\infty}^{\infty}\widehat{g}(k)\MT@delim@RIGHT\iver@star{k\text{ is
divisible by }t}=\sum_{k=-\infty}^{\infty}\widehat{g}(kt),$

so we can bound arithmetic progressions
$\MT@delim@RIGHT\abs@star{\sum_{k}\widehat{g}(kt)}\leqslant\lVert
f\rVert_{\mathrm{sup}}$, and this generalizes to other offsets, to bound
$\MT@delim@RIGHT\abs@star{\sum_{k}\widehat{g}(kt+o)}$ for some $o\in[t-1]$.
Notably, though, this approach does not say anything about sums like
$\sum_{k}a_{4k+1}$. The corresponding progression of Fourier coefficients
doesn’t give it, for example, since we pick up unwanted terms from the
negative Fourier coefficients.121212These sums are related to the Chebyshev
coefficients one gets from interpolating a function at Chebyshev points
[Tre19, Theorem 4.2].

$\displaystyle\sum_{k}\widehat{g}(4k+1)$
$\displaystyle=\MT@delim@RIGHT\parens@star{\widehat{g}(1)+\widehat{g}(5)+\widehat{g}(9)+\cdots}+\MT@delim@RIGHT\parens@star$
$\displaystyle=\frac{1}{2}\MT@delim@RIGHT\parens@star{a_{1}+a_{5}+a_{9}+\cdots}+\frac{1}{2}\MT@delim@RIGHT\parens@star{a_{3}+a_{7}+a_{11}+\cdots}=\sum_{k\geqslant
0}a_{2k+1}.$

In fact, by inspection of the distribution131313This is the functional to
integrate against to compute the sum,
$\tfrac{2}{\pi}\int_{-1}^{1}f(x)D(x)/\sqrt{1-x^{2}}=\sum a_{4k+1}$. The
distribution is not a function, but can be thought of as the limit object of
$D_{n}(x)=\sum_{k=0}^{n}T_{4k+1}(x)$ as $n\to\infty$, analogous to Dirichlet
kernels and the Dirac delta distribution.
$D(x)=\sum_{k=0}^{\infty}T_{4k+1}(x)$, it appears that this arithmetic
progression cannot be written as a linear combination of evaluations of
$f(x)$. Since the shape of the distribution appears to have $1/x$ behavior
near $x=0$, we conjecture that our analysis losing a log factor is, in some
respect, necessary.

######  Conjecture 7.9.

For any step size $t>1$ and offset $o\in[t-1]$ such that $o\neq t/2$, there
exists a function $f:[-1,1]\to\mathbb{R}$ such that $\lVert
f\rVert_{\mathrm{sup}}=1$ but
$\MT@delim@RIGHT\abs@star{\sum_{k=0}^{n}a_{tk+o}}=\Omega(\log(n))$.

##  8 Properties of the Clenshaw recursion

###  8.1 Deriving the Clenshaw recursions

Suppose we are given as input a degree-$d$ polynomial as a linear combination
of Chebyshev polynomials:

$\displaystyle p(x)=\sum_{k=0}^{d}a_{k}T_{k}(x).$ (26)

Then this can be computed with the _Clenshaw algorithm_ , which gives the
following recurrence.

$\displaystyle q_{d+1}$ $\displaystyle=q_{d+2}=0$ $\displaystyle q_{k}$
$\displaystyle=2xq_{k+1}-q_{k+2}+a_{k}$ (Clenshaw) $\displaystyle\tilde{p}$
$\displaystyle=\tfrac{1}{2}(a_{0}+q_{0}-q_{2})$

######  Lemma 8.1.

The recursion in Section 8.1 computes $p(x)$. That is, in exact arithmetic,
$\tilde{p}=p(x)$. In particular,

$\displaystyle q_{k}=\sum_{i=k}^{d}a_{i}U_{i-k}(x).$ (27)

###### Proof.

We show Eq. 27 by induction.

$\begin{split}q_{k}&=2xq_{k+1}-q_{k+2}+a_{k}\\\
&=2x\Big{(}\sum_{i=k+1}^{d}a_{i}U_{i-k-1}(x)\Big{)}-\Big{(}\sum_{i=k+2}^{d}a_{i}U_{i-k-2}(x)\Big{)}+a_{k}\\\
&=a_{k}+2xa_{k+1}U_{0}(x)+\sum_{i=k+2}^{d}a_{i}(2xU_{i-k-1}(x)-U_{i-k-2}(x))\\\
&=\sum_{i=k}^{d}a_{i}U_{i-k}(x).\end{split}$ (28)

Consequently, we have

$\begin{split}\frac{1}{2}(a_{0}+u_{0}-u_{2})&=\frac{1}{2}\Big{(}a_{0}+\sum_{i=0}^{d}a_{i}U_{i}(x)-\sum_{i=2}^{d}a_{i}U_{i-2}(x)\Big{)}\\\
&=a_{0}+a_{1}x+\sum_{i=2}^{d}\frac{a_{i}}{2}(U_{i}(x)-U_{i-2}(x))\\\
&=\sum_{i=0}^{d}a_{i}T_{i}(x).\end{split}$

∎

######  Remark 8.2.

Though the aforementioned discussion is specialized to the scalar setting, it
extends to the the matrix setting almost entirely syntactically: consider a
Hermitian $A\in\mathbb{C}^{n\times n}$ and $b\in\mathbb{C}^{n}$ with
$\|A\|,\|b\|\leqslant 1$. Then $p(A)b$ can be computed in the following way:

$\begin{split}u_{d+1}&=\vec{0}\\\ u_{d}&=a_{d}b\\\
u_{k}&=2Au_{k+1}-u_{k+2}+a_{k}b\\\
u:=p(A)b&=\frac{1}{2}(a_{0}b+u_{0}-u_{2})\end{split}$ (29)

The proof that this truly computes $p(A)b$ is the same as the proof of
correctness for Clenshaw’s algorithm shown above.

###  8.2 Evaluating even and odd polynomials

We will be considering evaluating odd and even polynomials. We again focus on
the scalar setting and note that this extends to the matrix setting in the
obvious way. The previous recurrence Section 8.1 can work in this setting, but
it’ll be helpful for our analysis if the recursion multiplies by $x^{2}$ each
time, instead of $x$ [MH02, Chapter 2, Problem 7]. So, in the case where the
degree-$(2d+1)$ polynomial $p(x)$ is _odd_ (so $a_{2k}=0$ for every $k$), it
can be computed with the iteration

$\displaystyle q_{d+1}$ $\displaystyle=q_{d+2}=0$ $\displaystyle q_{k}$
$\displaystyle=2T_{2}(x)q_{k+1}-q_{k+2}+a_{2k+1}U_{1}(x)$ (Odd Clenshaw)
$\displaystyle\tilde{p}$ $\displaystyle=\tfrac{1}{2}(q_{0}-q_{1})$

When $p(x)$ is a degree-$(2d)$ _even_ polynomial (so $a_{2k+1}=0$ for every
$k$), it can be computed via the same recurrence, replacing $a_{2k+1}U_{1}(x)$
with $a_{2k}$. However, we will use an alternative form that’s more convenient
for us (since we can reuse the analysis of the odd case).

$\displaystyle\tilde{a}_{2k}$ $\displaystyle\coloneqq
a_{2k}-a_{2k+2}+a_{2k+4}-\cdots\pm a_{2d}$ (30) $\displaystyle q_{d+1}$
$\displaystyle=q_{d+2}=0$ $\displaystyle q_{k}$
$\displaystyle=2T_{2}(x)q_{k+1}-q_{k+2}+\tilde{a}_{2k+2}U_{1}(x)^{2}$ (Even
Clenshaw) $\displaystyle\tilde{p}$
$\displaystyle=\tilde{a}_{0}+\tfrac{1}{2}(q_{0}-q_{1})$

These recurrences correctly compute $p$ follows from a similar analysis to the
standard Clenshaw algorithm, formalized below.

######  Lemma 8.3.

The recursions in Section 8.2 and Eq. Even Clenshaw correctly compute $p(x)$
for even and odd polynomials, respectively. That is, in exact arithmetic,
$\tilde{p}=p(x)$. In particular,

$\displaystyle q_{k}=\sum_{i=k}^{d}a_{i}U_{i-k}(x).$ (31)

###### Proof.

We can prove these statements by applying Eq. 27. In the odd case, Section 8.2
is identical to Section 8.1 except that $x$ is replaced by $T_{2}(x)$ and
$a_{k}$ is replaced by $a_{2k+1}U_{1}(x)$, so by making the corresponding
changes in the iterate, we get that

$\displaystyle q_{k}$
$\displaystyle=\sum_{i=k}^{d}a_{2i+1}U_{1}(x)U_{i-k}(T_{2}(x))=\sum_{i=k}^{d}a_{2i+1}U_{2(i-k)+1}(x)$
by Eq. 7 $\displaystyle\tilde{p}$
$\displaystyle=\tfrac{1}{2}(q_{0}-q_{1})=\sum_{i=0}^{d}\frac{a_{2i+1}}{2}\Big{(}U_{2i+1}(x)-U_{2i-1}(x)\Big{)}=p(x).$
by Eq. 5

Similarly, in the even case, Eq. Even Clenshaw is identical to Section 8.1
except that $x$ is replaced by $T_{2}(x)$ and $a_{k}$ is replaced by
$4\tilde{a}_{2k}x^{2}$ (see Definition 30), so that

$\displaystyle q_{k}$
$\displaystyle=\sum_{i=k}^{d}\tilde{a}_{2i+2}U_{1}(x)^{2}U_{i-k}(T_{2}(x))$
(32) $\displaystyle=\sum_{i=k}^{d}\tilde{a}_{2i+2}U_{1}(x)U_{2(i-k)+1}(x)$ by
Eq. 7
$\displaystyle=\sum_{i=k}^{d}\tilde{a}_{2i+2}(U_{2(i-k)}(x)+U_{2(i-k+1)}(x))$
by Eq. 4
$\displaystyle=\sum_{i=k}^{d+1}\tilde{a}_{2i+2}U_{2(i-k)}(x)+\sum_{i=k+1}^{d+1}\tilde{a}_{2i}U_{2(i-k)}(x)$
noticing that $\tilde{a}_{2d+2}=0$
$\displaystyle=\tilde{a}_{2k+2}+\sum_{i=k+1}^{d+1}(\tilde{a}_{2i}+\tilde{a}_{2i+2})U_{2(i-k)}(x)$
(33) $\displaystyle=\tilde{a}_{2k+2}+\sum_{i=k+1}^{d}a_{2i}U_{2(i-k)}(x)$ (34)
$\displaystyle=-\tilde{a}_{2k}+\sum_{i=k}^{d}a_{2i}U_{2(i-k)}(x)$ (35)

Finally, observe

$\displaystyle\tilde{a}_{0}+\frac{1}{2}(q_{0}-q_{1})$
$\displaystyle=\tilde{a}_{0}+\frac{1}{2}(a_{0}-\tilde{a}_{0}+\tilde{a}_{2})+\sum_{i=1}^{d}\frac{a_{2i}}{2}(U_{2i}(x)-U_{2i-2}(x))=p(x).$
(36)

∎

##  9 Stability of the scalar Clenshaw recursion

Before we move to the matrix setting, we warmup with a stability analysis of
the scalar Clenshaw recursion. Suppose we perform Section 8.1 to compute a
degree-$d$ polynomial $p$, except every addition, subtraction, and
multiplication incurs $\varepsilon$ relative error. Typically, this has been
analyzed in the finite precision setting, where the errors are caused by
truncation. These standard analyses show that this finite precision recursion
gives $p(x)$ to
$d^{2}(\sum\MT@delim@RIGHT\abs@star{a_{i}})\varepsilon=O(d^{3}\lVert
p\rVert_{\mathrm{sup}}\varepsilon)$ error. A use of Parseval’s formula [MH02,
Theorem 5.3] $\sum\MT@delim@RIGHT\abs@star{a_{i}}\leqslant\sqrt{d}\sqrt{\sum
a_{i}^{2}}=O(\sqrt{d}\lVert p\rVert_{\mathrm{sup}})$ reduces this by a factor
of $\sqrt{d}$.

However, this bound on $\sum\MT@delim@RIGHT\abs@star{a_{i}}$ is essentially
tight, and moreover, it is tight for a relevant class of polynomials:
truncations of well-behaved functions. Such polynomials generally have the
property that
$\MT@delim@RIGHT\abs@star{a_{k}}=\Theta((1-\log(1/\varepsilon)/d)^{-k})$
(Lemma 5.2), so most of the coefficients are constant-sized, eventually
decaying to $\varepsilon$.

We improve on prior stability analyses to show that the Clenshaw recurrence
for Chebyshev polynomials only incurs an error overhead of $d^{2}\log(d)$.
This is tight up to a logarithmic factor. This, for example, could be used to
improve the bound in [MMS18, Lemma 9] from $k^{3}$ to $k^{2}\log(k)$ (where in
that paper, $k$ denotes degree).

###  9.1 Analyzing error propagation

The following is a simple analysis of Clenshaw, with some rough resemblance to
an analysis of Oliver [Oli79].

######  Proposition 9.1 (Stability Analysis for Scalar Clenshaw).

Consider a degree-$d$ polynomial $p:[-1,1]\to\mathbb{R}$ with Chebyshev
coefficients $p(x)=\sum_{k=0}^{d}a_{k}T_{k}(x)$. Let
$\oplus,\ominus,\odot:\mathbb{C}\times\mathbb{C}\to\mathbb{C}$ be binary
operations representing addition, subtraction, and multiplication to
$\mu\varepsilon$ relative error, for $0<\varepsilon<1$:

$\displaystyle\MT@delim@RIGHT\abs@star{(x\oplus y)-(x+y)}$
$\displaystyle\leqslant\mu\varepsilon(\MT@delim@RIGHT\abs@star{x}+\MT@delim@RIGHT\abs@star{y})$
$\displaystyle\MT@delim@RIGHT\abs@star{(x\ominus y)-(x-y)}$
$\displaystyle\leqslant\mu\varepsilon(\MT@delim@RIGHT\abs@star{x}+\MT@delim@RIGHT\abs@star{y})$
$\displaystyle\MT@delim@RIGHT\abs@star{x\odot y-x\cdot y}$
$\displaystyle\leqslant\mu\varepsilon\MT@delim@RIGHT\abs@star{x}\MT@delim@RIGHT\abs@star{y}=\mu\varepsilon\MT@delim<EMAIL_ADDRESS>
Given an $x\in[-1,1]$, consider performing the Clenshaw recursion with these
noisy operations:

$\displaystyle\tilde{q}_{d+1}$ $\displaystyle=\tilde{q}_{d+2}=0$
$\displaystyle\tilde{q}_{k}$ $\displaystyle=(2\odot
x)\odot\tilde{q}_{k+1}\ominus(\tilde{q}_{k+2}\ominus a_{k})$ (Finite-Precision
Clenshaw) $\displaystyle\tilde{\tilde{p}}$
$\displaystyle=\tfrac{1}{2}\odot((a_{0}\oplus q_{0})\ominus q_{2})$

Then Proposition 9.1 outputs $p(x)$ up to $50\varepsilon\lVert
p\rVert_{\mathrm{sup}}$ error141414 We did not attempt to optimize the
constants for this analysis. , provided that $\mu>0$ satisfies the following
three criterion.

  1. (a)

$\mu\varepsilon\leqslant\frac{1}{50(d+2)^{2}}$;

  2. (b)

$\mu\sum_{i=0}^{d}\MT@delim@RIGHT\abs@star{a_{i}}\leqslant\lVert
p\rVert_{\mathrm{sup}}$;

  3. (c)

$\mu\MT@delim@RIGHT\abs@star{q_{k}}=\mu\MT@delim@RIGHT\abs@nostar{\sum_{i=k}^{d}a_{i}U_{i-k}(x)}\leqslant\frac{1}{d}\lVert
p\rVert_{\mathrm{sup}}$ for all $k\in\\{0,\ldots,d\\}$.

This analysis shows that arithmetic operations incurring $\mu\varepsilon$
error result in computing $p(x)$ to $\varepsilon$ error. In particular, the
stability of the scalar Clenshaw recurrence comes down to understanding how
small we can take $\mu$. Note that if we ignored coefficient sign,
$\MT@delim@RIGHT\abs@nostar{\sum_{i=k}^{d}a_{i}U_{i-k}(x)}\leqslant\MT@delim@RIGHT\abs@nostar{\sum_{i=k}^{d}\MT@delim@RIGHT\abs@star{a_{i}}U_{i-k}(x)}=\sum_{i=k}^{d}(i-k+1)\MT@delim@RIGHT\abs@star{a_{i}}$,
this would require setting $\mu=\Theta(1/d^{3})$. We show in Section 9.2 that
we can set $\mu=\Theta((d^{2}\log(d))^{-1})$ for all $x\in[-1,1]$ and
polynomials $p$.

######  Lemma 9.2.

In Proposition 9.1, it suffices to take $\mu=\Theta((d^{2}\log(d))^{-1})$.

###### Proof of Proposition 9.1.

We will expand out these finite precision arithmetic to get error intervals
for each iteration.

$\tilde{q}_{d+1}=\tilde{q}_{d+2}=0,$ (37)

and

$\begin{split}\tilde{q}_{k}&=(2\odot
x)\odot\tilde{q}_{k+1}\ominus(\tilde{q}_{k+2}\ominus a_{k})\\\ &=(2x\pm
2\mu\varepsilon\MT@delim@RIGHT\abs@star{x})\odot\tilde{q}_{k+1}\ominus(\tilde{q}_{k+2}-a_{k}\pm\mu\varepsilon(\MT@delim@RIGHT\abs@star{\tilde{q}_{k+2}}+\MT@delim@RIGHT\abs@star{a_{k}}))\\\
&=((2x\tilde{q}_{k+1}\pm
2\mu\varepsilon\MT@delim@RIGHT\abs@star{x}\tilde{q}_{k+1})\pm\mu\varepsilon\MT@delim@RIGHT\abs@star{(2x\pm
2\mu\varepsilon\MT@delim@RIGHT\abs@star{x})\tilde{q}_{k+1}})\ominus(\tilde{q}_{k+2}-a_{k}\pm\mu\varepsilon(\MT@delim@RIGHT\abs@star{\tilde{q}_{k+2}}+\MT@delim@RIGHT\abs@star{a_{k}}))\\\
&\in(2x\tilde{q}_{k+1}\pm(2\mu\varepsilon+\mu^{2}\varepsilon^{2})2\MT@delim@RIGHT\abs@star{x\tilde{q}_{k+1}})\ominus(\tilde{q}_{k+2}-a_{k}\pm\mu\varepsilon(\MT@delim@RIGHT\abs@star{\tilde{q}_{k+2}}+\MT@delim@RIGHT\abs@star{a_{k}}))\\\
&\in(2x\tilde{q}_{k+1}\pm
6\mu\varepsilon\MT@delim@RIGHT\abs@star{x\tilde{q}_{k+1}})\ominus(\tilde{q}_{k+2}-a_{k}\pm\mu\varepsilon(\MT@delim@RIGHT\abs@star{\tilde{q}_{k+2}}+\MT@delim@RIGHT\abs@star{a_{k}}))\\\
&=2x\tilde{q}_{k+1}-\tilde{q}_{k+2}+a_{k}\pm\mu\varepsilon(6\MT@delim@RIGHT\abs@star{x\tilde{q}_{k+1}}+\MT@delim@RIGHT\abs@star{\tilde{q}_{k+2}}+\MT@delim@RIGHT\abs@star{a_{k}})\\\
&\qquad+\mu\varepsilon\MT@delim@RIGHT\abs@star{2x\tilde{q}_{k+1}\pm
6\mu\varepsilon\MT@delim@RIGHT\abs@star{x\tilde{q}_{k+1}}}+\mu\varepsilon\MT@delim@RIGHT\abs@star{\tilde{q}_{k+2}-a_{k}\pm\mu\varepsilon(\MT@delim@RIGHT\abs@star{\tilde{q}_{k+2}}+\MT@delim@RIGHT\abs@star{a_{k}})}\\\
&\in
2x\tilde{q}_{k+1}-\tilde{q}_{k+2}+a_{k}\pm\mu\varepsilon(14\MT@delim@RIGHT\abs@star{x\tilde{q}_{k+1}}+3\MT@delim@RIGHT\abs@star{\tilde{q}_{k+2}}+3\MT@delim@RIGHT\abs@star{a_{k}}),\end{split}$
(38)

and,

$\begin{split}\tilde{\tilde{p}}&=\tfrac{1}{2}\odot((a_{0}\oplus q_{0})\ominus
q_{2})\\\
&=\tfrac{1}{2}\odot((a_{0}+q_{0}\pm\mu\varepsilon(\MT@delim@RIGHT\abs@star{a_{0}}+\MT@delim@RIGHT\abs@star{q_{0}}))\ominus
q_{2})\\\
&=\tfrac{1}{2}\odot((a_{0}+q_{0}-q_{2}\pm\mu\varepsilon(\MT@delim@RIGHT\abs@star{a_{0}}+\MT@delim@RIGHT\abs@star{q_{0}}))\pm\mu\varepsilon(\MT@delim@RIGHT\abs@star{a_{0}+q_{0}\pm\mu\varepsilon(\MT@delim@RIGHT\abs@star{a_{0}}+\MT@delim@RIGHT\abs@star{q_{0}})}+\MT@delim@RIGHT\abs@star{q_{2}}))\\\
&\in\tfrac{1}{2}\odot(a_{0}+q_{0}-q_{2}\pm\mu\varepsilon(3\MT@delim@RIGHT\abs@star{a_{0}}+3\MT@delim@RIGHT\abs@star{q_{0}}+\MT@delim@RIGHT\abs@star{q_{2}}))\\\
&=\tfrac{1}{2}(a_{0}+q_{0}-q_{2}\pm\mu\varepsilon(3\MT@delim@RIGHT\abs@star{a_{0}}+3\MT@delim@RIGHT\abs@star{q_{0}}+\MT@delim@RIGHT\abs@star{q_{2}}))\pm\mu\varepsilon\tfrac{1}{2}\MT@delim@RIGHT\abs@star{a_{0}+q_{0}-q_{2}\pm\mu\varepsilon(3\MT@delim@RIGHT\abs@star{a_{0}}+3\MT@delim@RIGHT\abs@star{q_{0}}+\MT@delim@RIGHT\abs@star{q_{2}})}\\\
&\in\tfrac{1}{2}(a_{0}+q_{0}-q_{2})\pm\tfrac{1}{2}\mu\varepsilon(7\MT@delim@RIGHT\abs@star{a_{0}}+7\MT@delim@RIGHT\abs@star{q_{0}}+3\MT@delim@RIGHT\abs@star{q_{2}}).\end{split}$
(39)

To summarize, we have

$\displaystyle\tilde{q}_{d+1}$ $\displaystyle=\tilde{q}_{d+2}=0$
$\displaystyle\tilde{q}_{k}$
$\displaystyle=2x\tilde{q}_{k+1}-\tilde{q}_{k+2}+a_{k}+\delta_{k},$
$\displaystyle\text{where }\MT@delim@RIGHT\abs@star{\delta_{k}}$
$\displaystyle\leqslant\mu\varepsilon(14\MT@delim@RIGHT\abs@star{x\tilde{q}_{k+1}}+3\MT@delim@RIGHT\abs@star{\tilde{q}_{k+2}}+3\MT@delim@RIGHT\abs@star{a_{k}})$
(40) $\displaystyle\tilde{\tilde{p}}$
$\displaystyle=\tfrac{1}{2}(a_{0}+q_{0}-q_{2})+\delta,$
$\displaystyle\text{where }\MT@delim@RIGHT\abs@star{\delta}$
$\displaystyle\leqslant\tfrac{1}{2}\mu\varepsilon(7\MT@delim@RIGHT\abs@star{a_{0}}+7\MT@delim@RIGHT\abs@star{q_{0}}+3\MT@delim@RIGHT\abs@star{q_{2}})$
(41)

By Lemma 8.1, this recurrence satisfies

$\displaystyle\tilde{q}_{k}$
$\displaystyle=\sum_{i=k}^{d}U_{i-k}(x)(a_{i}+\delta_{i})$ $\displaystyle
q_{k}-\tilde{q}_{k}$ $\displaystyle=\sum_{i=k}^{d}U_{i-k}(x)\delta_{i}$
$\displaystyle q-\tilde{q}$
$\displaystyle=\delta+\frac{1}{2}\Big{(}\sum_{i=0}^{d}U_{i}(x)\delta_{i}-\sum_{i=2}^{d}U_{i-2}(x)\delta_{i}\Big{)}$
$\displaystyle=\delta+\frac{1}{2}\delta_{0}+\sum_{i=1}^{d}T_{i}(x)\delta_{i}$
$\displaystyle\MT@delim@RIGHT\abs@star{q-\tilde{q}}$
$\displaystyle\leqslant\MT@delim@RIGHT\abs@star{\delta}+\frac{1}{2}\MT@delim@RIGHT\abs@star{\delta_{0}}+\sum_{i=1}^{d}\MT@delim@RIGHT\abs@star{T_{i}(x)\delta_{i}}\leqslant\MT@delim@RIGHT\abs@star{\delta}+\sum_{i=0}^{d}\MT@delim<EMAIL_ADDRESS>(42) This analysis so far has been fully standard. Let’s continue bounding.
$\displaystyle\leqslant\mu\varepsilon\MT@delim@RIGHT\parens@nostar[\Big{]}{\tfrac{7}{2}\MT@delim@RIGHT\abs@star{a_{0}}+\tfrac{7}{2}\MT@delim@RIGHT\abs@star{q_{0}}+\tfrac{3}{2}\MT@delim@RIGHT\abs@star{q_{2}}+\sum_{i=0}^{d}(14\MT@delim@RIGHT\abs@star{x\tilde{q}_{i+1}}+3\MT@delim@RIGHT\abs@star{\tilde{q}_{i+2}}+3\MT@delim@RIGHT\abs@star{a_{i}})}$
$\displaystyle\leqslant\mu\varepsilon\sum_{i=0}^{d}(20\MT@delim@RIGHT\abs@star{\tilde{q}_{i}}+10\MT@delim@RIGHT\abs@star{a_{i}}).$
(43)

Now, we will bound all of the $\delta_{k}$’s. Combining previous facts, we
have

$\displaystyle\MT@delim@RIGHT\abs@star{\tilde{q}_{k}}$
$\displaystyle=\MT@delim@RIGHT\abs@star{\sum_{i=k}^{d}U_{i-k}(x)(a_{i}+\delta_{i})}$
$\displaystyle\leqslant\MT@delim@RIGHT\abs@nostar[\Big{]}{\sum_{i=k}^{d}U_{i-k}(x)a_{i}}+\sum_{i=k}^{d}\MT@delim@RIGHT\abs@nostar[\Big{]}{U_{i-k}(x)\delta_{i}}$
$\displaystyle\leqslant\MT@delim@RIGHT\abs@nostar[\Big{]}{\sum_{i=k}^{d}U_{i-k}(x)a_{i}}+\sum_{i=k}^{d}(i-k+1)\MT@delim@RIGHT\abs@star{\delta_{i}}$
$\displaystyle\leqslant\MT@delim@RIGHT\abs@nostar[\Big{]}{\sum_{i=k}^{d}U_{i-k}(x)a_{i}}+\mu\varepsilon\sum_{i=k}^{d}(i-k+1)(14\MT@delim@RIGHT\abs@star{\tilde{q}_{i+1}}+3\MT@delim@RIGHT\abs@star{\tilde{q}_{i+2}}+3\MT@delim@RIGHT\abs@star{a_{i}})$
$\displaystyle\leqslant\MT@delim@RIGHT\parens@nostar[\Big{]}{\frac{1}{\mu
d}+3\mu\varepsilon\frac{d-k+1}{\mu}}\lVert
p\rVert_{\mathrm{sup}}+\mu\varepsilon\sum_{i=k}^{d}(i-k+1)(14\MT@delim@RIGHT\abs@star{\tilde{q}_{i+1}}+3\MT@delim@RIGHT\abs@star{\tilde{q}_{i+2}})$
$\displaystyle\leqslant\frac{1.5}{\mu d}\lVert
p\rVert_{\mathrm{sup}}+\mu\varepsilon\sum_{i=k}^{d}(i-k+1)(14\MT@delim@RIGHT\abs@star{\tilde{q}_{i+1}}+3\MT@delim@RIGHT\abs@star{\tilde{q}_{i+2}})$

Note that $\MT@delim@RIGHT\abs@star{\tilde{q}_{k}}\leqslant c_{k}$, where

$\displaystyle c_{d}$ $\displaystyle=0;$ $\displaystyle c_{k}$
$\displaystyle=\frac{1.5}{\mu d}\lVert
p\rVert_{\mathrm{sup}}+\mu\varepsilon\sum_{i=k}^{d}(i-k+1)(14c_{i+1}+3c_{i+2})$
(44) Solving this recurrence, we have that $c_{k}\leqslant\frac{2}{\mu
d}\lVert p\rVert_{\mathrm{sup}}$, since by strong induction, $\displaystyle
c_{k}$ $\displaystyle\leqslant\Big{(}\frac{1.5}{\mu
d}+\mu\varepsilon\sum_{i=k}^{d}(i-k+1)17\frac{2}{\mu d}\Big{)}\lVert
p\rVert_{\mathrm{sup}}$ $\displaystyle=\Big{(}\frac{1.5}{\mu
d}+17\mu\varepsilon\frac{1}{\mu d}(d-k+1)(d-k+2)\Big{)}\lVert
p\rVert_{\mathrm{sup}}\leqslant\frac{2}{\mu d}\lVert p\rVert_{\mathrm{sup}}$
(45)

Returning to Equation 43:

$\displaystyle\MT@delim@RIGHT\abs@star{q-\tilde{q}}$
$\displaystyle\leqslant\mu\varepsilon\sum_{i=0}^{d}(20\MT@delim@RIGHT\abs@star{\tilde{q}_{i}}+10\MT@delim@RIGHT\abs@star{a_{i}})\leqslant\mu\varepsilon\sum_{i=0}^{d}(20c_{i}+10\MT@delim@RIGHT\abs@star{a_{i}})$
(46) $\displaystyle\leqslant 40\varepsilon\lVert
p\rVert_{\mathrm{sup}}+10\mu\varepsilon\sum_{i=0}^{d}\MT@delim@RIGHT\abs@star{a_{i}}\leqslant
50\varepsilon\lVert p\rVert_{\mathrm{sup}}\qed$ (47)

###  9.2 Bounding the iterates of the Clenshaw recurrence

The goal of this section is to prove Lemma 9.2. In particular, we wish to show
that for $\mu=\Theta((d^{2}\log(d))^{-1})$, the following criteria hold:

  1. (a)

$\mu\varepsilon\leqslant\frac{1}{50(d+2)^{2}}$;

  2. (b)

$\mu\sum_{i=0}^{d}\MT@delim@RIGHT\abs@star{a_{i}}\leqslant\lVert
p\rVert_{\mathrm{sup}}$;

  3. (c)

$\mu\MT@delim@RIGHT\abs@star{q_{k}}=\mu\MT@delim@RIGHT\abs@nostar{\sum_{i=k}^{d}a_{i}U_{i-k}(x)}\leqslant\frac{1}{d}\lVert
p\rVert_{\mathrm{sup}}$ for all $k\in\\{0,\ldots,d\\}$.

For this choice of $\mu$, (a) is clearly satisfied, and since
$\MT@delim@RIGHT\abs@star{a_{i}}\leqslant 2\lVert p\rVert_{\mathrm{sup}}$
(Lemma 5.2), $\mu\sum_{i=0}^{d}\MT@delim@RIGHT\abs@star{a_{i}}\leqslant
2(d+1)\lVert p\rVert_{\mathrm{sup}}\leqslant\lVert p\rVert_{\mathrm{sup}}$, so
(b) is satisfied. In fact, both of these criterion are satisfied for
$\mu=O(1/d)$, provided $\varepsilon$ is sufficiently small.

Showing (c) requires bounding
$\lVert\sum_{\ell=k}^{d}a_{\ell}U_{\ell-k}(x)\rVert_{\mathrm{sup}}$ for all
$k\in[d]$. These expressions are also the iterates of the Clenshaw algorithm
(Lemma 8.1), so we are in fact trying to show that in the process of our
algorithm we never produce a value that’s much larger than the final value.
From testing computationally, we believe that the following holds true.

######  Conjecture 9.3.

Let $p(x)$ be a degree-$d$ polynomial with Chebyshev expansion
$p(x)=\sum_{\ell=0}^{d}a_{\ell}T_{\ell}(x)$. Then, for all $k$ from $0$ to
$d$,

$\displaystyle\Big{\lVert}\sum_{\ell=k}^{d}a_{\ell}U_{\ell-k}(x)\Big{\rVert}_{\mathrm{sup}}\leqslant(d-k+1)\lVert
p\rVert_{\mathrm{sup}},$

maximized for the Chebyshev polynomial $p(x)=T_{d}(x)$.

9.3 would imply that it suffices to take $\mu=\Theta(1/d^{2})$. We prove it up
to a log factor.

######  Theorem 9.4.

For a degree-$d$ polynomial $p$, consider the degree-$(d-k)$ polynomial
$q_{k}(x)=\sum_{\ell=k}^{d}a_{\ell}U_{\ell-k}(x)$. Then

$\displaystyle\lVert
q_{k}\rVert_{\mathrm{sup}}\leqslant(d-k+1)\Big{(}16+\frac{16}{\pi^{2}}\log(d)\Big{)}\lVert
p\rVert_{\mathrm{sup}}.$

###### Proof.

We proceed by carefully bounding the Chebyshev coefficients of $q_{k}$, which
turn out to be arithmetic progressions of the $a_{k}$’s which we bounded in
Section 7.

$\displaystyle q_{k}(x)$ $\displaystyle=\sum_{i}a_{i}U_{i-k}(x)$
$\displaystyle=\sum_{i}\sum_{j\geqslant
0}a_{i}T_{i-k-2j}(x)(1+\MT@delim@RIGHT\iver@star{i-k-2j\neq 0})$
$\displaystyle=\sum_{i}\sum_{j\geqslant
0}a_{i+k+2j}T_{i}(x)(1+\MT@delim@RIGHT\iver@star{i\neq 0})$
$\displaystyle=\sum_{i}T_{i}(x)(1+\MT@delim@RIGHT\iver@star{i\neq
0})\sum_{j\geqslant 0}a_{i+k+2j}$
$\displaystyle\MT@delim@RIGHT\abs@star{q_{k}(x)}$
$\displaystyle\leqslant\sum_{i}\MT@delim@RIGHT\iver@star{i\geqslant
0}(1+\MT@delim@RIGHT\iver@star{i\neq
0})\MT@delim@RIGHT\abs@nostar[\Big{]}{\sum_{j\geqslant 0}a_{i+k+2j}}$
$\displaystyle\leqslant
2\sum_{i=0}^{d-k}\MT@delim@RIGHT\abs@nostar[\Big{]}{\sum_{j\geqslant
0}a_{i+k+2j}}$ $\displaystyle\leqslant
4\sum_{i=0}^{d-k}\Big{(}4+\frac{4}{\pi^{2}}\log(i+k-1)\Big{)}\lVert
p\rVert_{\mathrm{sup}}$ by 7.3
$\displaystyle\leqslant(d-k+1)\Big{(}16+\frac{16}{\pi^{2}}\log(d)\Big{)}\lVert
p\rVert_{\mathrm{sup}}.\qed$

######  Remark 9.5.

We spent some time trying to prove 9.3, since its form is tantalizingly close
to that of the Markov brothers’ inequality [Sch41]
$\lVert\frac{d}{dx}p(x)\rVert_{\mathrm{sup}}=\lVert\sum_{\ell=0}^{d}a_{\ell}\ell
U_{\ell-1}(x)\rVert_{\mathrm{sup}}\leqslant d^{2}\lVert
p(x)\rVert_{\mathrm{sup}}$, except with the linear differential operator
$\frac{d}{dx}:T_{\ell}\mapsto\ell U_{\ell-1}$ replaced with the linear
operator $T_{\ell}\mapsto U_{\ell-k}$. However, calculations suggest that the
variational characterization of $\max_{\lVert
p\rVert_{\mathrm{sup}}=1}\MT@delim@RIGHT\abs@nostar{\frac{d}{dx}p(x)}$
underlying proofs of the Markov brothers’ inequality [Sha04] does not hold
here, and from our shallow understanding of these proofs, it seems that they
strongly use properties of the derivative.

##  10 Computing matrix polynomials

Our goal is to prove the following theorem:

######  Theorem 10.1.

Suppose we are given sampling and query access to $A\in\mathbb{C}^{m\times n}$
and $b\in\mathbb{C}^{n}$, and suppose we are given an even or odd degree-$d$
polynomial. Then we can output a description of a vector $y\in\mathbbm{C}^{n}$
such that
$\|y-p(A)b\|\leqslant\varepsilon\MT@delim@RIGHT\norm@star{p}_{\operatorname{Spec}(A)}$
with probability $\geqslant 0.9$ in time

$\displaystyle
O\MT@delim@RIGHT\parens@star{\min\MT@delim@RIGHT\braces@star{\operatorname{nnz}(A),\frac{d^{8}\|A\|_{F}^{4}}{\varepsilon^{4}}\log^{8}(d)\log^{2}(n)}+\frac{d^{9}\|A\|_{F}^{4}}{\varepsilon^{2}}\log^{4}(d)\log(n)}$

We can access the output description in the following way:

  1. (a)

Compute entries of $y$ in $O\MT@delim@RIGHT\Paren@star{\frac{d^{4}\lVert
A\rVert_{F}^{2}}{\varepsilon^{2}}\log^{4}(d)\log(n)}$ time;

  2. (b)

Sample $i\in[n]$ with probability
$\frac{\MT@delim@RIGHT\abs@star{y_{i}}^{2}}{\MT@delim@RIGHT\norm@star{y}^{2}}$
in $O\MT@delim@RIGHT\Paren@star{\frac{d^{6}\lVert
A\rVert_{F}^{4}}{(\varepsilon\MT@delim@RIGHT\norm@star{y})^{2}}\log^{8}(d)\log(n)}$
time;

  3. (c)

Estimate $\MT@delim@RIGHT\norm@star{y}^{2}$ to $\nu$ relative error in
$O\MT@delim@RIGHT\Paren@star{\frac{d^{6}\lVert
A\rVert_{F}^{4}}{(\nu\varepsilon\MT@delim@RIGHT\norm@star{y})^{2}}\log^{8}(d)\log(n)}$
time.

We now extend the error analysis in the scalar setting from Section 9 to the
matrix setting. For simplicity, we treat the odd case, and in Section 10.2,
describe what changes to prove the result for even polynomials.

###  10.1 Computing odd matrix polynomials

The main theorem we obtain is as follows:

######  Theorem 10.2.

Suppose we are given sampling and query access to $A\in\mathbb{C}^{m\times n}$
and $b\in\mathbb{C}^{n}$ with
$\MT@delim@RIGHT\norm@star{A},\MT@delim@RIGHT\norm@star{b}\leqslant 1$, and
suppose we are given a $(2d+1)$-degree odd polynomial, written in its
Chebyshev coefficients as

$\displaystyle p(x)=\sum_{i=0}^{d}a_{2i+1}T_{2i+1}(x).$

Then we can output a vector $x\in\mathbbm{C}^{n}$ such that
$\|Ax-p(A)b\|\leqslant\varepsilon\MT@delim@RIGHT\norm@star{p}_{\operatorname{Spec}(A)}$
with probability $\geqslant 0.9$ in time

$\displaystyle
O\MT@delim@RIGHT\parens@star{\min\MT@delim@RIGHT\braces@star{\operatorname{nnz}(A),\frac{\|A\|_{F}^{4}}{(\mu\varepsilon)^{4}}\log^{2}(n)}+\frac{d^{5}\|A\|_{F}^{4}}{(\mu\varepsilon)^{2}}\log(n)},$

assuming $\mu$ satisfies the following bounds:

  1. (a)

$\mu\varepsilon\leqslant\frac{1}{100d^{2}}$;

  2. (b)

$\mu\sum_{i=0}^{d}\MT@delim@RIGHT\abs@star{a_{2i+1}}\leqslant 1$;

  3. (c)

$\mu\lVert\sum_{i=k}^{d}a_{2i+1}U_{i-k}(T_{2}(x))\rVert_{\mathrm{sup}}\leqslant\frac{d-k+1}{d^{2}}$
for all $0\leqslant k\leqslant d$;

The output description has the additional properties

$\displaystyle\sum_{i}\MT@delim@RIGHT\norm@star{A_{i,*}}^{2}\MT@delim@RIGHT\abs@star{x_{i}}^{2}\lesssim\frac{\varepsilon^{2}}{\log(n)d^{2}}\qquad\lVert
x\rVert_{0}\lesssim\frac{\lVert
A\rVert_{F}^{2}}{(\mu\varepsilon)^{2}}\log(n),$

so that by Corollary 5.11, we can

  1. (a)

Compute entries of $Ax$ in $\lVert
x\rVert_{0}=O\MT@delim@RIGHT\Paren@star{\frac{\lVert
A\rVert_{F}^{2}}{(\mu\varepsilon)^{2}}\log(n)}$ time;

  2. (b)

Sample $i\in[n]$ with probability
$\frac{\MT@delim@RIGHT\abs@star{y_{i}}^{2}}{\MT@delim@RIGHT\norm@star{y}^{2}}$
in $O\MT@delim@RIGHT\Paren@star{\frac{\lVert
A\rVert_{F}}{\varepsilon^{2}\mu^{2}d^{2}\MT@delim@RIGHT\norm@star{y}^{2}}}$
time;

  3. (c)

Estimate $\MT@delim@RIGHT\norm@star{y}^{2}$ to $\nu$ relative error in
$O\MT@delim@RIGHT\Paren@star{\frac{\lVert
A\rVert_{F}^{4}\log(n)}{(\nu\mu^{2}\varepsilon
d\MT@delim@RIGHT\norm@star{y})^{2}}\log(n)}$ time.

The criterion for what $\mu$ need to be are somewhat non-trivial; essentially,
these bounds are what’s necessary for Section 8.2 to be numerically stable
when computing $p(x)/x$. This is necessary because we primarily work in the
“dual” space, maintaining our Clenshaw iterate $u_{k}$ as $Av_{k}$ where
$v_{k}$ is a sparse vector. For any bounded polynomial, though, we can always
take $\mu$ to be $\Omega((d^{2}\log^{2}(d))^{-1})$.

######  Proposition 10.3.

In Theorem 10.2, we can always take $1/\mu\lesssim
d^{2}\log^{2}(d)\MT@delim@RIGHT\norm@star{p}_{\operatorname{Spec}(A)}$.

We are now ready to dive into the proof of Theorem 10.2.

######  Algorithm 10.4 (Odd singular value transformation).

Input: A matrix $A\in\mathbb{C}^{m\times n}$, vector $b\in\mathbb{C}^{n}$,
$\varepsilon,\delta,\mu>0$, a degree $2d+1$ polynomial $p$, and the
corresponding coefficients in the Chebyshev basis, i.e.
$p(x)=\sum_{i=0}^{d}a_{2i+1}T_{2i+1}(x)$. Preprocessing sketches: Let
$s,t=\Theta\MT@delim@RIGHT\Paren@star{\frac{\lVert
A\rVert_{F}^{2}}{(\mu\varepsilon)^{2}}\log(\frac{n}{\delta})}$. P1. If
$\operatorname{SQ}(A^{\dagger})$ and $\operatorname{SQ}(b)$ are not given,
compute data structures to simulate them in $O(1)$ time; P2. Sample
$S\in\mathbbm{C}^{n\times s}$ to be
${\textsc{AMP}}\MT@delim@RIGHT\Paren@star{A,b,\\{\tfrac{1}{2}(\frac{\|A_{*,i}\|^{2}}{\lVert
A\rVert_{F}^{2}}+\frac{\MT@delim@RIGHT\abs@star{b_{i}}^{2}}{\MT@delim@RIGHT\norm@star{b}^{2}})\\}}$
(Theorem 6.5) with sketch size $s$; P3. Sample $T\in\mathbbm{C}^{t\times m}$
to be ${\textsc{AMP}}\MT@delim@RIGHT\Paren@star{S^{\dagger}A^{\dagger},AS}$
(Theorem 6.5) with sketch size $t$; P4. Compute a data structure that can
respond to $\operatorname{SQ}(TAS)$ queries in $O(1)$ time; Clenshaw
iteration: Let $r=\Theta(d^{4}\lVert
A\rVert_{F}^{2}(s+t))=\Theta\MT@delim@RIGHT\Paren@star{\frac{d^{4}\lVert
A\rVert_{F}^{4}}{(\mu\varepsilon)^{2}}\log\frac{n}{\delta}}$. Then starting
with $v_{d+1}=v_{d+2}=\vec{0}^{s}$ and going until $v_{0}$, I1. Let
$B^{(k)}={\textsc{best}}(TAS)$ and
$B_{\dagger}^{(k)}={\textsc{best}}((TAS)^{\dagger})$ with parameter $r$; I2.
Compute $v_{k}=2(2B_{\dagger}^{(k)}B^{(k)}-I)v_{k+1}-v_{k+2}+a_{2k+1}Sb$.
Output: Output $x=\tfrac{1}{2}S(v_{0}-v_{1})$ that satisfies
$\MT@delim@RIGHT\Norm@star{Ax-p(A)b}\leqslant\varepsilon$.

#### Preprocessing sketches.

Given a matrix $A\in C^{m\times n}$ and a vector $b\in\mathbbm{C}^{n}$, the
pre-processing phase of Algorithm 10.4 can be performed in
$O(\operatorname{nnz}(A))$ time. If $O(Q)$-time $\operatorname{SQ}(A)$ and
$\operatorname{SQ}(b)$ are already given, then this can be performed in
$O(Qst)$ time.

The way to do this is straightforward. Finally, we query all of $TAS$ in
$O(\min(\operatorname{nnz}(A),st))$ time. Using this, we can prepare the data
structure for $\operatorname{SQ}(TAS)$.