NumPy Linear Algebra — dot, matmul, linalg explained

# io.thecodeforge: Matrix multiplication — @ vs dot vs matmul
import numpy as np

# ── 2D ARRAYS: ALL THREE ARE EQUIVALENT ───────────────────────────────────────
A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])

result_at     = A @ B           # preferred — unambiguous, matches math notation
result_matmul = np.matmul(A, B) # explicit call to the same function as @
result_dot    = np.dot(A, B)    # identical for 2D; dangerous for 3D+

print(result_at)
# [[19 22]
#  [43 50]]
print(np.array_equal(result_at, result_matmul))  # True
print(np.array_equal(result_at, result_dot))     # True — only because these are 2D

# ── 1D VECTORS: DOT PRODUCT ────────────────────────────────────────────────────
u = np.array([1, 2, 3])
v = np.array([4, 5, 6])

print(u @ v)         # 32 — scalar dot product, clean and clear
print(np.dot(u, v))  # 32 — identical for 1D, np.dot is acceptable here

# ── 3D ARRAYS: WHERE @ AND DOT DIVERGE ────────────────────────────────────────
# Simulates batched embedding similarity: batch of user vectors against item matrices
batch_size = 10
batch_A = np.random.randn(batch_size, 3, 4)  # 10 matrices of shape (3, 4)
batch_B = np.random.randn(batch_size, 4, 5)  # 10 matrices of shape (4, 5)

# @ treats as stacked matrices — multiplies batch_A[i] @ batch_B[i] for each i
result_matmul = batch_A @ batch_B
print(f"@ shape:   {result_matmul.shape}")  # (10, 3, 5) — correct batched result

# np.dot contracts last axis of A with second-to-last of B — cross-batch product
result_dot = np.dot(batch_A, batch_B)
print(f"dot shape: {result_dot.shape}")     # (10, 3, 10, 5) — WRONG: cross-batch

# Verify that @ computes batch_A[i] @ batch_B[i] independently
for i in range(batch_size):
    assert np.allclose(result_matmul[i], batch_A[i] @ batch_B[i]), \
        f"Batch element {i} mismatch"
print("All batch elements verified: @ is correctly independent per element")

# ── SHAPE ASSERTION PATTERN: MANDATORY IN PRODUCTION ──────────────────────────
def batched_similarity(user_emb: np.ndarray, item_emb: np.ndarray) -> np.ndarray:
    """
    user_emb: (batch_size, embed_dim)
    item_emb: (batch_size, num_items, embed_dim)
    returns:  (batch_size, num_items)
    """
    # Transpose item_emb last two axes for matrix multiplication
    scores = user_emb[:, np.newaxis, :] @ item_emb.transpose(0, 2, 1)
    scores = scores.squeeze(1)  # (batch_size, num_items)
    assert scores.shape == (user_emb.shape[0], item_emb.shape[1]), \
        f"Similarity shape wrong: got {scores.shape}"
    return scores

# io.thecodeforge: Solving linear systems — solve vs inverse, condition numbers
import numpy as np
import time

# ── BASIC LINEAR SYSTEM ───────────────────────────────────────────────────────
# System: 2x + y = 5
#          x + 3y = 10
A = np.array([[2.0, 1.0],
              [1.0, 3.0]])
b = np.array([5.0, 10.0])

x = np.linalg.solve(A, b)
print(f"Solution: x={x[0]:.4f}, y={x[1]:.4f}")  # x=1.0, y=3.0
print(f"Residual ||Ax - b||: {np.linalg.norm(A @ x - b):.2e}")  # essentially 0
print(f"Verification A @ x ≈ b: {np.allclose(A @ x, b)}")       # True

# ── PERFORMANCE: solve vs inv ─────────────────────────────────────────────────
np.random.seed(42)
A_large = np.random.randn(1000, 1000)
# Make it positive definite to ensure it is well-conditioned
A_large = A_large @ A_large.T + 1000 * np.eye(1000)
b_large = np.random.randn(1000)

runs = 10
start = time.perf_counter()
for _ in range(runs):
    x_solve = np.linalg.solve(A_large, b_large)
time_solve = (time.perf_counter() - start) / runs

start = time.perf_counter()
for _ in range(runs):
    x_inv = np.linalg.inv(A_large) @ b_large
time_inv = (time.perf_counter() - start) / runs

print(f"\nPerformance on 1000×1000 system (mean of {runs} runs):")
print(f"  linalg.solve:  {time_solve*1000:.1f}ms")
print(f"  inv(A) @ b:    {time_inv*1000:.1f}ms")
print(f"  Speedup:       {time_inv/time_solve:.1f}x")
print(f"  Results match: {np.allclose(x_solve, x_inv, rtol=1e-10)}")

# ── CONDITION NUMBER AND ILL-CONDITIONING ─────────────────────────────────────
A_well = np.eye(3)                                    # condition number = 1.0
A_mild = np.array([[1.0, 0.99], [0.99, 1.0]])         # mildly ill-conditioned
A_bad  = np.array([[1.0, 1.0], [1.0, 1.0 + 1e-10]])  # nearly singular

print(f"\nCondition numbers:")
print(f"  Identity matrix:    {np.linalg.cond(A_well):.2e}")   # 1.00e+00
print(f"  Mildly ill-cond:    {np.linalg.cond(A_mild):.2e}")   # ~1.99e+02
print(f"  Nearly singular:    {np.linalg.cond(A_bad):.2e}")    # ~2.00e+10

