Quantum-inspired algorithms in practice

TL;DR

This paper evaluates the practical performance of quantum-inspired algorithms for recommendation systems and linear equations, revealing they work well under specific conditions but face challenges with high-rank, sparse data.

Contribution

The study provides a theoretical and empirical assessment of quantum-inspired algorithms, identifying their strengths and limitations in real-world scenarios.

Findings

Performance exceeds theoretical bounds in some cases.

Degrades with increasing matrix rank and condition number.

Effective for low-rank, large-dimension matrices in practice.

Abstract

We study the practical performance of quantum-inspired algorithms for recommendation systems and linear systems of equations. These algorithms were shown to have an exponential asymptotic speedup compared to previously known classical methods for problems involving low-rank matrices, but with complexity bounds that exhibit a hefty polynomial overhead compared to quantum algorithms. This raised the question of whether these methods were actually useful in practice. We conduct a theoretical analysis aimed at identifying their computational bottlenecks, then implement and benchmark the algorithms on a variety of problems, including applications to portfolio optimization and movie recommendations. On the one hand, our analysis reveals that the performance of these algorithms is better than the theoretical complexity bounds would suggest. On the other hand, their performance as seen in our implementation degrades noticeably as the rank and condition number of the input matrix are increased. Overall, our results indicate that quantum-inspired algorithms can perform well in practice provided that stringent conditions are met: low rank, low condition number, and very large dimension of the input matrix. By contrast, practical datasets are often sparse and high-rank, precisely the type that can be handled by quantum algorithms.

Tables

Table 1. Table 1: Errors in estimating singular values, singular vectors, coefficients, and solution vectors for input matrices of dimension n , m = 2 50 𝑛 𝑚 superscript 2 50 n,m=2^{50} . The coefficients were estimated using N = 10 4 𝑁 superscript 10 4 N=10^{4} samples and the errors in the vectors were calculated by computing their first L = 100 𝐿 100 L=100 entries. The error bars correspond to one standard deviation taken over ten repetitions of the algorithm. Runtimes correspond to a Python implementation of the algorithms running on 20 Intel Xeon CPUs operating at 2.4GHz with access to 252GB of shared memory.

	Errors
Case study	${\eta }_{\sigma }$	${\eta }_{𝒗}$	${\eta }_{\lambda }$	${\eta }_{𝒙}$	Runtime (hours)
$k=\kappa ={\kappa }_{\beta }=3$	$0.011\pm 0.012$	$0.124\pm 0.071$	$0.285\pm 0.276$	$0.414\pm 0.154$	4.3
$k=\kappa ={\kappa }_{\beta }=5$	$0.129\pm 0.016$	$0.212\pm 0.070$	$0.530\pm 0.616$	$1.235\pm 0.615$	8.6
$k=\kappa ={\kappa }_{\beta }=10$	$0.626\pm 0.018$	$1.619\pm 0.264$	$1.193\pm 1.720$	$4.138\pm 2.321$	29.0

Table 2. Table 2: Relative errors with corresponding standard deviations computed over 10 10 10 repetitions of the FKV algorithm of the quantum-inspired algorithm for linear systems A ​ 𝒙 = 𝒃 𝐴 𝒙 𝒃 A\bm{x}=\bm{b} . A 𝐴 A is a random matrix with dimension m = 40000 𝑚 40000 m=40000 , n = 20000 𝑛 20000 n=20000 , rank k = 5 𝑘 5 k=5 , and condition number κ = 5 𝜅 5 \kappa=5 . Approximated singular values and vectors were obtained by decomposing matrix C 𝐶 C with dimension r = 4250 𝑟 4250 r=4250 and c = 4250 𝑐 4250 c=4250 . The coefficients λ ~ ~ 𝜆 \widetilde{\lambda} have been estimated by performing N = 10 4 𝑁 superscript 10 4 N=10^{4} samples. Here, η σ subscript 𝜂 𝜎 \eta_{\sigma} is the error in approximating singular values σ ~ ~ 𝜎 \tilde{\sigma} , while η A subscript 𝜂 𝐴 \eta_{A} and η A + subscript 𝜂 superscript 𝐴 \eta_{A^{+}} are the errors of the reconstructed matrices A ~ ~ 𝐴 \widetilde{A} and A ~ + superscript ~ 𝐴 \widetilde{A}^{+} , respectively. η λ subscript 𝜂 𝜆 \eta_{\lambda} is the error in estimating coefficients, and η 𝒙 subscript 𝜂 𝒙 \eta_{\bm{x}} is the the median of the relative errors | x ~ i − x i | / | x i | subscript ~ 𝑥 𝑖 subscript 𝑥 𝑖 subscript 𝑥 𝑖 |\tilde{x}_{i}-x_{i}|/|x_{i}| with respect to all n = 20000 𝑛 20000 n=20000 entries of the approximate vector 𝒙 ~ ~ 𝒙 \tilde{\bm{x}} .

	Parameters			Error
Case study	$r$	$c$	$N$	${\eta }_{\sigma }$	${\eta }_A$	${\eta }_{A^{+}}$	${\eta }_{\lambda }$	${\eta }_{𝒙}$
Random matrix	$4250$	$4250$	${10}^4$	$0.010\pm 0.005$	$0.028\pm 0.004$	$0.101\pm 0.027$	$0.387\pm 0.191$	$0.087\pm 0.053$

Table 3. Table 3: Running times (in seconds) for the quantum-inspired algorithm applied to solve the system of linear equations A ​ 𝒙 = 𝒃 𝐴 𝒙 𝒃 A\bm{x}=\bm{b} with the same parameters as reported in Table 2 . The parameter t LS subscript 𝑡 LS t_{\mathrm{LS}} is the time required to construct length-squared (LS) probability distributions over rows and columns. t SVD C superscript subscript 𝑡 SVD 𝐶 t_{\mathrm{SVD}}^{C} accounts for the time spent sampling rows and columns from matrix A 𝐴 A and building and decomposing the sampled matrix C 𝐶 C . t λ subscript 𝑡 𝜆 t_{\lambda} is the running time of the estimation of the coefficients λ 𝜆 \lambda . Analogous quantities are defined for the direct calculation method. t 𝒙 subscript 𝑡 𝒙 t_{\bm{x}} is the elapsed time to sample 500 500 500 entries of the approximated solution vector in the quantum-inspired algorithm and the time spent to compute the exact solution vector in the direct calculation method. Runtimes correspond to a Python implementation of the algorithms running on two Intel Xeon CPUs operating at 2.4GHz with access to 252GB of shared memory.

	Quantum-inspired algorithm					Direct calculation
Case study	$t_{\mathrm{LS}}$	$t_{\mathrm{SVD}}^C$	$t_{\lambda }$	$t_{𝒙}$	$t_{\mathrm{total}}$	$t_{\mathrm{SVD}}^A$	$t_{\lambda }$	$t_{𝒙}$	$t_{\mathrm{total}}$
Random matrix	$1488.8$	$83.9$	$554.7$	$343$	$2470.4$	$5191.1$	$1.4$	$0.0003$	$5192.5$

Table 4. Table 4: Relative errors with corresponding standard deviations computed over 10 repetitions of the quantum-inspired algorithm applied to a portfolio optimization problem. Approximated singular values and vectors were obtained by decomposing matrix C 𝐶 C with dimension r = 340 𝑟 340 r=340 and c = 340 𝑐 340 c=340 . Coefficients λ 𝜆 \lambda have been estimated by performing N = 10 4 𝑁 superscript 10 4 N=10^{4} samples. Here, η σ subscript 𝜂 𝜎 \eta_{\sigma} is the error in approximating singular values σ ~ ~ 𝜎 \tilde{\sigma} , while η A subscript 𝜂 𝐴 \eta_{A} and η A + subscript 𝜂 superscript 𝐴 \eta_{A^{+}} are the errors of the reconstructed matrices A ~ ~ 𝐴 \widetilde{A} and A ~ + superscript ~ 𝐴 \widetilde{A}^{+} , respectively. η λ subscript 𝜂 𝜆 \eta_{\lambda} is the error in estimating coefficients, and η 𝒙 subscript 𝜂 𝒙 \eta_{\bm{x}} is the the median of the relative errors | x ~ i − x i | / | x i | subscript ~ 𝑥 𝑖 subscript 𝑥 𝑖 subscript 𝑥 𝑖 |\tilde{x}_{i}-x_{i}|/|x_{i}| with respect to all n = 475 𝑛 475 n=475 entries of the approximated vector 𝒙 ~ ~ 𝒙 \tilde{\bm{x}} .

	Parameters			Error
Case study	$r$	$c$	$N$	${\eta }_{\sigma }$	${\eta }_A$	${\eta }_A^{+}$	${\eta }_{\lambda }$	${\eta }_{𝒙}$
S&P 500 stock index	$340$	$340$	${10}^4$	$0.08\pm 0.02$	$0.16\pm 0.03$	$1.13\pm 0.12$	$1.58\pm 0.79$	$0.74\pm 0.19$

Table 5. Table 5: Running times (in seconds) for the quantum-inspired algorithm applied to the portfolio optimization problem. The algorithm was run using the parameters reported in Table 4 . The parameter t LS subscript 𝑡 LS t_{\mathrm{LS}} is the time required to construct length-squared (LS) probability distributions over rows and columns. t SVD C superscript subscript 𝑡 SVD 𝐶 t_{\mathrm{SVD}}^{C} accounts for the time spent sampling rows and columns from matrix A 𝐴 A and building and decomposing the sampled matrix C 𝐶 C . t λ subscript 𝑡 𝜆 t_{\lambda} is the running time of the estimation of the coefficients λ 𝜆 \lambda . Analogous quantities are defined for the direct calculation method. t 𝒙 subscript 𝑡 𝒙 t_{\bm{x}} is the elapsed time to sample 50 50 50 entries of the approximated solution vector in the quantum-inspired algorithm, and the time spent to compute the exact solution vector in the direct calculation method. Runtimes correspond to a python implementation of the algorithms running on two Intel Xeon CPUs operating at 2.4GHz with access to 252GB of shared memory.

	Quantum-inspired algorithm					Direct calculation
