Chapter 151 — Introduction to Matrices

# --- Visualization: Matrix as a transformation of the plane ---
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-whitegrid')

def draw_grid(ax, A, title, color='steelblue'):
    """Draw a unit grid before and after matrix transformation A."""
    # Grid lines in original space
    lines = np.linspace(-2, 2, 9)
    for v in lines:
        # Horizontal lines
        pts = np.array([[t, v] for t in np.linspace(-2, 2, 50)]).T
        transformed = A @ pts
        ax.plot(transformed[0], transformed[1], color=color, alpha=0.3, linewidth=0.8)
        # Vertical lines
        pts = np.array([[v, t] for t in np.linspace(-2, 2, 50)]).T
        transformed = A @ pts
        ax.plot(transformed[0], transformed[1], color=color, alpha=0.3, linewidth=0.8)
    # Basis vectors
    for vec, col, label in [(A[:,0], 'crimson', 'e₁'), (A[:,1], 'darkgreen', 'e₂')]:
        ax.annotate('', xy=vec, xytext=(0,0),
                    arrowprops=dict(arrowstyle='->', color=col, lw=2))
        ax.text(vec[0]*1.1, vec[1]*1.1, label, color=col, fontsize=11)
    ax.set_xlim(-3, 3); ax.set_ylim(-3, 3)
    ax.set_aspect('equal')
    ax.axhline(0, color='black', lw=0.5); ax.axvline(0, color='black', lw=0.5)
    ax.set_title(title, fontsize=11)

# Three characteristic matrices
I  = np.eye(2)                          # Identity: does nothing
R  = np.array([[0, -1], [1, 0]])        # 90-degree rotation
S  = np.array([[2, 0], [0, 0.5]])       # Scale x by 2, y by 0.5
Sh = np.array([[1, 1], [0, 1]])         # Shear

fig, axes = plt.subplots(1, 4, figsize=(16, 4))
configs = [
    (I,  'Identity\n[[1,0],[0,1]]'),
    (R,  'Rotation 90°\n[[0,-1],[1,0]]'),
    (S,  'Scaling\n[[2,0],[0,0.5]]'),
    (Sh, 'Shear\n[[1,1],[0,1]]'),
]
for ax, (mat, title) in zip(axes, configs):
    draw_grid(ax, mat, title)

plt.suptitle('Each Matrix = A Different Way of Warping the Plane', fontsize=13, y=1.02)
plt.tight_layout()
plt.show()

# --- Mathematical Formulation: Matrix-vector product two ways ---
import numpy as np

A = np.array([[2, 1],
              [0, 3],
              [1, -1]], dtype=float)  # shape: (3, 2) — maps R² → R³

v = np.array([4.0, 2.0])             # shape: (2,) — input vector in R²

# Method 1: Row interpretation (dot products)
result_rows = np.array([np.dot(A[i], v) for i in range(A.shape[0])])

# Method 2: Column interpretation (linear combination)
result_cols = v[0] * A[:, 0] + v[1] * A[:, 1]

# Method 3: NumPy @ operator
result_numpy = A @ v

print(f"A shape: {A.shape}  (maps R² → R³)")
print(f"v shape: {v.shape}")
print(f"Result via row dot products:  {result_rows}")
print(f"Result via column combo:      {result_cols}")
print(f"Result via A @ v:             {result_numpy}")
print(f"All equal: {np.allclose(result_rows, result_cols) and np.allclose(result_cols, result_numpy)}")

# --- Implementation: Matrix creation, inspection, and basic properties ---
import numpy as np

def matrix_info(A, name='A'):
    """
    Print key properties of a matrix.

    Args:
        A: 2D numpy array
        name: string label for display
    """
    m, n = A.shape
    print(f"--- Matrix {name} ---")
    print(f"  Shape:     {m} x {n}  ({'square' if m==n else 'rectangular'})")
    print(f"  dtype:     {A.dtype}")
    print(f"  Rank:      {np.linalg.matrix_rank(A)}  (max possible: {min(m,n)})")
    print(f"  Nonzeros:  {np.count_nonzero(A)} / {A.size}")
    print(f"  Frobenius norm: {np.linalg.norm(A, 'fro'):.4f}")
    print(f"  Trace:     {np.trace(A) if m==n else 'N/A (not square)'}")
    print()

# Common matrix constructors
A_zeros    = np.zeros((3, 4))           # all zeros
A_ones     = np.ones((2, 5))            # all ones
A_identity = np.eye(4)                  # identity matrix
A_diag     = np.diag([1, 2, 3, 4])     # diagonal matrix
A_random   = np.random.randn(3, 3)     # random Gaussian entries
A_range    = np.arange(12).reshape(3, 4)  # from a flat list

for mat, name in [(A_identity, 'Identity 4x4'),
                  (A_diag, 'Diagonal'),
                  (A_random, 'Random 3x3'),
                  (A_range, 'Arange 3x4')]:
    matrix_info(mat, name)

# Demonstrate rank deficiency
A_rank_deficient = np.array([[1, 2, 3],
                               [2, 4, 6],   # row 2 = 2 * row 1
                               [0, 1, 1]])
matrix_info(A_rank_deficient, 'Rank-deficient')

# --- Experiment 1: Where do the columns of A live after transformation? ---
# Hypothesis: The columns of A are exactly where the standard basis vectors land.
# Try changing: modify A and see where e1, e2 land.
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-whitegrid')

