# Split-and-Combine Singular Value Decomposition for Large-Scale Matrix

Authors: Jengnan Tzeng
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

The singular value decomposition (SVD) is a fundamental matrix decomposition in linear algebra. It is widely applied in many modern techniques, for example, high- dimensional data visualization, dimension reduction, data mining, latent semantic analysis, and so forth. Although the SVD plays an essential role in these fields, its apparent weakness is the order three computational cost. This order three computational cost makes many modern applications infeasible, especially when the scale of the data is huge and growing. Therefore, it is imperative to develop a fast SVD method in modern era. If the rank of matrix is much smaller than the matrix size, there are already some fast SVD approaches. In this paper, we focus on this case but with the additional condition that the data is considerably huge to be stored as a matrix form. We will demonstrate that this fast SVD result is sufficiently accurate, and most importantly it can be derived immediately. Using this fast method, many infeasible modern techniques based on the SVD will become viable.


---

## Body

## 1. Introduction

The singular value decomposition (SVD) and the principle component analysis (PCA) are fundamental in linear algebra and statistics. There are many modern applications based on these two tools, such as linear discriminate analysis [1], multidimensional scaling analysis [2], and feature extraction, high-dimensional data visualization. In recent years, digital information has been proliferating and many analytic methods based on the PCA and the SVD are facing the challenge of their significant computational cost. Thus, it is crucial to develop a fast approach to compute the PCA and the SVD.Currently there are some well-known methods for computing the SVD. For example, the GR-SVD is a two-step method which performs Householder transformations to reduce the matrix to bidiagonal form then performs theQR iteration to obtain the singular values [3, 4]. Since the off-diagonal regions are used to store the transform information, this approach is very efficient in saving the computational memory. If we only want to compute a few of the largest singular values and associated singular vectors of a large matrix, the Lanczos bidiagonalization is an important procedure for solving this problem [5–8]. However, the previously mentioned methods all require matrix multiplications for the SVD. One interesting problem is how do we compute the SVD for a matrix when the matrix size is huge and loading the whole matrix into the memory is not possible? The main purpose of this paper is to deal with this problem when the numerical rank of the huge matrix is small.The second purpose of this paper is to update the SVD when the matrix size is extended by new data updating. If the rank of matrix is much smaller than the matrix size, Matthew proposed a fast SVD updating method for the low-rank matrix in 2006 [9]. A rank-r thin SVD of an m×n matrix can be computed in O(mnr) time for r≤min(m,n). Although the matrix size remains unchanged in the Matthew’s method, we can still adopt the analysis of this paper for updating the SVD when the matrix size is changed.Multidimensional scaling (MDS) is a method of representing the high-dimensional data into the low-dimensional configuration [10–12]. One of the key approaches of the MDS is simply the SVD, that is, if we can find a fast approach of the MDS then it is possible to find a fast approach of the SVD. When the data configuration is Euclidean, the MDS is similar to the PCA, in that both can remove inherent noise with its compact representation of data. The order three computational complexity makes it infeasible to apply to huge data, for example, when the sample size is more than one million. In 2008, Tzeng et al. developed a fast multidimensional scaling method which turned the classical order three MDS method to be linear [13]. In this paper, we would like to implement SCMDS to the fast SVD approach, say SCSVD. The following subsections are reviews of the classical MDS and the SCSVD.


### 1.1. Review of the Classical MDS

Assume thatX is an n-by-r matrix with rank r, where the ith row of X represents the coordinates of the ith data point in an r-dimensional space. Let D=XXT be the product matrix of X. Torgerson [10] proposed the first classical MDS method to reconstruct a matrix X of Cartesian coordinates of points from a given symmetric matrix D. The key of the classical MDS is to apply the double centering operator and the singular value decomposition.Assume that1˙ is an n-by-1 matrix whose elements are all 1’s. We define a symmetric matrix B by
(1)B=(X-1n1˙1˙TX)(X-1n1˙1˙TX)T=XXT-1nXXT1˙1˙T-1n1˙1˙TXXT+1n21˙1˙TXXT1˙1˙T=D-1nD1˙1˙T-1n1˙1˙TD+1n21˙1˙TD1˙1˙T=D-DR-DC+DG,
where DR=(1/n)D1˙1˙T is the row mean matrix of D, DC=(1/n)1˙1˙TD is the column mean matrix of D, and DG=(1/n2)1˙1˙TXXT1˙1˙T is the ground mean matrix of D. The operator from D to D-DR-DC+DG is called double centering. If we define a matrix H by
(2)H=I-1n1˙1˙T,B can be simplified to B=HDH. Since matrix B is symmetric, the SVD decomposes B into B=VΣVT. Then we have
(3)B=VΣ1/2PT=X-1n1˙1˙TX,
for some unitary matrix P. In practice, we set P=I to obtain the MDS result, namely, B.AlthoughB≠X, B=VΣ1/2=(X-(1/n)1˙1˙TX)P for some unitary matrix P, it is the same as shifting all the rows of X such that the row’s mean is zero and then multiplying the result by P. Hence, the row of B preserves the pairwise relationship among the data points. B is derived from D by double centering, whose cost is lower than obtaining B by the SVD. When the data size is huge, the cost of the classical MDS is dominated by the cost of the SVD, whose O(n3) computational complexity makes it infeasible. Therefore, we need a fast MDS method to deal with the huge data problem.The MDS method is useful in the application of dimension reduction. If the matrixX is extended to the higher dimensional space and then shifted and/or rotated by an affine mapping, the double centering result of the corresponding product matrix stays the same. More generally, let X1∈Mn,r1 be derived from X∈Mn,r by an affine mapping such that X1=XQ+1˙bT, where Q is orthogonal and Q∈Mr,r1 and b is a nonzero vector in Rr1 and r≪min(n,r1). That is X1 is a big matrix but with a rank less than r. Assume that D1=X1X1T, the double centering of D1 obtains the same B as double centering of D.Theorem 1.
LetX1 be the dimensional extension from X by an affine mapping X1=XQ+1˙bT. Then the double centering of X1X1T is the same as the double centering of XXT.Proof.
LetD1=X1X1T=(XQ+1˙bT)(XQ+1˙bT)T and let D=XXT, we have
(4)D1=XQQTXT+XQb1˙T+1˙bTQTXT+∥b∥21˙1˙T=XXT+a1˙T+1˙aT+∥b∥21˙1˙T=D+a1˙T+1˙aT+∥b∥21˙1˙T,
where a=XQb is a vector in ℜm. By the definition of Dr, Dc and Dg, we have
(5)D1R=DR+a1˙T+∑air11˙1˙T+∥b∥21˙1˙T,D1C=DC+1˙aT+∑air11˙1˙T+∥b∥21˙1˙T,D1G=DG+2∑air11˙1˙T+∥b∥21˙1˙T,
where ai is the element of a. Then we have,
(6)D1-D1R-D1C+D1G=D-DR-DC+DG.Unlike the previous definition ofD, if D is a distance matrix with each element di,j=(xi-xj)T(xi-xj), the double centering of D2 is equivalent to -2B, provided that ∑i=1n‍xi=0. Hence, the MDS method performs double centering on D2, multiplies by -1/2, and then performs the SVD, which gives the configurations of the data set.