Case study	$t_{\mathrm{LS}}$	$t_{\mathrm{SVD}}^C$	$t_{\lambda }$	$t_{𝒙}$	$t_{\mathrm{total}}$	$t_{\mathrm{SVD}}^A$	$t_{\lambda }$	$t_{𝒙}$	$t_{\mathrm{total}}$
S&P 500 stock index	$0.46$	$0.36$	$88.7$	$24.63$	$114.15$	$0.15$	$0.0005$	$0.0006$	$0.15$

Table 6. Table 6: Relative errors with corresponding standard deviations computed over 10 repetitions of the FKV algorithm in the quantum-inspired algorithm for recommendation systems. Here, matrix A 𝐴 A is a preference matrix with dimension m = 611 𝑚 611 m=611 and n = 9274 𝑛 9274 n=9274 . Approximated singular values and vectors were obtained by decomposing matrix C 𝐶 C with dimension r = 450 𝑟 450 r=450 and c = 4500 𝑐 4500 c=4500 . The coefficients λ 𝜆 \lambda have been estimated by performing N = 10 4 𝑁 superscript 10 4 N=10^{4} samples. Here, η σ subscript 𝜂 𝜎 \eta_{\sigma} is the error in approximating singular values σ ~ ~ 𝜎 \tilde{\sigma} , while η A subscript 𝜂 𝐴 \eta_{A} and η A + subscript 𝜂 superscript 𝐴 \eta_{A^{+}} are the errors of the reconstructed matrices A ~ ~ 𝐴 \widetilde{A} and A ~ + superscript ~ 𝐴 \widetilde{A}^{+} , respectively. η λ subscript 𝜂 𝜆 \eta_{\lambda} is the error in estimating coefficients, and η 𝒙 subscript 𝜂 𝒙 \eta_{\bm{x}} is the the median of the relative errors | x ~ i − x i | / | x i | subscript ~ 𝑥 𝑖 subscript 𝑥 𝑖 subscript 𝑥 𝑖 |\tilde{x}_{i}-x_{i}|/|x_{i}| with respect to all n = 9274 𝑛 9274 n=9274 entries of the approximated solution vector 𝒙 𝒙 \bm{x} .

	Parameters			Error
Case study	$r$	$c$	$N$	${\eta }_{\sigma }$	${\eta }_A$	${\eta }_{A^{+}}$	${\eta }_{\lambda }$	${\eta }_{𝒙}$
MovieLens 100K	$450$	$4500$	${10}^4$	$0.06\pm 0.01$	$0.32\pm 0.02$	$0.66\pm 0.05$	$0.58\pm 0.18$	$0.71\pm 0.13$

Table 7. Table 7: Running times (in seconds) for the quantum-inspired algorithm for recommendation systems. The algorithm was run using the parameters reported in Table 6 . The parameter t LS subscript 𝑡 LS t_{\mathrm{LS}} is the time required to construct length-squared (LS) probability distributions over rows and columns. t SVD C superscript subscript 𝑡 SVD 𝐶 t_{\mathrm{SVD}}^{C} accounts for the time spent sampling rows and columns from matrix A 𝐴 A and building and decomposing the sampled matrix C 𝐶 C . t λ subscript 𝑡 𝜆 t_{\lambda} is the running time of the estimation of the coefficients λ 𝜆 \lambda . Analogous quantities are defined for the direct calculation method. t 𝒙 subscript 𝑡 𝒙 t_{\bm{x}} is the elapsed time to sample 500 500 500 entries of the approximated solution vector in the quantum-inspired algorithm and the time spent to compute the exact solution vector in the direct calculation method. Runtimes correspond to a python implementation of the algorithms running on two Intel Xeon CPUs operating at 2.4GHz with access to 252GB of shared memory.

	Quantum-inspired algorithm					Direct calculation
Case study	$t_{\mathrm{LS}}$	$t_{\mathrm{SVD}}^C$	$t_{\lambda }$	$t_{𝒙}$	$t_{\mathrm{total}}$	$t_{\mathrm{SVD}}^A$	$t_{\lambda }$	$t_{𝒙}$	$t_{\mathrm{total}}$
MovieLens 100K	$11.05$	$6.22$	$124.2$	$14.2$	$155.7$	$2.00$	$0.0003$	$0.001$	$2.00$

Equations

$$A=\sum _{\ell =1}^k{\sigma }_{\ell }{\mathbf{u}}^{(\ell )}{{\mathbf{v}}^{(\ell )}}^T,$$

$$\mathbf{x}=\sum _{\ell =1}^k{\lambda }_{\ell }{\mathbf{v}}^{(\ell )},$$

$${\lambda }_{\ell }=\frac{1}{{\sigma }_{\ell }^2}⟨{\mathbf{v}}^{(\ell )},A^T\mathbf{b}⟩.$$

$${\lambda }_{\ell }=⟨A_i^T,{\mathbf{v}}^{(\ell )}⟩,$$

$${\tilde{\mathbf{v}}}^{(\ell )}=\frac{1}{{\tilde{\sigma }}_{\ell }}{R}^T{\mathbf{\omega }}^{(\ell )},$$

$${\tilde{\mathbf{u}}}^{(\ell )}=A\left(\frac{1}{{\tilde{\sigma }}_{\ell }^2}{R}^T{\mathbf{\omega }}^{(\ell )}\right).$$

$$\lambda =⟨\mathbf{y},\mathbf{z}⟩=\sum _{i=1}^n{y}_i{z}_i,$$

$$\mathbb{E}(\chi )=\sum _{i=1}^n\frac{y_i{z}_i}{p_y(i)}{p}_y(i)=⟨\mathbf{y},\mathbf{z}⟩=\lambda .$$

$$\mathbb{E}({\chi }^2)=\sum _{i=1}^n{\left(\frac{y_i{z}_i}{p_y(i)}\right)}^2{p}_y(i)=\Vert \mathbf{y}{\Vert }^2\Vert \mathbf{z}{\Vert }^2,$$

$$N=$$

$$=$$

$${\lambda }_{\ell }=\frac{1}{{\sigma }_{\ell }^2}⟨{\mathbf{v}}^{(\ell )},A^T\mathbf{b}⟩=\frac{1}{{\sigma }_{\ell }}⟨{\mathbf{u}}^{(\ell )},\mathbf{b}⟩.$$

$$N$$

$$⟨A_i^T,{\mathbf{v}}^{(\ell )}⟩=\sum _{{\ell }^{\prime }=1}^k{\sigma }_{{\ell }^{\prime }}{u}_i^{({\ell }^{\prime })}⟨{\mathbf{v}}^{({\ell }^{\prime })},{\mathbf{v}}^{({\ell }^{\prime })}⟩={\sigma }_{\ell }{u}_i^{(\ell )}.$$

$$N$$

$$N=O\left(\frac{1}{{ϵ}^2}\frac{\Vert \mathbf{\nu }{\Vert }^2}{{\nu }_{\ell }^2}\right).$$

$$\frac{\Vert \mathbf{\sigma }{\Vert }^2}{{\sigma }_{\text{min}}^2}=\frac{\Vert \mathbf{\sigma }{\Vert }^2}{{\sigma }_{\text{max}}^2}\frac{{\sigma }_{\text{max}}^2}{{\sigma }_{\text{min}}^2}=\frac{\Vert \mathbf{\sigma }{\Vert }^2}{{\sigma }_{\text{max}}^2}{\kappa }^2.$$

$$\Vert {\mathbf{\sigma }}^{\prime }{\Vert }^2=({\sigma }_{\text{max}}^{\prime }{)}^2\left(k-1+\frac{1}{\kappa }\right)=1$$

$$\Rightarrow ({\sigma }_{\text{max}}^{\prime }{)}^2=\frac{1}{k-1+\frac{1}{\kappa }}=O(1/k),$$

$$\frac{\Vert \mathbf{\sigma }{\Vert }^2}{{\sigma }_{\text{max}}^2}{\kappa }^2=\frac{1}{({\sigma }_{\text{max}}^{\prime }{)}^2}{\kappa }^2=O(k{\kappa }^2).$$

$$N=O\left(\frac{k^2{\kappa }^2{\kappa }_{\beta }^2}{{ϵ}^2}\right),$$

$$N=O\left(\frac{k{\kappa }_{\nu }^2}{{ϵ}^2}\right).$$

$$\mathbf{w}:=\sum _{\ell =1}^k\frac{{\tilde{\lambda }}_{\ell }}{{\tilde{\sigma }}_{\ell }}{\mathbf{w}}^{(\ell )}.$$

$$V^{(\ell )}=\frac{1}{2^n}\sum _{\mathtt{y},\mathtt{z}\in \{0,1{\}}^n}(-1{)}^{{\mathtt{x}}^{(\ell )}\cdot (\mathtt{y}\oplus \mathtt{z})}{e}^{(\mathtt{y})}{e^{(\mathtt{z})}}^T,$$

$$A:=$$

$$p(\mathtt{z})=\frac{|a_{\mathtt{y},\mathtt{z}}{|}^2}{{\max }_{{\mathtt{z}}^{\prime }}|a_{\mathtt{y},{\mathtt{z}}^{\prime }}{|}^2}.$$

$$C_{s,t}=\frac{1}{\sqrt{2^{k}c}}\frac{\Vert \mathbf{a}\Vert }{\sqrt{{\sum }_{s^{\prime }=1}^r(a_{i_{s^{\prime }},j_t}{)}^2}}{a}_{i_s,j_t},$$

$${\chi }_{\mathtt{y},\mathtt{z}}=\frac{\Vert \mathbf{a}{\Vert }^2}{2^k}\frac{(-1{)}^{{\mathtt{x}}^{(\ell )}\cdot \mathtt{z}}}{{\sum }_{\ell =1}^k{\sigma }_{\ell }(-1{)}^{{\mathtt{x}}^{(\ell )}\cdot (\mathtt{y}\oplus \mathtt{z})}}\sum _{\ell =1}^k{\beta }_{\ell }(-1{)}^{{\mathtt{x}}^{(\ell )}\cdot \mathtt{y}},$$

