Title: SHIFTED RANDOMIZED SINGULAR VALUE DECOMPO-SITION

Abstract: We extend the randomized singular value decomposition (SVD) algorithm (Halko et al., 2011) to estimate the SVD of a shifted data matrix without explicitly constructing the matrix in the memory. With no loss in the accuracy of the original algorithm, the extended algorithm provides for a more efficient way of matrix factorization. The algorithm facilitates the low-rank approximation and principal component analysis (PCA) of off-center data matrices. When applied to different types of data matrices, our experimental results confirm the advantages of the extensions made to the original algorithm.

1. INTRODUCTION:

The singular value decomposition (SVD) is one of the most used matrix decompositions in many areas of science. Among the typical applications of SVD are the low-rank matrix approximation and principal component analysis (PCA) of data matrices (Jolliffe, 2002) . Using SVD to accurately estimate a low-rank factorization or the principal components of a data matrix, a mean-centering step should be carried out before performing SVD on the matrix. Despite its simplicity, the meancentering can be very costly if the data matrix is large and sparse. This cost is because the mean subtraction of a sparse matrix turns it to a dense matrix which requires a considerable amount of memory and CPU time to be analyzed. This motivates us to extend the randomized SVD algorithm introduced by (Halko et al., 2011) to estimate the singular value decomposition of a mean-centered matrix without explicitly forming the matrix in the memory.

More generally, we introduce a shifted randomized SVD algorithm that provides for the SVD estimation of a data matrix shifted by any vector in the eigenspace of its column vectors. The proposed algorithm facilitates the low-rank matrix approximation and the principal component analysis of a data matrix through merging the mean-centering and the SVD steps. The mean-centering is crucial to obtain the minimum PCA reconstruction error through a deterministic SVD. We experimentally show that it plays an essential role in case of using the randomized SVD algorithm too. Our experiments with different types of data matrices show that the extended algorithm performs better than the original algorithm when both are applied to a center-off data matrix.

In the followings, we briefly introduce the principal component analysis and its connection with the singular value decomposition. Then, we introduce the shifted randomized SVD algorithm and provide an analysis of its performance. Finally, we report our experimental results obtained from the principal component analysis of different types of data matrices using the extended and the original randomized SVD algorithms.

2. PRINCIPAL COMPONENT ANALYSIS:

Principal component analysis (PCA) is a method to study the variance of a random vector. PCA projects a random vector to the eigenspace of its covariance matrix. Let x be an m-dimensional random vector with the mean vector 0. PCA projects x to a latent random vector y as below:

y = A T x (1)

where the square matrix A is composed of the eigenvectors of the covariance matrix Σ x . The elements of y are called the principal components of x. In many use cases of PCA, A contains a subset of the eigenvectors of Σ x . The minimum PCA reconstruction error is obtained from the eigenvectors corresponding to the top eigenvalues of the covariance matrix (Jolliffe, 2002) .

The matrix of eigenvectors A in Equation 1 can be efficiently estimated from the singular value decomposition (SVD) of a sample matrix X. To this end, the sample matrix should be first centered around its mean vector:

X = X -µ x 1 T (2)

The fact that the matrix of left singular vectors of X is equal to the eigenvectors of the covariance matrix of X, the PCA projection of X is:

Y = U T X = SV T (3)

where X = U SV T .

The mean-centering step in Equation 2 can be very costly if X is a large sparse matrix and µ x = 0. In this case, X is a dense matrix that requires a vast amount of memory and cannot be processed in a reasonable time. In the next section, we introduce a randomized SVD algorithm to estimate the SVD of X without explicitly performing the mean subtraction step.

3. SHIFTED SINGULAR VALUE DECOMPOSITION:

Let X be an m × n (m ≤ n) matrix and µ be an m dimensional vector in the space of the column vectors of X. Algorithm 1, extends the randomized matrix factorization method introduced by Halko et al. (2011) to return a rank-k approximation of the singular value decomposition of the matrix X = X -µ1 T without explicitly forming X. The differences between the extended algorithm and the original one are in lines 6, 9, 10, and 12. In the followings, we explain the entire algorithm with a more in-depth focus on the modified parts.

The SHIFTED-RANDOMIZED-SVD algorithm works in three major steps:

1. Estimate a basis matrix for X 2. Project X to the space of the basis matrix 3. Estimate the SVD factors of X from its projection

