Linear Algebra II — Matrices, Transforms, Eigenvalues

Review boundary and corrected contracts

Small determinant magnitude is not a scale-invariant condition-number test: εI has tiny determinant and condition number one. solve avoids an explicit inverse but cannot remove an ill-conditioned problem. SVD complexity depends on sparsity, solver and accuracy. Eigenvectors keep an invariant line; a negative eigenvalue reverses direction.

The examples below teach the stated input contracts; fragments and deliberately incorrect counterexamples are exercises, not production implementations. Complexity counts elementary operations unless a bit-cost or allocation caveat is stated. Practice times are suggestions, not measured performance or hiring guarantees.

Learning objective: Read a matrix as a linear transformation, compose transformations by multiplying matrices, and know what eigenvalues/eigenvectors tell you about the transform's behaviour.

Why this session exists

A neural network is nothing but stacked matrix multiplications and non-linearities. A rendered 3-D scene is a sequence of matrix transformations. PCA, SVD, PageRank, JPEG, quantum states — all matrices. After choosing bases, every linear map between finite-dimensional vector spaces has a matrix representation. Abstract linear maps and infinite-dimensional operators need not be written as finite arrays. Get past "matrix = grid of numbers" to "matrix = function that stretches, rotates, and squashes space" and every downstream ML paper stops feeling like a wall of Greek letters.

You will be able to

• Read a matrix as a function that maps input vectors to output vectors.

• Multiply two matrices by hand for 2×2 and 3×3 cases (and know why row × column).

• Recognise identity, transpose, inverse, and diagonal — and say what each does geometrically.

• Explain eigenvalues and eigenvectors as the ‘axes the transform respects’ with a picture.

• Use NumPy to multiply matrices, solve linear systems, and find eigenvalues.

Prerequisites

• S017 · Vectors — dot product IS the atom of matrix multiplication.
• Comfortable with element-wise math and Σ notation.

Pre-read · skim before the session (35 min total)

Further reading: Mathematics for Machine Learning (Ch. 2.3–2.7 and 4.1–4.3)

Further reading: NumPy linalg reference

Further reading: Matrix multiplication

(a) Intuition · 5 min

A matrix is a function that transforms space

Real: Imagine a rubber sheet with a grid drawn on it. You pinch two corners and stretch it — the grid deforms. Lines map to lines or points, and the origin stays put; a singular transformation can collapse directions. That constrained stretch is a linear transformation, and it is entirely described by where the two basis vectors [1,0] and [0,1] end up. Those two ending points are the two columns of a 2×2 matrix. Multiply the matrix by any point on the original sheet and it tells you where that point moved to. Rotate the sheet 90°? Matrix. Stretch it 2× vertically? Matrix. Shear it? Matrix.

 
    In code, a matrix is a 2-D array. But its meaning is a function. A @ x in NumPy takes vector x and returns Ax — the transformed vector. Every operation you'll use — matrix-matrix multiply, inverse, eigendecomposition — is a fact about that function.
    The columns of A are literally "where the basis vectors go." That's the one sentence that unlocks the whole subject.
  

The three views of a matrix

A matrix is always all three

• A grid of numbers — m rows × n columns. This is the storage view.

• A linear function ℝⁿ → ℝᵐ — takes a length-n vector, returns a length-m vector. This is the transformation view.

• A composition of column vectors — the k-th column of A is where the k-th basis vector of the input space lands. This is the ‘how do I build one?’ view.

Matrix multiplication is not commutative

A · B ≠ B · A in general. Rotating then scaling ≠ scaling then rotating. If you write A @ B where you meant B @ A, the code still runs and produces wrong answers silently. Always check the shapes: (m×k) @ (k×n) = (m×n). If the inner dimensions don't match, you'll get an error — that's the friendly case. If the dimensions also work in the reversed order, the code may silently compute a different transformation. Some special matrices commute; the general rule does not promise inequality for every pair.

A quick history

In Arthur Cayley's 1858 memoir, the opening page connects matrix notation to systems of linear equations and distinguishes matrix addition from order-sensitive multiplication (called composition). That is the same bridge used here: a compact notation becomes an algebra of transformations. Modern numerical libraries add a different concern—how to compute those operations accurately with finite precision. NumPy's documented LAPACK-backed solvers illustrate that computational layer; a historical algebraic identity is not by itself a stable numerical algorithm.