$$\tilde{\mathbf{x}}=\frac{1}{\sqrt{2^n}}\sum _{\ell =1}^k{\tilde{\lambda }}_{\ell }{\tilde{\mathbf{v}}}^{(\ell )}.$$

$$\mathbf{r}=\frac{1}{n}\sum _{j=1}^n{\mathbf{r}}_j,\Sigma =\frac{1}{n}\sum _{j=1}^n{\mathbf{r}}_j{\mathbf{r}}_j^T.$$

$$\underset{\mathbf{w}}{\text{Minimize: }}{\mathbf{w}}^T\Sigma \mathbf{w}$$

Full text

1 Introduction

A driving force for studying quantum computing is the conviction that quantum algorithms can solve some problems more efficiently than classical methods. But the boundary between classical and quantum computing changes constantly: new classical and quantum algorithms disrupt the landscape, continuously updating the border that separates these two paradigms. This interplay is fruitful; it allows classical and quantum algorithm developers to feed off each other’s innovations and advance our understanding of the limits of computation. Indeed, several classical algorithms have been reported that rely at least partially on developments in quantum computing [20, 12, 6, 32, 26, 1, 3, 31, 5].

Quantum algorithms for linear algebra are a flagship application of quantum computing, particularly due to their relevance in machine learning [14, 17, 23, 16, 2]. These algorithms typically scale polylogarithmically with dimension, which, at the time they were reported, implied an asymptotic exponential speedup compared to state-of-the-art classical methods. For this reason, significant interest has been generated in the dequantization approach that led to breakthrough quantum-inspired classical algorithms for linear algebra problems with sublinear complexity [31, 30, 10, 4]. These dequantized algorithms work for general low-rank matrices, whereas quantum computers still exhibit an exponential speedup over all known classical algorithms for sparse, full-rank matrix problems, including the quantum Fourier transform, eigenvector and eigenvalue analysis, linear systems, and others. Dequantized algorithms for such problems – high-rank matrix inversion, for example – would imply that classical computers can efficiently simulate quantum computers, i.e., BQP = P, which is not currently considered to be likely.

The proof techniques used to determine the complexity of quantum-inspired algorithms lead to bounds that suggest that runtimes might be prohibitively large for practical applications: the proven complexity of the linear systems algorithm is $\tilde{O}(\kappa^{16}k^{6}\|A\|^{6}_{F}/\varepsilon^{6})$ [10], while for recommendation systems it is $\tilde{O}(k^{12}/\varepsilon^{12})$ [31]111We have taken $\sigma=O(\|A\|_{F}/\sqrt{k})$ as stated in [31] and ignored scaling with respect to the error parameter $\eta$.. Here, $A$ is the input matrix, $\|A\|_{F}$ is the Frobenius norm of $A$, $\kappa$ is the condition number, $k$ is the rank, and $\varepsilon$ is the precision of the algorithm. These results raise a few immediate questions: are these complexity bounds the consequence of proof techniques, or do they reflect a fundamental limitation on the capabilities of these algorithms? How do these algorithms actually perform in practice?

We address these questions by analyzing, implementing, and testing quantum-inspired algorithms for linear algebra. First, we perform a theoretical analysis aimed at identifying potential bottlenecks in their practical implementation. This allows us to anticipate the regimes where quantum-inspired algorithms may outperform previously-known classical methods. We then test the algorithms on artifical examples based on random matrices, where it is possible to control the dimension, rank, and condition number of the input matrix. We conclude by testing the algorithms on problems of practical interest. Specifically, we run the quantum-inspired algorithm for linear systems of equations applied to portfolio optimization on stocks from the S&P 500, and analyze the algorithm for recommendation systems on the MovieLens dataset [13]. Based on our analysis and tests, we find that, provided that the input matrices have small rank and condition number, quantum-inspired algorithms can perform well in practice: they provide good approximations in reasonable times, even for very large-dimensional problems. This shows that quantum-inspired algorithms can have significantly smaller runtimes than their worst-case complexity bounds would suggest, although we find evidence that the dependence on the error $\varepsilon$ may be tight for the linear systems algorithm. For problems whose input matrices have larger ranks and condition numbers, our implementation of the algorithms struggles noticeably to produce good results.

In the following sections, we begin by providing a brief description of quantum-inspired algorithms for linear algebra. This high-level summary may prove useful for those wishing to become more familiar with these techniques. We then analyze each of the main steps of the algorithm in further detail, focusing on identifying potential practical bottlenecks. We continue by implementing the algorithms and benchmarking their performance on artificially-generated problems and on real-world data. We conclude by summarizing the implications of our analysis.

2 Quantum-inspired algorithms for linear algebra

In this section, we give an overview of quantum-inspired algorithms for linear systems of equations and recommendation systems. More details can be found in Refs. [31, 30, 10, 4]. A unique feature of our description is that we view both algorithms as specific examples of a more general method to sample from vectors expressed in terms of the singular value decomposition (SVD) of an input matrix. Henceforth, capital letters are employed to denote matrices and bold letters to denote column vectors. For example, a linear system of equations is written as $A\,\bm{x}=\bm{b}$, where $A\in\mathbb{R}^{m\times n}$, $\bm{x}=(x_{1},x_{2},\ldots,x_{n})^{T}$, and $\bm{b}=(b_{1},b_{2},\ldots,b_{m})^{T}$.

We consider problems of the following form: Given an $m\times n$ matrix $A\in\mathbb{R}^{m\times n}$ with SVD

[TABLE]

the goal is to sample entries of the $n$-dimensional vector

[TABLE]

with respect to the length-square probability distribution $p_{x}(i)=x_{i}^{2}/\|\bm{x}\|^{2}$. The coefficients $\lambda_{\ell}$ depend on the matrix $A$ and possibly on other inputs. For linear systems of equations, $\bm{x}$ is the solution vector given by $\bm{x}=A^{+}\bm{b}$, with $A^{+}$ the Moore-Penrose pseudoinverse. As discussed in Ref. [10], the coefficients $\lambda_{\ell}$ of Eq. (2) are given by the inner product

[TABLE]

Similarly, for recommendation systems, $A$ is the preference matrix, whose entries $A_{ij}$ denote the rating that user $i$ has given to product $j$. In this case, the solution vector $\bm{x}$ is the $i$-th row of a low-rank approximation of $A$, which contains all the preferences of user $i$. The coefficients are given by

[TABLE]

where $A_{i}$ is the $i$-th row of $A$ [31].

Quantum-inspired algorithms assume that the entries of $A$ are given in a way that allows length-square sampling to be performed. This can be accomplished, for example, if the entries of $A$ are stored in a data structure of the form proposed by [16], which stores both entries of vectors and their sub-norms. Length-square sampling allows us to preferentially sample the rows of $A$ with large norm, and the large entries within each row. Given access to length-squared sampling, quantum-inspired algorithms solve low-rank linear algebra problems in three main steps:

Approximate SVD: The Frieze-Kannan-Vempala (FKV) algorithm [9] is used to compute approximate singular values $\tilde{\sigma}_{\ell}$ and approximate right singular vectors $\tilde{\bm{v}}^{(\ell)}$. The singular vectors are not calculated in their entirety; instead, a description is found allowing efficient computation of any given entry $v_{i}^{(\ell)}$. 2. 2.

Coefficient estimation: Based on the results of step 1, the coefficients $\tilde{\lambda}_{\ell}$ are approximated using Monte Carlo estimation techniques. 3. 3.

Sampling solution vectors: Using Eq. (2) and the results from steps 1 and 2, rejection sampling is employed to perform length-square sampling from the approximate vectors $\tilde{\bm{x}}=\sum_{\ell=1}^{k}\tilde{\lambda}_{\ell}\bm{\tilde{v}}^{(\ell)}$.

Before describing each of these steps in more detail, it is important to understand that the innovations of quantum-inspired algorithms are steps 2 and 3. In these steps, coefficient estimation and sampling from the solution vectors is performed in time $O(\text{poly}(k,\kappa,\varepsilon,\log n,\log m))$. By contrast, the strategy originally suggested in Ref. [9] is simply to compute coefficients and solution vectors directly from the approximate SVD of step 1, which requires time $O(kn)$. Asymptotically, quantum-inspired algorithms thus achieve an exponential asymptotic speedup in dimension, from $O(n,m)$ to $\text{polylog}(n,m)$, at the cost of a polynomial dependency on other parameters. In practice, the direct calculation can be done extremely fast since it scales only linearly with dimension, meaning that quantum-inspired algorithms require very large-dimensional problems before it becomes preferable to employ their sampling techniques.

2.1 Approximate SVD

The Frieze-Kannan-Vempala (FKV) algorithm [9] is a method to compute low-rank approximations of matrices. The FKV algorithm is a pioneering example of a randomized method for computing approximate matrix decompositions. These techniques have been studied extensively and several improvements have since been described [8, 27, 34, 24, 25, 19, 33, 7], see Ref. [11] for a survey of these techniques. As described in Ref. [11], these algorithms are variants of the same main strategy: (i) preprocess the matrix to calculate sampling probabilities, (ii) generate samples from the matrix, and (iii) use linear algebra techniques to postprocess the samples to compute a final approximation. It is an interesting question to understand to what extent other randomized methods for approximate matrix decompositions can be used to improve quantum-inspired algorithms. For concreteness and to allow a direct connection to Refs. [31, 30, 10, 4], in this work we focus on analyzing the FKV algorithm and implement quantum-inspired algorithms that rely on FKV. The conclusions and performance of the implementation are therefore subject to the specific properties of this algorithm.

The main idea behind FKV is the following: instead of performing a singular value decomposition of the large input matrix $A\in\mathbb{R}^{m\times n}$, a smaller matrix $C\in\mathbb{R}^{r\times c}$ is constructed by sampling $r$ rows and $c$ columns of $A$. The only restriction on $r$ and $c$ is that they have to be larger than the rank of the low-rank approximation of $A$: in particular, they can be independent of $m$ and $n$, i.e., not even $O(\log m)$, $O(\log n)$. The matrix $C$ captures the important features of $A$, making it possible to perform an SVD of the smaller matrix $C$ to retrieve information about the singular values and vectors of $A$. Since the rank of $C$ is bounded by its dimension, this approach is only possible when $A$ is low-rank or when a low-rank approximation of $A$ is sufficient.

