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

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?

• I wrote this FAQ-style question together with my own answer, because it is frequently being asked in various forms, but there is no canonical thread and so closing duplicates is difficult. Please provide meta comments in this accompanying meta thread. Commented Jan 22, 2015 at 11:25
• stats.stackexchange.com/questions/177102/… Commented Feb 3, 2016 at 10:45
• In addition to an excellent and detailed amoeba's answer with its further links I might recommend to check this, where PCA is considered side by side some other SVD-based techniques. The discussion there presents algebra almost identical to amoeba's with just minor difference that the speech there, in describing PCA, goes about svd decomposition of $\mathbf X/\sqrt{n}$ [or $\mathbf X/\sqrt{n-1}$] instead of $\bf X$ - which is simply convenient as it relates to the PCA done via the eigendecomposition of the covariance matrix. Commented Feb 3, 2016 at 12:18
• PCA is a special case of SVD. PCA needs the data normalized, ideally same unit. The matrix is nxn in PCA. Commented Oct 17, 2017 at 9:12
• @OrvarKorvar: What n x n matrix are you talking about ? Commented Mar 29, 2018 at 15:16

## 5 Answers

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, 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

• +1 for both Q&A. Thanks for sharing. I have one question: why do you have to assume that the data matrix is centered initially? Commented Aug 6, 2015 at 8:39
• @Antoine, covariance matrix is by definition equal to $\langle (\mathbf x_i - \bar{\mathbf x})(\mathbf x_i - \bar{\mathbf x})^\top \rangle$, where angle brackets denote average value. If all $\mathbf x_i$ are stacked as rows in one matrix $\mathbf X$, then this expression is equal to $(\mathbf X - \bar{\mathbf X})(\mathbf X - \bar{\mathbf X})^\top/(n-1)$. If $\mathbf X$ is centered then it simplifies to $\mathbf X \mathbf X^\top/(n-1)$. Think of variance; it's equal to $\langle (x_i-\bar x)^2 \rangle$. But if $\bar x=0$ (i.e. data are centered), then it's simply the average value of $x_i^2$. Commented Aug 6, 2015 at 9:43
• @amoeba for those less familiar with linear algebra and matrix operations, it might be nice to mention that $(A.B.C)^{T}=C^{T}.B^{T}.A^{T}$ and that $U^{T}.U=Id$ because $U$ is orthogonal Commented Feb 3, 2016 at 15:20
• A code sample for PCA by SVD: stackoverflow.com/questions/3181593/… Commented Mar 23, 2016 at 11:51
• @amoeba yes, but why use it? Also, is it possible to use the same denominator for $S$? The problem is that I see formulas where $\lambda_i = s_i^2$ and try to understand, how to use them?
– Dims
Commented Sep 5, 2017 at 22:49

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)


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.

• Thanks for your anser Andre. Just two small typos correction: 1. In the last paragraph youre confusing left and right. 2. In the (capital) formula for X, you're using v_j instead of v_i.
– Alon
Commented Sep 3, 2019 at 13:09

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)

[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.