Eigenvalue Problem

We've established that principal components are eigenvectors of the covariance matrix. But what exactly is an eigenvector? Why do eigenvectors of symmetric matrices have such beautiful properties? And how do we actually compute them for real-world datasets with thousands of dimensions?

This page dives deep into the eigenvalue problem—the mathematical structure that makes PCA tractable and elegant. Understanding eigendecomposition isn't just academic; it provides insight into why PCA works, when it fails, and how its results should be interpreted.

We'll explore the spectral theorem for symmetric matrices, understand the computational algorithms behind the scenes, and develop intuition for what eigenvalues and eigenvectors truly represent.

Let's start with the fundamental definitions and build intuition for what these objects represent.

The Definition

Given a square matrix $\mathbf{A} \in \mathbb{R}^{d \times d}$, a scalar $\lambda$ and a non-zero vector $\mathbf{v}$ are an eigenvalue-eigenvector pair if:

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

This seemingly simple equation says something profound: when $\mathbf{A}$ acts on $\mathbf{v}$, the result is just a scaled version of $\mathbf{v}$.

Most vectors get 'rotated' and 'stretched' by a matrix—their direction changes. Eigenvectors are special: they only get scaled, remaining along the same line through the origin.

Geometric Interpretation

Imagine a linear transformation represented by matrix $\mathbf{A}$. This transformation:

• Warps space: circles become ellipses, squares become parallelograms
• Rotates, shears, scales, and possibly reflects

But along eigenvector directions, something simple happens: the transformation is pure scaling.

• If $\lambda > 1$: the eigenvector direction is stretched
• If $0 < \lambda < 1$: the eigenvector direction is compressed
• If $\lambda < 0$: the eigenvector direction is reflected and scaled
• If $\lambda = 0$: vectors in that direction are collapsed to zero

Finding Eigenvalues: The Characteristic Equation

To find eigenvalues, we rearrange:

$$\mathbf{A}\mathbf{v} = \lambda\mathbf{v} \implies (\mathbf{A} - \lambda\mathbf{I})\mathbf{v} = \mathbf{0}$$

For a non-zero solution $\mathbf{v}$ to exist, the matrix $(\mathbf{A} - \lambda\mathbf{I})$ must be singular (not invertible). This happens when:

$$\det(\mathbf{A} - \lambda\mathbf{I}) = 0$$

This is the characteristic equation. It's a polynomial of degree $d$ in $\lambda$, so it has exactly $d$ roots (counting multiplicity, in the complex numbers).

For a $d \times d$ matrix, there are $d$ eigenvalues (some may be repeated, and for general matrices, some may be complex).

The covariance matrix has a special structure: it is symmetric ($\mathbf{S} = \mathbf{S}^T$). This symmetry grants us a powerful theorem that makes PCA particularly elegant.

The Spectral Theorem

Theorem: Let $\mathbf{S}$ be a real symmetric matrix. Then:

1. All eigenvalues are real (not complex)
2. Eigenvectors corresponding to distinct eigenvalues are orthogonal
3. $\mathbf{S}$ can be diagonalized by an orthogonal matrix

More precisely, there exists an orthogonal matrix $\mathbf{W}$ (meaning $\mathbf{W}^T = \mathbf{W}^{-1}$) such that:

$$\mathbf{S} = \mathbf{W} \mathbf{\Lambda} \mathbf{W}^T$$

where $\mathbf{\Lambda} = \text{diag}(\lambda_1, \lambda_2, \ldots, \lambda_d)$ is a diagonal matrix of eigenvalues.

Proof Sketch: Real Eigenvalues

Let $\mathbf{S}\mathbf{v} = \lambda\mathbf{v}$ where $\mathbf{v} eq \mathbf{0}$. Consider:

$$\mathbf{v}^T \mathbf{S} \mathbf{v} = \mathbf{v}^T (\lambda \mathbf{v}) = \lambda \|\mathbf{v}\|^2$$

Also, since $\mathbf{S}$ is symmetric and taking the transpose:

$$\mathbf{v}^T \mathbf{S} \mathbf{v} = (\mathbf{v}^T \mathbf{S} \mathbf{v})^T = \mathbf{v}^T \mathbf{S}^T \mathbf{v} = \mathbf{v}^T \mathbf{S} \mathbf{v}$$

This scalar equals its own transpose, so it's real. Since $\|\mathbf{v}\|^2 > 0$, we have $\lambda = \mathbf{v}^T \mathbf{S} \mathbf{v} / \|\mathbf{v}\|^2$ is real.

(For complex eigenvectors, a more careful argument using Hermitian adjoints is needed, but the conclusion holds: symmetric real matrices have only real eigenvalues.)

Proof Sketch: Orthogonal Eigenvectors

Let $\mathbf{S}\mathbf{v}_i = \lambda_i \mathbf{v}_i$ and $\mathbf{S}\mathbf{v}_j = \lambda_j \mathbf{v}_j$ with $\lambda_i eq \lambda_j$.

