Numerical Linear Algebra

import numpy as np

# Hilbert matrix: H_ij = 1/(i+j-1)
# Theoretical solution: x = [1, 1, ..., 1]
n = 12
H = np.array([[1/(i+j-1) for j in range(n)] for i in range(n)])
b = H @ np.ones(n)

# "Solve" Ax = b
x_numpy = np.linalg.solve(H, b)
print(f"Error (np.linalg.solve): {np.linalg.norm(x_numpy - np.ones(n)):.2e}")
# ~1e-4 (terrible for a 12x12 system!)

# The condition number tells you WHY
print(f"Condition number: {np.linalg.cond(H):.2e}")
# ~10^16 — the matrix is nearly singular in floating point

import numpy as np
import matplotlib.pyplot as plt

# Condition number growth for Hilbert matrices
ns = range(2, 15)
kappas = []
for n in ns:
    H = np.array([[1/(i+j+1) for j in range(n)] for i in range(n)])
    kappas.append(np.linalg.cond(H))

for n, k in zip(ns, kappas):
    print(f"H_{n}: κ = {k:.2e}")

# For n >= 12, κ exceeds 10^16 — beyond double precision capability

import numpy as np

def gaussian_elimination_no_pivot(A, b):
    """Unstable Gaussian elimination (no pivoting)."""
    n = len(b)
    Ab = np.hstack([A.astype(float), b.reshape(-1, 1)])
    for col in range(n):
        for row in range(col + 1, n):
            factor = Ab[row, col] / Ab[col, col]
            Ab[row, col:] -= factor * Ab[col, col:]
    # Back substitution
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (Ab[i, -1] - Ab[i, i+1:n] @ x[i+1:]) / Ab[i, i]
    return x

# Near-singular example
A = np.array([[1e-10, 1], [1, 1]], dtype=float)
b = A @ np.array([1.0, 1.0])

# numpy uses stable LAPACK (LU with partial pivoting)
x_stable = np.linalg.solve(A, b)
print(f"Stable: {x_stable} (error: {np.linalg.norm(x_stable - [1,1]):.2e})")

# No-pivot version may fail
try:
    x_unstable = gaussian_elimination_no_pivot(A, b)
    print(f"Unstable: {x_unstable} (error: {np.linalg.norm(x_unstable - [1,1]):.2e})")
except ZeroDivisionError:
    print("Unstable: Division by zero!")

import numpy as np

# Classic catastrophic cancellation
a = 1.0
b = 1.0 + 1e-15
print(f"Direct subtraction: {(b - a) - 1e-15}")  # May not be exactly 0

# Catastrophic cancellation in quadratic formula
# For b >> sqrt(b^2 - 4ac), the standard formula loses precision
a_coeff, b_coeff, c_coeff = 1, 1e8, 1
sqrt_disc = np.sqrt(b_coeff**2 - 4*a_coeff*c_coeff)

x1_standard = (-b_coeff + sqrt_disc) / (2 * a_coeff)
x2_standard = (-b_coeff - sqrt_disc) / (2 * a_coeff)
print(f"Standard: x1 = {x1_standard:.6e}, x2 = {x2_standard:.6e}")

# Stable alternative: use the other root to avoid cancellation
x2_stable = c_coeff / (a_coeff * x1_standard)  # = c/(a*x1)
x1_stable = (-b_coeff - np.sign(b_coeff) * sqrt_disc) / (2 * a_coeff)
print(f"Stable:   x1 = {x1_stable:.6e}, x2 = {x2_stable:.6e}")

def kahan_sum(values):
    """Kahan compensated summation — nearly exact for large sums."""
    s = 0.0
    c = 0.0  # compensation for lost low-order bits
    for x in values:
        y = x - c
        t = s + y
        c = (t - s) - y  # algebraically zero, but captures rounding error
        s = t
    return s

# Demonstration: sum 1e7 values of 1.0
n = 10_000_000
values = [1.0] * n
print(f"Naive sum:   {sum(values) - n}")       # May show rounding error
print(f"Kahan sum:   {kahan_sum(values) - n}")  # Exact (or nearly so)

# More dramatic: sum decreasing values
vals = [1e16, 1.0, -1e16]
print(f"Naive:   {sum(vals)}")          # 0.0 (correct by luck)
vals2 = [1e16, 1.0, -1e16, 1.0, -1e16]
print(f"Naive:   {sum(vals2)}")         # 0.0 (lost the 1.0!)
print(f"Kahan:   {kahan_sum(vals2)}")   # 1.0 (correct)

import numpy as np
from scipy.linalg import lu, solve_triangular

A = np.array([[2, 1, 1], [4, 3, 3], [8, 7, 9]], dtype=float)
b = np.array([1, 2, 3], dtype=float)

