Solved – Confirming an understanding of SVD

matrix decompositionsvd

I'm trying to sift through the concepts of low-rank approximations, matrix factorizations and SVD. There's a lot of info out there and the rabbit hole is deep so I just want to make sure my high level understanding is accurate:

A common application of SVD is to make low-rank approximations to a matrix, $A$
A low-rank approximation is akin to finding an approximation that retains the minimum number of rows (or columns) that are linearly independent (the rank of $A$). Loosely speaking, we think of these remaining rows as a minimum bound on the useful information contained in $A$. So finding the minimum number is good for making matrices smaller but useful.
One way to do this is to find another matrix, $B$, such that we minimize the Frobenius Norm, $\|A – B\|$
Matrix factorization is an application of low-rank approximation, wherein we can find some $RQ^{T} = B$
SVD is a common matrix factorization algorithm of the form $A = U \Sigma V^{T}$
In addition to factorization, SVD is used to solve equations of the form $Ax=b$. In this regard, under certain conditions, SVD is preferred to QR or Cholesky decomposition.
SVD is a general case (for $m \times n$ matrices) of an eigen-decomposition
The best way (or a very good way??) to find a low-rank approximation to $A$ is by performing SVD
SVD is very popular because it exposes eigenvalues and eigenvectors, which have a plethora of useful properties used in Machine Learning/Statistical applications, making it preferential to QR or Cholesky for factorization.
Using SVD, we can re-write $A$ as a linear combination of its singular values (from the diagonal of $\Sigma$) and the left and right singular vectors, i.e. the columns of $U$ and $V^{T}$
Calculating the SVD consists of finding the eigenvalues and eigenvectors of $AA^{T}$ ($U$) and $A^{T}A$ ($V$).
Calculating this product is very expensive, so other algorithms (like truncated SVD, randomized SVD) are being explored

Are the above points accurate? I understand there is a lot more to these subjects, but I hope to have captured their relationship at a high level. Please let me know if I'm missing something painfully obvious.

Best Answer

I have tried to dispel some of the popular myths regarding the SVD.

``A common application of SVD is to make low-rank approximations to a matrix, $A$''

This is one of the applications of the SVD but the SVD can do very many mathematical operations: Solve linear equations (including least-squares problems and minimum norm solutions), compute inverses (including the Moore-Penrose Pseudo inverse), compute null spaces, compute the numerical rank of a matrix, data compression, image compression. The SVD is the numerical backbone of the PCA (Principal Component Analysis). When I googled "SVD" today, I got 16.4 million hits (I do not know how many for this decomposition) and "Singular Value Decomposition" gave me 2.5 million. So there must be many applications.
``A low-rank approximation is akin to finding an approximation that retains the minimum number of rows (or columns) that are linearly independent (the rank of $A$). Loosely speaking, we think of these remaining rows as a minimum bound on the useful information contained in $A$. So finding the minimum number is good for making matrices smaller but useful.''

It depends on the definition you use to determine the numerical rank of the $m \times n$ matrix $A$. Usually, we compute the ratio $\sigma_{k}/\sigma_1$, where $\sigma_k$ is the $k$th singular value. If this ratio is less than some pre-determined tolerance then we can set the small singular values to zero. Any ratio less than the machine precision $\epsilon$ ( $2^{-53}$ in 8-byte variables ) is noise. A modest linear function of the form $f(m,n) \epsilon$ may be used for the tolerance.

A low rank approximation is of the form $$ \hat{A} = \hat{U} \hat{\Sigma} \hat{V}^t $$ where $\hat{U}$ is the first $k$ columns of $U$ and $\hat{\Sigma}$ is the principal $k \times k$ submatrix of $\Sigma$.

Note that the dimensions of $\hat{A}$ is same as $A$.

``One way to do this is to find another matrix, $B$, such that we minimise the Frobenius Norm, $\|A-B\|_{F}$.''

This is true for the Frobenius norm as well as the spectral norm, $\|A-B\|_2$, where the spectral norm is defined as the largest singular value of the matrix.

``Matrix factorisation is an application of low-rank approximation, wherein we can find some $RQ^t=B$.''

I do not understand this statement.

``SVD is a common matrix factorisation algorithm of the form $A=U\Sigma V^t$.''

The SVD is the most useful decomposition available for a rectangular matrix. It is also very useful for square matrices but there are other useful factorisations for square matrices.

``In addition to factorisation, SVD is used to solve equations of the form $Ax=b$. In this regard, under certain conditions, SVD is preferred to QR or Cholesky decomposition.''

If $A$ is a symmetric positive definite matrix, Cholesky factorisation is the preferable (inexpensive) method for linear solvers. If $A$ is a non-singular square matrix then the LU decomposition (Gaussian elimination with pivoting) can be used. When there are rank deficiencies or if the matrix is rectangular then the QR decomposition and the SVD are preferable. These cases include least-squares and minimum norm solutions. For difficult problems, the SVD is usually the best choice but there are experts who can use QR and LU decompositions in ingenious ways.

``SVD is a general case (for $m \times n$ matrices) of an eigen-decomposition.''

The SVD is based on two eigenvalue problems, $$ A^tA v_i = \sigma_i^2 v_i ,\;\;\; AA^t u_i = \sigma_i^2 u_i $$ where $u_i$ and the $v_i$ are the $i$th column of the singular vector matrices $U$ and $V$, respectively. However, numerical analysts very rarely compute $A^tA$ and $AA^t$. These two eigenvalue problems are solved implicitly without forming these numerically damaging products.