(b) Follow the basis vectors through a transformation · 15 min

The transformation view gives matrix arithmetic a concrete test: where does each input arrow end up? Follow the two basis vectors first, use their destinations to understand composition, then look for directions that the transformation only scales.

A matrix transforms space

The first two paths apply the same matrix to the two coordinate axes. Match each output pair to a column of the matrix; the third path explains why knowing those two destinations is enough for every other input vector.

The columns of A are exactly where the basis vectors land. Every other vector's destination is a weighted combination of those two column-destinations.

The matrix multiplication formula, unpacked

A second matrix supplies several input columns at once. Applying the first matrix to each of them produces the columns of the product; the row-by-column rule is how you compute each coordinate of those destinations.

Check shapes

Desc: (m × k) · (k × n) → (m × n). Inner k must match; outer m and n are the output shape.

Row × column

Desc: Entry (i, j) of the result = dot product of row i of A with column j of B.

For example, [[1,2],[3,4]] @ [[2,0],[1,2]] = [[4,4],[10,8]]: entry (2,1) is 3·2 + 4·1 = 10 using one-based labels. In a 3×3 case, [[1,2,3],[0,1,4],[5,6,0]] @ [[1,0,2],[0,1,3],[1,1,0]] = [[4,5,8],[4,5,3],[5,6,28]]. Work through each row-column pair; the output shape is not inferred from elementwise multiplication.

Result is a matrix of dot products

Desc: Every cell of the product is a scalar; it's the ‘answer’ to ‘what does row i of A × column j of B evaluate to?’

Alternative view · composition

Desc: A · B is the function ‘first apply B, then apply A.’ Order matters and reads right-to-left.

The special matrices worth memorising

Six matrices you'll see every day

Name: Identity Iₙ

Desc: 1s on the diagonal, 0s elsewhere. A · Iₙ = Iₘ · A = A for an m×n matrix. The ‘do nothing’ transform.

Name: Diagonal D

Desc: Non-zero only on the diagonal. Applying an n×n diagonal matrix to a vector scales each coordinate independently. Using its diagonal directly costs O(n), rather than O(n²) for a generic dense matrix-vector product.

Name: Transpose Aᵀ

Desc: Swap rows and columns. Preserves dot products in a sense: (Ax)·y = x·(Aᵀy). Central to gradient formulas.

Name: Inverse A⁻¹

Desc: The transform that undoes A. A · A⁻¹ = I. Only exists if A is square and det(A) ≠ 0.

Name: Orthogonal Q

Desc: Rows and columns are unit vectors, mutually perpendicular. Q · Qᵀ = I. Represents rotations/reflections — preserves lengths and angles.

Name: Symmetric S

Desc: For real S = Sᵀ, all eigenvalues are real and an orthonormal eigenbasis exists. Eigenvectors for distinct eigenvalues are orthogonal; arbitrary vectors chosen within a repeated-eigenvalue eigenspace need not be. Covariance matrices are real symmetric, which supports the eigenbasis used in PCA.

Solving Ax = b: three interpretations, one answer

So far the input has been known and the output computed. Solving a linear system reverses the question: given the destination, which input could have reached it? The row, column, and function views below describe the same search, with an inverse available only when the transformation is invertible.

Name: Row picture

Tagline: Intersection of hyperplanes

• Each row of A is a linear equation

• Each equation is a hyperplane in ℝⁿ

• Solution x is where they all meet

• Sensitive to picture-blowing-up in high dims

Name: Column picture

Tagline: Linear combination of columns

• b is a linear combination of columns of A

• Solve for the weights

• Cleaner geometric intuition

• Preferred view for ML

Name: Function picture

Tagline: Undo the transform

• A maps x → b

• We want x such that A applies to it gives b

• x = A⁻¹ b (when A is invertible)

• Numerical: prefer np.linalg.solve(A, b) for a nonsingular square system; use least squares when appropriate

Eigenvalues and eigenvectors — geometric meaning

Instead of fixing a destination, ask whether any nonzero input stays on its own line through the origin after the transformation. The diagram contrasts a vector that turns with two special vectors that are only scaled; identifying those special directions can make the transformation much easier to understand.

Eigenvector: a nonzero vector v with Av = λv. Its line maps into itself, possibly collapsing to zero. Eigenvalue λ: the signed scale factor over the reals.