In more detail, the first step of FKV is to perform a pass through the matrix $A$ and compute the Frobenius norm of each row $\|A_{i}\|_{F}$, as well as the norm of the matrix $\|A\|_{F}$. These quantitites define a length-square distribution over rows: $p(i)=\|A_{i}\|^{2}_{F}/\|A\|_{F}^{2}$. For each row $i$, a length-square distribution over its entries is also computed, $q_{i}(j)=\|A_{ij}\|^{2}_{F}/\|A_{i}\|_{F}^{2}$. To sample from these distributions, it is possible to employ a specialized data structure as suggested in Refs. [9, 16, 31] that allows sampling in time that is logarithmic in dimension. The algorithm proceeds as follows:

Sample $r$ row indices $i_{1},i_{2},\ldots,i_{r}$ from the row distribution $p(i)$. For each row index, select the row $A_{i_{s}}$ and renormalize it as $R_{i_{s}}=\frac{\|A\|_{F}}{\sqrt{r}\|A_{i_{s}}\|}A_{i_{s}}$. This defines an $r\times n$ matrix $R$ consisting of the rescaled rows $R_{i_{s}}$. 2. 2.

Select an index $s\in\{1,2,\ldots,r\}$ uniformly at random, then sample a column index $j$ from the column distribution $q_{i_{s}}(j)$. Repeat this a total of $c$ times to obtain column indices $j_{1},j_{2},\ldots,j_{c}$. For each column index, select the column $R_{\cdot,j_{t}}$ and renormalize it as $C_{\cdot,j_{t}}=\frac{\|A\|_{F}}{\sqrt{c}\|R_{\cdot,j_{t}}\|}R_{\cdot,j_{t}}$. This defines a matrix $C$ consisting of the rescaled columns $C_{\cdot,j_{t}}$. 3. 3.

Compute the singular value decomposition of $C$.

These steps are illustrated in Fig. 1. The FKV algorithm performs a useful tradeoff: it is possible to perform SVDs of a significantly smaller matrix $C$ at the cost of obtaining only approximations of the singular values and vectors of the original matrix $A$. We denote by $\tilde{\sigma}_{\ell}$ the singular values of $C$ and by $\bm{\omega}^{(\ell)}$ its left singular vectors. The FKV algorithm directly provides approximations of singular values, $\tilde{\sigma}_{\ell}\approx\sigma_{\ell}$ for $\ell=1,2,\ldots,k$. Approximate right singular vectors $\bm{\tilde{v}}^{(\ell)}$ of $A$ are obtained using the expression

[TABLE]

while approximate left singular vectors $\bm{\tilde{u}}^{(\ell)}$ are given by

[TABLE]

In existing classical methods, these vectors are computed directly, whereas in quantum-inspired algorithms it suffices to query their entries.

The quality of the FKV approximations depends on the number of sampled rows $r$ and columns $c$. The challenge is to find values of $r$ and $c$ such that sufficiently good approximations are obtained, while also ensuring that the complexity of computing the SVD of matrix $C$ is significantly less than for the original input matrix $A$, i.e., such that $r\ll m$ and $c\ll n$. Since $r,c$ can be much smaller than $m,n$, a substantial runtime reduction can potentially be achieved.

FKV is the core of quantum-inspired algorithms for linear algebra. Without the use of randomized methods, computing the SVD of an $m\times n$ matrix $A$ requires $O(\min\{m^{2}n,mn^{2}\})$ time using naive matrix multiplication. From that point onwards, computing coefficients $\lambda$ and the solution vector $\bm{x}$ takes only linear time in $m,n$. In FKV, we instead calculate the SVD of an $r\times c$ matrix, which requires $O(\min\{r^{2}c,rc^{2}\})$ time. When $r,c$ are much smaller than $m,n$, this is where the real savings of quantum-inspired algorithms happen. As explained previously, the subsequent steps of coefficient estimation and sampling from the solution vector are used to reduce a linear runtime to polylogarithmic runtime. As we discuss later, for practical problem sizes that are far from the asymptotic limit, it is preferable to compute coefficients and solution vectors explicitly by employing Eqs. (5) and (6) starting from the approximate SVD of matrix $C$ in the FKV algorithm.

2.2 Coefficient estimation

The coefficients $\lambda$ appearing in Eq. (2) are inner products between vectors, multiplied by a power of the singular values in the case of linear systems. Quantum-inspired algorithms compute approximations of these coefficients using Monte Carlo estimation. In general, a coefficient $\lambda$ is given by

[TABLE]

for some appropriate vectors $\bm{y},\bm{z}$, as in Eqs. (3) and (4). The strategy to estimate this inner product is to perform length-square sampling from one of the vectors. Without loss of generality, we take this vector to be $\bm{y}$. Define the random variable $\chi$ that takes values $\chi_{i}=y_{i}z_{i}/p_{x}(i)$, where the indices $i$ are sampled from the length-square distribution $p_{y}(i)=y_{i}^{2}/\|\bm{y}\|^{2}$. The expectation value of the random variable satisfies

[TABLE]

Similarly, the second moment is

[TABLE]

and $\sigma_{\chi}^{2}=\mathbb{E}(\chi^{2})-\mathbb{E}(\chi)^{2}$ is the variance of $\chi$. The strategy to estimate coefficients $\lambda$ is to draw $N$ samples $\chi^{(1)},\chi^{(2)},\ldots,\chi^{(N)}$ from $\chi$ and compute the unbiased estimator $\hat{\lambda}=\frac{1}{N}\sum_{j=1}^{N}\chi^{(j)}\approx\lambda$. This constitutes a form of importance sampling in Monte Carlo estimation: large entries of $\bm{y}$ are preferably selected since they contribute more significantly to the inner product. The error in the estimation is quantified by the ratio between the standard deviation and the mean of the estimator: $\epsilon:=\sqrt{\text{Var}(\hat{\lambda})}/|\lambda|$. The variance of the estimator is $\text{Var}(\hat{\lambda})=\sigma_{\chi}^{2}/N$, leading to a precision $\epsilon=\sigma_{\chi}/(|\lambda|\sqrt{N})$ in the estimator. This implies that the number of samples $N$ needed to, with high probability, achieve a precision $\epsilon$ is

[TABLE]

where $\theta$ is the angle between the vectors $\bm{y}$ and $\bm{z}$.

2.3 Complexity of coefficient estimation

Quantum-inspired algorithms differ from the FKV algorithm because they estimate the coefficients $\lambda$ rather than calculating them directly from the approximate SVD of the input matrix. Therefore, it is worthwhile to examine the complexity of coefficient estimation step to determine the regime where quantum-inspired algorithms become beneficial. For $n$-dimensional random vectors, the expectation value of the angle between two vectors satisfies $\mathbb{E}\left[\frac{1}{\cos^{2}\theta}\right]=n$, so we can anticipate that the number of samples needed may scale linearly with dimension. If $\frac{1}{\cos^{2}\theta}=\text{poly}(m,n)$, the algorithm no longer has the desired polylogarithmic complexity. In fact, as we show in Appendix A, regardless of what coefficient estimation strategy is employed, in a worst-case setting, inner products of vectors cannot be approximated in sublinear time. However, for matrices of low rank and condition number, as we show below, it is possible to achieve good estimation of coefficients.

For linear systems of equations, the coefficients $\lambda_{\ell}$ are given by

[TABLE]

In the quantum-inspired algorithm, we do not have direct access to the length-square distributions of either $A^{T}\bm{b}$ or $\bm{v}^{(\ell)}$. Instead, the strategy in Ref. [10] is to express the coefficients as $\lambda_{\ell}=\frac{1}{\sigma^{2}_{\ell}}\text{Tr}(A^{T}\bm{b}\,{\bm{v}^{(\ell)}}^{T})$ and sample the entries $A_{ij}$ of $A$. In this case the the number of samples needed in the estimation is [10]:

[TABLE]

where we have defined $\bm{\sigma}:=(\sigma_{1},\sigma_{2},\ldots,\sigma_{k})$, $\beta_{\ell}:=\langle\bm{u}^{(\ell)},\bm{b}\rangle$ such that $\bm{b}=\sum_{\ell=1}^{k}\beta_{\ell}\bm{u}^{(\ell)}$, and we have used Eq. (12) as well as the identity $\|A\|_{F}^{2}=\|\bm{\sigma}\|^{2}$. In the special cases of recommendation systems and portfolio optimization, the coefficients are proportional to an inner product $\langle A_{i}^{T},\bm{v}^{(\ell)}\rangle$. Writing $A_{i}=\sum_{\ell^{\prime}=1}^{k}\sigma_{\ell^{\prime}}u_{i}^{(\ell^{\prime})}{\bm{v}^{(\ell^{\prime})}}^{T}$ we have

[TABLE]

The number of samples then satisfies

[TABLE]

Define the $k$-dimensional vectors $\bm{\sigma}:=(\sigma_{1},\sigma_{2},\ldots,\sigma_{k})^{T}$, $\bm{\mu}:=(u_{i}^{(1)},u_{i}^{(2)},\ldots,u_{i}^{(k)})^{T}$ and their Hadamard product $\bm{\nu}=\bm{\sigma}\circ\bm{\mu}$ with entries $\nu_{\ell}=\sigma_{\ell}u_{i}^{(\ell)}$. It follows that

[TABLE]

In all cases, the number of samples depends on the ratio between the norm squared of a $k$-dimensional vector and the square of the $\ell$-th entry of that vector. We now determine the worst-case scaling for these ratios. We consider the ratio $\|\bm{\sigma}\|^{2}/\sigma_{\ell}^{2}$ for concreteness; the same will apply for $\bm{b}$ and $\bm{\nu}$. The largest ratio occurs for the smallest singular value $\sigma_{\text{min}}$ and we have

[TABLE]

