Chapter 140 — Span

# --- Visualization: Span of 1, 2, and 3 vectors in R^3 ---
# Shows how span grows from a line -> plane -> all of R^3

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

plt.style.use('seaborn-v0_8-whitegrid')
fig = plt.figure(figsize=(15, 5))

v1 = np.array([1.0, 0.5, 0.0])
v2 = np.array([0.0, 1.0, 0.5])
v3 = np.array([0.3, 0.0, 1.0])

# --- Panel 1: span({v1}) — a line ---
ax1 = fig.add_subplot(131, projection='3d')
t = np.linspace(-2, 2, 100)
line = np.outer(t, v1)
ax1.plot(line[:, 0], line[:, 1], line[:, 2], 'b-', lw=2)
ax1.quiver(0, 0, 0, *v1, color='red', lw=2.5, arrow_length_ratio=0.2)
ax1.set_title('span({v₁})\n= a line (1D)')
ax1.set_xlabel('x'); ax1.set_ylabel('y'); ax1.set_zlabel('z')

# --- Panel 2: span({v1, v2}) — a plane ---
ax2 = fig.add_subplot(132, projection='3d')
s_vals = np.linspace(-1.5, 1.5, 20)
t_vals = np.linspace(-1.5, 1.5, 20)
S, T = np.meshgrid(s_vals, t_vals)
X = S * v1[0] + T * v2[0]
Y = S * v1[1] + T * v2[1]
Z = S * v1[2] + T * v2[2]
ax2.plot_surface(X, Y, Z, alpha=0.25, color='blue')
ax2.quiver(0,0,0, *v1, color='red', lw=2.5, arrow_length_ratio=0.2)
ax2.quiver(0,0,0, *v2, color='green', lw=2.5, arrow_length_ratio=0.2)
ax2.set_title('span({v₁, v₂})\n= a plane (2D)')
ax2.set_xlabel('x'); ax2.set_ylabel('y'); ax2.set_zlabel('z')

# --- Panel 3: span({v1, v2, v3}) — all of R^3 (if independent) ---
ax3 = fig.add_subplot(133, projection='3d')
# Show that we can reach a point off the v1-v2 plane
B = np.column_stack([v1, v2, v3])
target = np.array([0.5, 0.5, 0.8])
c = np.linalg.solve(B, target)

ax3.plot_surface(X, Y, Z, alpha=0.1, color='blue')  # v1-v2 plane for reference
ax3.quiver(0,0,0,*v1,color='red',lw=2,arrow_length_ratio=0.2)
ax3.quiver(0,0,0,*v2,color='green',lw=2,arrow_length_ratio=0.2)
ax3.quiver(0,0,0,*v3,color='purple',lw=2,arrow_length_ratio=0.2)
ax3.quiver(0,0,0,*target,color='orange',lw=2.5,arrow_length_ratio=0.2)
ax3.text(*target+0.05, f'target\n(reachable!)', fontsize=8, color='orange')
ax3.set_title(f'span({{v₁,v₂,v₃}}) = ℝ³ (3D)\ncoords={c.round(2)}')
ax3.set_xlabel('x'); ax3.set_ylabel('y'); ax3.set_zlabel('z')

plt.suptitle('How span grows with each independent vector added', fontsize=12, y=1.01)
plt.tight_layout()
plt.show()

# --- Implementation: Span membership test and dimension ---

import numpy as np


def in_span(b, vectors, tol=1e-9):
    """
    Test whether b is in the span of the given vectors.

    Solves: A @ c = b (least squares) and checks the residual.

    Args:
        b: np.ndarray, shape (n,)
        vectors: list of np.ndarray, each shape (n,)
        tol: float — residual tolerance

    Returns:
        dict: {'in_span': bool, 'coefficients': np.ndarray or None, 'residual': float}
    """
    A = np.column_stack(vectors)
    c, _, _, _ = np.linalg.lstsq(A, b, rcond=None)
    residual = float(np.linalg.norm(A @ c - b))
    return {
        'in_span': residual < tol,
        'coefficients': c if residual < tol else None,
        'residual': residual
    }


def span_dimension(vectors, tol=1e-10):
    """
    Compute the dimension of span(vectors).

    Args:
        vectors: list of np.ndarray

    Returns:
        int — rank of the matrix formed by columns
    """
    A = np.column_stack(vectors)
    return np.linalg.matrix_rank(A, tol=tol)


def minimal_spanning_set(vectors, tol=1e-10):
    """
    Reduce a set of vectors to a minimal spanning set (a basis for their span).

    Uses column pivoting via QR decomposition to identify independent columns.

    Args:
        vectors: list of np.ndarray

    Returns:
        list of np.ndarray — a subset of the input that spans the same space
    """
    A = np.column_stack(vectors)
    # QR with column pivoting: first k columns (after pivoting) are independent
    Q, R, P = np.linalg.qr(A, mode='complete'), None, None
    # Use SVD to identify independent columns more robustly
    U, s, Vt = np.linalg.svd(A, full_matrices=False)
    rank = np.sum(s > tol)
    # Return U columns as the minimal spanning basis (orthonormalized)
    return [U[:, i] for i in range(rank)]