The graph illustrates possible cases, not a promise that every real matrix has two real eigenvectors. An eigenvector's line maps into itself: a positive eigenvalue preserves its orientation, a negative eigenvalue reverses it, and zero sends the vector to the origin. A scalar multiple of the identity has every nonzero vector as an eigenvector, so a randomly chosen vector need not change direction.

Why eigenvalues matter in ML

A real symmetric matrix has an orthonormal eigenbasis. PCA orders the covariance eigenvectors by decreasing variance; after selecting one, the next maximizes variance among directions orthogonal to those already selected. PageRank uses a stationary probability vector, with row/column conventions determining which side of its stochastic matrix it multiplies. Singular values quantify how a linear map stretches vector lengths: the largest bounds amplification and a very small relative smallest value signals sensitivity. Layer norms can bound a network Jacobian together with its nonlinearities, but individual weight eigenvalues alone do not guarantee stable training.

Basis-vector destinations, composition, and eigenvectors all describe what the matrix does, not merely how it is stored. That distinction also explains the next failure case: a transformation that flattens a dimension cannot uniquely reconstruct its input, however carefully you implement the arithmetic.

Myth: "A matrix is a 2D array of numbers — a grid, like a spreadsheet or a table. Matrix multiplication is a slightly odd rule you memorise: rows times columns."

