Chapter 152 — Matrix Representation

# --- Visualization: Building a matrix column by column ---
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-whitegrid')

# Define a linear transformation by specifying where basis vectors go
# e1 = [1,0] maps to [2, 1]  (first column)
# e2 = [0,1] maps to [-1, 2] (second column)
col1 = np.array([2.0, 1.0])
col2 = np.array([-1.0, 2.0])
A = np.column_stack([col1, col2])   # pack columns to get matrix

print(f"Matrix A:\n{A}")
print(f"A[:,0] = {A[:,0]}  (image of e1)")
print(f"A[:,1] = {A[:,1]}  (image of e2)")

# Verify: apply A to an arbitrary vector and compare with linear combination
v = np.array([3.0, 2.0])
via_matmul = A @ v
via_combination = v[0] * col1 + v[1] * col2
print(f"\nA @ [3,2] via matmul:      {via_matmul}")
print(f"3*col1 + 2*col2 via combo: {via_combination}")
print(f"Equal: {np.allclose(via_matmul, via_combination)}")

# Visualize
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Before transformation
ax = axes[0]
for vec, col, label in [([1,0], 'crimson', 'e₁=[1,0]'), ([0,1], 'darkgreen', 'e₂=[0,1]'), (v, 'steelblue', f'v=[{v[0]:.0f},{v[1]:.0f}]')]:
    vec = np.array(vec)
    ax.annotate('', xy=vec, xytext=(0,0), arrowprops=dict(arrowstyle='->', color=col, lw=2))
    ax.text(vec[0]+0.05, vec[1]+0.05, label, fontsize=9)
ax.set_xlim(-1, 4); ax.set_ylim(-1, 4); ax.set_aspect('equal')
ax.axhline(0, color='k', lw=0.5); ax.axvline(0, color='k', lw=0.5)
ax.set_title('Original Space', fontsize=11)

# After transformation
ax = axes[1]
Av = A @ v
for vec, col, label in [(col1, 'crimson', f'A·e₁={col1}'), (col2, 'darkgreen', f'A·e₂={col2}'), (Av, 'steelblue', f'A·v={Av}')]:
    ax.annotate('', xy=vec, xytext=(0,0), arrowprops=dict(arrowstyle='->', color=col, lw=2))
    ax.text(vec[0]+0.05, vec[1]+0.05, label, fontsize=9)
ax.set_xlim(-2, 8); ax.set_ylim(-1, 8); ax.set_aspect('equal')
ax.axhline(0, color='k', lw=0.5); ax.axvline(0, color='k', lw=0.5)
ax.set_title('Transformed Space (A applied)', fontsize=11)

plt.suptitle('A Matrix is Built from Where Basis Vectors Land', fontsize=12)
plt.tight_layout()
plt.show()

# --- Implementation: Indexing, slicing, and memory layout ---
import numpy as np

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

# Basic indexing
print("=== Indexing ===")
print(f"A[1, 2] = {A[1, 2]}             # row 1, col 2")
print(f"A[0, :] = {A[0, :]}          # first row")
print(f"A[:, 1] = {A[:, 1]}          # second column")
print(f"A[0:2, 0:2] =\n{A[0:2, 0:2]}  # top-left 2x2 block")
print()

# Memory layout
print("=== Memory Layout ===")
A_C = np.array(A, order='C')   # row-major (default)
A_F = np.array(A, order='F')   # column-major
print(f"C order raw buffer: {A_C.flatten(order='K')}")
print(f"F order raw buffer: {A_F.flatten(order='K')}")
print(f"A_C is C-contiguous: {A_C.flags['C_CONTIGUOUS']}")
print(f"A_F is F-contiguous: {A_F.flags['F_CONTIGUOUS']}")
print()

# Reconstructing a matrix from column vectors
def matrix_from_column_images(images):
    """
    Build a matrix from the list of where basis vectors map.

    Args:
        images: list of 1D arrays, each the image of one basis vector

    Returns:
        2D numpy array, shape (len(images[0]), len(images))
    """
    return np.column_stack(images)

# Define T: e1 -> [1,0,-1], e2 -> [2,1,0], e3 -> [0,3,1]
T_mat = matrix_from_column_images([[1,0,-1], [2,1,0], [0,3,1]])
print(f"Matrix from column images:\n{T_mat}")

# --- Implementation: Recover a transformation matrix from function evaluation ---
import numpy as np

