Linear Algebra: Eigenvalues and Eigenvectors Explained

Key Insights

Eigenvectors are special directions in space that remain unchanged (only scaled) when a matrix transformation is applied, with eigenvalues representing the scaling factor
Computing eigenvalues requires solving the characteristic equation det(A - λI) = 0, though production code should use optimized libraries like NumPy rather than manual calculation
Eigendecomposition powers critical algorithms including PCA for dimensionality reduction, Google’s PageRank, and stability analysis in control systems

What Are Eigenvalues and Eigenvectors?

When you apply a matrix transformation to most vectors, both their direction and magnitude change. Eigenvectors are the exceptional cases—vectors that maintain their direction under the transformation, only getting scaled by a factor. That scaling factor is the eigenvalue.

Mathematically, for a square matrix A, vector v is an eigenvector with eigenvalue λ if:

Av = λv

Think of stretching a piece of fabric. Most points move in complex ways, but there are special directions—the warp and weft threads—that only stretch or compress along their original direction. These are analogous to eigenvectors.

Here’s a simple example showing how eigenvectors preserve direction:

import numpy as np
import matplotlib.pyplot as plt

# Define a 2x2 matrix
A = np.array([[3, 1],
              [0, 2]])

# A vector that is NOT an eigenvector
v1 = np.array([1, 1])
Av1 = A @ v1

# An eigenvector (direction [1, 0])
v2 = np.array([1, 0])
Av2 = A @ v2

print(f"Non-eigenvector: {v1} -> {Av1}")
print(f"Direction changed: {Av1 / np.linalg.norm(Av1)}")
print(f"\nEigenvector: {v2} -> {Av2}")
print(f"Direction preserved: {Av2 / np.linalg.norm(Av2)}")

The eigenvector maintains its direction (1, 0), while the non-eigenvector changes direction entirely.

The Mathematical Foundation

To find eigenvalues, we rearrange the eigenvalue equation:

Av = λv
Av - λv = 0
(A - λI)v = 0

For non-trivial solutions (v ≠ 0), the matrix (A - λI) must be singular, meaning its determinant equals zero:

det(A - λI) = 0

This is the characteristic equation. Solving it yields the eigenvalues, which we then substitute back to find corresponding eigenvectors.

Let’s walk through a complete calculation:

import numpy as np
from numpy.linalg import eig

# Define matrix
A = np.array([[4, 2],
              [1, 3]])

print("Matrix A:")
print(A)

# Manual calculation of characteristic polynomial
# det(A - λI) = det([[4-λ, 2], [1, 3-λ]])
# = (4-λ)(3-λ) - 2*1
# = λ² - 7λ + 10 = 0
# Solutions: λ = 5, λ = 2

print("\nManual calculation:")
print("Characteristic polynomial: λ² - 7λ + 10 = 0")
print("Eigenvalues: λ₁ = 5, λ₂ = 2")

# For λ = 5: (A - 5I)v = 0
# [[-1, 2], [1, -2]]v = 0
# Eigenvector: v₁ = [2, 1]

# For λ = 2: (A - 2I)v = 0
# [[2, 2], [1, 1]]v = 0
# Eigenvector: v₂ = [1, -1]

# NumPy calculation
eigenvalues, eigenvectors = eig(A)
print("\nNumPy calculation:")
print(f"Eigenvalues: {eigenvalues}")
print(f"Eigenvectors:\n{eigenvectors}")

# Verify: Av = λv
for i in range(len(eigenvalues)):
    v = eigenvectors[:, i]
    lambda_i = eigenvalues[i]
    print(f"\nVerification for λ={lambda_i:.2f}:")
    print(f"Av = {A @ v}")
    print(f"λv = {lambda_i * v}")

Computing Eigenvalues and Eigenvectors

While manual calculation works for small matrices, real applications require efficient algorithms. The power iteration method is a fundamental technique for finding the dominant eigenvalue (largest in magnitude):