A = np.array([[2.0, -1.0],    # <-- modify this
              [1.0,  2.0]])

e1 = np.array([1.0, 0.0])    # standard basis vector 1
e2 = np.array([0.0, 1.0])    # standard basis vector 2

Ae1 = A @ e1   # should equal first column of A
Ae2 = A @ e2   # should equal second column of A

print(f"A[:,0] = {A[:,0]}   A @ e1 = {Ae1}   Equal: {np.allclose(A[:,0], Ae1)}")
print(f"A[:,1] = {A[:,1]}   A @ e2 = {Ae2}   Equal: {np.allclose(A[:,1], Ae2)}")

fig, axes = plt.subplots(1, 2, figsize=(10, 5))
for ax, (pts, transformed, title) in zip(axes, [
    (np.array([e1, e2]), np.array([Ae1, Ae2]), 'Basis vectors → columns of A'),
]):
    pass  # visualization already shown in Section 3

print("\nConclusion: The columns of A tell you where the coordinate axes land.")
print("A matrix is completely determined by where it sends the basis vectors.")

# --- Experiment 2: Rank and the dimension of the output ---
# Hypothesis: A rank-k matrix maps all of R^n into a k-dimensional subspace.
# Try changing: K — the number of linearly independent rows used.
import numpy as np

K = 2  # <-- change to 1, 2, or 3

# Build a matrix where only K rows are independent
base_rows = np.random.randn(K, 4)
coeffs = np.random.randn(4, K)
A = coeffs @ base_rows   # shape (4, 4) but rank K

print(f"Matrix shape: {A.shape}")
print(f"Constructed rank: {K}")
print(f"NumPy computed rank: {np.linalg.matrix_rank(A)}")

# Apply to 1000 random vectors and examine where they land
vecs = np.random.randn(4, 1000)
outputs = A @ vecs   # shape (4, 1000)

# Singular values show how many dimensions are 'alive'
_, s, _ = np.linalg.svd(outputs)
print(f"\nSingular values of output set (only {K} should be non-negligible):")
print(np.round(s[:8], 4))

# --- Experiment 3: Sparsity and computation cost ---
# Hypothesis: Sparse matrices can represent the same linear map with fewer stored numbers.
# Try changing: DENSITY — fraction of nonzero entries.
import numpy as np
import time

N = 1000
DENSITY = 0.01   # <-- modify: 0.001, 0.01, 0.1, 1.0

# Dense matrix
A_dense = np.random.randn(N, N)
mask = np.random.rand(N, N) > DENSITY
A_dense[mask] = 0.0

v = np.random.randn(N)

t0 = time.perf_counter()
for _ in range(100): _ = A_dense @ v
t_dense = (time.perf_counter() - t0) / 100

nonzero_frac = np.count_nonzero(A_dense) / A_dense.size
print(f"Density: {nonzero_frac:.4f}  ({np.count_nonzero(A_dense)} nonzeros out of {A_dense.size})")
print(f"Dense matmul time: {t_dense*1000:.3f} ms")
print("\nNote: dense numpy does not benefit from sparsity; use scipy.sparse for real sparse ops.")

# --- Mini Project: Dataset as Matrix ---
# Problem: A dataset of student scores across 5 subjects for 50 students.
# Represent it as a matrix, then compute per-subject means, per-student
# total scores, and identify which two subjects are most correlated.
# Task: Use only matrix operations — no explicit loops.

import numpy as np

np.random.seed(0)
SUBJECTS = ['Math', 'Physics', 'CS', 'English', 'History']
N_STUDENTS = 50

# Generate correlated scores: Math/Physics/CS share a 'STEM aptitude' factor
stem_factor = np.random.randn(N_STUDENTS)
humanities_factor = np.random.randn(N_STUDENTS)
noise = np.random.randn(N_STUDENTS, 5) * 0.5

scores = np.column_stack([
    70 + 10 * stem_factor + noise[:, 0],       # Math
    68 + 10 * stem_factor + noise[:, 1],       # Physics
    72 + 10 * stem_factor + noise[:, 2],       # CS
    65 + 10 * humanities_factor + noise[:, 3], # English
    63 + 10 * humanities_factor + noise[:, 4], # History
])
scores = np.clip(scores, 0, 100)

print(f"Data matrix shape: {scores.shape}  ({N_STUDENTS} students × {len(SUBJECTS)} subjects)")
print()

# TODO 1: Compute per-subject mean scores using axis argument
subject_means = # YOUR CODE
print("Subject means:", dict(zip(SUBJECTS, np.round(subject_means, 1))))

# TODO 2: Compute per-student total score (sum across subjects)
student_totals = # YOUR CODE
print(f"Top student score: {student_totals.max():.1f}")
print(f"Bottom student score: {student_totals.min():.1f}")

# TODO 3: Compute the 5×5 correlation matrix
# Hint: center the data first, then use matrix multiplication
centered = scores - subject_means  # broadcast subtract means
# corr[i,j] = (col_i · col_j) / (||col_i|| * ||col_j||)
corr_matrix = # YOUR CODE

print("\nCorrelation matrix (subjects x subjects):")
print(np.round(corr_matrix, 2))
print("\nMost correlated pair should be Math-Physics or English-History.")