In the first step (lines 2 to 11), a rank K basis matrix Q 1 (k < K m) that spans the column vectors of X is computed. In Line 2, a random matrix is drawn from the standard Gaussian distribution. This matrix is then used in Line 3 to form the sample matrix X 1 whose columns are independent random points in the range of X. This sample matrix is used to estimate a basis matrix for X in two steps. In Line 4, a basis matrix Q 1 is computed through QR-factorization of X 1 . Since X 1 is sampled from X, the basis matrix is considered as an approximation of the basis of X. Then in Line 6, the basis of X is estimated from the Q 1 by the QR-update algorithm proposed by Golub & Van Loan (1996, p. 607) . For a given QR factorization such as Q 1 R 1 = X 1 and two vectors u, and v, the QR-update algorithm computes the QR-factorization in Equation 4 by updating the already available factors

Q 1 and R 1 . QR = X 1 + uv T (4)

The computational complexity of the QR-update of the m × K matrix X 1 is O(m2 ).1 Replacing u with -µ and v with 1, the QR-update in Line 6 returns the basis matrix Q that span the range of the matrix X = X -µ1 T (5)

In other words, X can be approximated from Q:

X ≈ QQ T X (6)

Note that the basis matrix of X is computed without explicitly constructing the matrix X itself. The if statement in Line 4 controls the useless performance of the QR-update step with the null vector.

The for loop starting at Line 8, estimates a basis matrix for B = ( X XT ) q X using the basis of X, Q. The matrix B with a positive integer power has the same singular vectors as X, but with a sharper decay in its singular values since s j (B) = s j ( X) 2q+1 , where s j (.) returns the jth singular vectors of its input matrix. The sharp decay in singular values improves the accuracy of the randomized SVD when the singular values of X decay slowly. This effect is because the reconstruction error of the randomized SVD is directly proportional to the first top unused singular vector of X (see Equation 12).

The basis of ( X XT ) q X is computed via alternative applications of matrix product on XT and X.

For q = 1, in Line 9, a basis matrix of XT X is estimated through QR-factorization of XT Q. To avoid forming X explicitly, instead of direct multiplication XT Q, we use the distributive property of multiplication over addition:

XT Q = (X -µ1 T ) T Q = X T Q -1(µ T Q) (7

) where 1 is a vector of ones. The product 1(µ T Q) can be efficiently computed in O(nK) memory space if a higher priority is given to the parentheses. In Line 10, a basis matrix of X XT X is estimated through QR-factorization of XQ where Q is a basis matrix of XT Q. Similar to Equation 8, the product XQ is computed as:

XQ = (X -µ1 T )Q = XQ -µ(1 T Q ) (8)

with the same amount of memory space, O(nK). The matrix multiplication loop iterates q times. At this stage, we have the basis matrix Q that approximates a basis for X.

In the second major step of the algorithm, the matrix X is projected to the space spanned by Q: Y = Q T X (9) This step in done in Line 12 using the same trick as in Equation 8:

Y = Q T (X -µ1 T ) = Q T X -(Q T µ)1 T (10)

Finally, in the third step, the SVD factors of X are estimated from the K × n matrix Y in two steps. First, a rank-k approximation of Y is computed using a standard method of singular value decomposition, i.e., Y = U 1 ΣV T (Line 13). Then, the left singular vectors are updated by U = QU 1 resulting in U ΣV T = QY (Line 14). Replacing Y with Q T X and using Equation 6 ( X ≈ QQ T X), we have the rank-k approximation of X:

U ΣV T ≈ X (11)

Algorithm 1 The rank-k singular value decomposition of the m × n matrix X -µ1 T = U ΣV T with (m ≤ n) using the sampling parameter K (k < K m) and q ∈ {0, 1, 2, . . . }.

1: procedure SHIFTED-RANDOMIZED-SVD(X, µ, k, K, q) 2: Draw an n × K standard Gaussian matrix Ω 3: Form the sample matrix X 1 ← XΩ 4: Compute the QR factorization X 1 = Q 1 R 1 5: if µ = 0 then 6: Compute QR = Q 1 R 1 -µ1 T using the QR-update algorithm 7: end if 8: for i = 1, 2, . . . , q do 9: Compute the QR-factorization Q R = X T Q -1(µ T Q) 10: Compute the QR-factorization QR = XQ -µ(1 T Q ) 11: end for 12: Form Y ← Q T X -(Q T µ)1 T 13: Compute the singular value decomposition of Y = U 1 ΣV T 14: U ← QU 1 15: return (U, Σ, V ) 16: end procedure

The shifting vector µ in Algorithm 1 can be any vector in the space of the column vectors of X. If it is set to the null vector 0, then the algorithm reduces to the original randomized SVD algorithm of Halko et al. (2011) . If it is set to the mean vector of X, then the algorithm estimates the singular vectors of the mean-centered matrix X. In this case, the algorithm facilitates the principal component analysis of a data matrix X through merging the centering step in Equation 2 and the SVD step in Equation 3.

