--- title: "A Quantum-inspired Classical Algorithm for Separable Non-negative Matrix Factorization" authors: - "Zhihuai Chen" - "Yinan Li" - "Xiaoming Sun" - "Pei Yuan" - "Jialin Zhang" doi: "10.24963/ijcai.2019/627" url: "https://ijcai.org/proceedings/2019/627" pdf: "https://ijcai.org/proceedings/2019/0627.pdf" venue: "Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence" year: 2019 pages: "4511-4517" track: "Main track" keywords: - "Machine Learning Applications: Big data ; Scalability" - "Machine Learning Applications: Applications of Unsupervised Learning" --- # A Quantum-inspired Classical Algorithm for Separable Non-negative Matrix Factorization Zhihuai Chen, Yinan Li, Xiaoming Sun, Pei Yuan, Jialin Zhang *Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence*, pages 4511-4517. DOI: [10.24963/ijcai.2019/627](https://doi.org/10.24963/ijcai.2019/627) ## Abstract Non-negative Matrix Factorization (NMF) asks to decompose a (entry-wise) non-negative matrix into the product of two smaller-sized nonnegative matrices, which has been shown intractable in general. In order to overcome this issue, separability assumption is introduced which assumes all data points are in a conical hull. This assumption makes NMF tractable and widely used in text analysis and image processing, but still impractical for huge-scale datasets. In this paper, inspired by recent development on dequantizing techniques, we propose a new classical algorithm for separable NMF problem. Our new algorithm runs in polynomial time in the rank and logarithmic in the size of input matrices, which achieves an exponential speedup in the low-rank setting. ## Full text ## 1 Introduction Non-negative Matrix Factorization (NMF) aims to approximate a non-negative data matrix A 2 R m�n �0 by the product of two non-negative low rank factors, i.e., A � WH T , where W 2 R m�k �0 is called basis matrix, H 2 R n�k �0 is called encoding matrix and k � minfm; ng. In many applications, an NMF often results in more natural and interpretable part-based decomposition of data [Lee and Seung, 1999]. Therefore, NMF has been widely used in a number of practical applications, such as topic modeling in text, signal separation, social network, collaborative filtering, dimension reduction, sparse coding, feature selection and hyperspectral image analysis. Since computing an NMF is NP-hard [Vavasis, 2009], a series of heuristic algorithms have been proposed [Lee and Seung, 2001; Lin, 2007; Hsieh and Dhillon, 2011; Kim and Park, 2008; Ding et al., 2010; Guan et al., 2012]. All of the heuristic algorithms aim to minimize the reconstruction error, the formula which is a non-convex program and lack optimality guarantee: A natural assumption on the data called separability assumption, was observed in [Donoho and Stodden, 2004] . From a geometry perspective, the separable assumption means that all rows of A reside in a cone generated by a rather smaller number of rows. In particular, these generators are called anchors of A. To solve the Separable Non-Negative Matrix Factorizations (SNMF), it is su�cient to identify the anchors in the input matrices, which can be solved in polynomial time [Arora et al., 2012a; Arora et al., 2012b; Gillis and Vavasis, 2014; Esser et al., 2012; Elhamifar et al., 2012; Zhou et al., 2013; Zhou et al., 2014]. Separability assumption is favored by various practical applications. For example, in the unmixing task in hyperspectral imaging, separability implies the existence of 'pure' pixel [Gillis and Vavasis, 2014]. And in the topic detection task, it also means some words are associated with unique topic [Hofmann, 2017]. In huge datasets, it is useful to pick up some representative data points to stand for other points. Such 'self-expression' assumption helps to improve the data analysis procedure [Mahoney and Drineas, 2009; Elhamifar and Vidal, 2009]. ## 1.1 RelatedWork It is natural to assume all the rows of the input A has unit '1-norm, since '1-normalization translates the conical hull to convex hull while keeping the anchors unchanged. From this perspective, most algorithms essentially identify the extreme points in the convex hull of the ('1-normalized) data vectors. In [Arora et al., 2012a], the authors use m linear programs in O(m) variables to identify the anchors out of m data points, and it is therefore not suitable for dealing with large-scale realworld problems. Furthermore, [Recht et al., 2012] presents a single LP in n 2 variables for SNMF to deal with large-scale problems (but is still impractical for huge-scale problems). There is another class of algorithms based on greedy algorithms. The main idea is to opt a data point on the direction where the current residual decreases fast. The algorithms terminate with a su�ciently small error or a large iteration times. For example, Successful Projection Algorithm (SPA) [Gillis and Vavasis, 2014] derives from Gram-Schmidt orthogonalization with row or column pivoting. XRAY [Kumar et al., 2013] detects a new anchor referring to the residual of exterior data points and updates the residual matrix by solving a nonnegative least square regression. Both of these two algorithms based on greedy pursuit have smaller time complexity compared with LP-based methods. However, the time complexity is still too large for large-scaled data. [Zhou et al., 2013; Zhou et al., 2014] utilize a Divide-andConquer Anchoring (DCA) framework to tackle the SNMF. Namely, by projecting the data set into several low-dimension subspaces, and each projection can determines a small set of anchors. Moreover, it can be proven that all the k anchors can be identified by O(k log k) projections. Recently, a quantum algorithm for SNMF called Quantum Divide-and-Conquer Anchoring algorithm (QDCA), has been presented [Du et al., 2018], which uses the quantum technology to speed up the random projection step in [Zhou et al., 2013]. QDCA implements matrix-vector product (i.e., random projection) via quantum principal component analysis and then a quantum state encoding the projected data points could be prepared e�ciently. Moreover, there are also several papers utilizing dequantizing techniques to solve some low-rank matrix operations, such as recommendation systems [Tang, 2018] and matrix inversion [Gily´en et al., 2018; Chia et al., 2018]. Dequantizing techniques in those algorithms involve two technologies, the Monte-Carlo singular value decomposition and rejection sampling, which could e�ciently simulate some special operations on low-rank matrices. Inspired by QDCA and the dequantizing techniques , we propose a classical randomized algorithm which speeds up the random projection step in [Zhou et al., 2013] and thereby identifies all anchors e�ciently. Our algorithm takes time polynomial in rank k, condition number � and logarithm of the size of matrix. When rank k = O(log(mn)), our algorithm achieves exponentially speedup than any other classical algorithms for SNMF. ## 2 Preliminaries ## 2.1 Notations Let [n] := f1; 2; : : : ; ng. Let spanfxi 2 R n ji 2 [k]g := f P k i=1 �ixij�i 2 R ; i 2 [k]g denote the space spanned by xi for i 2 [k]. For a matrix A 2 R m�n , A(i) and A ( j) denote the ith row and the jth column of A for i 2 [m]; j 2 [n], respectively. Let AR = [A T (i1) ; A T (i2) ; : : : ; A T (ir) ] T where A 2 R m�n and R = fi1; i2; : : : ; irg � [m] (without loss of generality, assume i1 � i2 � � � � � ir). kAkF and kAk2 refer to Frobenius norm and spectral norm, respectively. For a vector v 2 R n , kvk denotes its '2-norm. For two probability distributions p; q (as density functions) over a discrete universe D, the total variation distance between them is defined as � � �p; q � � � TV := 1 2 P i2D jp(i)�q(i)j. �(A) := �max=�min denotes the condition number of A, where �max and �min are the maximal and minimal non-zero singular values of A. ## 2.2 Sample Model In query model, algorithms for SNMF problem require time which is at least linear in the number of nonzero elements of the matrix, since in the worst case, they have to read out all entries. However, we expect our algorithm to be e�cient even if the datasets are extremely large. Considering the QDCA in [Du et al., 2018], one of its advantage is that data is prepared in quantum state and can be access via 'quantum' way (like sampling). Thus, in quantum algorithm, quantum state is served to represent data implicitly which can be read out by measurement only. In order to avoiding reading the whole matrix, we introduce a new sample model other than the query model based on the idea of quantum state preparation assumption. Figure 1: Binary search tree for v = (v1; v2; v3; v4) T 2 R 4 . The leaf node stores v 2 i and interior node stores the sum of the children. In order to restore the original vector, we also store the sign of vi in leaf node. To sampling from Dv, we can start from top and randomly recurring on a child, with probability proportional to its weight. ![Image](https://ijcai.org/proceedings/2019/fulltext/627/image_000000_7677fcc09e8766d0c36d3ff119acecb961745518bdf7a00f07f1ca5c2a71d778.png) Definition 2.1 ('2-norm Sampling). Let Dv denote the distribution over [n] with density function Dv(i) = v 2 i =kvk 2 for v 2 R n . A sample from a distribution Dv is called a sample from v. Lemma 2.2 (Vector Sample Model). There is a data structure storing vector v 2 R n in O(n logn) space, and supporting following operations: - � Querying and updating a entry in O(log n) time; - � Sampling fromDv in O(logn) time; - � Findingkvk in O(1) time. Such a data structure can be easily implemented via Binary Search Tree (BST) (see Figure 1). Proposition 2.3 (Matrix Sample Model). Considering matrix A 2 R m�n , let ˜A and ˜A 0 be the vector whose entry is � � �A(i) � � � and � � �A ( j) � � �, respectively. There is a data structure storing matrix A 2 R m�n in O(mn) space and supporting following operations: - � Querying and updating an entry in O(logm+ log n) time; - � Sampling from A(i) for any i 2 [m] in time O(logn); - � Sampling from A ( j) for any j 2 [n] in time O(logm); - � FindingkAkF, � � �A(i) � � � and � � �A ( j) � � � in time O(1); - � Sampling ˜A and ˜A 0 in time O(logm) and O(logn), respectively. This data structure can be easily implemented via Lemma 2.2, we can just use two arrays of BST to store all rows and columns of A and use two extra BSTs store ˜A and ˜A 0 . ## 2.3 Low-rank Approximations in Sample Model FKV algorithm is a Monte-Carlo algorithm [Frieze et al., 2004] that returns approximate singular vectors of given matrix A in matrix sample model. The low-rank approximation of A can be reconstructed by approximate singular vectors. The query and sample complexity of FKV algorithm are independent of size of A. FKV algorithm outputs a short 'description' of ˆV, which is approximate to a right singular vectors V of matrix A. Similarly, FKV algorithm can output a description of approximate left singular vectors ˆU of A by inputting A T . Let FKV (A; k; �; �) denote the FKV algorithm, where A is a matrix given by sample model, k is the rank of approximate matrix of A, � is error parameter, and � is the failure probability. The FKV algorithm is described in Theorem 2.4. Theorem 2.4 (Low-rank Approximations, [Frieze et al., 2004]). Given matrix A 2 R m�n in matrix sample model, k 2 N and �; � 2 (0; 1), FKV algorithm outputs the description of the approximate right singular vectors ˆV 2 R n�k in O(poly(k; 1=�; log 1 � )) samples and queries of A with probability 1 � �, which satisfies Especially, if A is a matrix with rank k exactly, Theorem 2.4 also implies an inequality: � � �A ˆV ˆV T � A � � � F � p �kAkF . Description of ˆV. Note that FKV algorithm does not output the approximate right singular vectors ˆV directly since their lengths are linear of n. It returns a description of ˆV, which consists of three components: the row index sets T := fit 2 [m]jt 2 [p]g, a vector set U := fu ( j) 2 R p j j 2 [k]g which are singular vectors of a submatrix sampled from A , and its corresponding singular values � := f� ( j) j j 2 [k]g, where p = O(poly(k; 1 � )). In fact, ˆV (i) := ATu (i) =�i for i 2 [k]. Given a description of ˆV, we can sample from ˆV (i) in time O(poly(k; 1 � )) for i 2 [k] [Tang, 2018] and query its entry in time O(poly(k; 1 )). � Definition 2.5 (�-orthonormal). Given � > 0, ˆV 2 R n�k is called �-approximately orthonormal if 1 � �=k � � � � ˆV (i) � � � 2 � 1 + �=k for i 2 [k] and j ˆV (s) ˆV (t) j � �=k for s , t 2 [k]. The next lemma presents some properties of �-approximate orthonormal vectors. Lemma 2.6 (Properties of �-orthonormal Vectors, [Tang, 2018]). Given a set of k �-approximately orthonormal vectors ˆV 2 R n�k , then there exists a set of k orthonormal vectors V 2 R n�k spanning the columns of ˆV such that where �ˆV := VV T represents the orthonormal projector to image of ˆV and c1; c2 > 0 are constants. Lemma 2.7 ([Frieze et al., 2004]). The output vectors ˆV 2 R n�k of FKV (A; k; �; �) is �k=16-approximate orthonormal. ## 3 Fast Anchors Seeking Algorithm In this section, we present a randomized algorithm for SNMF which is called Fast Anchors Seeking (FAS) Algorithm. Especially, the input A 2 R m�n �0 of FAS is given by matrix sample model which is realized via a data structure described in Section 2. FAS returns the indices of anchors in time polynomial logarithmic to the size of matrix. ## 3.1 Description of Algorithm Recall that SNMF aims to factorize A = FAR where R is the index set of anchors. In this paper, an additional constraint is added: the sum of entries in any row of F is 1. Namely, any data point of A resides in convex hull which is the set of all convex combination of AR. In fact, normalizing each row of matrix A by '1-norm is valid, since the anchors remain unchanged. Moreover, Instead of storing '1-normalized matrix A, we can just maintain the '1-norms for all rows and columns. The Quantum Divide-and-Conquer Anchoring (QDCA) is a quantum algorithm for SNMF which achieves exponential speedup than any classical algorithms [Du et al., 2018]. After projecting any convex hull into an 1-dimensional space, the geometric information is partially preserved. Especially, the anchors in 1-dimensional projected subspace are still anchors in the original space. The main idea of QDCA is quantizing random projection step in DCA. It decomposes SNMF into several subproblems: projecting A onto a set of random unit vectors f�i 2 R n g s i=1 with s = O(k log k), i.e., computing A�i 2 R m . Such a matrix-vector product can be e�ciently implemented by Quantum Principle Component Analysis (QPCA). And then it returns a logm-qubits quantum state whose amplitudes are proportional to entries of A�i. Measurement of quantum state outcomes an index j 2 [m] which obeys distribution DA�i . Thus, we can prepare O(poly log m) copies of quantum states, measure each of them in computational basis and record the most frequent index. By repeating procedure above with s = O(k log k) times, we could successively identify all anchors with high probability. As discussed above, the core and most costly procedure is to simulate DA�i . At the first sight, traditional algorithms can not achieve exponential speedup on account of limits of computational model. In QDCA, vectors are encoded into quantum states and we can sample the entries with probability proportional to their magnitudes by measurements. This quantum state preparation overcomes the bottleneck of traditional computational model. Based on divided-and-conquer scheme and sample model (See Section 2.2), we present Fast Anchors Seeking (FAS) Algorithm inspired by QDCA. Designing FAS is quite hard and non-trivial although FAS and QDCA have the same scheme. Indeed, we can simulate DA�i directly by rejection sampling technology. However, the number of iterations of rejection sampling is unbounded. To overcome this di�culty, we translate matrix A into its approximation ˆU ˆU T A, where the columns ˆU 2 R m�k consists of k approximate left singular vectors of matrix A and k = rank(A). Next, it is obvious that y = ˆU T A�i 2 R k is a short vector and we can estimate its entries one by one (see Lemma 3.5) e�ciently. Now the problem becomes to simulate D ˆUy and it can be done by Lemma 3.4 . Given an error parameter �=2, the method described above will result in � � �A�i � ˆU ˆU T A�i � � � < �kAkF � � ��i � � � =2 via Theorem 2.4, which implies � � � �DA�i �D ˆU ˆU T A�i � � � � TV � �kAkF � � ��i � � � = � � �A�i � � �. Namely, the method above introduces an unbounded error in form �kAkF � � ��i � � � = � � �A�i � � � if �i is arbitrary vector in entire space R n . Fortunately, this issue can be solved by generating random vectors f�ig s i=1 lying in row space of A instead of those lying in entire space R n . To generate uniform random unit vectors on the row space of A, we need to find a basis of row space of A. If V 2 R n�k is a set of orthonormal basis of the row space of A (the space spanned by the right singular vectors), and xi is uniform random unit vector on S k�1 , then � = Vxi is a unit random vector in row space of A. Moreover, FKV algorithm will figure out approximate singular vectors ˆV for V, that can Figure 2: An illustration for how to approximate distribution DA�i . Oi represents the final distribution which approximate DA�i . frepresents 'approximate' and represents 'equal to' in a sense of total variation distance. To prove the upper bound for kDA�i ; OikTV , we introduce several medium distributions DA ˆVxi , D ˆU ˆUT A ˆVxi and D ˆUy. From left to right, (1) use a set of right singular vectors V 2 R n�k to generate the random unit vector lying in the row space of A; (2) however, since V cannot be gained e�ciently, use FKV algorithm to figure out an approximation ˆV given by a 'short desription'; (Lemma 2.6, Lemma 3.2) (3) translate A ˆVxi into ˆU ˆU T A ˆVxi since A � ˆU ˆ U T A (Theorem 2.4), where ˆU is the approximate left singular vectors generated by FKV given by a 'short desription'. (4) return an estimation ˜M of M = ˆU T A ˆV 2 R k�k by estimating each entry by Lemma 3.5 and then approximate Mxi (denoted as ˜yi); (5) finally, the rejection sampling works for ˆU˜yi by Lemma 3.4 since ˆU is approximately orthonormal. z ![Image](https://ijcai.org/proceedings/2019/fulltext/627/image_000001_46db4a3fa25ddb980c5109f1a69bbf7819a66eaeaa27f0f0a3c72d77a9f82804.png) x Figure 3: An illustration of finding an anchor of A with rank k = 3. The 3-dimensional space represents the row space of A and A's data points (blue and black points) lie on the '1-normalized plane P. The blue points also stand for the anchors of A. A random vector � is picked up from row space of A and then data points are projected on �. The anchors of the projected space on � are still the anchors of A, such that implies that the red point with the maximum absolute projection component on � is an anchor of A. help us make an approximate ˆ �i = ˆVxi for �i. Therefore, we will estimate distribution D ˆU ˆU T A ˆVxi instead of DA�i . Based on Corollary 3.5, ˆU T A ˆV can be estimated e�ciently. According to Lemma 3.4, ˆUy can be sampled e�ciently, thus we can treat ˜y as estimation of ˆU T A ˆVxi (see Figure 2). Once we can simulate distribution DAVxi , we can figure out the index of the largest component of vector AVxi by picking up O(poly logm) samples (Theorem 3.6). Moreover, according to [Zhou et al., 2013], by repeating this procedure with O(k log k) times, we can find all anchors of A with high probability (For single step of random projection, see Figure 3). ## 3.2 Analysis Now, we propose our main theorem and analyze the correctness and complexity of our algorithm FAS. Theorem 3.1 (Main Result). Given separable non-negative matrix A 2 R m�n �0 in matrix sample model, the rank k, condition number � and a constant � 2 (0; 1), Algorithm 1 returns the indices of anchors with probability at least 1 � � in time ## Correctness In this subsection, we will analyze the correctness of Algorithm 1. Firstly, we show that the columns of V defined in Lemma 2.6 form a basis of row space of matrix A, which is ## Algorithm 1 Fast Anchors Seeking Algorithm Input: Separable non-negative matrix A 2 R m�n �0 in matrix sample model, k = rank(A), condition number �, a constant � 2 (0; 1) and s = O(k log k). Output: The index set R of anchors for A. - 1: Initialize R = ;. q - 2: Set � < 2 2 log(4 log 2 (m)=�)= log 2 (m). p - 3: Set �V = O � min n �= k�; 1=k� 2 o� , �V = 1 � (1 � �) 1 4 . - 4: Run FKV (A; k; �V; �V) and output the description of approximate right singular vectors ˆV. - 5: Set �U = O � min n � k ; 1 k� 2 o� , �U = 1 � (1 � �) 1 4 . - 6: Run FKV (A T ; k; �U; �U) and output the description of approximate left singular vectors ˆU. - 7: By Lemma 3.5, estimate M := ˆU T A ˆV with relative error � = O � �=k 2 � � and failure probability � = 1 � (1 � �) 1 4 , and denote the result as ˜M. - 8: for i = 1 to s do - 9: Generate a unit random vector xi 2 R k . - 10: Directly compute ˜yi = ˜Mxi. - 11: By rejection sampling (Algorithm 2), simulate distribution D ˆU ˜y with failure probability � = 1 � (1 � �) 1 4s and pick up O(poly log m) samples. fLet Oi denotes the actual distribution which simulates - D ˆU˜yi g - 12: R R [ flg, where l is the most frequently index appearing in O(poly logm) samples. - 13: end for - 14: Return R necessary to generate unit vector in row space of A. The next, we prove that for each i 2 [s], distribution Oi is �-close to distribution DAVxi in total variant distance. Once again, we show how to gain the index of largest component of AVxi from distribution Oi. Finally, by O(k log k) random projection, it is enough for us to gain all anchors of matrix A. The following lemma tells us the approximate singular vectors outputted by FKV spans the row space of matrix A. And combining with Lemma 2.6, it gives us that V also spans the same space, i.e., V forms an orthonormal basis of row space of matrix A. Lemma 3.2. Let ˆV be the output of algorithmFKV (A; k; �; �). If � < 1 k� 2 , then with probability 1 � �, we obtain Proof. By contradiction, we assume that span n ˆV (i) � � �i 2 [k] o , span n A(i) � � �i 2 [m] o , which implies that there exists a unit vector x 2 span n A(i) � � �i 2 [m] o and x ? span n ˆV (i) � � �i 2 [k] o . Then we can obtain � � �Ax � A ˆV ˆV T x � � � = kAxk � �min(A) since ˆV T x = ~0. And according to Theorem 2.4, we have Thus �min(A) � p �k��min(A), which makes a contradiction if � < 1=k� 2 . � By Lemma 3.2, we can generate an approximate random vector in the row space of A with probability 1 � � in time O(poly(k; 1=�; log 1=�)) by FKV (A; k; �; �). Firstly, we obtain the description of approximate right singular vectors by FKV algorithm, where the error parameter � is bounded by rank k and condition number � (see in Lemma 3.2). Secondly, we generate a random unit vector xi 2 R k as a coordinate vector referring to a set of orthonormal vectors in Lemma 2.6. Let V denotes the matrix defined in Lemma 2.6, then it is obvious that its columns form the right singular vectors for matrix A. That is, ˆ �i = ˆVxi is an approximate vector of a random vector � = Vxi. Next, we show that total variant distance between Oi and DAVxi is bounded by constant �. For convenience, we assume that each step in Algorithm 1 succeeds and the final success probability will be given in next subsection. Lemma 3.3. For all i 2 [s], � � �Oi;DAVxi � � � TV � � holds simultaneously with probability 1 � �. In the rest, without ambiguity, we use notations O, x instead of Oi, xi. By applying triangle inequality, we divide the left part of inequality into four parts (the intuition idea please ref Figure 2): Thus, we only need to prove that 1 , 2 , 3 , 4 < � 4 , respectively. In addition, given u; v 2 R n , if ku � vk � � 2 kuk, thenkDu;DvkTV � �. For 1 , 2 and 3 , we only show their '2-norm version, i.e., - � kAVx � A ˆVxk � � 8 kAV xk; - � kA ˆVx � ˆU ˆU T A ˆVxk � � 8 kA ˆV xk; - � k ˆU ˆU T A ˆVx � ˆU ˜Mxk � � 8 k ˆU ˆU T A ˆVxk. For convenience, in the rest part, let �U = �Uk=16 and �V = �Vk=16 represent approximate ratio for orthonormality of ˆU and ˆV based on Lemma 2.7, respectively. Before we start our proof, we list two tools which are used to prove 3 and 4 , respectively. Based on rejection sampling, Lemma 3.4 shows that sampling from linear combination of �-approximately orthogonal vectors can be quickly realized without knowledge of norms of these vectors (see Algorithm 2). Algorithm 2 Rejection Sampling for D ˆUy Input: A set of approximately orthonormal vectors ˆU ( j) 2 R m (for j = 1; : : : ; k) in vector sample model and a vector y 2 R k . Output: A sample s subjecting to D ˆUy . - 1: Sample j according to probability proportional to k ˆU ( j) yjk. - 2: Sample s fromD ˆU ( j) . - 3: Compute rs = ( ˆUy) 2 s k P i(yi ˆUi j) . - 4: Accept s with probability rs, otherwise, restart. Lemma 3.4 ([Tang, 2018]). Given a set of �-approximately orthonormal vectors ˆV 2 R n�k in vector sample model, and an input vector w 2 R k , there exists an algorithm outputting a sample from a distribution � 1�� -close to DˆVw with probability 1 � � using O(k 2 log 1 (1 + O(�))) queries and samples. � Lemma 3.5. Given A 2 R m�n in matrix sample model and L 2 R k1�m and R 2 R n�k2 in query model, let M = LAR, then we can output a matrix ˜M 2 R k1�k2 , with probability 1 � �, such that Proof. Let Mi j = L(i)AR ( j) with i 2 [k1] and j 2 [k2]. In [Tang, 2018], there exists an algorithm that outputs an estimation of Mi j ( ˜Mi j) to precision �kAkF � � �L(i) � � � � � �R ( j) � � � with probability 1�� 0 in time O � 1 � 2 log 1 � 0 � . Let � 0 = 1�(1��) 1=(k1k2) . We can output ˜M with probability 1 � � utilizing O � k1k2 1 � 2 log 1 1�(1��) 1=k2 � = O � k1k2 1 � 2 log 1 � � queries and samples respectively where ˜M � proof of Lemma 3.3. Upper bound for 1 . By Lemma 3.2, Vx is a unit random vector sampled from the row space of A with probability 1 � �V := (1 � �) 1 4 if �V < 1 k� 2 . From Eq. (1) in Lemma 2.6, with probability 1 � �V Combing withkAkF � p k��min(A) � p k�kAVxk, we gain Upper bound for 2 . According to Lemma 3.2, the columns of ˆU span a space equal to the column space of A if �U � 1 k� 2 with probability 1��U := (1��) 1 4 . Let � ˆU denote the orthonormal projector to image of ˆU (column space of A). Similarly, � ? ˆU denotes the orthonormal projector to the orthogonal space of column space of ˆU. based on Eq. (2) in Lemma 2.6. If �U = O � min n � k ; 1 k� 2 o� , � � �A ˆVx � ˆU ˆU T A ˆVx � � � � � 8 � � �A ˆV x � � �. Upper bound for 3 . When �U and �V are discussed above, with probability 1 � � := (1 � �) 1 4 , we have According to Lemma 3.5, we obtain Combining Eq. (5) and Eq. (6), the following holds If � = O � � k 2 � � , then � � � ˆU ˆU T A ˆVx � ˆU ˜Mx � � � < � 8 � � � ˆU ˆU T A ˆVx � � � holds. Upper bound for 4 . Since �U = O � min n � k ; 1 k� 2 o� as discussed before, directly taking usage of Lemma 3.4, with probability 1 � � := (1 � �) 1 4s we have Hence, Algorithm 1 generates a distribution Oi which satisfies � � �Oi ; DAV xi � � � TV � � for s random unit vectors generated simultaneously with probability 1 � �. � The following theorem tells us how to find the largest component of AVxi from distribution Oi. Theorem 3.6 (Restatement of Theorem 1 in [Du et al., 2018]). Let D be a distribution over [m] and D 0 is another distribution simulating D with total variant error �. Let x1; : : : ; xN be examples independently sampled fromD 0 and Ni be the number of examples taking value of i. Let Dmax = maxfD1; : : : ;Dmg and Dsecmax = maxfD1; : : : ;DmgnDmax. If Dmax �Dsecmax > 2 p 2 log(4N=�)=N + �, then, for any � > 0, with a probability at least 1 � �, we have As mentioned in [Du et al., 2018], the assumption about the gap between Dmax and Dsecmax is easy to satisfy in practice. By choosing N = log 2 m and � < 2 q 2 log(4 log 2 m=�)= log 2 m, we have Dmax � Dsecmax > 4 q 2 log(4 log 2 m=�)= log 2 m, which will converge to zero as m goes to infinity. To estimate the number of random projections we need, we denote p � i the probability that after random projection �, a data point A(i) is identified as an anchor in subspace, i.e., In [Zhou et al., 2014], if p � i > k=� for a constant �, with s = 3 � k log k random projections, all anchors can be found with probability at least 1 � k exp(��s=3k). ## Complexity and Success Probability Note that Algorithm 1 involves operations that query and sample frommatrix A, ˆU and ˆV, but those operations can be implemented in O(log(mn)poly(k; �; 1=�)) time. Thus, in the following analysis, we just ignore the time complexity of those operations but multiple it to the final time complexity. The running time and failure probability mainly concentrates on lines 4, 6, 7 and 11 in Algorithm 1. The running time of lines 4 and 6 are O � poly(k; 1=�V ; log 1=�V ) � and O � poly(k; 1 �U ; log 1 �U ) � , respectively, according to Theorem 2.4. And line 7 takes O � k 2 1 � 2 log 1 � � to estimate matrix ˜M according to Lemma 3.5. And line 11 with s iterations totally spends O � sk 2 log 1 � poly log m � . In the perspective of failure probability, lines 4, 6 and 7 take the same failure probabilities (1 � �) 1 4 . And line 11 takes (1 � �) 1 4s for each iteration. Above all, the time complexity of FAS is O � poly � k; �; log 1 � ; logmn �� . The success probability is 1 � �. ## 4 Conclusion This paper presents a classical randomized algorithm FAS which dramatically reduces the running time to find anchors of low-rank matrix. Especially, we achieve exponential speedup when the rank is logarithmic of the input scale. Although our algorithm running in polynomial of logarithm of matrix dimension, it still has a bad dependence on rank k. In the future, we plan to improve its dependence on rank as well as analyze its noise tolerance. ## Acknowledgements Part of this work is done during Y. Li's visit at the Institute of Quantum Computing, Baidu Inc.. This work is supported in part by the National Natural Science Foundation of China Grants No. 61433014, 61832003, 61761136014, 61872334, 61502449, the Strategic Priority Research Program of Chinese Academy of Sciences Grant No. XDB28000000, and Anhui Initiative in Quantum Information Technologies, Grant No. AHY150100. Y. Li is supported by ERC Consolidator Grant 615307-QPROGRESS. Y. Li thanks Runyao Duan for hosting his visit at the Institute of Quantum Computing, Baidu Inc.. ## References 1. [Arora et al., 2012a] Sanjeev Arora, Rong Ge, Ravindran Kannan, and Ankur Moitra. Computing a nonnegative matrix factorization-provably. In Proceedings of the forty- fourth annual ACM symposium on Theory of computing, pages 145-162. ACM, 2012. 2. [Arora et al., 2012b] Sanjeev Arora, Rong Ge, and Ankur Moitra. Learning topic models-going beyond svd. In Foun- dations of Computer Science (FOCS), 2012 IEEE 53rd Annual Symposium on, pages 1-10. IEEE, 2012. 3. [Chia et al., 2018] Nai-Hui Chia, Han-Hsuan Lin, and Chun- hao Wang. Quantum-inspired sublinear classical algorithms for solving low-rank linear systems. arXiv preprint arX- iv:1811.04852, 2018. 4. [Ding et al., 2010] Chris HQ Ding, Tao Li, and Michael I Jordan. Convex and semi-nonnegative matrix factoriza- tions. IEEE transactions on pattern analysis and machine intelligence, 32(1):45-55, 2010. 5. [Donoho and Stodden, 2004] David Donoho and Victoria S- todden. When does non-negative matrix factorization give a correct decomposition into parts? In Advances in neural information processing systems, pages 1141-1148, 2004. 6. [Du et al., 2018] Yuxuan Du, Tongliang Liu, Yinan Li, Run- yao Duan, and Dacheng Tao. Quantum divide-and-conquer anchoring for separable non-negative matrix factorization. In Proceedings of the 27th International Joint Conference on Artificial Intelligence, pages 2093-2099. AAAI Press, 2018. 7. [Elhamifar and Vidal, 2009] Ehsan Elhamifar and René Vi- dal. Sparse subspace clustering. In 2009 IEEE Conference on Computer Vision and Pattern Recognition, pages 2790- 2797. IEEE, 2009. 8. [Elhamifar et al., 2012] Ehsan Elhamifar, Guillermo Sapiro, and Rene Vidal. See all by looking at a few: Sparse mod- eling for finding representative objects. In 2012 IEEE Conference on Computer Vision and Pattern Recognition, pages 1600-1607. IEEE, 2012. 9. [Esser et al., 2012] Ernie Esser, Michael Moller, Stanley Os- her, Guillermo Sapiro, and Jack Xin. A convex model for nonnegative matrix factorization and dimensionality reduction on physical space. IEEE Transactions on Image Processing, 21(7):3239-3252, 2012. 10. [Frieze et al., 2004] Alan Frieze, Ravi Kannan, and Santosh Vempala. Fast monte-carlo algorithms for finding low-rank approximations. Journal of the ACM (JACM), 51(6):1025- 1041, 2004. 11. [Gillis and Vavasis, 2014] Nicolas Gillis and Stephen A Vava- sis. Fast and robust recursive algorithmsfor separable non- negative matrix factorization. IEEE transactions on pattern analysis and machine intelligence, 36(4):698-714, 2014. 12. [Gilyén et al., 2018] András Gilyén, Seth Lloyd, and Ewin Tang. Quantum-inspired low-rank stochastic regression with logarithmic dependence on the dimension. arXiv preprint arXiv:1811.04909, 2018. 13. [Guan et al., 2012] Naiyang Guan, Dacheng Tao, Zhigang Luo, and John Shawe-Taylor. Mahnmf: Manhattan non- negative matrix factorization. stat, 1050:14, 2012. 14. [Hofmann, 2017] Thomas Hofmann. Probabilistic latent se- mantic indexing. In ACM SIGIR Forum, volume 51, pages 211-218. ACM, 2017. 15. [Hsieh and Dhillon, 2011] Cho-Jui Hsieh and Inderjit S D- hillon. Fast coordinate descent methods with variable selec- tion for non-negative matrix factorization. In Proceedings of the 17th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 1064-1072. ACM, 2011. 16. [Kim and Park, 2008] Jingu Kim and Haesun Park. Toward faster nonnegative matrix factorization: A new algorithm and comparisons. In Data Mining, 2008. ICDM'08. Eighth IEEE International Conference on, pages 353-362. IEEE, 2008. 17. [Kumar et al., 2013] Abhishek Kumar, Vikas Sindhwani, and Prabhanjan Kambadur. Fast conical hull algorithms for near-separable non-negative matrix factorization. In Inter- national Conference on Machine Learning, pages 231-239, 2013. 18. [Lee and Seung, 1999] Daniel D Lee and H Sebastian Seung. Learning the parts of objects by non-negative matrix factor- ization. Nature, 401(6755):788, 1999. 19. [Lee and Seung, 2001] Daniel D Lee and H Sebastian Seung. Algorithms for non-negative matrix factorization. In Ad- vances in neural information processing systems, pages 556-562, 2001. 20. [Lin, 2007] Chih-Jen Lin. Projected gradient methods for nonnegative matrix factorization. Neural computation, 19(10):2756-2779, 2007. 21. [Mahoney and Drineas, 2009] Michael W Mahoney and Pet- ros Drineas. Cur matrix decompositions for improved data analysis. Proceedings of the National Academy of Sciences, pages pnas-0803205106, 2009. 22. [Recht et al., 2012] Ben Recht, Christopher Re, Joel Tropp, and Victor Bittorf. Factoring nonnegative matrices with linear programs. In Advances in Neural Information Pro- cessing Systems, pages 1214-1222, 2012. 23. [Tang, 2018] Ewin Tang. A quantum-inspired classical algo- rithm for recommendation systems. arXiv preprint arX- iv:1807.04271, 2018. 24. [Vavasis, 2009] Stephen A Vavasis. On the complexity of nonnegative matrix factorization. SIAM Journal on Opti- mization, 20(3):1364-1377, 2009. 25. [Zhou et al., 2013] Tianyi Zhou, Wei Bian, and Dacheng Tao. Divide-and-conquer anchoring for near-separable nonnega- tive matrix factorization and completion in high dimensions. In 2013 IEEE 13th International Conference on Data Min- ing, pages 917-926. IEEE, 2013. 26. [Zhou et al., 2014] Tianyi Zhou, Jeff A Bilmes, and Carlos Guestrin. Divide-and-conquer learning by anchoring a conical hull. In Advances in Neural Information Processing Systems, pages 1242-1250, 2014.