Year: 2021

Abstract:
Principal component analysis (PCA) is an important method for dimensionality reduction in data science and machine learning. But, it is expensive for large matrices when only a few principal components are needed. Existing fast PCA algorithms typically assume the user will supply the number of components needed, but in practice, they may not know this number beforehand. Thus, it is important to have fast PCA algorithms depending on a tolerance. For m×n matrices where a few principal components explain most of the variance in the data, we develop one such algorithm that runs in O(mnl) time, where l ≪ min(m,n) is a small multiple of the number of principal components. We provide approximation error bounds that are within a constant factor away from optimal and demonstrate its utility with data from a variety of applications.

heading: 1 INTRODUCTION

heading: 1.1 THE TRUNCATED SVD AND PCA
text: Let A be an m×n real matrix and A = UΣV T its singular value decomposition (SVD). The rank-k truncated SVD (rank-k TSVD) of A is the matrix Ak := UkΣkV TK , where Uk and Vk are the first k columns of U and V , respectively, and Σk is the leading k × k block of Σ. The columns uj and vj of Uk and Vk are the left and right singular vectors of A, respectively, and the diagonal entries σ1(A) ≥ · · · ≥ σk(A) ≥ 0 are the singular values of A. One common use of TSVD is in principal component analysis (PCA). This is a dimensionality reduction technique that aims to find the directions in which the data varies the most. It turns out that these directions are given by the top k right singular vectors vj , 1 ≤ j ≤ k. Projecting the original data onto these directions then transforms the data from a high-dimensional space into a lower-dimensional one. Typically, k is chosen so that these principal directions explain a certain amount of variance in the data. For example, if the user wanted to explain 99% of the variance, they would choose k so that ∑k i=1 σi(A) 2/ ∑n i=1 σi(A) 2 ≥ 0.99. PCA is used to compress data in a variety of settings such as images, training data for machine learning algorithms, and even neural network weight matrices (Xue et al., 2013).
In this work, we will be interested in a slightly different version of PCA. Let ε be a user-prescribed tolerance, and instead choose k so that σk(A) ≥ ε ≥ σk+1(A). Since the ith principal component explains σi(A)2/ ∑n i=1 σi(A)
2 fraction of the total variance, we are essentially ignoring any components that explain less than ε2/ ∑n i=1 σi(A)
2 fraction of the total variance. This can be interpreted as discarding principal components corresponding to noise, where ε describes the size of the noise.
Despite its utility and importance, the main drawback of using TSVD for PCA is that it is expensive, especially when the user needs only the top few singular values/vectors. Thus, a large body of research has been devoted to finding faster ways of computing it without sacrificing too much accuracy.

heading: 1.2 PRIOR WORK
text: The literature on fast, approximate TSVD algorithms typically assume the user knows what rank k to use. Recent work uses randomization to reduce the run time while still maintaining a high level of
accuracy. See, for example, Rokhlin et al. (2010) and Halko et al. (2010) and the references therein for typical examples of these types of algorithms. In Musco & Musco (2015), the authors present a randomized algorithm based on block Krylov subspace methods to compute an approximate TSVD. For a matrix A, rank k, and tolerance ε, the algorithm produces a matrix Z whose columns approximate the top k left singular vectors of A and such that ∥∥A− ZZTA∥∥ 2 ≤ (1 + ε) ∥A−Ak∥2. They also prove stronger bounds on the quality of the singular vectors, which is important for PCA. This algorithm is especially suited to sparse matrices, which can be multiplied quickly.
In some cases, the user may not know ahead of time what k to use, so it is useful to consider algorithms which accept a desired precision ε as input rather than rank. Algorithms in this vein incrementally build a matrix Q with orthogonal columns and another matrix B until ∥A−QB∥ < ε. Typically, the number of columns of Q (or the number of rows of B) is quite small so that QB is a compact approximation of the original data A. An approximate TSVD of A can then be produced from the SVD of B. Recent work again uses randomization to reduce the run time. We present a prototypical example of this style of algorithm in Algorithm 1 (Yu et al., 2018). See Halko et al. (2010) and Martinsson & Voronin (2016) for more examples. While these algorithms guarantee a small approximation error, there are no guarantees on the accuracy of the singular values or vectors. We will compare the accuracy for Algorithm 1 to the proposed algorithm in Experiments.
Algorithm 1 The randQB EI algorithm for the fixed-precision problem Input: an m× n matrix A; desired accuracy tolerance ε; block size b Output: Q, B such that ∥A−QB∥F < ε Q = [ ]; B = [ ]; E = ∥A∥2F for i = 1, 2, 3, . . . do Ωi = randn(n, b) Qi = orth(AΩi −Q(BΩi)) Qi = orth(Qi −Q(QTQi)) Bi = Q T i A
Q = [Q,Qi]
B = [ B Bi ] E = E − ∥Bi∥2F if E < ε2 then stop
end for