def power_iteration(A, num_iterations=100):
    """Find dominant eigenvalue and eigenvector using power iteration."""
    n = A.shape[0]
    # Start with random vector
    v = np.random.rand(n)
    v = v / np.linalg.norm(v)
    
    for _ in range(num_iterations):
        # Multiply by A
        v_new = A @ v
        # Normalize
        v_new = v_new / np.linalg.norm(v_new)
        v = v_new
    
    # Compute eigenvalue using Rayleigh quotient
    eigenvalue = (v.T @ A @ v) / (v.T @ v)
    return eigenvalue, v

# Test on a matrix
A = np.array([[4, 2],
              [1, 3]])

lambda_power, v_power = power_iteration(A)
print(f"Power iteration result:")
print(f"Eigenvalue: {lambda_power:.4f}")
print(f"Eigenvector: {v_power}")

# Compare with NumPy
eigenvalues, eigenvectors = np.linalg.eig(A)
dominant_idx = np.argmax(np.abs(eigenvalues))
print(f"\nNumPy result (dominant):")
print(f"Eigenvalue: {eigenvalues[dominant_idx]:.4f}")
print(f"Eigenvector: {eigenvectors[:, dominant_idx]}")

For production code, always use optimized libraries. NumPy’s eig() uses LAPACK routines that handle numerical stability, complex eigenvalues, and edge cases properly.

Practical Applications

Principal Component Analysis (PCA)

PCA uses eigendecomposition of the covariance matrix to find principal components—directions of maximum variance in data. This is the foundation of dimensionality reduction:

from sklearn.datasets import load_iris
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt

# Load dataset
iris = load_iris()
X = iris.data
y = iris.target

# Standardize features
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# Compute covariance matrix
cov_matrix = np.cov(X_scaled.T)

# Eigendecomposition
eigenvalues, eigenvectors = np.linalg.eig(cov_matrix)

# Sort by eigenvalue magnitude
idx = eigenvalues.argsort()[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]

print("Eigenvalues (variance explained):")
print(eigenvalues)
print(f"\nVariance explained ratio: {eigenvalues / eigenvalues.sum()}")

# Project onto first 2 principal components
X_pca = X_scaled @ eigenvectors[:, :2]

# Visualize
plt.figure(figsize=(10, 4))

plt.subplot(1, 2, 1)
plt.scatter(X_scaled[:, 0], X_scaled[:, 1], c=y, cmap='viridis', alpha=0.6)
plt.xlabel('Feature 1')
plt.ylabel('Feature 2')
plt.title('Original Features')

plt.subplot(1, 2, 2)
plt.scatter(X_pca[:, 0], X_pca[:, 1], c=y, cmap='viridis', alpha=0.6)
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.title('Principal Components')
plt.tight_layout()
plt.savefig('pca_comparison.png', dpi=150, bbox_inches='tight')
print("\nVisualization saved as 'pca_comparison.png'")

The first two principal components capture most of the variance, allowing us to reduce from 4 dimensions to 2 while preserving data structure.

PageRank Algorithm

Google’s PageRank treats the web as a giant matrix where eigenvalues determine page importance. The dominant eigenvector of the link matrix gives page rankings.

def pagerank(adjacency_matrix, damping=0.85, iterations=100):
    """Simple PageRank implementation using power iteration."""
    n = adjacency_matrix.shape[0]
    
    # Normalize columns (outgoing links)
    col_sums = adjacency_matrix.sum(axis=0)
    col_sums[col_sums == 0] = 1  # Avoid division by zero
    M = adjacency_matrix / col_sums
    
    # Add damping factor
    M = damping * M + (1 - damping) / n * np.ones((n, n))
    
    # Power iteration
    v = np.ones(n) / n
    for _ in range(iterations):
        v = M @ v
        v = v / v.sum()
    
    return v

# Example web graph (4 pages)
# Page connections: 0->1, 0->2, 1->2, 2->0, 2->3, 3->2
adjacency = np.array([[0, 1, 1, 0],
                      [0, 0, 1, 0],
                      [1, 0, 0, 1],
                      [0, 0, 1, 0]], dtype=float)

