PCA: Eigendecomposition

The diagonal direction $[1/\sqrt{2},1/\sqrt{2}]$ from post 03 captures 5.67 variance — beating either original axis. But 5.67 is not the maximum. The true optimal direction captures 5.98 — and finding it requires a tool we haven't used yet.

That tool is the eigendecomposition of the covariance matrix. An eigenvector of C is a direction that C only stretches, never rotates. The amount of stretch (the eigenvalue) equals the variance captured. The eigenvector with the largest eigenvalue is PC1 — the single direction that captures the most possible variance. This post computes both eigenvalues and eigenvectors by hand on the same 4-sample dataset, then connects them to the SVD that sklearn uses under the hood.

Anchor: Continuing the 4-sample 2D dataset from post 03 (hours_studied vs exam_score).

import numpy as np

# Centered data from post 03
X_c = np.array([
    [-1.5, -2.5],
    [-0.5, -0.5],
    [ 0.5,  0.5],
    [ 1.5,  2.5],
])

# Covariance matrix from post 03
C = np.array([[1.667, 2.667],
              [2.667, 4.333]])

print("Covariance matrix C:")
print(C)
print(f"trace(C) = {np.trace(C):.3f}")
print(f"det(C)   = {np.linalg.det(C):.4f}")

Covariance matrix C:
[[1.667 2.667]
 [2.667 4.333]]
trace(C) = 6.000
det(C)   = 0.1112

How an eigenvector captures the direction of maximum stretch

A vector v is an eigenvector of matrix C if multiplying by C only stretches it — it doesn't rotate:

$C\mathbf{v}=\lambda \mathbf{v}$

$\lambda$ is the eigenvalue: the scalar stretch factor. For PCA:

• Eigenvectors of C give the principal component directions
• Eigenvalues give the variance captured by each PC

The eigenvector with the largest eigenvalue = PC1 (most variance). Smallest = last PC (least variance).

The Plan — Four Steps to PCA via Eigendecomposition

We'll compute PCA from the covariance matrix in four steps. First we find the eigenvalues (how much variance each PC captures) and eigenvectors (the PC directions themselves). Then we verify they satisfy the eigen-equation. Finally we project the centered data onto the new axes to get PC scores.

Step 1: Find Eigenvalues from the Characteristic Equation

Now we compute the eigenvalues by solving the characteristic equation:

$\det (C-\lambda I)=0$

Substituting:

$C-\lambda I=\left[\begin{array}{cc}1.667-\lambda & 2.667\\ 2.667& 4.333-\lambda \end{array}\right]$

$\det (C-\lambda I)=(1.667-\lambda )(4.333-\lambda )-2.667^2$

Expanding:

$={\lambda }^2-6\lambda +(1.667\times 4.333-7.112)$ $={\lambda }^2-6\lambda +(7.222-7.112)$ $={\lambda }^2-6\lambda +0.111=0$

Quadratic formula:

$\lambda =\frac{6\pm \sqrt{36-4\times 0.111}}{2}=\frac{6\pm \sqrt{35.556}}{2}=\frac{6\pm 5.963}{2}$

${\lambda }_1=\frac{6+5.963}{2}=5.981{\lambda }_2=\frac{6-5.963}{2}=0.019$

Sanity checks:

${\lambda }_1+{\lambda }_2=5.981+0.019=6.000=\text{trace}(C)=1.667+4.333✓$

${\lambda }_1\times {\lambda }_2=5.981\times 0.019=0.114\approx \det (C)=0.111✓$

The sum of eigenvalues equals the total variance in the data. The eigenvalues partition this total.

✓ Step 1 complete. λ₁ = 5.981 captures 99.68% of total variance; λ₂ = 0.019 captures the remaining 0.32%.

Step 2: Find Eigenvectors

Now we find the eigenvector for each eigenvalue, starting with the one that captures the most variance.

PC1 direction (λ₁ = 5.981)

