Relationship between SVD and PCA. How to use SVD to perform PCA?

Question

Principal component analysis (PCA) is usually explained via an eigen-decomposition of the covariance matrix. However, it can also be performed via singular value decomposition (SVD) of the data matrix $\mathbf X$. How does it work? What is the connection between these two approaches? What is the relationship between SVD and PCA?

Or in other words, how to use SVD of the data matrix to perform dimensionality reduction?

Answer 1

Let the real values data matrix $\mathbf X$ be of $n \times p$ size, where $n$ is the number of samples and $p$ is the number of variables. Let us assume that it is centered, i.e. column means have been subtracted and are now equal to zero.

Then the $p \times p$ covariance matrix $\mathbf C$ is given by $\mathbf C = \mathbf X^\top \mathbf X/(n-1)$. It is a symmetric matrix and so it can be diagonalized: $$\mathbf C = \mathbf V \mathbf L \mathbf V^\top,$$ where $\mathbf V$ is a matrix of eigenvectors (each column is an eigenvector) and $\mathbf L$ is a diagonal matrix with eigenvalues $\lambda_i$ in the decreasing order on the diagonal. The eigenvectors are called principal axes or principal directions of the data. Projections of the data on the principal axes are called principal components, also known as PC scores; these can be seen as new, transformed, variables. The $j$-th principal component is given by $j$-th column of $\mathbf {XV}$. The coordinates of the $i$-th data point in the new PC space are given by the $i$-th row of $\mathbf{XV}$.

If we now perform singular value decomposition of $\mathbf X$, we obtain a decomposition $$\mathbf X = \mathbf U \mathbf S \mathbf V^\top,$$ where $\mathbf U$ is a unitary matrix (with columns called left singular vectors), $\mathbf S$ is the diagonal matrix of singular values $s_i$ and $\mathbf V$ columns are called right singular vectors. From here one can easily see that $$\mathbf C = \mathbf V \mathbf S \mathbf U^\top \mathbf U \mathbf S \mathbf V^\top /(n-1) = \mathbf V \frac{\mathbf S^2}{n-1}\mathbf V^\top,$$ meaning that right singular vectors $\mathbf V$ are principal directions (eigenvectors) and that singular values are related to the eigenvalues of covariance matrix via $\lambda_i = s_i^2/(n-1)$. Principal components are given by $\mathbf X \mathbf V = \mathbf U \mathbf S \mathbf V^\top \mathbf V = \mathbf U \mathbf S$.

To summarize:

1. If $\mathbf X = \mathbf U \mathbf S \mathbf V^\top$, then the columns of $\mathbf V$ are principal directions/axes (eigenvectors).
2. Columns of $\mathbf {US}$ are principal components ("scores").
3. Singular values are related to the eigenvalues of covariance matrix via $\lambda_i = s_i^2/(n-1)$. Eigenvalues $\lambda_i$ show variances of the respective PCs.
4. Standardized scores are given by columns of $\sqrt{n-1}\mathbf U$ and loadings are given by columns of $\mathbf V \mathbf S/\sqrt{n-1}$. See e.g. here and here for why "loadings" should not be confused with principal directions.
5. The above is correct only if $\mathbf X$ is centered. Only then is covariance matrix equal to $\mathbf X^\top \mathbf X/(n-1)$.
6. The above is correct only for $\mathbf X$ having samples in rows and variables in columns. If variables are in rows and samples in columns, then $\mathbf U$ and $\mathbf V$ exchange interpretations.
7. If one wants to perform PCA on a correlation matrix (instead of a covariance matrix), then columns of $\mathbf X$ should not only be centered, but standardized as well, i.e. divided by their standard deviations.
8. To reduce the dimensionality of the data from $p$ to $k<p$, select $k$ first columns of $\mathbf U$, and $k\times k$ upper-left part of $\mathbf S$. Their product $\mathbf U_k \mathbf S_k$ is the required $n \times k$ matrix containing first $k$ PCs.
9. Further multiplying the first $k$ PCs by the corresponding principal axes $\mathbf V_k^\top$ yields $\mathbf X_k = \mathbf U_k^\vphantom \top \mathbf S_k^\vphantom \top \mathbf V_k^\top$ matrix that has the original $n \times p$ size but is of lower rank (of rank $k$). This matrix $\mathbf X_k$ provides a reconstruction of the original data from the first $k$ PCs. It has the lowest possible reconstruction error, see my answer here.
10. Strictly speaking, $\mathbf U$ is of $n\times n$ size and $\mathbf V$ is of $p \times p$ size. However, if $n>p$ then the last $n-p$ columns of $\mathbf U$ are arbitrary (and corresponding rows of $\mathbf S$ are constant zero); one should therefore use an economy size (or thin) SVD that returns $\mathbf U$ of $n\times p$ size, dropping the useless columns. For large $n\gg p$ the matrix $\mathbf U$ would otherwise be unnecessarily huge. The same applies for an opposite situation of $n\ll p$.

