Use this Skill for SVD, PCA, eigendecomposition, Cholesky, iterative solvers, sparse matrix formats, and condition number analysis with numpy and scipy.
Use this Skill for SVD, PCA, eigendecomposition, Cholesky, iterative solvers, sparse matrix formats, and condition number analysis with numpy and scipy.
One-line summary: Apply SVD, PCA, eigendecomposition, Cholesky factorization, iterative solvers (CG, GMRES), and sparse matrix techniques for scientific computing.
When to Use This Skill
When decomposing matrices for dimensionality reduction or latent factor models
When solving large linear systems $Ax = b$ efficiently (direct or iterative)
When analyzing matrix condition numbers to diagnose ill-conditioning
When working with sparse matrices from finite element or graph problems
When implementing PCA or truncated SVD for data analysis
When factorizing positive definite matrices (Cholesky for fast solves)
Trigger keywords: SVD, PCA, eigendecomposition, Cholesky, sparse matrix, GMRES, conjugate gradient, condition number, matrix factorization, linear system
Background & Key Concepts
Singular Value Decomposition (SVD)
Any $m \times n$ matrix $A$ decomposes as:
$$
A = U \Sigma V^T
$$
where $U \in \mathbb{R}^{m \times m}$, $\Sigma \in \mathbb{R}^{m \times n}$ (diagonal, non-negative), $V \in \mathbb{R}^{n \times n}$ are orthogonal. SVD underpins PCA, pseudoinverse, and low-rank approximations.
Eigendecomposition
For a square matrix $A$: $Av = \lambda v$ where $\lambda$ is an eigenvalue and $v$ the eigenvector. For symmetric $A = Q \Lambda Q^T$ (spectral theorem).
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_digits
# Load example data
digits = load_digits()
X = digits.data.astype(float) # shape (1797, 64)
X -= X.mean(axis=0) # center# Covariance matrix approach (small p)
C = X.T @ X / (len(X) - 1) # 64×64 covariance
eigenvalues, eigenvectors = np.linalg.eigh(C)
# eigh returns in ascending order; reverse for descending
eigenvalues = eigenvalues[::-1]
eigenvectors = eigenvectors[:, ::-1]
# Explained variance ratio
total_var = eigenvalues.sum()
exp_var_ratio = eigenvalues / total_var
cumulative_var = np.cumsum(exp_var_ratio)
print(f"n_components for 90% variance: {np.argmax(cumulative_var >= 0.90) + 1}")
print(f"Top 5 eigenvalues: {eigenvalues[:5].round(2)}")
# Project to first 2 PCs
X_pca = X @ eigenvectors[:, :2]
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
axes[0].plot(cumulative_var[:30], 'o-')
axes[0].axhline(0.90, color='r', linestyle='--', label='90% threshold')
axes[0].set_xlabel("Number of components")
axes[0].set_ylabel("Cumulative explained variance")
axes[0].legend(); axes[0].set_title("PCA Explained Variance")
scatter = axes[1].scatter(X_pca[:, 0], X_pca[:, 1],
c=digits.target, cmap='tab10', s=5, alpha=0.7)
plt.colorbar(scatter, ax=axes[1], label='Digit')
axes[1].set_title("PCA Projection (first 2 PCs)")
plt.tight_layout()
plt.savefig("pca_digits.png", dpi=150)
plt.show()
Step 3: Sparse Matrices and Iterative Solvers
import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as spla
import matplotlib.pyplot as plt
defbuild_1d_laplacian(n):
"""Build tridiagonal Laplacian (finite difference) as sparse CSR matrix."""
diags = [-np.ones(n-1), 2*np.ones(n), -np.ones(n-1)]
return sp.diags(diags, [-1, 0, 1], shape=(n, n), format='csr')
n = 1000
A = build_1d_laplacian(n)
b = np.ones(n)
print(f"Matrix: {A.shape}, nnz={A.nnz}, density={A.nnz/(n**2):.6f}")
print(f"Condition number estimate: {spla.norm(A) * spla.norm(spla.inv(A.tocsc())):.2e}")
# Conjugate Gradient solver (works for symmetric positive definite)
residuals_cg = []
defcg_callback(xk):
r = b - A @ xk
residuals_cg.append(np.linalg.norm(r))
x_cg, info_cg = spla.cg(A, b, tol=1e-10, maxiter=5000, callback=cg_callback)
print(f"\nCG converged: {info_cg == 0}, iterations: {len(residuals_cg)}")
print(f"Residual: {np.linalg.norm(b - A @ x_cg):.2e}")
# GMRES solver (general non-symmetric)
residuals_gmres = []
defgmres_callback(rk):
residuals_gmres.append(rk)
x_gmres, info_gmres = spla.gmres(A, b, tol=1e-10, maxiter=500, callback=gmres_callback)
print(f"GMRES converged: {info_gmres == 0}, iterations: {len(residuals_gmres)}")
# ILU preconditioner for GMRES
ilu = spla.spilu(A.tocsc(), fill_factor=2.0)
M = spla.LinearOperator(A.shape, ilu.solve)
x_prec, info_prec = spla.gmres(A, b, M=M, tol=1e-10, maxiter=100)
print(f"Preconditioned GMRES converged: {info_prec == 0}")
fig, ax = plt.subplots(figsize=(8, 5))
ax.semilogy(residuals_cg, label="CG")
if residuals_gmres:
ax.semilogy(residuals_gmres, label="GMRES")
ax.set_xlabel("Iteration"); ax.set_ylabel("Residual norm")
ax.legend(); ax.set_title("Solver Convergence")
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("solver_convergence.png", dpi=150)
plt.show()
Advanced Usage
Cholesky Factorization for Fast Linear Solves
import numpy as np
import scipy.linalg as la
import time
defgenerate_spd(n, seed=42):
"""Generate a random symmetric positive definite matrix."""
rng = np.random.default_rng(seed)
A = rng.standard_normal((n, n))
return A @ A.T + n * np.eye(n) # add n*I to ensure positive definiteness
n = 500
A = generate_spd(n)
b = np.ones(n)
# Cholesky factorization
t0 = time.time()
L, lower = la.cho_factor(A)
x = la.cho_solve((L, lower), b)
print(f"Cholesky solve: {time.time()-t0:.4f}s, residual={np.linalg.norm(A@x-b):.2e}")
# LU (for comparison)
t0 = time.time()
x_lu = la.solve(A, b)
print(f"LU solve: {time.time()-t0:.4f}s, residual={np.linalg.norm(A@x_lu-b):.2e}")
# Multiple RHS — Cholesky amortizes cost
B = np.random.randn(n, 50)
t0 = time.time()
X = la.cho_solve((L, lower), B)
print(f"Cholesky (50 RHS): {time.time()-t0:.4f}s")
Condition Number and Preconditioning
import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as spla
defanalyze_conditioning(A):
"""Report condition number and recommend preconditioning strategy."""if sp.issparse(A):
# Estimate via power iteration (full SVD too expensive)
s_max = spla.norm(A)
# Smallest singular value via inverse iterationtry:
s_min = 1.0 / spla.norm(spla.inv(A.tocsc()))
kappa = s_max / s_min
except Exception:
kappa = float('inf')
else:
s = np.linalg.svd(A, compute_uv=False)
kappa = s[0] / s[-1]
print(f"Condition number κ(A) = {kappa:.2e}")
if kappa < 1e4:
print(" Well-conditioned: direct solver sufficient")
elif kappa < 1e10:
print(" Moderately ill-conditioned: consider Jacobi/ILU preconditioning")
else:
print(" Severely ill-conditioned: use regularization (Tikhonov, truncated SVD)")
return kappa
# Test with progressively ill-conditioned matricesfor cond_target in [1e2, 1e6, 1e10]:
# Build matrix with prescribed condition number
n = 50
U, _ = np.linalg.qr(np.random.randn(n, n))
s = np.logspace(0, -np.log10(cond_target), n)
A = U @ np.diag(s) @ U.T
kappa = analyze_conditioning(A)
Troubleshooting
Error: numpy.linalg.LinAlgError: Matrix is singular
Cause: Matrix is rank-deficient or numerically singular.
Fix:
# Use pseudoinverse for rank-deficient systems
x = np.linalg.lstsq(A, b, rcond=None)[0]
# Or add regularization (Tikhonov)
lambda_reg = 1e-6
x = np.linalg.solve(A.T @ A + lambda_reg * np.eye(A.shape[1]), A.T @ b)
Issue: CG doesn't converge
Cause: Matrix is not symmetric positive definite.
Fix:
# Check symmetryprint(f"Max asymmetry: {np.abs(A - A.T).max():.2e}") # should be ~0# Check positive definiteness
eigenvalues = np.linalg.eigvalsh(A)
print(f"Min eigenvalue: {eigenvalues.min():.2e}") # must be > 0