The SHIFTED-RANDOMIZED-SVD algorithm explained in the previous section approximates the SVD of a shifted data matrix X = X -µ1 T without explicitly constructing the matrix in the memory. In this section, we study the performance of the algorithm based on the accuracy and the efficiency of the original randomized SVD algorithm of Halko et al. (2011) .

To estimate the singular value decomposition of a sifted matrix X using the original randomized SVD algorithm, X should be explicitly formed and passed to the algorithm. Since SHIFTED-RANDOMIZED-SVD adds no extra randomness to the original algorithm, we have the same reconstruction error bound as if the original algorithm factorized the shifted matrix X:

EQUATION

)

where σ k+1 is the (k + 1)th singular value of the m × n matrix X with m ≤ n, 2 ≤ k ≤ m 2 is the decomposition rank, and q ∈ Z + is a power value as explained in Algorithm 1.

In the followings, we study the computational complexity of the SVD factorization of X using the original randomized SVD algorithm and its extended version in Algorithm 1. For an m × n matrix X, the computational complexity of the original randomized SVD algorithm is:

EQUATION

)

where α is the cost of the matrix-vector multiplication with the input matrix X. If X is a dense matrix then α = mn , and if X is a sparse matrix then α = T , a small constant value.

Algorithm 1 adds a QR-update step (Line 6) and three matrix-matrix multiplications (lines 9, 10, and 12) to the original algorithm. The matrix multiplications do not affect the computational complexity of the original algorithm since their computational complexity is equal to the complexity of computing Q T X in the original algorithm. The QR-update step, running in O(m 2 ), however, can affect the computational complexity of the algorithm.

Assuming that µ = 0, if both X and X are dense matrices then both algorithms have the same computational complexity as:

O(mnk + (m + n)k 2 ) (14)

This equality is because the computational complexity of the QR-update step, O(m 2 ), is dominated by the complexity of the original algorithm (i.e., m 2 ≤ mn for every m ≤ n where m, n ∈ N). In addition, the construction of X to be use by the original algorithm takes O(mn) time which is greater than or equal to the complexity of the QR-update step. Hence, the added operations do not affect the computational complexity of the original algorithm. Halko et al. (2011) show that for a dense input matrix, the randomized SVD algorithm can be performed in O(mn log k + (m + n)k 2 ) if instead of the random normal matrix Ω in Line 2, a structured random matrix such as the sub-sampled random Fourier transform is used. This improvement can be considered for the SHIFTED-RANDOMIZED-SVD algorithm too.

If the input matrix X is sparse, then X is a dense matrix for every µ = 0. In this case, the computational complexity of the SHIFTED-RANDOMIZED-SVD algorithm is:

O(T k + m 2 + (m + n)k 2 ) (15)

where the constant T is the cost of multiplying a sparse matrix to a vector, and the parameter m 2 is related to the complexity of the QR-update step. On the other hand, since X is a dense matrix, the complexity of the original algorithm is O(mnk + (m + n)k 2 ) which is higher than the complexity of the extended algorithm.

In a special case where X is a dense matrix and X is a sparse matrix, the original algorithm can factorize X in O(T k +(m+n)k 2 ), but SHIFTED-RANDOMIZED-SVD needs O(mnk +(m+n)k 2 ) time. In this case, if Algorithm 1 is applied to X with µ = 0, the factorization can be performed in the same processing time as the original algorithm. As a summary, we showed that the SHIFTED-RANDOMIZED-SVD algorithm as illustrated in Algorithm 1 is as efficient as the randomized SVD algorithm proposed by Halko et al. (2011) if the input matrix is dense, and more efficient than it if the input matrix is sparse.

5. EXPERIMENTS:

We experimentally study the difference between performing PCA with the randomized SVD algorithm (RSVD) proposed by Halko et al. (2011) and its extended version in Algorithm 1 (S-RSVD).

The fact that a minimum PCA reconstruction error is obtained from the deterministic SVD of a mean-centered data matrix, the performance of S-RSVD on an off-center data matrix with its mean vector as the shifting vector is expected to be more accurate than the performance of RSVD on the same data matrix. Since the two algorithms are randomized in nature, it is important to test if this expectation is valid for different types of data matrices.

We experimentally compare the two algorithms based on the mean of the squared L 2 norm of PCA reconstruction error (MSE). The same parameters K = 2k and q = 0 are used for both S-RSVD and RSVD unless it is clearly mentioned. The shifting vector µ for S-RSVD is set to the mean vector of data matrices. The experiments are carried out on different types of data matrices including randomly generated data, a set of images, and word co-occurrence matrices. The characteristics of the data matrices are illustrated in the corresponding sections.