---

Further links

• What is the intuitive relationship between SVD and PCA -- a very popular and very similar thread on math.SE.

• Why PCA of data by means of SVD of the data? -- a discussion of what are the benefits of performing PCA via SVD [short answer: numerical stability].

• PCA and Correspondence analysis in their relation to Biplot -- PCA in the context of some congeneric techniques, all based on SVD.

• Is there any advantage of SVD over PCA? -- a question asking if there any benefits in using SVD instead of PCA [short answer: ill-posed question].

• Making sense of principal component analysis, eigenvectors & eigenvalues -- my answer giving a non-technical explanation of PCA. To draw attention, I reproduce one figure here:

Answer 2

I wrote a Python & Numpy snippet that accompanies @amoeba's answer and I leave it here in case it is useful for someone. The comments are mostly taken from @amoeba's answer.

import numpy as np
from numpy import linalg as la
np.random.seed(42)


def flip_signs(A, B):
    """
    utility function for resolving the sign ambiguity in SVD
    http://stats.stackexchange.com/q/34396/115202
    """
    signs = np.sign(A) * np.sign(B)
    return A, B * signs


# Let the data matrix X be of n x p size,
# where n is the number of samples and p is the number of variables
n, p = 5, 3
X = np.random.rand(n, p)
# Let us assume that it is centered
X -= np.mean(X, axis=0)

# the p x p covariance matrix
C = np.cov(X, rowvar=False)
print("C = \n", C)
# C is a symmetric matrix and so it can be diagonalized:
l, principal_axes = la.eig(C)
# sort results wrt. eigenvalues
idx = l.argsort()[::-1]
l, principal_axes = l[idx], principal_axes[:, idx]
# the eigenvalues in decreasing order
print("l = \n", l)
# a matrix of eigenvectors (each column is an eigenvector)
print("V = \n", principal_axes)
# projections of X on the principal axes are called principal components
principal_components = X.dot(principal_axes)
print("Y = \n", principal_components)

# we now perform singular value decomposition of X
# "economy size" (or "thin") SVD
U, s, Vt = la.svd(X, full_matrices=False)
V = Vt.T
S = np.diag(s)

# 1) then columns of V are principal directions/axes.
assert np.allclose(*flip_signs(V, principal_axes))

# 2) columns of US are principal components
assert np.allclose(*flip_signs(U.dot(S), principal_components))

# 3) singular values are related to the eigenvalues of covariance matrix
assert np.allclose((s ** 2) / (n - 1), l)

# 8) dimensionality reduction
k = 2
PC_k = principal_components[:, 0:k]
US_k = U[:, 0:k].dot(S[0:k, 0:k])
assert np.allclose(*flip_signs(PC_k, US_k))

# 10) we used "economy size" (or "thin") SVD
assert U.shape == (n, p)
assert S.shape == (p, p)
assert V.shape == (p, p)

Answer 3

Let me start with PCA. Suppose that you have $n$ data points comprised of $d$ numbers (or dimensions) each. If you center this data (subtract the mean data point $\mu$ from each data vector $x_i$) you can stack the data to make a matrix

$$ X = \left( \begin{array}{ccccc} && x_1^T - \mu^T && \\ \hline && x_2^T - \mu^T && \\ \hline && \vdots && \\ \hline && x_n^T - \mu^T && \end{array} \right)\,. $$

The covariance matrix

$$ S = \frac{1}{n-1} \sum_{i=1}^n (x_i-\mu)(x_i-\mu)^T = \frac{1}{n-1} X^T X $$

measures to which degree the different coordinates in which your data is given vary together. So, it's maybe not surprising that PCA -- which is designed to capture the variation of your data -- can be given in terms of the covariance matrix. In particular, the eigenvalue decomposition of $S$ turns out to be

$$ S = V \Lambda V^T = \sum_{i = 1}^r \lambda_i v_i v_i^T \,, $$

