Quantum-inspired algorithms in practice

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 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 but only 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.

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 [18]. 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 ε may be tight for the linear systems algorithm. For problems whose input matrices have larger ranks and condition numbers, the algorithms struggle 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.

II. 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. [8,[15][16][17]. 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 x = b, where A ∈ R m×n , x = (x 1 , x 2 , . . . , x n ) T , and b = (b 1 , b 2 , . . . , b m ) T .
We consider problems of the following form: Given an m × n matrix A ∈ R m×n with SVD the goal is to sample entries of the n-dimensional vector with respect to the length-square probability distribution p x (i) = x 2 i / x 2 . The coefficients λ depend on the matrix A and possibly on other inputs. For linear systems of equations, x is the solution vector given by x = A + b, with A + the Moore-Penrose pseudoinverse. As discussed in Ref. [16], the coefficients λ of Eq. (2) are given by the inner product 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 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 where A i is the i-th row of A [8].
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 [13], 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: 1. Approximate SVD: The Frieze-Kannan-Vempala (FKV) algorithm [19] is used to compute approximate singular valuesσ and approximate right singular vectorsṽ ( ) . The singular vectors are not calculated in their entirety; instead, a description is found allowing efficient computation of any given entry v ( ) i .

Coefficient estimation:
Based on the results of step 1, the coefficientsλ are approximated using Monte Carlo estimation techniques.
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 vectorsx = k =1λ ṽ ( ) . 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(poly(k, κ, ε, log n, log m)). By contrast, the strategy originally suggested in Ref. [19] 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 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.

A. Approximate SVD
The Frieze-Kannan-Vempala (FKV) algorithm [19] is a method to compute low-rank approximations of matrices. The FKV algorithm is just a particular example among several randomized methods for computing approximate matrix decompositions [20][21][22][23][24][25], see Ref. [26] for a survey of these techniques. It is an interesting question whether other randomized methods for approximate matrix decompositions can be used to improve quantum-inspired algorithms. For concreteness, in this work we focus on analyzing the FKV algorithm.
The main idea behind FKV is the following: instead of performing a singular value decomposition of the large input matrix A ∈ R m×n , a smaller matrix C ∈ R r×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 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: For each row i, a length-square distribution over its entries is also computed, To sample from these distributions, it is possible to employ a specialized data structure as suggested in Refs. [8,13,19] that allows sampling in time that is logarithmic in dimension. The algorithm proceeds as follows: 1. Sample r row indices i 1 , i 2 , . . . , i r from the row distribution p(i). For each row index, select the row A is and renormalize it as This defines an r × n matrix R consisting of the rescaled rows R is .
2. Select an index s ∈ {1, 2, . . . , r} uniformly at random, then sample a column index j from the column distribution q is (j). Repeat this a total of c times to obtain column indices j 1 , j 2 , . . . , j c . For each column index, select the column R ·,jt and renormalize it as C ·,jt = A F √ c R·,j t R ·,jt . This defines a matrix C consisting of the rescaled columns C ·,jt .
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σ the singular values of C and by ω ( ) its left singular vectors. The FKV FIG. 1: Schematic representation of the FKV algorithm. The input is an m × n matrix A, which in this example has ten rows and columns. After performing one pass through the matrix, it is possible to perform length-square sampling of rows and columns. A total of r rows are sampled from this length-square distribution, then renormalized and stacked together to construct a new r × n matrix R. In this example, we set r = 4 and obtain the row indices (i1, i2, i3, i4) = (2,8,3,6). Once the matrix R has been built, c column indices are drawn by repeatedly choosing rows of R uniformly at random, then sampling from the corresponding length-square distribution of column indices. In this example, we set c = 4 and obtain the column indices (j1, j2, j3, j4) = (3,8,1,3). Note that repeated indices are allowed in the algorithm. The final matrix C is built by renormalizing these columns of R and grouping them together. The singular values of C approximate those of A, while the singular vectors of C can be used to reconstruct approximations to the singular vectors of A.
algorithm directly provides approximations of singular values,σ ≈ σ for = 1, 2, . . . , k. Approximate right singular vectorsṽ ( ) of A are obtained using the expressioñ while approximate left singular vectorsũ ( ) are given bỹ 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 m and c 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 × n matrix A requires O(min{m 2 n, mn 2 }) time using naive matrix multiplication. From that point onwards, computing coefficients λ and the solution vector x takes only linear time in m, n. In FKV, we instead calculate the SVD of an r × c matrix, which requires O(min{r 2 c, rc 2 }) time. When r, c are much smaller than m, n, this 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.

B. Coefficient estimation
The coefficients λ 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 λ is given by for some appropriate vectors y, 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 y. Define the random variable χ that takes values χ i = y i z i /p x (i), where the indices i are sampled from the length-square distribution p y (i) = y 2 i / y 2 . The expectation value of the random variable satisfies Similarly, the second moment is and σ 2 χ = E(χ 2 ) − E(χ) 2 is the variance of χ. The strategy to estimate coefficients λ is to draw N samples χ (1) , χ (2) , . . . , χ (N ) from χ and compute the unbiased estimatorλ = 1 N N j=1 χ (j) ≈ λ. This constitutes a form of importance sampling in Monte Carlo estimation: large entries of 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: := Var(λ)/|λ|. The variance of the estimator is Var(λ) = σ 2 χ /N , leading to a precision = σ χ /(|λ| √ N ) in the estimator. This implies that the number of samples N needed to, with high probability, achieve a precision is where θ is the angle between the vectors y and z.

C. Complexity of coefficient estimation
Quantum-inspired algorithms differ from the FKV algorithm because they estimate the coefficients λ 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 E 1 cos 2 θ = n, so we can anticipate that the number of samples needed may scale linearly with dimension. If 1 cos 2 θ = 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 λ are given by In the quantum-inspired algorithm, we do not have direct access to the length-square distributions of either A T b or v ( ) . Instead, the strategy in Ref. [16] is to express the coefficients as λ = 1 σ 2 Tr(A T b v ( ) T ) and sample the entries A ij of A. In this case the the number of samples needed in the estimation is [16]: where we have defined σ : and we have used Eq. (12) as well as the identity A 2 F = σ 2 . In the special cases of recommendation systems and portfolio optimization, the coefficients are proportional to an inner product The number of samples then satisfies Define the k-dimensional vectors σ := (σ 1 , σ 2 , . . . , σ k ) T , µ := (u i ) T and their Hadamard product 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 -th entry of that vector. We now determine the worst-case scaling for these ratios. We consider the ratio σ 2 /σ 2 for concreteness; the same will apply for b and ν. The largest ratio occurs for the smallest singular value σ min and we have Similarly, this ratio is largest when σ max is as small as possible. This occurs when the largest k − 1 singular values are equal to each other. In that case, defining σ := σ / σ such that σ = 1 we have and therefore Defining κ β := β max /β min and κ ν := ν max /ν min , where β max/min and ν max/min are the largest and smallest entries of their corresponding vectors, we conclude that the number of samples required to estimate coefficients in quantuminspired algorithms for linear systems scales as and for recommendation systems as 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 = κ = κ β = 100 means that an order of N = 10 16 samples are needed to estimate the coefficients with an error of = 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.

D. Sampling solution vectors
The approximate SVD and coefficient estimation steps respectively provide approximate singular valuesσ and approximate coefficientsλ . The approximate right singular vectors of A are not constructed explicitly, but are implicitly defined by Eq. (5), namelyṽ ( ) = 1 σ R T w ( ) . This allows an implicit calculation of the approximate The challenge is to sample from these vectors using only query access to the entries of w and R. Quantum-inspired algorithms achieve this using a well-known technique known as rejection sampling. The steps are as follows [8]: 1. Sample a row index i uniformly at random.
2. Sample a column index j from the length-square distribution q i (j) = |R ij | 2 / R i 2 .
3. Output j with probability | w,R·,j | 2 R·,j 2 w 2 , or sample j again. The expected number of repetitions before outputting a column index j is r w 2 /( n j=1 | w, R ·,j | 2 ) = r w 2 / x 2 . The steps of rejection sampling can be performed quickly. The overhead arises from calculating the inner product w, R ·,j , 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.

III. 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. We examine the performance of each step of the algorithm individually, revealing the ones that contribute most to the overall runtime and which ones lead to the most significant source of errors. 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.

A. 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 where x, y, z are n-bit strings, the vectors e (y) form a basis of the 2 n -dimensional vector space, and the x ( ) are distinct n-bit strings for = 1, 2, . . . , k. We define the input matrix as where σ 1 , σ 2 , . . . , σ k are the singular values and we have implicitly defined a y,z = k =1 σ (−1) x ( ) ·(y⊕z) . The right and left singular vectors of A are the same, given by . 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 y, we select a column index z uniformly at random and accept it with probability The optimization is trivial since max z |a y,z | 2 = k =1 (σ ) 2 . Once row indices i 1 , . . . , i r and column indices j 1 , . . . , j c have been sampled, the entries of the r × c matrix C in FKV can be expressed as where a = (a 1 , a 2 , . . . , a 2 k ) T is a vector containing all the possible values of the entries a y,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 and use its sample mean as an estimator for the coefficients. Finally, the approximate solution vectorx is given bỹ Any entry of the approximate solution vector can be computed efficiently from the entries of the approximate right singular vectorsṽ ( ) , 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 ≈ 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 x ( ) that define each singular vector are chosen uniformly at random. We study three examples with k = κ = κ β = 3, 5, 10, employing N = 10 4 samples in the estimation of the coefficients. We report the average and standard deviation of the errors |λ −λ | |λ | 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 |xi−xi| |xi| . The results are summarized in Table I.  These results show that for problems with the appropriate properties -namely very low rank and condition numberquantum-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 O(κ 16 k 6 A 6 F /ε 6 ). Even setting A F = 1 and ε = 1, for k = κ = κ β = 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.

B. Random matrices
We study the quantum-inspired algorithm for linear systems of equations Ax = b for randomly chosen A and b.
To generate a random matrix A of dimension m × n, rank k, and condition number κ, instead of creating A directly, we build the components of its singular value decomposition in matrix form, A = U ΣV . Here, U is an m × k matrix whose columns are the left singular vectors of A. Similarly, V is a k × n matrix whose rows are the right singular vectors, and Σ is a k × 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 Σ are explained in Appendix C. The vector b has been chosen to have the form b = k =1 β u ( ) , with coefficients β drawn from the standard normal distribution 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 store the entire distribution in memory and use fast built-in numerical sampling algorithms in Python. In Figs. 2(a)-2(c) we monitor, respectively, the mean relative error η σ = 1 in estimating the singular values as well as the errors η A = Ã − A F / A F and η + 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/ √ r for these random matrices, which is in agreement with the matrix Chernhoff bound appearing in Theorem 3 of Ref. [16]. We also observe that the reconstructed matrix A has smaller errors compared to the approximated pseudo-inverse matrix A + . This is because in reconstructing In all cases, error bars denote the standard deviation for 10 repetitions of the algorithm. The main results are shown in Fig. 3, where we plot the error η x of the approximated solution vectorx calculated as the median of the relative errors |x i − x i |/|x i | with respect to all n entries of vectorx. First, we investigate in Fig. 3(a) how the error η 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 λ defined in Eq. (12). As expected, η 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. 3(b) we set r = c = 4250 and κ = 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 η 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. 3(c) for the solution vector error as a function of the condition number κ. Small errors of the order of 10% are only obtained for small condition numbers, while large values of η x ≥ 100% already appear for κ ≥ 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.
Finally, Tables II and III 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 κ = 5. In Table II, we show the average relative errors of singular values η σ and coefficients η λ as defined in the previous section. We also report the relative errors η A = A − A F / A F and η A + = A + − A + F / A + F of the reconstructed matrices A and A + . Sampling 4250 rows and columns from the original 40000 × 20000 matrix A is enough to have a good approximation full matrix. The largest error in Table II corresponds to the estimated coefficients. Table III 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.   Table II. The parameter tLS is the time required to construct length-squared (LS) probability distributions over rows and columns. t C SVD accounts for the time spent sampling rows and columns from matrix A and building and decomposing the sampled matrix C. t λ is the running time of the estimation of the coefficients λ. Analogous quantities are defined for the direct calculation method. tx is the elapsed time to sample 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 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 [27], 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 r j = (r 1,j , r 2,j , . . . , r n,j ) T , where r i,j is the return of asset i on the j-th day. The vector of expected returns r and the correlation matrix Σ are defined as The vector r captures the expected returns for each asset, while the correlation matrix Σ 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 Minimize: subject to: where w is a portfolio allocation vector that determines how much wealth is invested in each asset. We allow w to have negative entries, which is interpreted as an indication to short-sell the corresponding asset. As discussed in Ref. [28], this problem can be equivalently cast as a linear system of equations: where ν is a Lagrange multiplier. In this case, the linear system Ax = b is given by 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 The solution vector x is then where we have identified the coefficients λ = v ( ) , A T b /σ 2 as in Eq. (3). The vector b has only one non-zero entry, so it holds that 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 [29]. 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 × 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 σ max = 40.2, σ min = 1.8 × 10 −4 , and its condition number is κ = 2.23 × 10 4 . More details can be found in Appendix D. We set a target return µ equal to the average return over all 474 stocks in the index.
In Tables IV 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 λ . The relative errors characterizing the approximate solution vector are considerable, with values of η x = 0.74. In this example, we have calculated the approximate solution vector by using Eq. (34) within a low-rank approximation with k = 10. Notice from Table IV that the large errors η + A in the reconstructed pseudo-inverse matrix and estimated coefficients η λ translate into large relative errors of the solution vector.
Finally, in Table V 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. IV: Relative errors with corresponding standard deviations computed over 10 repetitions of the quantum-inspired algorithm applied to a portfolio optimization problem. Approximated singular vales and vectors were obtained by decomposing matrix C with dimension r = 340 and c = 340. Coefficients λ have been estimated by performing N = 10 4 samples. Here, ησ is the error in approximating singular valuesσ, while ηA and η A + are the errors of the reconstructed matrices A and A + , respectively. η λ is the error in estimating coefficients, and ηx is the the median of the relative errors |xi − xi|/|xi| with respect to all n = 475 entries of the approximated vectorx.

Parameters Error
Case study r c N ησ ηA η + A η λ ηx S&P 500 stock index 340 340 10 4 0.08 ± 0.02 0.16 ± 0.03 1.13 ± 0.12 1.58 ± 0.79 0.74 ± 0.19 TABLE V: 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 IV. The parameter tLS is the time required to construct lengthsquared (LS) probability distributions over rows and columns. t C SVD accounts for the time spent sampling rows and columns from matrix A and building and decomposing the sampled matrix C. t λ is the running time of the estimation of the coefficients λ. Analogous quantities are defined for the direct calculation method. tx is the elapsed time to sample 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 We analyze the performance of the algorithm for recommendation systems on the MovieLens 100K database [18], 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 σ max = 534.4, σ min = 2.95, and its condition number is κ = 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 of the preference matrix A, which can be written in terms of the SVD of A as A = k =1 σ u ( ) v ( ) 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 and therefore the i-th row of A is given by Expressing A i as a column vector, we recognize the target vector as where λ = A T i , v ( ) as in Eq. (4).
The matrix A 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, 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 VI and VII. We observe from Table VI 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 η σ = 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 η x = 0.71. In Table VII 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.    Table VI. The parameter tLS is the time required to construct length-squared (LS) probability distributions over rows and columns. t C SVD accounts for the time spent sampling rows and columns from matrix A and building and decomposing the sampled matrix C. t λ is the running time of the estimation of the coefficients λ. Analogous quantities are defined for the direct calculation method. tx is the elapsed time to sample 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 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 Chernhoff bound of Ref. [16] indicates that the error in the FKV algorithm scales as ε = O(1/ √ 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. III B. Since the complexity of computing the SVD of an r × 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/ε 6 ). This is compatible with the complexity bound from Ref. [16], suggesting that complexity with respect to error may be tight.
It was shown in Ref. [19] that the FKV algorithm can produce an approximate SVD with error ε by asymptotically setting r = c = O(max{k 4 /ε 2 , k 2 /ε 4 }). This was improved in Ref. [16], where it was shown that an error ε could be obtained when r = c by choosing r = O(c) = O(k 2 /ε 2 ). This implies a complexity of O(k 6 /ε 6 ) when all other parameters are fixed. Importantly, such a choice of r only makes sense if the resulting r × c matrix is smaller than the original input matrix, i.e., if r ≤ m and c ≤ n. In practice, it is always possible to set r ≤ m and c ≤ 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/ε 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. [16] and [8]. Numerical tests support the theoretical calculations. The linear dependency shown in Fig. 3 between error, rank, and condition number for a fixed number of samples N indicates that N must be a polynomial in (kκ/ε). Our bound states that N = O(k 2 κ 2 /ε 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.

IV. 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Õ(κ 16 k 6 A 6 F /ε 6 ) [16], while for recommendation systems it isÕ(k 12 /(ε 12 η 6 )) [8], where η is a failure probability. 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, such 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. A direct calculation, as done in FKV, requires linear time in the dimension of the input matrix. Sampling methods scale polylogarithmically with dimension, but they incur additional polynomial costs in the rank, condition number, and error in the estimation. We performed a precise analysis of the complexity of sampling methods for coefficient estimation, revealing that typically a very large number of samples is required to obtain low errors in the estimation. This implies that quantum-inspired techniques only become advantageous for problems of extremely large dimension. 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. Our results suggest that breakthroughs in computational complexity theory such as the dequantization approach that led to the quantum-inspired algorithms of Refs. [8,[15][16][17] are unlikely to alter the set of problems where quantum algorithms may have a practical impact.
where a = (a 1 , a 2 , . . . , a 2 k ) T is the vector containing all the possible values of the entries a y,z and we used Finally, to estimate the coefficient λ = 1 This random variable is given by

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 III B 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 λ. To generate U , we first sample an m × k Gaussian random matrix U with entries drawn independently from the standard normal distribution N (0, 1). The columns of this matrix are generally not orthogonal, so we perform a QR decomposition U = QR, where Q is an m × k orthogonal matrix and R is a k × k upper triangular matrix. We then simply set U := Q. An analogous method is used to generate V . Finally, given a target condition number κ, we select the largest singular value σ max uniformly at random in the interval [1,500]. This fixes the smallest singular value as σ min = σ max /κ. Other singular values are sampled from the interval (σ min , σ max ) by using the quadrant law for singular values [32].
In Figs. 4(a)-4(c), we benchmark the approximate SVD of random matrices with the same dimensions as the one in Fig. 2, but with increasing values of the rank k. In this case the condition number of these matrices has been fixed to κ = 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 η A and η 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. 5(a)-5(c) we benchmark the approximate SVD of low-rank random matrices with k = 5 while increasing their condition number. We observe from Figs 5(a)-5(b) almost no dependence of the errors η σ and η A with condition number. On the contrary, Fig. 5(c) demonstrates that the relative error η A + of the reconstructed pseudo-inverse matrix grows linearly with condition number, exhibiting values greater than 100% already for κ > 100. The latter is somehow expected since the error of the approximated singular vectors are propagated to matrix A + roughly as 1/σ , i.e., the smaller the singular values, the larger the error η A + due to the approximate SVD of matrix A. In Figs. 6(a)-(c) we investigate the error η λ = k =1 |λ − λ |/|λ | of the estimated coefficients λ . In all cases, the estimation of the coefficients was performed by taking N = 10 4 samples of the corresponding random variable. In Fig. 6(a) we show that the values of η λ 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 η λ > 100 % are observed for matrices with rank k > 20 as shown in Fig. 6(b). Furthermore, Fig. 6(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.