5.1. RANDOM DATA:

In this section, we examine how the two SVD algorithms are affected by the parameters such as the number of principal components, the size and the distribution of a data matrix, and the power iteration scheme. Our experiments are based on two comparison metrics, 1) an MSE value obtained from a fixed number of principal components, and 2) the sum of MSE values obtained from different number of principal components ranging from 1 to 100.

Figure 1a represents the effect of the number of principal components on the MSE values obtained from a 100 × 1000 matrix sampled from a 100-dimensional random vector uniformly distribution in the range 0, 1. The results show that mean-centering leads to substantial reduction to the reconstruction error when the number of principal components is small. This observation is in line with the fact that the contribution of the mean-centering is mostly to the accuracy of top principal components.

The effect of the sample size on the two factorization algorithms is represented in Figure 1b in which the X-axis is the sample size and the Y-axis is the sum of the MSE values obtained from different number of principal components ranging from 1 to 100. The data matrices are generated by a 100-dimensional uniform random vector in the range 0, 1. The results show that S-RSVD is more accurate and stable than RSVD. Despite the fact that both algorithms are randomized, the stability of S-RSVD is less sensitive to the sample size.

Figure 1c compares the performance of the two algorithms with respect to the data distribution. The Y-axis is the sum of MSE values.Regardless of the data distribution, S-RSVD is more accurate than RSVD. This observation is in line with the fact that PCA does not make any assumption about the data distribution.

To examine whether both algorithms are equally accurate for the factorization of a mean-centered data matrix X, a comparison between the two algorithms is provided in Figure 1d . The Y-axis is the sum of the MSE values.In the experiments with RSVD, the matrix X is explicitly constructed and factorized, but in the S-RSVD experiments the singular factors of X are implicitly estimated from X. The results show that S-RSVD is as accurate as RSVD applied to an already centered data matrix X. This observation supports Equation 11.

An important parameter of the randomized SVD algorithm is the power value q that has a positive effect on the accuracy of the algorithm, . Figure 1e shows the MSE values obtained from each of the factorization algorithms with different values of q. The data matrix in this experiment is sampled from a 100-dimensional uniform distribution. The Y-axis is the sum of MSE values and the X-axis represents q. The importance of mean-centering is clear when the value of q is small. The accuracy of RSVD is significantly improved as the values of q increases, while the accuracy of S-RSVD is only slightly improved. This observation on a set of uniformly distributed vectors suggest that RSVD with a positive value of q (1 or 2 as suggested by Halko et al. (2011) ) can be as accurate as S-RSVD.

To test whether RSVD with a large value of q can be as accurate as S-RSVD, we run the same experiment as above but on data with different distributions. Figure 1f shows the difference between q, and (f) the difference of their accuracies with respect to the power value q and the data distribution.

the sum of MSE values obtained from each of the algorithms (i.e., Y-axis is MSE-SUM(S-RSVD) -MSE-SUM(S-RSVD)) with respect to the parameter q. Being all the results negative means that S-RSVD is more accurate than RSVD. Except for the data with Zipfian distribution, the difference between the accuracy of the two algorithms approaches to zero as the value of q increases. The Zipfian graph fluctuates widely for small values of q, but it becomes flat as q becomes larger. In the best case, the difference between the two algorithms on Zipfian data is -64 at q = 200. This indicates that the power iteration scheme cannot fully recover the reconstruction loss of an offcenter data matrix, but depending on the data distribution, it can be helpful. The power iteration in Algorithm 1 is a computationally heavy step which can negatively affect the efficiency of the algorithm when the value of q is large.

5.2. IMAGE DATA:

In this section, we experiment with handwritten digits and facial image matrices. The handwritten digits are a copy of the test set of the UCI ML hand-written digits datasets consisting of 1979 images of size 8 × 8.2 We vectorize individual image matrices and stack them into a single 16 × 1979 data matrix. The facial images consisting of 13233 images each of size 250 × 250 are downloaded from Labeled Faces in the Wild (LFW). 3 The facial image matrix after vectorizing and stacking all of the images matrices is a 62500 × 13233 matrix.

The left side of Table 1 summarizes the results obtained from 10-dimensional PCA of the image matrices.The MSE values represented in the first two rows of the table show that S-RSVD is more accurate than RSVD. To ensure that the results are not due to chance, we run the experiment 30 times and perform two t-tests with the following null hypotheses:

