Connection Between SVD and Principal Component Analysis

import matplotlib.pyplot as plt
%matplotlib inline
import matplotlib_inline
matplotlib_inline.backend_inline.set_matplotlib_formats('svg')
import seaborn as sns

Connection Between SVD and Principal Component Analysis#

PCA is a linear dimensionality-reduction technique. We typically use PCA in supervised or unsupervised learning to reduce the number of features in the dataset. The main idea behind PCA is to find a new set of features that are uncorrelated and ordered by the amount of variance they explain.

Let \(\mathbf{X}\) be an \(n \times m\) matrix, where \(n\) is the number of samples and \(m\) is the number of features. We will use \(\mathbf{x}_i\) to denote the \(i\)-th row of \(\mathbf{X}\). So, write:

\[\begin{split} \mathbf{X} = \begin{bmatrix} -\mathbf{x}_1- \\ -\mathbf{x}_2-\\ \vdots \\ -\mathbf{x}_n- \end{bmatrix}, \end{split}\]

where:

\[ -\mathbf{x}_j- \equiv \mathbf{x}_j^T. \]

When we do PCA, we always want to center the data. To this end, we calculate the empirical mean of the data:

\[ \bar{\mathbf{x}} = \langle \mathbf{x}_i \rangle = \frac{1}{n} \sum_{i=1}^n \mathbf{x}_i, \]

and we make the centered data matrix:

\[\begin{split} \mathbf{B} = \mathbf{X} - \bar{\mathbf{x}} = \begin{bmatrix} -\mathbf{x}_1-\bar{\mathbf{x}}- \\ -\mathbf{x}_2-\bar{\mathbf{x}}- \\ \vdots \\ -\mathbf{x}_n-\bar{\mathbf{x}}- \end{bmatrix}. \end{split}\]

Now take an arbitrary unit vector \(\mathbf{v}\).

Hide code cell source

 
import numpy as np
A = np.linalg.cholesky([[1, 0.9], [0.9, 1]])
X = (np.array([3, 3])[:, None] + A @ np.random.randn(2, 500)).T
bar_x = X.mean(axis=0)
v = np.array([np.cos(0.5), np.sin(0.5)])
x_i = np.array([4, 5])
b_i = x_i - bar_x
fig, ax = plt.subplots(figsize=FIGURE_SIZES["half_standard"])
ax.plot(X[:, 0], X[:, 1], '.', alpha=0.2)
ax.plot(bar_x[0], bar_x[1], 'o', color='red')
ax.quiver(0, 0, x_i[0], x_i[1], angles='xy', scale_units='xy', scale=1)
ax.text(3, 4, '$x_i$', fontsize=12)
ax.quiver(0, 0, bar_x[0], bar_x[1], angles='xy', scale_units='xy', scale=1, color='red')
ax.text(bar_x[0], bar_x[1], '$\\bar{x}$', fontsize=12, color='red')
ax.quiver(bar_x[0], bar_x[1], v[0], v[1], angles='xy', scale_units='xy', scale=1, color='green')
ax.text(bar_x[0] + v[0], bar_x[1] + v[1], '$v$', fontsize=12, color='green')
ax.quiver(bar_x[0], bar_x[1], b_i[0], b_i[1], angles='xy', scale_units='xy', scale=1, color='blue')
ax.text(bar_x[0] + b_i[0], bar_x[1] + b_i[1], '$b_i$', fontsize=12, color='blue')
finalize_axes(keep_box=False)

array([<Axes: >], dtype=object)
A two-dimensional data cloud with the sample mean, one observation vector, its centered vector, and a candidate projection direction drawn as labeled arrows.

Consider the projection of the centered sample \(\mathbf{b}_i=\mathbf{x}_i-\bar{\mathbf{x}}\) (the \(i\)-th row of \(\mathbf{B}\)) on \(\mathbf{v}\). It is:

\[ \text{Proj}_i = \mathbf{b}_i^T \mathbf{v}. \]

What is the empirical mean of the projections?

\[ \langle \text{Proj}_i \rangle = \langle \mathbf{b}_i^T \mathbf{v} \rangle = \langle \mathbf{b}_i\rangle^T \mathbf{v} = 0, \]

because \(\mathbf{b}_i\) is centered.

What is the variance of the projections?

\[ \langle \text{Proj}_i^2\rangle = \langle \mathbf{v}^T\mathbf{b}_i \mathbf{b}_i^T \mathbf{v} \rangle = \mathbf{v}^T \langle \mathbf{b}_i \mathbf{b}_i^T \rangle \mathbf{v} = \mathbf{v}^T \mathbf{C} \mathbf{v}, \]

where \(\mathbf{C}\) is the covariance matrix of the centered data:

\[ \mathbf{C} = \frac{1}{n} \mathbf{B}^T \mathbf{B}. \]

PCA finds \(\mathbf{v}\) by solving the following problem:

\[ \mathbf{v}_1 = \arg\max \mathbf{v}^T \mathbf{C} \mathbf{v}, \]

subject to the constraint that \(\mathbf{v}\) is a unit vector.

Using the method of Lagrange multipliers, we can show that the solution to this problem is the eigenvector of \(\mathbf{C}\) with the largest eigenvalue.

In a similar way, we can find the second principal direction, \(\mathbf{v}_2\). And so on. We always get a sequence of orthogonal unit vectors, \(\mathbf{v}_1, \mathbf{v}_2, \ldots, \mathbf{v}_m\), corresponding to the eigenvalues of \(\mathbf{C}\) in decreasing order.

And now, here is the connection between PCA and SVD. Do the SVD of the centered matrix \(\mathbf{B}\):

