Chapter 138 — Subspaces

# --- Visualization: Subspaces of R^2 and R^3 ---
# Shows which geometric objects are subspaces (pass through origin, flat)
# and which are not (shifted, curved)

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

plt.style.use('seaborn-v0_8-whitegrid')
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# --- Plot 1: R^2 subspaces ---
ax = axes[0]
t = np.linspace(-2, 2, 100)

# Subspace: line through origin (y = 2x)
ax.plot(t, 2*t, 'b-', linewidth=2, label='Line through origin (subspace)')

# Non-subspace: line NOT through origin (y = 2x + 1)
ax.plot(t, 2*t + 1, 'r--', linewidth=2, label='Shifted line (NOT subspace)')

# Non-subspace: circle
theta = np.linspace(0, 2*np.pi, 100)
ax.plot(np.cos(theta), np.sin(theta), 'g:', linewidth=2, label='Circle (NOT subspace)')

ax.plot(0, 0, 'ko', markersize=8, label='Origin')
ax.set_xlim(-2.5, 2.5); ax.set_ylim(-2.5, 2.5)
ax.set_aspect('equal')
ax.set_title('Subspaces of ℝ²', fontsize=13)
ax.set_xlabel('x'); ax.set_ylabel('y')
ax.legend(fontsize=8)

# --- Plot 2: Closure test visualization ---
ax2 = axes[1]
ax2.set_xlim(-3, 3); ax2.set_ylim(-3, 3)
ax2.set_aspect('equal')
ax2.axhline(0, color='gray', linewidth=0.5)
ax2.axvline(0, color='gray', linewidth=0.5)

# Subspace: the y-axis (all vectors [0, y])
ax2.plot([0, 0], [-3, 3], 'b-', linewidth=3, label='y-axis (subspace)')

# Demonstrate closure: two vectors in y-axis
v1 = np.array([0, 1.5])
v2 = np.array([0, -0.5])
v_sum = v1 + v2
ax2.annotate('', xy=v1, xytext=[0,0], arrowprops=dict(arrowstyle='->', color='blue', lw=2))
ax2.annotate('', xy=v2, xytext=[0,0], arrowprops=dict(arrowstyle='->', color='cyan', lw=2))
ax2.annotate('', xy=v_sum, xytext=[0,0], arrowprops=dict(arrowstyle='->', color='navy', lw=2, linestyle='dashed'))
ax2.text(0.1, 1.5, 'v₁', fontsize=11); ax2.text(0.1, -0.5, 'v₂', fontsize=11)
ax2.text(0.1, v_sum[1], 'v₁+v₂ ∈ subspace', fontsize=9)
ax2.set_title('Closure under addition', fontsize=13)
ax2.set_xlabel('x'); ax2.set_ylabel('y')
ax2.legend(fontsize=9)

# --- Plot 3: Failure of closure ---
ax3 = axes[2]
ax3.set_xlim(-3, 3); ax3.set_ylim(-3, 3)
ax3.set_aspect('equal')

# Non-subspace: first quadrant (not closed under negation)
x_fill = np.array([0, 3, 3, 0])
y_fill = np.array([0, 0, 3, 3])
ax3.fill(x_fill, y_fill, alpha=0.2, color='red', label='First quadrant (NOT subspace)')
ax3.plot([0, 3], [0, 0], 'r-', lw=2)
ax3.plot([0, 0], [0, 3], 'r-', lw=2)

v = np.array([1.5, 1.0])
neg_v = -v
ax3.annotate('', xy=v, xytext=[0,0], arrowprops=dict(arrowstyle='->', color='blue', lw=2))
ax3.annotate('', xy=neg_v, xytext=[0,0], arrowprops=dict(arrowstyle='->', color='red', lw=2))
ax3.text(v[0]+0.1, v[1], 'v ∈ set', fontsize=10)
ax3.text(neg_v[0]+0.1, neg_v[1], '-v ∉ set!', fontsize=10, color='red')
ax3.plot(0, 0, 'ko', markersize=8)
ax3.set_title('Closure failure (negation)', fontsize=13)
ax3.set_xlabel('x'); ax3.set_ylabel('y')
ax3.legend(fontsize=9)

