Eigenvalues & Eigenvectors

import numpy as np

print("=" * 55)
print("EIGENVALUES & EIGENVECTORS")
print("=" * 55)

# Simple 2×2 example
A = np.array([[4, 2], [1, 3]])
eigenvalues, eigenvectors = np.linalg.eig(A)

print(f"\nA =\n{A}")
print(f"\nEigenvalues: λ₁ = {eigenvalues[0]:.4f}, λ₂ = {eigenvalues[1]:.4f}")
print(f"Eigenvectors:")
for i in range(2):
    v = eigenvectors[:, i]
    print(f"  v{i+1} = {np.round(v, 4)}")

# Verify: Av = λv
for i in range(2):
    v = eigenvectors[:, i]
    lam = eigenvalues[i]
    av = A @ v
    lv = lam * v
    print(f"\n  Av{i+1} = {np.round(av, 4)}")
    print(f"  λ{i+1}v{i+1} = {np.round(lv, 4)}")
    print(f"  Equal? {np.allclose(av, lv)} ✓")

# Trace and Determinant
print(f"\nTrace(A) = {np.trace(A)}, Sum(λ) = {sum(eigenvalues):.4f}")
print(f"Det(A) = {np.linalg.det(A):.4f}, Prod(λ) = {np.prod(eigenvalues):.4f}")

# Symmetric matrix: real eigenvalues, orthogonal eigenvectors
print("\n" + "=" * 55)
print("SYMMETRIC MATRIX: SPECTRAL THEOREM")
print("=" * 55)

S = np.array([[5, 2, 1],
              [2, 6, 3],
              [1, 3, 7]])

eigenvalues_s, eigenvectors_s = np.linalg.eigh(S)  # Use eigh for symmetric
print(f"\nSymmetric S =\n{S}")
print(f"Eigenvalues (sorted): {np.round(eigenvalues_s, 4)}")
print(f"All real? {np.all(np.isreal(eigenvalues_s))} ✓")

# Orthogonality check
Q = eigenvectors_s
print(f"\nQᵀQ (should be I):")
print(np.round(Q.T @ Q, 4))

# Reconstruction: A = QΛQᵀ
Lambda = np.diag(eigenvalues_s)
S_reconstructed = Q @ Lambda @ Q.T
print(f"\nReconstruction error: {np.linalg.norm(S - S_reconstructed):.2e} ✓")

# PCA Connection
print("\n" + "=" * 55)
print("📊 PCA = EIGENDECOMPOSITION OF COVARIANCE")
print("=" * 55)

np.random.seed(42)
# 2D correlated data
x1 = np.random.normal(0, 3, 200)
x2 = 0.7 * x1 + np.random.normal(0, 1, 200)
data = np.column_stack([x1, x2])

# Center data
data_centered = data - data.mean(axis=0)
cov_matrix = np.cov(data_centered, rowvar=False)

eigenvalues_pca, eigenvectors_pca = np.linalg.eigh(cov_matrix)
# Sort descending
idx = np.argsort(eigenvalues_pca)[::-1]
eigenvalues_pca = eigenvalues_pca[idx]
eigenvectors_pca = eigenvectors_pca[:, idx]

print(f"\nCovariance matrix:\n{np.round(cov_matrix, 3)}")
print(f"\nPrincipal Components:")
for i in range(2):
    var_explained = eigenvalues_pca[i] / sum(eigenvalues_pca) * 100
    print(f"  PC{i+1}: direction={np.round(eigenvectors_pca[:,i], 3)}, variance={eigenvalues_pca[i]:.3f} ({var_explained:.1f}%)")
print(f"\n💡 PC1 captures {eigenvalues_pca[0]/sum(eigenvalues_pca)*100:.1f}% of variance!")