Similarly, this ratio is largest when $\sigma_{\text{max}}$ is as small as possible. This occurs when the largest $k-1$ singular values are equal to each other. In that case, defining $\sigma^{\prime}_{\ell}:=\sigma_{\ell}/\|\bm{\sigma}\|$ such that $\|\bm{\sigma^{\prime}}\|=1$ we have

[TABLE]

and therefore

[TABLE]

Defining $\kappa_{\beta}:=\beta_{\text{max}}/\beta_{\text{min}}$ and $\kappa_{\nu}:=\nu_{\text{max}}/\nu_{\text{min}}$, where $\beta_{\text{max/min}}$ and $\nu_{\text{max/min}}$ are the largest and smallest entries of their corresponding vectors, we conclude that the number of samples required to estimate coefficients in quantum-inspired algorithms for linear systems scales as

[TABLE]

and for recommendation systems as

[TABLE]

These asymptotic formulas – or more precisely Eqs. (13) and (15) – can be used to estimate the problem sizes for which the complexity of coefficient estimation is smaller than a direct calculation. As an illustration, in linear systems, setting $k=\kappa=\kappa_{\beta}=100$ means that an order of $N=10^{16}$ samples are needed to estimate the coefficients with an error of $\epsilon=10^{-2}$. Finally, it is important to understand that a large number of samples can improve the precision of the estimation, but not its accuracy: if the singular values and vectors obtained from FKV have large errors, the expectation value of the estimator will not coincide with the actual value of the coefficient.

2.4 Sampling solution vectors

The approximate SVD and coefficient estimation steps respectively provide approximate singular values $\tilde{\sigma}_{\ell}$ and approximate coefficients $\tilde{\lambda}_{\ell}$. The approximate right singular vectors of $A$ are not constructed explicitly, but are implicitly defined by Eq. (5), namely $\bm{\tilde{v}}^{(\ell)}=\frac{1}{\tilde{\sigma}_{\ell}}R^{T}\bm{w}^{(\ell)}$. This allows an implicit calculation of the approximate solution vector $\tilde{\bm{x}}=R^{T}\bm{w}$, where

[TABLE]

The challenge is to sample from these vectors using only query access to the entries of $\bm{w}$ and $R$. Quantum-inspired algorithms achieve this using a well-known technique known as rejection sampling. The steps are as follows [31]:

Sample a row index $i$ uniformly at random. 2. 2.

Sample a column index $j$ from the length-square distribution $q_{i}(j)=|R_{ij}|^{2}/\|R_{i}\|^{2}$. 3. 3.

Output $j$ with probability $\frac{|\langle\bm{w},R_{\cdot,j}\rangle|^{2}}{\|R_{\cdot,j}\|^{2}\|\bm{w}\|^{2}}$, or sample $j$ again.

The expected number of repetitions before outputting a column index $j$ is $r\|\bm{w}\|^{2}/(\sum_{j=1}^{n}|\langle\bm{w},R_{\cdot,j}\rangle|^{2})=r\|\bm{w}\|^{2}/\|\bm{\tilde{x}}\|^{2}$. The steps of rejection sampling can be performed quickly. The overhead arises from calculating the inner product $\langle\bm{w},R_{\cdot,j}\rangle$, which requires $O(rk)$ time; significantly less than computing the approximate SVD. As we illustrate in the following sections, the number of expected repetitions is usually moderate in practice, so sampling from the approximate solution vector can be performed in comparatively less time than other steps in the algorithm.

3 Numerical benchmarking

In the previous section we performed a theoretical analysis of quantum-inspired algorithms. Here, we complement that analysis by benchmarking the algorithms on specific problems. The main goal of this numerical implementation is to individually examine the performance of each step of the algorithm and use this information to reveal which steps contribute most to the overall runtime as well as which ones lead to the most significant source of errors. A secondary goal is to gain insights on the actual scaling of the algorithm with respect to its various parameters.

As an initial benchmark, we study the algorithm for linear systems of equations on artificial problems of extremely large dimension. We design matrices such that performing length-square sampling from its rows and columns, as well as all other steps of the quantum-inspired algorithm, can be done using only query access to its entries. This allows to study the performance of the algorithms in their intended asymptotic regime, without handling large datasets directly. We then test the algorithm for linear systems on Gaussian random matrices. This choice is made to allow full control over properties of the input matrices: dimension, rank, and condition number, while collecting large amounts of statistics. Afterwards, we apply the algorithms to portfolio optimization and movie recommendations. All algorithms were implemented in Python and the source code is available at https://github.com/XanaduAI/quantum-inspired-algorithms.

3.1 Design principles

In conducting numerical experiments, choices must be made regarding how the algorithms are implemented. Here we describe the principles that guide the numerical experiments performed to benchmark quantum-inspired algorithms. The first choice is the programming language. To implement resource-intensive algorithms, it is common to employ compiler-based languages such as C++. However, our main goal is not to develop high-performance implementation of quantum-inspired algorithms. Rather, the goal is to identify practical bottlenecks. Crucially, it is also important to share code that is easily-reproducible by the quantum information community. For these reasons, the algorithm is in Python. The runtimes reported in this work can therefore be likely improved by using a language more suitable for high-performance computation, and by employing specialized hardware.

In implementing the algorithms, whenever possible, we have been faithful to the descriptions in Refs. [10] and [31]. These algorithms require the use of a specialized data structure to ensure polylogarithmic scaling of length-square sampling. As shown in Fig. 2, for the matrix sizes that we study in this work, which have dimensions smaller than $10^{5}$, using built-in sampling functions is orders-of-magnitude faster than employing our implementation of the data structure. Therefore, to speed up the practical runtime of the algorithm for these problems, we employ direct sampling from the distribution using built-in functions.

Finally, although it is natural to consider the total variation distance as a way to quantify errors, we identify that the mean relative error, defined as $\eta=\sum_{i}\frac{|x_{i}-\tilde{x_{i}}|}{|x_{i}|}$ for two vectors $x$ and $\tilde{x}$, is more informative in determining the quality of approximations. Moreover, it can be simultaneously used to measure errors in approximating singular values, coefficients, and solution vectors.

3.2 High-dimensional problems

Quantum-inspired algorithms for linear algebra are aimed at tackling extremely large datasets, but it is very challenging to test them in such scenarios. Instead, we design a class of problems where it is possible to run quantum-inspired algorithms using only black-box access to the entries of the input matrix and vector. We focus on the algorithm for linear systems and report the main formulas, whose derivations can be found in Appendix B.

Define the matrices

[TABLE]

where $\mathtt{x,y,z}$ are $n$-bit strings, the vectors $e^{(\mathtt{y})}$ form a basis of the $2^{n}$-dimensional vector space, and the $\mathtt{x}^{(\ell)}$ are distinct $n$-bit strings for $\ell=1,2,\ldots,k$. We define the input matrix as

[TABLE]

where $\sigma_{1},\sigma_{2},\ldots,\sigma_{k}$ are the singular values and we have implicitly defined $a_{\mathtt{y},\mathtt{z}}=\sum_{\ell=1}^{k}\sigma_{\ell}(-1)^{x^{(\ell)}\cdot(y\oplus z)}$. The right and left singular vectors of $A$ are the same, given by $v^{(\ell)}=\frac{1}{\sqrt{2^{n}}}\sum_{\mathtt{y}\in\{0,1\}^{n}}(-1)^{\mathtt{x}^{(\ell)}\cdot\mathtt{y}}e^{(\mathtt{y})}$. As we show in the Appendix B, by construction, all rows of $A$ have the same norm, therefore length-square sampling can be done by choosing rows uniformly at random. For length-square sampling of columns, we use rejection sampling: given a row $\mathtt{y}$, we select a column index $\mathtt{z}$ uniformly at random and accept it with probability

[TABLE]

The optimization is trivial since $\max_{\mathtt{z}^{\prime}}|a_{\mathtt{y,z^{\prime}}}|^{2}=\sum_{\ell=1}^{k}(\sigma_{\ell})^{2}$. Once row indices $i_{1},\ldots,i_{r}$ and column indices $j_{1},\ldots,j_{c}$ have been sampled, the entries of the $r\times c$ matrix $C$ in FKV can be expressed as

[TABLE]

where $\bm{a}=(a_{1},a_{2},\ldots,a_{2^{k}})^{T}$ is a vector containing all the possible values of the entries $a_{\mathtt{y},\mathtt{z}}$. The approximate singular values and vectors of $A$ can be reproduced from the SVD of $C$ as detailed in previous sections. To estimate the coefficient $\lambda_{\ell}=\frac{1}{\sigma_{\ell}^{2}}\langle\bm{v}^{(\ell)},A^{T}\bm{b}\rangle$, we sample the random variable

[TABLE]

and use its sample mean as an estimator for the coefficients. Finally, the approximate solution vector $\bm{\tilde{x}}$ is given by

[TABLE]

Any entry of the approximate solution vector can be computed efficiently from the entries of the approximate right singular vectors $\bm{\tilde{v}}^{(\ell)}$, which can be obtained from the approximate SVD of $C$.

The crucial property of this construction is that the only dependency on dimension comes from the sampling of rows and columns, which can be done efficiently even for extremely high-dimensional problems. The limitations arise in computing the matrix $C$, which takes $O(r^{2}c)$ time, and from computing the $2^{k}$ different matrix elements. This limits the values of $r,c$ and $k$ that can be explored in this approach. We fix the input matrices to have an extremely large dimension of $2^{50}\approx 10^{15}$. We set $r=c=150$, the largest that could be handled in a reasonable amount of time. This means that we are aiming to approximate the SVD of the input matrix by using a smaller matrix whose dimension differs by thirteen orders of magnitude. The bit strings $\mathtt{x}^{(\ell)}$ that define each singular vector are chosen uniformly at random. We study three examples with $k=\kappa=\kappa_{\beta}=3,5,10$, employing $N=10^{4}$ samples in the estimation of the coefficients. We report the average and standard deviation of the errors $\eta_{\sigma}=\frac{1}{k}\sum_{\ell=1}^{k}\frac{|\tilde{\sigma}_{\ell}-\sigma_{\ell}|}{\sigma_{\ell}}$ and $\eta_{\lambda}=\frac{1}{k}\sum_{\ell=1}^{k}\frac{|\tilde{\lambda}_{\ell}-\lambda_{\ell}|}{|\lambda_{\ell}|}$ over ten repetitions of the algorithms. For vectors, we evaluate the first $L$ entries and similarly report the average and standard deviation of the errors $\eta_{\bm{v}}=\frac{1}{kL}\sum_{i=1}^{L}\sum_{\ell=1}^{k}\frac{|v^{(\ell)}_{i}-\tilde{v}^{(\ell)}_{i}|}{|v^{(\ell)}_{i}|}$ and $\eta_{\bm{x}}=\frac{1}{kL}\sum_{i=1}^{L}\sum_{\ell=1}^{k}\frac{|x_{i}-\tilde{x}_{i}|}{|x_{i}|}$. The results are summarized in Table 1.

