Principal Component Analysis & SVD | AI & Machine Learning Engineering

Real-world datasets often possess high ambient dimensionality while lying predominantly on low-dimensional linear or affine subspaces. High dimensionality inflates sample complexity, exacerbates the curse of dimensionality, and renders direct visualization impossible.

Principal Component Analysis (PCA) is the foundational linear dimensionality reduction technique. It rotates the coordinate frame along directions of maximal variance, decorrelating features and yielding optimal low-rank reconstructions.

Spectral IdentityLinear Algebra

Sigma * u = lambda * u

The optimal orthogonal projection axes correspond exactly to the eigenvectors of the empirical covariance matrix.

1. Variance Maximization Derivation#

Let $X \in \mathbb{R}^{N \times d}$ denote a zero-centered dataset where sample mean $\bar{x} = \frac{1}{N}\sum_{i=1}^N x_i = 0$ .

The empirical covariance matrix $\Sigma \in \mathbb{R}^{d \times d}$ is symmetric and positive semi-definite:

\Sigma = \frac{1}{N} X^T X = \frac{1}{N} \sum_{i=1}^N x_i x_i^T

Architecture Diagram

Rendering diagram...

We seek a unit projection vector $u_1 \in \mathbb{R}^d$ ( $\|u_1\|_2^2 = u_1^T u_1 = 1$ ) such that projecting data points onto $u_1$ maximizes the sample variance.

Step 1: Formulate the Projected Variance

The scalar projection of sample $x_i$ onto unit vector $u_1$ is $y_{i,1} = x_i^T u_1$ .

Because $X$ is zero-centered, the mean of the projected coordinates is also zero:

\bar{y}_1 = \frac{1}{N} \sum_{i=1}^N x_i^T u_1 = \left( \frac{1}{N} \sum_{i=1}^N x_i^T \right) u_1 = 0

The empirical variance $\sigma_1^2$ along direction $u_1$ is:

\sigma_1^2 = \frac{1}{N} \sum_{i=1}^N (x_i^T u_1)^2 = \frac{1}{N} \sum_{i=1}^N (u_1^T x_i)(x_i^T u_1) = u_1^T \left( \frac{1}{N} \sum_{i=1}^N x_i x_i^T \right) u_1 = u_1^T \Sigma u_1

Step 2: Formulate Constrained Optimization via Lagrange Multipliers

We maximize $u_1^T \Sigma u_1$ subject to the unit-norm constraint $u_1^T u_1 = 1$ :

\mathcal{L}(u_1, \lambda_1) = u_1^T \Sigma u_1 - \lambda_1 (u_1^T u_1 - 1)

Step 3: Differentiate and Solve for Stationarity

Differentiate $\mathcal{L}$ with respect to vector $u_1$ :

\nabla_{u_1} \mathcal{L} = 2 \Sigma u_1 - 2 \lambda_1 u_1 = 0

Dividing by 2 yields the fundamental Eigenvalue Equation:

\Sigma u_1 = \lambda_1 u_1

Step 4: Maximize the Objective

Premultiply both sides by $u_1^T$ :

u_1^T \Sigma u_1 = \lambda_1 u_1^T u_1 = \lambda_1

Because $u_1^T \Sigma u_1$ is the projected variance, the variance along direction $u_1$ equals the eigenvalue $\lambda_1$ . To maximize variance, $u_1$ must be the eigenvector of covariance matrix $\Sigma$ corresponding to the largest eigenvalue $\lambda_{\max}$ .

Subsequent principal components $u_2, \dots, u_k$ are obtained by repeating this optimization under the additional mutual orthogonality constraints $u_j^T u_m = 0$ for all $m < j$ .

2. Equivalence to Reconstruction Error Minimization#

An alternative derivation of PCA seeks a $k$ -dimensional orthogonal projection subspace that minimizes the mean squared reconstruction error between original samples $x_i$ and reconstructed vectors $\tilde{x}_i = \sum_{j=1}^k (x_i^T u_j) u_j$ .

By the Pythagorean theorem in Euclidean metric space:

\|x_i\|_2^2 = \|\tilde{x}_i\|_2^2 + \|x_i - \tilde{x}_i\|_2^2

Summing over all $N$ data points and dividing by $N$ :

\frac{1}{N} \sum_{i=1}^N \|x_i\|_2^2 = \frac{1}{N} \sum_{i=1}^N \|\tilde{x}_i\|_2^2 + \frac{1}{N} \sum_{i=1}^N \|x_i - \tilde{x}_i\|_2^2

Observe that the left side $\frac{1}{N} \sum_{i=1}^N \|x_i\|_2^2 = \text{trace}(\Sigma)$ is a constant independent of projection basis $U$ .

Therefore, minimizing reconstruction error $\frac{1}{N} \sum_{i=1}^N \|x_i - \tilde{x}_i\|_2^2$ is mathematically identical to maximizing projected variance $\frac{1}{N} \sum_{i=1}^N \|\tilde{x}_i\|_2^2$ .

3. Singular Value Decomposition (SVD) Connection#

