An Efficient and Reliable Tolerance-Based Algorithm for Principal Component Analysis
Authors: Michael Yeh
Ming Gu
Venue: ICLR 2022 Submitted (2022)
Keywords: principal component analysis
dimensionality reduction
data compression

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\times 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\ll \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.

Under review as a conference paper at ICLR 2022
AN EFFICIENT AND RELIABLE TOLERANCE-BASED ALGORITHM FOR PRINCIPAL COMPONENT ANALYSIS
Anonymous authors Paper under double-blind review
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 few principal components are needed. Existing fast PCA al- gorithms 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 X n matrices where few principal components explain most of the variance in the data, we develop one such algorithm that runs in ((mll) time, where I < 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.
INTRODUCTION
1.1 THE TRUNCATED SVD AND PCA
real matrix and A = UEVT its singular value decomposition (SVD). The rank-k truncated SVD (rank-k TSVD) of A is the matrix A := 22 V where U* and V are the first k columns of U and V, respectively, and Σg is the leading k x k block of Σ. The columns Uy and 2y of Ug and V& are the left and right singular vectors of A, respectively, and the diagonal entries are the singular values of A.
Let A be
(A) 2 a (A)
One common use of TSVD is in principal component analysis (PCA). This is a dimensionality reduction technique that aims find the directions in which the data varies the most. It turns that these directions are given by the top right singular vectors By, 1 I j I k. Projecting the original data onto these directions then transforms the data from high-dimensional space into lower-dimensional one. Typically, & is chosen so that these principal directions explain a certain amount of variance in the data. For example, if the user wanted explain 99% of the variance, they would choose so that 1 (A)) Σi= (A)2 > 0.99. PCA is used to compress data 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, will be interested in slightly different version of PCA. Let be user- eser-prrseribed tolerance, and instead choose is so that (A) 2 E 2 (A). Since the ith principal component explains o, (A)2 / o; (A)2 fraction of the total variance, ponents that explain less than i=1 o, (A) fraction of the total variance. This can be interpreted as discarding principal components corresponding to noise, where E describes the size of the noise.
essentially ignoring any com-
Despite its utility and importance, the main drawback of using TSVD for PCA is that it is expen- sive, 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.
1.2
PRIOR WORK
The literature on fast, approximate TSVD algorithms typically assume the user knows what rank use. Recent work uses randomization reduce the time while still maintaining a high level of

---

Under review as a conference paper at ICLR 2022
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 e, the algorithm produces a matrix Z whose columns approxi- mate the top & left singular vectors of A and such that A Zz" A|l, S (1 + E) ||4 - Akll2- 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 desired precision as input rather than rank. Algorithms in this vein incrementally build matrix Q with orthogonal columns and another matrix B until || - QBI < E. Typically, the number of columns of Q (or the number of rows of B) quite small so that QB is 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 reduce the run time. We present prototypical example of this style of algorithm in Algorithm Yu et al.) 2018). See Halko et al. (010)) and Martinsson & Voronin (2016) for more examples. While these algorithms guarantee small approximation error, there no guarantees the accuracy of the singular values vectors. We will compare the accuracy for Algorithm to the proposed algorithm Experiments.
Algorithm The randQB El algorithm for the fixed- sxeddprecision problem
Input: an X matrix A; desired accuracy tolerance e; block size 6
Output: Q, B such that || QB <
Q []: B [];
E |Al|
for i 1,2,3....00 do
2 randn(n, b)
Q; orth( AΩ; Q (B2)))
Q: orth(Q; Q; QQQ Q.))
B, Q A
Q = [Q. Qil
B B.
E = E |
if E < 2 then stop
end for
1.3 OUR WORK
In this work, we propose algorithm that, for a matrix A, accuracy tolerance 8, and singular value tolerance e, produces an approximate TSVD Åg satisfying:
1. The rank & of Å does not exceed the rank of A, determined by the tolerance as described above,
2. a, (As) 2 (11-8)o,(4) for 1 S S k.
3. A Av Ael1, S 1 ~ (1 + c.) and
4. If l coincides with the true rank, then | Ag S (1)))
= (1+5) ||4 |A- Asllz-
i.e. the truncation error factor of + 5 from optimal.
These properties are verified in the Appendix. The algorithm thus yields high-quality approxi- mation 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.
2