These results show that for problems with the appropriate properties – namely very low rank and condition number – quantum-inspired algorithms can provide good estimates even for extremely high-dimensional problems. Importantly, these results are obtained in considerably less time than would be expected from the theoretical complexity bound $\tilde{O}(\kappa^{16}k^{6}\|A\|^{6}_{F}/\varepsilon^{6})$. Even setting $\|A\|_{F}=1$ and $\varepsilon=1$, for $k=\kappa=\kappa_{\beta}=10$, the bound would suggest that an order of $10^{22}$ operations would be needed; this is not compatible with the runtimes we experience. This suggests that these theoretical bounds do not in fact reflect the practical runtime of the algorithms and therefore should not be used to evaluate their performance. Nevertheless, our findings do reflect that the quality of the estimates worsens considerably as rank and condition number are increased.

3.3 Random matrices

We study the quantum-inspired algorithm for linear systems of equations $A\bm{x}=\bm{b}$ for randomly chosen $A$ and $\bm{b}$. To generate a random matrix $A$ of dimension $m\times n$, rank $k$, and condition number $\kappa$, instead of creating $A$ directly, we build the components of its singular value decomposition in matrix form, $A=U\Sigma V$. Here, $U$ is an $m\times k$ matrix whose columns are the left singular vectors of $A$. Similarly, $V$ is a $k\times n$ matrix whose rows are the right singular vectors, and $\Sigma$ is a $k\times k$ diagonal matrix whose entries are the corresponding singular values of $A$, which has rank $k$ by construction. The details of numerical methods used to generate matrices $U$, $V$ and $\Sigma$ are explained in Appendix C. The vector $\bm{b}$ has been chosen to have the form $\bm{b}=\sum_{\ell=1}^{k}\beta_{\ell}\bm{u}^{(\ell)}$, with coefficients $\beta_{\ell}$ drawn from the standard normal distribution $\mathcal{N}(0,1)$.

We consider matrices consisting of $m=40,000$ rows and $n=20,000$ columns. These numbers are chosen to deal with the general case of rectangular matrices, while working at the limit of dimensions for which it would still be possible to calculate SVDs in a reasonable amount of time. This means, for simplicity, that in this setting we are not required to use specialized data structures to perform length-square sampling: we can instead use fast built-in numerical sampling algorithms in Python which are very fast in practice. More precisely, generating a single sample from a distribution of dimension $40,000$ using built-in algorithms takes approximately $2\times 10^{-3}$ seconds on a standard desktop computer. This is even more efficient when several samples are taken: generating one million samples from such vectors requires roughly $0.2$ seconds.

In Figs. 3(a)-3(c) we monitor, respectively, the mean relative error $\eta_{\sigma}=\frac{1}{k}\sum_{\ell=1}^{k}\frac{|\tilde{\sigma}_{\ell}-\sigma_{\ell}|}{\sigma_{\ell}}$ in estimating the singular values as well as the errors $\eta_{A}=\|\tilde{A}-A\|_{F}/\|A\|_{F}$ and $\eta_{A}^{+}=\|\tilde{A}^{+}-A^{+}\|_{F}/\|A^{+}\|_{F}$ of the reconstructed matrices $A$ and $A^{+}$. The approximation is accurate across all singular values even for a relatively small number of sampled rows and columns.

Empirically, the error in the approximation scales roughly as $1/\sqrt{r}$ for these random matrices, which is in agreement with the matrix Chernoff bound appearing in Theorem 3 of Ref. [10]. We also observe that the reconstructed matrix $\tilde{A}$ has smaller errors compared to the approximated pseudo-inverse matrix $A^{+}$. This is because in reconstructing matrix $\tilde{A}$, the largest singular values – which carry the dominant contributions – are approximated better by the FKV algorithm.

The main results are shown in Fig. 4, where we plot the error $\eta_{\bm{x}}$ of the approximated solution vector $\tilde{\bm{x}}$ calculated as the median of the relative errors $|\tilde{x}_{i}-x_{i}|/|x_{i}|$ with respect to all $n$ entries of vector $\tilde{\bm{x}}$. First, we investigate in Fig. 4(a) how the error $\eta_{\bm{x}}$ changes as we increase the number of sampled rows $r$ and columns $c=r$. All results have been obtained by taking $N=10^{4}$ samples to estimate the coefficients $\lambda_{\ell}$ defined in Eq. (12). As expected, $\eta_{\bm{x}}$ decreases as the number of sampled rows and columns is increased. We have also superimposed error bars to show the standard deviation over 10 independent repetitions of the FKV algorithm.

In Fig. 4(b) we set $r=c=4250$ and $\kappa=5$ and investigate the error of the solution vector as a function of the rank $k$ of matrix $A$. This plot reveals that the error $\eta_{\bm{x}}$ increases linearly with rank. For $k=5$, a remarkably small relative error of $8.7\%$ is found for the entries of the approximated vector. However, we observe that the accuracy of the algorithm deteriorates rapidly for larger values of the rank. A similar trend is observed in Fig. 4(c) for the solution vector error as a function of the condition number $\kappa$. Small errors of the order of $10\%$ are only obtained for small condition numbers, while large values of $\eta_{\bm{x}}\geq 100\%$ already appear for $\kappa\geq 100$. Note that this is the behaviour we would expect from Eq. (21) if the dominant source of errors originates from coefficient estimation. Similar plots for the error in the approximate SVD of matrix $A$ and in coefficients estimation are shown in Appendix C. Overall, we conclude that for random matrices, errors decrease as $1/\sqrt{r}$, and increase linearly with both rank and condition number.

Finally, Tables 2 and 3 respectively summarize the errors and running time of the quantum-inspired algorithm for the case where the matrix $A$ has rank $k=5$ and condition number $\kappa=5$. In Table 2, we show the average relative errors of singular values $\eta_{\sigma}$ and coefficients $\eta_{\lambda}$ as defined in the previous section. We also report the relative errors $\eta_{A}=\|A-\widetilde{A}\|_{F}/\|A\|_{F}$ and $\eta_{A^{+}}=\|A^{+}-\widetilde{A}^{+}\|_{F}/\|A^{+}\|_{F}$ of the reconstructed matrices $\widetilde{A}$ and $\widetilde{A}^{+}$. Sampling $4250$ rows and columns from the original $40000$ $\times$ $20000$ matrix $A$ is enough to have a good approximation full matrix. The largest error in Table 2 corresponds to the estimated coefficients.

Table 3 shows that the time to compute the approximate SVD using the FKV algorithm is orders-of-magnitude smaller than a direct calculation. Most of the overhead in running the quantum-inspired algorithms is associated with constructing the length-square distributions and estimating the coefficients. However, it is important to keep in mind that length-square distributions can in principle be constructed on-the-fly as the data from the matrix is generated, so it is best to separate this from the runtime of the actual algorithm. For this particular case of a matrix $A$ with low rank and small condition number, we find that the quantum-inspired algorithm outperforms the direct calculation.

We can place these runtimes into perspective by attempting to identify whether they are due to our implementation or are inherent to the algorithm. There are $8\times 10^{8}$ entries in this matrix, so our implementation is taking roughly $\sim 1.8\times 10^{-6}$ seconds per entry to compute the length-squared distribution. Similarly, to estimate the coefficients, we generate $10^{4}$ samples, repeated 10 times to calculate the median. This is done for each of the five coefficients. For each such sample, the entries of the right singular vectors of the matrix $R$ must be estimated, which requires summing over the $r=4,250$ sampled rows. Multiplying these numbers gives a total of $~{}2\times 10^{9}$ operations, for roughly $2.6\times 10^{-7}$ seconds per operation. In both cases, runtimes per operation are reasonably low, showing that total runtimes are inherently large given the number of required operations.

3.4 Portfolio optimization

We study an application of the quantum-inspired algorithm for linear systems of equations to a canonical financial problem: portfolio optimization. In its simplest version, commonly known as the Markowitz mean-variance model [18], the goal is to optimally invest wealth across $n$ possible assets that are modeled only by their expected returns and correlations. Consider a vector of returns $\bm{r}_{j}=(r_{1,j},r_{2,j},\ldots,r_{n,j})^{T}$, where $r_{i,j}$ is the return of asset $i$ on the $j$-th day. The vector of expected returns $\bm{r}$ and the correlation matrix $\Sigma$ are defined as

[TABLE]

The vector $\bm{r}$ captures the expected returns for each asset, while the correlation matrix $\Sigma$ expresses fluctuations around mean values, which in the model is interpreted as risk. An optimal portfolio in this setting is one that minimizes risk for any given target expected return. This corresponds to the optimization problem

[TABLE]

where $\bm{w}$ is a portfolio allocation vector that determines how much wealth is invested in each asset. We allow $\bm{w}$ to have negative entries, which is interpreted as an indication to short-sell the corresponding asset. As discussed in Ref. [22], this problem can be equivalently cast as a linear system of equations:

[TABLE]

where $\nu$ is a Lagrange multiplier. In this case, the linear system $A\bm{x}=\bm{b}$ is given by $A=\begin{bmatrix}0&\bm{r}^{T}\\ \bm{r}&\Sigma\end{bmatrix}$, $\bm{b}=\begin{bmatrix}\mu\\ \bm{0}\end{bmatrix}$, and $\bm{x}=\begin{bmatrix}\nu\\ \bm{w}\end{bmatrix}$. To compute a solution to the linear system, we employ the Moore-Penrose pseudo-inverse $A^{+}$, which can be expressed in terms of the singular values and vectors of $A$ as $A^{+}=\sum_{\ell=1}^{k}\frac{1}{\sigma_{\ell}}\bm{v}^{(\ell)}{\bm{u}^{(\ell)}}^{T}$. The solution vector $\bm{x}$ is then