heading: 1.3 OUR WORK
text: In this work, we propose an algorithm that, for a matrix A, accuracy tolerance δ, and singular value tolerance ε, produces an approximate TSVD Ãk satisfying:
1. The rank k of Ãk does not exceed the true rank of A, determined by the tolerance ε as described above,
2. σj(Ãk) ≥ (1− δ)σj(A) for 1 ≤ j ≤ k, 3. ∥∥∥A− Ãk∥∥∥
2 ≤ 1+δ1−δ ε ≈ (1 + 2δ)ε, and
4. If k coincides with the true rank, then ∥∥∥A− Ãk∥∥∥
2 ≤ (1+δ)σk+1(A) = (1+δ) ∥A−Ak∥2,
i.e. the truncation error is a factor of 1 + δ from optimal.
These properties are verified in the Appendix. The algorithm thus yields a high-quality approximation to TSVD and can be used in applications as an approximate PCA. The algorithm is fast for matrices whose singular values decay quickly when ε is set so that k will be relatively small.

heading: 2 PRELIMINARIES
text: Unless otherwise stated, we will consider matrices with more rows than columns. For a matrix with more columns than rows, apply the algorithm to its transpose.

heading: 2.1 FLIP-FLOP SPECTRUM REVEALING QR
text: The proposed algorithm is essentially a tolerance-based version of Flip-Flop Spectrum Revealing QR (FFQR) (Feng et al., 2019). For a matrix A and integers k ≤ l, FFQR produces an approximation to the rank-k TSVD Ak whose accuracy depends on the ratio σk+1(A)/σl+1(A). Thus, if A has rapidly decaying singular values, FFQR will be close to TSVD.
FFQR is computed as follows. Let A be an m × n matrix (m ≥ n) and k ≤ l. Perform l steps of Randomized QR with Column Pivoting (RQRCP) (Duersch & Gu, 2017) to get the factorization
AΠ = QR = Q ( R11 R12 0 R22 ) ,
where Π is an n× n permutation matrix, Q is an m×m orthogonal matrix, R is m× n, and R11 is an l × l upper triangular matrix. The next phase of FFQR involves performing extra “spectrum-revealing” column swaps on R and using Givens rotations to restore its upper trapezoidal form. These swaps ensure ∥R22∥2 = O(σl(A)). See Xiao et al. (2017) for more details.
Next, perform l steps of QR on RT to get
RT = PLT = (P1 P2) ( L11 0 L21 L22 )T ,
where P is an n× n orthogonal matrix, P1 is its leading l columns, L is an m× n matrix, and L11 is l × l lower triangular. Putting the above together yields
A = QRΠT = Q ( L11 0 L21 L22 ) PTΠT .
Discard L22 (as in truncated QRCP) and approximate ( L11 L21 ) with its rank-k TSVD ÛkΣ̂kV̂ Tk :
Q ( L11 0 L21 L22 ) PTΠT ≈ Q ( L11 L21 ) PT1 Π T ≈ Q(ÛkΣ̂kV̂ Tk )PT1 ΠT
Setting Ũk := QÛk, Σ̃k := Σ̂k, Ṽk := ΠP1V̂k gives the rank-k approximation A ≈ ŨkΣ̃kṼ Tk . In Feng et al. (2019), the authors prove the following bounds for FFQR. Given ε > 0 and g > 1,
there are matrix-dependent quantities g1 ≤ √ 1+ε 1−ε , g2 ≤ g, τ ≤ g1g2 √ (l + 1)(n− l), and τ̂ ≤
g1g2 √ l(n− l) such that for 1 ≤ j ≤ k,
σj(Σk) ≥ σj(A)
4 √ 1 + min ( 2τ̂4, τ4(2 + 4τ̂4) ( σl+1(A) σj(A) )4) and ∥∥∥A− ŨkΣ̃kṼ Tk ∥∥∥
2 ≤ σk+1(A)
4 √ 1 + 2τ4 ( σl+1(A)
σk+1(A)
)4 .