Truth: A matrix is a linear transformation written in coordinates. Its columns are where the basis vectors land. Matrix multiplication is function composition — AB means "do B, then A" — and every property that looks arbitrary (why the inner dimensions must match, why it isn't commutative, why the identity looks like that) falls out of that single fact.

Why: Because the grid model is sufficient for storage and for the mechanics of computing a product, and because in data work you genuinely do use matrices as tables of records. The model fails the moment you need to reason rather than compute. Under the grid view, AB ≠ BA is a strange quirk to remember; under the transformation view it is obvious, because rotating-then-stretching is visibly not the same as stretching-then-rotating. Similarly "singular" sounds like a numerical accident, when it actually means the transformation collapsed a dimension and destroyed information — which is why it cannot be inverted, and why a poorly conditioned design can make regression coefficients sensitive to small data changes. Everything mysterious about linear algebra becomes mechanical once you stop seeing grids.

Test: The columns literally are the images of the basis vectors — check it:import numpy as np A = np.array([[0., -1.], [1., 0.]]) # 90 degree rotation print(A @ np.array([1., 0.])) # [0, 1] <- column 1 of A print(A @ np.array([0., 1.])) # [-1, 0] <- column 2 of A

B = np.array([[2., 0.], [0., 1.]]) # stretch x by 2 print(A @ B) print(B @ A) # different: order is composition order

Question: Why does a matrix fail to be invertible exactly when its determinant is zero? These look like unrelated facts — one about solving equations, one about a strange alternating sum.

Premise: A matrix maps input vectors to output vectors. Inverting it means recovering the unique input that produced a given output.

Because: an inverse function must be well defined, so each output needs exactly one preimage

Premise: The determinant measures how the transformation scales volume: the unit cube maps to a parallelepiped whose signed volume is det(A).

Because: linearity means the volume scaling factor is the same everywhere, so one number suffices

Premise: If det(A) = 0, the output volume is zero — the entire input space has been flattened onto a lower-dimensional subspace (a plane, a line, a point).

Because: only a degenerate shape has zero volume, which means the columns became linearly dependent

Premise: Flattening is many-to-one: infinitely many distinct inputs map to the same output, because an entire direction was crushed to zero.

Because: rank-nullity gives a nonzero null vector v, so A(x + tv) = Ax for every real t

Premise: Therefore no inverse function can exist. There is no rule that recovers which of infinitely many preimages you started from; the information was genuinely destroyed, not merely hidden.

Because: an inverse must map each output back to a single input, and there isn't one

The magnitude of the determinant is scale-dependent. For example εI has determinant εⁿ but condition number one, so determinant magnitude alone is not a stability test. Analyze singular values/condition number. For λ>0, XᵀX+λI is positive definite, but regularization changes the estimation problem.

Name: Matrices are verbs, not tables

Picture: Read every matrix as an action on space. Its columns tell you where the basis vectors go, and because the transformation is linear, that's enough to determine where everything goes. Multiplying Ax is applying the verb to a vector; multiplying AB is chaining two verbs, applied right to left.Once matrices are verbs, the vocabulary decodes itself: the identity is "do nothing", the inverse is "undo", the determinant is "how much does this squash or expand", singular is "collapsed a dimension", transpose swaps the roles of input and output space, and real eigenvectors span lines that the map sends into themselves, including sign reversal or collapse.

• Inner dimensions must match because the output space of the first transformation must be the input space of the second. Shape errors are type errors about which space you're in.

• AB ≠ BA because composition order matters. This is not a quirk; it is why the order of operations in a neural network layer stack is meaningful.

• Prefer a solver for Ax = b rather than computing an inverse solely to multiply it by b. This avoids unnecessary work and rounding steps. Neither route removes the sensitivity of an ill-conditioned problem; an inverse is appropriate when the inverse itself is the required result.

• Eigen/SVD decompositions answer "what does this transformation do, in its own natural coordinates?" — which is why SVD underlies PCA, low-rank approximation, and recommender factorisation alike.

Trigger: Fire this model the moment you see: a shape mismatch error · LinAlgError: Singular matrix · regression coefficients that are huge or flip sign between runs · a condition number in a warning · PCA, SVD, or embedding factorisation · a transformer's attention as QK'.

Decision: You need to solve a large least-squares problem Ax ≈ b. Normal equations, QR/SVD decomposition, or iterative methods?

Name: Normal equations (A'A x = A'b)

Gain: simplest to derive and implement; A'A is small when features are far fewer than rows, and it can be accumulated in one streaming pass over the data

Cost: for full-column-rank A, κ₂(AᵀA) = κ₂(A)². This can substantially worsen rounding-error sensitivity; it does not mean every fit loses exactly half its digits or that QR cures ill-conditioning.

Pick: features are few, well-conditioned, and you need a single pass over data that doesn't fit in memory

Name: QR or SVD decomposition

Gain: numerically stable — it works directly on A, so the condition number is never squared; SVD additionally reveals rank and provides the pseudo-inverse for genuinely rank-deficient problems

Cost: needs the full matrix (or a blocked variant), higher cost than the normal equations, and SVD in particular is expensive on very large or very wide matrices

Pick: accuracy matters, the matrix might be ill-conditioned, or you need to diagnose rank deficiency rather than merely survive it — the default for any serious numerical work

Name: Iterative methods (conjugate gradient, LSQR, SGD)

Gain: methods such as LSQR use products with A and Aᵀ without a dense factorization. Conjugate gradient requires a symmetric positive-definite system, so it is not directly interchangeable with LSQR on arbitrary rectangular A. SGD uses sampled loss gradients; these methods have different convergence contracts.

Cost: convergence depends on conditioning and often requires a preconditioner to be practical; the answer is approximate with an error you must monitor; more knobs to get wrong

Pick: the matrix is huge and sparse, or exists only as an operator, or an approximate solution reached quickly beats an exact one reached slowly — which is the entire premise of training by gradient descent

Verdict: Use a library's least-squares routine (which will pick QR or SVD) unless scale forces you elsewhere, and check the condition number before trusting any coefficient. The recurring production mistake is not choosing the wrong solver — it's reporting coefficients from an ill-conditioned fit as though they were meaningful, when collinearity has made their individual values essentially arbitrary while the predictions remain fine.That distinction is worth holding onto: a model can predict well on similar inputs while individual coefficients are unstable. Prediction can also fail under perturbations or extrapolation, so neither is automatically protected.

(c) Hands-on · 25 min

You'll implement matrix multiplication from scratch, compare to NumPy, solve a toy linear system, and find principal directions of a synthetic 2-D dataset. Save as matrices_demo.py. The teaching helpers expect nonempty rectangular numeric matrices with compatible dimensions; assert checks are diagnostics, not production input validation (optimized Python can remove them).

"""matrices_demo.py — matrices from scratch + NumPy + a mini-PCA."""
from __future__ import annotations
import numpy as np
 
Matrix = list[list[float]]
Vector = list[float]
 
# ── 1. Matrix-vector and matrix-matrix multiplication from scratch ───────
def matvec(A: Matrix, x: Vector) -> Vector:
    """Ax  where A is m×n and x is length n."""
    m, n = len(A), len(A[0])
    assert len(x) == n, f"shape mismatch: A is {m}x{n}, x length is {len(x)}"
    return [sum(A[i][j] * x[j] for j in range(n)) for i in range(m)]
 
def matmul(A: Matrix, B: Matrix) -> Matrix:
    """AB  where A is m×k and B is k×n → result m×n."""
    m, k = len(A), len(A[0])
    k2, n = len(B), len(B[0])
    assert k == k2, f"inner dims must match: {k} != {k2}"
    return [
        [sum(A[i][p] * B[p][j] for p in range(k)) for j in range(n)]
        for i in range(m)
    ]
 
def transpose(A: Matrix) -> Matrix:
    return [list(row) for row in zip(*A)]
 
def identity(n: int) -> Matrix:
    return [[1.0 if i == j else 0.0 for j in range(n)] for i in range(n)]
 
# ── 2. Sanity vs NumPy on 200 random cases ──────────────────────────────
def compare_to_numpy() -> None:
    rng = np.random.default_rng(42)
    for _ in range(200):
        m, k, n = rng.integers(1, 6, size=3)
        A = rng.normal(size=(m, k))
        B = rng.normal(size=(k, n))
        ours = matmul(A.tolist(), B.tolist())
        np.testing.assert_allclose(ours, A @ B, atol=1e-9)
    print("200 random matrix products match NumPy within the stated tolerance ✅")
 
# ── 3. Composition: (rotate 90°) followed by (scale x by 2) ─────────────
def demo_composition() -> None:
    R = [[0, -1], [1, 0]]     # 90° CCW rotation
    S = [[2,  0], [0, 1]]     # scale x by 2
    v = [1, 0]
    Rv = matvec(R, v);   print("Rv     =", Rv)     # → (0, 1)
    SRv = matvec(S, Rv); print("S(Rv)  =", SRv)     # → (0, 1)
    # Compose: apply R first, then S  ⇔  (S · R) v
    SR = matmul(S, R)
    print("(S·R)v =", matvec(SR, v))                # same as above
 
# ── 4. Solving a real linear system: 3 equations, 3 unknowns ────────────
def demo_solve() -> None:
    """
    A restaurant sells apples, bananas, cherries.
    3 receipts:
       2 apples + 1 banana + 3 cherries = $8.60
       1 apple  + 4 bananas + 2 cherries = $9.30
       3 apples + 2 bananas + 1 cherry   = $7.70
    Find the per-fruit price.
    """
    A = np.array([[2, 1, 3],
                  [1, 4, 2],
                  [3, 2, 1]], dtype=float)
    b = np.array([8.60, 9.30, 7.70])
    x = np.linalg.solve(A, b)
    print(f"prices  apple={x[0]:.2f}  banana={x[1]:.2f}  cherry={x[2]:.2f}")
    # Verify: A @ x should reproduce b
    np.testing.assert_allclose(A @ x, b, atol=1e-9)
 
# ── 5. A tiny PCA — find the principal direction of a 2D cloud ──────────
def demo_pca() -> None:
    rng = np.random.default_rng(0)
    # Cloud stretched along the (1,2)/sqrt(5) direction
    n = 500
    t = rng.normal(scale=3.0, size=n)      # along principal direction
    noise = rng.normal(scale=0.4, size=n)  # perpendicular
    X = np.column_stack([t, 2*t]) + rng.normal(scale=0.4, size=(n, 2))
 
    X_centred = X - X.mean(axis=0)
    cov = (X_centred.T @ X_centred) / (n - 1)
    eigvals, eigvecs = np.linalg.eigh(cov)   # symmetric → eigh
    order = np.argsort(eigvals)[::-1]
    eigvals, eigvecs = eigvals[order], eigvecs[:, order]
 
    print(f"eigenvalues (variances): {eigvals.round(3)}")
    print(f"principal direction: {eigvecs[:, 0].round(3)}")
    # Expected: principal direction ≈ ±(1, 2)/sqrt(5) ≈ (0.447, 0.894)
 
if __name__ == "__main__":
    compare_to_numpy()
    print("\n── COMPOSITION ─────")
    demo_composition()
    print("\n── SOLVE Ax = b ────")
    demo_solve()
    print("\n── MINI PCA ────────")
    demo_pca()

Run it:

uv run --with numpy python matrices_demo.py

Observed output with the seeded example and this review's NumPy runtime (rounding and the eigenvector's overall sign may differ across implementations):