# --- Tests ---
print("=== Span membership tests ===")

v1 = np.array([1.0, 0.0, 0.0])
v2 = np.array([0.0, 1.0, 0.0])

# b in the xy-plane: should be in span
b1 = np.array([3.0, -2.0, 0.0])
# b off the xy-plane: should NOT be in span
b2 = np.array([1.0, 1.0, 1.0])

r1 = in_span(b1, [v1, v2])
r2 = in_span(b2, [v1, v2])

print(f"b1={b1} in span({{e1,e2}}): {r1['in_span']} (coeffs: {r1['coefficients']})")
print(f"b2={b2} in span({{e1,e2}}): {r2['in_span']} (residual: {r2['residual']:.4f})")

print()
print("=== Span dimension ===")
# Three vectors in R^3 — but two are parallel
redundant = [
    np.array([1.0, 2.0, 3.0]),
    np.array([2.0, 4.0, 6.0]),  # = 2 * first
    np.array([0.0, 1.0, 0.0])
]
print(f"3 vectors (one redundant): span dim = {span_dimension(redundant)}  (expected 2)")

independent = [
    np.array([1.0, 0.0, 0.0]),
    np.array([0.0, 1.0, 0.0]),
    np.array([0.0, 0.0, 1.0])
]
print(f"3 independent vectors: span dim = {span_dimension(independent)}  (expected 3)")

# --- Visualization: how span dimension changes as you add vectors ---

import numpy as np
import matplotlib.pyplot as plt

plt.style.use('seaborn-v0_8-whitegrid')

rng = np.random.default_rng(7)

AMBIENT_DIM = 5   # <-- try 3, 8, 10
MAX_VECTORS = 10  # number of vectors to add one by one

vectors = []
span_dims = []

for i in range(MAX_VECTORS):
    # Alternate between truly new vectors and dependent ones
    if i < AMBIENT_DIM:
        # Independent: random vector unlikely to be in current span
        v = rng.standard_normal(AMBIENT_DIM)
    else:
        # Dependent: linear combination of existing vectors
        coeffs = rng.standard_normal(len(vectors))
        v = sum(c * vec for c, vec in zip(coeffs, vectors))
    vectors.append(v)
    if len(vectors) == 1:
        dim = 1
    else:
        dim = span_dimension(vectors)
    span_dims.append(dim)

fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(range(1, MAX_VECTORS+1), span_dims, 'bo-', markersize=8, linewidth=2)
ax.axhline(AMBIENT_DIM, color='red', linestyle='--', label=f'Max dim = {AMBIENT_DIM}')
ax.fill_between(range(1, MAX_VECTORS+1), span_dims,
                 alpha=0.15, color='blue')
ax.set_xlabel('Number of vectors added')
ax.set_ylabel('dim(span)')
ax.set_title(f'Span dimension grows until saturation at dim={AMBIENT_DIM}\n'
              f'Adding dependent vectors (after index {AMBIENT_DIM}) changes nothing')
ax.set_xticks(range(1, MAX_VECTORS+1))
ax.legend()
plt.tight_layout()
plt.show()

# --- Experiment 1: When does a set of vectors span R^n? ---
# Hypothesis: k vectors can span R^n only if k >= n AND they are independent.
# Try changing: n and k to observe the threshold.

import numpy as np

rng = np.random.default_rng(42)

n = 4       # dimension of ambient space  <-- modify
k_values = range(1, 8)  # number of vectors to try  <-- modify

print(f"Testing random vectors in R^{n}:")
print(f"{'k vectors':>12}  {'dim(span)':>10}  {'spans R^n?':>12}")

for k in k_values:
    A = rng.standard_normal((n, k))
    d = np.linalg.matrix_rank(A)
    spans_all = (d == n)
    print(f"  k={k:2d}          {d:10d}  {str(spans_all):>12}")

print(f"\nConclusion: need at least k={n} linearly independent vectors to span R^{n}.")
print(f"With k > {n}, extra vectors are redundant (rank still capped at {n}).")

# --- Experiment 2: Span as the column space of a matrix ---
# Hypothesis: b in span({v1,...,vk}) iff Ax=b has a solution.
# Try changing: the target b and observe when solutions exist.

import numpy as np

# Matrix A: columns are the spanning vectors
A = np.array([
    [1, 0, 1],
    [0, 1, 1],
    [0, 0, 0]   # third row zero: rank < 3
], dtype=float)

print(f"Matrix A (columns = spanning vectors):")
print(A)
print(f"Rank = {np.linalg.matrix_rank(A)} (spans a 2D subspace of R^3)")
print()

# b vectors to test
targets = [
    (np.array([2.0, 3.0, 0.0]), "In column space (z=0)"),
    (np.array([1.0, 1.0, 1.0]), "NOT in column space (z=1)"),
    (np.array([0.5, -0.5, 0.0]), "In column space (z=0)"),  # <-- modify
]