heading: 2.2 THE QLP DECOMPOSITION
text: The basis of FFQR is the QLP decomposition (Stewart, 1999). Let A be an m× n matrix. Perform QRCP on A to obtain AΠ = QR and then perform QRCP on RT to get RTΠ1 = PLT , where
L is lower triangular. Putting these together yields A = QΠ1LPTΠT . This is the pivoted QLP decomposition of A. Stewart observed that the diagonal entries Lii of L closely track the singular values of A.
For the proposed algorithm, we choose not to pivot when factoring RT . In this case, Lii will not track σi(A) as well. We will discuss a partial remedy for this below. The advantage of not pivoting is that, just as in FFQR, we do not have to finish computing R before computing L. Once we have performed l steps of QRCP on A, we can compute the first l rows of L and thus have access to its first l diagonal entries.

heading: 3 A FAST, APPROXIMATE, TOLERANCE-BASED PCA ALGORITHM

heading: 3.1 BLOCKED FFQR
text: Since we do not know what l is beforehand, we compute R and L incrementally in blocks. Select a block size b and perform b steps of RQRCP to get
AΠ1 = Q1
( R
(b) 11 R (b) 12 0 R (b) 22
) ,
where R(b)11 is b × b upper triangular. The first b rows of R are essentially done since subsequent steps of RQRCP will only permute the columns of R(b)12 . Perform QR on them (to keep the notation simple, we write this as an LQ factorization):(
R (b) 11 R (b) 12
) = (L11 0)P T 1 ,
where L11 is b × b lower triangular, and P1 is orthogonal. We have just computed the first b rows of L and know the first b diagonal entries. For the next block, continue RQRCP for another b steps. The permutation matrix Π2 in this block will affect only columns b + 1 through n, leaving the first
b columns untouched. Thus, Π2 can be written in block form as Π2 = ( Ib 0
0 Π̃2
) , where Ib is the
b× b identity matrix and Π̃2 is an (n− b)× (n− b) permutation matrix. We now have
AΠ1Π2 = Q2Q1  R (b) 11 R (b) 12 Π̃2
0 R (2b) 11 R (2b) 12 0 0 R (2b) 22  , where R(2b)11 is b× b upper triangular. Since the first b rows have changed, we must account for this in the previous LQ: (
R (b) 11 R (b) 12 Π̃2
) = ( R
(b) 11 R (b) 12 ) Π2 = (L11 0)P T 1 Π2.
Now apply the matrix ΠT2 P1 to the newly completed rows ( 0 R (2b) 11 R (2b) 12 ) and perform LQ on the last n− b columns to get( 0 R
(2b) 11 R (2b) 12 ) ΠT2 P1 = (L21 L22 0)P T 2 ,
where L22 is b× b lower triangular. The orthogonal matrix P2 affects only the last n − b columns and can therefore be written in block
form as P2 = ( Ib 0
0 P̃2 ) . Hence, (L11 0) = (L11 0)PT2 and(
R (b) 11 R (b) 12 Π̃2
0 R (2b) 11 R (2b) 12
) = ( L11 0 0 L21 L22 0 ) PT2 P T 1 Π2,
showing that we have computed the first 2b rows of L. We can continue this procedure, computing b rows of L at a time. Once we decide to stop, we finish the remaining rows of L by applying the
orthogonal matrices from all previous LQ factorizations to the last rows of R. For example, if we wanted to stop after 2 blocks, apply ΠT2 P1P2 to ( 0 0 R (2b) 22 ) to get(
0 0 R (2b) 22 ) ΠT2 P1P2 = (L31 L32 L33)
and the partial QLP decomposition
AΠ1Π2 = Q2Q1  R (b) 11 R (b) 12 Π̃2
0 R (2b) 11 R (2b) 12 0 0 R (2b) 22
 = Q2Q1(L11 0 0L21 L22 0 L31 L32 L33 ) PT2 P T 1 Π2.