# LU decomposition via scipy
P, L, U = lu(A)

# Verify: PA = LU
print("PA = LU check:", np.allclose(P @ A, L @ U))

# Solve using LU factors
y = solve_triangular(L, P @ b, lower=True)
x = solve_triangular(U, y)
print(f"Solution: {x}")
print(f"Verify:   {A @ x}")  # Should equal b

# numpy's solve uses LAPACK's LU implementation internally
x_numpy = np.linalg.solve(A, b)
print(f"numpy:    {x_numpy}")

import numpy as np
from scipy.linalg import cho_factor, cho_solve

# Create SPD matrix
A = np.array([[4, 2, 1], [2, 5, 3], [1, 3, 6]], dtype=float)
b = np.array([1, 2, 3], dtype=float)

# Cholesky decomposition
L = np.linalg.cholesky(A)
print("A = LL^T check:", np.allclose(L @ L.T, A))

# Solve using Cholesky
c, low = cho_factor(A)
x = cho_solve((c, low), b)
print(f"Solution: {x}")
print(f"Verify:   {A @ x}")

import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import cg, gmres

# Create a large SPD tridiagonal system
n = 10000
A = diags([-1, 2, -1], [-1, 0, 1], shape=(n, n), format='csr')
b = np.ones(n)

# Conjugate Gradient
x_cg, info = cg(A, b, tol=1e-10, maxiter=1000)
print(f"CG: {np.linalg.norm(A @ x_cg - b):.2e} residual, info={info}")

# Jacobi iteration (manual)
def jacobi(A, b, x0=None, tol=1e-10, maxiter=1000):
    x = x0 if x0 is not None else np.zeros_like(b)
    D = A.diagonal()
    for k in range(maxiter):
        r = b - A @ x
        if np.linalg.norm(r) < tol:
            return x, k
        x = x + r / D  # Jacobi step: x_new = D^{-1}(b - (A-D)x)
    return x, maxiter

x_jac, iters = jacobi(A, b)
print(f"Jacobi: {np.linalg.norm(A @ x_jac - b):.2e} residual, {iters} iterations")

import numpy as np

# Overdetermined system: 5 equations, 3 unknowns
np.random.seed(42)
A = np.random.randn(5, 3)
b = np.random.randn(5)

# Method 1: Normal equations (UNSTABLE)
x_normal = np.linalg.solve(A.T @ A, A.T @ b)
print(f"Normal equations: residual = {np.linalg.norm(A @ x_normal - b):.2e}")
print(f"  κ(A^T A) = {np.linalg.cond(A.T @ A):.2e}")

# Method 2: QR decomposition (STABLE)
x_qr, res_qr, rank_qr, sv_qr = np.linalg.lstsq(A, b, rcond=None)
print(f"QR (lstsq):       residual = {np.linalg.norm(A @ x_qr - b):.2e}")

# Method 3: SVD (MOST STABLE)
U, s, Vt = np.linalg.svd(A, full_matrices=False)
x_svd = Vt.T @ np.diag(1/s) @ U.T @ b
print(f"SVD:              residual = {np.linalg.norm(A @ x_svd - b):.2e}")

# Rank-deficient example
A_rankdef = np.array([[1, 2], [2, 4], [3, 6]], dtype=float)  # rank 1
b_rankdef = np.array([1, 2, 3.1], dtype=float)
x_rd = np.linalg.lstsq(A_rankdef, b_rankdef, rcond=None)[0]
print(f"\nRank-deficient least squares: x = {x_rd}")
print(f"  Residual: {np.linalg.norm(A_rankdef @ x_rd - b_rankdef):.2e}")

import numpy as np
from scipy.sparse import diags, random
from scipy.sparse.linalg import cg, spilu, LinearOperator

# Create a difficult SPD system (convection-diffusion)
n = 500
A = diags([-1, 4, -1], [-1, 0, 1], shape=(n, n), format='csr')
A = A + 0.1 * random(n, n, density=0.01, format='csr')  # Add some non-symmetric perturbation
A = A @ A.T  # Make SPD
b = np.random.randn(n)

# Without preconditioning
x_no_precon, info1 = cg(A, b, tol=1e-10, maxiter=1000)
print(f"No preconditioner: {np.linalg.norm(A @ x_no_precon - b):.2e}")

# With Jacobi (diagonal) preconditioner
M_inv = 1.0 / A.diagonal()
M = LinearOperator(A.shape, matvec=lambda x: M_inv * x)
x_jac, info2 = cg(A, b, tol=1e-10, maxiter=1000, M=M)
print(f"Jacobi precon:     {np.linalg.norm(A @ x_jac - b):.2e}")