Solve $(C-{\lambda }_{1}I)\mathbf{v}=0$:

$C-5.981I=\left[\begin{array}{cc}1.667-5.981& 2.667\\ 2.667& 4.333-5.981\end{array}\right]=\left[\begin{array}{cc}-4.314& 2.667\\ 2.667& -1.648\end{array}\right]$

From row 1: $-4.314v_1+2.667v_2=0⟹v_2=\frac{4.314}{2.667}{v}_1=1.618v_1$

Set $v_1=1$: $\mathbf{v}=[1,1.618]$. Normalize: $\Vert \mathbf{v}\Vert =\sqrt{1^2+1.618^2}=\sqrt{3.618}=1.902$

${\mathbf{v}}_1=\left[\frac{1}{1.902},\frac{1.618}{1.902}\right]=[0.526,0.851]$

PC2 direction (λ₂ = 0.019)

From row 1 of $(C-0.019I)$: $1.648v_1+2.667v_2=0⟹v_2=-0.618v_1$

Normalized: ${\mathbf{v}}_2=[0.851,-0.526]$

✓ Step 2 complete. v₁ = [0.526, 0.851] is the direction of maximum variance (PC1); v₂ = [0.851, −0.526] is the orthogonal second component (PC2).

Step 3: Verify

# Verify: C × v1 = λ1 × v1
v1 = np.array([0.526, 0.851])
v2 = np.array([0.851, -0.526])
lam1, lam2 = 5.981, 0.019

print("C × v1:", (C @ v1).round(4))
print("λ1 × v1:", (lam1 * v1).round(4))
print()
print("C × v2:", (C @ v2).round(4))
print("λ2 × v2:", (lam2 * v2).round(4))
print()
print("v1 · v2 (should be 0):", np.dot(v1, v2).round(6))
print("|v1| =", np.linalg.norm(v1).round(4), "  |v2| =", np.linalg.norm(v2).round(4))

C × v1: [3.1498 5.0918]
λ1 × v1: [3.1462 5.0905]

C × v2: [0.0161 -0.0099]
λ2 × v2: [0.0162 -0.0100]

v1 · v2 (should be 0): 0.0
|v1| = 1.0   |v2| = 1.0

$C{\mathbf{v}}_1\approx {\lambda }_1{\mathbf{v}}_1$ — confirmed within rounding. The two eigenvectors are orthogonal (${\mathbf{v}}_1\cdot {\mathbf{v}}_2=0$) and unit length — guaranteed for symmetric matrices by the Spectral Theorem.

✓ Step 3 complete. Cv₁ ≈ λ₁v₁ and Cv₂ ≈ λ₂v₂ confirmed within rounding. Eigenvectors satisfy both the eigen-equation and the orthogonality check.

Step 4: Project Data onto Principal Components

Scores $=X_c\cdot V$ where $V=[{\mathbf{v}}_1{\mathbf{v}}_2]$ (columns are eigenvectors):

V = np.column_stack([v1, v2])  # (2,2) matrix
scores = X_c @ V               # (4,2) scores matrix
print("PC1 scores:", scores[:, 0].round(4))
print("PC2 scores:", scores[:, 1].round(4))
print()
print("Var(PC1 scores):", scores[:, 0].var(ddof=1).round(4), " ≈ λ₁ =", lam1)
print("Var(PC2 scores):", scores[:, 1].var(ddof=1).round(4), " ≈ λ₂ =", lam2)

PC1 scores: [-2.917  -0.6888  0.6888  2.917]
PC2 scores: [ 0.0379 -0.1624  0.1624 -0.0379]

Var(PC1 scores): 5.9827  ≈ λ₁ = 5.981
Var(PC2 scores): 0.0173  ≈ λ₂ = 0.019

The variance of each set of scores equals the corresponding eigenvalue — this is not a coincidence, it's the definition. Eigenvalues ARE the variances captured by each PC.