Afterwards, spectrum-revealing swaps can be performed if desired. For each swap and uppertrapezoidal restoration, some nonzero entries will appear above the diagonal in L. These are easily eliminated with Givens rotations.
3.2 DETERMINING l
We now derive a criterion to determine when to stop factoring in blocked FFQR and to find l. Let ε be the tolerance parameter, and define the rank k by σk+1(A) ≤ ε ≤ σk(A). One could use the bounds derived in Feng et al. (2019), using the diagonal entries of L to estimate the singular values of A and stopping when σl+1(A)/σk+1(A) is sufficiently small. However, the dimension-dependent bounds for τ and τ̂ are impractical, so we will use a different bound.
In Feng et al. (2019), the authors prove that σj(A)4 ≤ σj(Σk)4+2 ∥R22∥42 , 1 ≤ j ≤ k. Rearranging this inequality gives
σj(Σk) ≥ σj(A) 4 √ 1− 2
∥R22∥42 σj(A)4 , 1 ≤ j ≤ k.
They also prove the following bound on the truncation error:∥∥∥A− ŨkΣ̃kṼ Tk ∥∥∥ 2 ≤ σk+1(A) 4 √ 1 + 2 ∥R22∥42 σk+1(A)4 . (1)
These bounds hold even without spectrum-revealing swaps. So, if ∥R22∥2 /σk+1(A), is small, then the leading k singular values of A will be revealed up to a certain number of digits and ŨkΣ̃kṼ Tk will be a nearly optimal rank-k approximation. In practice, the above two bounds are sufficient because ∥R22∥2 = O(σl(A)) already, without extra swaps. The earlier bounds still have theoretical value in that they show the algorithm works well when A has rapidly decaying singular values.
The factors 4 √ 1− 2∥R22∥ 4 2
σj(A)4 and 4
√ 1 + 2
∥R22∥42 σk+1(A)4 are equal to 1 − 12 ∥R22∥42 σj(A)4 and 1 + 12 ∥R22∥42 σk+1(A)4 ,
respectively, up to first order. Introduce an accuracy parameter δ, and say we have 12 ∥R22∥42 σk+1(A)4 ≤ δ. Then we have σj(Σ̃k) ≥ σj(A)(1− δ), 1 ≤ j ≤ k, and ∥∥∥A− ŨkΣ̃kṼ Tk ∥∥∥
2 ≤ σk+1(A)(1+ δ) up to
first order. This means that ≈ − log δ digits of the top k singular values of A and optimal truncation error have been computed correctly. We can rewrite 12 ∥R22∥42 σk+1(A)4 ≤ δ as ∥R22∥2 ≤ σk+1(A) 4 √ 2δ. This is the tolerance-based criterion to determine l.
3.2.1 ESTIMATING σk+1 AND ∥R22∥2 To estimate ∥R22∥2 and σk+1(A) accurately, we use Stewart’s observation that the diagonal entries Lii of L closely track the singular values σi(A) of A. As stated above, Lii will not track σi(A) as well because we are not pivoting when factoring RT . A partial remedy is simply to sort the Lii’s. We show below that the resulting tracking behavior is similar in quality to that of fully pivoted QLP.
Let L(j) be the j-th largest diagonal entry of L in magnitude, i.e. ∣∣L(1)∣∣ ≥ ∣∣L(2)∣∣ ≥ · · · ≥ ∣∣L(n)∣∣. In
light of Stewart’s observation, we will assume that there are constants α and β such that α ∣∣L(j)∣∣ ≤
σj(A) ≤ β ∣∣L(j)∣∣, 1 ≤ j ≤ n. The values of α and β will be estimated empirically below. A simple
way to interpret these inequalities is that for each diagonal entry Ljj , there is a singular value of A in the interval [α |Ljj | , β |Ljj |]. Consider {Ljj : β |Ljj | ≤ ε}. This just corresponds to all the intervals [α |Ljj | , β |Ljj |] contained in (−∞, ε]. For each Ljj in the set, there is a singular value σi(A) in the corresponding interval. Thus σi(A) ≤ ε. Since σk+1(A) is the largest singular value of A less than or equal to ε, we must have σi(A) ≤ σk+1(A), which then implies α |Ljj | ≤ σk+1(A). This yields a lower bound on σk+1(A), namely max{α |Ljj | : β |Ljj | ≤ ε}. Since we will not know all the Ljj’s, we can obtain only a sub-optimal lower bound sk+1 on σk+1(A). Initialize sk+1 = 0. After i blocks of blocked FFQR, update sk+1 = max{α |Ljj | : β |Ljj | ≤ ε and j ≤ ib}. To estimate ∥R22∥2 = σ1(R22), we will use the first diagonal entry L11 of L in the fully pivoted QLP factorization. First, consider a general matrix A, and perform QRCP: AΠ = QR. Then in fully pivoted QLP, |L11| is just the largest row norm of R. As noted in Stewart (1999), the largest row of R is usually among the first few rows. Thus we can estimate ∥A∥2 using max1≤i≤q ∥R(i, :)∥2, for some small integer q. This requires only q steps of QRCP. We can apply this idea to estimate ∥R22∥2. After i steps of QRCP on A, the R factor has the form( R
(i) 1
0 R (i) 22
) , where R(i)1 is i×n upper triangular. The R factor in the QRCP factorization of R (i) 22 is
just R(n)1 (i+1 : m, i+1 : n). Thus, QRCP-factoring A automatically yields the QRCP factorizations of all the trailing blocks R(i)22 . Using the 2-norm estimation scheme in the previous paragraph, after
j steps of QRCP, we have the upper bound ∥∥∥R(i)22 ∥∥∥
2 ≤ β max i≤ι≤i+q−1 ∥∥∥R(j)1 (ι, :)∥∥∥ 2 := β ∥∥∥R(i)22 ∥∥∥ j,q for
1 ≤ i ≤ j − q + 1. Putting these estimates together will give us the final stopping criterion. After each block, we first update sk+1 with the newly computed Ljj’s and then check if ∥∥∥R(i)22 ∥∥∥ j,q ≤ 1β sk+1 4 √ 2δ for some i. The smallest i for which this inequality holds will be l + 1.
Algorithm 2 Approximate, tolerance-based PCA Inputs: A, tolerance ε, accuracy δ, block size b, number of rows q, oversampling size p for RQRCP Outputs: Rank k, Ũk, Σ̃k, Ṽk c← 0, sk+1 ← 0 while c < n do
Perform steps c+ 1 to c+ b of RQRCP on A; update Q, R, and Π Compute rows c+ 1 to c+ b of L; update P for j = c+ 1 : c+ b do
if |Ljj | ≤ ε/β and α |Ljj | ≥ sk+1 then sk+1 ← α |Ljj |
end if end for for i = 1 : c+ b− q + 1 do
if ∥∥∥R(i)22 ∥∥∥
c+b,q ≤ 1β sk+1
4 √ 2δ then
l← i− 1 exit while loop
end if end for c← c+ b
end while Compute rows c+ b+ 1 to m of L. Compute TSVD ÛkΣ̂kV̂k of L(:, 1 : l), where k satisfies σk(L(:, 1 : l)) ≥ ε ≥ σk+1(L(:, 1 : l)). Ũk ← QÛk, Σ̃k ← Σ̂k, Ṽk ← ΠP1V̂k