200 random matrix products match NumPy within the stated tolerance ✅

── COMPOSITION ─────
Rv     = [0, 1]
S(Rv)  = [0, 1]
(S·R)v = [0, 1]

── SOLVE Ax = b ────
prices  apple=1.22  banana=1.19  cherry=1.66

── MINI PCA ────────
eigenvalues (variances): [46.409  0.178]
principal direction: [0.451 0.892]

What each block does

Anatomy of the code

Name: matvec + matmul from scratch

Desc: Nested comprehensions expose each dot product. NumPy can use optimized native kernels; speed depends on shapes, dtype, hardware and library build, so no fixed ratio is implied.

Name: compare_to_numpy

Desc: Differential test on 200 seeded shape/value cases. assert_allclose checks abs(actual - desired) <= atol + rtol * abs(desired); here atol is 1e-9 and the NumPy 2.4 default rtol is 1e-7. Passing these cases is evidence, not a proof for all matrices.

Name: demo_composition

Desc: R first (rotate 90°), then S (scale x), is the composition S∘R and the matrix product S @ R. Function composition and matrix multiplication use the same order convention: the rightmost operation acts first.

Name: np.linalg.solve(A, b)

Desc: For a dense nonsingular n×n system, solve calls LAPACK's gesv solver (LU with pivoting), with O(n³) factorization work. Prefer it when you need x rather than A⁻¹; inspect conditioning and residuals because a stable algorithm can still solve a sensitive problem.