where $v_i$ is the $i$-th Principal Component, or PC, and $\lambda_i$ is the $i$-th eigenvalue of $S$ and is also equal to the variance of the data along the $i$-th PC. This decomposition comes from a general theorem in linear algebra, and some work does have to be done to motivate the relation to PCA.

SVD is a general way to understand a matrix in terms of its column-space and row-space. (It's a way to rewrite any matrix in terms of other matrices with an intuitive relation to the row and column space.) For example, for the matrix $A = \left( \begin{array}{cc}1&2\\0&1\end{array} \right)$ we can find directions $u_i$ and $v_i$ in the domain and range so that

You can find these by considering how $A$ as a linear transformation morphs a unit sphere $\mathbb S$ in its domain to an ellipse: the principal semi-axes of the ellipse align with the $u_i$ and the $v_i$ are their preimages.

In any case, for the data matrix $X$ above (really, just set $A = X$), SVD lets us write

$$ X = \sum_{i=1}^r \sigma_i u_i v_i^T\,, $$

where $\{ u_i \}$ and $\{ v_i \}$ are orthonormal sets of vectors. A comparison with the eigenvalue decomposition of $S$ reveals that the "right singular vectors" $v_i$ are equal to the PCs, the "right singular vectors" are

$$ u_i = \frac{1}{\sqrt{(n-1)\lambda_i}} Xv_i\,, $$

and the "singular values" $\sigma_i$ are related to the data matrix via

$$ \sigma_i^2 = (n-1) \lambda_i\,. $$

It's a general fact that the left singular vectors $u_i$ span the column space of $X$. In this specific case, $u_i$ gives us a scaled projection of the data $X$ onto the direction of the $i$-th principal component. The right singular vectors $v_i$ in general span the row space of $X$, which gives us a set of orthonormal vectors that spans the data much like PCs.

I go into some more details and benefits of the relationship between PCA and SVD in this longer article.

Answer 4

Let's try to understand using a data matrix $X$ of dimension $n \times d$, where $d \gg n$ and $rank(X)=n$.

Then $\underset{n \times d}{X}=\underset{n \times n}{U}\underset{n \times n}{\Sigma} \underset{n \times d}{V^T}$ (reduced SVD)

with $X^TX=V\Sigma^TU^TU\Sigma V^T=V\Sigma^2V^T$ (since unitary / orthonormal $U$, $V$ and diagonal $\Sigma$)

and the covariance matrix (assuming $X$ is already mean centered, i.e., the columns of $X$ have $0$ means)

$C=E[X^TX]-E[X]^TE[X]=\frac{X^TX}{n-1}-0=V\frac{\Sigma^2}{n-1}V^T=\tilde{V}\Lambda \tilde{V}^T$ (PCA, by spectral decomposition)

$\implies \Lambda = \frac{\Sigma^2}{n-1}$ and $V=\tilde{V}$ upto sign flip.

Let's validate the above with eigenfaces (i.e., the principal components / eigenvectors of the covariance matrix for such a face dataset) using the following face dataset:

import numpy as np
from sklearn.decomposition import PCA
from sklearn.datasets import fetch_olivetti_faces

X = fetch_olivetti_faces().data
X.shape # 400 face images of size 64×64 flattened
# (400,4096)
n = len(X)

# z-score normalize
X = X - np.mean(X, axis=0) # mean-centering
# X = X / np.std(X, axis=0)  # scaling to have sd=1

# choose first k eigenvalues / eigenvectors for dimensionality reduction
k = 25

# SVD
U, Σ, Vt = np.linalg.svd(X, full_matrices=False) 

# PCA
pca = PCA(k).fit(X)
PC = pca.components_.T
#Vt.shape, PC.shape

Now let's compare the eigenvalues and eigenvectors computed:

# first k eigenvalues of Λ = Σ^2/(n-1) # from SVD
print(Σ[:k]**2/(n-1))              # from SVD
# [18.840178   11.071763    6.304614    3.9545844   2.8560426   2.49771
# 1.9200633   1.611159    1.5492224   1.3229507   1.2621089   1.1369102
# 0.98639774  0.90758985  0.84092826  0.77355367  0.7271429   0.64526594
# 0.59645116  0.5910001   0.55270135  0.48628208  0.4619924   0.45075357
# 0.4321357 ]
print(pca.explained_variance_[:k]) # from PCA
# [18.840164   11.07176     6.3046117   3.9545813   2.8560433   2.4977121
# 1.9200654   1.6111585   1.549223    1.3229507   1.2621082   1.1369106
# 0.98639697  0.9075892   0.84092826  0.773553    0.72714305  0.64526534
# 0.59645087  0.5909973   0.55269724  0.4862703   0.461944    0.45075053 0.43211046]

# plot PC the k dominant principal components / eigenvectors as images 
# 1. using PC obtained with PCA
# 2. using Vt[:k,:].T obtained from SVD

Here the differences in the eigenvectors are due to sign ambiguity (refer to https://www.osti.gov/servlets/purl/920802)

Answer 5

[This is just to mention how matrix ${ U }$ in the SVD of a data matrix ${ X }$ (also determined by eigenvectors of ${ X X ^{T} }$) gives the principal directions in PCA.]

Consider data ${ x _1, \ldots, x _n \in \mathbb{R} ^{d} }.$

Goal is to find, for every ${ k ,}$ an affine subspace ${ \mathscr{A} _k }$ of dimension ${ k }$ minimising the mean squared loss/error $${ \frac{1}{n} \sum _{i=1} ^{n} \text{dist}( x _i, \mathscr{A} _k) ^2 }.$$

Best 0 dimensional fit:

We want a ${ p \in \mathbb{R} ^{d} }$ such that the loss ${ \frac{1}{n} \sum _{i=1} ^{n} \lVert x _i - p \rVert ^2 }$ is minimised.
Note loss $${\begin{align*} \frac{1}{n} \sum _{i=1} ^{n} \lVert x _i - p \rVert ^2 &= \frac{1}{n} \sum _{i=1} ^{n} \lVert (x _i - q) - (p - q) \rVert ^2 \\ &= \frac{1}{n} \sum _{i=1} ^{n} \lVert x _i - q \rVert ^2 - \frac{1}{n} \sum _{i=1} ^{n} 2 (p-q) ^{T} (x _i - q) + \lVert p - q \rVert ^2 .\end{align*} }$$
Setting ${ q = \overline{x} = \frac{1}{n} \sum _{i=1} ^{n} x _i }$ gives the inequality $${ \begin{align*} \frac{1}{n} \sum _{i=1} ^{n} \lVert x _i - p \rVert ^2 &= \frac{1}{n} \sum _{i=1} ^{n} \lVert x _i - \overline{x} \rVert ^2 + \lVert p - \overline{x} \rVert ^2 \\ &\geq \frac{1}{n} \sum _{i=1} ^{n} \lVert x _i - \overline{x} \rVert ^2 , \end{align*} }$$ with equality holding when ${ p = \overline{x} }.$
Hence the loss ${ \frac{1}{n} \sum _{i=1} ^{n} \lVert x _i - p \rVert ^2 }$ is minimised when ${ p = \overline{x} = \frac{1}{n} \sum _{i=1} ^{n} x _i }.$

Best k dimensional fit:

We want orthonormal vectors ${ [\hat{w} _1, \ldots, \hat{w} _k ] }$ in ${ \mathbb{R} ^d }$ so that ${ \overline{x} + \text{span}(\hat{w} _1, \ldots, \hat{w} _k) }$ is minimiser of the loss $${ \frac{1}{n} \sum _{i=1} ^{n} \text{dist}( x _i, \mathscr{A} _k) ^2 }$$ taken over all ${ k}$ dimensional affine subspaces ${ \mathscr{A} _{k} = \overline{x} + \text{span}(w _1, \ldots, w _k) }.$

Firstly, for any orthonormal vectors ${ W = [w _1, \ldots, w _k] }$ in ${ \mathbb{R} ^d }$ we have $${\begin{align*} & \text{dist}(x _i, \overline{x} + \text{span}(w _1, \ldots, w _k)) ^2 \\ &= \text{dist}(x _i - \overline{x}, \text{span}(w _1, \ldots, w _k)) ^2 \\ &= \lVert x _i - \overline{x} \rVert ^2 - \sum _{j=1} ^{k} ( w _j ^{T} (x _i - \overline{x}) ) ^2 \\ &= \lVert x _i - \overline{x} \rVert ^2 - \sum _{j=1} ^{k} w _j ^{T} \left( (x _i - \overline{x}) (x _i - \overline{x}) ^{T} \right) w _j \end{align*} }$$ and hence the loss $${ \begin{align*} &\frac{1}{n} \sum _{i=1} ^{n} \text{dist}(x _i, \overline{x} + \text{span}(w _1, \ldots, w _k)) ^2 \\ &= \frac{1}{n} \sum _{i=1} ^{n} \lVert x _i - \overline{x} \rVert ^2 - \sum _{j=1} ^{k} w _j ^{T} \left( {\color{red}{\frac{1}{n} \sum _{i=1} ^{n} (x _i - \overline{x}) (x _i - \overline{x}) ^{T}}} \right) w _j. \end{align*} }$$

Writing ${ X = [x _1 - \overline{x}, \ldots, x _n - \overline{x}] \in \mathbb{R} ^{d \times n} }$ we have the matrix $${ \frac{1}{n} \sum _{i=1} ^{n} (x _i - \overline{x}) (x _i - \overline{x}) ^{T} = \frac{1}{n} X X ^{T} }.$$
The problem of minimising the loss is reduced to the problem of finding orthonormal ${ [w _1, \ldots, w _k] \in \mathbb{R} ^{d \times k} }$ which maximise $${ \sum _{j=1} ^{k} w _j ^{T} \left( \frac{1}{n} X X ^{T} \right) w _j .}$$

Writing the SVD of ${ X }$ as ${ X = U \Sigma V ^{T} },$ we have $${ \begin{align*} \frac{1}{n} X X ^{T} &= \frac{1}{n} U \Sigma \Sigma ^T U ^{T} \\ &=\frac{1}{n} U \underbrace{\begin{pmatrix} \sigma _1 ^{2} & & & \\ &\ddots & & \\ & &\sigma _r ^2 & \\ & & & \end{pmatrix}} _{d \times d} U ^{T} \\ &= \frac{1}{n} U \Sigma _d ^2 U ^{T} . \end{align*} }$$

Now for any orthonormal ${ [w _1, \ldots, w _k] },$ we have $${ \begin{align*} &\sum _{j=1} ^{k} w _j ^{T} \left( \frac{1}{n} X X ^{T} \right) w _j \\ &= \frac{1}{n} \sum _{j=1} ^{k} w _j ^{T} U \Sigma _d ^2 U ^{T} w _j \\ &= \frac{1}{n} \sum _{j=1} ^{k} \tilde{w} _j ^{T} \Sigma _d ^2 \tilde{w} _j \end{align*} }$$ where orthonormal set ${ [\tilde{w} _1, \ldots, \tilde{w} _k] = U ^{T} [w _1, \ldots, w _k] }.$

For ${ k \leq r }$ we have $${ \frac{1}{n} \sum _{j=1} ^{k} \tilde{w} _j ^{T} \Sigma _d ^2 \tilde{w} _j \leq \frac{1}{n} (\sigma _1 ^2 + \ldots + \sigma _k ^2) },$$ with equality attained for eg when $${ [\tilde{w} _1, \ldots, \tilde{w} _k] = [e _1, \ldots, e _k] }$$ that is ${ U ^{T} [w _1, \ldots, w _k] = [ e _1, \ldots, e _k] }.$ So the required orthonormal set is $${ [w _1, \ldots, w _k] = [U _1, \ldots, U _k ] \text{ for } k \leq r }.$$

Finally, for ${ k \leq r },$ the set ${ \overline{x} + \text{span}(U _1, \ldots, U _k) }$ is the best fit ${ k }$ dimensional affine subspace for the data.

The corresponding minimum loss is $${ \begin{align*} &\frac{1}{n} \sum _{i=1} ^{n} \text{dist}(x _i, \overline{x} + \text{span}(U _1, \ldots, U _k)) ^2 \\ &= \frac{1}{n} \sum _{i=1} ^{n} \lVert x _i - \overline{x} \rVert ^2 - \frac{1}{n} (\sigma _1 ^2 + \ldots + \sigma _k ^2 ) \\ &= \frac{1}{n} \text{tr}(X ^{T} X) - \frac{1}{n} ( \sigma _1 ^2 + \ldots + \sigma _k ^2) \\ &= \frac{1}{n} (\sigma _{k+1} ^2 + \ldots + \sigma _r ^2). \end{align*} }$$

One can also consider the best rank ${ k }$ approximation of ${ X = [x _1 - \overline{x}, \ldots, x _{n} - \overline{x}] \in \mathbb{R} ^{d \times n} }$ to get the same result.