heading: 4 EXPERIMENTS
text: We use the following test matrices for our experiments:
1. A random 3000×3000 matrix with singular values decaying geometrically from 1 down to 10−12. We generate a 3000×3000 matrix with entries from a standard normal distribution, compute its SVD UΣV T , and replace the diagonal of Σ with the desired singular value distribution.
2. A 2003×2003 data matrix (bcsstk13) from the SuiteSparse matrix collection (Davis & Hu, 2011). This matrix arises from a computational fluid dynamics problem. Its singular values decay from ≈ 1012 down to ≈ 102.
3. A 19200×5322 matrix generated from a video from the UCF-Crime dataset (Sultani et al., 2018). The original video was a 240 × 320 RGB video consisting of 5322 frames. We resized the video by half to 120× 160, converted it to grayscale, flattened each frame into a 19200× 1 column vector, and then stacked these horizontally to form the final matrix.
4. A 5000×5000 kernel matrix generated from 5000 data points from the MNIST handwritten digits dataset (Lecun et al., 1998). We used the kernel function k(x, x′) = e−γ∥x−x
′∥2 , where γ = 1/(median of pairwise distances between data points)2.
The singular values of each test matrix are plotted in Figure 1.
4.1 ESTIMATING α AND β
For each matrix A, we tested three singular value estimation schemes. See Table 1. For the first two columns (“Unpivoted”), we ran RQRCP to get AΠ = QR and then QR-factored RT = PL. We recorded the minimum (α) and maximum (β) values of σi(A)/ |Lii|. For the second two (“Sorted”), we recorded the minimum and maximum values of σi(A)/
∣∣L(i)∣∣, where L(i) is the ith largest diagonal entry of L in magnitude. For the last two columns (“Pivoted”), we QRCP-factored RTΠ1 = PL and then recorded the minimum and maximum values of σi(A)/ |Lii|. We observed that the tracking behavior can break down when σi(A) is smaller than machine precision and thus ignored ratios corresponding to such σi(A) when computing the minimum and maximum values. Therefore, it is recommended that the tolerance ε be set at least a small factor above machine epsilon.
“Sorted” and “Pivoted” have similar α and β values, with the latter slightly better overall, while “Unpivoted” tends to be worse than the other two. Based on the middle two columns, it seems that α ≈ 0.7 and β ≈ 2 are reasonable values.