Name: np.linalg.eigh(cov)

Desc: eigh is for real symmetric or complex Hermitian matrices and returns real eigenvalues. The sample covariance here is real symmetric. It reads one selected triangle rather than validating symmetry of an arbitrary input.

Name: PCA in 5 lines

Desc: For this centered real dataset, form its sample covariance, diagonalize it, and sort by decreasing eigenvalue. The first eigenvector gives the maximum-variance axis; later axes are constrained to be orthogonal. The unused noise draw in this retained program advances the seeded generator; the noise actually added to X is the independent two-coordinate draw on the next line.

What to observe

ClassName: ml-4 list-decimal space-y-1

All 200 random products meet the combined absolute/relative tolerance above; this is not a pure absolute 1e-9 guarantee. The composition demo shows: applying two transformations is the same as multiplying their matrices — but READ RIGHT-TO-LEFT. The restaurant solver returns approximately $1.22 /$1.19 / $1.66. The exact decimal solution is 1.22, 1.192 and 1.656, which satisfies all three receipts; plugging rounded display prices back in introduces rounding error. The seeded sample's principal direction is approximately (0.451, 0.892), close to the population direction (1, 2)/√5. Its negative represents the same axis. The observed eigenvalues are approximately 46.409 and 0.178, so the first component captures about 99.62% of this sample's variance. That supports a low reconstruction error for this generated cloud, not a universal guarantee for a two-dimensional dataset.

Goal: Feel that rotation is orthogonal

Hint: Verify by computing R @ R.T — you should get the identity to machine precision. Orthogonal matrices satisfy Q · Qᵀ = I. This is the defining property. Now try it with S (the scale matrix) — S @ S.T is NOT identity; S is not orthogonal.

import numpy as np
R = np.array([[0.0, -1.0], [1.0, 0.0]])   # 90° rotation
print("R @ R.T =")
print(R @ R.T)                             # should be identity
print("det(R) =", np.linalg.det(R))        # should be 1
 
theta = np.pi / 4                          # 45°
R45 = np.array([[np.cos(theta), -np.sin(theta)],
                [np.sin(theta),  np.cos(theta)]])
v = np.array([1.0, 0.0])
print("R45 v =", R45 @ v)                  # should have length 1
print("length =", np.linalg.norm(R45 @ v)) # approximately 1 — floating arithmetic rounds

(d) Production reality · 15 min

Where this shows up in the rest of the plan

Root: Matrices are the atom of ML infrastructure

Label: S019/S020 · Calculus & gradient descent

Desc: Gradients are vectors; Jacobians and Hessians are matrices; chain rule is matrix multiplication.

Label: S078 · PCA & dimensionality reduction

Desc: SVD can produce a rank-50 approximation of 10 000-dimensional data. How much information it retains depends on the singular-value spectrum and the downstream task.

Label: S085 · Convolutional nets

