ch192 — Rank, Nullity, and the Fundamental Theorem of Linear Algebra

# --- Visualization: Rank and null space via SVD ---
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-whitegrid')

rng = np.random.default_rng(42)

# Build matrices with known rank
def low_rank_matrix(m, n, rank, noise=0.0, rng=rng):
    """Construct an m x n matrix of given rank with optional noise."""
    U = rng.normal(0, 1, (m, rank))
    V = rng.normal(0, 1, (n, rank))
    A = U @ V.T
    if noise > 0:
        A += noise * rng.normal(0, 1, (m, n))
    return A

# Several matrices with different ranks
configs = [
    (6, 8, 1, 'Rank 1 (6x8)'),
    (6, 8, 3, 'Rank 3 (6x8)'),
    (6, 8, 6, 'Rank 6 — full row rank (6x8)'),
    (6, 8, 6, 'Rank 6 + noise (numerical rank)'),
]

fig, axes = plt.subplots(1, 4, figsize=(14, 4))
for ax, (m, n, r, title) in zip(axes, configs):
    noise = 0.1 if 'noise' in title else 0.0
    A = low_rank_matrix(m, n, r, noise=noise)
    _, s, _ = np.linalg.svd(A, full_matrices=False)
    ax.bar(range(1, len(s)+1), s, color='steelblue', alpha=0.8)
    ax.axhline(1e-10, color='crimson', ls='--', lw=1, label='ε threshold')
    ax.set_title(title, fontsize=8)
    ax.set_xlabel('Singular value index')
    ax.set_ylabel('Singular value')
    ax.legend(fontsize=6)

plt.suptitle('Singular Value Spectra for Matrices of Different Ranks', fontsize=10)
plt.tight_layout()
plt.show()
print("Rank = number of singular values above threshold.")
print("With noise, 'numerical rank' depends on threshold choice.")

# --- Implementation: Four Fundamental Subspaces via SVD ---

def four_fundamental_subspaces(A, tol=1e-10):
    """
    Compute the four fundamental subspaces of A via SVD.
    
    Args:
        A:   (m, n) matrix
        tol: threshold below which singular values are treated as zero
    
    Returns:
        col_space:       (m, r) basis for column space
        null_space:      (n, n-r) basis for null space
        row_space:       (n, r) basis for row space
        left_null_space: (m, m-r) basis for left null space
        rank:            r
    """
    U, s, Vt = np.linalg.svd(A, full_matrices=True)
    
    r = int(np.sum(s > tol))   # numerical rank
    
    col_space       = U[:, :r]        # first r columns of U
    left_null_space = U[:, r:]        # remaining columns of U
    row_space       = Vt[:r, :].T     # first r rows of Vt transposed
    null_space      = Vt[r:, :].T     # remaining rows of Vt transposed
    
    return col_space, null_space, row_space, left_null_space, r

def numerical_rank(A, tol=1e-10):
    """Numerical rank of A: number of singular values above tol."""
    _, s, _ = np.linalg.svd(A, full_matrices=False)
    return int(np.sum(s > tol))

def nullity(A, tol=1e-10):
    """Nullity = n - rank."""
    return A.shape[1] - numerical_rank(A, tol)

# --- Validate on known matrices ---
# Test 1: rank-2 matrix (outer product of two vectors)
u = np.array([1, 2, 3, 4], dtype=float)
v = np.array([1, 1, 0, 0, 1], dtype=float)
A1 = np.outer(u, v)   # rank 1

# Test 2: general 4x5 rank-3 matrix
A2 = low_rank_matrix(4, 5, 3)

# Test 3: 5x3 full column rank
A3 = rng.normal(0, 1, (5, 3))

for A, name in [(A1, 'Outer product (expected rank=1)'),
                (A2, '4x5 rank-3 matrix'),
                (A3, '5x3 full col rank')]:
    cs, ns, rs, lns, r = four_fundamental_subspaces(A)
    m, n = A.shape
    print(f"\n{name}:")
    print(f"  Shape: {m}x{n},  Rank={r},  Nullity={n-r}")
    print(f"  Rank + Nullity = {r} + {n-r} = {r + (n-r)} (should equal n={n})")
    print(f"  col_space dim: {cs.shape[1]}, null_space dim: {ns.shape[1] if ns.size > 0 else 0}")
    
    # Orthogonality check: col_space ⊥ left_null_space
    if lns.size > 0 and cs.size > 0:
        orth_err = np.abs(cs.T @ lns).max()
        print(f"  Orthogonality (col ⊥ left_null): max error = {orth_err:.2e}")

# --- Experiment 1: Solvability of Ax = b ---
# Hypothesis: Ax=b is solvable iff b is in the column space of A.
# Try changing: RANK, the RHS vector b

RANK = 2   # <-- try 1, 2, 3
m, n = 4, 5
A = low_rank_matrix(m, n, RANK)
cs, ns, rs, lns, r = four_fundamental_subspaces(A)