plt.suptitle('Subspaces: What Qualifies and What Doesn\'t', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

# --- Implementation: Subspace verification and fundamental subspaces ---
import numpy as np

def is_closed_under_addition(vectors, test_pairs=50, atol=1e-10):
    """
    Test closure under addition for a set of vectors.
    
    Strategy: sample random linear combinations and check membership
    (membership approximated by checking if the projection residual is near zero).
    
    Args:
        vectors: array of shape (n_vecs, dim) spanning the candidate subspace
        test_pairs: number of random pairs to test
        atol: tolerance for membership check
    
    Returns:
        bool: True if closure appears to hold
    """
    n = len(vectors)
    rng = np.random.default_rng(42)
    
    for _ in range(test_pairs):
        # Pick two random vectors from the span
        c1 = rng.standard_normal(n)
        c2 = rng.standard_normal(n)
        v1 = vectors.T @ c1
        v2 = vectors.T @ c2
        v_sum = v1 + v2
        
        # Check if v_sum is in the span of vectors
        # Project v_sum onto span(vectors) and check residual
        V = vectors.T  # columns are basis vectors
        # Solve least squares: V @ x ≈ v_sum
        x, residuals, rank, sv = np.linalg.lstsq(V, v_sum, rcond=None)
        projection = V @ x
        residual = np.linalg.norm(v_sum - projection)
        
        if residual > atol:
            return False
    return True


def column_space_basis(A, tol=1e-10):
    """
    Compute an orthonormal basis for the column space of A.
    Uses QR decomposition with column pivoting.
    
    Args:
        A: matrix of shape (m, n)
        tol: threshold for detecting zero singular values
    
    Returns:
        basis: array of shape (m, r) where r = rank(A)
    """
    # SVD gives the most numerically stable basis
    U, s, Vt = np.linalg.svd(A, full_matrices=False)
    rank = np.sum(s > tol)
    return U[:, :rank]  # columns of U span the column space


def null_space_basis(A, tol=1e-10):
    """
    Compute an orthonormal basis for the null space of A.
    The null space: all x such that Ax = 0.
    
    Args:
        A: matrix of shape (m, n)
        tol: threshold
    
    Returns:
        basis: array of shape (n, (n - rank)) or empty if trivial
    """
    U, s, Vt = np.linalg.svd(A, full_matrices=True)
    rank = np.sum(s > tol)
    # The last (n - rank) rows of Vt span the null space
    null_basis = Vt[rank:].T  # shape (n, n-rank)
    return null_basis


def verify_null_space(A, null_basis, atol=1e-10):
    """
    Verify: every column of null_basis satisfies Ax = 0.
    """
    for i in range(null_basis.shape[1]):
        x = null_basis[:, i]
        residual = np.linalg.norm(A @ x)
        if residual > atol:
            return False, i, residual
    return True, -1, 0.0


# --- Example: Find subspaces of a 3x4 matrix ---
np.random.seed(7)
A = np.array([
    [1, 2, 0, 1],
    [2, 4, 1, 3],
    [0, 0, 1, 1]
], dtype=float)

print(f"Matrix A ({A.shape[0]}x{A.shape[1]}):")
print(A)
print()

col_basis = column_space_basis(A)
nul_basis = null_space_basis(A)
rank = col_basis.shape[1]
nullity = nul_basis.shape[1] if nul_basis.size > 0 else 0

print(f"Rank (dimension of column space): {rank}")
print(f"Nullity (dimension of null space): {nullity}")
print(f"Rank + Nullity = {rank + nullity} = n = {A.shape[1]}  ← Rank-Nullity Theorem")
print()
print("Column space basis (orthonormal):")
print(np.round(col_basis, 4))
print()
print("Null space basis:")
print(np.round(nul_basis, 4))

ok, idx, res = verify_null_space(A, nul_basis)
print(f"\nNull space verification: A @ null_basis ≈ 0? {ok}")

# --- Experiment 1: Which sets are subspaces? ---
# Hypothesis: A set defined by a homogeneous linear condition is a subspace;
# a set defined by an inhomogeneous condition is not.

import numpy as np

def sample_from_hyperplane(a, b, n=200):
    """
    Sample vectors from {x in R^2 : a^T x = b}.
    When b=0 this is a subspace; when b≠0 it is an affine subspace (not a subspace).
    """
    # Parametrize: for x1 free, x2 = (b - a[0]*x1) / a[1]
    x1 = np.linspace(-3, 3, n)
    x2 = (b - a[0] * x1) / a[1]
    return np.column_stack([x1, x2])

a = np.array([1.0, 2.0])
B_values = [0, 1, -2]  # <-- modify: try different values of b

fig, axes = plt.subplots(1, 3, figsize=(14, 4))
for ax, b in zip(axes, B_values):
    pts = sample_from_hyperplane(a, b)
    ax.plot(pts[:, 0], pts[:, 1], 'b-', lw=2)
    ax.plot(0, 0, 'ko', markersize=8)
    contains_origin = abs(a @ np.array([0, 0]) - b) < 1e-10
    is_sub = "SUBSPACE" if contains_origin else "NOT a subspace"
    ax.set_title(f'x + 2y = {b}\n{is_sub}', fontsize=11)
    ax.set_xlim(-3, 3); ax.set_ylim(-4, 4)
    ax.set_xlabel('x'); ax.set_ylabel('y')
    ax.axhline(0, color='gray', lw=0.5); ax.axvline(0, color='gray', lw=0.5)

plt.suptitle('Homogeneous (b=0) vs Inhomogeneous (b≠0) Linear Constraints', fontsize=12)
plt.tight_layout(); plt.show()
print("Observation: only the line through the origin is a subspace.")

# --- Experiment 2: Rank-Nullity Theorem ---
# Hypothesis: For any m×n matrix, rank + nullity = n always.
# Try changing: matrix dimensions, rank, entries

import numpy as np

def check_rank_nullity(m, n, target_rank=None):
    rng = np.random.default_rng(0)
    if target_rank is None:
        target_rank = min(m, n)
    # Build a rank-deficient matrix: A = U @ V where U is m×r, V is r×n
    U = rng.standard_normal((m, target_rank))
    V = rng.standard_normal((target_rank, n))
    A = U @ V
    
    col_basis = column_space_basis(A)
    nul_basis = null_space_basis(A)
    rank = col_basis.shape[1]
    nullity = nul_basis.shape[1] if nul_basis.size > 0 else 0
    
    print(f"  Matrix {m}×{n}, target_rank={target_rank}: rank={rank}, nullity={nullity}, "
          f"rank+nullity={rank+nullity} (n={n}) ✓" if rank+nullity==n else " MISMATCH")

print("Rank-Nullity checks:")
for m, n, r in [(3, 5, 2), (4, 4, 3), (2, 6, 1), (5, 3, 3)]:  # <-- modify these
    check_rank_nullity(m, n, target_rank=r)

# --- Experiment 3: Intersection of subspaces ---
# Hypothesis: The intersection of two subspaces is always a subspace.
# Two planes through the origin in R^3 intersect in a line (also through origin).

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Plane 1: normal n1 = [1, 0, 0]  → the y-z plane (x=0)
# Plane 2: normal n2 = [0, 1, 0]  → the x-z plane (y=0)
# Intersection: x=0 AND y=0 → the z-axis

n1 = np.array([1.0, 0, 0])  # <-- modify: change the normal vectors
n2 = np.array([0.0, 1, 0])

# Intersection direction = n1 × n2
intersection_dir = np.cross(n1, n2)
if np.linalg.norm(intersection_dir) < 1e-10:
    print("Planes are parallel (or identical) — intersection is a plane or empty")
else:
    intersection_dir = intersection_dir / np.linalg.norm(intersection_dir)
    print(f"Plane 1 normal: {n1}")
    print(f"Plane 2 normal: {n2}")
    print(f"Intersection direction (line through origin): {np.round(intersection_dir, 4)}")
    print("Verify: n1 · dir =", np.dot(n1, intersection_dir).round(10), "(should be 0)")
    print("Verify: n2 · dir =", np.dot(n2, intersection_dir).round(10), "(should be 0)")

# --- Mini Project: Constraint Satisfier via Null Space ---
#
# Problem: In robotics and physics, you often have a system with constraints.
# Given a robot arm with n joints, subject to m linear constraints (Ax = b),
# the space of solutions is an affine subspace: one particular solution + null space of A.
#
# This is also how linear systems work: the general solution to Ax = b
# is x = x_particular + x_null, where x_null ∈ N(A).
#
# Task: Given a constraint matrix A and target b, find the general solution space
# and visualize the solution set.

import numpy as np
import matplotlib.pyplot as plt

# Constraint: 2x + y = 5 (1 equation, 2 unknowns → 1D solution space)
A = np.array([[2.0, 1.0]])
b = np.array([5.0])

# Step 1: Find a particular solution (minimum norm solution via pseudoinverse)
x_particular = np.linalg.lstsq(A, b, rcond=None)[0]
print(f"Particular solution: x_p = {x_particular}")
print(f"Verify Ax_p = b: {A @ x_particular} (should be {b})")

# Step 2: Find null space basis
null_basis = null_space_basis(A)
print(f"\nNull space basis:\n{null_basis}")
print(f"Null space dimension: {null_basis.shape[1]} (rank-nullity: n=2, rank=1)")

# Step 3: General solution = x_particular + t * null_basis[:,0]
t_values = np.linspace(-4, 4, 100)
solutions = x_particular[:, None] + null_basis @ (t_values[None, :])  # shape (2, 100)

# Verify all solutions satisfy the constraint
residuals = A @ solutions - b[:, None]
print(f"\nMax constraint violation over all solutions: {np.max(np.abs(residuals)):.2e}")

# Step 4: Visualize
fig, ax = plt.subplots(figsize=(7, 6))
ax.plot(solutions[0], solutions[1], 'b-', lw=2, label='Solution set (affine subspace)')
ax.plot(*x_particular, 'ro', markersize=10, label='Particular solution (min norm)')
ax.annotate('', xy=x_particular + null_basis[:, 0], xytext=x_particular,
            arrowprops=dict(arrowstyle='->', color='green', lw=2))
ax.text(*(x_particular + null_basis[:, 0] + 0.1), 'null space\ndirection', fontsize=10, color='green')
ax.plot(0, 0, 'ks', markersize=8, label='Origin')
ax.set_xlim(-4, 6); ax.set_ylim(-4, 6)
ax.set_xlabel('x₁'); ax.set_ylabel('x₂')
ax.set_title('General solution to 2x₁ + x₂ = 5\n= particular solution + null space', fontsize=12)
ax.legend(); ax.set_aspect('equal')
plt.tight_layout(); plt.show()

# TODO: Try a 2-equation, 3-unknown system (solution is a line in R^3)
# TODO: Visualize the null space shifting with different b values