heading: 4.2 COMPARISON TO TSVD AND RANDQB EI
text: Here we compare the proposed algorithm to TSVD and randQB EI. Tests were coded in Fortran and run on a laptop with a 2.00 GHz Intel i7-4510U CPU with 16.0 GB of RAM.
First, we compare the proposed algorithm to TSVD. To compute the latter, the LAPACK routine dgesdd routine is used to compute the full SVD, which is then truncated based on the tolerance ε. For the random, data, and kernel matrices, tolerances corresponding to 99% explained variance are chosen. For the video matrix, we choose one corresponding to 99.9% explained variance because the first principal component already accounts for 99% of the variance.
It is not so simple to translate a Frobenius norm tolerance to a corresponding 2-norm tolerance, but we find that for matrices with geometrically decaying singular values, 99% explained variance roughly corresponds to a tolerance of 0.1 ∥A∥2. The 2-norm of A can be estimated after the first block of FFQR. For experimental purposes, we computed the SVD of each matrix and then selected the tolerance.
We set δ = 10−4 for all test matrices because machine learning algorithms typically only need a few digits of accuracy. But the larger singular values are computed with more accuracy because the accuracy of the jth singular value depends on ∥R22∥2 /σj(A). Finally, for all matrices, we set the block size b = 64, number of rows q = 5, and oversampling size p = 5.
The results are listed in Table 2. The proposed algorithm detects the rank k correctly for each test matrix and is much faster than dgesdd. The column REOTE contains the relative error in the optimal truncation error ∣∣∣∣∥A−Ãk∥2∥A−Ak∥2 − 1 ∣∣∣∣. This is always bounded above by δ, but we see that in practice this relative error is much smaller.
Figure 2 plots the relative errors in the approximate singular values ∣∣∣1− σj(Σ̃k)σj(A) ∣∣∣, 1 ≤ j ≤ k for each of the test matrices. These are again bounded by δ. As expected, the larger singular values are computed more accurately. We also plot the relative errors in the singular values computed by randQB EI, using block size 64; ε = √ 0.01 ∥A∥F for the random, data, and kernel matrices; and
ε = √ 0.001 ∥A∥F for the video matrix. The authors of randQB EI also include a power parameter P in their implementation. We set P = 1 as in their paper.
Although we did not analyze the error in the singular vectors (in general, the error depends on the gap between the singular values), we compute the angles θ(vj , ṽj) between the right singular vectors and their approximations for both our algorithm and randQB EI, and plot them in Figure 3. For our algorithm, the angles are all quite small, so the proposed algorithm finds good-quality approximations to the true singular vectors/principal directions.
Accuracy-wise, our algorithm performs better than randQB EI, except on the kernel matrix. We found that setting P = 0 causes the accuracy of randQB EI to drop below ours. The first few
singular values of the kernel matrix are extremely large compared to the rest; thus, performing even just one power iteration effectively enhances the accuracy of randQB EI on this matrix.