Compute $\mathbf{v}_j^T \mathbf{S} \mathbf{v}_i$ two ways:

$$\mathbf{v}_j^T \mathbf{S} \mathbf{v}_i = \mathbf{v}_j^T (\lambda_i \mathbf{v}_i) = \lambda_i (\mathbf{v}_j^T \mathbf{v}_i)$$

$$\mathbf{v}_j^T \mathbf{S} \mathbf{v}_i = (\mathbf{S} \mathbf{v}_j)^T \mathbf{v}_i = (\lambda_j \mathbf{v}_j)^T \mathbf{v}_i = \lambda_j (\mathbf{v}_j^T \mathbf{v}_i)$$

Equating: $\lambda_i (\mathbf{v}_j^T \mathbf{v}_i) = \lambda_j (\mathbf{v}_j^T \mathbf{v}_i)$

Since $\lambda_i eq \lambda_j$, we must have $\mathbf{v}_j^T \mathbf{v}_i = 0$. QED.

The spectral theorem guarantees that the covariance matrix can be decomposed as:

$$\mathbf{S} = \mathbf{W} \mathbf{\Lambda} \mathbf{W}^T = \sum_{j=1}^{d} \lambda_j \mathbf{w}_j \mathbf{w}_j^T$$

Let's unpack this decomposition and understand what it tells us.

The Outer Product Form

The sum form is particularly illuminating:

$$\mathbf{S} = \lambda_1 \mathbf{w}_1 \mathbf{w}_1^T + \lambda_2 \mathbf{w}_2 \mathbf{w}_2^T + \cdots + \lambda_d \mathbf{w}_d \mathbf{w}_d^T$$

Each term $\mathbf{w}_j \mathbf{w}_j^T$ is a projection matrix onto the line spanned by $\mathbf{w}_j$. The covariance matrix is a weighted sum of projection matrices, where the weights are the eigenvalues.

Low-Rank Approximation

If we truncate the sum after $k$ terms:

$$\hat{\mathbf{S}}k = \sum{j=1}^{k} \lambda_j \mathbf{w}_j \mathbf{w}_j^T$$

This is the best rank-$k$ approximation to $\mathbf{S}$ in the Frobenius norm. The approximation error is:

$$\|\mathbf{S} - \hat{\mathbf{S}}_k\|F = \sqrt{\sum{j=k+1}^{d} \lambda_j^2}$$

This is the Eckart-Young-Mirsky theorem applied to the covariance matrix.

The Inverse and Pseudoinverse

The eigendecomposition also reveals matrix inverses:

$$\mathbf{S}^{-1} = \mathbf{W} \mathbf{\Lambda}^{-1} \mathbf{W}^T = \sum_{j=1}^{d} \frac{1}{\lambda_j} \mathbf{w}_j \mathbf{w}_j^T$$

This only works if all $\lambda_j > 0$. If some eigenvalues are zero (rank-deficient case), the matrix is singular and has no inverse.

The pseudoinverse handles this by only inverting non-zero eigenvalues:

$$\mathbf{S}^{+} = \sum_{j: \lambda_j > 0} \frac{1}{\lambda_j} \mathbf{w}_j \mathbf{w}_j^T$$

This is useful in high-dimensional settings where $n < d$.

The covariance matrix has an additional special property beyond symmetry: it is positive semi-definite (PSD). This has important implications for its eigenvalues and for PCA.

Definition of PSD

A symmetric matrix $\mathbf{S}$ is positive semi-definite if:

$$\mathbf{v}^T \mathbf{S} \mathbf{v} \geq 0 \quad \text{for all } \mathbf{v} \in \mathbb{R}^d$$

For the covariance matrix, this is easily verified:

$$\mathbf{v}^T \mathbf{S} \mathbf{v} = \mathbf{v}^T \left(\frac{1}{n}\mathbf{X}^T\mathbf{X}\right) \mathbf{v} = \frac{1}{n}\|\mathbf{X}\mathbf{v}\|^2 \geq 0$$

The squared norm is always non-negative, so the covariance matrix is always PSD.

Implications for Eigenvalues

Positive semi-definiteness is equivalent to all eigenvalues being non-negative:

$$\lambda_j \geq 0 \quad \text{for all } j$$

Proof: If $\mathbf{S}\mathbf{v} = \lambda\mathbf{v}$ with $\|\mathbf{v}\| = 1$, then:

$$\lambda = \mathbf{v}^T \mathbf{S} \mathbf{v} \geq 0$$

by the PSD property.

This is crucial for PCA: eigenvalues represent variances, and variances cannot be negative. The mathematics confirms what intuition demands.

Computing the characteristic polynomial and finding its roots is mathematically exact but computationally impractical for large matrices. Real-world eigensolvers use iterative methods that are more stable and efficient.