---

Under review as a conference paper at ICLR 2022
2
PRELIMINARIES
Unless otherwise stated, we will consider matrices with more rows than columns. For a matrix with more columns than rows, apply the algorithm to transpose.
2.1 FLIP-FLOP SPECTRUM REVEALING QR
The proposed algorithm is essentially folerance-based version of Flip-Flop Spectrum Revealing QR (FFQR) (Feng et al., 2019). For a matrix and integers S I, FFQR produces an approximation to the rank-k TSVD Ag whose accuracy depends the ratio 0 (A+1(4))44((A)) Thus, if A has rapidly decaying singular values, FFQR will be close to TSVD.
FFQR is computed follows. Let A be Randomized QR with Column Pivoting (RQRCP) (Duersch & Gu 2017) to get the factorization
matrix (m 2 n) and & s l. Perform / steps of
R1i
R12 R22
AlI = =R
where an n x n permutation matrix, Q is an m m orthogonal matrix, is m n, and Ru is an / x upper triangular matrix.
The next phase of FFQR involves performing extra "spectrum-revealing" column swaps on and using Givens rotations to restore its upper trapezoidal form. These swaps ensure || 1222l22 (A)) See Xiao et al. (2017) for more details.
Next, perform / steps of QR on Rd to get
Lu 0
T
RT = P2)
L21
L22
where an n x n orthogonal matrix, Pi is its leading columns, an m x n matrix, and In I x / lower triangular. Putting the above together yields
A= Q= Q
L11
0
PTTTT.
L21
L22
L11 Lzn
Discard L22 (as in truncated QRCP) and approximate
with its rank-k TSVD ÚgE.V?
o o(tn (xn) PTTITT = Q
(1l L2n
P ~ Q
Setting Ug = QŪx, En = Ex. V/2 = IP, gives the rank approximation A = V.. In Feng et al. 2019)) the authors prove the following bounds for FFQR. Given 0 and g 1, there are matrix-dependent quantities S 9192 (( I) such that for S S.S.A k.
V I- 92 9. S 9122 (1+ 1)(n l), and 7 <
(A(A)
o, 2
| 1 min 274 74(2 + 474)
(441
and
A (0mm
(A))
ok+1
2.2 THE QLP DECOMPOSITION
The basis of FFQR is the QLP decomposition (Stewart 1999). Let be an m X n matrix. Perform QRCP to obtain ATT QR and then perform QRCP RT to get RT PTT, where
3

---

Under review as a conference paper at ICLR 2022
L is lower triangular. Putting these together yields A LP"T.. This is the pivoted QLP decomposition of A. Stewart observed that the diagonal entries Lu of L closely track the singular values of A.
For the proposed algorithm, we choose not to pivot when factoring R In this case, La will not track o; (A) well. We will discuss a partial remedy for this below. The advantage of not pivoting that, just as in FFQR, we do not have to finish computing R before computing L. Once we have performed steps of QRCP A, can compute the first / rows of L and thus have access to its first diagonal entries.
3 A FAST, APPROXIMATE, TOLERANCE-BAN BASED PCA ALGORITHM
3.1 BLOCKED FFQR
Since we do not know what is beforehand, we compute R and L incrementally in blocks. Select block size b and perform b steps of RQRCP to get
AlL = Q1
12
R
66)
where 12 is b b upper triangular. The first b rows of are essentially done since subsequent steps of RQRCP will only permute the columns of F/22 Perform QR on them (to keep the notation simple, we write this an LQ factorization):
(R(
112)) (L 0)
where Ln is b x lower triangular, and P1 is orthogonal. We have just computed the first rows of L and know the first diagonal entries. For the next block, continue RQRCP for another b steps. The permutation matrix II, in this block will affect only columns b + 1 through n, leaving the first
Tb 0
b columns untouched. Thus, 1l can be written in block form as ll.
ii2
where 4n the
b x identity matrix and πl is an (n - b) (n b) permutation matrix. We now have
R22
Tl2
AI, = Q2Q1
0
R
11
(2))
R
(2))
R 26)
where (2(26) is b b upper triangular. Since the first b rows have changed, account for this in the previous LQ:
R = -n1
1E I12 = (bn1 0) 0'I''I1.
Now apply the matrix I1' Pi to the newly completed rows (o R!AM the last b columns get
(2b) and perform LQ on
(o 22
R
(2b)
IT Pi
(L21
Lm 0) 0)
where Lyg is x 6 lower triangular.
The orthogonal matrix P2 affects only the last columns and can therefore be written in block
Ib 0
form as P2
Hence, (In
0) = (Tm
0) Pf
and
P
(1
R
(b)
Il
L1 Lai
0
0
PT
I
R
L2z
0
showing that we have computed the first 2b rows of L. We can continue this procedure, computing b rows of L. at time. Once decide to stop, finish the remaining of T. by applying the
4

---

Under review as a conference paper at ICLR 2022
orthogonal matrices from all previous LQ factorizations the last rows of R. For example, if we wanted to stop after blocks, apply 11 Pi to 0
(2h)
0
R22
to get
(o
22)
IT PiP = (L31
L32 L33)
and the partial QLP decomposition
R1
R TT,
L11
0
AIIII2 Q2Q1
R1
(6))
R
R(26)
Q2Q1
L21
L22
0
Ps
Il.
R (b)
Ly1
Ly2
Lgg//
0
22
Afterwards, spectrum-revealing swaps can be performed if desired. For each swap and upper- trapezoidal restoration, some nonzero entries will appear above the diagonal in L. These are easily climinated with Givens rotations.
3.2 DETERMINING
We derive criterion determine when stop factoring in blocked FFQR and to find 1. Let e be the tolerance parameter, and define the rank &or by ( < E S (A). One could use the bounds derived in Feng al. (219)), using the diagonal entries of estimate the singular values of A and stopping when 0++1 A)/ok+ (A) is sufficiently small. However, the dimension-deeenddet bounds for 5 and & are impractical, so we will use different bound.
In Feng et al. 2019). the authors prove that oj (A)² S o, (Ex) R22l 11 S S k. Rearranging this inqquality gives
Ruz
o, 2 o,(4)
isisk.
(A)
They also prove the following bound on the truncation error:
11 Ro2ll
A-- S (a)+ 1 +2 Dx+1(A))
(1)
These bounds hold even without spectrum-revealieee swaps. So, if [| /or+1 (A), is small, then the leading l singular values of A will be revealed up to certain number of digits and Ü, E. V/T will be a nearly optimal rank-k approximation. In practice, the above two bounds are sufficient because || ( (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 - 2
(A) and V 1 +
Rea
Reall (A)+ are equal to 1 1 || and 1 + 2 (A)) (A)*
||alll
respectively, up to first order. Introduce an accuracy parameter 8. and say we have 2 (A)
RRalll
< 6.
Then we have (āx) 2 (A))(1 6). S < k, and A -- V/ S 5o++ (A)(1 8) up to first order. This means that log digits the top k singular values of and optimal truncation error have been computed correctly. We can rewrite 1 || it This is the tolerance-bar based criterion determine l.
2ax 1 S d |R22 (A+1(A) 225.
3.2.1 ESTIMATING da+1 AND ||22 l
To estimate Reall2 2 and 0x++ (A) accurately,
Stewart's observation that the diagonal entries
Li of closely track the singular values (A) of A. As stated above, I, will not track a; (A) as well because we not pivoting when factoring RT. partial remedy is simply sort the L 's.
We show below that the resulting tracking behavior is similar in quality to that of fully pivoted QLP.
Let L be the j-th largest diagonal entry L in magnitude. i.c. | > | 2 ... 2 [(()) In
light of Stewart's observation, we will assume that there are constants & and such that [ <
(A) < 3 [,,0] S The values of a and will be estimated empirically below. A simple
5

---

Under review as a conference paper at ICLR 2022
way to interpret these inequalities is that for each diagonal entry Lyj. there is singular value of A in the interval [a 1, Lnyl]- Consider {Lo 3 | S ε}. This just corresponds to all the intervals [a ا , اأأرا contained in (-oo, For each Ljj in the set, there is singular value a, (A) in the corresponding interval. Thus (A) S e. Since 1 is the largest singular value of A less than or equal to E, must have (A(A) Σ (A+)(A), which then implies |T & ak+1(A).
This yields a lower bound on (A)) namely max(a |رمبا : B L EJ..
Since we will know all the Lis's, we can obtain only a sub-optimal lower bound 8 1 on (A)) Initialize 81+1 O. After i blocks of blocked FFQR, update 8k+11 max aa Laal : 8 8|4 and < ib}.
To estimate || la = 01 (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: All = QR. Then in fully pivoted QLP, |L11| is just the largest row norm of R. As noted in Stewart (1999), the largest row of usually among the first few rows. Thus we can estimate
||Al| using max [R(i,)) 1si50
for some small integer q. This requires only steps of QRCP.
We apply this idea to estimate || After steps of QRCP on A, the R factor has the form
R(1)
R(0) 122 where R( is upper triangular. The R factor in the QRCP factorization of Ro2 is (
0
just (i+1 : m, in n.). Thus, QRCP- QRCP-factoring A automatically yields the QRCP factorizations
of all the trailing blocks R(g Using the 2-norm estimation scheme in the previous paragraph, after H for j.q
j steps of QRCP, we have the upper bound
3
max iSeSit9
[ (4, 3
1 isis 1
Putting these estimates together will give us the final stopping criterion. After each block, we first update S6+1 with the newly computed Ljj's and then check if
R2 I #88+1 20 for some i.
1.
The smallest i for which this inequality holds will be l + 1.
Algorithm 2 Approximate, tolerance-bait based PCA
Inputs: A, tolerance E, accuracy 8, block size number of rows q, oversampling size for RQRCP
Outputs: Rank k; UK. Σg, Vi c e 0, 8k+1 +-
while do
Perform steps + to c of RQRCP on A; update Q, R, and II Compute rows + 1 c b of L; update P
for j c 1 : do
if رزر S E/3 and أززاا 2 s then
8k+1 + Lzjl
end if end for for if
q + do
3 Sh- 288 then
cbb
exit while loop
end if end for
c f c
end while
Compute rows c b+b+++1 to m of L.
Compute TSVD U, Eg of L(;11: (), where k satisfies on (b(L 1: :1) 2 oeni Onti(L(; (L(;1 :)).
ÜR + QŪ. Eg <- Er V/ < ITP,

---

Under review as a conference paper at ICLR 2022
Randor
10
Date
Video
Herrel
10:9
12"
10
102
104
10"
10"
Siggua
ar vikuo
Figure 1: Singular value distributions of the test matrices
Table 1: a and values for the test matrices under the three schemes.
Unpivoted
Sorted
Pivoted
8
Q
8
Random
0.710
1.39
0.742
1.33
0.745
1.32
Data
0.523
2.75
0.793
1.17
0.817
1.17
Video
0.321
2.85
0.838
1.85
0.840
1.58
Kernel
0.635
1.60
0.802
1.32
0.817
1.31
EXPERIMENTS
We use the following test matrices for our experiments:
1. A random 3000 3000 matrix with singular values decaying geometrically from down to 10-12. We generate a 3000 x 3000 matrix with entries from a standard normal distribution, compute its SVD UEVT, and replace the diagonal of with the desired singular value distribution.
2. A 2003 x 2003 data matrix (besstk13) 13) from the SuiteSparse matrix collection (Davis & Hu 201T). This matrix arises from computational fluid dynamics problem. Its singular values decay from R 1012 down to Rz 102.
3. A 19200 x 5322 matrix generated from a video from the UCF-Crine dataset Sultani al. 2018). The original video was 240 x 320 RGB video consisting of 5322 frames. We resized the video by half to 120 x 160, converted to grayscale, flattened each frame into 19200 x 1 column vector, and then stacked these horizontally to form the final matrix.
4. A 5000 x 5000 kernel matrix generated from 5000 data points from the MNIST handwritten
digits dataset (Lecun et al., 1998). We used the kernel function k(r, x') = e >||==2/11 where 7 1/ 1/mmeiian of pairwise distances between data points)?.
The singular values of each test matrix are plotted in Figure 1
4.1 ESTIMATING AND B
For each matrix A. we tested three singular value estimation schemes. See Table For the first two columns ("Unpivoted"), we ran RQRCP to get AIT = Q and then fR-factored RT PL. We recorded the minimum (a) and maximum (3) values of (A)) L.dl For the second two ("Sorted")), we recorded the minimum and maximum values of o,(AA) [ where (() is the ith largest diago- nal entry of L in magnitude. For the last two columns ("Pivoted"), we QRCP-factored R² 111 PL and then recorded the minimum and maximum values of (A)) [Lall-
We observed that the tracking behavior can break down when (AA) is smaller than machine preci- sion and thus ignored ratios corresponding to such o, (A) when computing the minimum and max- imum values. Therefore, it is recommended that the tolerance & be set least a small factor above machine epsilon.
7

---

Under review as a conference paper at ICLR 2022
127
Randwe
10P
Data
120
Vi380
12
Kamel
124
13"
929
40.55
1T 1
17th
132
300
200
100
150
10
G.ggua ralue
Sirgu
Sirga
Figure 2: Relative singular value errors for the test matrices. Our algorithm is the bold line, randQB El is the thin one.
"Sorted" and "Pivoted" have similar and f values, with the latter slightly better overall, while "Unpivoted" tends be worse than the other two. Based on the middle two columns, it seems that 0.7 and / x 2 are reasonable values.
4.2 COMPARISON TSVD AND RANDQB
Here we compare the proposed algorithm to TSVD and randQB_a _E.. Tests were coded in Fortran and on laptop with 2.00 GHz Intel 17-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 s. For the random, data, and kernel matrices, tolerances corresponding to 99% explained variance are chosen. For the video matrix, choose corresponding to 99.9% explained variance because the first principal component already accounts for 99% of the variance.
It is not simple translate Frobenius norm tolerance
corresponding 2-norm tolerance, but we find that for matrices with geometrically decaying singular values, 99% explained variance roughly corresponds
tolerance of 0.1 |- 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 5 10-1 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 | (A). Finally, for all matrices, we set the block size b 64, number of rows g 5, and oversampling size p =
The results are listed in Table 2 The proposed algorithm detects the rank i correctly for each test matrix and is much faster than dgesdd. The column REOTE contains the relative error in the optimal
A A.
This is always bounded above by 8, but see that in practice this relative error is much smaller.
truncation error
A
Figure plots the relative errors in the approximate singular values | (Sك) ,
S j S h for each of the test matrices. These are again bounded by 8. As expected, the larger singular values are computed more accurately. We also plot the relative errors in the singular values computed by randQB_El El, using block size 64; = v0.0I ||Al| F for the random, data, and kernel matrices; and v0.001 r| for the video matrix. The authors of randQB EEI also include a power parameter P in their implementation. We set P = 1 as in their paper.
Although 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 (Vv, 0;) between the right singular vectors and their approximations for both 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 the true singular vectors/principal directions.
Accuracy-wise our algorithm performs better than randQB El, except on the kernel matrix. We found that setting P 0 causes the accuracy of randQB EI to drop below ours. The first few
8

---

Under review as a conference paper at ICLR 2022
Randon
Oata
Video
10-
Hersel
A.M.
137
12"
X70 (14)
12C
203
300
150
Singu 193ec
Degali
Figure 3: Angle between the top & right singular vectors v and their approximations ūn for the test matrices. Our algorithm is the bold line, randQB El is the thin one.
Table 2: Comparison of the proposed algorithm to TSVD.
dgesdd
Proposed algorithm
Matrix
S
k
Time(s)
1
Time (s)
REOTE
Speed-up
Random
1 10
250
14.9
804
250
2.76
1.11 X 10 16
5.39 x
Data
x 10
128
4.41
759
128
1.38
2.22 10-16
3.19x
Video
5.6 103
36
161
1692
36
50.2
9.55 X 10-15
3.21x
Kernel
7.68 x 10
72.7
370
2.75
1.11 x 10-15 15
26.4x
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.
CONCLUSION
In this work, developed efficient algorithm for computing approximate truncated SVD. In contrast to much of the literature, this algorithm truncates according to tolerance rather than fixed rank. We have also demonstrated that it provides high-quality approximations both the singular values and vectors of the original matrix, thus making it suitable for use in applications as an approximate PCA.
ACKNOWLEDGMENTS
would like to thank Jed Duersch for providing his code for RQRCP.
REFERENCES
Timothy A. Davis and Yifan Hu. The university of florida sparse matrix collection. ACM Transac- fions on Mathematical Software, 38(1):1-25, 2011.
Jed A. Duersch and Ming Gu. Randomized qr with column pivoting. SIAM Journal on Scientific Computing, C264)C223,, 2017.
Yuehua Feng, Jianwei Xiao, and Ming Gu. Low-rank matrix approximations with flip-flop spectrum- revealing qr factorization. Electronic Transactions on Numerical Analysis, 51:469-494, 2019.
Nathan Halko, Per-Gunnar Martinsson, and Joel A. Tropp. Finding structure with randomness:
Probabilistic algorithms for constructing approximate matrix decompositions, 2010.
Y. Lecun, L. Bottou, Y. Bengio, and P. Haffner. Gradient based learning applied to document recog- nition. Proceedings of the IEEE, 86(1 1):227)-22333 1998. doi: 10. 10.1109/5.726711.
Per- Per-Gunnar Martinsson and Sergey Voronin. A randomized blocked algorithm for efficiently com- puting rank-revealing factorizations matrices. SIAM Journal on Scientific Computing, 38(5): $485-$507, 2016.

---

Under review as a conference paper at ICLR 2022
Cameron Musco and Christopher Musco. Randomized block krylov methods for stronger and faster approximate singular value decomposition. In Corinna Cortes, Neil D. Lawrence, Daniel D. Lee. Masashi Sugiyama, and Roman Gareett (eds.), Advances in Neural Information Processing Systems 28: Annual Conference on Neural Information Processing Systems 2015, December 7- 12, 2015, Montreal, Quebec, Canada, pp. 1396-1404, 2015.
Vladimir Rokhlin, Arthur Szlam, and Mark Tygert. A randomized algorithm for principal compo- nent analysis. SIAM Journal on Matrix Analysis and Applications, 31(3):1100-1122 2010.
G. W. Stewart. The qlp upproximation the singular value decomposition. SIAM Journal on Scientific Computing, (144))336-1348, 1999.
Waqas Sultani, Chen Chen, and Mubarak Shah. Real-world anomaly detection in surveillance videos. In Proceedings of the IEEE Conference on Computer Vision and Puttern Recognition (CVPR), June 2018
Jianwei Xiao, Ming Gu, and Julien Langou. Fast parallel randomized qr with column pivoting algo- rithms for reliable low-rank matrix approximations. In 2017 IEEE 24th International Conference on High Performance Computing (HiPC), pp. 233-242, 2017. doi: 10. 1010///P...201..000555
Jian Xue, Jinyu Li, and Yifan Gong. Restructuring of deep neural network acoustic models with singular value decomposition. In Interspeech, January 2013.
Wenjian Yu, Yu Gu, and Yaohang Li. Efficient randomized algorithms for the fixed-precisioe low- rank matrix approximation. SIAM Journal on Matrix Analysis and Applications, 3933):1339- 1359, 2018.
A APPENDIX
A.1 VERIFICATION OF PROPERTIES 1-4
In this subsection, we denote the true rank as kirres the rank detected by our algorithm as k and
the matrix output by our algorithm Aς- Recall ktrue is defined by I S E Tkeeu (A). The detected rank is determined as follows. We run Blocked FFQR until the trailing block of R small enough and then take be the smallest integer such that || + 1: m, + 1 n)]ll S Tkrru +1 (1A) 28. Denote R(I + 1: m, I ++ 1 :n) by Rzz for short. Afterwards, compute the SVD of L = L(: :)) and define ie by
S & E (In))..
First, note that i < ktrue- To see this, first observe that ktrue = # (A) > o]. By the Cauchy Interlacing Theorem, o,(L,)) I (A), 1 < j I I. 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 2019) that or (A() S (L)) +2 Reall
or
Rizll r,(Ln) 2 od) 1- (A)
~ (A( 1 || 1
1 15i5L. 51
2 (A)4
Plugging in | S O (A)) 28 gives o, (L1) 2 (A)( - 6) for 1 S j S ktrue + 1 and in particular for S j < k. This is Property 2.
For the last two properties, we refer Equation In the notation for this section, reads:
A A Axl12 < (A)+ 1-2 Rez Tk+1 (A)*
= (A) (1 1+ 20x+1(A)
Rez 11
Again, plugging in [|ezell S Okrec+1 ( 2/2 and using the fact that & + 1 S kerue + 1 gives A Ag S 50411(1) (A)(1 + 6). We see that if ī kirues then this inequality Property 4.
10

---

Under review as a conference paper at ICLR 2022
Finally, from the proof of Property above, we have (+1((11) 2 (++ Thus, (i41(A). 2
1-8%k+1 (L1) and A A S 1 1+t1 (L1) - H which is Property 3.
Note that k < kerue only when there are singular values slightly above the tolerance. The tolerance- based criterion ensures that up to first order (اار 2 ()({ 8). So only singular values satisfying (A) 2 2 (A)(1 - 8) can be perturbed below E and decrease the rank.


Reviews:
[{"review_id": "pjH3TWh8Lc7", "reviewer": "Reviewer_ps4E", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Strengths: The paper is well written. The problem, all notations, and the derivation of the idea are all clearly presented, so the paper is easy to follow.\n\nWeaknesses: The proposed methods seems to be a non-essential variation of the FFQR algorithm in Feng et al., (2019). While specifying $\\varepsilon$ is more practical than specifying $k$ in practice and modifying the FFQR algorithm for accepting $\\varepsilon$ requires tedious deviation, the overall contribution of the paper seems incremental. The main theoretical results are essentially the same as in Feng et al., (2019).\n\nMinor suggestions:\n- In Algorithm 1, it is better to make the meaning of 'randn()' and 'orth' more clear.\n- Add a period at the ends of the third last displayed equation of page 3.\n- Line 2 of page 4, change \"Stewart\" to \"Stewart (1999)\". \n", "summary_of_the_paper": "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.\n", "summary_of_the_review": "It is a well written paper, but the contribution seems to be incremental.\n", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "iQmFiTg6Btg", "reviewer": "Reviewer_oomW", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "\nThere’s extensive work on scaling up PCA such as the randomized algorithm Halko et., al. (2010). In the prior work section, the authors argue that the existing randomized algorithms only guarantee a small approximation error, but not the accuracy of the computed singular values or space. I think the claim is not true as the accuracy of the singular values/space follows from classic results, e.g., the \nHoffman-Wielandt bound, in matrix perturbation theory. \n\nThe threshold parameter epsilon in the proposed algorithm is introduced to specify the noise level. However, this parameter may not be very useful as data is typically standardized. I think it is helpful to establish the connection between the parameter with the rank of the algorithm output.\n\nThe presentation and the experiments could be improved. This paper lacks necessary evaluation on real-life datasets. The authors report numerical results of the proposed algorithm on synthetic data. For a fair evaluation, I think it would be helpful to include the state-of-the-art randomized algorithms as well as the robust PCA as baselines.", "summary_of_the_paper": "In this paper, the authors propose an approximate algorithm for PCA on large data. The proposed algorithm computes only the principle components associated with the singular values larger than a given threshold.", "summary_of_the_review": "I think the originality and evaluation of the paper need to be strengthened, including comprehensive comparison with state-of-the-art randomized algorithms as well as empirical evaluation on real-life data.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "1: The contributions are neither significant nor novel."}, "extra_text": {}}, {"review_id": "QPYo2aFZevM", "reviewer": "Reviewer_ELuj", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "## Strengths\n- PCA is an important problem that has attracted the attention of theoreticians and practitioners for many years. Accordingly, the paper is well motivated and certainly of interest to many. \n- Parameter tuning (for sparsity level or regularization) is taken for granted but often relies on heuristics or solving a similar problem many times which is far from ideal. As such, a tolerance based approach with low time complexity is appealing and a valuable contribution. \n- Their method seems to perform quite well compared to Algorithm 1 on realistic problems (examples 2-4) encountered in the world. \n- The authors also do a nice job summarizing several of the methods their work is based on which lays a good foundation for their paper. \n\n## Weaknesses\n- There is a claim concerning the time complexity in the abstract but no further discussion. It would strengthen the paper if the authors presented their analysis. \n- It is also unclear if the complexity analysis assumes access to the singular values of $A$, or to the matrix that tracks the singular values, i.e., $L$. \n- In general, the paper provides mathematical reasoning, however the informal structure makes it difficult to follow arguments and whether or not the implemented algorithm satisfies the properties in 1.3 and what conditions are necessary for the stated time complexity. \n- Although the method eliminates the need to provide a target rank, the proposed algorithm requires several other parameters, in particular $\\alpha$, $\\beta$, $b$, and $\\epsilon$. From the discussion in 3.2.1, 4.1, and 4.2 these values must be considered well before choosing. It seems that we have traded one unknown parameter for several others. \n\n\n## Questions/comments:\n- Section 1.1\n    - Should be lower case K for definition of TSVD\n- Section 1.3\n    - The \"true rank\" as defined in the appendix does not appear in the main text body (rank $k$ appears at the beginning of 3.2). This is confusing since true rank suggests number of singular values greater than zero...not greater than some threshold. I recommend using a different name.\n    - Bullet point 4: if $k$ is the true rank of $k$, then $\\sigma_{k+1}(A) = 0$, correct? Should this be $\\sigma_{k+1}(\\tilde A_k) = 0$? This seems to contradict bullet point 1 (maybe this is based on the alternative definition of rank?)\n- Section 3.2.1\n    - Do we need to assume the there exist such $\\alpha$ and $\\beta$? This would be true for a very small $\\alpha$ and very large \n- Section 4.2\n    - How does runtime compare between randQB and the proposed method? A column in table 2 would be helpful in evaluating the methods merit. \n", "summary_of_the_paper": "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.", "summary_of_the_review": "The problem is well motivated and of considerable interest to the community, the numerical experiments appear promising as well. On the other hand, formal proofs are missing as is any discussion of time complexity results which seems to be a central contribution. At the moment, I believe the paper is below the acceptance threshold. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Meta-review: The reviewers recommended rejection. There was no reply from the authors. The main weaknesses are:
- No experiment on real-life dataset (only simulated)
- Unsubstantiated claims about the literature 
- No discussion on the time complexity
- Incremental contribution
Decision: Reject