heading: 5 CONCLUSION
text: In this work, we developed an efficient algorithm for computing an approximate truncated SVD. In contrast to much of the literature, this algorithm truncates according to a tolerance rather than a fixed rank. We have also demonstrated that it provides high-quality approximations to both the singular values and vectors of the original matrix, thus making it suitable for use in applications as an approximate PCA.

heading: ACKNOWLEDGMENTS
text: We would like to thank Jed Duersch for providing his code for RQRCP.

heading: A APPENDIX
text: A.1 VERIFICATION OF PROPERTIES 1-4
In this subsection, we denote the true rank as ktrue, the rank detected by our algorithm as k̃, and the matrix output by our algorithm Ãk̃. Recall ktrue is defined by σktrue+1(A) ≤ ε ≤ σktrue(A). The detected rank k̃ is determined as follows. We run Blocked FFQR until the trailing block of R is small enough and then take l to be the smallest integer such that ∥R(l + 1 : m, l + 1 : n)∥2 ≤ σktrue+1(A) 4 √ 2δ. Denote R(l + 1 : m, l + 1 : n) by R22 for short. Afterwards, compute the SVD of L1 := L(:, 1 : l) and define k̃ by σk̃+1(L1) ≤ ε ≤ σk̃(L1).
First, note that k̃ ≤ ktrue. To see this, first observe that ktrue = #{j : σj(A) > ε}. By the Cauchy Interlacing Theorem, σj(L1) ≤ σj(A), 1 ≤ j ≤ l. Thus we can only shift the singular values of A downward, which will not increase the size of the above set. This proves Property 1.
For Property 2, it follows from computations in Feng et al. (2019) that σj(A)4 ≤ σ4j (L1)+2 ∥R22∥ 4 2, or
σj(L1) ≥ σj(A) 4 √ 1− 2
∥R22∥42 σj(A)4 ≈ σj(A)
( 1− 1
2 ∥R22∥42 σj(A)4
) , 1 ≤ j ≤ l.
Plugging in ∥R22∥2 ≤ σktrue+1(A) 4 √ 2δ gives σj(L1) ≥ σj(A)(1 − δ) for 1 ≤ j ≤ ktrue + 1 and in particular for 1 ≤ j ≤ k̃. This is Property 2. For the last two properties, we refer to Equation 1. In the notation for this section, it reads:∥∥∥A− Ãk̃∥∥∥
2 ≤ σk̃+1(A) 4
√ 1 + 2
∥R22∥42 σk̃+1(A) 4 ≈ σk̃+1(A)
( 1 + 1
2 ∥R22∥42 σk̃+1(A) 4
) .
Again, plugging in ∥R22∥2 ≤ σktrue+1(A) 4 √ 2δ and using the fact that k̃ + 1 ≤ ktrue + 1 gives∥∥∥A− Ãk̃∥∥∥
2 ≤ σk̃+1(A)(1 + δ). We see that if k̃ = ktrue, then this inequality is Property 4.
Finally, from the proof of Property 2 above, we have σk̃+1(L1) ≥ σk̃+1(A)(1−δ). Thus, σk̃+1(A) ≤ 1 1−δσk̃+1(L1) and ∥∥∥A− Ãk̃∥∥∥
2 ≤ 1+δ1−δσk̃+1(L1) ≤ 1+δ 1−δ ε, which is Property 3.
Note that k̃ < ktrue only when there are singular values slightly above the tolerance. The tolerancebased criterion ensures that up to first order σj(L1) ≥ σj(A)(1 − δ). So only singular values satisfying σj(A) ≥ ε ≥ σj(A)(1− δ) can be perturbed below ε and decrease the rank.

Review summaries:
This paper presents a PCA algorithm that terminates after approximate singular values fall below a user provided threshold. The proposed method is based on the FFQR algorithm of Feng et al, but includes a tolerance based stopping criteria. It is stated that the algorithm runs $O(mnl)$ time. Experimental results are provided for comparison between TSVD and randQB_EI of Yue et al.

This paper develops an algorithm that approximates the principal components corresponding to the top singular values. The algorithm requires the user specified tolerance instead of the number of top singular values, which is often the case in practice. It presents the detailed algorithm and evaluate its performance using numerical experiments.