Power Iteration

The simplest iterative method for finding the largest eigenvalue and its eigenvector:

1. Start with a random vector $\mathbf{v}^{(0)}$
2. Repeat: $\mathbf{v}^{(k+1)} = \frac{\mathbf{S}\mathbf{v}^{(k)}}{\|\mathbf{S}\mathbf{v}^{(k)}\|}$
3. Under mild conditions, $\mathbf{v}^{(k)} \to \mathbf{w}_1$

The eigenvalue is then $\lambda_1 = (\mathbf{v}^{(k)})^T \mathbf{S} \mathbf{v}^{(k)}$.

Why it works: Decompose $\mathbf{v}^{(0)} = \sum_j c_j \mathbf{w}_j$. After $k$ iterations:

$$\mathbf{S}^k \mathbf{v}^{(0)} = \sum_j c_j \lambda_j^k \mathbf{w}_j = \lambda_1^k \left(c_1 \mathbf{w}1 + \sum{j>1} c_j \left(\frac{\lambda_j}{\lambda_1}\right)^k \mathbf{w}_j\right)$$

Since $|\lambda_j / \lambda_1| < 1$ for $j > 1$, the other components decay exponentially.

QR Algorithm

The workhorse algorithm for computing all eigenvalues simultaneously:

1. Initialize $\mathbf{A}_0 = \mathbf{S}$
2. Repeat:
   • Compute QR decomposition: $\mathbf{A}_k = \mathbf{Q}_k \mathbf{R}_k$
   • Update: $\mathbf{A}_{k+1} = \mathbf{R}_k \mathbf{Q}_k$
3. $\mathbf{A}_k$ converges to diagonal (or block-diagonal) form

The diagonal entries converge to eigenvalues. For symmetric matrices, convergence is particularly fast.

With shifts: Adding shifts $\mathbf{A}_k - \mu_k \mathbf{I}$ before QR decomposition accelerates convergence dramatically. The Wilkinson shift strategy achieves cubic convergence.

Singular Value Decomposition (SVD) is intimately related to PCA and often preferred computationally. Understanding this connection is essential for practical PCA implementation.

SVD Definition

Any matrix $\mathbf{X} \in \mathbb{R}^{n \times d}$ can be decomposed as:

$$\mathbf{X} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^T$$

where:

• $\mathbf{U}$ is $n \times n$ orthogonal (left singular vectors)
• $\mathbf{\Sigma}$ is $n \times d$ diagonal (singular values $\sigma_1 \geq \sigma_2 \geq \cdots \geq 0$)
• $\mathbf{V}$ is $d \times d$ orthogonal (right singular vectors)

Connecting SVD to Eigendecomposition

Consider the covariance matrix (assuming centered $\mathbf{X}$):

$$\mathbf{S} = \frac{1}{n}\mathbf{X}^T\mathbf{X} = \frac{1}{n}(\mathbf{U}\mathbf{\Sigma}\mathbf{V}^T)^T(\mathbf{U}\mathbf{\Sigma}\mathbf{V}^T)$$

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

Comparing with $\mathbf{S} = \mathbf{W}\mathbf{\Lambda}\mathbf{W}^T$:

• Principal component directions: $\mathbf{W} = \mathbf{V}$ (right singular vectors)
• Eigenvalues: $\lambda_j = \sigma_j^2 / n$ (squared singular values divided by $n$)

Why SVD is Preferred

1. Numerical Stability: Forming $\mathbf{X}^T\mathbf{X}$ squares the condition number, amplifying numerical errors. SVD works directly on $\mathbf{X}$.

2. Memory Efficiency: For $n < d$, we never need to form the $d \times d$ covariance matrix.

3. Works on Rectangular Matrices: SVD is defined for any matrix, not just square ones.

4. Truncated Computation: Efficient algorithms compute only the top $k$ singular values/vectors without computing all $d$.

Real-world PCA implementation must handle various numerical challenges. Understanding these helps you interpret results correctly and avoid common mistakes.

Condition Number

The condition number $\kappa(\mathbf{S}) = \lambda_{\max} / \lambda_{\min}$ measures sensitivity to numerical errors.

• Small $\kappa$ (close to 1): Well-conditioned, stable computation
• Large $\kappa$ (> 10^6): Ill-conditioned, results may be unreliable
• $\kappa = \infty$: Matrix is singular (has zero eigenvalues)

Ill-conditioning often occurs when features have vastly different scales or when features are highly correlated.

Handling Near-Singular Matrices

When eigenvalues are very small, practical approaches include:

1. Regularization: Add a small value to diagonal: $\mathbf{S} + \epsilon \mathbf{I}$
2. Truncation: Only keep components with $\lambda > \tau$ for some threshold $\tau$
3. Variance Threshold: Keep components explaining at least $x\%$ of variance
4. Pseudoinverse: Use Moore-Penrose pseudoinverse for downstream computations