✓ Step 4 complete. Var(PC1 scores) = 5.983 and Var(PC2 scores) = 0.017, matching λ₁ and λ₂. The eigendecomposition is verified end-to-end on the anchor data.

Explained Variance Ratio

${\text{EVR}}_1=\frac{{\lambda }_1}{{\lambda }_1+{\lambda }_2}=\frac{5.981}{6.000}=0.9968=99.68\%$

${\text{EVR}}_2=\frac{{\lambda }_2}{{\lambda }_1+{\lambda }_2}=\frac{0.019}{6.000}=0.0032=0.32\%$

PC1 accounts for 99.68% of variance. PC2 is almost purely noise. We can project to 1D with near-zero information loss.

Reconstruction Quality

# Reconstruct using only PC1
pc1_scores = scores[:, 0]          # (4,)
reconstructed = np.outer(pc1_scores, v1)  # (4,2)

errors = np.sqrt(((X_c - reconstructed)**2).sum(axis=1))
print("Reconstruction from PC1 only:")
for i, (true, recon, err) in enumerate(zip(X_c, reconstructed, errors)):
    print(f"  Sample {i+1}: true={true}, recon={recon.round(3)}, error={err:.4f}")

Reconstruction from PC1 only:
  Sample 1: true=[-1.5 -2.5], recon=[-1.534 -2.483], error=0.0380
  Sample 2: true=[-0.5 -0.5], recon=[-0.362 -0.586], error=0.1623
  Sample 3: true=[ 0.5  0.5], recon=[ 0.362  0.586], error=0.1623
  Sample 4: true=[ 1.5  2.5], recon=[ 1.534  2.483], error=0.0380

Compare to post 03: reconstruction error with the diagonal direction u₃ was 0.707. With the true eigenvector v₁, it's only 0.038 — 18× smaller. Eigenvectors genuinely are the optimal projection directions.

SVD — The Numerically Stable Alternative

In practice, sklearn doesn't compute the covariance matrix at all. It uses Singular Value Decomposition (SVD) directly on $X_c$:

$X_c=U\Sigma V^T$

• $U$ (n×n): left singular vectors (one per sample)
• $\Sigma$ (n×p): diagonal matrix of singular values ${\sigma }_i$
• $V$ (p×p): right singular vectors = principal component directions

The connection to eigendecomposition: since $C=X_c^T{X}_c/(n-1)=V({\Sigma }^2/(n-1))V^T$, the right singular vectors of $X_c$ are exactly the eigenvectors of $C$, and:

${\lambda }_i=\frac{{\sigma }_i^2}{n-1}$

U, singular_values, Vt = np.linalg.svd(X_c, full_matrices=False)
eigenvalues_from_svd = singular_values**2 / (len(X_c) - 1)

print("Eigenvalues from characteristic equation: [5.981, 0.019]")
print(f"Eigenvalues from SVD:                      {eigenvalues_from_svd.round(4)}")
print()
print("PC1 from hand computation: [0.526, 0.851]")
print(f"PC1 from SVD (Vt[0]):      {np.abs(Vt[0]).round(4)}")

Eigenvalues from characteristic equation: [5.981, 0.019]
Eigenvalues from SVD:                      [5.9815 0.0185]

PC1 from hand computation: [0.526, 0.851]
PC1 from SVD (Vt[0]):      [0.5245 0.8514]

Identical within rounding. Why use SVD over direct eigendecomposition? Computing $X_c^T{X}_c$ squares the condition number, amplifying numerical errors for nearly-degenerate matrices. SVD avoids this by working directly on $X_c$.

Why Eigenvectors of C are Orthogonal

The covariance matrix $C$ is symmetric ($C^T=C$). The Spectral Theorem guarantees:

1. All eigenvalues are real (no complex numbers)
2. Eigenvectors for distinct eigenvalues are orthogonal
3. $C$ is diagonalizable: $C=V\text{diag}(\lambda )V^T$