### 1.2. Review of the SCMDS

In 2008, we adapted the classical MDS so as to reduce the originalO(n3) complexity to O(n) [13], in which we have proved that when the data dimension is significantly smaller than the number of data entries, there is a fast linear approach for the classical MDS.The main idea of fast MDS is using statistical resampling to split data into overlapping subsets. We perform the classical MDS on each subset and get the compact Euclidean configuration. Then we use the overlapping information to combine each configuration of subsets to recover the configuration of the whole data. Hence, we named this fast MDS method by the split-and-combine MDS (SCMDS).LetD be the n-by-n product matrix for some X and let ℐ be a random permutation of 1,2,…,n. Assume that ℐk is a nonempty subset of ℐ for k=1,…,K, where ℐk⋂‍ℐk+1≠∅ for k=1,…,K-1 and ⋃‍ℐk=ℐ. ℐk is considered as the index of overlapping subgroups. Dk is the submatrix of D that is extracted from the rows and columns of D by the index ℐk. Then Dk is the product matrix of the submatrix of X by index ℐk. Since the size of Dk is much smaller than D, we can perform the classical MDS easily on each Dk to obtain the MDS result, say Yk. After obtaining the representation coordinates of Yk for each subgroup, we have to combine these representations into one. Without loss of generality, we only demonstrate how to combine the coordinates of the first two subgroups.Assume thatX1 and X2 are matrices whose columns are the two coordinate sets of the overlapping points with respect to Y1 and Y2, respectively. Then there exists an affine mapping that maps X1 to X2. Let X-1 and X-2 be the means of columns of X1 and X2, respectively. In order to obtain the affine mapping, we apply QR factorization to both X1-X-11˙T and X2-X-21˙T. Then we have X1-X-11˙T=Q1R1 and X2-X-21˙T=Q2R2. It is clear that the mean of the center of column vectors of X1-X-11˙T has been shifted to zero. Because X1 and X2 come from the same data set, the difference between X1-X-11˙T and X2-X-21˙T is a rotation. Therefore, the triangular matrices R1 and R2 should be identical when there is no noise in X1 and X2. Due to randomness of the sign of columns of Qi in QR factorization, the sign of columns of Qi might need to be adjusted according to the corresponding diagonal elements of Ri, so that the signs of tridiagonal elements of R1 and R2 are the same.After necessary modification to the sign of columns ofQi, we conclude that
(7)Q1T(X1-X-11˙T)=Q2T(X2-X-21˙T).
Furthermore, we have
(8)X1=Q1Q2TX2-Q1Q2T(X-21˙T)+X-11˙T.
Here, the unitary operator is U=Q1Q2T and the shifting is b=-Q1Q2TX-2+X-1. Since the key step of finding this affine mapping is QR decomposition, the computational cost is O((nI)3), where nI is the number of columns of X1 and X2. Therefore, the cost O((nI)3) complexity is limited by the number of samples in each overlapping region. The proof of the computational cost of SCMDS is given later.Assume that there aren points in a data set. We divide these n samples into K overlapping subgroups, and let nG be the number of points in each subgroup and let  nI be the number of points in each intersection region. Then we have the relationship
(9)KnG-(K-1)nI=n
or
(10)K=(N-nI)(nG-nI).
For each subgroup, we apply the classical MDS to compute the configuration of each group data, which costs O((nG)3). In each overlapping region, we apply QR factorization to compute the affine transformation, which costs O((nI)3). Assume that the true data dimension is r, and the lower bound of nI is r+1. For convenience, we take nG=αp for some constant α>2. Then the total computational cost is about
(11)N-r(α-1)rO(α3r3)+N-αr(α-1)rO(r3)≈O(r2N).
When r≪n, the computational cost O(r2n) is much smaller than O(nn), which is the computation time of the fast MDS method proposed by Morrison et al., 2003 [14]. The key idea of our fast MDS method is to split data into subgroups, then combine the configurations to recover the whole one. Since all the order three complexities are restricted in the small number of data entries, we can therefore speed up MDS.Proposition 2.
When the number of samplesn is much greater than the dimension of data r, the computational complexity of SCMDS is O(r2n), which is linear.