# ── HANDLING ILL-CONDITIONED SYSTEMS ──────────────────────────────────────────
# lstsq returns the minimum-norm least-squares solution — no exception on singular input
A_singular = np.array([[1.0, 2.0], [2.0, 4.0]])  # rank 1 — singular
b_singular  = np.array([3.0, 6.0])

try:
    np.linalg.solve(A_singular, b_singular)  # raises LinAlgError
except np.linalg.LinAlgError as e:
    print(f"\nlinalg.solve raised: {e}")

x_ls, residuals, rank, sv = np.linalg.lstsq(A_singular, b_singular, rcond=None)
print(f"lstsq solution:      {x_ls}")  # minimum-norm solution
print(f"Matrix rank:         {rank}")   # 1 — correctly identified as rank-deficient
print(f"Singular values:     {sv}")

# ── MULTIPLE RIGHT-HAND SIDES: FACTORIZE ONCE ─────────────────────────────────
from scipy.linalg import lu_factor, lu_solve

A_reuse = np.random.randn(500, 500)
A_reuse = A_reuse @ A_reuse.T + 500 * np.eye(500)  # well-conditioned SPD
lu, piv = lu_factor(A_reuse)  # factorize once

for i in range(100):           # solve 100 different right-hand sides
    b_i = np.random.randn(500)
    x_i = lu_solve((lu, piv), b_i)  # reuse factorization — fast
print("\n100 right-hand sides solved with one LU factorization.")

# io.thecodeforge: Eigenvalues, eigenvectors, SVD, and PCA from scratch
import numpy as np

# ── EIGENDECOMPOSITION: GENERAL MATRIX ────────────────────────────────────────
A = np.array([[4.0, 2.0],
              [1.0, 3.0]])

values, vectors = np.linalg.eig(A)
print(f"Eigenvalues:  {values}")     # [5. 2.]
print(f"Eigenvectors:\n{vectors}")   # columns are eigenvectors

# Verify: A @ v = lambda * v for each eigenvector
for i in range(len(values)):
    lhs = A @ vectors[:, i]
    rhs = values[i] * vectors[:, i]
    print(f"  λ={values[i]:.1f}: A@v matches λv: {np.allclose(lhs, rhs)}")

# ── EIGENDECOMPOSITION: SYMMETRIC MATRIX (use eigh) ───────────────────────────
S = np.array([[4.0, 1.0],
              [1.0, 3.0]])

# eigh: faster, real output, eigenvalues in ascending order
vals_h, vecs_h = np.linalg.eigh(S)
print(f"\neigh eigenvalues (ascending): {vals_h}")  # [2.76 4.24] — always real

# Compare: eig on the same symmetric matrix may show tiny imaginary parts
vals_g, _ = np.linalg.eig(S)
print(f"eig  eigenvalues: {vals_g}")   # may show 0j imaginary parts
print(f"Symmetric check: {np.allclose(S, S.T)}")  # True — safe to use eigh

# ── SVD ────────────────────────────────────────────────────────────────────────
np.random.seed(42)
M = np.random.randn(100, 20)  # 100 observations, 20 features

# full_matrices=False: thin SVD — U is (100, 20) not (100, 100)
U, s, Vt = np.linalg.svd(M, full_matrices=False)
print(f"\nThin SVD shapes: U={U.shape}, s={s.shape}, Vt={Vt.shape}")
# U=(100,20), s=(20,), Vt=(20,20)

# Reconstruct and verify
M_reconstructed = U @ np.diag(s) @ Vt
print(f"Reconstruction error: {np.linalg.norm(M - M_reconstructed):.2e}")  # ~1e-13

# ── LOW-RANK APPROXIMATION ─────────────────────────────────────────────────────
# Eckart-Young theorem: rank-k truncation minimizes Frobenius norm error
for k in [1, 5, 10, 20]:
    M_k = U[:, :k] @ np.diag(s[:k]) @ Vt[:k, :]
    error = np.linalg.norm(M - M_k, 'fro')
    variance_explained = (s[:k]**2).sum() / (s**2).sum()
    print(f"  Rank-{k:2d}: error={error:.3f}, variance explained={variance_explained:.1%}")

# ── PCA FROM SCRATCH USING SVD ────────────────────────────────────────────────
np.random.seed(0)
X = np.random.randn(200, 5)  # 200 samples, 5 features
X[:, 1] += 2 * X[:, 0]      # introduce correlation to make PCA interesting

# Step 1: center the data
X_centered = X - X.mean(axis=0)

# Step 2: thin SVD of centered data
U_pca, s_pca, Vt_pca = np.linalg.svd(X_centered, full_matrices=False)

# Step 3: principal components are the rows of Vt_pca
print(f"\nPCA principal components shape: {Vt_pca.shape}")  # (5, 5)

# Step 4: explained variance ratio
explained_variance_ratio = (s_pca**2) / (s_pca**2).sum()
print("Explained variance ratio per component:")
for i, evr in enumerate(explained_variance_ratio):
    print(f"  PC{i+1}: {evr:.1%}")

# Step 5: project to top-2 principal components
X_pca = X_centered @ Vt_pca[:2].T  # (200, 2)
print(f"Projected data shape: {X_pca.shape}")  # (200, 2)

# ── MATRIX NORM ────────────────────────────────────────────────────────────────
print(f"\nFrobenius norm of A: {np.linalg.norm(A, 'fro'):.4f}")
print(f"Spectral norm of A:  {np.linalg.norm(A, 2):.4f}")   # largest singular value
print(f"Nuclear norm of M:   {s.sum():.4f}")                  # sum of singular values