Verified numerically above: ${\mathbf{v}}_1\cdot {\mathbf{v}}_2=0.526\times 0.851+0.851\times (-0.526)=0.4477-0.4477=0$.

This orthogonality is why PC scores are uncorrelated — the new axes are constructed to be independent directions.

Sign Convention

Eigenvectors satisfy $C\mathbf{v}=\lambda \mathbf{v}$ but also $C(-\mathbf{v})=\lambda (-\mathbf{v})$. Both $+{\mathbf{v}}_1$ and $-{\mathbf{v}}_1$ are valid eigenvectors. Sklearn normalizes so the component with the largest absolute loading is positive — but you may see PC1 flip sign between runs or sklearn versions. The scores flip sign accordingly; distances and EVR are unaffected.

Trace Table

Related Concepts

Backward: This post depends on the geometric intuition from the previous post (PCA as maximizing variance along orthogonal directions). The covariance matrix computation, eigen-decomposition basics, and the concept of explained variance ratio are all used here. Understanding SVD as an alternative to eigendecomposition for numerical stability is helpful but not required.

Forward: The PCA implementation post builds the full sklearn pipeline using these eigen-decomposition concepts. The comparison with LDA shows how supervised dimensionality reduction differs (maximizing between-class scatter instead of total variance). For non-negative data, NMF provides an alternative decomposition with parts-based interpretation.

Honest Limitations

1. Eigendecomposition is numerically unstable for near-singular matrices. When features are highly correlated or one feature is a linear combination of others (rank-deficient), the covariance matrix is near-singular and eigenvalues become unreliable. Sklearn's PCA uses SVD instead for precisely this reason — SVD operates on the data matrix directly without computing $X^{T}X$.
2. Eigenvectors are not uniquely signed. Both $+{\mathbf{v}}_1$ and $-{\mathbf{v}}_1$ are valid eigenvectors. Sklearn normalizes so the largest absolute loading is positive, but PC1 scores can flip sign between runs or versions. Distances and EVR are unaffected, but interpretations of "positive loading" require caution.
3. PCA is not feature selection. Each principal component is a linear combination of all original features. The loadings show contribution weight, not feature importance. If interpretability of individual features matters (e.g., medical diagnosis where each feature has a clear clinical meaning), L1-penalized methods (sparse PCA, Lasso) are more appropriate.

Test Your Understanding

1. The eigenvalue ${\lambda }_1=5.981$ and ${\lambda }_2=0.019$ sum to 6.000 = trace(C). The diagonal entries of C are the individual feature variances (1.667 and 4.333). What does it mean that trace is preserved under eigendecomposition — specifically, what geometric quantity is conserved when you rotate the coordinate system?

2. We normalized the eigenvector by dividing by $\Vert \mathbf{v}\Vert =1.902$. What would happen to the PC1 scores if we had used the unnormalized eigenvector [1, 1.618] instead? Would the explained variance ratio change?

3. The reconstruction error for sample 1 using the true eigenvector v₁ was 0.038, compared to 0.707 with the diagonal direction u₃ from post 03. Calculate: what fraction of the total reconstruction error (across all 4 samples) is explained by the variance captured by PC2 (${\lambda }_2=0.019$)?

4. For a 3×3 covariance matrix, the characteristic equation $\det (C-\lambda I)=0$ becomes a cubic — three eigenvalues. The eigenvalues of a PSD (positive semi-definite) matrix are non-negative. Why must all eigenvalues of a covariance matrix be $\ge 0$ — what geometric interpretation forbids a negative eigenvalue?

5. Sklearn uses SVD instead of eigendecomposition because computing $X_c^T{X}_c$ amplifies numerical errors. Construct a simple example where this matters: what would happen if you added a third feature to X that is exactly equal to 0.1 × feature_1? Would the covariance matrix become degenerate, and how would eigendecomposition fail while SVD would not?