## 1.1. Review of the Classical MDS

Assume thatX is an n-by-r matrix with rank r, where the ith row of X represents the coordinates of the ith data point in an r-dimensional space. Let D=XXT be the product matrix of X. Torgerson [10] proposed the first classical MDS method to reconstruct a matrix X of Cartesian coordinates of points from a given symmetric matrix D. The key of the classical MDS is to apply the double centering operator and the singular value decomposition.Assume that1˙ is an n-by-1 matrix whose elements are all 1’s. We define a symmetric matrix B by
(1)B=(X-1n1˙1˙TX)(X-1n1˙1˙TX)T=XXT-1nXXT1˙1˙T-1n1˙1˙TXXT+1n21˙1˙TXXT1˙1˙T=D-1nD1˙1˙T-1n1˙1˙TD+1n21˙1˙TD1˙1˙T=D-DR-DC+DG,
where DR=(1/n)D1˙1˙T is the row mean matrix of D, DC=(1/n)1˙1˙TD is the column mean matrix of D, and DG=(1/n2)1˙1˙TXXT1˙1˙T is the ground mean matrix of D. The operator from D to D-DR-DC+DG is called double centering. If we define a matrix H by
(2)H=I-1n1˙1˙T,B can be simplified to B=HDH. Since matrix B is symmetric, the SVD decomposes B into B=VΣVT. Then we have
(3)B=VΣ1/2PT=X-1n1˙1˙TX,
for some unitary matrix P. In practice, we set P=I to obtain the MDS result, namely, B.AlthoughB≠X, B=VΣ1/2=(X-(1/n)1˙1˙TX)P for some unitary matrix P, it is the same as shifting all the rows of X such that the row’s mean is zero and then multiplying the result by P. Hence, the row of B preserves the pairwise relationship among the data points. B is derived from D by double centering, whose cost is lower than obtaining B by the SVD. When the data size is huge, the cost of the classical MDS is dominated by the cost of the SVD, whose O(n3) computational complexity makes it infeasible. Therefore, we need a fast MDS method to deal with the huge data problem.The MDS method is useful in the application of dimension reduction. If the matrixX is extended to the higher dimensional space and then shifted and/or rotated by an affine mapping, the double centering result of the corresponding product matrix stays the same. More generally, let X1∈Mn,r1 be derived from X∈Mn,r by an affine mapping such that X1=XQ+1˙bT, where Q is orthogonal and Q∈Mr,r1 and b is a nonzero vector in Rr1 and r≪min(n,r1). That is X1 is a big matrix but with a rank less than r. Assume that D1=X1X1T, the double centering of D1 obtains the same B as double centering of D.Theorem 1.
LetX1 be the dimensional extension from X by an affine mapping X1=XQ+1˙bT. Then the double centering of X1X1T is the same as the double centering of XXT.Proof.
LetD1=X1X1T=(XQ+1˙bT)(XQ+1˙bT)T and let D=XXT, we have
(4)D1=XQQTXT+XQb1˙T+1˙bTQTXT+∥b∥21˙1˙T=XXT+a1˙T+1˙aT+∥b∥21˙1˙T=D+a1˙T+1˙aT+∥b∥21˙1˙T,
where a=XQb is a vector in ℜm. By the definition of Dr, Dc and Dg, we have
(5)D1R=DR+a1˙T+∑air11˙1˙T+∥b∥21˙1˙T,D1C=DC+1˙aT+∑air11˙1˙T+∥b∥21˙1˙T,D1G=DG+2∑air11˙1˙T+∥b∥21˙1˙T,
where ai is the element of a. Then we have,
(6)D1-D1R-D1C+D1G=D-DR-DC+DG.Unlike the previous definition ofD, if D is a distance matrix with each element di,j=(xi-xj)T(xi-xj), the double centering of D2 is equivalent to -2B, provided that ∑i=1n‍xi=0. Hence, the MDS method performs double centering on D2, multiplies by -1/2, and then performs the SVD, which gives the configurations of the data set.


## 1.2. Review of the SCMDS