for b, description in targets:
    c, residuals, rank, _ = np.linalg.lstsq(A, b, rcond=None)
    residual = np.linalg.norm(A @ c - b)
    in_cs = residual < 1e-9
    print(f"b = {b}  ({description})")
    print(f"  Residual: {residual:.2e}  ->  in span: {in_cs}")
    if in_cs:
        print(f"  Coefficients: {c.round(4)}")
    print()

# --- Experiment 3: Removing vectors from a spanning set ---
# Hypothesis: A spanning set can be reduced to a basis without losing span.
# Approach: greedily remove vectors while maintaining the same span dimension.

import numpy as np

rng = np.random.default_rng(5)

N_DIM = 3  # ambient space  <-- modify
N_VECTORS = 7  # start with an overcomplete set  <-- modify

# Generate an overcomplete spanning set
base_vectors = [rng.standard_normal(N_DIM) for _ in range(N_DIM)]
extra_vectors = [sum(rng.random(N_DIM) * v for v in base_vectors)
                 for _ in range(N_VECTORS - N_DIM)]
all_vectors = base_vectors + extra_vectors

target_dim = np.linalg.matrix_rank(np.column_stack(all_vectors))
print(f"Starting with {len(all_vectors)} vectors, span dim = {target_dim}")

# Greedy reduction: try removing each vector; keep it removed if span dim unchanged
current_set = list(all_vectors)
removed = []

i = 0
while i < len(current_set):
    trial = current_set[:i] + current_set[i+1:]
    if len(trial) == 0:
        break
    d = np.linalg.matrix_rank(np.column_stack(trial))
    if d == target_dim:
        removed.append(i)
        current_set = trial  # remove it
    else:
        i += 1

final_dim = np.linalg.matrix_rank(np.column_stack(current_set))
print(f"After removal: {len(current_set)} vectors remain, span dim = {final_dim}")
print(f"This is a basis for the span: {len(current_set) == target_dim}")

# --- Mini Project: Feature Reachability Analyzer ---
#
# Problem: In ML, a model can only represent targets that lie in the span
# of its features (for linear models). Given a dataset of features and targets,
# determine which targets are representable (in the column space of the
# feature matrix) and which are not.
#
# This is the fundamental question of linear regression: is y reachable?

import numpy as np
import matplotlib.pyplot as plt

plt.style.use('seaborn-v0_8-whitegrid')
rng = np.random.default_rng(0)

# --- Dataset ---
# Features: x1 (income), x2 (age) — 2 features, N samples
N = 50
x1 = rng.uniform(20, 100, N)   # income (thousands)
x2 = rng.uniform(20, 65, N)    # age
X = np.column_stack([x1, x2])  # feature matrix: shape (N, 2)

# Three target vectors:
# y1: perfectly representable (linear combo of features)
y1 = 0.5*x1 - 0.3*x2 + 0.0   # exact linear combo — in span
# y2: representable with noise
y2 = 0.5*x1 - 0.3*x2 + rng.normal(0, 2, N)  # noisy — approximately in span
# y3: not representable (depends on x1^2)
y3 = 0.01 * x1**2 + rng.normal(0, 1, N)  # nonlinear — NOT in span of linear features

def fit_and_residual(X, y):
    """
    Fit y to span of X (least squares) and return residual norm.
    Large residual = y not in span(X).
    """
    coeffs, _, _, _ = np.linalg.lstsq(X, y, rcond=None)
    y_hat = X @ coeffs
    residual = np.linalg.norm(y - y_hat)
    r2 = 1 - np.var(y - y_hat) / np.var(y)
    return coeffs, y_hat, residual, r2

# --- TODO: run analysis for each target ---
targets = [('y1 (exact linear)', y1), ('y2 (linear + noise)', y2), ('y3 (nonlinear)', y3)]

fig, axes = plt.subplots(1, 3, figsize=(13, 4))

for ax, (name, y) in zip(axes, targets):
    coeffs, y_hat, residual, r2 = fit_and_residual(X, y)
    ax.scatter(y, y_hat, alpha=0.5, s=20)
    mn, mx = min(y.min(), y_hat.min()), max(y.max(), y_hat.max())
    ax.plot([mn, mx], [mn, mx], 'r--', lw=1.5, label='Perfect fit')
    ax.set_xlabel('Actual y')
    ax.set_ylabel('Predicted ŷ')
    ax.set_title(f'{name}\nR²={r2:.3f}  residual={residual:.2f}')
    ax.legend(fontsize=8)
    print(f"{name}: R²={r2:.4f}, residual norm={residual:.4f}")
    print(f"  Coefficients: {coeffs.round(4)}")
    in_span = r2 > 0.99
    print(f"  In span of X: {in_span}")
    print()

plt.suptitle('Feature Reachability: Is the target in span(X)?', fontsize=12)
plt.tight_layout()
plt.show()