How to Calculate Eigenvalues and Eigenvectors in Python

Key Insights

NumPy’s eig() and eigh() functions handle eigenvalue computation efficiently, with eigh() offering better performance and numerical stability for symmetric matrices—always prefer it when applicable.
For large sparse matrices, SciPy’s scipy.sparse.linalg.eigs() can compute only the top-k eigenvalues, reducing computation time from O(n³) to near-linear complexity in many cases.
Eigenvectors form the mathematical foundation of PCA by identifying directions of maximum variance in data—understanding this connection transforms eigendecomposition from abstract math into a practical dimensionality reduction tool.

Introduction to Eigenvalues and Eigenvectors

Eigenvalues and eigenvectors reveal fundamental properties of linear transformations. When you multiply a matrix A by its eigenvector v, the result is simply a scaled version of that same vector: Av = λv, where λ (lambda) is the eigenvalue representing the scaling factor.

This seemingly simple relationship has profound practical applications. In Principal Component Analysis (PCA), eigenvectors identify directions of maximum variance in high-dimensional data. In image compression, eigendecomposition enables algorithms like Singular Value Decomposition (SVD) to represent images with fewer parameters. Stability analysis in dynamical systems relies on eigenvalues to determine whether systems converge or diverge over time.

The key insight: eigenvectors represent directions that remain unchanged under transformation (only stretched or compressed), while eigenvalues quantify that scaling.

Manual Calculation with NumPy

NumPy provides numpy.linalg.eig() for general matrices and numpy.linalg.eigh() for Hermitian (symmetric real) matrices. Always use eigh() when possible—it’s faster and more numerically stable.

import numpy as np

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

eigenvalues, eigenvectors = np.linalg.eig(A)

print("Eigenvalues:", eigenvalues)
print("Eigenvectors:\n", eigenvectors)

Output shows two eigenvalues (approximately 5 and 2) and corresponding eigenvectors as columns in the eigenvector matrix.

Let’s verify the fundamental equation Av = λv:

# Verify for the first eigenvector
v1 = eigenvectors[:, 0]
lambda1 = eigenvalues[0]

left_side = A @ v1
right_side = lambda1 * v1

print("Av =", left_side)
print("λv =", right_side)
print("Match:", np.allclose(left_side, right_side))

For symmetric matrices, use eigh():

# Symmetric matrix
S = np.array([[6, 2, 1],
              [2, 5, 2],
              [1, 2, 4]])

eigenvalues_sym, eigenvectors_sym = np.linalg.eigh(S)

print("Eigenvalues (sorted):", eigenvalues_sym)

Note that eigh() returns eigenvalues in ascending order, while eig() doesn’t guarantee ordering. You can reconstruct the original matrix using eigendecomposition:

# Reconstruct: A = Q * Λ * Q^T (for symmetric matrices)
Q = eigenvectors_sym
Lambda = np.diag(eigenvalues_sym)

S_reconstructed = Q @ Lambda @ Q.T
print("Reconstruction error:", np.max(np.abs(S - S_reconstructed)))

The reconstruction error should be near machine epsilon (around 1e-15), confirming numerical accuracy.

Using SciPy for Advanced Cases

SciPy extends NumPy’s capabilities for specialized scenarios. For sparse matrices—common in graph analysis and large-scale machine learning—computing all eigenvalues is wasteful and slow.

from scipy.sparse import csr_matrix
from scipy.sparse.linalg import eigs
import time

# Create a large sparse matrix
n = 1000
density = 0.01
np.random.seed(42)

# Generate sparse symmetric matrix
A_dense = np.random.randn(n, n)
A_dense = (A_dense + A_dense.T) / 2  # Make symmetric
mask = np.random.rand(n, n) > density
A_dense[mask] = 0

A_sparse = csr_matrix(A_dense)

print(f"Matrix size: {n}x{n}, Non-zeros: {A_sparse.nnz}")

# Compute only top 5 eigenvalues
start = time.time()
eigenvalues_sparse, eigenvectors_sparse = eigs(A_sparse, k=5, which='LM')
sparse_time = time.time() - start

print(f"Top 5 eigenvalues: {eigenvalues_sparse.real}")
print(f"Computation time: {sparse_time:.4f}s")

The which='LM' parameter specifies “Largest Magnitude” eigenvalues. Other options include ‘SM’ (smallest), ‘LR’ (largest real), and ‘SR’ (smallest real).

For comparison, computing all eigenvalues on the dense version:

# Dense computation (slower)
start = time.time()
eigenvalues_dense, _ = np.linalg.eigh(A_dense)
dense_time = time.time() - start

print(f"Dense computation time: {dense_time:.4f}s")
print(f"Speedup: {dense_time / sparse_time:.1f}x")

For a 1000×1000 matrix with 1% density, sparse methods typically achieve 10-50x speedup when computing only top eigenvalues.

Practical Application: Principal Component Analysis (PCA)

PCA uses eigendecomposition to find orthogonal directions of maximum variance. Let’s implement it from scratch and compare with scikit-learn.

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

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

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# Manual PCA implementation
def pca_manual(X, n_components=2):
    # Compute covariance matrix
    cov_matrix = np.cov(X.T)
    
    # Compute eigenvalues and eigenvectors
    eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix)
    
    # Sort in descending order
    idx = eigenvalues.argsort()[::-1]
    eigenvalues = eigenvalues[idx]
    eigenvectors = eigenvectors[:, idx]
    
    # Select top components
    components = eigenvectors[:, :n_components]
    
    # Project data
    X_transformed = X @ components
    
    return X_transformed, eigenvalues, components

