Quantum machine learning (QML) is a field of study with a rapidly growing number of proposals for how quantum computers could significantly speed up machine learning tasks [DW20, CHI + 18]. If any of these proposals yield substantial practical speedups, it could be the killer application motivating the development of scalable quantum computers [Pre18]. At first glance, many applications of QML seem to admit exponential speedups. However, these exponential speedups are less likely to manifest in practice compared to, say, Shor's algorithm for factoring [Sho97], because unlike their classical counterparts, QML algorithms must make strong input assumptions and learn relatively little from their output [Aar15]. These caveats arise because both loading input data into a quantum computer and extracting amplitude data from an output quantum state are hard in their most generic forms.
A recent line of research analyzes the speedups of QML algorithms by developing classical counterparts that carefully exploit these restrictive input and output assumptions. This began with a breakthrough 2018 paper by Tang [Tan19] showing that the quantum recommendation systems algorithm [KP17], previously believed to be one of the strongest candidates for a practical exponential speedup in QML, does not give an exponential speedup. Specifically, Tang describes a "dequantized" algorithm that solves the same problem as the quantum algorithm and only suffers a polynomial slowdown. Tang's algorithm crucially exploits the input data structure assumed by the quantum algorithm, which is used for efficiently preparing states. Subsequent work relies on similar techniques to dequantize a wide range of QML algorithms, including those for principal component analysis and supervised clustering [Tan21], low-rank linear system solving [CGL + 20], low-rank semidefinite program solving [CLLW20], support vector machines [DBH21], nonnegative matrix factorization [CLS + 19], and minimal conical hull [DHLT20]. These results show that the advertised exponential speedups of many QML algorithms disappear when compared to classical algorithms with input assumptions analogous to the state preparation assumptions of the quantum algorithms, drastically changing our understanding of the landscape of potential QML algorithm speedups.
A recent line of work in quantum algorithms has worked to unify many quantum algorithms ranging from quantum walks to QML, under a quantum linear algebra framework called quantum singular value transformation (QSVT) [LC17,CGJ19,GSLW19]. Since this framework captures essentially all known linear algebraic QML techniques [MRTC21], including all prior dequantized QML algorithms (up to minor technical details), a natural question is whether this framework can be dequantized. One cannot hope to dequantize all of QSVT, because with sparse block-encodings of input data, QSVT can simulate Harrow, Hassidim, and Lloyd's pioneering poly-logarithmic time algorithm (HHL) for the BQP-complete problem of sampling from the solution of a sparse system of linear equations [HHL09]. However, one could hope to dequantize QSVT provided that the input data comes in the state preparation data structure used commonly for quantum linear algebra. This data structure allows for efficient QML when the input is close-to-low-rank, but as dequantized algorithms show, it also gives significant power to classical algorithms. Prior work [Tan21, CGL + 20] has made similar speculations that these techniques could feasibly dequantize wide swathes of quantum linear algebra. In this work, we give evidence for these hopes by presenting a classical analogue of the QSVT framework and applying it to dequantize QML algorithms.
We describe a simple framework for quantum-inspired classical algorithms with wide applicability, grasping the capabilities and limitations of these techniques. We use this framework to dequantize many quantum linear algebra algorithms. We also prove QSVT-like extensibility properties of our framework, giving evidence that it can dequantize any QSVT algorithms in the QRAM input model.
Input model: Oversampling and query access. Our framework assumes a specific input model called oversampling and query access, which can be thought of as a classical analogue to quantum state preparation assumptions, i.e., the ability to prepare a quantum state |v⟩ proportional to some input vector v. Our conceptual contribution is to define this generalization of sampling and query access, because it has better closure properties.
We have sampling and query access to a vector v ∈ C n , denoted SQ(v), if we can efficiently make the following kinds of queries (Definition 3.2): (1) given an index i ∈ [n], output the corresponding entry v(i); (2) sample an index j ∈ [n] with probability |v(j)| 2 /∥v∥ 2 ; and (3) output the vector's ℓ 2 -norm ∥v∥. We have sampling and query access to a matrix A ∈ C m×n , denoted SQ(A), if we have SQ(A(i, •)) for all rows of A, A(i, •), and also SQ(a) for a the vector of row norms (i.e., a(i) := ∥A(i, •)∥). We have ϕ-oversampling and query access to a vector v, denoted SQ ϕ (v), if (1) we can query for entries of v and (2) we have sampling and query access to an "entry-wise upper bound" vector v ˜satisfying ∥v ˜∥2 = ϕ∥v∥ 2 and |v ˜(i)| ≥ |v(i)| for all indices i; the definition for a matrix is analogous (Definition 3.4).
The parameter ϕ should be seen as a form of overhead that comes out in the runtime of algorithms: through rejection sampling, SQ ϕ (v) can do approximate versions of all the queries of SQ(v) with a factor ϕ of overhead (Lemma 3.5). In this paper, we most often think of ϕ as being independent of input size.
To motivate this definition, we make the following observations about this input model. First, as far as we know, if input data is given classically,foot_0 classical algorithms in the sampling and query model can be run whenever the corresponding algorithms in the quantum model can (Remark 3.11). For example, if input is loaded in the QRAM data structure, as commonly assumed in QML in order to satisfy state preparation assumptions [Pra14, CHI + 18], then we have log-time sampling and query access to it. So, a fast classical algorithm for a problem in this classical model implies lack of quantum speedup for the problem, at least in the usual settings explored in the QML literature. In particular, a polynomial-time classical algorithm in this model implies lack of exponential quantum speedup. Second, oversampling and query access has many similarities to the notion of quantum block-encodings in quantum singular value transformation [GSLW19]. The commonly used datastructures that enable oversampling and query access to a matrix A also enable implementing an efficient quantum circuit whose unitary is a block-encoding of A. Further, in both input models one can perform efficient matrix arithmetic.
Matrix arithmetic. The thrust of our main results is to demonstrate that oversampling and query access is approximately closed under arithmetic operations. We argue that the essential power of quantum-inspired algorithms lies in their ability to leverage sampling and query access of the input matrices to provide oversampling and query access to complex arithmetic expressions as output (possibly with some approximation error), without paying the (at least) linear time necessary to compute such expressions in conventional ways. While the proven closure properties are important from a complexity theoretic point of view, some of them don't explicitly come into play in the demonstrated applications of our framework for dequantizing QML, since they often come with undesirable polynomial overhead. This is in contrast with quantum block-encodings, which generally compose with minimal overhead.
We now list the closure properties that we show, along with the corresponding closure properties proven for block-encodings in [GSLW19]. For all of these, the query time for access to the output is just polynomial in the query times for access to the input, so in particular, these procedures run in time independent of input dimension. More specifically, we will compare what we call (sub)normalization "overhead" between the two, which is the value ϕ in the classical setting and what [GSLW19] denotes as α in the quantum setting. The two quantities are analogous, and roughly correspond to overheads in rejection sampling and post-selection, inducing a multiplicative factor in sampling times for both.
• Given access to a constant number of vectors v 1 , . . . , v τ , we have access to linear combinations ∑︁ τ t=1 λ t v t , and analogously with linear combinations of matrices (Lemmas 3.6 and 3.9). This is a classical analogue to the "linear combinations of unitaries" technique for block-encodings [GSLW19, Lemma 52]. In the quantum setting, there is less overhead, 2 allowing for efficient block encodings for linear combinations of arbitrarily many matrices in certain settings.
• Given access to two matrices A, B with Frobenius norm at most one, we have access to a matrix Z ε-close to the product A † B in Frobenius norm (Lemma 4.6 and Remark 4.7). In the quantum setting, closure of block-encodings under products is almost immediate [GSLW19, Lemma 53] and is not approximate. In both cases the individual input overheads of A and B are multiplied. With the same overheads one can also form Kronecker products A ⊗ B exactly-this is immediate both in the classical and quantum case [CVB20]. In particular, given access to two vectors u and v, we have access to their outer product uv † (Lemma 3.8).
• Given access to a matrix A with Frobenius norm at most one and a Lipschitz function f , we have access to a matrix Z ε-close to f (A † A) in Frobenius norm (Theorem 5.1) 3 . In the quantum setting, block-encodings are closed under even and odd polynomial singular value transformations [GSLW19, Lemmas 8, 10] without approximation, provided the polynomial is low-degree and bounded. This block-encoding closure property can be viewed as a corollary of the above two properties, but can also be achieved directly with some more efficient technical machinery.
An even polynomial singular value transformation of A is precisely f (A † A) for f a low-degree polynomial (which means f is Lipschitz), and odd polynomials can be decomposed into a product of an even polynomial with A, so our closure property is as strong as the quantum one. The full details are derived in Section 6.1.
To summarize, every arithmetic operation of matrices with block-encodings in [GSLW19] (with the possible exception of long linear combinations) can be mimicked by matrices with oversampling and query access, up to Frobenius norm error, provided that an input matrix in a block-encoding corresponds to having 4 SQ(A) and SQ(A † ). The "linear combinations of vectors" and "outer products" classical closure properties have been used in prior work [Tan19, CGL + 20]. However, without our new definition of oversampling and query access, it was not clear that these algorithms could be chained indefinitely as we show with these closure properties.
2 In fact, already √ ϕ ≥ α/ ⃦ ⃦ ∑︁ τ t=1 λtvt ⃦ ⃦ which can be seen using that the root mean-square is at least the average. 3 For a Hermitian matrix H and a function f : R ↦ → C, f (H) denotes applying f to the eigenvalues of H. That is, f (H) := ∑︁ n i=1 f (λi)viv † i , for λi and vi the eigenvalues and eigenvectors of H. 4 We take some care here to distinguish whether we have oversampling and query access to A or A † . We don't need to: we show that having either one of them implies having the other, up to approximation (Remark 6.4). However, the accesses assumed in our closure properties are in some sense the most natural choices and require the least overhead.
Implications for quantum singular value transformation. Our results give compelling evidence that there is indeed no exponential speedup for QRAM-based QSVT, and show that oversampling and query access can be thought of as a classical analogue to block-encodings in the bounded Frobenius norm regime. Nevertheless, we do not rule out the possibility for large polynomial speedups, as the classical runtimes tend to have impractically large polynomial exponents. To elaborate more on this connection, we now recall the QSVT framework in more detail.
The QSVT framework of Gilyén, Su, Low, and Wiebe [GSLW19] assumes that the input matrix A is given by a block-encoding, which is a quantum circuit implementing a unitary transformation whose top-left block contains (up to scaling) A itself [LC17]. Given a block-encoding of A, one can apply it to a quantum state or form block-encodings of other expressions using the closure properties mentioned above. One can get a block-encoding of an input matrix A through various methods. If A is s-sparse with efficiently computable elements and ∥A∥ ≤ 1, then one can directly get a block-encoding of A/s [GSLW19, Lemma 48]. If A is in the QRAM data structure (used for efficient state preparation for QML algorithms [Pra14]), one can directly get a block-encoding of A/∥A∥ F [GSLW19, Lemma 50]. We will use the term QRAM-based QSVT to refer to the family of quantum algorithms possible in the QSVT framework when all input matrices & vectors are given in the QRAM data structure.
The normalization in QRAM-based QSVT means that it has an implicit dependence on the Frobenius norm ∥A∥ F . Since ∥A∥ F is also the key parameter in the complexity of our corresponding classical algorithms, this suggests that QRAM-based QSVT does not give inherent exponential quantum speedups (though, if input preparation/output analysis protocols have no classical analogues, they can act as a subroutine in an algorithm that does give an exponential quantum speedup). Our closure results confirm this: if input matrices and vectors are given in QRAM data structure, then on one (quantum) hand we can construct block-encodings of these matrices normalized to have Frobenius norm one and on the other (classical) hand we have sampling and query access to the input. The conclusion is that, up to some controllable approximation error, an algorithm using the block-encoding framework has a classical analogue in the oversampling and query access model. The classical algorithm's runtime is only polynomially slower than the corresponding quantum algorithm, except in the ε parameter. 5 One can argue similarly that there should be no exponential speedup for QSVT for block-encodings derived from (purifications of) density operators [GSLW19, Lemma 45] that come from some well-structured classical data (see Section 6.1). This stands in contrast to, for example, block-encodings that come from sparsity assumptions [GSLW19,Lemma 48], where the matrix in the block-encoding can have Frobenius norm as large as √ sn (where we take A to be n × n), and so the classical techniques cannot be applied without incurring dependence on n in the runtime.
We now illustrate the flavor of the algorithmic ideas underlying our main results, by showing why the "oversampling" input model is closed under approximate matrix products. Suppose we are given sampling and query access to two matrices A ∈ C m×n and B ∈ C m×p , and desire (over)sampling and query access to A † B. A † B is a sum of outer products of rows of A with rows of B (that is,
, so a natural idea is to use the outer product closure property to get access to each outer product individually, and then use the linear combination closure property to get access to their sum, which is A † B as desired. However, there are m terms in the sum, which is too large: we can't even compute entries of A † B in time independent of m. So, we use sampling to approximate this sum of m terms by a linear combination over far fewer terms, allowing us to get access to Z for Z ≈ A † B. This type of matrix product approximation is well-known in the classical literature [DKM06]. Given SQ(A), we can pull samples i 1 , . . . , i s according to the row norms of A, a distribution we will denote p (so p(i) = ∥A(i, •)∥ 2 /∥A∥ 2 F ). Consider Z := 1 s ∑︁ s k=1 1
Further, the variance of this estimator is small. In the following computation, we consider s = 1, because the variance for general s decreases as 1/s.
By Chebyshev's inequality, we can choose
)︁ to get that ∥Z -A † B∥ F < ε∥A∥ F ∥B∥ F with probability 0.99. Since Z is a linear combination of s outer products, this gives us oversampling and query access to Z as desired. In our applications we would keep Z as an outer product A ′ † B ′ for convenience. We have just sketched the proof of our key lemma: an approximate matrix product protocol.
Key lemma [DKM06] (informal version of Lemma 4.6). Suppose we are given SQ(X) ∈ C m×n and SQ(Y ) ∈ C m×p . Then we can find normalized submatrices of X and Y ,
We subsequently have O(s)-time SQ(X ′ ), SQ(X ′ † ), SQ(Y ′ ), SQ(Y ′ † ).
Prior quantum-inspired algorithms [Tan19, Tan21, CLW18, CLLW20] indirectly used this lemma by using [FKV04], which finds a low-rank approximation to the input matrix in the form of an approximate low-rank SVD and relies heavily on this lemma in the analysis.
One of our main results, mentioned earlier as our singular value transformation closure property, is that, given SQ(A) ∈ C m×n , in time independent of m and n, we can access an approximation of f (A † A) for a Lipschitz-function f that, without loss of generality, satisfies f (0) = 0 (Theorem 5.1). One could use [FKV04] to give a classical algorithm for SVT, but a more efficient approach is to directly apply the key lemma twice to get an approximate decomposition of f (A † A):
r×n is a subset of rows of the input matrix and U is a matrix with size independent of input dimension (R corresponds to the 'R' of the RUR decomposition, and f ¯(CC † ) ∈ C r×r corresponds to the 'U'). In other words, an RUR decomposition expresses a desired matrix as a linear combination of r 2 outer products of rows of the input matrix ( ∑︁ i,j [f ¯(CC † )](i, j)R(i, •) † R(j, •), for example). 6 We want our output in the form of an RUR decomposition, since we can describe such a decomposition implicitly just as a list of row indices and some additional coefficients, which avoids picking up a dependence on m or n in our runtimes. Further, having SQ(A) gives us SQ ϕ (R † U R) via closure properties, enabling efficient access to matrix-vector expressions like R † U Rb. More general results follow as corollaries of our main result on even SVT (Theorems 7.1 and 7.2). However, using only our main theorem about even SVT, we can directly recover most existing quantum-inspired machine learning algorithms without using these general results, yielding faster dequantization for QML algorithms. We now outline our results recovering such applications.
We use the results above to recover existing quantum-inspired algorithms for recommendation systems [Tan19], principal component analysis [Tan21], supervised clustering [Tan21], support vector machines [DBH21], low-rank matrix inversion [CGL + 20], and semidefinite program solving [CLLW20]. We also propose new quantum-inspired algorithms for low-rank Hamiltonian simulation and discriminant analysis (dequantizing the quantum algorithm of Cong & Duan [CD16]).
Fig. 1 has a summary of our results, along with a comparison of runtimes to the corresponding quantum algorithms and prior quantum-inspired work, where it exists. All our results match or improve on prior dequantized algorithms apart from that for matrix inversion, where prior work gives an incomparable runtime that only holds for strictly low-rank matrices of rank k. Our results for matrix inversion and semidefinite program solving solve the problem in greater generality than prior work, without the restriction that the input matrices are strictly rank-k. 7We do not claim any meaningful breakthroughs for these problems in the classical literature: the problems that these QML algorithms solve differ substantially from their usual classical counterparts. For example, the quantum recommendation systems algorithm of Kerenidis and Prakash [KP17] performs sampling from a low-rank approximation of the input instead of low-rank matrix completion, which is the typical formalization of the recommendation systems problem [Tan19]. Evaluating these quantum algorithms' justifications for their versions of problems is outside the scope of this work: instead, we argue that these algorithms would likely not give exponential speedups when implemented, regardless of whether such implementations would be useful. The goal of our framework is to demonstrate what can be done classically and establish a classical frontier for quantum algorithms to push past.
The proofs for these dequantization results follow the same general structure: consider the quantum algorithm and formulate the problem that this algorithm solves, and in particular, the linear algebra expression that the quantum algorithm computes. From there, repeatedly use the SVT result and key lemma to approximate this expression by something like an RUR decomposition. Finally, use closure properties to gain oversampling and query access to that output decomposition. This procedure is relatively straightfoward and flexible. Also, unlike previous work [CGL + 20, CLLW20], our results need not assume that the input is strictly low-rank. Instead, following [Tan19, GSLW19], our algorithms work on close-to-low-rank matrices by doing SVTs that smoothly threshold to effectively only operate on large-enough singular values.
Quantum-inspired algorithms. Our approach and analysis is much simpler than that of Frieze, Kannan, and Vempala [FKV04], while it also gives improved results in our applications, and has several other advantages. For example, the reduction to [FKV04] first given by Tang to get an quantum algorithm prior work this work simple QSVT [GSLW19], §6.1
discriminant analysis [CD16], §6.9
Figure 1: The time complexity for our algorithms, the quantum algorithms they are based on, and prior quantum-inspired algorithms (where they exist). We assume our sampling and query accesses to the input takes O(1) time. There are data structures that can support such queries (Remark 3.11), and if the input is in QRAM, the runtime only increases by at most a factor of log of input size.
We list the runtime of the algorithm, not including the time it takes to access the output (denoted with ˜︂ sq). The runtimes as listed ignore polylog terms, particularly those in error parameters (ε and δ) and dimension parameters (m and n). The matrices and vectors referenced in these runtimes are always the input, σ refers to a singular value threshold of the input matrices, λ refers to an eigenvalue threshold (which can be thought of here as σ 2 ), and η > ε is a (dimensionless) gap parameter.
(♣) indicates that the error analyses of the corresponding results are incomplete; we list the runtime they achieve for completeness.
(♢) indicates that the corresponding results only hold in the restricted setting where the input matrices are strictly rank k. For the quantum algorithms with this tag, they allow for general matrices, but only have an informal error analysis arguing that singular values outside the range considered don't affect the final result.
SVT-based low-rank approximation bound from the standard notion of low-rank approximation [Tan19, Theorem 4.7] induces a quadratic loss in precision, which appears to be only an artifact of the analysis. Also, [FKV04] gives Frobenius norm error bounds, though for applications we often only need spectral norm bounds; our main theorem can get improved runtimes by taking advantage of the weaker spectral norm bounds. Finally, we take a reduced number of rows compared to columns, whereas [FKV04] approximates the input by taking the same number of rows and columns.
Randomized numerical linear algebra. All of the results presented here are more or less randomized linear algebra algorithms [Mah11,Woo14]. The kind of sampling we get from sampling and query access is called importance sampling or length-square sampling in that body of work: see the survey by Kannan and Vempala [KV17] for more on importance sampling. Importance sampling, and specifically, its approximate matrix product property, is the core primitive of this work. In addition to the low-rank approximation algorithms [FKV04] used in the quantum-inspired literature, others have used importance sampling for, e.g., orthogonal tensor decomposition [DM07, MMD08, SWZ16] (generalizing low-rank approximation [FKV04]) and support vector machines [HKS11].
The fundamental difference between quantum-inspired algorithms and traditional sketching algorithms is that we assume "we can perform quantum measurements" of states corresponding to input in time independent of input dimension (that is, we have efficient sampling and query access to input), and in exchange want algorithms that run in time independent of dimension and provide only (over)sampling and query access to the output. This quantum-inspired model is weaker than the standard sketching algorithm model (Remark 3.11): an algorithm taking T time in the quantum-inspired model for an input matrix A can be converted to a standard algorithm that runs in time O(nnz(A) + T ), where nnz(A) is the number of nonzero entries of A. So, we can also think about an O(T )-time quantum-inspired algorithm as an O(nnz(A) + T )-time sketching algorithm, where the nnz(A) portion of the runtime can only be used to facilitate importance sampling. 8 This restriction makes for algorithms that may perform worse in generic sketching settings, but work in more settings, and so demonstrate lack of exponential quantum speedup for a wider range of problems.
A natural question is whether more modern sketching techniques can be used in our model. After all, importance sampling is only one of many sketching techniques studied in the large literature on sketching algorithms. Notably, though, other types of sketches seem to fail in the input regimes where quantum machine learning succeeds: assuming sampling and query access to input, importance sampling takes time independent of dimension, whereas other randomized linear algebra methods such as Count-Sketch and Johnson-Lindenstrauss still take time linear in input-sparsity.
Subsequent work by Chepurko, Clarkson, Horesh, Lin, and Woodruff [CCH + 22] notes that importance sampling oversamples leverage score sampling, so usual analyses for leverage score sampling also hold for importance sampling, up to some small overhead. It is reasonable to suspect that this connection could lead to significant improvements over the results presented here. However, exploiting this connection for improved runtimes seems nontrivial, since most approaches using leverage score sampling requires performing O(nnz(A))-time pre-processing operations, even if one has SQ(A). As a simple example, for low-rank approximation of an input matrix A, if we wish to adapt the algorithm of Clarkson and Woodruff [CW17,Woo14], importance sampling of A can replace CountSketch for one of the sketches [Woo14, Lemma 4.2], but importance sampling of SA cannot replace the other [Woo14,Theorem 4.3]. Versions of importance sampling may work for the second sketch (for example, sampling from a low-rank approximation of SA), but we are aware of none that can be obtained easily from SQ(A). This may be a manifestation of the difficulty of achieving relative-error estimates to quantities like leverage scores and residuals in this model.
An alternative approach is to use a projection-cost preserving sketch (PCP) like ridge leverage score sampling to sketch A on both sides [CMM17]. The importance sampling from SQ(A) λ/2oversamples ridge leverage score sampling, where λ := ∥A∥ 2 F /∥A -A k ∥ 2 F ≥ ∥A∥ 2 F /σ 2 k+1 , so using importance sampling in place of ridge leverage score sampling can give algorithms. This is how [CCH + 22] gets their algorithms for low-rank sampling (recommendation systems) and quantuminspired linear regression. This does appear to be a promising approach, with possibility to extend to dequantizing all of QSVT, but to the authors' knowledge, it is still an open question how to improve the algorithms presented in this work, with the exception of [CCH + 22] improving over our recommendation systems algorithm. Their use of PCPs significantly improves the runtime for this low-rank sampling task down below ∥A∥ 6 F σ 6 ε 6 , but getting a similarly good runtime for all functions seems nontrivial. For example, their linear regression algorithm requires that the input matrix is strictly rank-k (or is regularized).
The quantum-like closure properties of importance sampling shown here may be useful in the context of classical sketching algorithms. This insight unlocks surprising power in importance sampling. For example, it reveals that Frieze, Kannan, and Vempala's low-rank approximation algorithm [FKV04], which, as stated, requires O(kmn) time to output the desired matrix, actually can produce useful results (samples and entries) in time independent of input dimension. To use the language in [FKV04], if Assumptions 1 and 2 hold for the input matrix, they also hold for the output matrix! Classical algorithms for quantum problems. We are aware of two important prior results from before Tang's first paper [Tan19] that connect quantum algorithms to randomized numerical linear algebra. The first is Van den Nest's work on using probabilistic methods for quantum simulation [VdN11], which defines a notion of "computationally tractable" (CT) state equivalent to our notion of sampling and query access and then uses it to simulate restricted classes of quantum circuits. We share some essential ideas with this work, such as the simple sampling lemmas Lemmas 3.6 and 4.12, but also differ greatly since we focus on low-rank matrices relevant for QML, whereas [VdN11] focuses on simulating potentially large quantum circuits that correspond to high-rank matrices. The second is a paper by Rudi, Wossnig, Ciliberto, Rocchetto, Pontil, and Severini [RWC + 20] that uses the Nyström method to simulate a sparse Hamiltonian H on a sparse input state in time poly-logarithmic in dimension and polynomial in ∥H∥ F , assuming sampling and query access to H. Our Hamiltonian simulation results do not require a sparsity assumption and still achieve a dimension-independent runtime, but get slightly larger exponents in exchange.
Practical implementation. A work by Arrazola, Delgado, Bardhan, and Lloyd [ADBL20] implements and benchmarks quantum-inspired algorithms for regression and recommendation systems. The aforementioned paper of Chepurko, Clarkson, Horesh, Lin, and Woodruff [CCH + 22] does the same for the quantum-inspired algorithms they introduce. The former work makes various conclusions, including that the ε 2 scaling in the number of rows/columns taken in our recommendation systems algorithm is inherent and that the quantum-inspired algorithms performed slower and worse than direct computation for practical datasets. The latter work finds that their algorithms perform faster than direct algorithms, with an accompanying increase in error comparable to that of other sketching algorithms [DKW18]. This improvement appears to come from both a better-performing implementation as well as an algorithm with better asymptotic runtime. Nevertheless, it is difficult to draw definitive conclusions about the practicality of quantum-inspired algorithms as a whole from these experimental results. Since quantum-inspired algorithms are a restricted, weaker form of computation than classical randomized numerical linear algebra algorithms (see the comparison made above), it seems possible that they perform worse than standard sketching algorithms, despite seemingly having exponentially improved runtime in theory.
Modern sketching algorithms use similar techniques to quantum-inspired algorithms, but are more natural to run on a classical computer and are likely to be faster. For example, Dahiya, Konomis, and Woodruff [DKW18] conducted an empirical study of sketching algorithms for lowrank approximation on both synthetic datasets and the movielens dataset, reporting that their implementation "finds a solution with cost at most 10 times the optimal one . . . but does so 10 times faster." Sketching algorithms like those in [DKW18] may become a relevant point of reference for benchmarking quantum linear algebra, when the implementation of these quantum algorithms on actual quantum hardware becomes possible. In a sense, our work shows using asymptotic runtime bounds that in many scenarios sketching and sampling techniques give similar computational power to quantum linear algebra, which is a counterintuitive point since the former typically leads to linear runtimes and the latter leads to poly-logarithmic ones.
Quantum machine learning. Our work has major implications for the landscape of quantum machine learning. Since we have presented many dequantized versions of QML algorithms, the question remains of what QML algorithms don't have such versions. In other words, what algorithms still have the potential to give exponential speedups?
There are two general paradigms for employing quantum linear algebra techniques in quantum machine learning: the low-rank approach and the high-rank approach. In both, we need to turn classical input data vectors to quantum states (via what's called an amplitude encoding), perform linear algebra operations on those vectors, and extract information about the output via sampling [Aar15]. Since preparing generic quantum states require a number of quantum gates that is proportional to the dimension of the vectors, we need some state preparation assumptions (like having QRAM with an appropriate data structure, etc., cf. Remark 3.11) in order to achieve sublinear runtimes. The main difference between the two paradigms is how the input matrices are given: in the low-rank approach, they are also given in QRAM, and so must be rescaled to have Frobenius norm one. QSVT with block-encodings coming from the QRAM data structurefoot_5 or density operators is an example framework in this vein. The restrictions of this block-encoding method means that the rank needs to be small, but assuming the hardware needed for QRAM can be realized, this is still a flexible setting, one that QML researchers find interesting. In the high-rank approach, which is used in the HHL algorithm [HHL09] and its derivatives, the matrix needs to be represented by a concise quantum circuit and have a small (poly-logarithmic in input dimension) condition number in order to gain an exponential speedup over classical algorithms. This doesn't happen in typical datasets. The collection of these demanding requirements hamstrings most attempts to find applications of HHL [HHL09] with the potential for practical super-polynomial speedups.
Our results give evidence for the lack of exponential speedup for the former, low-rank approach. It is important to note however, that our results do not rule out the possibility of large polynomial quantum speedups. In order to assess the potential usefulness of QML algorithms in this regime it is important to improve on classical upper and lower bounds for these problems, which we leave as an open question. On the other hand, high-rank block-encodings, such as those coming from sparsity assumptions in the original HHL algorithm [HHL09], remain impervious to our techniques. This suggests that the most promising way to get exponential quantum speedups for QML algorithms is by assuming sparse matrices as input, or utilizing other efficiently implementable high-rank quantum operation such as the Quantum Fourier Transform.
Works from Zhao, Fitzsimons, and Fitzsimons on Gaussian process regression [ZFF19]; from Lloyd, Garnerone, and Zanardi on topological data analysis [LGZ16,GCD20]; and from Yamasaki, Subramanian, Sonoda, and Koashi [YSSK20] on learning random features attempt to address these issues to get a super-polynomial quantum speedup. Though these works avoid the dequantization barrier to large quantum speedups, it remains to be seen how broad will be their impact on QML and whether these speedups manifest for data seen in practice.
Related independent work. Independently from our work, Jethwani, Le Gall, and Singh simultaneously derived similar results [JLGS20]. They implicitly derive a version of our even SVT result, and use it to achieve generic SVT (approximate SQ(b † f (SV) (A)) for a vector b) by writing
x and then using sampling subroutines to get the solution from the resulting expression b † AR † U R. It is difficult to directly compare the main SVT results, because the parameters that appear in their runtime bounds are somewhat non-standard, but one can see that for typical choices of f , their results require a strictly low-rank A. In comparison our results apply to general A, and we also demonstrate how to apply them to (re)derive dequantized algorithms.
Our framework recovers recent dequantization results, and we hope that it will be used for dequantizing more quantum algorithms. In the meantime, our work leaves several natural open questions:
(a) Is there an approach to QML that does not go through HHL (whose demanding assumptions make exponential speedups difficult to demonstrate even in theory) or a low-rank assumption (which, as we demonstrate, makes the tasks "easy" for classical computers) and yields a provable superpolynomial speedup for a practically relevant ML problem? (b) Our algorithms still have significant slowdown as compared to their quantum counterparts.
Can we shave condition number factors to get runtimes of the form ˜︁ O (︂
(for the recommendation systems application, for instance), without introducing additional assumptions? Can we get even better runtimes by somehow avoiding SVD computation? (c) Do the matrix arithmetic closure properties we showed for ℓ 2 -norm importance sampling hold for other kinds of sampling and sketching distributions, like leverage score or ℓ p -norm sampling?
(d) In the quantum setting, linear algebra algorithms [GSLW19] can achieve logarithmic dependence on the precision ε. Can classical algorithms also achieve such exponentially improved dependence, when the goal is restricted to sampling from the output (i.e., without the requirement to query elements of the output)? If not, is there a mildly stronger classical model that can achieve this? Can one prove that this exponential advantage for sampling problems cannot be conferred to estimation/decision problems?
The paper proceeds as follows. Section 3 introduces the notion of (over)sampling and query access and some of its closure properties. Section 4 gives the fundamental idea of using sampling and query access to sketch matrices used for the approximation results in Section 4.1 and singular value transformation results in Section 5. These results form the framework that is used to dequantize QSVT in Section 6.1 and recover all the quantum-inspired results in Section 6. These applications of our framework contain various tricks and patterns that we consider to be "best practice" for coercing problems into our framework, since they have given us the best complexities and generality. More general results of SVT are shown in Section 7.
To begin with, we define notation to be used throughout this paper. For n ∈ N, [n] := {1, . . . , n}.
For z ∈ C, its absolute value is |z| = √ z * z, where z * is the complex conjugate of z. f ≲ g denotes the ordering f = O(g) (and respectively for ≳ and ≂). ˜︁ O(g) is shorthand for O(g poly(log g)). log refers to the natural logarithm. Finally, we assume that arithmetic operations (e.g., addition and multiplication of real numbers) and function evaluation oracles (computing f (x) from x) take unit time, and that queries to oracles (like the queries to input discussed in Section 3) are at least unit time cost.
In this paper, we consider complex matrices A ∈ C m×n for m, n ∈ N. For i ∈ [m], j ∈ [n], we let A(i, •) denote the i-th row of A, A(•, j) denote the j-th column of A, and A(i, j) denote the (i, j)-th element of A. (A | B) denotes the concatenation of matrices A and B and vec(A) ∈ C mn denotes the vector formed by concatenating the rows of A. For vectors v ∈ C n , ∥v∥ denotes standard Euclidean norm (so ∥v∥ := (
and the spectral norm of A is ∥A∥ := ∥A∥ Op := sup x∈C n ,∥x∥=1 ∥Ax∥. We say that U is an isometry if ∥U x∥ = ∥x∥ for all x, or equivalently, if U is a subset of columns of a unitary.
A singular value decomposition (SVD) of A is a representation A = U DV † , where for N := min(m, n), U ∈ C m×N and V ∈ C n×N are isometries and D ∈ R N ×N is diagonal with σ i := D(i, i) and σ 1 ≥ σ 2 ≥ • • • ≥ σ N ≥ 0. We can also write this decomposition as A = ∑︁ N i=1 σ i u i v † i , where u i := U (•, i) and v i := V (•, i). For Hermitian A, an (unitary) eigendecomposition of A is a singular value decomposition where U = V , except the entries of D are allowed to be negative.
Using SVD, we can define the rank-k approximation of A to be
We now formally define singular value transformation:
Definition 2.1. For a function f : [0, ∞) → C such that f (0) = 0 and a matrix A ∈ C m×n , we define the singular value transform of A via a singular value decomposition A = ∑︁ min(m,n)
The requirement that f (0) = 0 ensures that the definition is independent of the (not necessarily unique) choice of SVD.
Definition 2.2. For a function f : R → C and a Hermitian matrix A ∈ C n×n , we define the eigenvalue transform of A via a unitary eigendecomposition
Since we only consider eigenvalue transformations of Hermitian matrices, where singular vectors/values and eigenvectors/values (roughly) coincide, the key difference between singular value transformation and eigenvalue transformation is that the latter can distinguish eigenvalue sign. As eigenvalue transformation is the standard notion of a matrix function, we will usually drop the superscript in notation: f (A) := f (EV) (A).
We will use the following standard definition of a Lipschitz function.
Definition 2.3. We say f :
We define approximate isometry as follows:
If V is an α-approximate isometry, among other things, it implies that |∥V ∥ 2 -1| ≤ α and that there exists an isometry U ∈ C m×n with im(U ) = im(V ) such that ∥U -V ∥ ≤ α. We show this and other basic facts in the following lemma, whose proof is deferred to Appendix B.
Lemma 2.5. If X ˆ∈ C m×n is an α-approximate isometry, then there is an exact isometry X ∈ C m×n with the same columnspace as X ˆsuch that ∥X ˆ-X∥ ≤ α. Furthermore, for any matrix
Since we want our algorithms to run in time sublinear in input size, we must carefully define our access model. The sampling and query oracle we present below is unconventional, being designed as a reasonable classical analogue for the input model of some quantum algorithms. It will also be used heavily to move between intermediate steps of these quantum-inspired algorithms. First, as a warmup, we define a simple query oracle: Definition 3.1 (Query access). For a vector v ∈ C n , we have Q(v), query access to v, if for all i ∈ [n], we can query for v(i). Likewise, for a matrix A ∈ C m×n , we have Q(A) if for all (i, j) ∈ [m] × [n], we can query for A(i, j). Let q(v) (respectively q(A)) denote the (time) cost of such a query. For example, in the typical RAM access model, we are given our input v ∈ C n as Q(v) with q(v) = 1. For brevity, we will sometimes abuse this notation (and other access notations) and, for example, abbreviate "Q(A) for A ∈ C m×n " as "Q(A) ∈ C m×n ". We will also sometimes abuse complexity notation like q to refer to known bounds on the complexity, instead of the complexity itself.
Definition 3.2 (Sampling and query access to a vector). For a vector v ∈ C n , we have SQ(v), sampling and query access to v, if we can:
1. query for entries of v as in Q(v); 2. obtain independent samples i ∈ [n] following the distribution D v ∈ R n , where D v (i) := |v(i)| 2 /∥v∥ 2 ; 3. query for ∥v∥.
Let q(v), s(v), and n(v) denote the cost of querying entries, sampling indices, and querying the norm respectively. Further define sq(v) := max(q(v), s(v), n(v)).
We will refer to these samples as importance samples from v, though one can view them as measurements of the quantum state |v⟩ := 1 ∥v∥ ∑︁ v i |i⟩ in the computational basis. Quantum-inspired algorithms typically don't give exact sampling and query access to the output vector. Instead, we get a more general version of sampling and query access, which assumes we can only access a sampling distribution that oversamples the correct distribution.foot_7 Definition 3.3. For p, q ∈ R n ≥0 that are distributions, meaning
The motivation for this definition is the following: if p ϕ-oversamples q, then we can convert a sample from p to a sample from q with probability 1/ϕ using rejection sampling: sample an i distributed as p, then accept the sample with probability q(i)/(ϕp(i)) (which is ≤ 1 by definition). Definition 3.4 (Oversampling and query access). For v ∈ C n and ϕ ≥ 1, we have SQ ϕ (v), ϕoversampling and query access to v, if we have Q(v) and SQ(v ˜) for v ˜∈ C n a vector satisfying ∥v ˜∥2 = ϕ∥v∥ 2 and |v ˜(i)| 2 ≥ |v(i)| 2 for all i ∈ [n]. Denote s ϕ (v) := s(v ˜), q ϕ (v) := q(v ˜), n ϕ (v) := n(v ˜), and sq ϕ (v) := max(s ϕ (v), q ϕ (v), q(v), n ϕ (v)).
The distribution D v ˜ϕ-oversamples D v , since for all i ∈ [n],
For this reason, we call D v ˜a ϕ-oversampled importance sampling distribution of v. SQ(v) is the same as SQ 1 (v), by taking v ˜= v. Note that we do not assume knowledge of ϕ (though it can be estimated, (though it can be estimated as shown in Lemma 3.5). However, we do need to know ∥v ˜∥ (even if ∥v∥ is known), as it cannot be deduced from a small number of queries, samples, or probability computations. So, we will be choosing v ˜(and, correspondingly, ϕ) such that ∥v ˜∥2 remains computable, even if potentially some cv ˜satisfies all our other requirements for some c < 1 (giving a smaller value of ϕ).
Intuitively speaking, estimators that use D v can also use D v ˜via rejection sampling at the expense of a factor ϕ 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 3.5. Suppose we are given SQ ϕ (v) and some δ ∈ (0, 1]. Denote ˜︂ sq(v) := ϕ sq ϕ (v) log 1 δ . We can sample from D v with probability ≥ 1 -δ in O(˜︂ sq(v)) time. We can also estimate ∥v∥ to ν multiplicative error for ν ∈ (0, 1] with probability
Proof. Consider the following rejection sampling algorithm to generate samples: sample an index i from v ˜, and output it as the desired sample with probability r(i)
)︁ time and r(i) ≤ 1 since v ˜bounds v. The probability of accepting a sample in a round is
and, conditioned on a sample being accepted, the probability of it being i is |v(i)| 2 /∥v∥ 2 , so the output distribution is D v as desired. So, to get a sample with ≥ 1 -δ probability, run rejection sampling for at most 2ϕ log 1 δ rounds. To estimate ∥v∥ 2 , notice that we know ∥v ˜∥2 , so it suffices to estimate ∥v∥ 2 /∥v ˜∥2 which is ϕ -1 . The probability of accepting the rejection sampling routine is ϕ -1 , so we run 3ν -2 ϕ log 2 δ rounds of it for estimating ϕ -1 . Let Z denote the fraction of them which end in acceptance. Then, by a Chernoff bound we have
so Z∥v ˜∥2 is a good multiplicative approximation to ∥v∥ 2 with probability ≥ 1 -δ.
Generally, compared to a quantum algorithm that can output (and measure) a desired vector |v⟩, our algorithms will output SQ ϕ (u) such that ∥u -v∥ is small. So, ˜︂ sq(u) is the relevant complexity measure that we will analyze and bound: if we wish to mimic samples from the output of the quantum algorithm we dequantize, we will pay a one-time cost to run our quantum-inspired algorithm for "obtaining" SQ ϕ (u), and then pay ˜︂ sq(u) cost per additional measurement. 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.
) one-time pre-processing cost to query for norms).
Proof. Denote u := ∑︁ λ t v t . To compute u(s) for some s ∈ [n], we just need to query v t (s) for all t ∈ [τ ], paying O( ∑︁ q(v t )) cost. So, it suffices to get SQ(u ˜) for an appropriate bound u ˜. We choose
We have SQ(u ˜): we can compute ∥u ˜∥2 by querying for all norms ∥v ˜t∥, compute u ˜(s) by querying v ˜t(s) for all t ∈ [τ ]. We can sample from u ˜by first sampling t ∈ [τ ] with probability ∥λtv ˜t∥ 2 ∑︁ ℓ ∥λ ℓ v ˜ℓ∥ 2 , and then taking our sample to be j ∈ [n] from v ˜t. The probability of sampling j ∈ [n] is correct:
If we pre-process by querying all the norms ∥v ˜ℓ∥ in advance, we can sample from the distribution over i's in O(1) time, using an alias sampling data structure for the distribution (Remark 3.11), and we can sample from v ˜t using our assumed access to it, SQ φt (v t ).
So, our general goal will be to express our output vector as a linear combination of a small number of input vectors that we have sampling and query access to. Then, we can get an approximate SQ access to our output using Lemma 3.5, where we pay an additional "cancellation constant" factor of ϕ = τ ∑︁ φt∥λtvt∥ 2 ∥ ∑︁ λtvt∥ 2 . 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 α 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.
We also define oversampling and query access for a matrix. The same model (under an alternative definition) is also discussed in prior work [FKV04,DKR02] and is the right notion for the sampling procedures we will use. Definition 3.7 (Oversampling and query access to a matrix). For a matrix A ∈ C m×n , we have SQ(A) if we have SQ(A(i, •)) for all i ∈ [m] and SQ(a) for a ∈ R m the vector of row norms (a(i) := ∥A(i, •)∥).
We have SQ ϕ (A) if we have Q(A) and SQ(A ˜) for
The complexity of (over)sampling and querying from the matrix A is denoted by
, and n ϕ (A) := n(a ˜) respectively. We also denote sq ϕ (A) := max(s ϕ (A), q ϕ (A), q(A), n ϕ (A)). We omit subscripts if ϕ = 1.
Observe that access to a matrix, SQ ϕ (A), implies access to its vectorized version, SQ ϕ (vec(A)):
we can take ˜︂ vec(A) = vec(A ˜), and the distribution for vec(A ˜) is sampled by sampling i from D a ˜, and then sampling j from D A ˜(i,•) . This gives the output (i, j) with probability |A ˜(i, j)| 2 /∥A ˜∥2 F . Therefore, one can think of SQ ϕ (A) as SQ ϕ (vec(A)), with the addition of having access to samples (i, j) from vec(A), conditioned on fixing a particular row i and also knowing the probabilities of these conditional samples. Now we prove that oversampling and query access is closed under taking outer products. The same idea also extends to taking Kronecker products of matrices.
Proof. We can query an entry A(i, j) = u(i)v(j) † by querying once from u and v. Our choice of upper bound is A ˜= u ˜v ˜ †. Clearly, this is an upper bound on uv † and ∥A ˜∥2
after querying for u ˜(i), and a ˜= ∥v ˜∥2 u ˜, so we have SQ(a ˜) from SQ(u ˜) after querying for ∥v ˜∥.
Using the same ideas as in Lemma 3.6, we can extend sampling and query access of input matrices to linear combinations of those matrices.
)︁ cost. So, it suffices to get SQ(A ˜) for an appropriate bound A ˜. We choose
That |A ˜(i, j)| ≥ |A(i, j)| follows from Cauchy-Schwarz, and we get the desired value of ϕ:
We have SQ(A ˜): we can compute ∥A ˜∥F by querying for all norms ∥A
)︂ time.
Remark 3.10. With the lemmas we've introduced, we can already get oversampling and query access to some modest expressions. For example, consider RUR decompositions, which show up frequently in our results: suppose we have SQ(A) for A ∈ C m×n , R ∈ C r×n a (possibly normalized) subset of rows of A, and a matrix U ∈ C r×r . Then
which is a linear combination of r 2 outer products involving rows of A. So, by Lemma 3.8 and Lemma 3.9, we have SQ ϕ (R † U R).
For us, the most interesting scenario is when our sampling and query oracles take poly-logarithmic time, since this corresponds to the scenarios where quantum state preparation procedures can run in time polylog(n). In these scenarios, quantum machine learning have the potential to achieve exponential speedups. We can provide such classical access in various ways.
Remark 3.11. Below, we list settings where we have sampling and query access to input matrices and vectors, and whenever relevant, we compare the resulting runtimes to the time to prepare analogous quantum states. Note that because we do not analyze classical algorithms in the bit model, i.e., we do not count each operation bitwise, their runtimes may be missing log factors that should be counted for a fair comparison between classical and quantum. If the input vector (resp. matrix) is given as a list of nnz(v) (resp. nnz(A)) of its non-zero entries, then the pre-processing time is linear in that number of entries. 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 O(T sq(A)) for T independent of input size) implies an algorithm in the standard computational model (running in O(nnz(A) + T ) time).
(b) (Dynamic data structure) QML algorithms often assume that their input is in a data structure with a certain kind of quantum access [Pra14, KP20, GLM08, WZP18, RSW + 19, CGJ19]. They argue that, since this data structure allows for circuits preparing input states with linear gate count but polylog depth, hardware called QRAM might be able to parallelize these circuits enough so that they run in effectively polylog time. In the interest of considering the best of all possible worlds for QML, we will treat circuit depth as runtime for QRAM and ignore technicalities.
This data structure (see Fig. 2) admits sampling and query access to the data it stores with just-as-good runtimes: specifically, for a matrix A ∈ C m×n , we get SQ(A) with q(A) = O(1), s(A) = O(log mn), and n(A) = O(1). So, quantum-inspired algorithms can be used whenever QML algorithms assume this form of input.
Further, unlike the alias method stated above, this data structure supports updating entries in O(log mn) time, which is used in applications of QML where data accumulates over time [KP17]. , simply by using the uniform distribution over non-zero entries for the oversampling and query oracles. For example, for SQ(a ˜) we can set a ˜(i) := c √ s, and for A ˜(i, •) we use the vector with entries c at the non-zeros of A(i, •) (potentially adding some "dummy" zero locations to have exactly s non-zeroes).
Note that similar sparse-access assumptions are often seen in the QML and Hamiltonian simulation literature [HHL09]. Also, if A is not much smaller than we expect, then ϕ can be independent of dimension. For example, if A has exactly s non-zero entries per row and
(f) (CT states) In 2009, Van den Nest defined the notion of a "computationally tractable" (CT) state [VdN11]. Using our notation, |ψ⟩ ∈ C n is a CT state if we have SQ(ψ) with sq(ψ) = polylog(n). Van den Nest's paper identifies several classes of CT states, including product states, quantum Fourier transforms of product states, matrix product states of polynomial bond dimension, stabilizer states, and states from matchgate circuits. For more details on how can one get efficient sampling and query access to such vectors we direct the reader to [VdN11].
We now introduce the workhorse of our algorithms: the matrix sketch. Using sampling and query access, we can generate these sketches efficiently, and these allow one to reduce the dimensionality of a problem, up to some approximation. Most of the results presented in this section are known in the classical sketching literature: we present them here for completeness, and to restate them in the context of sampling and query access.
Definition 4.1. For a distribution p ∈ R m , we say that a matrix S ∈ R s×m is sampled according to p if each row of S is independently chosen to be e i / √︁ s • p(i) with probability p(i), where e i is the vector that is one in the ith position and zero elsewhere. If p is an ℓ 2 -norm sampling distribution D v as defined in Definition 3.2, then we also say S is sampled according to v.
We call S an importance sampling sketch for A ∈ C m×n if it is sampled according to A's row norms a, and we call S a ϕ-oversampled importance sampling sketch if it is sampled according to the bounding row norms from SQ ϕ (A), a ˜(or, more generally, from a ϕ-oversampled importance sampling distribution of a).
One should think of S as a description of how to sketch A down to SA. The following lemma shows that ∥SA∥ F approximates ∥A∥ F , giving a simple example of the phenomenon that SA approximates A in certain senses: it shows that ∥SA∥ F = Θ(∥A∥ F ) with probability ≥ 0.9 when S has Ω( 1 ϕ 2 ) rows. We show later (Lemma 4.8) that a similar statement holds for spectral norm: ∥SA∥ = Θ(∥A∥) with probability ≥ 0.9 when S has Ω ˜(ϕ 2 ∥A∥ 2 F /∥A∥ 2 ) rows.
Lemma 4.2 (Frobenius norm bounds for matrix sketches). Let S ∈ C r×m be a ϕ-oversampled importance sampling sketch of A ∈ C m×n . Then ∥[SA](i, •)∥ ≤ √︁ ϕ/r∥A∥ F for all i ∈ [r], so ∥SA∥ 2 F ≤ ϕ∥A∥ 2 F (unconditionally). Equality holds when ϕ = 1. Further,
Proof. Let p be the distribution used to create S, and let s i be the sample from p used for row i of S. Then ∥SA∥ 2 F is the sum of the row norms ∥[SA](i, •)∥ 2 over all i ∈ [r], and
The first equation shows the unconditional bounds on ∥SA∥ F . When ϕ = 1, p(i) = ∥A(i, •)∥ 2 /∥A∥ 2 F so the inequality becomes an equality. By the second equation, ∥SA∥ 2 F -∥A∥ 2 F has expected value zero and is the sum of independent random variables bounded in [-∥A∥ 2 F , (ϕ -1)∥A∥ 2 F ], so the probabilistic bound follows immediately from Hoeffding's inequality.
In the standard algorithm setting, computing an importance sampling sketch requires reading all of A, since we need to sample from D a . If we have SQ ϕ (A), though, we can efficiently create a
)︁ time: for each row of S, we pull a sample from p, and then compute √︁ p(i). After finding this sketch S, we have an implicit description of SA: it is a normalized multiset of rows of A, so we can describe it with the row indices and corresponding normalization, (i 1 , c 1 ), . . . , (i s , c s ).
SA can be used to approximate matrix expressions involving A. Further, we can chain sketches using the lemma below, which shows that from SQ ϕ (A), we have SQ ≤2ϕ ((SA) † ), under a mild assumption on the size of the sketch S. This can be used to find a sketch T † of (SA) † . The resulting expression SAT is small enough that we can compute functions of it in time independent of dimension, and so will be used extensively. When we discuss sketching A down to SAT , we are referring to the below lemma for the method of sampling T . Lemma 4.3. Consider SQ φ (A) ∈ C m×n and S ∈ R r×m sampled according to a ˜, described as pairs (i 1 , c 1 ), . . . , (i r , c r ). If r ≥ 2φ 2 log 2 δ , then with probability ≥ 1-δ, we have SQ ϕ (SA) and SQ ϕ ((SA) † ) for some ϕ satisfying ϕ ≤ 2φ. If φ = 1, then for all r, we have SQ(SA) and SQ((SA) † ).
The runtimes for SQ ϕ (SA) are q(SA) = q(A), s ϕ (SA) = s φ (A), q ϕ (SA) = q φ (A), and
Proof. By Lemma 4.2, ∥SA∥ 2 F ≥ ∥A∥ 2 F /2 with probability ≥ 1 -δ. Suppose this bound holds. To get SQ ϕ (SA), we take ˜︂ SA = SA ˜, which bounds SA by inspection. Further, ∥SA ˜∥2
Analogously, (SA ˜) † works as a bound for SQ ϕ ((SA) † ). We can query an entry of SA by querying the corresponding entry of A, so all that suffices is to show that we have SQ(SA ˜) and SQ((SA ˜) † ) from SQ(A ˜). (When φ = 1, we can ignore the above argument: the rest of the proof will show that we have SQ(SA) and SQ((SA) † ) from SQ(A).)
We have SQ(SA ˜). Because the rows of SA ˜are rescaled rows of A ˜, we have SQ access to them from SQ access to A ˜. Because ∥SA ˜∥2
∥A ˜∥2 F , we have SQ access to the vector of row norms of SA ˜(pulling samples simply by pulling samples from the uniform distribution).
We have SQ((SA ˜) † ). (This proof is similar to one from [FKV04].) Since the rows of (SA ˜) † are length r, we can respond to SQ queries to them by reading all entries of the row and performing some linear-time computation. ∥(SA ˜) † ∥ 2 F = ∥A ˜∥2 F , so we can respond to a norm query by querying the norm of A ˜. Finally, we can sample according to the row norms of (SA ˜) † by first querying an index i ∈ [r] uniformly at random, then outputting the index j ∈ [n] sampled from [SA ˜](i, •) (which we can sample from because it is a row of A ˜). The distribution of the samples output by this procedure is correct: the probability of outputting j is
Here, we present approximation results on sketched matrices that we will use heavily throughout our results. We begin with a fundamental observation: given sampling and query access to a matrix A, we can approximate the matrix product A † B by a sum of rank-one outer products. We formalize this with two variance bounds, which we can use together with Chebyshev's inequality. Consider X ∈ C m×n , Y ∈ C m×p , and take S ∈ R r×m to be sampled according to p ∈ R m a ϕ-oversampled importance sampling distribution from X or Y . Then,
Proof. To show the first equation, we use that
] is a sum of variances, one for each entry (i, j), since E[X † S † SY -XY ] is zero in every entry. Furthermore, for every entry (i, j), the matrix expression is the sum of r independent, mean-zero terms, one for each row of S:
So, we can use standard properties of variancesfoot_8 to conclude that
The second other inequality follows by the same computation:
The above result shows that, given SQ(X), X † Y can be approximated by a sketch with constant failure probability. If we have SQ(X) and SQ(Y ), we can make the failure probability exponential small. To show this tighter error bound, we use an argument of Drineas, Kannan, and Mahoney for approximating matrix multiplication. We state their result in a slightly stronger form, which is actually proved in their paper. For completeness, a proof of this statement is in the appendix. Lemma 4.5 (Matrix multiplication by subsampling [DKM06, Theorem 1]). Suppose we are given X ∈ C n×m , Y ∈ C n×p , r ∈ N and a distribution p ∈ R n satisfying the oversampling condition that, for some ϕ ≥ 1,
Let S ∈ R r×n be sampled according to p. Then X † S † SY is an unbiased estimator for X † Y and
From a simple application of Lemma 4.5, we get a key lemma used frequently in Section 6.
Lemma 4.6 (Approximating matrix multiplication to Frobenius norm error; corollary of [DKM06, Theorem 1]). Consider X ∈ C m×n , Y ∈ C m×p , and take S ∈ R r×m to be sampled according to q := q 1 +q 2 2 , where q 1 , q 2 ∈ R m are ϕ 1 , ϕ 2 -oversampled importance sampling distributions from x, y, the vector of row norms for X, Y , respectively. Then S is a 2ϕ 1 , 2ϕ 2 -oversampled importance sampling sketch of X, Y , respectively. Further,
Proof. First, notice that 2q(i) ≥ q 1 (i) and 2q(i) ≥ q 2 (i), so q oversamples the importance sampling distributions for X and Y with constants 2ϕ 1 and 2ϕ 2 , respectively. We get the bound by using Lemma 4.5; q satisfies the oversampling condition with ϕ = √ ϕ 1 ϕ 2 ∥X∥ F ∥Y ∥ F ∑︁ ℓ ∥X(ℓ,•)∥∥Y (ℓ,•)∥ , using the inequality of arithmetic and geometric means:
Remark 4.7. Lemma 4.6 implies that, given SQ ϕ 1 (X) and SQ ϕ 2 (Y ), we can get SQ ϕ (M ) for M a sufficiently good approximation to
. This is an approximate closure property for oversampling and query access under matrix products.
Given the above types of accesses, we can compute the sketch S necessary for Lemma 4.6 by taking p = D x ˜and q = D y ˜), thereby finding a desired M := X † S † SY . We can compute entries of M with only r queries each to X and Y , so all we need is to get SQ(M ˜) for M ˜the appropriate bound. We choose
showing that we have SQ(M ) follows from the proofs of Lemmas 3.8 and 3.9, since M is simply a linear combination of outer products of rows of X ˜with rows of Y ˜. Finally, this bound has the appropriate norm. Notating the rows sampled by the sketch as s 1 , . . . , s r , we have
If X = Y , we can get an improved spectral norm bound: instead of depending on ∥X∥ 2 F , error depends on ∥X∥∥X∥ F . Lemma 4.8 (Approximating matrix multiplication to spectral norm error [RV07, Theorem 3.1]). Suppose we are given A ∈ R m×n , ε > 0, δ ∈ [0, 1], and S ∈ R r×n a ϕ-oversampled importance sampling sketch of A. Then
The above results can be used to approximate singular values, simply by directly translating the bounds on matrix product error to bounds on singular value error. Lemma 4.9 (Approximating singular values). Given SQ ϕ (A) ∈ C m×n and ε ∈ (0, 1], we can form importance sampling sketches S ∈ R r×m and
)︁ time satisfying the following property. Take r, c ≥ s for some sufficiently large
. Then, if σ i and σ ˆi are the singular values of A and SAT , respectively (where σ ˆi = 0 for i > min(r, c)), we have with probability
If we additionally assume that ε ≲ ∥A∥/∥A∥ F , we can conclude
This result follows from results bounding the error between singular values by errors of matrix products. For notation, let σ i (M ) be the ith largest singular value of M . We will use the following inequalities relating norm error of matrices to error in their singular values:
When A, B are Hermitian, the same bound holds for their eigenvalues. 13Proof of Lemma 4.9. We use known theorems, plugging in the values of r and c. Using Lemma 4.6 for the sketch S, we know that
by Lemma 4.3, T † is an ≤ 2ϕ-oversampled importance sampling sketch of (SA) † , so by Lemma 4.6 for
By rescaling δ and union bounding, we can have all events happen with probability ≥ 1 -δ. Then, from triangle inequality followed by Lemma 4.10,
The analogous result holds for spectral norm via Lemma 4.8 and Lemma 4.11; the only additional complication is that we need to assert that ∥SA∥ ≲ ∥A∥. We use the following argument, using the upper bound on ε:
Finally, if we wish to approximate a vector inner product u † v, a special case of matrix product, we can do so with only sampling and query access to one of the vectors while still getting log 1 δ dependence on failure probability. The proof of this statement is in Appendix B.
Lemma 4.12 (Inner product estimation, [Tan19, Proposition 4.2]).
)︁ .
Remark 4.13. Lemma 4.12 also applies to higher-order tensor inner products:
(a) (Trace inner products, [GLT18, Lemma 11]) Given SQ ϕ (A) ∈ C n×n and Q(B) ∈ C n×n , we can estimate Tr[AB † ] to additive error ε with probability at least 1 -δ by using
time. To do this, note that SQ ϕ (A) and Q(B) imply SQ ϕ (vec(A)) and Q(vec(B)). Tr[AB] = ⟨vec(B), vec(A)⟩, so we can just apply Lemma 4.12 to conclude.
time. To do this, observe that x † Ay = Tr(x † Ay) = Tr(Ayx † ) and that Q(yx † ) can be simulated with Q(x), Q(y). So, we just apply the trace inner product procedure.
Finally, we observe a simple technique to convert importance sampling sketches into approximate isometries, by inserting the appropriate pseudoinverse. This will be used in some of the more involved applications.
Lemma 4.14. Given A ∈ C m×n , S ∈ C r×m sampled from a ϕ-oversampled importance sampling distribution of A, and T † ∈ C n×c sampled from an ≤ ϕ-oversampled importance sampling distribution of (SA) † , let R := SA and C := SAT . Let σ k be the kth singular value of A. If, for α ∈ (0, 1],
Proof. The following occurs with probability ≥ 1 -δ. By Lemma 4.9,
. By Lemma 4.8,
We get the computation for the α-approximate isometry by restricting attention to the span of U :
One can also observe that, for a sufficiently good sketch C, R ≈ C k (C k ) + R in spectral norm, giving a generic way to approximate a sketch R by a product of a small matrix with an approximate projective isometry. We do not need it in our proofs, so this computation is not included.
Our main result is that, given SQ ϕ (A) and a smooth function f , we can approximate f (A † A) by a decomposition R † U R + f (0)I. This primitive is based on the even singular value transformation used by Gilyén, Su, Low, and Wiebe [GSLW19].
Theorem 5.1 (Even singular value transformation). Let A ∈ C m×n and f : R + → C be such that f (x) and f ¯(x) := (f (x)-f (0))/x are L-Lipschitz and L ¯-Lipschitz, respectively, on
i +d] for some d > 0. Take parameters ε and δ such that 0 < ε ≲ min(L∥A∥ 2 * , L ¯∥A∥ 2 * ∥A∥ 2 ) and δ ∈ (0, 1]. Choose a norm * ∈ {F, Op}.
Suppose we have SQ ϕ (A). Consider the importance sampling sketch S ∈ R r×m corresponding to SQ ϕ (A) and the importance sampling sketch T † ∈ R c×n corresponding to SQ ≤2ϕ ((SA) † ) (which we have by Lemma 4.3). Then, for R := SA and C := SAT , we can achieve the bound
First, we make some technical remarks. The assumption that ε ≲ L∥A∥ 2 * is for non-degeneracy: if ε ≥ L∥A∥ 2 , then the naive approximation f (0)I of f (A † A) would suffice, since ∥f (0)I -f (A † A)∥ ≤ L∥A∥ 2 ≤ ε as desired. 14 The parameter d (or, rather, the parameter ε ¯) specifies the domain where 14 The choice f (0)I assumes that f is Lipschitz on {0, ∥A∥ 2 }. More generally, we can choose f (x)I for any
in order to get a sufficiently good naive approximation.
f (x) and f ¯(x) should be smooth: the condition in the theorem is that they should be Lipschitz on the spectrum of A † A, with ε ¯room for approximation. This will not come into play often, though, since we can often design our singular value transforms such that we can take d = ∞. For example, if our desired transform f becomes non-smooth outside the relevant interval [0, ∥A∥ 2 ], we can apply Theorem 5.1 with d = ∞ and the function g such that g(x) = g(∥A∥ 2 ) for x ≥ ∥A∥ 2 and g(x) = f (x) otherwise. Then g(A † A) = f (A † A) and g is smooth everywhere, so we do not need to worry about the d parameter. Finally, we note that no additional log terms are necessary (i.e., Ω ˜becomes Ω) when the Frobenius norm is used.
By our discussion in Section 4, finding the sketches S and T for Theorem 5.1 takes time
, querying for all of the entries of C takes additional time O(rc q(A)), and computing f ¯(CC † ) takes additional time O (︁ min(r 2 c, rc 2 ) )︁ (if done naively). For our applications, this final matrix function computation will dominate the runtime, and the rest of the cost we will treat as O (︁ rc sq ϕ (A) )︁ . For some intuition on error bounds and time complexity, we consider how the parameters in our main theorem behave in a restricted setting: suppose we have SQ(A) with minimum singular value σ and such that ∥A∥ F /σ is dimension-independent. 15 This condition simultaneously bounds the rank and condition number of A. Further supposefoot_11 that f is L-Lipschitz on the interval [0, ∥A∥ 2 ] and satisfies L∥A∥ 2 < ΓD where
for some dimension-independent Γ. Γ must be at least one, so we can think about such an f as being at most Γ times "steeper" compared to the least possible "steepness". Under these assumptions, we can get a decomposition satisfying
The time to compute the decomposition is
These quantities are all dimensionless. Dependence on σ arises because we bound L ¯≤ L/σ 2 : our algorithm's dependence on L ¯implicitly enforces a low-rank constraint in this case. All of our analyses give qualitatively similar results to this, albeit in more general settings allowing approximately low-rank input.
To perform error analyses, we will need bounds on the norms of the matrices in our decomposition. The following lemma gives the bounds we need for Section 6.
Lemma 5.2 (Norm bounds for even singular value transformation). Suppose the assumptions from Theorem 5.1 hold. Then with probability at least 1 -δ, the event in Eq. (3) occurs (that is,
and moreover, the following bounds also hold:
Eq. ( 7) is typically a better bound than combining Eqs. ( 5) and (6). For intuition, notice this is true if ε, ε ¯= 0: the left-hand and right-hand sides of the following inequality are the two ways to bound ∥R † √︁ f ¯(CC † )∥ 2 , up to constant factors (σ below runs over the singular values of A):
The rest of this section will be devoted to proving Theorem 5.1 and Lemma 5.
Proof of Theorem 5.1 and Lemma 5.2. Since g(A † A) = f (A † A) + g(0)I for f (x) := g(x) -g(0), we can assume without loss of generality that f (0) = 0. As a reminder, in the statement of Theorem 5.1 we take
These values are chosen such that the following holds with probability ≥ 1 -δ simultaneously.
1. The ith singular value of CC † does not differ from the ith singular value of A † A by more than ε ¯. This follows from Lemma 4.9 with error parameter ε∥A∥ -1
). This immediately implies Eq. (6).
. This is the spectral norm bound in Eq. (5) (the Frobenius norm bound follows from Lemma 4.2). We use Lemma 4.8:
We need the polylog factors in our number of samples to deal with the log r that arises from Lemma 5.4 in the spectral norm case.
(plugging in value for r)
. This follows similarly to the above point.
When all of the above bounds hold, we can conclude:
This gives Eq. (3) after rescaling ε by an appropriate constant factor. When * = Op, we also have Eq. ( 7), since
We remark here that the log term in Lemma 5.4 unfortunately cannot be removed (because some Lipschitz functions are not operator Lipschitz). However, several bounds hold under various mild assumptions, and for particular functions, the log term can be improved to log log or completely removed. For example, the QSVT literature [GSLW19] cites the following result:
.
Though we will not use it, we can extend these results on eigenvalue transformation of Hermitian matrices to singular value transformation of general matrices via the reduction from [GSLW19, Corollary 21]. For example, Lemma 5.4 implies the following: Lemma 5.6. Let A, B ∈ C m×n be matrices and let f : [0, ∞) → C be L-Lipschitz continuous on the singular values of A and B such that f (0) = 0. Then
In Section 7, we prove results on generic singular value transformation and eigenvalue transformation by bootstrapping Theorem 5.1. Since these are slower, though, we will use primarily the even singular value transformation results that we just proved to recover "dequantized QML" results. This will be the focus of next section.
6 Applying the framework to dequantizing QML algorithms Now, with our framework, we can recover previous dequantization results: recommendation systems (Section 6.2), supervised clustering (Section 6.3), principal component analysis (Section 6.4), lowrank matrix inversion (Section 6.5), support-vector machines (Section 6.6), and low-rank semidefinite programs (Section 6.8). We also propose new quantum-inspired algorithm for other applications, including QSVT (Section 6.1), Hamiltonian simulation (Section 6.7), and discriminant analysis (Section 6.9). We give applications in roughly chronological order; this also happens to be a rough difficulty curve, with applications that follow more easily from our main results being first.
Everywhere it occurs, K := ∥A∥ 2 F /σ 2 , where A is the input matrix. κ := ∥A∥ 2 2 /σ 2 . For simplicity, we will often describe our runtimes as if we know spectral norms of input matrices (so, for example, we know κ). If we do not know the spectral norm, we can run Lemma 4.9 repeatedly with multiplicatively decreasing ε until we find a constant factor upper bound on the spectral norm, which suffices for our purposes. Alternatively, we can bound the spectral norm by the Frobenius norm, which we know from sampling and query access to input.
We begin by dequantizing the quantum singular value transformation described by Gilyén, Su, Low, and Wiebe [GSLW19] for close-to-low-rank matrices. Definition 6.1. For a matrix A ∈ C m×n and p(x) ∈ C[x] degree-d polynomial of parity-d (i.e., even if d is even and odd if d is odd), we define the notation p (QV) (A) in the following way:
1. If p is even, meaning that we can express p(x) = q(x 2 ) for some polynomial q(x), then
2. If p is odd, meaning that we can express p(x) = x • q(x 2 ) for some polynomial q(x), then Then with probability ≥ 1 -δ, for ε a sufficiently small constant, we can get
)︂ time. Specifically, for p even, the runtime is
and for p odd, the runtime is
.
From this result it follows that QSVT, as described in [GSLW19, Theorem 17], has no exponential speedup when the block-encoding of A comes from a quantum-accessible "QRAM" data structure as in [GSLW19, Lemma 50]. In the setting of QSVT, given A and b in QRAM, one can prepare |b⟩ and construct a block-encoding for A/∥A∥ F = A in polylog(mn) time. Then one can apply (quantum) SVT by a degree-d polynomial on A and apply the resulting map to |b⟩ with d•polylog(mn) gates and finally project down to get the state |p (QV) (A)b⟩ with probability ≥ 1 -δ after Θ (︁ ∥p (QV) (A)b∥ , and log(mn). We are exponentially slower in ε and δ (these errors are conflated for the quantum algorithm). However, this exponential advantage vanishes if the desired output is not a quantum state but some fixed value (or an estimate of one). In that case, the quantum algorithm must also pay 1 ε during the sampling or tomography procedures and the classical algorithm can boost a constant success probability to ≥ 1 -δ, only paying a log 1 δ factor. Note that, unlike in the quantum output, we can query entries of the output, which a quantum algorithm cannot do without paying at least a 1 ε factor. Theorem 6.2 also dequantizes QSVT for block-encodings of density operators when the density operator comes from some well-structured classical data. Indeed, [GSLW19, Lemma 45] assumes we can efficiently prepare a purification of the density operator ρ. The rough classical analogue is the assumption that we have sampling and query access to some A ∈ C m×n with ρ = A † A. Since Tr(ρ) = 1, we have ∥A∥ F = 1. Then, p (QV) (ρ) = r (QV) (A) for r(x) = p(x 2 ) and ∥ρ∥ = ∥A∥ 2 , so we can repeat the above argument to show the lack of exponential speedup for this input model too.
We can mimic the quantum algorithm with our techniques because low-degree polynomials are smooth, in the sense that we formalize with the following lemma (proven in Appendix B).
These bounds are tight for Chebyshev polynomials. In general, these bounds can be loose, so for any particular QML application we recommend using our main results for faster algorithms.
Proof of Theorem 6.2. Consider the even case: take p(x) = q(x 2 ) for q a degree-d/2 polynomial, so p (QV) (A) = q(A † A), and we have the correct form to apply Theorem 5.1. q is uncontrolled outside of [-1, 1], so we instead apply the singular value transformation which is constant outside of [-1, 1]:
We can do this because the singular values of A lie in [0, 1], so q(A † A) = f (A † A). Then, by Lemma 6.3, f and f ¯are Lipschitz with
(We will later rescale ε; note that ε∥p (QV) (A)b∥ ≲ L∥A∥ 2 ∥b∥, so ε is small enough for the theorem assumption.) This reduces the problem to approximating
, using that ∥R∥ F ≲ ∥A∥ F (Eq. ( 5)) and sq(R † ) = O(r sq(A)) (Lemma 4.3). This suffices to maintain the error bound because (using Eqs. ( 6) and ( 7) and Lemma 6.3),
In the last step, we use that ε ≲ ∥p (QV) (A)b∥; if we do not have that assumption, ˜︂ sq(v) ≲ r 2 d 4 ∥v∥ 2 log(mn) log 1 δ . We rescale ε ← ε∥p (QV) (A)b∥ to get the desired bound. The runtime is dominated by finding
We also need to compute the matrix-vector product f ¯(CC † )u, but this can be done in O(rc) time by instead multiplying through with the expression
where U ∈ C r×c comes from the SVD of C. Now for the odd case: we could use Theorem 7.1 here, but we will continue to use Theorem 5.1 here. Similarly to the even case, we take g(x) to be q(x) in [-1, 1] and held constant ouside it, so p (QV) (A) = A • g(A † A). Then, by plugging in the smoothness parameters from Lemma 6.3, we get R, C such that ∥R † g ¯(CC † )R + g(0)I -g(A † A)∥ < ε ∥A∥ with probability ≥ 1 -δ where
We now use the approximating matrix product lemmas Lemmas 4.4 and 4.6 three times.
. We can do this since we have SQ(A † ) (by assumption) and SQ
samples, with each sample costing O(r log(mn)) time (by Lemma 4.3). Then using Lemma 5.2 and the bounds on |q(x)| and |q ¯(x)| in Lemma 6.3,
2. We use Lemma 4.6 to approximate Rb ≈ u such that ∥Rb -u∥ ≤ ε d 2 ∥A∥ , where we use
Above, we used that ε ≲ ∥p (QV) (A)b∥ ≤ ∥A∥∥q(A † A)∥∥b∥ ≤ d∥A∥. Now, we rescale ε ← ε∥p (QV) (A)b∥ to get the desired statement. The runtime is dominated by the sampling for C in O(rc log(mn)) time, the computation of g ¯(CC † ) in O (︁ r 2 c )︁ time, and the approximation of
Remark 6.4. Here, we make a brief remark about a technical detail we previously elided. Technically, QSVT can use A † in QRAM instead of A (cf. [GSLW19, Lemma 50]), leaving open the possibility that there is a quantum algorithm that does not give an exponential speedup when A is in QRAM, but does when A † is in QRAM. We sketch an argument why this is impossible by showing that, given SQ(A), we can simulate SQ ϕ (B) (and SQ ϕ (B † )) for B such that ∥B -A † ∥ ≤ ε∥A∥ with probability ≥ 1 -δ. Unfortunately, this argument is fairly involved, so we defer it to Appendix A.
Our framework gives a simpler and faster variant of Tang's dequantization [Tan19] of Kerenidis and Prakash's quantum recommendation systems [KP17]. Tang's result is notable for being the first result in this line of work and for dequantizing what was previously believed to be the strongest candidate for practical exponential quantum speedups for a machine learning problem [Pre18].
We want to find a product j ∈ [n] that is a good recommendation for a particular user i ∈ [m], given incomplete data on user-product preferences. If we store this data in a matrix A ∈ R m×n with sampling and query access, in the strong model described by Kerenidis and Prakash [KP17], finding good recommendations reduces to the following: Problem 6.5. For a matrix A ∈ R m×n , given SQ(A) and a row index i ∈ [m], sample from A ˆ(i, •) up to δ error in total variation distance, where ∥A ˆ-A σ,η ∥ F ≤ ε∥A∥ F .
Here, A σ,η is a certain type of low-rank approximation to A. The standard notion of low-rank approximation is that of
, which is the rank-r matrix closest to A in spectral and Frobenius norms. Using singular value transformation, we define an analogous notion thresholding singular values instead of rank. Definition 6.6 (A σ,η ). We define A σ,η as a singular value transform of A satisfying:
Note that P σ,η is not fully specified in the range [σ(1 -η), σ(1 + η)), so A σ,η is any of a family of matrices with error η.
σ,η (A) is A for (right) singular vectors with value ≥ σ(1 + η), zero for those with value < σ(1 -η), and something in between for the rest. Our analysis simplifies the original algorithm, which passes through the low-rank approximation guarantee of Frieze, Kannan, and Vempala [FKV04].
Our algorithm uses that we can rewrite our target low-rank approximation as A • t(A † A), where t is a smoothened projector. So, we can use our main theorem, Theorem 5.1, to approximate t(A † A) by some R † U R. Then, the i th row of our low-rank approximation is A(i, •)R † U R, which is a product of a vector with an RUR decomposition. Thus, using the sampling techniques described in Section 4.1, we have SQ ϕ (A(i, •)R † U R), so we can get the sample from this row as desired.
Corollary 6.7. Suppose 0 < ε ≲ ∥A∥/∥A∥ F and η ≤ 0.99. A classical algorithm can solve Problem 6.5 in time
The assumption on ε is a weak non-degeneracy condition in the low-rank regime. For reference, η = 1/6 in the application of this algorithm to recommendation systems. So, supposing the first term of the runtime dominates, the runtime is ˜︁ O (︂ ∥A∥ 6 F ∥A∥ 10 σ 16 ε 6 Proof. Note that A σ,η = A • t(A † A), where t is the thresholding function shown below.
We will apply Theorem 5.1 with error parameter ε to get matrices R, C such that
Since t(x) is (4ησ 2 ) -1 -Lipschitz and t(x)/x is (4η(1 -η) 2 σ 4 ) -1 -Lipschitz, the sizes of r and c are
So, it suffices to compute the SVD of an r × c matrix, which has a runtime of
Next, we want to approximate AR † ≈ A ′ R ′ † . If we had SQ(A † ) (in particular, if we could compute column norms ∥A(•, j)∥), we could do this via Lemma 4.6, and if we were okay with paying factors of 1 δ , we could do this via Lemma 4.4. Here, we will instead implicitly define an approximation by approximating each row [AR † ](i, •) = A(i, •)R † via Lemma 4.6, since we then have SQ(A(i, •) † ) and SQ(R † ). With this proposition, we can estimate [AR † ](i,
Using this and the observation that max x t ¯(x) = (1 + η) -2 σ -2 ≤ σ -2 , we can bound the quality of our final approximation as
by Lemma 5.2 and Eq. ( 9)
We can sample from A ˆ(i,
taking O(r ′ r + rc) time. Then, we use Lemmas 3.5 and 3.6 to get a sample from xR with probability
, where sq ϕ (xR) = O(r) and
Then, using previously established bounds and bounds from Lemma 5.2, we have
This sampling procedure and the SVD dominate the runtime. Since the sampling is exact, the only error in total variation distance is the probability of failure.
Remark 6.8. This algorithm implicitly assumes that the important singular values are ≥ σ.
Without such an assumption, we can take σ = ε∥A∥ F and η = 1/2, and have meaningful bounds on the output matrix A ˆ. Observe that, for p(x) = x(t( √ x) -1),
So, our low-rank approximation output A ˆsatisfies ∥A ˆ-A∥ ≲ ε∥A∥ F , with no assumptions on A,
This can be subsequently used to get SQ ϕ (A ˆ(i, •)) = SQ ϕ (e i A ˆ) where
in a myopic sense, solving the same problem as Problem 6.5), or more generally, any product of A ˆwith a vector, in time independent of dimension.
The 2013 paper of Lloyd, Mohseni, and Rebentrost [LMR13] gives two algorithms for the machine learning problem of clustering. The first algorithm is a simple swap test procedure that was dequantized by Tang [Tan21] (the second is an application of the quantum adiabatic algorithm with no proven runtime guarantees). We will reproduce the algorithm from [Tan21] here: since the dequantization just uses the inner product protocol, so it rather trivially fits into our framework.
We have a dataset of points in R d grouped into clusters, and we wish to classify a new data point by assigning it to the cluster with the nearest average, aka centroid. We do this by estimating the distance between the new point p ∈ R d to the centroid of a cluster of points q 1 , . . . , q n-1 ∈ R d , namely, ∥p -1 n-1 (q 1 + • • • + q n-1 )∥ 2 . This is equal to ∥wM ∥ 2 , where
. . .
Because the quantum algorithm assumes input in quantum states, we can assume sampling and query access to the data points, giving the problem Problem 6.9. Given SQ(M ) ∈ R n×d , Q(w) ∈ R n , approximate (wM )(wM ) T to additive ε error with probability at least 1 -δ.
There is a classical algorithm to solve Problem 6.9 in
Note that ∥M ∥ 2 F = 2 and ∥w∥ 2 = ∥p∥ 2 + 1 n-1
∑︁ n-1 i=1 ∥q i ∥ 2 . The quantum algorithm has a quadratically faster runtime of O
Proof. Recall our notation for the vector of row norms m := [∥M (1, •)∥, . . . , ∥M (n, •)∥] coming from Definition 3.7. We can rewrite (wM )(wM ) T as an inner product ⟨u, v⟩ where
where u and v are three-dimension tensors. By flattening u and v, we can represent them as two vectors in R (n•d•n)×1 . We clearly have Q(v) from queries to M and w. As for getting SQ(u) from SQ(M ): to sample, we first sample i according to m, sample j according to M (i, •), and sample k according to m; to query, compute u i,j,k = M (i, j)m(k). Finally, we can apply Lemma 4.12 to estimate ⟨u, v⟩. ∥u∥ = ∥M ∥ 2 F and ∥v∥ = ∥w∥ 2 , so estimating ⟨u, v⟩ to ε additive error with probability at least 1 -δ requires O(∥M ∥ 4 F ∥w∥ 4 ε -2 log 1 δ ) samples.
Principal component analysis (PCA) is an important data analysis tool, first proposed to be feasible via quantum computation by Lloyd, Mohseni, and Rebentrost [LMR14]. Given copies of states with density matrix ρ = X † X, the quantum PCA algorithm can prepare the state ∑︁ λ i |v i ⟩⟨v i | ⊗ |λ ˆi⟩⟨λ ˆi|, where λ i and v i are the eigenvalues and eigenvectors of X † X, and λ ˆi are eigenvalue estimates (up to additive error). See Prakash's PhD thesis [Pra14, Section 3.2] for a full analysis and Chakraborty, Gilyén, and Jeffery for a faster version of this algorithm in the block-encoding model [CGJ19]. Directly measuring the eigenvalue register is called spectral sampling, but such sampling is not directly useful for machine learning applications.
Though we do not know how to dequantize this protocol exactly, we can dequantize it in the low-rank setting, which is the only useful poly-logarithmic time application that Lloyd, Mohseni, and Rebentrost [LMR14] suggests for quantum PCA. Problem 6.11 (PCA for low-rank matrices). Given a matrix SQ(X) ∈ C m×n such that X † X has top k eigenvalues {λ i } k i=1 and eigenvectors {v i } k i=1 , with probability ≥ 1 -δ, compute eigenvalue estimates {λ ˆi} k i=1 such that
Note that we should think of λ i as σ 2 i , where σ i is the ith largest singular value of X. To robustly avoid degeneracy conditions, our runtime must depend on parameters for condition number and spectral gap:
We also denote κ := ∥X∥ 2 /λ k . Dependence on K and η are necessary to reduce Problem 6.11 to spectral sampling. If K = poly(n), then λ k = Tr(X † X)/ poly(n), so distinguishing λ k from λ k+1 necessarily takes poly(n) samples, and even sampling λ k once takes poly(n) samples. As a result, learning v k is also impossible. A straightforward coupon collector argument (given e.g. by Tang [Tan21]) shows that Problem 6.11 can be solved by a quantum algorithm performing spectral samplingfoot_13 , with runtime depending polynomially on K and 1 η . We omit this argument for brevity.
Classically, we can solve this PCA problem with quantum-inspired techniques, as first noted in [Tan21].
Corollary 6.12. For 0 < ε ≲ η∥X∥ 2 /∥X∥ 2 F , we can solve Problem 6.11 in ˜︁ O (︂
This improves significantly over prior work [Tan21,Theorem 8], which achieves the runtime of
. 19 The best quantum algorithm for this problem runs in ˜︁ O (︂
We approach the problem as follows. First, we use that an importance-sampled submatrix of X has approximately the same singular values as X itself (Lemma 4.9) to get our estimates {λ ˆi} k i=1 . With these estimates, we can define smoothened step functions
We can then use our main theorem to find an RUR decomposition for f i (X † X). We use additional properties of the RUR description to argue that it is indeed a rank-1 outer product v ˆ † i v ˆi, which is our desired approximation for the eigenvector. We have sampling and query access to v ˆi because it is R † x for some vector x. Our runtime is quite good because these piecewise linear step functions have relatively tame derivatives, as opposed to the thresholded inverse function, whose Lipschitz constants must incur quadratic and cubic overheads in terms of condition number.
Proof. We will assume that we know λ k and η. If both are unknown, then we can estimate them with the singular value estimation procedure described below (Lemma 4.9).
Notice that η∥X∥ 2 ≤ λ k follows from our definition of η. The algorithm will proceed as follows: first, consider C := SXT ∈ C c×r as described in Theorem 5.1, with parameters
Consider computing the eigenvalues of CC † ; denote the i th eigenvalue λ ˆi. Since r, c = Ω(
, by Lemma 4.9 with error parameter ε
These λ ˆi's for i ∈ [k] have the desired property for eigenvalue estimates:
19 This runtime comes from taking εσ = εv = ε and changing the normalization of the gap parameter η = η∥X∥ 2 /∥X∥ 2 F to correspond to the problem as formulated here. 20 Given X in QRAM, this follows from applying Theorem 27 to a quantum state with density matrix of X † X with α = ∥X∥F and ∆ = ε∥X∥ 2 F ∥X∥ ≲ η∥X∥ ∥X∥ F . The output is some estimate of √ λi to ∆ error, which when squared is an estimate of λi to ∆∥X∥ = ε∥X∥ 2 F error as desired. Then, the density matrix is a probability distribution over eigenvectors with their corresponding eigenvalue estimate (which is enough to identify the eigenvector). The coupon collector argument mentioned above gives us access to all the top k eigenvalues and eigenvectors by running this algorithm ∥X∥ 2 F /λ k times [Tan21].
This bound also implies that, for all i, |λ ˆi -
This is a function that is one when |x-λ ˆi| ≤ 1 8 η∥X∥ 2 , zero when |x-λ ˆi| ≥ 1 4 η∥X∥ 2 , and interpolates between them otherwise. From the eigenvalue gap and the aforementioned bound
i exactly. Further, by Theorem 5.1, we can conclude that
i , with C, R the exact approximations used to estimate singular values. The conditions of Theorem 5.1 are satisfied because ε ≲ 8 ≤ 8 η = L∥X∥ 2 for L the Lipschitz constant of f i . The values of r, c are chosen so that
Further, f i is chosen with respect to λ ˆi such that R † f ¯i(CC † )R is rank one, since CC † has one eigenvalue between λ ˆi -1 4 η∥X∥ 2 and λ ˆi + 1 4 η∥X∥ 2 . Thus, this approximation is an outer product, R † f ¯i(CC † )R = v ˆiv ˆ † i , and we take the corresponding vector to be our eigenvector estimate:
which is the desired bound. By choosing failure probability δ/k, the bound can hold true for all k with probability ≥ 1 -δ. Finally, we can get access to
i , using Lemmas 3.5 and 3.6, we have SQ ϕ (v ˆi) with
The low-rank matrix inversion algorithms given by Gilyén, Lloyd, and Tang [GLT18] and Chia, Lin, and Wang [CLW18] dequantize Harrow, Hassidim, and Lloyd's quantum matrix inversion algorithm (HHL) [HHL09] in the regime where the input matrix is low-rank instead of sparse. The corresponding quantum algorithm in this regime is given by Chakraborty, Gilyén, and Jeffery [CGJ19], among others. Since sparse matrix inversion is BQP-complete, it is unlikely that one can efficiently dequantize it. However, the variant of low-rank (non-sparse) matrix inversion appears often in quantum machine learning [Pra14, WZP18, RML14, CD16, RL18], making it an influential primitive in its own right.
Using our framework, we can elegantly derive the low-rank matrix inversion algorithm in a manner similar to prior quantum-inspired work [GLT18,CLW18]. Moreover, we can also handle the approximately low-rank regime and only invert the matrix on a well-conditioned subspace, solving principal component regression-for more discussion see [GSLW19]. Namely, we can find a thresholded pseudoinverse of an input matrix: Definition 6. 13 (A + σ,η ). We define A + σ,η to be any singular value transform of A satisfying:
This definition is analogous to A σ,η in Section 6.2: it is A + for singular vectors with value ≥ σ, zero for singular vectors with value ≤ σ(1 -η), and a linear interpolation between the two in between.
Problem 6.14. Given SQ φ (A) ∈ C m×n , Q(b) ∈ C m , with probability ≥ 1 -δ, get SQ ϕ (x ˆ) such that ∥x ˆ-x * ∥ ≤ ε∥A∥ -1 ∥b∥, where x * := A + σ,η b.
Corollary 6.15. For 0 < ε ≲ ∥A∥ 2 σ 2 and η ≤ 0.99, we can solve Problem 6.14 in ˜︁ O (︂
This should be compared to [GLT18], which applies only to strictly rank-k A with φ = 1 and gets the incomparable runtime of ˜︁ O (︂
. The corresponding quantum algorithm using block-encodings takes O(∥A∥ F /σ) time, up to polylog(m, n) factors, to get this result for constant η [GSLW19, Theorem 41].
If we further assume that ε < 0.99 and b is in the image of A, then ˜︂ sq ϕ (x ˆ) can be simplified,
However, this algorithm also works for larger ε; namely, if we only require that ∥x ˆ-x * ∥ ≤ εσ -1 ∥b∥ (a "worst-case" error bound), then this algorithm works with runtime smaller by a factor of κ 3 (and ˜︂ sq ϕ (x ˆ) smaller by a factor of κ). The algorithm comes from rewriting A + σ,η b = ι(A † A)A † b for ι a function encoding a thresholded inverse. Namely, ι(x) = 1/x for x ≥ σ 2 , ι(x) = 0 for x ≤ (1 -η) 2 σ 2 , and is a linear interpolation between the endpoints for x ∈ [(1 -η) 2 σ 2 , σ 2 ]. By our main theorem, we can find an RUR decomposition for ι(A † A), from which we can then get SQ(R † U RA † b) via sampling techniques.
Proof. We will solve our problem for
By Theorem 5.1 with L = 1 (2η-η 2 )σ 4 and L ¯= 1 (1-η) 2 (2η-η 2 )σ 6 , we can find such R and C with
Computing the SVD of a matrix of this size dominates the runtime, giving the complexity in the theorem statement. Next, we would like to further approximate R † ι ¯(CC † )RA † b. We will do this by estimating RA † b by some vector u to εσ 3 ∥A∥ -1 ∥b∥ = ε∥A∥ 2 F ∥b∥K -1 κ -1 2 error, since then, using the bounds from Lemma 5.2,
We use Remark 4.13 to estimate
samples for each of the r entries. This implies that x ˆ:= R † ι ¯(CC † )u has the desired error and failure probability. Finally, we can use Lemmas 3.5 and 3.6 with matrix R † and vector ι ¯(CC † )u to get SQ ϕ (x ˆ) for
In this section, we use our framework to dequantize Rebentrost, Mohseni, and Lloyd's quantum support vector machine [RML14], which was previously noted to be possible by Ding, Bao, and
Huang [DBH21]. Mathematically, the support vector machine is a simple machine learning model attempting to label points in R m as +1 or -1. Given input data points x 1 , . . . , x m ∈ R n and their corresponding labels y ∈ {±1} m . Let w ∈ R n and b ∈ R be the specification of hyperplanes separating these points. It is possible that no such hyperplane satisfies all the constraints. To resolve this, we add a slack vector e ∈ R m such that e(j) ≥ 0 for j ∈ [m]. We want to minimize the squared norm of the residuals:
The dual of this problem is to maximize over the Karush-Kuhn-Tucker multipliers of a Lagrange function, taking partial derivatives of which yields the linear system
where 1 ⃗ is the all-ones vector and X = {x 1 , . . . , x m } ∈ C m×n . Call the above m + 1 × m + 1 matrix F , and F ˆ:= F/ Tr(F ).
The quantum algorithm, given X and y in QRAM, outputs a quantum state |F ˆ+ λ,0.01
)︁ time. The quantum-inspired analogue is as follows.
Problem 6.16. Given SQ(X) ∈ R m×n and SQ(y
Note that we must assume ∥F ˆ∥ ≤ 1; the quantum algorithm makes the same assumptionfoot_14 . Another dequantization was reported in [DBH21], which, assuming X is strictly low-rank (with minimum singular value σ), outputs a description of (XX T ) + y that can be used to classify points. This can be done neatly in our framework: express (XX T ) + (or, more generally, (XX T ) + σ,η ) as Xf (X T X)X T for the appropriate choice of f . Then, use Theorem 5.1 to approximate f (X T X) ≈ R T ZR and use Lemma 4.4 to approximate XR T ≈ CW T . This gives an approximate "CUC" decomposition of the desired matrix, since Xf (X T X)X T ≈ XR T ZRX T ≈ CW T ZW C T , which we can use for whatever purpose we like.
For our solution to Problem 6.16, though, we simply reduce to matrix inversion as described in Section 6.5: we first get SQ ϕ (F ˆ), and then we apply Corollary 6.15 to complete. Section VI.C of [DBH21] claims to dequantize this version, but gives no correctness boundsfoot_15 or runtime bounds (beyond arguing it is polynomial in the desired parameters). Corollary 6.17. For 0 < ε ≲ 1 and η ≤ 0.99, we can solve Problem 6.16 in
time, where we get
)︁ .
The runtimes in the statement are not particularly tight, but we chose the form to mirror the runtime of the QSVM algorithm, which similarly depends polynomially on 1 λ and 1 η .
Proof. Consider constructing SQ φ (K) ∈ C m×m as follows. To query an entry K(i, j), we estimate X(i, •)X(j, •) T to ε∥X(i, •)∥∥X(j, •)∥ error. We define K(i, j) to be this estimate. Using Lemma 4.12, we can do this in O (︁ 1 ε 2 log q δ )︁ time. q here refers to the number of times the query oracle is used, so in total the subsequent algorithm will only have an errant query with probability ≥ 1 -δ. (q will not appear in the runtime because it's folded into a polylog term.) Then, we can take K ˜:= xx T , where x ∈ R m is the vector of row norms of X, since by Cauchy-Schwarz,
Since we have SQ(x) from SQ(X), we have SQ(K ˜) with sq(K ˜) = O(1) by Lemma 3.8. ∥K ˜∥2
. We can trivially get SQ(L) for L := [︁ 0 1
with sq(L) = O(1). Our approximation to F ˆis
Using Lemma 3.9, we have SQ φ ′ (M ) with
where the last inequality uses that Tr(F ) ≥ √ m, which follows from ∥F ˆ∥ ≤ 1:
.
Note that we can compute Tr(F ) given SQ(X). So, applying Corollary 6.15, we can get the desired SQ ϕ (v) in runtime
Here, we used that ∥M ∥ ≤ ∥M ∥ F ≲ 1, which we know since φ ′ ≥ 1 (by our definition of oversampling and query access). That Q
does not affect the runtime, since the dominating cost is still the SVD. On the other hand, this does come into play for the runtime for sampling:
We take q = m to guarantee that all future queries will be correct with probability ≥ 1 -δ.
The normalization used by the quantum and quantum-inspired SVM algorithms means that these algorithms fail when X has too small Frobenius norm, since then the singular values from XX T are all filtered out. In Appendix B, we describe an alternative method that relies less on normalization assumptions, instead simply computing F + . This is possible if we depend on ∥X∥ 2 F γ in the runtime. Recall from Eq. ( 12) that we regularize by adding γ -1 I, so γ -1 acts as a singular value lower bound and ∥X∥ 2 F γ implicitly constrains. Corollary 6.18. Given SQ(X T ) and SQ(y), with probability ≥ 1 -δ, we can output a real number b ˆsuch that |b -b ˆ| ≤ ε(1 + b) and SQ ϕ (α ˆ) such that ∥α ˆ-α∥ ≤ εγ∥y∥, where α and b come from Eq. (12). Our algorithm runs in ˜︁
. Note that when γ -1/2 is chosen to be sufficiently large (e.g. O(∥X∥ F )) and ∥α∥ = Ω(γ∥y∥), this runtime is dimension-independent.
Notice that εγ∥y∥ is the right notion, since γ is an upper bound on the spectral norm of the inverse of the matrix in Eq. ( 12). We assume SQ(X T ) instead of SQ(X) for convenience, though both are possible via the observation that f (XX T ) = Xf ¯(X T X)X T .
The problem of simulating the dynamics of quantum systems was the original motivation for quantum computers proposed by Feynman [Fey82]. Specifically, given a Hamiltonian H, a quantum state |ψ⟩, a time t > 0, and a desired error ε > 0, we ask to prepare a quantum state |ψ t ⟩ such that
This problem, known as Hamiltonian simulation, sees wide application, including in quantum physics and quantum chemistry. A rich literature has developed on quantum algorithms for Hamiltonian simulation [Llo96, ATS03, BCK15], with an optimal quantum algorithm for simulating sparse Hamiltonians given in [LC17]. In this subsection, we apply our framework to develop classical algorithms for Hamiltonian simulation. Specifically, we ask: Problem 6.19. Consider a Hermitian matrix H ∈ C n×n , a unit vector b ∈ C n , and error parameters ε, δ > 0. Given SQ(H) and SQ(b), output SQ ϕ (b ˆ) with probability ≥ 1 -δ for some b ˆ∈ C n satisfying ∥b ˆe iH b∥ ≤ ε.
We give two algorithms that are fundamentally the same, but operate in different regimes: the first works for low-rank H, and the second for arbitrary H. Corollary 6.20. Suppose H has minimum singular value σ and ε < min(0.5, σ). We can solve Problem 6.19 in ˜︁ O (︂
This runtime is dimensionless in a certain sense. The natural error bound to require is that ∥b ˆe iH b∥ ≤ ε∥H∥, since
, which is dimensionless. The runtime of the algorithm in the following corollary does not have this property, so its scaling with ∥H∥ is worse, despite being faster for, say, ∥H∥ = 1. Corollary 6.21. For ε < min(0.5, ∥H∥ 3 ), we can solve Problem 6.19 in ˜︁ O (︁
Our strategy proceeds as follows: consider a generic function f (x) and Hermitian H. We can write f (x) as a sum of an even function a(x) := 1 2 (f (x) + f (-x)) and an odd function b(x) := 1 2 (f (x) -f (-x)). For the even function, we can use Theorem 5.1 to approximate it via the function f a (x) := a( √ x); the odd function can be written as H times an even function, which we approximate using Theorem 5.1 for
the Lipschitz constants don't blow up by splitting f into even and odd parts. Now, we specialize to Hamiltonian simulation. We first rewrite the problem, using the function sinc(x) := sin(x)/x.
where f cos (λ) := cos( √ λ) and f sinc (λ) := i • sinc( √ λ). When applying Theorem 5.1 on f cos and f sinc , we will use the following bounds on the smoothness of f cos and f sinc .
We separate these bounds into the case where x ≥ 1, which we use when we assume H has a minimum singular value, and the case where x < 1, which we use for arbitrary H.
Proof of Corollary 6.21. Using the Lipschitz bounds above with Theorem 5.1, we can find
where, using that our Lipschitz constants are all bounded by constants,
As a consequence,
Note that, by Lemma 5.2, ∥R cos ∥ ≲ ∥H∥, ∥f ¯cos (C cos C † cos )∥ ≲ 1, and ∥R † cos √︂ f ¯cos (C cos C † cos )∥ ≲ 1; the same bounds hold for the sinc analogues. We now approximate using Lemma 4.6 four times.
4. We approximate Hb ≈ R † w to ε∥b∥ accuracy, requiring r
Our output will be b
ˆas a linear combination of a small number of vectors, all of which we have sampling and query access to. We can complete using Lemmas 3.5 and 3.6, where the matrix is the concatenation
, and the vector is the concatenation
. The length of this vector is r cos + 1 + r sinc + r ≲ r sinc . We get SQ ϕ (b ˆ) where
In the second inequality, we use the same bounds for proving ∥b ˆe iH b∥ ≤ ε, repurposed to argue that all approximations are sufficiently close to the values they are estimating, up to relative error.
Proof of Corollary 6.20. Our approach is the same, though with different parameters. For Theorem 5.1, we use that in the interval [σ 2 /2, ∞), f cos has Lipschitz constants of L = O(1/σ) and L ¯= O(1/σ 3 ) and f sinc has L = O(1/σ 2 ) and L ¯= O(1/σ 4 ). So, if we take
all the conditions of Theorem 5.1 are satisfied: in particular, σ 2 /2 > ε ¯in both cases, up to rescaling ε by a constant factor:
Here, we used our initial assumption that ε ≤ σ. So, the bounds Eqs. ( 13) and ( 14) hold. Note that, by Lemma 5.2, ∥R cos ∥ ≲ ∥H∥, ∥f ¯cos (C cos C † cos )∥ ≲ σ -2 , and ∥R † cos √︂ f ¯cos (C cos C † cos )∥ ≤ 1; the same bounds hold for the sinc analogues. We now approximate using Lemma 4.6 four times.
4. We approximate Hb ≈ R † w to ε∥b∥ accuracy, requiring r
Our output will be
which is close to e iH b by the argument
Now, we have expressed b ˆas a linear combination of a small number of vectors, all of which we have sampling and query access to. We can complete using Lemmas 3.5 and 3.6, where the matrix is the concatenation
, and the vector is the concatenation
. The length of this vector is r cos + 1 + r sinc + r ≲ r sinc . We get SQ ϕ (b ˆ) where
. Since ε < σ, the r 2 sinc term dominates.
Remark 6.22. In the case where H is not low-rank, we could still run a modified version of Corollary 6.20 to compute a modified "exp σ,η (iH)" where singular values below σ are smoothly thresholded away. Following the same logic as Definition 6.13, we could redefine f cos such that f cos (x) = 1 for x < σ 2 (1 -η), f cos (x) = cos( √ λ) for x ≥ σ 2 , and is a linear interpolation between the endpoints for the x in between (and f sinc similarly). These functions have the same Lipschitz constants as their originals, up to factors of 1 η , and give the desired behavior of "smoothing away" small singular values (though we do keep the 0th and 1st order terms of the exponential). Remark 6.23. Our result generalizes those of Ref. [RWC + 20], which achieves essentially the same result only in the much easier regime where H and b are sparse. They achieve a significant speedup due to these assumptions: note that when H is sparse, and a subsample of rows R is taken, RR † can be computed in time independent of dimension; so, we only need to take a subsample of rows, and not of columns. More corners can be cut from our algorithm in this fashion. In summary, though our algorithm is significantly slower, their sparsity assumptions are essential for their fast runtime, and our framework can identify where these tradeoffs occur.
A recent line of inquiry in quantum computing [BS17, vAGGdW20, BKL + 19, vAG19] focuses on finding quantum speedups for semidefinite programs (SDPs), a central topic in the theory of convex optimization with applications in algorithms design, operations research, and approximation algorithms. Chia, Li, Lin, and Wang [CLLW20] first noticed that quantum-inspired algorithms could dequantize these quantum algorithms in certain regimes. We improve on their result, giving an algorithm which is as general as the quantum algorithms, if the input is given classically (e.g., in a data-structure in RAM). Our goal is to solve the ε-feasibility problem; solving an SDP reduces by binary search to solving log(1/ε) instances of this feasibility problem. Problem 6.24 (SDP ε-feasibility). Given an ε > 0, m real numbers b 1 , . . . , b m ∈ R, and Hermitian n × n matrices SQ(A (1) ), . . . , SQ(A (m) ) such that -I ⪯ A (i) ⪯ I for all i ∈ [m], we define S ε as the set of all X satisfyingfoot_16
If S ε = ∅, output "infeasible". If S 0 ̸ = ∅, output an X ∈ S ε . (If neither condition holds, either output is acceptable.) Corollary 6.25. Let F ≥ max j∈[m] (∥A (j) ∥ F ), and supposefoot_17 F = Ω(1). Then we can solve Problem 6.24 with success probability ≥ 1 -δ in cost
providing sampling and query access to a solution.
. For the same feasibility problem, the previous quantum-inspired SDP solver [CLLW20] proved a complexity bound
, assuming that the constraint matrices have rank at most r. Since the rank constraint implies that
, under this assumption our algorithm has complexity
. So, our new algorithm both solves a more general problem and also greatly improves the runtime. The paper with the current best runtime for SDP solving does not discuss this precise model, but if we use the runtime they achieve in quantum state input model, making reasonable substitutions of γ → 1 ε and B → F 2 , the corresponding quantum runtime is ˜︁ O (︂
)︂ , up to polylog(n) factors. Like prior work on quantum algorithms for SDP-solving, we use the matrix multiplicative weights (MMW) framework [AK16,Kal07] to solve Problem 6.24. Corollary 6.25 immediately follows from running the algorithm this framework admits (Algorithm 1), where we solve an instance of the problem described in Lemma 6.26 with precision θ = ε/4 in each of the O (︁ log(n)/ε 2 )︁ iterations. MMW works as a zero-sum game with two players, where the first player wants to provide an X ∈ S ε , and the second player wants to find a violation for any proposed X, i.e., a j ∈ [m] such that Tr[A (j) X] > b j + ε. At the t th round of the game, if the second player points out a violation j t for the current solution X t , the first player proposes a new solution Our task is to execute Lines 3 and 4 of Algorithm 1, for an implicitly defined matrix with the form given in Line 6.
Algorithm 1: MMW based feasibility testing algorithm for SDPs 1 Set X 1 := In n , and the number of iterations T := 16 log n ε 2 ; 2 for t = 1, . . . , T do
else conclude that X t ∈ S ε 8 return X t 9 end 10 If no solution found, conclude that the SDP is infeasible and terminate the algorithm Lemma 6.26 ("Efficient" trace estimation). Consider the setting described in Corollary 6.25. Given θ ∈ (0, 1], t ≤ log(n)
where sq(A) = max j∈[m] sq(A (j) ), and s(A), q(A), n(A) are defined analogously.
To estimate Tr[A (i) H], we first notice that we have SQ ϕ (θ ∑︁ t i=1 A (j i ) ), since it is a linear combination of matrices that we have sampling and query access to (Lemma 3.9). Then, we can find approximations of the Gibbs state by applying eigenvalue transformation (Theorem 7.2) according to the exponential function to get exp[-θ ∑︁ t i=1 A (j i ) ] as an RUR decomposition. Then the estimation of Tr[A (i) H] can be performed by usual techniques (namely, Remark 4.13).
In order to understand how precisely we need to approximate the matrix in Line 6 we prove the following lemmas. Our first lemma will show that, to estimate Tr(A (i) H)/ Tr(H) to θ precision, it suffices to estimate both Tr(A (i) H) and Tr(H) to 1 3 θ Tr(H) precision.
Lemma 6.27. Suppose that θ ∈ [0, 1] and a, a ˜, Z, Z ˜are such that |a| ≤ Z, |a -a ˜| ≤ θ 3 Z, and
Next, we will prove that the approximations we will use to Tr(A (i) H) and Tr(H) suffice. We introduce some useful properties of matrix norms. For a matrix A ∈ C m×n and p ∈ [1, ∞], we denote by ∥A∥ p the Schatten p-norm, which is the ℓ p -norm of the singular values ( ∑︁ i σ p i (A)) 1/p . In particular, ∥A∥ F = ∥A∥ 2 and ∥A∥ Op = ∥A∥ ∞ . We recall some useful inequalities [Bha97, Section IV.2]. Hölder's inequality states that for all B ∈ C n×k and r, p, q ∈ (0, ∞] such that 1 p + 1 q = 1 r , we have ∥AB∥ r ≤ ∥A∥ p ∥B∥ q . The trace-norm inequality states that if n = m, then |Tr(A)| ≤ ∥A∥ 1 . Proof. First, recall that, by Lemma 2.5, there is unitary U such that ∥ ˜︁ U -U ∥ ≤ δ. Consequently, also using facts from Lemma 2.5, along with bounds on δ, Now we are ready to devise our upper bound on the trace estimation subroutine.
Proof of Lemma 6.26. By Lemma 6.27, it suffices to find estimates of Tr(e H ) and Tr(Ae H ) for all A = A (i) , to θ 3 Tr(e H ) additive precision. Recall from the statement that H := -θ ∑︁ t i=1 A (j i ) . By triangle inequality, ∥H∥ F ≤ F θ log(n). Because H is a linear combination of matrices, by Lemma 3.9, after paying log(n) θ 2 n(A) cost, we can obtain SQ ϕ (H) for ϕ ≤ F 2 log 2 (n)
θ 2 q(A) and s ϕ (H) = s(A).
√ n, so Tr(A (i) e H )/ Tr(e H ) ≤ e 2 F/ √ n ≤ θ, and outputting 0 as estimates is acceptable. We use Theorem 7.2 (with f (x) = x, so that L = 1, and choosing ε := Θ(θ)) to find a diagonal matrix
together with an approximate isome-
If every diagonal element is less than 3/4, then we conclude that ∥H∥ ≤ 1, and return 0. Otherwise we have ∥H∥ ≥ 1/2 and thus by Theorem 7.2 we have Tr(e H )). Note that by Lemma 6.29 and Eq. (15) we have
and since s = ˜︁ O (︁ F 6 θ -12 log 6 (n) log(1/δ) )︁ and q(H) ≤ log(n) θ 2 q(A), we also have
Therefore, Remark 4.13 tells us that given SQ(A), a ( θ 9 Tr(e H ))-approximation of Tr(A ˜︁ U † (e D -I) ˜︁ U ) can be computed with success probability at least 1 -δ 2m in time
)︃ .
Since we do this for all i ∈ [m], the overall complexity of obtaining the desired estimates Tr(A (i) e H ) with success probability at least 1 -δ 2 is m times
)︃ .
Discriminant analysis is used for dimensionality reduction and classification over large data sets.
Cong and Duan introduced a quantum algorithm to perform both with Fisher's linear discriminant analysis [CD16], a generalization of principal component analysis to data separated into classes. The problem is as follows: given classified data, we wish to project our data onto a subspace that best explains between-class variance, while minimizing within-class variance. Suppose there are M input data points {x i ∈ R N : 1 ≤ i ≤ M } each belonging to one of k classes. Let µ c denote the centroid (mean) of class c ∈ [k], and x ¯denote the centroid of all data points. Following the notation of [CD16], let
denote the between-class and within-class scatter matrices of the dataset respectively. The original goal is to solve the generalized eigenvalue problem S B v i = λ i S W v i and output the top eigenvalues and eigenvectors; for dimensionality reduction using linear discriminant analysis, we would project onto these top eigenvectors. If S W would be full-rank, this problem would be equivalent to finding the eigenvalues of S -1 W S B . However, this does not happen in general, and therefore various relaxations are considered in the literature [BHK97,Wel09]. For example, Welling [Wel09] considers the eigenvalue problem of
Cong and Duan further relax the problem, as they ignore small eigenvalues of S W and S B , and only compute approximate eigenvalues of Eq. (19) (after truncating eigenvalues), leading to inexact eigenvectors. We construct a classical analogue of their quantum algorithm. 26 Cong and Duan also describe a quantum algorithm for discriminant analysis classification; this algorithm does a matrix inversion procedure very similar to those described in Section 6.5 and Section 6.6, so for brevity we will skip dequantizing this algorithm.
To formally analyze this algorithm, we could, as in Section 6.4, assume the existence of an eigenvalue gap, so the eigenvectors are well-conditioned. However, let us instead use a different convention: if we can find diagonal D and an approximate isometry U such that S
B U ≈ U D, then we say we have found approximate eigenvalues and eigenvectors of S + W S B . Problem 6.30 (Linear discriminant analysis). Consider the functions
We prove this by using Theorem 5.1 to approximate sqrt(
Finding eigenvectors from an RUR decomposition follows from an observation (Lemma 4.14): for a matrix C W formed by sampling columns from R W (using SQ(W )), and [C W ] k the rank-k approximation to C W (which can be computed because C W has size independent of dimension), (([C W ] k ) + R W ) † has singular values either close to zero or close to one. This roughly formalizes the intuition of C W preserving the left singular vectors and singular values of R W . We can rewrite
which holds by choosing k sufficiently large and choosing C to be the same sketch used for U . Then, we can compute the eigendecomposition of the center
V is an approximate isometry, so we choose its columns to be our eigenvectors, and our eigenvalues are the diagonal entries of D. We show that this has the approximation properties analogous to the quantum algorithm.
. These approximations suffice for us:
each of which is bounded by ε∥B∥ 2 /σ 2 . Next, we approximate
and so
By Lemma 4.14, ([
† is an εσ/∥B∥-approximate projective isometry 28 onto the image of
(where we use that ε < σ/∥B∥). To turn this approximate projective isometry into an isometry, we compute the eigendecomposition
-and thus is an εσ/∥B∥-approximate isometry. So, our eigenvectors are U and our eigenvalues are D := Σ. This satisfies the desired bounds because
28 We get more than we need here: an ε-approximate projective isometry would suffice for the subsequent arguments.
The eigenvalues are correct because, by the approximate isometry condition, ∥U -U ˜∥ ≲ ε σ ∥B∥ for U ãn isometry, and so we can use Lemma 2.5 to conclude
U ˜DU ˜ † is an eigendecomposition. Furthermore, this is an approximation of a Hermitian PSD matrices, where singular value error bounds align with eigenvalue error bounds. So, Weyl's inequality (Lemma 4.11) implies the desired bound |D ii -λ i | ≲ ε ∥B∥ 2 σ 2 for λ i the true eigenvalues. We have SQ ϕ (U (•, i)) by Lemmas 3.5 and 3.6, since U (
This gives the stated runtime.
In this section, we present more general versions of our algorithm for even SVT to get results for generic SVT (Theorem 7.1) and eigenvalue transformation (Theorem 7.2). In applications we mainly use even SVT to allow for more fine-tuned control over runtime, but we do use eigenvalue transformation in Section 6.8. For generic SVT: consider a matrix A ∈ C m×n and a function f : [0, ∞) → C satisfying f (0) = 0 (so the singular value transformation f (SV) (A) is well-defined as in Definition 2.1). Given SQ(A) and SQ(A † ), we give an algorithm to output a CUR decomposition approximating f (SV) (A).
Theorem 7.1 (Generic singular value transformation). Let A ∈ C m×n be given with both SQ ϕ (A) and SQ ϕ (A † ) and let f : [0, ∞) → C be a function such that f (0) = 0, g(x) := f ( √ x)/ √ x is L-Lipschitz, and g ¯(x) := (g(x) -g(0))/x is L ¯-Lipschitz. Then, for 0 < ε ≤ min(L∥A∥ 3 , L ¯∥A∥ 5 ), we can output sketches R := SA ∈ C r×n and C := AT ∈ C m×c , along with M ∈ C r×c such that
. Finding S, M , and T takes time
If we only wish to assume SQ ϕ (A), we can do so by using Lemma 4.4 instead of Lemma 4.6 in our proof, paying an additional factor of 1 δ . Note that if sq ϕ (A), sq ϕ (A † ) = O(1), then this runtime is dominated by the last term. Moreover, if A is strictly low-rank, with minimum singular value σ, or essentially equivalently, if f (x) = 0 for x ≤ σ and so g(x) = 0 for x ≤ σ 2 , then L ≤ ℓ/σ 2 and L ¯= 2ℓ/σ 4 for ℓ the Lipschitz constant of f . In this case the complexity is
)︂(︂ ℓ∥A∥ ε
Importantly, when ε = O(ℓ∥A∥) (that is, if we want relative error), this runtime is independent of dimension. If one desires greater generality, where we only need to depend on the Lipschitz constant of f , we can use a simple trick: as we aim for a spectral norm bound, we can essentially treat A as if it had strictly low rank. Consider the variant of f , f ≥σ , which is zero below σ/2, f above σ, and is a linear interpolation in between.
Further, the Lipschitz constant of f ≥ε/ℓ is at most 2ℓ: the slope of the linear interpolation is 2f (σ)/σ ≤ 2ℓσ/σ. So, we can run our algorithm for arbitrary ℓ-Lipschitz f in the time given by Eq. ( 20), with σ = ε/ℓ.
Our proof strategy is to apply our main result Theorem 5.1 to g(A † A), for g(x) := f ( √ x)/ √ x, and subsequently approximate matrix products with Lemma 4.6 to get an approximation of the form A ′ R ′ † U R + g(0)A:
Here, A ′ R ′ † U R is a CUR decomposition as desired, since A ′ is a normalized subset of columns of A. One could further approximate g(0)A by a CUR decomposition if necessary (e.g. by adapting the eigenvalue transformation result below).
We do not use this theorem in our applications. Sometimes we implicitly use a similar strategy (e.g. in Section 6.5), but because we apply our matrix to a vector (f (A † )b) we can use Lemma 4.12 instead of Lemma 4.6 when approximating. This allows for the algorithm to work with only SQ ϕ (A) and still achieve a poly-logarithmic dependence on 1 δ .
Proof. If we want to compute f ˆ(SV) (A), we can work with f (x) := f ˆ(x) -g(0)x, so that g(0) = 0 without loss of generality. Notice that f ˆ(SV) (A) = f (SV) (A)+g(0)A, so if we get a CUR decomposition for f (SV) (A) we can add g(0)A after to get the decomposition in the theorem statement. Consider the SVT g(x)
Second, use Lemma 4.6 to get a sketch T ′ † ∈ C c ′ ×n such that, with probability ≥ 1 -δ/4,
The choices of parameters necessary are as follows (using that ∥SA∥ F = O(∥A∥ F ) by Eq. ( 5) and we have a 2ϕ-oversampled distribution for (SA) † by Lemma 4.3):
This implies the desired bound through the following sequence of approximations:
This gives us a CUR decomposition of f (SV) (A). These two approximations only incur O(ε) error in spectral norm; for the first, notice that
(by ( 21))
For the second approximation observe that |g(x)| ≤ L|x| (and, by corollary, g ¯(x) ≤ L) due to g being L-Lipschitz and g(0) = 0, therefore
The time complexity of this procedure is
which comes from producing sketches, querying all the relevant entries of SAT and SAT ′ , the singular value transformation of SAT , and the matrix multiplication in M . We get r factors in the latter two terms because we can separate g ¯((SAT )(SAT )
As for eigenvalue transformation, consider a function f : R → C and a Hermitian matrix
, so we can use Theorem 5.1 to compute the eigenvalue transform f (A). For non-even f , we cannot use this result, and present the following algorithm to compute it. Theorem 7.2 (Eigenvalue transformation). Suppose we are given a Hermitian SQ ϕ (A) ∈ C n×n with eigenvalues λ
L , and some ε ∈ (0, L∥A∥/2]. Then we can output matrices S ∈ C r×n , N ∈ C s×r , and
Moreover, this decomposition satisfies the following two properties. First, N SA is an approximate isometry:
Under the reasonable assumptionsfoot_21 that sq(A) and ϕ are small (O(1), say) and ε ≤ L∥A∥ ∥A∥ ∥A∥ F , the complexity of this theorem is ˜︁ O (︁ L 22 ε -22 ∥A∥ 16 ∥A∥ 6 F log 3 1 δ )︁ . We now outline our proof. Our strategy is similar to the one used for quantum-inspired semidefinite programming [CLLW20]: first we find the eigenvectors and eigenvalues of A and then apply f to the eigenvalues. Let π(x) be a (smoothened) step function designed so it can zeroes out small eigenvalues of A (in particular, eigenvalues smaller than ε/ √ 2L). Then
Here, C σ is the low-rank approximation of C formed by transforming C according to the "filter" function on x that is 0 for x < σ and x otherwise. U ˆ:= C + σ R ∈ C c×n is an approximate isometry by Lemma 4.14. We are nearly done now: since the rest of the matrix expression,
consists of submatrices of A of size independent of n, we can directly compute its unitary eigendecomposition U DU † . This gives the approximate decomposition A ≈ (U ˆU )D(U ˆU ) † , with U ˆU and D acting as approximate eigenvectors and eigenvalues of A, respectively. An application of Lemma 5.4 shows that f (A) ≈ (U ˆU )f (D)(U ˆU ) † in the desired sense. Therefore, our output approximation of f (A) comes in the form of an RUR decomposition that can be rewritten in the form of an approximate eigendecomposition. The only major difference between this proof sketch and the proof below is that we perform our manipulations on the SVD of C σ , to save on computation time: note that the SVD can be made small in dimension, using that the rank of C σ is bounded by ∥C∥ 2 F /σ 2 .
Proof. Throughout this proof ε is not dimensionless; if choices of parameters are confusing, try replacing ε with ε∥A∥. We will take f (0) = 0 without loss of generality. First, consider the "smooth projection" singular value transformation
Since π is a projector onto the large eigenvectors of A, we can add these projectors to our expression without incurring too much spectral norm error.
The necessary sizes for these bounds to hold are as follows (Lipschitz constants for π are 2L 2 /ε 2 and 4L 4 /ε 4 , ∥SA∥ F = O(∥A∥ F ) by Eq. ( 5), and we have a 2ϕ-oversampled distribution for (SA) † by Lemma 4.3):
This approximation does not incur too much error:
Third, apply Remark 4.13(b) r 2 times to approximate each entry of RAR † : pull t samples from SQ ϕ (A) for t
)︂ such that, given some Q(x), Q(y), with probability ≥ 1 -δ r 2 , one can output an estimate of x † Ay up to ε 3 ∥x∥∥y∥ L 3 ∥A∥ 2 F additive error with no additional queries to SQ ϕ (A). Then, by union bound, with probability ≥ 1 -δ, using the same t samples from A each time, one can output an estimate of R(i, •)AR(j, •) † up to ε 3 ∥R(i,•)∥∥R(j,•)∥ L 3 ∥A∥ 2 F error for all i, j ∈ [r] such that i ≤ j. Let M be the matrix of these estimates. Then, using that ∥R∥ F = O(∥A∥ F ) from Eq. (5),
From Eqs. ( 6) and (7),
So far, we have shown that we can find an RUR approximation to A, with
However, if we wish to apply an eigenvalue transformation to A, we need to access the eigenvalues of A as well. To do this, we will express this decomposition as an approximate unitary eigendecomposition.
Using that π ¯zeroes out singular values that are smaller than ε 2 2L 2 , we can write our expression as U ˆD ˆU ˆ †, for U ˆ∈ C n×s and D ˆ∈ C s×s :
Here, we are using an SVD of C truncated to ignore singular values smaller than ε √ 2L , where
approximate isometry; we rescale ε until this is at most 1 2 . The rest of this error analysis will show that, since U ˆis an approximate isometry, f (A) ≈ U ˆf (D ˆ)U ˆ † in the senses required for the theorem statement. Though D ˆis not diagonal, since it is s×s, we can compute its unitary eigendecomposition U (D ˆ)D (D ˆ)(U (D ˆ)) † ; so, we can take
) † to get the decomposition in the theorem statement. (Including the isometry (U (D ˆ)) † in our expression for U ˆdoes not change the value of α).
First, consider the eigenvalues of D ˆ. Note that ∥U ˆ †U ˆ-I∥ ≤ α since U ˆis an α-approximate isometry, and by Lemma 2.5, there exists an isometry U such that ∥U -U ˆ∥ ≤ α. We first observe that, using Lemma 2.5 and our bound on α,
Consequently, by Weyl's inequality (Lemma 4.11), the eigenvalues of
and by assumption, f is L-Lipschitz on the spectrums of A and U D ˆU . From this, we can conclude that we can compute estimates for the eigenvalues of f (A), since the eigenvalues of U D ˆU † are the eigenvalues of D ˆ(padded with zero eigenvalues) which we know from our eigendecomposition of D ˆFurther, our estimates f (λ ˆi) satisfy the desired bound |f (λ ˆi) -f (λ i )| ≤ ε. Finally, since f is Lipschitz on our spectrums of concern, the desired error bound for our approximation holds by the following computation (which uses Lemma 2.5 extensively): The complexity analysis takes some care: we want to compute our matrix expressions in the correct order. First, we will sample to get S and T , and then compute the truncated singular value decomposition of C := SAT , which we have denoted
∈ C c×s . Then, we will perform the inner product estimation protocol r 2 times to get our estimate M ∈ C r×r , and compute the eigendecomposition of
) 2 )(U
) 2 )D ) + (U
By evaluating the expression for D ˆfrom left-to-right, we only need to perform matrix multiplications that (naively) take s 3 or sr 2 time. The only cost of c we incur is in computing the SVD of C. The runtime is O (︁ (r + c + t)(s ϕ (A) + q ϕ (A)) + (rc + r 2 t) q(A) + s 3 + r 2 s + r 2 c + n ϕ (A) )︁ .
and has the same columnspace as X ˆ, and
, and consequently Let S ∈ R r×n be sampled according to p. Then X † S † SY is an unbiased estimator for X † Y and
Proof. Using that the rows of S are selected independently, we can conclude the following:
To prove concentration, we use McDiarmid's "independent bounded difference inequality" [McD89]. )︂ .
To use Lemma B.1, we think about this expression as a function of the indices that are randomly chosen from p. That is, let f be the function [n] r → R defined to be
Then, by Jensen's inequality, we have
Now suppose that the index sequences i ⃗ and i ⃗ ′ only differ at the s-th position. Then by the triangle inequality, So, with probability ≥ 1 -δ, So, we just need to boost the quality of this random variable. Consider taking Z ¯to be the mean of x := 8ϕ∥u∥ 2 ∥v∥ 2 1 ε 2 independent copies of Z. Then, by Chebyshev's inequality (stated here for complex-valued random variables),
Next, we take the (component-wise) median of y := 8 log 1 δ independent copies of Z ¯, which we call Z ˜, to decrease failure probability. Consider the median of the real parts of Z ¯. The key observation is that if ℜ(Z ˜-E[Z]) ≥ ε/ √ 2, then at least half of the Z ¯'s satisfy ℜ(Z ¯-
2) be the characteristic function for this event for a particular mean. The above argument implies that Pr[E i ] ≤ 1 4 . So, by Hoeffding's inequality, Pr
]︄ ≤ exp(-q/8) ≤ δ 2 .
With this combined with our key observation, we can conclude that Pr[ℜ(Z ˜-⟨u, v⟩) ≥ ε/ √ 2] ≤ δ/2. From a union bound together with the analogous argument for the imaginary component, we have Pr[|Z ˜-⟨u, v⟩| ≥ ε] ≤ δ as desired. The time complexity is the number of samples multiplied by the time to create one instance of the random variable Z, which is O(sq(u) + q(v)). Lemma 6.3. Consider p(x) a degree-d polynomial of parity-d such that |p(x)| ≤ 1 for x ∈ [-1, 1]. Recall that, for a function f : C → C, we define f ¯(x) := (f (x) -f (0))/x (and f ¯(0) = f ′ (0) when f is differentiable at zero).
• If p is even, then max
where ≲ hides a constant depending on k. Note that by replacing x in the above equation with 2y -1, we get that max y∈[0,1] |p (k) (y)| ≲ d 2k max y∈[0,1] |p(y)| (paying an additional 2 k constant factor). We make a couple observations about r ¯(x) using Taylor expansions, where r(x) is any degree-d polynomial. . Note that when γ -1/2 is chosen to be sufficiently large (e.g. O(∥X∥ F )) and ∥α∥ = Ω(γ∥y∥), this runtime is dimension-independent.
Proof. Denote σ 2 := γ -1 , and redefine X ← X T (so we have SQ(X) instead of SQ(X T )). By the block matrix inversion formulafoot_23 we know that
So, we have reduced the problem of inverting the modified matrix to just inverting M -1 where M = X T X + σ -2 I. M is invertible because M ⪰ σ 2 I. Note that M -1 = f (X T X), where
So, by Theorem 5.1, we can find R † f ¯(CC † )R such that ∥R † f ¯(CC † )R + 1 σ 2 I -f (X T X)∥ ≤ εσ -2 , where (because L = σ -4 , L ¯= σ -6 )
So, the runtime for estimating this is ˜︁ O (︂
. We further approximate using Lemma 4.6:
we find r 1 ≈ R † 1 ⃗ , r y ≈ R † y ⃗, and γ ≈ 1 ⃗ † y in O(r K ε 2 log 1 δ ) time (for the first two) and O( 1 ε 2 log 1 δ ) time (for the last one) such that the following bounds hold:
Via Lemma 5.2, we observe the following additional bounds:
Now, we compute what the subsequent errors are for replacing M -1 with N := R † ZR + 1 σ 2 I, With these bounds, we can conclude that (continuing to use Eqs. ( 25) and ( 26))
by Eq. ( 27)
So, by rescaling ε down by ∥X∥ σ , it suffices to sample from α ˆ:= R † Z (︂ r y -r † 1 Zr y + σ -2 γ r †
1 Zr 1 + σ -2 m r 1 )︂ -σ -2 (︂ y -r † 1 Zr y + σ -2 γ r †
1 Zr 1 + σ -2 m 1 ⃗ )︂ .
To gain sampling and query access to the output, we consider this as a matrix-vector product, where the matrix is (R † | y | 1 ⃗ ) and the vector is the corresponding coefficients in the linear combination. Then, by Lemmas 3.5 and 3.6, we can get SQ ϕ (α ˆ) for ϕ = (r + 2)
Now we consider z(t) := g I (t) = Tr (︁ e H+t(H ˜-H) )︁ . From Eq. ( 29) we have z ′ (t) ≤ ∥H ˜-H∥z(t). Using Grönwall's differential inequality, we can conclude that z(t) ≤ z(0)e t∥H ˜-H∥ for every t ∈ [0, ∞).
Finally, we use the fact that there exists a matrix A of operator norm at most 1 such that ⃦ ⃦ e H ˜e H ⃦ ⃦ 1 = Tr(A(e H ˜e H )) (take, e.g., sgn(e H ˜e H )). We finish the proof by observing that for such an A,
This assumption is important. When input data is quantum (say, it is coming directly from a quantum system), a classical computer has little hope of performing linear algebra on it efficiently, see for example[ACQ22,HKP21].
The QML algorithms we discuss generally only incur polylog( 1 ε ) terms, but need to eventually pay poly(1/ε) to extract information from output quantum states. So, we believe this exponential speedup is artificial. See the open questions section for more discussion of this error parameter.
This is the relevant variant of the notion of a CUR decomposition from the randomized numerical linear algebra and theoretical computer science communities[DMM08].
For semidefinite program solving, ∥A (•) ∥F ≤ √ k, which makes the runtimes comparable.
The same holds for quantum algorithms using the QRAM data structure input model: the data structure itself can be built during an O(nnz(A))-time (classical) preprocessing phase.
The QRAM data structure used here has alternatives which in some sense vary the "norm" in which one stores the input matrix [KP20, Theorem IV.4], [CGJ19, Lemma 25]. We do not expect that these alternatives would give sampling and query access to the matrix they store or could be otherwise dequantized, since these datastructures generalize and strengthen the sparse-access input model, which is known to be BQP-complete[HHL09].
This is the notion of approximate orthonormality as given by the first arXiv version of[Tan19].
Oversampling turns out to be the "natural" form of approximation in this setting; other forms of error do not propagate through quantum-inspired algorithms well.
See the proof of Lemma 4.5 in Appendix B for this kind of computation done with more detail.
[Bha97, Corollary III.2.2] actually proves the Hermitian version. The result about singular values is an easy consequence, see for example the blog of Terence Tao [Tao10, Exercise 22(iv)].
By a dimension-independent or dimensionless quantity, we mean a quantity that is both independent of the size of the input matrix and is scale-invariant, i.e., does not change under scaling A ← αA.
This criterion is fairly reasonable. For example, the polynomials used in QSVT satisfy it.
Formally, to get a true approximation AR ≈ A ′ R, we need to union bound the failure probability for each row, paying a log m factor in runtime. However, we will ignore this consideration: our goal is to sample from one row, so we only need to succeed in our particular row.
The quantum analogue to SQ(X) is efficient state preparation of X, a purification of ρ.
The algorithm as written in[RML14] assumes that ∥F ∥ ≤ 1; we confirmed with an author that this is a typo.
The correctness of this dequantization is unclear, since the approximations performed in this section incur significant errors.
For simplicity, we assume here that X is normalized to have trace 1. This can be relaxed; for an example, see[vAGGdW20].
Because of the normalization assumption that ∥A (•) ∥ ≤ 1, F is effectively a dimensionless "stable rank"-type constant, normalized by maxi ∥A (i) ∥.
In case applying Theorem 7.2 would result in s > n, we instead directly diagonalize H ensuring s ≤ n.
Analyzing whether or not the particular relaxation used in this and other quantum machine learning papers provides a meaningful output is unfortunately beyond the scope of our paper.
This is the runtime of Step 2 of Algorithm 1. The normalization factor of max(∥B∥F, ∥S∥F) is implicit there,κ ef f corresponds to max(∥B∥ F ,∥S∥ F ) σ 2, and the error bound the algorithm achieves is the one we describe here, since the authors must implicitly rescale the inverse and square root function by a cumulative factor of ∥B∥ 2 /σ 2 to apply their Theorem 1.
The correct way to think about ε is as some constant fraction of L∥A∥. If ε > L∥A∥ then f (0)I is a satisfactory approximation. The bound we give says that we want an at least ∥A∥F/∥A∥ improvement over trivial, which is modest in the close-to-low-rank regime that we care about. Similar assumptions appear in Section 6.
The constraint on the size of ε from Theorem 5.1 here is ε(L∥A∥) -1 ≲ min(L 2 ∥A∥ 2 /ε 2 , L 4 ∥A∥ 4 /ε 4 ), which is true since ε ≤ L∥A∥/2.
In a more general setting, we would use the Sherman-Morrison inversion formula, or the analogous formula for functions of matrices subject to rank-one perturbations.
Note that in case A = I, by the cyclicity of trace, this equation implies that d dt Tr(e H(t) ) = Tr(e H(t) d dt H(t)).