# SVD by QR and Choleski decomposition - What is going on?

Here's an algorithm I found that performs Singular Value Decomposition. I preferred this algorithm because it can be parallelized, and I don't have to calculate the huge $AA^T$ matrix when the number of rows are very large:

Let A be the matrix for which SVD has to be performed.

$A = U{\Sigma}V^T$

1) First perform $QR$ decomposition of $A$.

$A=QR$

We calculate this by performing the 2nd and the 3rd step.

2) To calculate the $R$ matrix, calculate $A^TA$ and find the Cholesky decomposition of $A^TA$.

$R = Cholesky(A^TA)$

3) Calculate the $Q$ matrix:

$Q = AR^{-1}$

4) Also, perform SVD on the $R$ matrix (a small matrix when compared to A).

$U_R{\Sigma}V^T = R$

The $\Sigma$ and $V$ matrix for $R$ will be the same as for the SVD for $A$. Now calculate the $U$ matrix:

$U = QU_R$

Now, I know the concept of SVD through eigenvectors and eigenvalues of $AA^T$ and $A^TA$, but what is going on here? Is there an intuitive explanation?

• Why do you perform steps 2 and 3? If you have the QR decomposition of $A$, you have already the $Q$ factor you recompute at step 3 (in a very inaccurate way) and $R$ is already the Cholesky factor of $A^TA$. The algorithm would work even if you completely removed the these two steps. Oct 24, 2013 at 14:21
• Oh, I don't calculate the QR decomposition in the first step. What I meant was that I perform QR decomposition by doing step 2 and 3. I'll edit the question accordingly. Oct 24, 2013 at 17:15
• Well, that's not a very good way how to compute the QR factorisation. It is even worse than the classical implementation of the Gram-Schmidt which can be pretty well parallelized. Oct 24, 2013 at 20:21
• The problem is that I'm doing this on MapReduce, and to calculate QR using Gram-Schmidt would require me to calculate the orthogonal vectors of Q iteratively and then the R matrix, which would take another MapReduce job. With the above method, I can calculate $A^TA$ in one go, then calculate R, S and V locally and then calculate U in the next MapReduce job. Anyway, I validated the output of the above algorithm with standard implementations in R and the outputs were matching. How are you saying that it is worse than anything if it gives me the required output AND I'm not worried about its time? Oct 24, 2013 at 21:19
• Yes, taking A=rand(100)^3 (in Matlab notation). Minimal singular value computed by: 1) CGS - 2 valid digits, Cholesky of $A^TA$ - 0 valid digits. Oct 24, 2013 at 22:15

The reason calculating R by the Cholesky factorization is worse than other methods is because it works on sums of squared and crossproducts of floating point numbers, which has been known since the 1980's to be numerically less accurate than the Householder and Givens QR algorithms. In the 1990's during a short course on PLS, Sevante Wold, explained that many multivariate analysis software packages used the Householder algorithm to calculate PC's (for PCA) because it is the numerically most stable of the available algorithms.The Householder algorithm uses planar reflections and the Givens algorithm uses planar rotations. The advantage of the Cholesky algorithm is memory usage which is much less than the other algorithms. However, this is less of an issue now, but was a serious concern in the 1960s and into the early 1990s.

I have generated random matrices for which the above mentioned algorithm does not give the same results as calculating the SVD of A directly from A. This is likely to be true when A contains large values, which when squared, have more digits than can be represented by the the computer being used for the calculations. (Such large numbers are 'rounded' or represented by the nearest floating point number.) Keep in mind that intermediate steps in many data analyses will often generate numbers which push the limits of floating point number systems and that floating point number systems are not the same as the theoretical number systems we learned about in school.

One nice feature of the Givens algorithm (pointed out to me by a friend who worked for SAS many years ago) when used in least-squares problems is that one can add more rows to the end of A and merely spend the time to update Q and R by continuing to do more planar rotations. Also, the Givens algorithms yields what are called least-squares residuals which are statistically independent/uncorrelated. For some situations, having independent residuals is very important and the Givens algorithm calculates them. However, this independence comes at a cost: if A has p columns, then p of the least-squares residuals will be 0.