# b in column space (solvable)
b_solvable = cs @ rng.normal(0, 1, r)      # linear combo of col space basis

# b not in column space (no solution)
if lns.size > 0:
    b_impossible = lns[:, 0]   # lives in left null space — perpendicular to col space
else:
    b_impossible = rng.normal(0, 1, m)   # random (likely not in col space)

# Check: project b onto col space; if residual is zero, b is in col space
def in_column_space(b, col_space, tol=1e-9):
    """Check if b lies in the span of col_space columns."""
    proj = col_space @ col_space.T @ b
    residual = np.linalg.norm(b - proj)
    return residual < tol, residual

in_cs, res_s = in_column_space(b_solvable, cs)
in_ci, res_i = in_column_space(b_impossible, cs)

print(f"Rank-{RANK} matrix ({m}x{n}):")
print(f"  b_solvable  in col space: {in_cs}  (residual={res_s:.2e})")
print(f"  b_impossible in col space: {in_ci}  (residual={res_i:.2e})")

# --- Experiment 2: Null Space and Solution Families ---
# Hypothesis: if nullity > 0, there are infinitely many solutions to Ax = b.
# Try changing: RANK_EXP

RANK_EXP = 2   # <-- try 1, 2, 3, 4
m_e, n_e = 4, 5
A_e = low_rank_matrix(m_e, n_e, RANK_EXP)
cs_e, ns_e, _, _, r_e = four_fundamental_subspaces(A_e)

b_e = cs_e @ rng.normal(0, 1, r_e)   # solvable b

# Find a particular solution via least-norm pseudoinverse
# Pseudoinverse via SVD: A+ = V S+ U.T
U_e, s_e, Vt_e = np.linalg.svd(A_e, full_matrices=True)
s_inv = np.where(s_e > 1e-10, 1/s_e, 0)
# Build S+ (m x n)
S_plus = np.zeros((n_e, m_e))
for i in range(min(m_e, n_e)):
    if i < len(s_inv):
        S_plus[i, i] = s_inv[i]
A_pinv = Vt_e.T @ S_plus @ U_e.T
x_particular = A_pinv @ b_e

# Verify: Ax_particular ≈ b
print(f"Ax_particular - b residual: {np.linalg.norm(A_e @ x_particular - b_e):.2e}")

# General solution family: x_particular + alpha * null_vec for each null vector
if ns_e.size > 0:
    ALPHA_VALS = np.linspace(-2, 2, 20)
    null_vec = ns_e[:, 0]
    residuals = [np.linalg.norm(A_e @ (x_particular + a * null_vec) - b_e)
                 for a in ALPHA_VALS]
    
    fig, ax = plt.subplots(figsize=(8, 3))
    ax.plot(ALPHA_VALS, residuals, 'o-', ms=4, color='steelblue')
    ax.set_xlabel('Alpha (weight of null space component)')
    ax.set_ylabel('||Ax - b|| (should be 0)')
    ax.set_title(f'All solutions Ax=b: x_particular + alpha * null_vec (nullity={n_e - r_e})')
    ax.set_ylim(-0.001, max(residuals) * 1.2 + 0.001)
    plt.tight_layout()
    plt.show()
    print("Residual is ~0 for all alpha: infinitely many exact solutions.")

# --- Mini Project: Network Constraint Analysis ---
# Problem: A water distribution network has nodes (junctions) and edges (pipes).
#          The incidence matrix B[e, v] = +1 if edge e leaves node v, -1 if it enters.
#          Flow conservation: B.T @ flow = demand. The rank of B.T determines
#          whether the system is solvable and how many degrees of freedom exist.
# Task: Build a network incidence matrix, analyze its four subspaces,
#       and interpret each subspace in terms of the network.

# Network: 5 nodes, 6 edges
#   0-1, 1-2, 2-3, 3-0, 1-3, 0-2
nodes = 5
edges = [(0,1), (1,2), (2,3), (3,0), (1,3), (0,2), (2,4), (3,4)]
n_edges = len(edges)

# Build incidence matrix B (n_edges x nodes)
B = np.zeros((n_edges, nodes))
for e_idx, (i, j) in enumerate(edges):
    B[e_idx, i] = +1
    B[e_idx, j] = -1

print(f"Incidence matrix B: {B.shape}  ({n_edges} edges x {nodes} nodes)")
print(f"B =\n{B.astype(int)}")

cs, ns, rs, lns, r = four_fundamental_subspaces(B)
print(f"\nRank(B):   {r}")
print(f"Nullity(B): {n_edges - r}  (= number of independent cycles in the network)")
print(f"Left Null Space dim: {lns.shape[1]}  (= number of connected components - 1)")
print(f"\nNull space of B (cycle basis):")
print(f"  Each null space vector represents a circulation (flow around a cycle).")
if ns.size > 0:
    print(ns.round(3))
print(f"\nLeft null space of B:")
print(f"  The vector [1,1,1,...,1]/sqrt(n) is always here (sum of all rows = 0).")
if lns.size > 0:
    print(lns.round(3))