In 2008, we adapted the classical MDS so as to reduce the originalO(n3) complexity to O(n) [13], in which we have proved that when the data dimension is significantly smaller than the number of data entries, there is a fast linear approach for the classical MDS.The main idea of fast MDS is using statistical resampling to split data into overlapping subsets. We perform the classical MDS on each subset and get the compact Euclidean configuration. Then we use the overlapping information to combine each configuration of subsets to recover the configuration of the whole data. Hence, we named this fast MDS method by the split-and-combine MDS (SCMDS).LetD be the n-by-n product matrix for some X and let ℐ be a random permutation of 1,2,…,n. Assume that ℐk is a nonempty subset of ℐ for k=1,…,K, where ℐk⋂‍ℐk+1≠∅ for k=1,…,K-1 and ⋃‍ℐk=ℐ. ℐk is considered as the index of overlapping subgroups. Dk is the submatrix of D that is extracted from the rows and columns of D by the index ℐk. Then Dk is the product matrix of the submatrix of X by index ℐk. Since the size of Dk is much smaller than D, we can perform the classical MDS easily on each Dk to obtain the MDS result, say Yk. After obtaining the representation coordinates of Yk for each subgroup, we have to combine these representations into one. Without loss of generality, we only demonstrate how to combine the coordinates of the first two subgroups.Assume thatX1 and X2 are matrices whose columns are the two coordinate sets of the overlapping points with respect to Y1 and Y2, respectively. Then there exists an affine mapping that maps X1 to X2. Let X-1 and X-2 be the means of columns of X1 and X2, respectively. In order to obtain the affine mapping, we apply QR factorization to both X1-X-11˙T and X2-X-21˙T. Then we have X1-X-11˙T=Q1R1 and X2-X-21˙T=Q2R2. It is clear that the mean of the center of column vectors of X1-X-11˙T has been shifted to zero. Because X1 and X2 come from the same data set, the difference between X1-X-11˙T and X2-X-21˙T is a rotation. Therefore, the triangular matrices R1 and R2 should be identical when there is no noise in X1 and X2. Due to randomness of the sign of columns of Qi in QR factorization, the sign of columns of Qi might need to be adjusted according to the corresponding diagonal elements of Ri, so that the signs of tridiagonal elements of R1 and R2 are the same.After necessary modification to the sign of columns ofQi, we conclude that
(7)Q1T(X1-X-11˙T)=Q2T(X2-X-21˙T).
Furthermore, we have
(8)X1=Q1Q2TX2-Q1Q2T(X-21˙T)+X-11˙T.
Here, the unitary operator is U=Q1Q2T and the shifting is b=-Q1Q2TX-2+X-1. Since the key step of finding this affine mapping is QR decomposition, the computational cost is O((nI)3), where nI is the number of columns of X1 and X2. Therefore, the cost O((nI)3) complexity is limited by the number of samples in each overlapping region. The proof of the computational cost of SCMDS is given later.Assume that there aren points in a data set. We divide these n samples into K overlapping subgroups, and let nG be the number of points in each subgroup and let  nI be the number of points in each intersection region. Then we have the relationship
(9)KnG-(K-1)nI=n
or
(10)K=(N-nI)(nG-nI).
For each subgroup, we apply the classical MDS to compute the configuration of each group data, which costs O((nG)3). In each overlapping region, we apply QR factorization to compute the affine transformation, which costs O((nI)3). Assume that the true data dimension is r, and the lower bound of nI is r+1. For convenience, we take nG=αp for some constant α>2. Then the total computational cost is about
(11)N-r(α-1)rO(α3r3)+N-αr(α-1)rO(r3)≈O(r2N).
When r≪n, the computational cost O(r2n) is much smaller than O(nn), which is the computation time of the fast MDS method proposed by Morrison et al., 2003 [14]. The key idea of our fast MDS method is to split data into subgroups, then combine the configurations to recover the whole one. Since all the order three complexities are restricted in the small number of data entries, we can therefore speed up MDS.Proposition 2.
When the number of samplesn is much greater than the dimension of data r, the computational complexity of SCMDS is O(r2n), which is linear.


## 2. Methodology

After reviewing the MDS and the SCMDS, We will now demonstrate how to adapt the SCMDS method to become the fast PCA and how the fast PCA can become the fast SVD with further modification.


### 2.1. From the SCMDS to the SCPCA

Because the MDS is similar to the PCA when the data configuration is Euclidean, we can adapt the SCMDS method to obtain the fast PCA with the similar constraint when the rankr is much smaller than the size m and n.Equation (1) shows how to use double centering to convert the product matrix D to B and the B=X-1˙1˙TX. We note that the MDS result B is equal to the matrix obtained from shifting the coordinates of X with the mean of data points to the original point followed by a transformation by some unitary matrix P.LetA∈Mn,r; the purpose of the PCA is to find the eigenvectors of the covariance matrix of A and to obtain the corresponding scores. Since the covariance matrix of A defined by (1/(r-1))(A-A(1/r)1˙1˙T)(A-(1/r)1˙1˙T)T is symmetric and finding the eigenvectors of the covariance matrix of A is similar to obtaining V of B in (1), we set A-A(1/r)1˙1˙T=X-(1/n)1˙1˙TX by omitting the constant 1/r-1. We note that such omission of the constant will not change the direction of the eigenvectors. When we want to carry out the PCA on A, the related matrix X can be obtained by X=A-A(1/r)1˙1˙T+1˙μ, where μ∈Rr is defined by
(12)μ=12n1˙T(A-A1r1˙1˙T).
In other words, carrying out the PCA on A is equivalent to carrying out the MDS on D=XXT. That is
(13)(X-1n1˙1˙TX)(X-1n1˙1˙TX)T=(A-A1r1˙1˙T)(A-1r1˙1˙T)T=VΣVT,
and the result of PCA will be V and (1/r-1)Σ1/2. From the result of the MDS (B=VΣ1/2), we can easily obtain V and Σ by normalizing the column vectors of B.In practice, we do not want to produce the product matrixD=XXT when n is huge. Rather, after we obtain matrix X from A, we will randomly produce the index permutation ℐ and use the index of overlapping subgroup ℐk to compute Dk. And then we follow the step of the SCMDS to obtain the result, B.However, the result of the SCMDS is of the formBQ+1˙bT, where Q is an unitary matrix in Mr and b∈ℜr is usually a nonzero vector. Removal of the factor b can be easily done by computing the average of the row vectors of the SCMDS result. After removing the component b from each row vector, we obtain BQ.SinceBQ∈Mm,r and r≪m, we have
(14)(BQ)TBQ=QTΣ1/2VTVΣ1/2Q=QTΣQ,
which is a r×r small matrix. It is easy to solve this small eigenvalue problem to obtain Q, Σ, and V. The computational cost from BQ+1˙bT to obtain V and Σ1/2 is O(m). Hence, the SCMDS approach to computing the SVD is still linear.Notice that there are some linear SVD methods for thin matrix (m≪n or m≫n). Even in the case of big matrix with small rank, the traditional SVD method, for example the GR-SVD, can be implemented to be linear (because of the dependency, many components of the matrix become zero in the GR-SVD algorithm). However, if the rank of big matrix is almost full and the numerical rank is small or the matrix size is considerably huge that loading the whole matrix is impossible, the SCMDS has the advantage in computing speed. And we call this approach to obtaining V and Σ1/2 by the SCMDS to be SCPCA.