[TABLE]

where we have identified the coefficients $\lambda_{\ell}=\langle\bm{v}^{(\ell)},A^{T}\bm{b}\rangle/\sigma_{\ell}^{2}$ as in Eq. (3). The vector $\bm{b}$ has only one non-zero entry, so it holds that $A^{T}\bm{b}=\mu A_{1}$, where $A_{1}$ is the first row of $A$.

To test the algorithm in this setting, we employ publicly-available pricing data for the stocks comprising the S&P 500 stock index during the five-year period 2013-2018 [21]. Since not all stocks remain in the index during this time, we restrict to the 474 stocks that do, leading to a matrix $A$ of dimension $475\times 475$. Returns are calculated on a daily basis based on opening price. The matrix has full rank of $k=475$, its largest and smallest singular values are respectively $\sigma_{max}=40.2$, $\sigma_{min}=1.8\times 10^{-4}$, and its condition number is $\kappa=2.23\times 10^{4}$. More details can be found in Appendix D. We set a target return $\mu$ equal to the average return over all 474 stocks in the index.

In Tables 4 we report the errors of the quantum-inspired algorithm. The SVD of matrix $A$ was approximated by sampling $340$ rows and columns from the full matrix. As in the previous example, we took $N=10^{4}$ samples to estimate the coefficients $\lambda_{\ell}$. The relative errors characterizing the approximate solution vector are considerable, with values of $\eta_{\bm{x}}=0.74$. In this example, we have calculated the approximate solution vector by using Eq. (3.4) within a low-rank approximation with $k=10$. Notice from Table 4 that the large errors $\eta_{A}^{+}$ in the reconstructed pseudo-inverse matrix and estimated coefficients $\eta_{\lambda}$ translate into large relative errors of the solution vector.

Finally, in Table 5 we report the running-times for the algorithm. Since the matrices involved in this example are small, exact and approximate SVDs can be computed quickly, but for the quantum-inspired algorithm, there is significant overhead due to coefficient estimation and sampling from the solution vector.

3.5 Movie recommendations

We analyze the performance of the algorithm for recommendation systems on the MovieLens 100K database [13], which consists of a preference matrix with 100,000 ratings from 611 users across 9,724 movies. Ratings are specified on a half-star scale in the range $[0.5,5]$. For example, an entry of 4.5 corresponds to a rating of four-and-a-half stars. The lowest possible rating is 0.5 and the zero entries in the preference matrix correspond to movies that have not been rated by the user. Since most users watch only a small fraction of available movies, the matrix is sparse. The preference matrix has full rank of $k=611$, its largest and smallest singular values are respectively $\sigma_{max}=534.4$, $\sigma_{min}=2.95$, and its condition number is $\kappa=181.2$. More details can be found in Appendix D.

The goal of the algorithm is to predict missing ratings and subsequently recommend movies that have a high predicted rating. The approach is to consider a low-rank approximation $A^{\prime}$ of the preference matrix $A$, which can be written in terms of the SVD of $A$ as $A^{\prime}=\sum_{\ell=1}^{k}\sigma_{\ell}\bm{u}^{(\ell)}{\bm{v}^{(\ell)}}^{T}$. Here $k$ is a parameter of choice in the algorithm, not the actual rank of the original preference matrix. This low-rank approximation can be equivalently written as $A^{\prime}=A\sum_{\ell=1}^{k}\bm{v}^{(\ell)}{\bm{v}^{(\ell)}}^{T}$ and therefore the $i$-th row of $A^{\prime}$ is given by

[TABLE]

Expressing $A^{\prime}_{i}$ as a column vector, we recognize the target vector as

[TABLE]

where $\lambda_{\ell}=\langle A_{i}^{T},\bm{v}^{(\ell)}\rangle$ as in Eq. (4).

The matrix $A^{\prime}$ is typically no longer sparse, so the previously missing ratings can be predicted based on the non-zero entries of this low-rank approximation. In this sense, the algorithm can be interpreted as an approximate reconstruction of a preference matrix $P$ that contains ratings for all users and movies, of which we are only given a few sample ratings in the form of the sparse matrix $A$. If $P$ is well approximated by a low-rank matrix and $A$ has sufficiently many entries of $P$, then a low-rank approximation of $A$ is also close to $P$.

As in the previous section, we summarize the results of applying the quantum-inspired algorithm for recommendation systems in Tables 6 and 7. We observe from Table 6 that sampling $450$ rows and $4500$ columns from the full preference matrix is enough to describe the first ten singular values with a small relative error of $\eta_{\sigma}=0.06$. On the other hand, the errors of the reconstructed matrices $A$, $A^{+}$, and estimated coefficients show larger values of $0.32$, $0.66$ and $0.58$, respectively. Finally, we find that the entries of the approximated solution vector, computed using Eq. (36) within the low-rank approximation with $k=10$ show typical relative errors of $\eta_{\bm{x}}=0.71$. In Table 7 we compare the running times between the quantum-inspired algorithm and the direct calculation method. As in the portfolio optimization example, the matrices are small enough that computing their SVD can be done in a short time. The runtime of the quantum-inspired algorithm is again dominated by the coefficient estimation step, even when the resulting errors in estimating these coefficients is relatively large.

The numerical examples we have studied indicate that when applied to moderately-sized real-life data sets, because they rely on a more intricate procedure, the quantum-inspired algorithms take more time than exact diagonalization, and because they rely on sampling for coefficient estimation, lead to higher inaccuracies. These results suggest that in order to provide a speedup over preexisting classical algorithms, the quantum-inspired algorithms must be applied to very large data sets where exact diagonalization is impossible and where even the linear scaling of FKV prevents its direct application.

3.6 Practical complexity of quantum-inspired algorithms

In this section, we comment on the complexity of quantum-inspired algorithms based on the numerical tests and theoretical analysis performed thus far. The matrix Chernoff bound of Ref. [10] indicates that the error in the FKV algorithm scales as $\varepsilon=O(1/\sqrt{r})$ with the dimension $r$ of the matrix $C$. This is consistent with the results of empirical tests on Gaussian random matrices as reported in Sec. 3.3. Since the complexity of computing the SVD of an $r\times c$ matrix with $c=O(r)$ scales as $O(r^{3})$, this shows that the complexity of approximate SVD calculations with the FKV algorithm grows as $O(1/\varepsilon^{6})$. This is compatible with the complexity bound from Ref. [10], suggesting that complexity with respect to error may be tight.

It was shown in Ref. [9] that the FKV algorithm can produce an approximate SVD with error $\varepsilon$ by asymptotically setting $r=c=O(\max\{k^{4}/\varepsilon^{2},k^{2}/\varepsilon^{4}\})$. This was improved in Ref. [10], where it was shown that an error $\varepsilon$ could be obtained when $r\neq c$ by choosing $r=O(c)=O(k^{2}/\varepsilon^{2})$. This implies a complexity of $O(k^{6}/\varepsilon^{6})$ when all other parameters are fixed. Importantly, such a choice of $r$ only makes sense if the resulting $r\times c$ matrix is smaller than the original input matrix, i.e., if $r\leq m$ and $c\leq n$. In practice, it is always possible to set $r\leq m$ and $c\leq n$, resulting in a guaranteed reduced runtime in computing the SVD at the expense of an error in the approximation. This is the strategy adopted in our experiments. It is important to understand that, when implemented properly, the runtime of quantum-inspired algorithms is never larger than a direct calculation. The crucial point is whether the resulting errors are sufficiently low, since it becomes prohibitively expensive to reduce them to arbitrarily small values due to the $O(1/\varepsilon^{6})$ scaling.

In regimes where the largest source of errors originates from coefficient estimation, it is possible that the runtime of quantum-inspired algorithms is dominated by the complexity of this step. This occurred in several examples considered in this work, including high-dimensional problems with low-rank matrices. The complexity of coefficient estimation is captured in Eqs. (13) and (15), which have smaller exponents than the complexity bounds of Refs. [10] and [31]. Numerical tests support the theoretical calculations. The linear dependency shown in Fig. 4 between error, rank, and condition number for a fixed number of samples $N$ indicates that $N$ must be a polynomial in $(k\kappa/\varepsilon)$. Our bound states that $N=O(k^{2}\kappa^{2}/\varepsilon^{2})$ in accordance with this behavior. Depending on the properties of the input problem, different complexity regimes are possible, each determined by the steps of the algorithms that dominate the runtime.

4 Conclusion

In terms of asymptotic complexity, quantum-inspired algorithms constitute a breakthrough in our understanding of the boundary between classical and quantum algorithms: they imply that certain linear algebra problems can be performed in sublinear classical time. Our results show that the proven complexity bounds for these algorithms do not actually reflect their practical runtimes. The proven complexity of the linear systems algorithm is $\tilde{O}(\kappa^{16}k^{6}\|A\|^{6}_{F}/\varepsilon^{6})$ [10], while for recommendation systems it is $\tilde{O}(k^{12}/(\varepsilon^{12}\eta^{6}))$ [31], where $\eta$ is an error parameter. In the implementation of the algorithm, we observe a significantly faster runtime than these bounds would suggest. This indicates that care must be taken when employing these bounds to make statements about the performance of the algorithms. Our results are also encouraging for complexity theorists aiming to improve these bounds.