ranks = pagerank(adjacency)
print("PageRank scores:")
for i, rank in enumerate(ranks):
    print(f"Page {i}: {rank:.4f}")

Visualization and Interpretation

Visualizing how matrices transform space helps build intuition about eigenvectors:

def visualize_transformation(A):
    """Visualize matrix transformation and eigenvectors."""
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
    
    # Create unit circle
    theta = np.linspace(0, 2*np.pi, 100)
    circle = np.array([np.cos(theta), np.sin(theta)])
    
    # Transform circle
    ellipse = A @ circle
    
    # Get eigenvectors
    eigenvalues, eigenvectors = np.linalg.eig(A)
    
    # Plot original space
    ax1.plot(circle[0], circle[1], 'b-', alpha=0.3, label='Unit circle')
    ax1.arrow(0, 0, eigenvectors[0, 0], eigenvectors[1, 0], 
              head_width=0.1, color='red', label=f'v₁ (λ={eigenvalues[0]:.2f})')
    ax1.arrow(0, 0, eigenvectors[0, 1], eigenvectors[1, 1],
              head_width=0.1, color='green', label=f'v₂ (λ={eigenvalues[1]:.2f})')
    ax1.set_aspect('equal')
    ax1.grid(True, alpha=0.3)
    ax1.legend()
    ax1.set_title('Original Space')
    
    # Plot transformed space
    ax2.plot(ellipse[0], ellipse[1], 'b-', alpha=0.3, label='Transformed')
    ax2.arrow(0, 0, eigenvalues[0]*eigenvectors[0, 0], 
              eigenvalues[0]*eigenvectors[1, 0],
              head_width=0.1, color='red', label=f'λ₁v₁')
    ax2.arrow(0, 0, eigenvalues[1]*eigenvectors[0, 1],
              eigenvalues[1]*eigenvectors[1, 1],
              head_width=0.1, color='green', label=f'λ₂v₂')
    ax2.set_aspect('equal')
    ax2.grid(True, alpha=0.3)
    ax2.legend()
    ax2.set_title('Transformed Space')
    
    plt.tight_layout()
    plt.savefig('transformation_viz.png', dpi=150, bbox_inches='tight')

# Example: rotation + scaling
A = np.array([[2, 0.5],
              [0.5, 1]])
visualize_transformation(A)
print("Transformation visualization saved")

Common Pitfalls and Best Practices

Numerical Stability: Ill-conditioned matrices can produce inaccurate eigenvalues. Always check condition numbers:

# Ill-conditioned matrix example
A_bad = np.array([[1, 1],
                  [1, 1.0001]])

print(f"Condition number: {np.linalg.cond(A_bad):.2e}")

# Small perturbations cause large eigenvalue changes
eigenvalues_1 = np.linalg.eigvals(A_bad)
A_perturbed = A_bad + np.random.randn(2, 2) * 1e-6
eigenvalues_2 = np.linalg.eigvals(A_perturbed)

print(f"Original eigenvalues: {eigenvalues_1}")
print(f"Perturbed eigenvalues: {eigenvalues_2}")
print(f"Difference: {np.abs(eigenvalues_1 - eigenvalues_2)}")

Non-Diagonalizable Matrices: Not all matrices have a complete set of eigenvectors. Defective matrices require generalized eigenvectors or Jordan normal form.

Computational Complexity: Full eigendecomposition is O(n³). For large sparse matrices, use iterative methods targeting specific eigenvalues (e.g., scipy.sparse.linalg.eigs).

Conclusion

Eigenvalues and eigenvectors reveal the fundamental structure of linear transformations. They’re not just mathematical curiosities—they’re computational workhorses powering machine learning, network analysis, and control systems.

The key is understanding the geometric intuition: eigenvectors are invariant directions under transformation, and eigenvalues measure the scaling along those directions. Once you grasp this, applications from PCA to PageRank become natural extensions of the same principle.

For deeper study, explore singular value decomposition (SVD), which generalizes eigendecomposition to non-square matrices, and spectral graph theory, which uses eigenvalues to analyze network structure. The linear algebra underlying these concepts is identical—you’re already equipped with the foundation.