### 2.2. From the SCPCA to the SCSVD

The concepts of the SVD and the PCA are very similar. Since the PCA starts from decomposing the covariance matrix of a data set, it can be considered as adjusting the center of mass of a row vector to zero. On the other hand, the SVD operates directly on the product matrix without shifting. If the mean of the matrix rows is zero, the eigenvectors derived by the SVD are equal to the eigenvectors derived by the PCA. We are looking for a method which will give a fast approach to produce the SVD result without recomputing the eigenvectors of the whole data set, when the PCA result is given. The following is the mathematical analysis for this process.LetX be a column matrix of data set: X~=X-X-·1˙T, where X- is the mean of columns of X. Hence, the row mean of X~ is zero. Assume that we have the PCA result of X, that is, X~X~T=UΣ2UT. Then we have X~=UΣWT for some orthogonal matrix W. Assume that the rank of X~ is r and r is much smaller than the matrix size. We observe that rank(X)=r or rank(X)=r+1, depending on whether X- is spanned by columns of U. If X- is spanned by U, then
(15)X=X~+X-·1˙T=UΣWT+U·c·1˙T=U(ΣWT+c·1˙T),
where c is the coefficient vector of X- when represented by U, that is, X-=U·c.If the singular value decomposition ofΣWT+c·1˙T is U2Σ1V1T, we have
(16)X=U(U2Σ1V1T)=(UU2)Σ1V1T=U1Σ1V1T.
Because the matrix U2 is unitary, U1=UU2 is automatically an orthogonal matrix as well. Then we have the SVD of X.Checking the matrix size ofΣWT+c·1˙T, we can see that to compute the SVD of ΣWT+c·1˙T is not a big task. This is because ΣWT+c·1˙T is a r-by-n matrix, and under our assumption, r is much smaller than n, so we can apply the economic SVD to obtain the decomposition of ΣWT+c·1˙T.On the other hand, ifX- is not spanned by U, the analysis becomes
(17)X=X~+X~·1˙T=[U∣ur+1]([Σ000][VT0]+c·1˙T),
where ur+1 is a unit vector defined by
(18)ur+1=(I-UUT)X-∥(I-UUT)X-∥.
Using the same concept of diagonalization in the case when X- is spanned by columns of U, we find the SVD of
(19)([Σ000][VT0]+c·1˙T)=U2Σ1V1T.
Then X=[U∣ur+1]U2Σ1V1T=U1Σ1V1T, where U1=[U∣ur+1]U2 is another orthogonal matrix, and hence the SVD of X is completed.Proposition 3.
LetA∈Mm,n(ℜ) be a big size matrix and let its rank r be much smaller than min{m,n}. Then the SCSVD or the SCPCA has linear computational complexity.From the above analysis, we can have a fast PCA approach by computing the SCMDS first and then adapt the MDS result to obtain the PCA. We named this approach SCPCA. Similarly, the fast SVD approach, which computes the SCMDS first, then adapts the MDS result to obtain the PCA, and finally adapts the PCA result to the SVD, is called the SCSVD. These two new approaches work when the rank ofX is much smaller than the number of samples and the number of variables. To obtain the exact solution, the parameter nI must be greater than the rank of X. In the SCPCA or SCMDS method, if nI≤r, we only get the approximated solutions of the PCA and the SVD. Under the necessary criterion, we can reduce the computational complexity from min{O(rmn),O(mn2),O(nm2} to min{O(rm),O(rn)}. If the numerical rank is not small, for example, r≈min(m,n), the computational complexity becomes almost the same as the original PCA and SVD. Our approach has no advantage in the latter case.


## 2.1. From the SCMDS to the SCPCA

Because the MDS is similar to the PCA when the data configuration is Euclidean, we can adapt the SCMDS method to obtain the fast PCA with the similar constraint when the rankr is much smaller than the size m and n.Equation (1) shows how to use double centering to convert the product matrix D to B and the B=X-1˙1˙TX. We note that the MDS result B is equal to the matrix obtained from shifting the coordinates of X with the mean of data points to the original point followed by a transformation by some unitary matrix P.LetA∈Mn,r; the purpose of the PCA is to find the eigenvectors of the covariance matrix of A and to obtain the corresponding scores. Since the covariance matrix of A defined by (1/(r-1))(A-A(1/r)1˙1˙T)(A-(1/r)1˙1˙T)T is symmetric and finding the eigenvectors of the covariance matrix of A is similar to obtaining V of B in (1), we set A-A(1/r)1˙1˙T=X-(1/n)1˙1˙TX by omitting the constant 1/r-1. We note that such omission of the constant will not change the direction of the eigenvectors. When we want to carry out the PCA on A, the related matrix X can be obtained by X=A-A(1/r)1˙1˙T+1˙μ, where μ∈Rr is defined by
(12)μ=12n1˙T(A-A1r1˙1˙T).
In other words, carrying out the PCA on A is equivalent to carrying out the MDS on D=XXT. That is
(13)(X-1n1˙1˙TX)(X-1n1˙1˙TX)T=(A-A1r1˙1˙T)(A-1r1˙1˙T)T=VΣVT,
and the result of PCA will be V and (1/r-1)Σ1/2. From the result of the MDS (B=VΣ1/2), we can easily obtain V and Σ by normalizing the column vectors of B.In practice, we do not want to produce the product matrixD=XXT when n is huge. Rather, after we obtain matrix X from A, we will randomly produce the index permutation ℐ and use the index of overlapping subgroup ℐk to compute Dk. And then we follow the step of the SCMDS to obtain the result, B.However, the result of the SCMDS is of the formBQ+1˙bT, where Q is an unitary matrix in Mr and b∈ℜr is usually a nonzero vector. Removal of the factor b can be easily done by computing the average of the row vectors of the SCMDS result. After removing the component b from each row vector, we obtain BQ.SinceBQ∈Mm,r and r≪m, we have
(14)(BQ)TBQ=QTΣ1/2VTVΣ1/2Q=QTΣQ,
which is a r×r small matrix. It is easy to solve this small eigenvalue problem to obtain Q, Σ, and V. The computational cost from BQ+1˙bT to obtain V and Σ1/2 is O(m). Hence, the SCMDS approach to computing the SVD is still linear.Notice that there are some linear SVD methods for thin matrix (m≪n or m≫n). Even in the case of big matrix with small rank, the traditional SVD method, for example the GR-SVD, can be implemented to be linear (because of the dependency, many components of the matrix become zero in the GR-SVD algorithm). However, if the rank of big matrix is almost full and the numerical rank is small or the matrix size is considerably huge that loading the whole matrix is impossible, the SCMDS has the advantage in computing speed. And we call this approach to obtaining V and Σ1/2 by the SCMDS to be SCPCA.


## 2.2. From the SCPCA to the SCSVD

The concepts of the SVD and the PCA are very similar. Since the PCA starts from decomposing the covariance matrix of a data set, it can be considered as adjusting the center of mass of a row vector to zero. On the other hand, the SVD operates directly on the product matrix without shifting. If the mean of the matrix rows is zero, the eigenvectors derived by the SVD are equal to the eigenvectors derived by the PCA. We are looking for a method which will give a fast approach to produce the SVD result without recomputing the eigenvectors of the whole data set, when the PCA result is given. The following is the mathematical analysis for this process.LetX be a column matrix of data set: X~=X-X-·1˙T, where X- is the mean of columns of X. Hence, the row mean of X~ is zero. Assume that we have the PCA result of X, that is, X~X~T=UΣ2UT. Then we have X~=UΣWT for some orthogonal matrix W. Assume that the rank of X~ is r and r is much smaller than the matrix size. We observe that rank(X)=r or rank(X)=r+1, depending on whether X- is spanned by columns of U. If X- is spanned by U, then
(15)X=X~+X-·1˙T=UΣWT+U·c·1˙T=U(ΣWT+c·1˙T),
where c is the coefficient vector of X- when represented by U, that is, X-=U·c.If the singular value decomposition ofΣWT+c·1˙T is U2Σ1V1T, we have
(16)X=U(U2Σ1V1T)=(UU2)Σ1V1T=U1Σ1V1T.
Because the matrix U2 is unitary, U1=UU2 is automatically an orthogonal matrix as well. Then we have the SVD of X.Checking the matrix size ofΣWT+c·1˙T, we can see that to compute the SVD of ΣWT+c·1˙T is not a big task. This is because ΣWT+c·1˙T is a r-by-n matrix, and under our assumption, r is much smaller than n, so we can apply the economic SVD to obtain the decomposition of ΣWT+c·1˙T.On the other hand, ifX- is not spanned by U, the analysis becomes
(17)X=X~+X~·1˙T=[U∣ur+1]([Σ000][VT0]+c·1˙T),
where ur+1 is a unit vector defined by
(18)ur+1=(I-UUT)X-∥(I-UUT)X-∥.
Using the same concept of diagonalization in the case when X- is spanned by columns of U, we find the SVD of
(19)([Σ000][VT0]+c·1˙T)=U2Σ1V1T.
Then X=[U∣ur+1]U2Σ1V1T=U1Σ1V1T, where U1=[U∣ur+1]U2 is another orthogonal matrix, and hence the SVD of X is completed.Proposition 3.
LetA∈Mm,n(ℜ) be a big size matrix and let its rank r be much smaller than min{m,n}. Then the SCSVD or the SCPCA has linear computational complexity.From the above analysis, we can have a fast PCA approach by computing the SCMDS first and then adapt the MDS result to obtain the PCA. We named this approach SCPCA. Similarly, the fast SVD approach, which computes the SCMDS first, then adapts the MDS result to obtain the PCA, and finally adapts the PCA result to the SVD, is called the SCSVD. These two new approaches work when the rank ofX is much smaller than the number of samples and the number of variables. To obtain the exact solution, the parameter nI must be greater than the rank of X. In the SCPCA or SCMDS method, if nI≤r, we only get the approximated solutions of the PCA and the SVD. Under the necessary criterion, we can reduce the computational complexity from min{O(rmn),O(mn2),O(nm2} to min{O(rm),O(rn)}. If the numerical rank is not small, for example, r≈min(m,n), the computational complexity becomes almost the same as the original PCA and SVD. Our approach has no advantage in the latter case.


## 3. SVD for Continuously Growing Data

In this section, we look for the solution when the data is updated constantly and we need to compute the SVD continuously. Instead of scanning all the data again, we try to use the previous result of the SCSVD together with the new updated data to compute the next SVD. Before introducing our updating method, we review the general updating methods first.LetA be an m-by-n matrix, where m is the number of variables and n is the number of samples. And we assume that both m and n are huge. When new data comes in, we collect these new data to form a column matrix which is denoted by E. Assume that we have the singular value decomposition of A, that is,
(20)A=UΣVT,
where U∈Mm,r(ℜ), V∈Mn,r(ℜ) are orthogonal and Σ∈Mr(ℜ) is a diagonal. Since the data gets updated, the data matrix becomes
(21)A1=[A∣E].
To compute the singular value decomposition of A1, we need to compute the eigenvalue and eigenvector of A1A1T.IfE is spanned by the column space of  U, we can represent the column matrix E by E=UC, where C is the coefficient matrix of U with columns of S as the basis. Since U is orthogonal, the coefficient matrix C can be computed easily by C=UTE. Then we have
(22)A1A1T=[A∣UC][A∣UC]T=AAT+UC(UC)T=U(Σ2+CCT)UT=UQΣ12QTUT=U1Σ12U1T.Note that the matrixΣ2+CCT is positive symmetric and it is a small size matrix. Using the spectrum theorem, we can decompose this matrix into QΣ12QT. Because the matrix Q is unitary, U1 is U rotated by Q. In this case, the computational cost is dominated by the rotation from U to U1, and it is O(mr2).In general, the extended matrixE is not spanned by the column space of U and the decomposition of A1A1T should be modified by
(23)A1A1T=[A∣UC+E1][A∣UC+E1]T=AAT+(UC+E1)(UC+E1)T=[U∣Q1][Σ2+CCTCR1TR1CTR1R1T][UTQ1T]=[U∣Q1]Q2Σ12Q2T[UTQ1T]=U1Σ12U1T,
where E=UC+E1, U1=[U∣Q1]Q2 and E1=Q1R1 is the QR decomposition of E1. The cost to obtain C and E1 is O(mkr), the computational cost of QR decomposition of E1 is O(mk2), the cost to obtain Σ12 is O((r+k)3), and the cost of rotation U to U1 is O(m(r+k)2). When m≫r,k the cost of update process is linear O(m). However, if the rank is not small (r≈m or r≈n), the computational cost is dominated by O((r+k)3).Now, we discuss how to update the SVD in the SCSVD approach. If the updated column vectorsE are spanned by the column vectors of some subgroup of the original data, say Xi, the dimension will not increase by the extended matrix E. We just reset the distance matrix Di in this subgroup by
(24)Di=([Xi∣Ei]-1nG,i+k[Xi∣Ei]1˙1˙T)T×([Xi∣Ei]-1ni+k[Xi∣Ei]1˙1˙T),
where nG,i is the number of vectors in the original i-th subgroup. To compute the MDS result of this updated subgroup, it costs O((nG,i+k)3). Then we will combine this nG,i+k MDS representation points with the remaining points. The cost of obtaining the combination affine mapping is still O((nI+k)3) which is increased by the new update. Combining this new subgroup with the whole data set needs O(r2(nG,i-nI)). Then transfering the MDS result to the SVD needs O(mr2), which is linear too.IfE is not spanned by any subgroup of the previous SCMDS process, then we have to update E to the largest group of the subgroups. Because the dimension of this updated group will increase, we have to make some modifications in the combined approach. If the column matrix of the MDS configuration of the new updated group is X1∈Mr1,k1 and the column matrix of the configuration of the other elements is X2∈Mr,k2, where r1>r, then we will extend the matrix X2 to an r1-by-k2 matrix by filling zeros in the bottom of the original X2. After adjusting X1 and X2 to the same dimension of rows, we perform the same combination step as before. Here, the width of matrix X2 is almost the number of whole data, and this combination step takes O(mr12).There is one important difference between the general approach and the SCSVD approach. To obtain the SVD from the MDS result, we need the column mean vector of the matrix. The column mean vector of the original matrix must be computed, and then the column mean vector is updated by the simple weighted average when the new data comes in.Notice that no matter whetherE is spanned by U or not, the update method in the SCSVD has no advantage over the general update method although the computational costs are of the same order. This is because the true computational cost of the SCSVD updating method is always more than the general updating method, even though they are of the same order. Thus, when we have to update the new data to the matrix on which the SVD will be performed, the general method is recommended. However, when the matrix size is so huge that to expand Σ2+CCT becomes impossible, it is better to recompute the SVD result. In this case and with the condition that the true rank is much smaller than the matrix size, the SCSVD approach is recommended.


## 4. Experimental Result

In this section, we show that our fast PCA and fast SVD methods work well for big-sized matrices with small ranks. The simulated matrix is created by the product of two slender matrices and one small diagonal matrix, that is,Am,r×Λr×Br,n, where m,n≫r. The size of the first matrix is m×r, the second matrix is r×r, and the third matrix is r×n. Then the product of these three matrixes is of size m×n and its rank is smaller than r. The absolute values of the diagonal elements of Λ are decreasing to simulate the situation as data normally comes in. When m and n are large and r is much smaller than m and n, the simulated matrix satisfies our SCSVD condition. We pick m=4000, n=4000, and r=50 as our first example. Each element of the simulated matrix is generated from the normal distribution 𝒩(0,1) and then the diagonal terms of Λ are adjusted. The average elapsed time of the SCSVD is 3.98 seconds, while the economical SVD takes 16.14 seconds. Here, each average value is derived from 100 repeats.If we increase the matrix tom=20000, n=20000 and fix the same rank r=50, the elapsed time of the economical SVD is 1209.92 seconds, but the SCSVD is only 195.85 seconds. We observe that our SCSVD method demonstrates significant improvement. Figure 1 shows the speed comparison between the economical SVD (solid line) and the SCSVD (dashed line) with the square matrix size from 500 to 4000 by fixed rank 50. We also use fixed parameter NI=51 and NG=2NI in each simulation test. We can see that the computational cost of the SVD follows the order 3 increase, compared with linear increase of the SCSVD.Figure 1
Comparison of the elapsed time between economical SVD (the solid line) and SCSVD (the dashed line).Note that when the estimated rank used in the SCSVD is greater than the real rank of data matrix, there is almost no error (except rounding error) between the economic SVD and the SCSVD. The error between the economical SVD and the SCSVD is computed by comparing the orthogonality. Assume thatA=UΣVT is derived by the SCSVD and A=U1Σ1V1T is derived from the economical SVD. Since the output of the column vectors of U and V might reverse, the error is defined by
(25)∥U20TU1,20-I20∥,
where U20 is the first 20th column of U, U1,20 is the first 20th column of U1, and I20 is the identity matrix of size 20. The error is about 10-13 only. The average error in the same experiment in Figure 1 is 9.7430×10-13 and the standard deviation is 5.2305×10-13. Thus, when the estimated rank of the SCSVD is greater than the true rank, the accuracy of the SCSVD is pretty much the same as the SVD in the case of a small rank matrix.We would like to explore what happens if the estimated rank is smaller than the true rank. According to the experiment of the SCMDS [13], if the estimated rank is smaller than the true rank, the performance of the SCMDS will deteriorate. In our experiment, we set the true rank of the simulated matrix to be 50 and then observe the estimated rank from 49 to 39. The relationship between the error and the estimated rank with different matrix size is shown in Figure 2.Figure 2
The relationship between errors and estimated dimension of matrix size from 500 to 4000 with step size 500. The true rank of these matrices is 50. From the bottom to top is matrix of sizes 500, 1000, 1500,…, respectively.We can see that when the estimated rank decreases, the error arises rapidly. Lines in Figure2 from the bottom to the top are the matrix size with 500×500 to 4000×4000, respectively. The error increases slowly when the matrix size increases. We observe that making the estimated rank greater than the true dimension is essential for our method.The purpose of the second simulation experiment is to observe the approximation performance of applying the SCPCA to a big full rank matrix. We generate a random matrix with a fixed number of columns and rows, say 1000. The square matrix is created by the form,A1000,r×Λ×Br,1000+αE1000,1000, where r=50 is the essential rank, E is the perturbation, and α is a small coefficient for adjusting the influence to the previous matrix. Such a matrix can be considered as a big-sized matrix with a small rank added by a full rank perturbation matrix. We will show that our method works well for this type of matrices.Figure3 shows the error versus estimated rank, where the error is defined as (25), which is a comparison of the orthogonality between U and U1. All the elements of matrices A, B, E are randomly generated from the normal distribution 𝒩(0,1) and Λ is a diagonal matrix in which the diagonal terms decay by the order, where α=0.01 and the essential rank r=50. We can see that when the estimated rank increases, the composition error decreases. Especially when the estimated rank is greater than the essential rank r, there is almost no error. Thus, it is important to make sure that the estimated rank is greater than the essential rank. In other words, when the estimated rank of the SCSVD is smaller than the essential rank, our SCSVD result can be used as the approximated solution of the SVD.Figure 3
The effect of estimated rank to the error. The matrix size is from 500-by-500 to 4000-by-4000 (from the bottom to the top, resp.) and its essential rank is 50 (α=0.01). When the estimated rank is greater than 50, there is almost no composition error.In the last experimental result, we let the matrix be growing and observe the performance between the general and the SCSVD update approaches. We start from a matrix that is formed byA4000,50×Λ50×B50,4000+0.01·E4000,4000, where Λ is a diagonal matrix such that the diagonal terms are positive and decayed rapidly like natural data and each element in A and B is a random number from the normal distribution 𝒩(0,1). Then we update a 100×4000 random matrix into the matrix 10 times of the size. The element of the updated matrix is also from the normal distribution 𝒩(0,1). When the matrix gets updated, we compute the SVD decomposition to obtain the first 20 columns of U. The updated result of the general approach is U1×Σ1×V1, and the updated result of the SCMDS approach is U2×Σ2×V2. We compute the error as (25) for U1 and U2. The reason that we concern ourselves only with the first 20 bases of the column space is that it is rare to use such a high dimension in the dimension reduction applications.The computational time of the SCMDS is more than the general approach; however, the difference is within one second. The error of the SCMDS approach to update100 rows compared with recomputing the SVD decomposition is 4.0151×10-5 with standard deviation 3.5775×10-5. The error of general approach is 1.4371×10-12 with standard deviation 1.5215×10-12. The update time of the SCMDS and general approach is shown in Figure 4. The red solid line is the time for the SCMDS and the blue line is that for general approach. We can see that the SCMDS approach is about 8 times of general approach. Hence the SCSVD update approach is not recommended for either saving time or controlling error.Figure 4
Compare the update time for SCMDS approach (red line) and the general approach (blue line).


## 5. Conclusion

We proposed the fast PCA and fast SVD methods derived from the technique of the SCMDS method. The new PCA and the SVD have the same accuracy as the traditional PCA and the SVD method when the rank of a matrix is much smaller than the matrix size. The results of applying the SCPCA and the SCSVD to a full rank matrix are also quite reliable when the essential rank of matrix is much smaller than the matrix size. Thus, we can use the SCSVD in a huge data application to obtain a good approximated initial value. The updating algorithm of the SCMDS approach is discussed and compared with the general update approach. The performances of the SCMDS approach both in the computational time and error are worse than the general approach. Hence, the SCSVD method is only recommended for computing the SVD for a large matrix but is not recommended for the update approach.



---

*Source: 683053-2013-03-13.xml*