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

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.

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:

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

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}$

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.

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:

=== 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:

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.