• H 1 0 :there is no difference between the MSE of S-RSVD and RSVD. • H 2 0 :there is no difference between the individual column reconstruction errors of S-RSVD and RSVD. The former hypothesis is validated on the 30 MSE pairs obtained from the SVD methods, but the later hypothesis is validate on the pairs of the reconstruction error of individual images. The p-values represented in Table 1 reject both hypotheses and confirm that the better results obtained from S-RSVD are not by chance. The rejection of H 2 0 indicates that S-RSVD results in not only lower MSE for the entire image matrices, but also for individual images.

To provide a better picture of how well the SVD algorithms perform on individual images, we plot the first 10 images of each data matrix and estimate the win-rates of the algorithms. Figure 2 shows the examples of the original handwritten and facial images (the top rows), and their reconstructions using S-RSVD (the middle rows) and RSVD (the bottom rows) with the reconstruction error values on top of each image. For most images, S-RSVD is more accurate than RSVD. To generalize this observation, we estimate the win-rate (WR) of the algorithms (i.e., the number of images for which one algorithm is more accurate than the other algorithm out of the total number of images). The results shown in Table 1 indicate that 66% of the handwritten images and 82% of the facial images are reconstructed more accurately by S-RSVD than RSVD.

:

The computational complexity of the QR-update of an m × n matrix in O(N

) where N = max(m, n).

https://archive.ics.uci.edu/ml/datasets/Optical+Recognition+of+ Handwritten+Digits

http://vis-www.cs.umass.edu/lfw/lfw.tgz

Figures and Tables:

Figure1: A comparison between S-RSVD and RSVD based on (a) number of principal components, (b) sample size, (c) data distribution, (d) explicit versus implicit mean-centering, (e) the power value q, and (f) the difference of their accuracies with respect to the power value q and the data distribution.

1e3 n = 1e4 n = 1e5 n = 3e5 MSE of S-RSVD 415.7 15.3e7 195e-5 235e-5 763e-5 994e-5 MSE of RSVD 430.6 16.1e7 200e-5 236e-5 765e-5 998e-5 p 1 -value 0

Figure 2: The effect of mean-centering on (a) handwritten digits and (b) facial data. In each subfigure, the top row is the original image, the rows in the middle and bottom are the S-RSVD and RSVD reconstructed images, respectively. The reconstruction errors are shown on top of each image.

The reconstruction error statistics of image and word data.

Table content: <table><tr><td/><td>.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>p 2 -value</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>WR of S-RSVD</td><td>66%</td><td>82%</td><td>71%</td><td>73%</td><td>77%</td><td>70%</td></tr><tr><td>WR of RSVD</td><td>34%</td><td>18%</td><td>29%</td><td>27%</td><td>23%</td><td>30%</td></tr></table>

Peer review (iclr2020):
This paper extends the existing randomized SVD algorithm by Halko et al. to propose a shifted randomized SVD algorithm. The proposed algorithm performs randomized SVD on _x0008_ar{X}=X-mu 1^T with a given constant vector mu without explicitly constructing _x0008_ar{X} from X. The proposed algorithm seems to be a straightforward modification of the randomized SVD algorithm, so that the contribution of this paper would be incremental. Also, the claimed benefits of the proposed method are better accuracy, less memory usage, and less computational time. On memory usage and computational time, however, I could not find any experimental comparison in Section 5. At least in Section 5.3 where the data matrix is indeed sparse, there should be some comparison in these respects. Because of these reasons, I would not be able to recommend acceptance of this paper. In Algorithm 1, line 6, I think that one has to subtract not mu 1^T but mu 1^T Omega. I guess that equation (12) comes from Theorem 1.2 in Halko et al. (2011). An explicit citation should be needed here to put an appropriate credit to that paper. Also, the conditions under which equation (12) holds should be explicitly stated. In fact, Theorem 1.2 in Halko et al. is about the expected reconstruction error with rank-2k factorization obtained by randomized SVD, which is not consistent with the authors claim in page 3, lines 4-5: The first top unused singular value is not sigma_{k+1} but sigma_{2k+1}. Page 1, line 25: that provides (for) the SVD estimation Page 2, line 38: that span(s) the range Page 4, line 6: a s(h)ifted matrix Page 4, line 29: to be use(d) by the Page 5, line 22: uniformly distribut(ion -> ed) Page 5, lines 33, 38: Space is missing after the period. Page 5, line 43: There is an extra comma. Page 5, line 48: suggest(s) that Page 6, line 4: approaches (to) zero Page 6, line 14: 16 	imes 1979 should read 64 	imes 1979. Page 7, line 2: is validate(d)