PCA: Eigendecomposition

Post 03 showed that PCA finds the direction u that maximizes ${\mathbf{u}}^{T}C\mathbf{u}$, and that a diagonal direction outperforms the original axes. This post computes the exact optimal direction by solving the eigendecomposition of the covariance matrix — by hand, step by step.

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}")

What is an Eigenvector?

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).

Step 1: Find Eigenvalues from 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 2: Find Eigenvectors

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 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{\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 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)

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.

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}")

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)}")

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

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?