---
title: PCA From Scratch in Python: The Math of Eigenvectors and Dimensionality Reduction Explained
source: https://app.sythra.ai/learn/machine-learning/pca-from-scratch-eigenvectors-dimensionality-reduction-python
topic: Machine Learning
updated: 2026-09-09
publisher: Sythra (https://app.sythra.ai)
---

# PCA From Scratch in Python: The Math of Eigenvectors and Dimensionality Reduction Explained

Principal Component Analysis (PCA) is an unsupervised linear dimensionality reduction technique that transforms correlated features into a set of linearly uncorrelated orthogonal axes called principal components. These components align with the directions of maximum variance in the data, derived mathematically as the eigenvectors of the feature covariance matrix. The corresponding eigenvalues quantify the exact variance preserved along each axis, allowing high-dimensional data to be compressed into fewer dimensions with minimal reconstruction loss.

_Source: [https://app.sythra.ai/learn/machine-learning/pca-from-scratch-eigenvectors-dimensionality-reduction-python](https://app.sythra.ai/learn/machine-learning/pca-from-scratch-eigenvectors-dimensionality-reduction-python) — free to read on Sythra._

## Key points

- Explains PCA geometrically using the teapot camera angle and 2D cigar scatter plot analogies.
- Derives the mathematical objective via Lagrange multipliers: proving why maximizing variance leads directly to the eigenvector equation $C v = \lambda v$.
- Provides a complete, step-by-step numerical hand trace on a 3-point dataset: calculating the covariance matrix $C$, solving $\det(C - \lambda I) = 0$, finding $\lambda_1 \approx 4.3028$ (86.06% variance), and projecting points.
- Explains why production libraries use Singular Value Decomposition (SVD) on $X$ directly rather than computing $X^T X$, avoiding condition number squaring $\kappa(X^T X) = \kappa(X)^2$.
- Implements a complete object-oriented PCA class from scratch in NumPy with fit, transform, and inverse_transform, verified side-by-side against Scikit-Learn.
- Covers stopping criteria (Scree plot elbow, 95% cumulative variance, Kaiser criterion) and contrasts PCA against LDA, Kernel PCA, t-SNE, and UMAP.

**Principal Component Analysis (PCA)** is the foundational workhorse of unsupervised machine learning and multivariate statistics. Whenever tabular tables grow to hundreds of columns, computer vision models ingest thousands of raw pixel intensities, or genomics pipelines process tens of thousands of gene expression levels, data scientists face the notorious **curse of dimensionality**. Distances between points collapse, predictive models overfit, training times balloon, and human visualization becomes impossible beyond three dimensions.

> **THE TEAPOT AND THE CAMERA ANGLE:** Imagine placing a complex, three-dimensional ceramic teapot on a table. Your job is to take a single two-dimensional photograph to display on an e-commerce website. A photograph is 2D, but the teapot is 3D — so your choice of camera angle dictates how much information survives. If you point your camera straight down from the ceiling directly above the lid, the teapot collapses into an uninformative circle; the spout, the handle, and the body contours are completely hidden. But if you walk around the table and position the camera at an angle showing the long profile, a single 2D snapshot captures almost the entire geometry. **PCA is an algorithm that computes the optimal camera angle for high-dimensional data clouds** — the vantage point that maximizes the visible spread (variance) while discarding redundant noise.

## 1. Key Concepts & Mathematical Notation Glossary

Before establishing the mathematical derivations, review the standard notation used throughout multivariate statistics, linear algebra, and machine learning literature:

| Symbol | Mathematical Concept | Dimensionality | Interpretation / Role in PCA |
| --- | --- | --- | --- |
| $X$ | Raw Data Matrix | $n \times d$ | $n$ observations (samples) across $d$ original input features. |
| $\mu$ | Empirical Feature Mean Vector | $1 \times d$ | Mean vector $\mu = \frac{1}{n} \sum_{i=1}^n x_i$ used to center the feature cloud at the origin. |
| $X_{\text{centered}}$ | Centered Data Matrix | $n \times d$ | Zero-mean data $X - \mu$; ensures variance is measured relative to the centroid. |
| $C$ | Sample Covariance Matrix | $d \times d$ | Symmetric matrix $C = \frac{1}{n-1} X_{\text{centered}}^T X_{\text{centered}}$ capturing pairwise feature covariances. |
| $v_i$ | Eigenvector / Principal Axis | $d \times 1$ | Orthogonal unit direction vector ($\\|v_i\\| = 1$) pointing along an axis of data variance. |
| $\lambda_i$ | Eigenvalue | Scalar | Variance magnitude of the data projected onto eigenvector $v_i$. |
| $W$ | Projection Matrix | $d \times k$ | Matrix formed by concatenating the top $k$ principal eigenvectors as columns. |
| $Z$ or $X_{\text{proj}}$ | Transformed / Latent Matrix | $n \times k$ | Compressed low-dimensional coordinates: $Z = X_{\text{centered}} W$. |
| $\text{EVR}_i$ | Explained Variance Ratio | Scalar $\in [0, 1]$ | Fraction of total variance accounted for by component $i$: $\lambda_i / \sum_{j=1}^d \lambda_j$. |
| $\Sigma, U, V^T$ | Singular Value Decomposition (SVD) | Matrices | Decomposition $X_{\text{centered}} = U \Sigma V^T$ used in production PCA for numerical stability. |

## 2. Why PCA Is Needed: Overcoming the Curse of Dimensionality

In modern data science, more features do not always yield better models. As the feature count $d$ grows relative to the sample count $n$, three critical problems emerge:

- **Multicollinearity and Redundant Information:** In financial and customer analytics, variables such as monthly income, annual tax bracket, and home value are highly correlated. Retaining all of them destabilizes regression coefficients, inflates standard errors, and adds computational weight without introducing novel signal.
- **Distance Concentration in High Dimensions:** In high-dimensional Euclidean space, the ratio between the distance to the nearest neighbor and the distance to the farthest neighbor approaches 1. Distance-based algorithms (such as KNN, K-Means, and SVM with RBF kernels) degrade because every point appears virtually equidistant from every other point.
- **Storage, Throughput, and Human Interpretability:** Processing raw 1024-dimensional embeddings or 4096-pixel images incurs severe GPU memory and latency penalties. Compressing features down to an informative low-dimensional subspace accelerates downstream gradient descent and enables 2D/3D scatter visualization.

## 3. Geometric Intuition: The Long Axis of the Data Cloud

To visualize the geometry of PCA, picture a two-dimensional scatter plot displaying two correlated student metrics: _hours studied_ on the horizontal axis and _final exam score_ on the vertical axis. Because students who study more generally score higher, the cloud of points forms an elongated, tilted ellipse — reminiscent of a cigar angled at roughly 45 degrees.

Now ask the central question of dimensionality reduction: _If you are allowed to draw only ONE straight line through this 2D cloud, which line preserves the most information?_

- If you project every dot horizontally onto the exam score axis, you discard how many hours they studied.
- If you project every dot vertically onto the study hours axis, you discard their exam score.
- However, if you draw a line straight through the **longest axis of the ellipse** and drop perpendicular shadows (orthogonal projections) onto it, the points that were far apart remain far apart. The variation along that diagonal captures the underlying academic commitment story.

That longest axis is the **first principal component (PC1)**. It is not an original feature; it is a linear combination of both features ($v_1 = w_1 \cdot \text{hours} + w_2 \cdot \text{score}$). The **second principal component (PC2)** is then constrained to be strictly orthogonal (perpendicular) to PC1, capturing whatever residual dispersion remains. In $d$ dimensions, PCA discovers $d$ orthogonal axes ordered strictly by the variance of the data along them.

## 4. Complete Mathematical Derivation of PCA

PCA can be derived from two equivalent mathematical perspectives: **maximizing the projected variance** or **minimizing the orthogonal reconstruction error**. Here, we present the classic variance maximization derivation using the method of Lagrange multipliers.

### Step 1: Centering and Standardizing the Data

Let $X \in \mathbb{R}^{n \times d}$ be the data matrix with $n$ rows and $d$ columns. Before computing directional spread, we must center each column around zero so the coordinate origin coincides with the data centroid:

$$\mu_j = \frac{1}{n} \sum_{i=1}^n X_{ij}, \qquad X_{\text{centered}} = X - \mathbf{1}\mu^T$$

When features have disparate units (for example, age in years ranging from 18 to 80, alongside annual income in dollars ranging from 20,000 to 500,000), features with massive numerical ranges would artificially dominate the variance. In such cases, we standardize each feature by dividing by its sample standard deviation: $Z_{ij} = (X_{ij} - \mu_j) / s_j$.

### Step 2: The Sample Covariance Matrix

The empirical sample covariance matrix $C \in \mathbb{R}^{d \times d}$ summarizes the joint dispersion between all pairs of features:

$$C = \frac{1}{n - 1} X_{\text{centered}}^T X_{\text{centered}}$$

Notice key structural properties of $C$:

- **Symmetry:** $C_{ij} = C_{ji}$ because $\text{Cov}(X_i, X_j) = \text{Cov}(X_j, X_i)$, so $C^T = C$.
- **Positive Semi-Definite:** For any vector $v \in \mathbb{R}^d$, $v^T C v = \frac{1}{n-1} (X_{\text{centered}} v)^T (X_{\text{centered}} v) = \frac{1}{n-1} \|X_{\text{centered}} v\|^2 \ge 0$. Consequently, all eigenvalues of $C$ are real and non-negative ($\lambda_i \ge 0$).
- **Diagonal Entries:** $C_{jj}$ represents the sample variance of feature $j$ alone.

### Step 3: Variance Maximization via Lagrange Multipliers

Let $v \in \mathbb{R}^d$ be a candidate unit direction vector ($\|v\|^2 = v^T v = 1$). The projection of a centered sample point $x_i \in \mathbb{R}^d$ onto $v$ is the scalar dot product $x_i^T v$. The sample variance of these projected scalars across all $n$ points is:

$$\sigma_v^2 = \frac{1}{n - 1} \sum_{i=1}^n (x_i^T v)^2 = \frac{1}{n - 1} (X_{\text{centered}} v)^T (X_{\text{centered}} v) = v^T \left( \frac{1}{n - 1} X_{\text{centered}}^T X_{\text{centered}} \right) v = v^T C v$$

We wish to find the direction $v$ that maximizes the projected variance $v^T C v$, subject to the constraint that $v$ is a unit vector ($v^T v = 1$):

$$\max_{v} \; v^T C v \quad \text{subject to} \quad v^T v = 1$$

We formulate the Lagrangian function with Lagrange multiplier $\lambda$:

$$\mathcal{L}(v, \lambda) = v^T C v - \lambda (v^T v - 1)$$

Taking the vector derivative of $\mathcal{L}$ with respect to $v$ and setting it to zero:

$$\nabla_v \mathcal{L} = 2 C v - 2 \lambda v = 0 \implies C v = \lambda v$$

> **THE FUNDAMENTAL REVELATION:** The condition $C v = \lambda v$ is the textbook definition of an **eigenvector equation**! It proves that any stationary point of the projected variance must be an eigenvector of the covariance matrix $C$. Furthermore, substituting $C v = \lambda v$ back into the variance expression yields: $$\text{Var}(\text{projection}) = v^T C v = v^T (\lambda v) = \lambda (v^T v) = \lambda$$ Thus, **the eigenvalue $\lambda$ directly equals the amount of variance captured along that eigenvector**. To maximize variance, we simply pick the eigenvector associated with the largest eigenvalue $\lambda_1$.

### Step 4: Subsequent Components and the Projection Matrix

For the second principal component $v_2$, we maximize $v_2^T C v_2$ subject to two constraints: $v_2^T v_2 = 1$ and orthogonality to the first component ($v_2^T v_1 = 0$). By the spectral theorem for symmetric matrices, this yields the eigenvector corresponding to the second largest eigenvalue $\lambda_2$.

We sort the eigenvalues in descending order: $\lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_d \ge 0$. Selecting the top $k$ eigenvectors forms the projection matrix $W \in \mathbb{R}^{d \times k}$:

$$W = \begin{bmatrix} | & | & & | \\ v_1 & v_2 & \dots & v_k \\ | & | & & | \end{bmatrix}$$

The dimensionally reduced coordinates $Z \in \mathbb{R}^{n \times k}$ and the approximate reconstruction $\hat{X} \in \mathbb{R}^{n \times d}$ are given by:

$$Z = X_{\text{centered}} W, \qquad \hat{X} = Z W^T + \mu$$

### Step 5: Explained Variance Ratio and Cumulative Variance

Because the covariance matrix is symmetric, the sum of its diagonal elements (the total feature variance) equals the sum of its eigenvalues (the trace invariance property: $\text{Tr}(C) = \sum_{j=1}^d C_{jj} = \sum_{j=1}^d \lambda_j$). The fraction of total information preserved by component $i$ is:

$$\text{Explained Variance Ratio}_i = \frac{\lambda_i}{\sum_{j=1}^d \lambda_j}$$

## 5. Grounded Step-by-Step Hand-Worked Numerical Walkthrough

To eliminate all abstraction, let us execute the entire PCA pipeline on a verified $3$-sample, $2$-feature dataset with zero hand-waving. Consider the following pre-centered matrix $X_{\text{centered}}$:

$$X_{\text{centered}} = \begin{bmatrix} 2 & 1 \\ 0 & -1 \\ -2 & 0 \end{bmatrix}, \quad (n = 3, \; d = 2)$$

**Verification of zero mean:** Feature 1 mean is $\frac{2 + 0 + (-2)}{3} = 0$. Feature 2 mean is $\frac{1 + (-1) + 0}{3} = 0$.

### Sub-step 5.1: Compute the Sample Covariance Matrix C

With $n = 3$, degrees of freedom $n - 1 = 2$:

$$X^T X = \begin{bmatrix} 2 & 0 & -2 \\ 1 & -1 & 0 \end{bmatrix} \begin{bmatrix} 2 & 1 \\ 0 & -1 \\ -2 & 0 \end{bmatrix} = \begin{bmatrix} 2(2)+0(0)+(-2)(-2) & 2(1)+0(-1)+(-2)(0) \\ 1(2)+(-1)(0)+0(-2) & 1(1)+(-1)(-1)+0(0) \end{bmatrix} = \begin{bmatrix} 8 & 2 \\ 2 & 2 \end{bmatrix}$$

$$C = \frac{1}{n - 1} X^T X = \frac{1}{2} \begin{bmatrix} 8 & 2 \\ 2 & 2 \end{bmatrix} = \begin{bmatrix} 4 & 1 \\ 1 & 1 \end{bmatrix}$$

### Sub-step 5.2: Calculate the Eigenvalues

We find scalar values $\lambda$ where the characteristic polynomial $\det(C - \lambda I) = 0$:

$$\det \begin{bmatrix} 4 - \lambda & 1 \\ 1 & 1 - \lambda \end{bmatrix} = (4 - \lambda)(1 - \lambda) - (1)(1) = \lambda^2 - 5\lambda + 3 = 0$$

Applying the quadratic formula $\lambda = \frac{-b \pm \sqrt{b^2 - 4ac}}{2a}$:

$$\lambda = \frac{5 \pm \sqrt{25 - 4(1)(3)}}{2} = \frac{5 \pm \sqrt{13}}{2} \approx \frac{5 \pm 3.60555}{2}$$

$$\lambda_1 = \frac{5 + 3.60555}{2} \approx 4.3028, \qquad \lambda_2 = \frac{5 - 3.60555}{2} \approx 0.6972$$

**Trace Check:** Total variance $= \lambda_1 + \lambda_2 = 4.3028 + 0.6972 = 5.0000 = C_{11} + C_{22} = 4 + 1 = 5.0000$.

$$\text{EVR}_1 = \frac{4.3028}{5.0000} = 0.86056 \; (86.06\%), \qquad \text{EVR}_2 = \frac{0.6972}{5.0000} = 0.13944 \; (13.94\%)$$  A single principal component retains over 86% of the dataset's total variance!

### Sub-step 5.3: Compute the First Eigenvector v1

We solve $(C - \lambda_1 I) v_1 = 0$ for vector $v_1 = [v_{11}, v_{12}]^T$:

$$\begin{bmatrix} 4 - 4.3028 & 1 \\ 1 & 1 - 4.3028 \end{bmatrix} \begin{bmatrix} v_{11} \\ v_{12} \end{bmatrix} = \begin{bmatrix} -0.3028 & 1 \\ 1 & -3.3028 \end{bmatrix} \begin{bmatrix} v_{11} \\ v_{12} \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix}$$

From row 1: $-0.3028 v_{11} + v_{12} = 0 \implies v_{12} = 0.3028 v_{11}$. Imposing unit norm $\|v_1\|^2 = v_{11}^2 + v_{12}^2 = 1$:

$$v_{11}^2 + (0.3028 v_{11})^2 = v_{11}^2(1 + 0.09169) = 1.09169 v_{11}^2 = 1 \implies v_{11} = \frac{1}{\sqrt{1.09169}} \approx 0.9571$$

$$v_{12} = 0.3028 \times 0.9571 \approx 0.2898 \implies v_1 = \begin{bmatrix} 0.9571 \\ 0.2898 \end{bmatrix}$$

### Sub-step 5.4: Project the Data Points onto PC1

Projecting each 2D observation $x_i$ onto $v_1$ ($z_i = x_{i1} v_{11} + x_{i2} v_{12}$):

- **Sample 1 $[2, 1]$:** $z_1 = 2(0.9571) + 1(0.2898) = 1.9142 + 0.2898 = \mathbf{+2.2040}$
- **Sample 2 $[0, -1]$:** $z_2 = 0(0.9571) + (-1)(0.2898) = \mathbf{-0.2898}$
- **Sample 3 $[-2, 0]$:** $z_3 = -2(0.9571) + 0(0.2898) = \mathbf{-1.9142}$

> **SIGN INDETERMINACY IN EIGENVECTORS:** If you negate an eigenvector, $(-v)$ is equally valid: $C(-v) = -(Cv) = -(\lambda v) = \lambda (-v)$. Both $v_1$ and $-v_1$ define the exact same one-dimensional geometric axis line. Consequently, NumPy or Scikit-Learn might return projected values $[+2.204, -0.290, -1.914]$ or $[-2.204, +0.290, +1.914]$. This sign flip is normal and mathematically indistinguishable.

## 6. Production Engineering: Why Scikit-Learn Uses SVD Instead of Covariance

If PCA is defined by the eigenvectors of $C = \frac{1}{n-1} X^T X$, why does production code in `sklearn.decomposition.PCA` bypass computing $C$ entirely and run **Singular Value Decomposition (SVD)** on $X_{\text{centered}}$?

- **Condition Number Squaring:** The condition number $\kappa(M)$ measures numerical sensitivity to floating-point rounding errors. Forming the product $X^T X$ squares the condition number: $\kappa(X^T X) = (\kappa(X))^2$. If $X$ has an ill-conditioned ratio of singular values $\kappa(X) = 10^5$, $X^T X$ has $\kappa = 10^{10}$, causing severe loss of numerical precision.
- **Memory Footprint when $d \gg n$:** In genomic or NLP applications where $d = 50{,}000$ features but $n = 200$ patients, forming the covariance matrix requires allocating a massive $50{,}000 \times 50{,}000$ float matrix (10 GB RAM). SVD decomposes $X$ directly without allocating $C$.

Recall the thin SVD decomposition of centered data: $X_{\text{centered}} = U \Sigma V^T$, where $U \in \mathbb{R}^{n \times k}$, $\Sigma \in \mathbb{R}^{k \times k}$ is diagonal with singular values $\sigma_1 \ge \sigma_2 \ge \dots \ge \sigma_k$, and $V \in \mathbb{R}^{d \times k}$. Substituting this into the covariance matrix formulation:

$$C = \frac{1}{n - 1} X^T X = \frac{1}{n - 1} (U \Sigma V^T)^T (U \Sigma V^T) = \frac{1}{n - 1} V \Sigma U^T U \Sigma V^T = V \left( \frac{\Sigma^2}{n - 1} \right) V^T$$

Because $U^T U = I$, the right singular vectors $V$ of $X$ are **identically the eigenvectors of $C$**, and the eigenvalues relate directly to the singular values by: $$\lambda_i = \frac{\sigma_i^2}{n - 1}$$ Modern libraries compute SVD on $X$ directly using LAPACK routines (`gesdd`), achieving $O(n d \min(n, d))$ runtime with superior numerical stability.

## 7. Python Implementation From Scratch

Here is a complete, production-grade object-oriented implementation of PCA built using pure NumPy. It includes zero-mean centering, covariance computation, symmetric eigendecomposition via `np.linalg.eigh`, variance ratios, projection, and inverse reconstruction.

```python
import numpy as np

class PCAFromScratch:
    """
    Principal Component Analysis (PCA) implemented from first principles.
    
    Parameters
    ----------
    n_components : int or float
        If int, the number of top principal components to retain.
        If float between 0.0 and 1.0, the fraction of cumulative variance to preserve.
    """
    def __init__(self, n_components=2):
        self.n_components = n_components
        self.components_ = None            # Shape: (d, k) projection vectors
        self.explained_variance_ = None    # Retained eigenvalues
        self.explained_variance_ratio_ = None # Fraction of variance per component
        self.mean_ = None                  # Empirical feature means
        self.n_features_in_ = None

    def fit(self, X):
        X = np.asarray(X, dtype=np.float64)
        n_samples, n_features = X.shape
        self.n_features_in_ = n_features

        # Step 1: Center the data (mean = 0)
        self.mean_ = np.mean(X, axis=0)
        X_centered = X - self.mean_

        # Step 2: Sample covariance matrix C (d x d)
        # rowvar=False denotes columns are features, rows are samples
        cov_matrix = np.cov(X_centered, rowvar=False)
        if cov_matrix.ndim == 0:  # scalar fallback for 1D edge-case
            cov_matrix = np.array([[cov_matrix]])

        # Step 3: Eigendecomposition of symmetric matrix C
        # np.linalg.eigh is optimized for Hermitian/symmetric matrices
        eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix)

        # Step 4: Sort eigenvalues and eigenvectors in descending order
        idx_desc = np.argsort(eigenvalues)[::-1]
        eigenvalues = eigenvalues[idx_desc]
        eigenvectors = eigenvectors[:, idx_desc]

        # Zero out tiny negative eigenvalues caused by floating point precision
        eigenvalues = np.maximum(eigenvalues, 0.0)
        total_variance = np.sum(eigenvalues)
        variance_ratio = eigenvalues / total_variance if total_variance > 0 else np.zeros_like(eigenvalues)

        # Step 5: Select number of components
        if isinstance(self.n_components, float) and 0.0 < self.n_components <= 1.0:
            cumulative = np.cumsum(variance_ratio)
            k = int(np.searchsorted(cumulative, self.n_components) + 1)
        else:
            k = int(min(self.n_components, n_features))

        self.components_ = eigenvectors[:, :k]
        self.explained_variance_ = eigenvalues[:k]
        self.explained_variance_ratio_ = variance_ratio[:k]
        return self

    def transform(self, X):
        X = np.asarray(X, dtype=np.float64)
        X_centered = X - self.mean_
        return np.dot(X_centered, self.components_)

    def fit_transform(self, X):
        return self.fit(X).transform(X)

    def inverse_transform(self, Z):
        """Reconstruct high-dimensional approximations from latent scores: Z @ W.T + mean."""
        Z = np.asarray(Z, dtype=np.float64)
        return np.dot(Z, self.components_.T) + self.mean_


# --- Verification on the 3-sample hand-worked dataset ---
if __name__ == "__main__":
    X_demo = np.array([
        [ 2.0,  1.0],
        [ 0.0, -1.0],
        [-2.0,  0.0]
    ])

    # 1. Run From-Scratch PCA
    pca_scratch = PCAFromScratch(n_components=1)
    Z_scratch = pca_scratch.fit_transform(X_demo)
    X_reconstructed = pca_scratch.inverse_transform(Z_scratch)

    print("=== PCA FROM SCRATCH VERIFICATION ===")
    print(f"Eigenvalue (Variance) lambda_1: {pca_scratch.explained_variance_[0]:.4f}")
    print(f"Explained Variance Ratio:      {pca_scratch.explained_variance_ratio_[0]:.4%}")
    print(f"Principal Axis Direction (v1): {pca_scratch.components_[:, 0]}")
    print("Transformed 1D Coordinates (Z):")
    print(np.round(Z_scratch, 4))
    print("Reconstructed Approximation (X_hat):")
    print(np.round(X_reconstructed, 4))

    # 2. Compare Against Scikit-Learn
    from sklearn.decomposition import PCA
    pca_sklearn = PCA(n_components=1)
    Z_sklearn = pca_sklearn.fit_transform(X_demo)

    print("\n=== SCIKIT-LEARN VALIDATION ===")
    print(f"Sklearn Explained Variance:    {pca_sklearn.explained_variance_[0]:.4f}")
    print(f"Sklearn Variance Ratio:        {pca_sklearn.explained_variance_ratio_[0]:.4%}")
    print("Sklearn Transformed (Z):")
    print(np.round(Z_sklearn, 4))
    # Check absolute value congruence
    np.testing.assert_allclose(np.abs(Z_scratch), np.abs(Z_sklearn), rtol=1e-5)
    print("\n[SUCCESS] From-scratch implementation matches scikit-learn within 1e-5 precision!")
```

Executing this test harness confirms exact mathematical congruence with our manual calculations:

```python
=== PCA FROM SCRATCH VERIFICATION ===
Eigenvalue (Variance) lambda_1: 4.3028
Explained Variance Ratio:      86.0555%
Principal Axis Direction (v1): [-0.95709203 -0.28978415]
Transformed 1D Coordinates (Z):
[[-2.204 ]
 [ 0.2898]
 [ 1.9142]]
Reconstructed Approximation (X_hat):
[[ 2.1094  0.6387]
 [-0.2774 -0.084 ]
 [-1.832  -0.5547]]

=== SCIKIT-LEARN VALIDATION ===
Sklearn Explained Variance:    4.3028
Sklearn Variance Ratio:        86.0555%
Sklearn Transformed (Z):
[[ 2.204 ]
 [-0.2898]
 [-1.9142]]

[SUCCESS] From-scratch implementation matches scikit-learn within 1e-5 precision!
```

## 8. Selecting the Optimal Number of Principal Components

In exploratory data analysis and ML preprocessing, how do practitioners systematically decide how many components $k$ to retain? Three standard criteria are used:

1. **Cumulative Explained Variance Threshold (e.g., 90% or 95%):** Plot cumulative explained variance $\sum_{i=1}^k \text{EVR}_i$ against $k$. Choose the smallest $k$ where the curve surpasses 0.90 or 0.95. This guarantees that 95% of information is preserved while discarding the noise tail.
2. **The Scree Plot Elbow Method:** Plot individual eigenvalues $\lambda_i$ in descending order. The graph typically plunges steeply and then flattens out into an 'elbow.' The components before the elbow represent primary structural signal, while points after the elbow represent random experimental noise.
3. **The Kaiser-Guttman Criterion:** When PCA is applied to standardized data ($Z$-scores, where each variable has variance 1), any principal component with an eigenvalue $\lambda_i < 1.0$ contains less variance than a single original feature. The Kaiser rule states: retain only components where $\lambda_i \ge 1.0$.

## 9. Common Pitfalls, Edge Cases, and Failure Modes

- **Failing to Standardize Varied Scales:** If a real estate dataset pairs square footage ($800 \text{ to } 5{,}000$) with bedroom count ($1 \text{ to } 5$), the raw variance of square footage is tens of thousands of times greater. Without standard scaling, PC1 will align 99.9% with square footage simply due to units of measurement.
- **Data Leakage Across Train/Test Splits:** Never call `fit_transform()` on the entire dataset prior to splitting. Centering and computing covariance across test samples leaks future variance distributions into training. Always run `pca.fit(X_train)` and then `pca.transform(X_test)` using the training mean and components.
- **Failure on Non-Linear Manifolds:** PCA is strictly a _linear_ projection technique. If data lies on a curved manifold — such as concentric spheres, a Swiss roll, or non-linear parabolas — PCA flattens the geometry across chords, collapsing separate clusters on top of each other. In such domains, non-linear techniques (Kernel PCA, t-SNE, or UMAP) must be utilized.
- **Sensitivity to Extreme Outliers:** Because variance minimizes squared deviations ($d^2$), a single distant outlier can pull the principal component axis heavily toward itself, warping the coordinate frame for all normal points.
- **Sparse Text Data Collapse:** In NLP TF-IDF matrices (which are 99% sparse), centering ($X - \mu$) fills millions of zeros with non-zero mean values, immediately converting a sparse matrix into a dense matrix that consumes terabytes of memory. For sparse data, use `TruncatedSVD`, which decomposes without zero-centering.

## 10. Dimensionality Reduction Comparison Matrix

Understanding where PCA sits in relation to supervised, non-linear, and manifold learning methods:

| Algorithm | Supervised / Unsupervised | Linear / Non-linear | Primary Objective | Best Use Case |
| --- | --- | --- | --- | --- |
| **PCA** | Unsupervised | Linear | Maximizes global variance, minimizes reconstruction error | Fast preprocessing, feature compression, noise reduction, multicollinearity removal |
| **Kernel PCA** | Unsupervised | Non-linear | Projects data into high-dimensional Hilbert space via kernel trick ($k(x, y)$) | Separating non-linear geometric structures (circles, curves) for downstream linear models |
| **LDA** | Supervised (uses labels $y$) | Linear | Maximizes between-class variance relative to within-class variance ($S_B / S_W$) | Dimensionality reduction for linear classification (e.g., face recognition, cancer subtyping) |
| **t-SNE** | Unsupervised | Non-linear | Preserves local neighborhood probabilities using Student-t distributions | 2D/3D visualization of complex multi-cluster datasets (e.g., single-cell RNA, image embeddings) |
| **UMAP** | Unsupervised (or Semi-supervised) | Non-linear | Preserves local and global fuzzy simplicial set topology | High-performance non-linear visualization and general clustering preprocessing |
| **Autoencoder** | Self-supervised Neural Net | Non-linear (with non-linear activations) | Minimizes non-linear bottleneck reconstruction loss | Complex multimodal compression (images, audio, latent generative spaces) |

## 11. Hands-On Practice & Curriculum Roadmap

To consolidate your mastery of PCA and linear algebra in machine learning, tackle these practical challenges:

1. **Implement PCA with SVD:** Refactor the from-scratch class to compute `np.linalg.svd(X_centered, full_matrices=False)` and demonstrate that right singular vectors $V^T$ match the eigenvectors of the covariance matrix within numerical tolerance.
2. **The 95% Variance MNIST Compression:** Load Scikit-Learn's digits dataset ($8 \times 8 = 64$ dimensions). Plot the cumulative explained variance curve. Determine how many components are required to capture 95% of the variance, reconstruct the images, and visually compare the compressed digits against the originals.
3. **Explore Non-Linear Manifold Breakdown:** Generate a 3D Swiss roll using `sklearn.datasets.make_swiss_roll`. Apply standard PCA to project it to 2D, then compare the output against t-SNE and UMAP to visually witness why linear PCA fails on folded manifolds.

> **WHAT TO LEARN NEXT:** Now that you understand linear global dimensionality reduction, advance to non-linear manifold visualization. In the next guide, **t-SNE vs UMAP: The Math Behind High-Dimensional Visualization**, we deconstruct how stochastic neighbor embedding and Riemannian geometry resolve the limitations of PCA.

---

Written by Sythra — Learn machine learning by building. Practice this topic with Sythra's AI tutor: https://app.sythra.ai/pricing