5.9. SVD and PCA#
In this notebook, we will learn how SVD and PCA are related. We can use SVD to find the principal axes (also called loadings) of a dataset directly, without explicitly computing the covariance matrix: for a centered data matrix \(A = U\Sigma V^T\), the columns of \(V\) are exactly the eigenvectors of the covariance matrix of \(A\). The principal components themselves — the data projected onto these axes — are then given by \(AV = U\Sigma\).
To see this, we need to recall how to diagonalize a matrix with its eigenvalues and eigenvectors.
A symmetric matrix \(A \in \mathbb{R}^{n \times n}\) can be diagonalized using its eigenvalues and eigenvectors. If \(A\) has eigenvalues \(\lambda_1, \lambda_2, \ldots, \lambda_n\) and corresponding eigenvectors \(\mathbf{v}_1, \mathbf{v}_2, \ldots, \mathbf{v}_n\), then we can write
where \(V = [\mathbf{v}_1 \mid \mathbf{v}_2 \mid \cdots \mid \mathbf{v}_n]\) is the matrix whose columns are the eigenvectors, and \(\Lambda = \text{diag}(\lambda_1, \lambda_2, \ldots, \lambda_n)\) is the diagonal matrix of eigenvalues.
For a symmetric matrix, the eigenvectors are orthogonal, so \(V^T V = I\) and \(V\) is an orthogonal matrix. This decomposition is called the eigendecomposition or spectral decomposition.
Setup
For PCA, we have a data matrix \(X \in \mathbb{R}^{n \times p}\) where \(n\) is the number of observations (data points) and \(p\) is the number of features (dimensions). Each row represents one observation.
PCA requires centering the data by subtracting the mean of each feature
where \(\bar{X}\) is a matrix where each row contains the mean of all observations.
PCA finds the principal components by analyzing the covariance matrix. By definition, the covariance matrix is
Let us examine the SVD decomposition of the centered data
The right singular vectors \(V\) are exactly the principal components (eigenvectors of \(C\)). To see this, we can verify that \(V\) diagonalizes the covariance matrix
where we used \(U^T U = I\) in the last step. This is the eigendecomposition of \(C\). Therefore, the columns of \(V\) are the principal components, and the variances along these principal components are \(\frac{\sigma_i^2}{n-1}\).
Practical Implementation of SVD
Traditional PCA involves centering the data, computing the covariance matrix \(C = \frac{1}{n-1}X^T X\), computing the eigendecomposition of \(C\), and extracting the principal components as the eigenvectors.
PCA via SVD involves centering the data, computing the SVD \(X_{\text{centered}} = U\Sigma V^T\), and extracting the principal components as the columns of \(V\). The variance explained by each component is \(\frac{\sigma_i^2}{n-1}\).
There are several advantages to computing PCA via SVD:
Numerical stability: SVD is more numerically stable than computing \(X^T X\) followed by an eigendecomposition
Efficiency: For tall matrices (\(n \gg p\)), SVD can be more efficient
Avoids precision loss: Does not require forming \(X^T X\) explicitly, which can lose precision
Professional-grade implementations of PCA (such as sklearn.decomposition.PCA) use SVD internally for these reasons.
Example Exercise
Let us verify the SVD approach to PCA by comparing the SVD method to the eigenvalue method and sklearn results. We generate data from a bivariate normal distribution with a known mean and covariance matrix determined by \(\sigma_x\), \(\sigma_y\), and \(\rho\).
Bivariate Normal Distribution
A bivariate normal distribution has a probability density function \(\phi\) given by
where \(\mathbf{x} = (x_1, x_2)\) is a two-component vector of outcomes, \(\boldsymbol{\mu}\) is a vector denoting the mean along the \(x_1\) and \(x_2\) directions, respectively, and \(\Sigma\) is the covariance matrix. The probability density function \(\phi\) represents the probability density to draw an outcome \(\mathbf{x} = (x_1, x_2)\). More specifically, \(\phi(\mathbf{x})dx_1 dx_2\) is the probability to find \(\mathbf{x} = (x_1, x_2)\) in the infinitesimal rectangle \([x_1, x_1 + dx_1] \times [x_2, x_2 + dx_2]\).
The covariance matrix \(\Sigma\) is a symmetric, positive definite matrix that represents the variance and covariance between multidimensional random variables. For the bivariate case, the covariance matrix can be constructed from the correlation coefficient \(\rho\)
where \(\sigma_x\) and \(\sigma_y\) are the standard deviations in the \(x\) and \(y\) directions. Note that the covariance matrix is symmetric: \(\Sigma^T = \Sigma\).
The correlation coefficient \(\rho\) is a number between \(-1\) and \(1\). When \(\rho = 0\), the distributions along each axis are uncorrelated. When \(\rho = 1\), the distributions are maximally positively correlated, and when \(\rho = -1\), they are maximally negatively correlated.
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
# ============================================================================
# Generate Data from Bivariate Normal Distribution
# ============================================================================
# Set random seed for reproducibility
np.random.seed(42)
# Define parameters for bivariate normal distribution
n = 500 # Number of observations
mu_x, mu_y = 3.0, 4.0 # Means
sigma_x, sigma_y = 2.5, 4.0 # Standard deviations
rho = 0.6 # Correlation coefficient
# Construct mean vector and covariance matrix
mean = np.array([mu_x, mu_y])
cov_matrix = np.array([
[sigma_x**2, rho * sigma_x * sigma_y],
[rho * sigma_x * sigma_y, sigma_y**2]
])
print("True Covariance Matrix:")
print(cov_matrix)
print()
# Generate data
X = np.random.multivariate_normal(mean, cov_matrix, size=n)
print(f"Data shape: {X.shape}")
print(f"Data mean: {X.mean(axis=0)}")
print(f"Sample covariance:\n{np.cov(X.T)}")
True Covariance Matrix:
[[ 6.25 6. ]
[ 6. 16. ]]
Data shape: (500, 2)
Data mean: [2.93495196 4.01364237]
Sample covariance:
[[ 5.8228266 5.41899649]
[ 5.41899649 15.15016765]]
# ============================================================================
# Visualize Generated Data
# ============================================================================
fig, ax = plt.subplots(1, 1, figsize=(8, 6))
ax.scatter(X[:, 0], X[:, 1], alpha=0.6, s=30, edgecolors='k', linewidths=0.5)
ax.set_xlabel('$x_1$', fontsize=14)
ax.set_ylabel('$x_2$', fontsize=14)
ax.set_title(f'Bivariate Normal Data (n={n}, ρ={rho})', fontsize=16, fontweight='bold')
ax.grid(True, alpha=0.3)
ax.axhline(y=mu_y, color='red', linestyle='--', alpha=0.5, label=f'Mean: ({mu_x:.1f}, {mu_y:.1f})')
ax.axvline(x=mu_x, color='red', linestyle='--', alpha=0.5)
ax.legend(fontsize=12)
plt.tight_layout()
plt.show()
print(f"\nCorrelation coefficient: {rho}")
print(f"Standard deviations: σ_x = {sigma_x}, σ_y = {sigma_y}")
Correlation coefficient: 0.6
Standard deviations: σ_x = 2.5, σ_y = 4.0
# ============================================================================
# Method 1: Traditional PCA via Eigendecomposition
# ============================================================================
print("Method 1: Traditional PCA (Eigendecomposition)")
print("-"*60)
# Center the data
X_centered = X - X.mean(axis=0)
# Compute sample covariance matrix
C = (X_centered.T @ X_centered) / (n - 1)
print("\nSample Covariance Matrix:")
print(C)
# Compute eigendecomposition
eigenvalues, eigenvectors = np.linalg.eig(C)
# Sort by eigenvalues in descending order
idx = eigenvalues.argsort()[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]
print("\nEigenvalues (variances):")
print(eigenvalues)
print("\nEigenvectors (principal axes):")
print(eigenvectors)
print("\nVariance explained (%):")
print(eigenvalues / eigenvalues.sum() * 100)
Method 1: Traditional PCA (Eigendecomposition)
------------------------------------------------------------
Sample Covariance Matrix:
[[ 5.8228266 5.41899649]
[ 5.41899649 15.15016765]]
Eigenvalues (variances):
[17.63599681 3.33699744]
Eigenvectors (principal axes):
[[-0.41694888 -0.90892994]
[-0.90892994 0.41694888]]
Variance explained (%):
[84.08907475 15.91092525]
# ============================================================================
# Method 2: PCA via SVD
# ============================================================================
print("\nMethod 2: PCA via SVD")
print("-"*60)
# Center the data (already done above, but included for clarity)
X_centered = X - X.mean(axis=0)
# Compute SVD
U, sigma, Vt = np.linalg.svd(X_centered, full_matrices=False)
# Extract principal components (columns of V)
principal_axes_svd = Vt.T
# Compute variances from singular values
variances_svd = (sigma**2) / (n - 1)
print("\nSingular values:")
print(sigma)
print("\nVariances (σ² / (n-1)):")
print(variances_svd)
print("\nPrincipal axes (right singular vectors):")
print(principal_axes_svd)
print("\nVariance explained (%):")
print(variances_svd / variances_svd.sum() * 100)
Method 2: PCA via SVD
------------------------------------------------------------
Singular values:
[93.81024682 40.80639314]
Variances (σ² / (n-1)):
[17.63599681 3.33699744]
Principal axes (right singular vectors):
[[-0.41694888 -0.90892994]
[-0.90892994 0.41694888]]
Variance explained (%):
[84.08907475 15.91092525]
# ============================================================================
# Method 3: sklearn PCA
# ============================================================================
print("\nMethod 3: sklearn PCA")
print("-"*60)
# Create PCA object and fit to data
pca = PCA(n_components=2)
pca.fit(X)
print("\nVariances (explained variance):")
print(pca.explained_variance_)
print("\nPrincipal axes:")
print(pca.components_.T) # Transpose to match our convention
print("\nVariance explained (%):")
print(pca.explained_variance_ratio_ * 100)
Method 3: sklearn PCA
------------------------------------------------------------
Variances (explained variance):
[17.63599681 3.33699744]
Principal axes:
[[ 0.41694888 0.90892994]
[ 0.90892994 -0.41694888]]
Variance explained (%):
[84.08907475 15.91092525]
As another example we will compute PCA on the breast cancer dataset. Instead of using sklearn’s PCA class, we will compute it directly from the singular value decomposition, and plot the first two principal components.
Recall the two related but distinct quantities from the SVD, \(A = U\Sigma V^T\), where \(A\) is the centered data matrix. The columns of \(V\) are the called the loadings or principle axes, they inform us on the directions in the original feature space along which the data varies most. The principal components are the data itself, projected coordinates along those directions, computed as \(AV = U\Sigma\).
from sklearn.datasets import load_breast_cancer
data = load_breast_cancer()
X = data.data
y = data.target
target_names = data.target_names
X_centered = X - X.mean(axis=0)
U, S, Vt = np.linalg.svd(X_centered, full_matrices=False)
loadings = Vt.T # principal axes: columns of V
principal_components = U * S # equivalent to X_centered @ loadings
explained_variance = S**2 / (X.shape[0] - 1)
explained_variance_ratio = explained_variance / explained_variance.sum()
fig, ax = plt.subplots(figsize=(7, 5.5))
for label, name, color in zip([0, 1], target_names, ['tab:red', 'tab:blue']):
mask = y == label
ax.scatter(principal_components[mask, 0], principal_components[mask, 1],
s=15, alpha=0.7, label=name, color=color)
ax.set_xlabel(f"PC1 ({explained_variance_ratio[0]*100:.1f}% variance)")
ax.set_ylabel(f"PC2 ({explained_variance_ratio[1]*100:.1f}% variance)")
ax.legend()
plt.tight_layout()
5.9.1. More on the PCA Theory#
Let \(X_{\text{centered}}\) be the \(N \times n\) data matrix whose rows are the centered data vectors, so that \(\mathbf{x}_i \in \mathbb{R}^n\) denotes the \(i\)-th row of \(X_{\text{centered}}\), written as a column vector. The sample covariance matrix is
where the second equality holds because \((X_{\text{centered}}^T X_{\text{centered}})_{jk} = \sum_i (\mathbf{x}_i)_j (\mathbf{x}_i)_k = \left(\sum_i \mathbf{x}_i \mathbf{x}_i^T\right)_{jk}\).
The principal axes \(\mathbf{v}_1, \mathbf{v}_2, \ldots, \mathbf{v}_{n}\) create an orthonormal basis for \(\mathbb{R}^{n}\).
The component of a data vector \(\mathbf{x}_i\) in the direction of principal axis \(\mathbf{v}_k\) is \((\mathbf{v}_k^T\mathbf{x}_i)\mathbf{v}_k\), so \(\mathbf{x}_i\) can be represented exactly as a linear combination of the principal axes as $\(\mathbf{x}_i = \sum_{k = 1}^{n}(\mathbf{v}_k^T\mathbf{x}_i) \mathbf{v}_k.\)$
We can find a one-dimensional approximation to \(\mathbf{x}_i\) using the component along the first principal axis \(\mathbf{x}_i \approx (\mathbf{v}_1^T\mathbf{x}_i) \mathbf{v}_1\). More generally, we can find a \(D\)-dimensional approximation as \(\mathbf{x}_i \approx \sum_{k = 1}^{D}(\mathbf{v}_k^T\mathbf{x}_i) \mathbf{v}_k\).
In what sense is this a good approximation?
We can show that finding the direction of maximum variance $\(\max_{\mathbf{u}} \mathbf{u}^TC\mathbf{u}\)\( \)\(\text{subject to } \mathbf{u}^T\mathbf{u} = 1\)$ is equivalent to finding the one-dimensional approximation that minimizes the squared error
To see the equivalence, we will expand the error to isolate the term that depends on \(\mathbf{u}\):
Since \(\mathbf{x}_i^T\mathbf{x}_i\) does not depend on \(\mathbf{u}\), minimizing \(\frac{1}{N-1}\sum_{i = 1}^N \left(\mathbf{x}_i^T\mathbf{x}_i - (\mathbf{u}^T\mathbf{x}_i)^2\right)\) over \(\mathbf{u}\) is equivalent to maximizing \(\frac{1}{N-1}\sum_{i=1}^N (\mathbf{u}^T\mathbf{x}_i)^2\), so
So, the best one-dimensional approximation to \(\mathbf{x}_i\) is along the first principal axis \(\mathbf{x}_i \approx (\mathbf{v}_1^T\mathbf{x}_i) \mathbf{v}_1\).
Extending this argument shows that the best \(D\)-dimensional approximation is \(\mathbf{x}_i \approx \sum_{k = 1}^{D}(\mathbf{v}_k^T\mathbf{x}_i) \mathbf{v}_k\).