\[ \mathbf{B} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^T. \]

Here \(\mathbf{\Sigma}\) is \(n\times m\), so \(\mathbf{\Sigma}^T\mathbf{\Sigma}\) is an \(m\times m\) diagonal matrix. If \(m>n\), its last \(m-n\) entries are zero. Form the covariance:

\[\begin{split} \begin{aligned} \mathbf{C} &= \frac{1}{n}\mathbf{B}^T\mathbf{B} \\ &= \frac{1}{n}\mathbf{V}\mathbf{\Sigma}^T(\mathbf{U}^T\mathbf{U})\mathbf{\Sigma}\mathbf{V}^T \\ &= \frac{1}{n}\mathbf{V}\mathbf{\Sigma}^T\mathbf{\Sigma}\mathbf{V}^T. \end{aligned} \end{split}\]

So, we get that the SVD diagonalizes the covariance matrix. We can read off the eigenvalues and eigenvectors. The \(j\)-th column of \(\mathbf{V}\) is the \(j\)-th principal direction. The \(j\)-th singular value squared and divided by \(n\) is the variance in the \(j\)-th principal direction.

The total variance of the data is:

\[ \text{Tr}(\mathbf{C}) = \frac{1}{n} \text{Tr}(\mathbf{\Sigma}^T\mathbf{\Sigma}) = \sum_{j=1}^{\min(n,m)} \frac{\sigma_j^2}{n}. \]

The projection coefficients are known as the principal components. They are:

\[ \mathbf{Z} = \mathbf{B} \mathbf{V} = \mathbf{U} \mathbf{\Sigma}. \]

The principal components are uncorrelated#

Let \(\mathbf{z}_i\) be the \(i\)-th row of \(\mathbf{Z}\). Then:

\[ \langle z_{ik}z_{il}\rangle = \langle u_{ik}\sigma_k u_{il}\sigma_l\rangle = \sigma_k \sigma_l \langle u_{ik} u_{il}\rangle = \frac{\sigma_k^2}{n} \delta_{kl}. \]

The displayed calculation uses \(k,l\leq\min(n,m)\); any remaining principal components are identically zero. Here \(\langle\cdot\rangle\) denotes the average over the sample index \(i\) (there is no sum over \(k\) or \(l\)), and we used the fact that the columns of \(\mathbf{U}\) are orthonormal.

Example: PCA on the MNIST dataset#

We use the MNIST dataset (LeCun et al., n.d.) to illustrate PCA. The notebook loads OpenML dataset 554 (mnist_784, version 1).

from sklearn.datasets import fetch_openml
mnist = fetch_openml('mnist_784', version=1)

Let’s do the steps described above:

import numpy as np
# 1. Extract the relevant data
X = mnist.data
# 2. Find the empirical mean
x_bar = np.mean(X, axis=0)
# 3. Center the data
B = X - x_bar
# 4. Do the SVD
U, S, Vt = np.linalg.svd(B, full_matrices=False)
# 5. Project the data
Z = U @ np.diag(S)

Let’s look at the explained variance:

Hide code cell source

 
cum_var = np.cumsum(S**2)/np.sum(S**2)
best_k = np.argmax(cum_var > 0.9)

fig, ax = plt.subplots(figsize=FIGURE_SIZES["half_standard"])
ax.plot(cum_var)
ax.axhline(0.9, color='red', linestyle='--')
ax.axvline(best_k, color='red', linestyle='--')
ax.text(best_k, 0.5, f'90% at k={best_k}', verticalalignment='center', color='red')
ax.set(yscale='log', xlabel="Singular value index", ylabel="Cumulative fraction of total variance")
finalize_axes(keep_box=False)

array([<Axes: xlabel='Singular value index', ylabel='Cumulative fraction of total variance'>],
      dtype=object)
Cumulative explained variance of MNIST principal components versus component index, with the first rank exceeding 90 percent marked.

Let’s look at the first 10 principal directions:

Hide code cell source

 
fig, axes = plt.subplots(5, 5, figsize=FIGURE_SIZES["half_tall"])
for i, ax in enumerate(axes.ravel()):
    ax.imshow(Vt[i].reshape(28, 28), cmap='gray')
    ax.axis('off')
A five-by-five grid of leading MNIST principal directions rendered as grayscale 28-by-28 pixel patterns resembling combinations of digit strokes.

We fit PCA to all 70,000 images and display a fixed random subset of 100 images per digit. Marker shapes identify the digits; subsampling keeps the overlapping groups visible.

y = np.asarray(mnist.target, dtype=int)
display_rng = np.random.default_rng(697)
markers = ["o", "s", "^", "v", "D", "<", ">", "p", "h", "8"]
fig, ax = plt.subplots(
    figsize=FIGURE_SIZES["full_standard"], constrained_layout=True
)
for digit, marker in enumerate(markers):
    mask = display_rng.choice(np.flatnonzero(y == digit), size=100, replace=False)
    ax.scatter(
        Z[mask, 0],
        Z[mask, 1],
        label=str(digit),
        marker=marker,
        s=16,
        facecolors="none",
        edgecolors="black",
        linewidths=0.55,
        alpha=0.65,
    )
ax.legend(
    loc="upper center", bbox_to_anchor=(0.5, -0.14),
    ncol=5, title="Digit", handletextpad=0.35, columnspacing=0.9,
)
ax.set(xlabel="First principal component", ylabel="Second principal component")
finalize_axes(keep_box=False)

array([<Axes: xlabel='First principal component', ylabel='Second principal component'>],
      dtype=object)
MNIST digits projected onto the first two principal components, with digit classes forming overlapping but structured clusters.