def linear_map_to_matrix(T_func, input_dim, output_dim=None):
    """
    Given a linear function T: R^n -> R^m, compute its matrix.
    Uses the fundamental theorem: column j = T(e_j).

    Args:
        T_func: callable, accepts 1D array of length input_dim
        input_dim: n, dimension of input space
        output_dim: m (inferred if None)

    Returns:
        2D numpy array of shape (output_dim, input_dim)
    """
    columns = []
    for j in range(input_dim):
        e_j = np.zeros(input_dim)
        e_j[j] = 1.0                 # j-th standard basis vector
        columns.append(T_func(e_j))  # T(e_j) is the j-th column
    return np.column_stack(columns)

# Test: define a transformation as a Python function
def my_transform(v):
    """Rotate and scale: x' = 2x - y, y' = x + y, z' = 3z"""
    return np.array([
        2*v[0] - v[1],
        v[0] + v[1],
        3*v[2]
    ])

M = linear_map_to_matrix(my_transform, input_dim=3)
print(f"Recovered matrix:\n{M}")

# Verify on a test vector
test_v = np.array([1.0, 2.0, 3.0])
print(f"\nDirect function:  {my_transform(test_v)}")
print(f"Via matrix:       {M @ test_v}")
print(f"Match: {np.allclose(my_transform(test_v), M @ test_v)}")

# --- Experiment: Row access vs column access speed ---
# Hypothesis: In row-major storage, row access is faster than column access.
# Try changing: N to see how the gap scales.
import numpy as np
import time

N = 2000   # <-- try 500, 2000, 5000
A = np.random.randn(N, N)
A_F = np.asfortranarray(A)

REPS = 200

# Row access on C-order
t0 = time.perf_counter()
for _ in range(REPS):
    for i in range(0, N, N//10):
        _ = A[i, :].sum()
t_row_C = (time.perf_counter() - t0) / REPS

# Column access on C-order
t0 = time.perf_counter()
for _ in range(REPS):
    for j in range(0, N, N//10):
        _ = A[:, j].sum()
t_col_C = (time.perf_counter() - t0) / REPS

print(f"N = {N}")
print(f"Row access (C-order):    {t_row_C*1000:.3f} ms")
print(f"Column access (C-order): {t_col_C*1000:.3f} ms")
print(f"Column/Row ratio: {t_col_C/t_row_C:.2f}x")
print("\nFor small N the gap is negligible; for large N, memory layout matters.")

# --- Mini Project: Image as a Matrix ---
# Problem: A grayscale image is literally a 2D matrix of pixel intensities.
# Task: Generate a synthetic image, then access specific rows/columns/patches
# and show how matrix indexing maps to spatial image operations.

import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-whitegrid')

# Generate a synthetic 64x64 'image' with structure
SIZE = 64
x = np.linspace(0, 4*np.pi, SIZE)
y = np.linspace(0, 4*np.pi, SIZE)
X, Y = np.meshgrid(x, y)
image = (np.sin(X) * np.cos(Y) + 1) / 2   # range [0, 1]

fig, axes = plt.subplots(2, 3, figsize=(13, 8))

axes[0,0].imshow(image, cmap='gray', vmin=0, vmax=1)
axes[0,0].set_title('Full Image (64x64 matrix)')

# Row 20: a horizontal slice of the image
axes[0,1].plot(image[20, :], color='steelblue')
axes[0,1].set_title('Row 20: A[20, :]')
axes[0,1].set_xlabel('Column index'); axes[0,1].set_ylabel('Intensity')

# Column 30: a vertical slice
axes[0,2].plot(image[:, 30], color='crimson')
axes[0,2].set_title('Column 30: A[:, 30]')
axes[0,2].set_xlabel('Row index'); axes[0,2].set_ylabel('Intensity')

# Top-left 16x16 patch
patch = image[:16, :16]
axes[1,0].imshow(patch, cmap='gray', vmin=0, vmax=1)
axes[1,0].set_title('Top-left patch: A[:16, :16]')

# Every other row and column (downsampling)
downsampled = image[::2, ::2]
axes[1,1].imshow(downsampled, cmap='gray', vmin=0, vmax=1)
axes[1,1].set_title(f'Downsampled 2x: A[::2, ::2]\n({downsampled.shape[0]}x{downsampled.shape[1]})')

# Flipped image (reverse rows)
flipped = image[::-1, :]
axes[1,2].imshow(flipped, cmap='gray', vmin=0, vmax=1)
axes[1,2].set_title('Flipped vertically: A[::-1, :]')

for ax in axes.flat: ax.axis('off') if ax.get_images() else None
plt.suptitle('Matrix Indexing = Image Spatial Operations', fontsize=13)
plt.tight_layout()
plt.show()