Desc: A fixed convolution is a linear operator representable by a structured matrix. Implementations may use direct kernels, transformed algorithms or matrix-multiply lowering rather than materializing that matrix.

Label: S095 · Word embeddings

Desc: E[word_id] selects an embedding row. E · Eᵀ gives all pairwise dot products; for cosine similarities, normalize nonzero rows first.

Label: S115 · Transformers

Desc: Scaled dot-product attention is softmax(Q · Kᵀ / sqrt(d)) · V: two matrix products after Q, K and V have been formed. Their learned projections require additional work.

Label: S125 · Distributed training

Desc: Model parallelism partitions model computation or parameters across devices. All-reduce combines corresponding tensor entries across participants and returns the reduction to each; it is not matrix multiplication.

(e) Recall + stretch · 10 min

Self-check (answer follows each question).

Recall — click each to reveal

Q: In one sentence, what does a matrix ‘do’ to a vector?

A: It applies a linear transformation — it maps the input vector to a new vector by stretching, rotating, shearing, and/or projecting. Lines map to lines or points; the origin stays put.

Q: What are the columns of an m×n matrix A, geometrically?

A: Each column of A is where a basis vector of the input space lands after A is applied. Column j = A · eⱼ, where eⱼ is the j-th standard basis vector.

Q: State the shape rule for matrix multiplication.

A: (m × k) · (k × n) = (m × n). Inner dimensions must match; the outer dimensions form the result shape.

Q: Why is A · B ≠ B · A in general?

A: Because matrix multiplication is function composition. Applying A after B is not the same as applying B after A — rotating and then stretching gives a different final shape than stretching and then rotating.

Q: What is the identity matrix and what does it do?

A: Iₙ has 1s on the diagonal and 0s elsewhere. For m×n A, A · Iₙ = Iₘ · A = A. Use the identity with the correct dimension on each side.

Q: When does A⁻¹ exist?

A: Only when A is square (n × n) AND det(A) ≠ 0 (equivalently: A has full rank; equivalently: all eigenvalues are non-zero; equivalently: A · x = 0 only when x = 0). Rectangular matrices don't have true inverses (use pseudoinverse).

Q: You need to solve A x = b in code. What's the wrong way and the right way?

A: Prefer np.linalg.solve for a nonsingular square system and np.linalg.lstsq for least squares. Neither eliminates sensitivity inherent in an ill-conditioned problem; avoid explicitly forming an inverse unless the inverse itself is needed.

Q: What is an eigenvector of A? What is the eigenvalue?

A: An eigenvector v is nonzero and satisfies A · v = λ · v. Its line maps into itself; a negative real λ reverses orientation and λ=0 collapses it. The scalar λ is its eigenvalue.

Q: Why does PCA compute the eigenvectors of the covariance matrix?

A: The leading unit eigenvector maximizes projected variance; each later one does so subject to orthogonality to the earlier choices. Eigenvalues are the variances along those axes. Projection onto the top k minimizes squared reconstruction error among k-dimensional linear subspaces of centered data. This objective is not a guarantee of downstream predictive quality.

Q: Stretch — you have a 10M × 10k user × item ratings matrix. Compute a rank-50 approximation. What's the operation, what does it cost, and why?

Stretch: true

A: A rank-k approximation uses truncated SVD or an iterative/randomized method. Cost depends on nonzeros, matrix products, iterations and requested accuracy. Sparse solvers do not universally cost O(mnk), and matrix completion from missing ratings is not simply SVD of a zero-filled matrix.

Explain-out-loud test

1. What does a matrix do to a vector, in plain English?
2. Why is A @ B not the same as B @ A?
3. What is an eigenvector, and why does that word appear in every PCA paper?

What comes next

Hub: The 6-Month Learning Plan

Part of a 130-session evergreen learning series. Session structure: intuition → visual → hands-on → production war stories → recall. Duration: 90 minutes.

Corrective source notes

• numpy.linalg.lstsq — NumPy v2.5 Manual — numpy.linalg.lstsq minimizes the residual Euclidean norm and returns a minimum-norm solution if minimizers are nonunique. Checked 2026-10-07.
• numpy.linalg.cond — NumPy v2.5 Manual — The condition number is a numerical sensitivity measure distinct from raw determinant magnitude. Checked 2026-10-07.

These citations support the named claims, not blanket certification of the lesson. Worked proofs are derived in the text; execution scope and finite-test limits are recorded separately.