Note that $$ A v_i = \sigma_i u_i, \;\;\; A^tu_i = \sigma_i v_i $$ are the two primary singular vector expressions for the SVD. These two equations can be amalgamated to obtain an eigenvalue problem $$ \left[ \begin{array}{cc} 0 & A \\ A^t & 0 \end{array} \right] \left[ \begin{array}{c} u_i \\ v_i \end{array} \right] = \sigma_i \left[ \begin{array}{c} u_i \\ v_i \end{array} \right] . $$

The best way (or a very good way??) to find a low-rank approximation to $A$ is by performing SVD

Yes.

``SVD is very popular because it exposes eigenvalues and eigenvectors, which have a plethora of useful properties used in Machine Learning/Statistical applications, making it preferential to QR or Cholesky for factorisation.''

In statistics, $A^tA$ (or $AA^t$) are sample covariance matrices if $A$ is a data matrix (with the means removed).

``Using SVD, we can re-write $A$ as a linear combination of its singular values (from the diagonal of $\Sigma$) and the left and right singular vectors, i.e. the columns of $U$ and $V^t$''

The matrix may be expressed as a sum of rank one outer products $$ A = \Sigma_{i=1}^{r} \sigma_i u_i v_i^t, \;\;\; r=\min(m,n) $$ For reduced rank approximations, use $ r < \min(m,n)$.

Calculating the SVD consists of finding the eigenvalues and eigenvectors of $AA^t(U)$ and $A^tA (V)$.

No (as explained in 7.)

Calculating this product is very expensive, so other algorithms (like truncated SVD, randomised SVD) are being explored

It is not expensive unless the matrix is extremely large. There are different algorithms (Lanczos methods) for very large (sparse) matrices.

For $m \ge n$, (if not transpose the matrix), the nominal cost is $O(mn^2)$ as measured by multiplication counts. However, in modern computing architectures, some of the arithmetic can be done in parallel. Furthermore, the cache management using level 3 BLAS (Basic Linear Algebra Subroutines) increases the efficiency. Thus, $O(mn^2)$ complexity is often a crude measure.

The first part of the SVD algorithm is to reduce the matrix to an $n \times n$ triangular matrix using the QR decomposition and the second part is to transform that to an $n \times n$ bidiagonal matrix. The first two steps are not iterative algorithms; they are straight forward steps which use orthogonal transformations. The iterative part is finding the SVD of the bidiagonal matrix which is performed using the QR-algorithm (not to be confused with the QR decomposition) or a divide and conquer algorithm. If only few singular values/vectors are required, then usually a bisection algorithm is used for the required singular values followed by inverse iteration to get the corresponding singular vectors.

Most of the details can be found in "Matrix Computations" by Gene Golub and Charles Van Loan.

Related Solutions

SVD – Singular Value Decomposition of a Three-Dimensional Array Explained

There are several notions of decomposition of such a tensor. Last year I asked essentially the same question on the MaplePrimes site, answered it myself by referring to wikipedia, and provided an implementation for one of those notions (the CANDECOMP/PARAFAC decomposition) in a follow-up post (applied to decomposing the $3\times m \times n$ tensor given by the R,G,B entries of an image).

Solved – Finding matrix eigenvectors using QR decomposition

The (basic) algorithm with QR decomposition is as follows.

Let $X$ by a symmetric matrix.
Let $X_1 = X$, and iterate the following:
Given $X_k$, write a QR decomposition $X_k = Q_k R_k$, and let $X_{k+1} = R_k Q_k$;
The matrices sequence $X_n$ converges to some diagonal matrix $D$ with the eigenvalues on the diagonal; you retrieve the corresponding eigenvectors as the columns of $\prod_i Q_i$.

Here is an example code in R.

# some symmetric matrix
A <- matrix( sample(1:30,16), ncol=4)
A <- A + t(A);

# initialize
X <- A;
pQ <- diag(1, dim(A)[1]);

# iterate 
for(i in 1:30)
{
  d <- qr(X);
  Q <- qr.Q(d);
  pQ <- pQ %*% Q;
  X <- qr.R(d) %*% Q;
}

Now we have a look on the result

> A
     [,1] [,2] [,3] [,4]
[1,]   52   30   49   28
[2,]   30   50    8   44
[3,]   49    8   46   16
[4,]   28   44   16   22

The matrix X contains the eigenvalues on the diagonal:

> round(X,5)
         [,1]    [,2]      [,3]     [,4]
[1,] 132.6279  0.0000   0.00000  0.00000
[2,]   0.0000 52.4423   0.00000  0.00000
[3,]   0.0000  0.0000 -11.54113  0.00000
[4,]   0.0000  0.0000   0.00000 -3.52904

And the product of all Q contains the eigenvectors:

> round(pQ,5)
        [,1]     [,2]     [,3]     [,4]
[1,] 0.60946 -0.29992 -0.09988 -0.72707
[2,] 0.48785  0.65200  0.57725  0.06069
[3,] 0.46658 -0.60196  0.22156  0.60898
[4,] 0.41577  0.35013 -0.77956  0.31117

We can compare to the result of eigen(A) :

> eigen(A)
$values
[1] 132.627875  52.442300  -3.529045 -11.541131

$vectors
           [,1]       [,2]        [,3]        [,4]
[1,] -0.6094595 -0.2999194  0.72707077  0.09987744
[2,] -0.4878528  0.6519967 -0.06068999 -0.57724915
[3,] -0.4665778 -0.6019623 -0.60897966 -0.22156327
[4,] -0.4157690  0.3501285 -0.31117293  0.77956246

Of there is room for lots of improvements, but basically here it is. I once read lots of papers on the subject but my memory is leacking :(

Note that, as your problem is to perform PCA, you will find easily many PCA programs on the internet, you may prefer to do than rather than program it yourself.

Best Answer

Related Solutions

SVD – Singular Value Decomposition of a Three-Dimensional Array Explained

Solved – Finding matrix eigenvectors using QR decomposition

Related Question