Directly computing covariance matrix $\Sigma = \frac{1}{N} X^T X$ incurs $O(N \cdot d^2)$ operations and risks squaring the condition number of ill-conditioned data matrices.

Modern PCA implementations compute principal components directly via Singular Value Decomposition (SVD) of centered data matrix $X \in \mathbb{R}^{N \times d}$ :

X = U S V^T

Where:

$U \in \mathbb{R}^{N \times N}$ is orthogonal ( $U^T U = I$ ) containing left singular vectors.
$S \in \mathbb{R}^{N \times d}$ is rectangular diagonal containing singular values $s_1 \ge s_2 \ge \dots \ge 0$ .
$V \in \mathbb{R}^{d \times d}$ is orthogonal ( $V^T V = I$ ) containing right singular vectors.

Evaluating $X^T X$ :

X^T X = (V S^T U^T)(U S V^T) = V (S^T S) V^T

Therefore, the right singular vectors $V$ of $X$ are identically the eigenvectors of covariance matrix $X^T X$ , and the eigenvalues relate to singular values by:

\lambda_k = \frac{s_k^2}{N}

The low-dimensional projection $Y = X V_k$ evaluates without explicit covariance matrix construction.

4. From-Scratch Implementation: Zero Dependencies#

The following implementations compute PCA from scratch in Python and Rust using Power Iteration to extract dominant eigenvectors without relying on third-party linear algebra packages.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
import math
from typing import List, Tuple

class PCAFromScratch:
    def __init__(self, n_components: int = 2, max_iter: int = 300, tol: float = 1e-6):
        self.n_components = n_components
        self.max_iter = max_iter
        self.tol = tol
        self.components: List[List[float]] = []
        self.mean: List[float] = []

    @staticmethod
    def _dot(a: List[float], b: List[float]) -> float:
        return sum(x * y for x, y in zip(a, b))

    @staticmethod
    def _norm(a: List[float]) -> float:
        return math.sqrt(sum(x ** 2 for x in a))

    def _power_iteration(self, cov: List[List[float]], num_features: int) -> Tuple[float, List[float]]:
        # Initialize candidate vector
        b_k = [1.0 / math.sqrt(num_features)] * num_features

        for _ in range(self.max_iter):
            # Matrix-vector product: cov * b_k
            b_next = [self._dot(cov[i], b_k) for i in range(num_features)]
            norm = self._norm(b_next)
            if norm == 0.0:
                break
            b_next = [x / norm for x in b_next]

            # Check convergence
            diff = self._norm([x - y for x, y in zip(b_next, b_k)])
            b_k = b_next
            if diff < self.tol:
                break

        # Rayleigh quotient: lambda = b^T * cov * b
        eigenvalue = self._dot(b_k, [self._dot(cov[i], b_k) for i in range(num_features)])
        return eigenvalue, b_k

    def fit(self, X: List[List[float]]) -> None:
        num_samples = len(X)
        num_features = len(X[0])

        # 1. Compute empirical mean and center data
        self.mean = [sum(row[j] for row in X) / num_samples for j in range(num_features)]
        X_centered = [[row[j] - self.mean[j] for j in range(num_features)] for row in X]

        # 2. Compute covariance matrix: 1/N * X^T * X
        cov = [[0.0] * num_features for _ in range(num_features)]
        for i in range(num_features):
            for j in range(i, num_features):
                val = sum(X_centered[k][i] * X_centered[k][j] for k in range(num_samples)) / num_samples
                cov[i][j] = val
                cov[j][i] = val

        # 3. Extract top components using power iteration and deflation
        self.components = []
        for _ in range(self.n_components):
            eigenval, eigenvector = self._power_iteration(cov, num_features)
            self.components.append(eigenvector)

            # Deflate covariance matrix: cov = cov - lambda * (v * v^T)
            for i in range(num_features):
                for j in range(num_features):
                    cov[i][j] -= eigenval * eigenvector[i] * eigenvector[j]

    def transform(self, X: List[List[float]]) -> List[List[float]]:
        num_features = len(self.mean)
        projected = []
        for row in X:
            centered = [row[j] - self.mean[j] for j in range(num_features)]
            coords = [self._dot(centered, comp) for comp in self.components]
            projected.append(coords)
        return projected

if __name__ == "__main__":
    X = [
        [2.5, 2.4], [0.5, 0.7], [2.2, 2.9],
        [1.9, 2.2], [3.1, 3.0], [2.3, 2.7],
        [2.0, 1.6], [1.0, 1.1], [1.5, 1.6], [1.1, 0.9]
    ]

    pca = PCAFromScratch(n_components=1)
    pca.fit(X)
    X_proj = pca.transform(X)
    print("Principal Component Vector:", pca.components[0])
    print("Projected 1D Coordinates (first 3):", X_proj[:3])

5. Active Knowledge Checks#

Concept check

In the constrained optimization derivation of PCA via Lagrange multipliers, why does the stationarity condition nabla_u L = 0 resolve directly to the eigenvalue equation Sigma * u = lambda * u?

Concept check

How does the Singular Value Decomposition X = U * S * V^T relate to the Principal Components of covariance matrix 1/N * X^T * X?