X_pca_manual, eigenvals, components = pca_manual(X_scaled, n_components=2)

# Compare with sklearn
pca_sklearn = SklearnPCA(n_components=2)
X_pca_sklearn = pca_sklearn.fit_transform(X_scaled)

print("Manual eigenvalues:", eigenvals[:2])
print("Sklearn explained variance:", pca_sklearn.explained_variance_)
print("Results match:", np.allclose(np.abs(X_pca_manual), np.abs(X_pca_sklearn)))

The variance explained by each principal component equals its corresponding eigenvalue:

# Explained variance ratio
explained_variance_ratio = eigenvals / eigenvals.sum()
print("\nExplained variance ratio:")
for i, ratio in enumerate(explained_variance_ratio):
    print(f"PC{i+1}: {ratio:.3f} ({ratio*100:.1f}%)")

For the Iris dataset, the first two principal components typically explain 95%+ of total variance, justifying dimensionality reduction from 4D to 2D.

Visualization and Interpretation

Visualizing eigenvectors helps build intuition about their geometric meaning:

# Plot PCA results
plt.figure(figsize=(12, 5))

# Subplot 1: Transformed data
plt.subplot(1, 2, 1)
for label in np.unique(y):
    mask = y == label
    plt.scatter(X_pca_manual[mask, 0], X_pca_manual[mask, 1], 
                label=iris.target_names[label], alpha=0.7)
plt.xlabel('First Principal Component')
plt.ylabel('Second Principal Component')
plt.title('PCA Projection')
plt.legend()
plt.grid(True, alpha=0.3)

# Subplot 2: Scree plot
plt.subplot(1, 2, 2)
plt.bar(range(1, len(eigenvals) + 1), eigenvals)
plt.xlabel('Principal Component')
plt.ylabel('Eigenvalue (Variance)')
plt.title('Scree Plot')
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('pca_visualization.png', dpi=150, bbox_inches='tight')

To visualize how eigenvectors behave under transformation:

# 2D transformation visualization
A_2d = np.array([[2, 1],
                 [1, 2]])

eigenvalues_2d, eigenvectors_2d = np.linalg.eig(A_2d)

# Create grid of points
theta = np.linspace(0, 2*np.pi, 50)
circle = np.array([np.cos(theta), np.sin(theta)])

# Transform circle
ellipse = A_2d @ circle

plt.figure(figsize=(10, 5))

# Original space
plt.subplot(1, 2, 1)
plt.plot(circle[0], circle[1], 'b-', label='Original')
plt.quiver(0, 0, eigenvectors_2d[0, 0], eigenvectors_2d[1, 0], 
           angles='xy', scale_units='xy', scale=1, color='r', width=0.006)
plt.quiver(0, 0, eigenvectors_2d[0, 1], eigenvectors_2d[1, 1], 
           angles='xy', scale_units='xy', scale=1, color='g', width=0.006)
plt.axis('equal')
plt.grid(True, alpha=0.3)
plt.title('Before Transformation')
plt.legend()

# Transformed space
plt.subplot(1, 2, 2)
plt.plot(ellipse[0], ellipse[1], 'b-', label='Transformed')
plt.quiver(0, 0, eigenvalues_2d[0]*eigenvectors_2d[0, 0], 
           eigenvalues_2d[0]*eigenvectors_2d[1, 0],
           angles='xy', scale_units='xy', scale=1, color='r', width=0.006)
plt.quiver(0, 0, eigenvalues_2d[1]*eigenvectors_2d[0, 1], 
           eigenvalues_2d[1]*eigenvectors_2d[1, 1],
           angles='xy', scale_units='xy', scale=1, color='g', width=0.006)
plt.axis('equal')
plt.grid(True, alpha=0.3)
plt.title('After Transformation')
plt.legend()

plt.tight_layout()

The eigenvectors (red and green arrows) maintain their direction but scale by their eigenvalues.

Common Pitfalls and Best Practices

Numerical Stability: Near-singular matrices produce unreliable eigenvalues. Always check the condition number:

# Check matrix conditioning
A_ill = np.array([[1, 1],
                  [1, 1.0000001]])

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

if cond > 1e10:
    print("Warning: Matrix is ill-conditioned")
    
eigenvalues_ill, _ = np.linalg.eig(A_ill)
print("Eigenvalues:", eigenvalues_ill)

Complex Eigenvalues: Non-symmetric matrices can have complex eigenvalues. Handle them explicitly:

# Rotation matrix (complex eigenvalues)
theta = np.pi / 4
R = np.array([[np.cos(theta), -np.sin(theta)],
              [np.sin(theta), np.cos(theta)]])

eigenvalues_complex, _ = np.linalg.eig(R)
print("Complex eigenvalues:", eigenvalues_complex)
print("Magnitudes:", np.abs(eigenvalues_complex))

Eigenvector Normalization: NumPy returns unit eigenvectors, but verify this when implementing custom algorithms:

# Check normalization
for i in range(eigenvectors.shape[1]):
    norm = np.linalg.norm(eigenvectors[:, i])
    print(f"Eigenvector {i} norm: {norm:.6f}")

Performance Guidelines: Use eigh() for symmetric matrices, sparse methods for large matrices with low density, and compute only required eigenvalues when possible. For matrices larger than 10,000×10,000, consider iterative methods or randomized algorithms.

Eigenvalue computation is fundamental to modern data science. Master these techniques, understand when to apply each method, and you’ll have powerful tools for dimensionality reduction, feature extraction, and system analysis.