In our analysis, we showed that quantum-inspired algorithms can provide reasonably low errors in relatively short times even in the regime of extremely large-dimensional problems. However, in our implementation, the performance requires matrices of very low rank and condition number; the errors in the outputs grow noticeably when rank and condition number are increased. Compared to previously-known methods such as the Frieze-Kannan-Vempala (FKV) algorithm, quantum-inspired algorithms differ in their use of sampling techniques to estimate coefficients and sample from the solution vectors. On the one hand, a direct calculation, as done in FKV, requires linear time in the dimension of the input matrix, which can be done extremely fast even for problems of large size. On the other hand, sampling methods scale polylogarithmically with dimension, but they incur additional polynomial costs in the rank, condition number, and error in the estimation. Therefore, one conclusion from our work is that despite their asymptotic scaling, sampling techniques for coefficient estimation do not lead to practical improvements compared to the direct computation used in previously-existing classical algorithms: quantum-inspired techniques only become advantageous for problems of extremely large dimension.

Additionally, we have shown that when employing the FKV algorithm for approximate SVD, a scaling of $O(k^{6}/\varepsilon^{6})$ with respect to the error $\varepsilon$ and rank $k$ follows directly from the matrix Chernoff bound and the error bounds from Ref. [10]. Therefore, these exponents are likely to be fundamental and not the result of proof techniques. Overall, our results indicate that quantum-inspired algorithms perform well in practice, but only under the restrictive conditions of large-dimensional input matrices with very low rank and condition number. It remains unclear whether datasets with these properties actually appear in practice.

By contrast, the complexity of quantum algorithms does not depend on matrix rank, and therefore they work properly even for full-rank problems. Matrices originating in practical applications are often sparse and they typically have large effective ranks; these are the problems that can be tackled by quantum algorithms. Dealing with matrices of large condition number remains challenging for all techniques.

Several open questions remain. The numerical tests on random matrices performed in Section 3.3 indicate that when all other parameters are fixed, errors in the approximation of the output vector increase linearly with rank and condition number. It is important to understand whether this is a general feature of the algorithm. Additionally, it is still of interest to determine the empirical scaling of quantum-inspired algorithms with respect to all relevant parameters. Our results indicate that for linear systems, the dependency on rank and error is likely tight, but the findings of section 3.2 show that not all exponents are tight. Finally, it would be of great interest to implement highly-optimized versions of our code on powerful supercomputers capable of operating much larger datasets than those tested in this work.

Acknowledgements

We thank Christian Weedbrook, Nathan Killoran, Nicolás Quesada, and Iordanis Kerenidis for useful discussions. We are grateful to András Gilyén for providing valuable feedback on an early draft of this manuscript. S. Lloyd was funded by AFOSR under a MURI on Optimal Quantum Measurements and State Verification, by IARPA under the QEO program, by ARO, and by NSF.

Appendix

Appendix A Worst-case hardness of coefficient estimation

Here we show that the coefficient estimation cannot generally done in sublinear time. The main strategy is to encode an NP-Hard problem – approximating the partition function of Ising models – into the coefficients to be estimated. The exponential time hypothesis then implies that this coefficients cannot be estimated in sublinear time.

We begin by noting that, for any sum of the form $\mathcal{Z}=\sum_{i=1}^{n}w_{i}$, it is possible to find vectors $\bm{y},\bm{z}$ such that $w_{i}=y_{i}z_{i}$ for all $i=1,2,\ldots,n$. This means that any sum $\mathcal{Z}$ can be expressed as an inner product $\mathcal{Z}=\langle\bm{y},\bm{z}\rangle$. If estimating $\mathcal{Z}$ is NP-Hard, then so is calculating the coefficient $\lambda=\langle\bm{x},\bm{y}\rangle$. As a concrete example, consider the partition function of the Ising model

[TABLE]

where $s=(s_{1},\ldots,s_{m})\in\{-1,1\}^{m}$ is the configuration of $m$ spins, $\beta$ is the inverse temperature, and $E(s)=-\sum_{i<j}J_{ij}s_{i}s_{j}-\mu\sum_{j}h_{j}s_{j}$ is the energy of the configuration. It is known that approximating this partition function is NP-Hard in a general setting [29]. Since there are $n:=2^{m}$ different configurations $s$, each can be labeled by a number $k=1,\ldots,n$, also expressing energies as a function of $k$. As discussed above, we can also define $n$-dimensional vectors $\bm{x},\bm{y}$ such that $\mathcal{Z}_{I}=\sum_{k=1}^{n}x_{k}y_{k}=\langle\bm{x},\bm{y}\rangle$, for instance by setting $x_{k}=y_{k}=e^{-\beta E(k)/2}$. Therefore, there exist vectors such that computing their inner product is as difficult as computing the NP-Hard partition function $\mathcal{Z}_{I}$. The exponential time hypothesis [15] states that the best possible algorithms for computing NP-complete problems require $2^{\delta m}$ time for some $\delta>0$. Thus, unless this hypothesis is false, estimating coefficients requires $2^{\delta m}=O(n)$ time in a worst-case setting.

Appendix B Formulas for high-dimensional problems

We begin by showing that all rows of the high-dimensional matrices considered in this paper have equal norm. Denote by $A_{y}$ the $\mathtt{y}$-th row of $A$. This can be written as

[TABLE]

where $\mathtt{z}^{\prime}=\mathtt{z}\oplus\mathtt{y}$ and we have defined a new basis $\{\tilde{e}^{(\mathtt{z}^{\prime})}\}=\{e^{(\mathtt{z}^{\prime}\oplus\mathtt{y})}\}$. Therefore the norm of the row is

[TABLE]

which does not depend on $\mathtt{y}$, meaning it is the same for all rows. Indeed, the rows of $A$ have the same elements; they simply appear in different positions depending on the corresponding row label $\mathtt{y}$.

The entries of matrix $C$ can be expressed in terms of known quantities as

[TABLE]

where $\bm{a}=(a_{1},a_{2},\ldots,a_{2^{k}})^{T}$ is the vector containing all the possible values of the entries $a_{\mathtt{y},\mathtt{z}}$ and we used $\|R_{\cdot,j_{t}}\|=\sqrt{\frac{2^{n}}{r}}\sqrt{\sum_{s^{\prime}=1}^{r}A_{i_{s^{\prime}},j_{t}}^{2}}$.

Finally, to estimate the coefficient $\lambda_{\ell}=\frac{1}{\sigma_{\ell}^{2}}\langle\bm{v}^{(\ell)},A^{T}\bm{b}\rangle$, we sample the random variable $X_{\mathtt{y},\mathtt{z}}=\frac{\|A\|_{F}^{2}}{A_{\mathtt{y},\mathtt{z}}}b_{\mathtt{y}}v_{\mathtt{z}}^{(\ell)}$ whose average is equal to $\langle\bm{v}^{(\ell)},A^{T}\bm{b}\rangle$. This random variable is given by

[TABLE]

Appendix C Complementary results for random matrices

In this section we explain the methodology to generate the low-rank Gaussian random matrices used in Section 3.3 and show further numerical results to benchmark the accuracy of the algorithms used to approximate the SVD of matrix $A$ and to estimate the coefficients $\lambda$. To generate $U$, we first sample an $m\times k$ Gaussian random matrix $U^{\prime}$ with entries drawn independently from the standard normal distribution $\mathcal{N}(0,1)$. The columns of this matrix are generally not orthogonal, so we perform a QR decomposition $U^{\prime}=QR$, where $Q$ is an $m\times k$ orthogonal matrix and $R$ is a $k\times k$ upper triangular matrix. We then simply set $U:=Q$. An analogous method is used to generate $V$. Finally, given a target condition number $\kappa$, we select the largest singular value $\sigma_{\mathrm{max}}$ uniformly at random in the interval $[1,500]$. This fixes the smallest singular value as $\sigma_{\mathrm{min}}=\sigma_{\mathrm{max}}/\kappa$. Other singular values are sampled from the interval $(\sigma_{\mathrm{min}},\sigma_{\mathrm{max}})$ by using the quadrant law for singular values [28].

In Figs. 5(a)-5(c), we benchmark the approximate SVD of random matrices with the same dimensions as the one in Fig. 3, but with increasing values of the rank $k$. In this case the condition number of these matrices has been fixed to $\kappa=5$. Although the errors in the estimation of singular values does not appear to be strongly influenced by the rank, errors in the approximated singular vector become apparent in the relative errors $\eta_{A}$ and $\eta_{A^{+}}$ for the reconstructed matrices as $k$ increases. In particular, we notice that the error of the pseudo-inverse matrix, which is the one required to compute the solution vector, scales linearly with rank showing relatively large errors larger than $45$ % for $k=100$.

In Figs. 6(a)-6(c) we benchmark the approximate SVD of low-rank random matrices with $k=5$ while increasing their condition number. We observe from Figs 6(a)-6(b) almost no dependence of the errors $\eta_{\sigma}$ and $\eta_{A}$ with condition number. On the contrary, Fig. 6(c) demonstrates that the relative error $\eta_{A^{+}}$ of the reconstructed pseudo-inverse matrix grows linearly with condition number, exhibiting values greater than $100\%$ already for $\kappa>100$. The latter is somehow expected since the error of the approximated singular vectors are propagated to matrix $A^{+}$ roughly as $1/\sigma_{\ell}$, i.e., the smaller the singular values, the larger the error $\eta_{A^{+}}$ due to the approximate SVD of matrix $A$.

In Figs. 7(a)-(c) we investigate the error $\eta_{\lambda}=\sum_{\ell=1}^{k}|\tilde{\lambda}_{\ell}-\lambda_{\ell}|/|\lambda_{\ell}|$ of the estimated coefficients $\lambda_{\ell}$. In all cases, the estimation of the coefficients was performed by taking $N=10^{4}$ samples of the corresponding random variable. In Fig. 7(a) we show that the values of $\eta_{\lambda}$ display a weak dependence on the number of sampled rows and columns. We notice, however, that this is not the case as we increase the rank and condition number of the matrix. For example, large values of $\eta_{\lambda}>100$ % are observed for matrices with rank $k>20$ as shown in Fig. 7(b). Furthermore, Fig. 7(c) shows clearly that even for a matrix with rank as low as $k=5$, when the condition number is increased, the error of the estimated coefficients ramps up rapidly to large values.

Appendix D Singular values and coefficients for portfolio optimization and recommendation systems

In Fig. 8, we show the first ten singular values of the portfolio optimization matrix as well as the exact and approximate coefficients. Similarly, Fig. 9 shows analogous plots for the case of movie recommendations.