# With incomplete Cholesky
ilu = spilu(A.tocsc(), drop_tol=1e-4)
M_ilu = LinearOperator(A.shape, matvec=ilu.solve)
x_ilu, info3 = cg(A, b, tol=1e-10, maxiter=1000, M=M_ilu)
print(f"ILU precon:        {np.linalg.norm(A @ x_ilu - b):.2e}")

import numpy as np

# === Core Operations ===

# Solve Ax = b
A = np.array([[3, 1], [1, 2]], dtype=float)
b = np.array([9, 8], dtype=float)
x = np.linalg.solve(A, b)
print(f"Solve: x = {x}")

# Determinant and inverse (AVOID inverting in production!)
print(f"det(A) = {np.linalg.det(A):.4f}")
print(f"inv(A) = \n{np.linalg.inv(A)}")

# Matrix factorizations
Q, R = np.linalg.qr(A)
print(f"\nQR: Q = \n{Q}\nR = \n{R}")

L = np.linalg.cholesky(A @ A.T + np.eye(2))  # Ensure SPD
print(f"Cholesky: L = \n{L}")

U, s, Vt = np.linalg.svd(A)
print(f"\nSVD: singular values = {s}")

# Eigenvalues
eigvals, eigvecs = np.linalg.eig(A)
print(f"Eigenvalues: {eigvals}")

# Matrix norms
print(f"\n||A||_F = {np.linalg.norm(A, 'fro'):.4f}")
print(f"||A||_2 = {np.linalg.norm(A, 2):.4f}")
print(f"||A||_1 = {np.linalg.norm(A, 1):.4f}")
print(f"||A||_inf = {np.linalg.norm(A, np.inf):.4f}")

# Condition number
print(f"κ(A) = {np.linalg.cond(A):.4f}")
print(f"κ(A^T A) = {np.linalg.cond(A.T @ A):.4f}")  # κ^2

from scipy import sparse
from scipy.sparse import linalg as spla
import numpy as np

# Create sparse matrix (5-point Laplacian)
n = 100
A = sparse.diags([-1, 4, -1], [-1, 0, 1], shape=(n, n), format='csr')
b = np.random.randn(n)

# Solve sparse system
x, info = spla.cg(A, b, tol=1e-10)
print(f"CG solve: residual = {np.linalg.norm(A @ x - b):.2e}")

# Sparse eigenvalues (largest/smallest)
vals, vecs = spla.eigs(A, k=5, which='LM')
print(f"Largest eigenvalues: {np.sort(vals.real)}")

vals_sm, _ = spla.eigs(A, k=3, which='SM')
print(f"Smallest eigenvalues: {np.sort(vals_sm.real)}")

# Sparse SVD (useful for large-scale PCA)
U, s, Vt = spla.svds(A, k=5)
print(f"Top 5 singular values: {s}")

import numpy as np

# Simulating mixed-precision issues
def forward_pass_fp32(W, x):
    """Full precision forward pass."""
    return W @ x

def forward_pass_fp16(W, x):
    """Simulated mixed-precision forward pass."""
    W_16 = W.astype(np.float16)
    x_16 = x.astype(np.float16)
    result = W_16 @ x_16
    return result.astype(np.float32)

np.random.seed(42)
W = np.random.randn(100, 100) * 0.1
x = np.random.randn(100)

y_fp32 = forward_pass_fp32(W, x)
y_fp16 = forward_pass_fp16(W, x)

print(f"FP32 result norm: {np.linalg.norm(y_fp32):.6f}")
print(f"FP16 result norm: {np.linalg.norm(y_fp16):.6f}")
print(f"Relative error:   {np.linalg.norm(y_fp32 - y_fp16) / np.linalg.norm(y_fp32):.2e}")

# Condition number of a typical layer's weight matrix
W_layer = np.random.randn(512, 512) * (2 / 512)**0.5  # Xavier init
print(f"\nLayer condition number: {np.linalg.cond(W_layer):.2f}")
# Usually ~10-100 with good initialization, much worse after training

import numpy as np

# BAD: Using inv(A) @ b
A = np.array([[1, 2], [3, 4]], dtype=float)
b = np.array([5, 6], dtype=float)

# Wrong
x_bad = np.linalg.inv(A) @ b  # Slower, less stable

# Correct
x_good = np.linalg.solve(A, b)  # Direct, stable

print(f"inv @ b: {x_bad}")
print(f"solve:  {x_good}")
print(f"Equal:  {np.allclose(x_bad, x_good)}")

# BAD: Comparing floats with ==
a = 0.1 + 0.2
print(f"\n0.1 + 0.2 == 0.3? {a == 0.3}")          # False!
print(f"0.1 + 0.2 ≈ 0.3?  {np.isclose(a, 0.3)}")  # True