Linear Algebra for ML

Solution: (x, y, z) = (1, 2, 3). Back-substitution starts from the bottom row (simplest equation) and works up, plugging known values into earlier equations.

Row 2 is exactly 2 × Row 1 — it carries no new information. R2 ← R2 − 2R1 gives [0 0 | 0]. The system reduces to x + y = 3, x − y = 1, which has the unique solution x = 2, y = 1. The "extra" equation was linearly dependent — its rank contribution was zero. In ML, this is like having a duplicate feature.

Verify: A · A⁻¹ = [[2×3+1×(−5), 2×(−1)+1×2],[5×3+3×(−5), 5×(−1)+3×2]] = [[1,0],[0,1]]. The fact that det = 1 is special — integer matrices with det = ±1 have integer inverses.

Line 5 is the bug. The inverse formula for x0 should use A[1][1]*b[0] − A[0][1]*b[1], not +. Recall that A⁻¹ = (1/det)[[d, −b],[−c, a]]. The off-diagonal entries are negated. With the + sign, you're multiplying by [[d, b],[−c, a]] instead — wrong sign on the (1,2) entry.

[[2,4],[3,6]] has det = 0. Notice that row 2 = 1.5 × row 1 — the rows are linearly dependent. The matrix maps all of 2D space onto a line. You can't invert that — information is destroyed. In ML, this is like having one feature that's a perfect linear function of another.

Check: λ₁ + λ₂ = 7 = tr(A). λ₁ × λ₂ = 10 = det(A). Both match.

So v = [1, 1] (or any scalar multiple). The eigenvector direction is the 45-degree line. A stretches vectors along this line by a factor of 5.

Av = λv confirmed. The matrix A applied to [1,1] produces 5×[1,1] — pure scaling, no rotation. That's exactly what "eigenvector" means.

So v = [1, −2] (or any scalar multiple). This eigenvector points in a different direction from the λ=5 eigenvector. Together, these two eigenvectors form a basis for R² — any vector can be decomposed into these two directions, each scaled by its eigenvalue.

javascript
function eigenvalues2x2(A) {
  const tr = A[0][0] + A[1][1];
  const det = A[0][0]*A[1][1] - A[0][1]*A[1][0];
  const disc = Math.sqrt(tr*tr - 4*det);
  return [(tr + disc)/2, (tr - disc)/2];
}

The Spectral Theorem guarantees that real symmetric matrices have: (1) all real eigenvalues, and (2) orthogonal eigenvectors. This is why covariance matrices (which are symmetric) always have real eigenvalues and orthogonal principal components. The eigenvalues can be negative (just not for PSD matrices) and the matrix can be singular (if an eigenvalue is 0).

For a diagonal matrix, the singular values are just the absolute values of the diagonal entries. U = V = I (identity), Σ = A itself.

Since A is symmetric with eigenvalues 3 and 1, the singular values are |3| = 3 and |1| = 1. For symmetric PSD matrices, singular values equal eigenvalues — the two decompositions coincide.

That's 6003 vs 1,000,000 — a 166× compression. And by the Eckart-Young theorem, this is the best possible rank-3 approximation in both Frobenius and spectral norm. The approximation error is σ₄ = 1 (in spectral norm).

Rank-20 needs 512×20 + 20 + 512×20 = 20,500 numbers vs the original 512×512 = 262,144. Ratio: 262,144/20,500 ≈ 12.8, so roughly 13:1. This is the principle behind truncated SVD for image compression and dimensionality reduction.

The rank of a matrix equals the number of non-zero singular values. Three non-zero singular values (10, 5, 2) means rank 3. The matrix maps 5D space into a 3D subspace — two dimensions of information are lost. In ML, this means only 3 of the 5 features are linearly independent.

C₁₁ = 4 is the variance of feature 1. C₁₂ = C₂₁ = 2 is the covariance between features — they're positively correlated. The covariance matrix is always symmetric.

The first principal component captures variance 6, the second captures 2. The eigenvector for λ=6 is [1,1]/√2 (the 45-degree direction). Data varies most along this diagonal — projecting onto it keeps 75% of the total variance.

75% of the variance is captured by a single direction. In practice, you pick k such that the cumulative explained variance exceeds a threshold (commonly 95%). With just one component here, you'd lose 25% of the information.

PCA finds the direction of maximum variance, which is generally a linear combination of the original features. With covariance 35 between features 1 and 2, a rotated axis capturing both features jointly explains more variance than either feature alone. Feature selection picks individual features; PCA finds optimal linear combinations. They're fundamentally different approaches.

The Jacobian of a linear transformation y = Wx is simply the matrix W itself. This makes sense: y₁ = 2x₁ + x₂, so ∂y₁/∂x₁ = 2. This is why backpropagation through a linear layer just multiplies by W^T.

Quadratic forms are scalar-valued: they map a vector to a single number. The loss function in ridge regression is x^TAx + bx, a quadratic form plus a linear term.

Since A is symmetric (A = A^T), A + A^T = 2A. So the gradient simplifies to 2Ax = 2[4, 7] = [8, 14]. Symmetric A is the common case in ML (covariance matrices, Hessians).

Negative gradient means x should increase to reduce the loss. The residual Ax − b = [0, −1, −1] tells us the predictions are too low for data points 2 and 3. Gradient descent would step x in the direction [2, 4] (negating the gradient).

The chain unfolds as: ∂L/∂W₁ = ∂L/∂y · ∂y/∂h₂ · ∂h₂/∂h₁ · ∂h₁/∂z₁ · ∂z₁/∂W₁. That's 5 Jacobian matrices multiplied together. Each layer contributes two factors: one for the linear part (W), one for the activation (σ). The loss contributes one more. This is why deep networks are prone to vanishing/exploding gradients — you're multiplying many matrices together.

Line 7 is the bug. The gradient is 2A^T(Ax − b), not 2A(Ax − b). The code uses matmul(A, r) instead of matmul(transpose(A), r). Missing the transpose is one of the most common bugs in implementing gradient descent for linear models. Dimensionally: if A is [m,n] and r is [m,1], then A·r would be an invalid shape — only A^T·r = [n,m]×[m,1] = [n,1] has the right shape to update x.

(A^TA)₁₁ = 1² + 2² + 3² = 14. This is the sum of squares of the x-coordinates. The full A^TA is a 2×2 symmetric matrix encoding the structure of the feature data.

The best-fit line is y = 1.5x + 0.33. At x=1: y=1.83 (error 0.17), x=2: y=3.33 (error 0.33), x=3: y=4.83 (error −0.17). The residuals sum to approximately zero, as they must for any least squares solution with an intercept term.

Since x₂ = 1.8x₁ + 32, column 2 of A is a linear combination of column 1 and the intercept column. A^TA is singular (det = 0), and the inverse doesn't exist. There are infinitely many solutions: you can trade off weight on x₁ vs x₂ arbitrarily. The SVD pseudoinverse gives the minimum-norm solution, and ridge regression picks a unique solution by penalizing large weights.

javascript
function normalEquation(A, b) {
  const n = A.length;
  // A^T A (2x2)
  let a00=0,a01=0,a11=0, b0=0,b1=0;
  for (let i=0; i<n; i++) {
    a00 += A[i][0]*A[i][0];
    a01 += A[i][0]*A[i][1];
    a11 += A[i][1]*A[i][1];
    b0  += A[i][0]*b[i];
    b1  += A[i][1]*b[i];
  }
  const det = a00*a11 - a01*a01;
  return [(a11*b0 - a01*b1)/det, (-a01*b0 + a00*b1)/det];
}

The L∞ norm only cares about the single largest entry. It's used in adversarial robustness: an L∞ adversarial perturbation can change every pixel, as long as no single pixel changes by more than ε.

The Frobenius norm treats the matrix as a long vector and takes its L2 norm. It equals √(∑ σ_i²) — the square root of the sum of squared singular values. Weight decay in neural networks penalizes the squared Frobenius norm of the weight matrices.

The nuclear norm is the "L1 norm of the singular values." Just as L1 on vectors promotes sparsity, the nuclear norm on matrices promotes low rank. Minimizing ||W||_* pushes small singular values to exactly zero, effectively reducing the matrix rank. Used in matrix completion (recommender systems) and low-rank matrix recovery.

Always: ||v||_∞ ≤ ||v||₂ ≤ ||v||₁. The max absolute entry is at most the Euclidean length (since the other terms are non-negative), and the Euclidean length is at most the sum of absolute values (by Cauchy-Schwarz). Verify with v = [3, −4, 0, 5]: ||v||_∞ = 5 ≤ ||v||₂ = 7.07 ≤ ||v||₁ = 12.

Both eigenvalues are positive, so A is positive definite. This means x^TAx > 0 for all non-zero x, and the quadratic form defines a bowl-shaped surface with a unique minimum at x = 0.

Verify: LL^T = [[2,0],[1,√2]] × [[2,1],[0,√2]] = [[4,2],[2,3]] = A. Cholesky is 2× faster than LU decomposition for solving Ax = b when A is PD, and it's the standard method for sampling from multivariate Gaussians: x = L · z + μ where z ~ N(0, I).

tr(B) = 2, det(B) = 1 − 4 = −3. Eigenvalues: λ² − 2λ − 3 = 0 → λ = 3, −1. Since −1 < 0, B is not positive definite. In fact, take x = [1, −1]: x^TBx = 1 − 2 − 2 + 1 = −2 < 0. Positive diagonal entries are necessary but not sufficient for PD — the off-diagonal correlations can make the matrix indefinite.

Adding εI shifts every eigenvalue by +ε. If ε = 0.01, the problematic −0.001 becomes 0.009 > 0, making the matrix PD. This is called jitter or ridge regularization and is standard practice in Gaussian processes, numerical Cholesky, and any code that inverts near-singular matrices. In GPyTorch, JAX, and TensorFlow, you'll see + 1e-6 * I everywhere.

Forward substitution through L is O(n²), much cheaper than computing A⁻¹ which is O(n³). The Cholesky factorization itself is O(n³/3), but once you have L, each new right-hand side b costs only O(n²).

For a matmul [m,k]×[k,n], FLOPs = 2mkn. Here m=n=50, k=1000. On a modern CPU doing ~10 GFLOPs, this takes about 0.5 ms. Trivial. But if d=50000 instead of 50, the cost jumps to 5 trillion FLOPs — 500 seconds on the same CPU.

30× less than the covariance matrix computation. The eigendecomposition is cheap because the covariance matrix is only [50,50]. But the cost grows as d³ — at d=5000 it would be 167 billion FLOPs, 33,000× more expensive.

Cheaper than the covariance computation because we only project onto 10 directions, not all 50. The output Z [1000, 10] has 80% fewer features, which makes subsequent training much faster.

Training cost is dominated by Z^TZ at 0.2M FLOPs — 25× less than the covariance step. PCA reduced d from 50 to 10, which made training 25× cheaper (cost scales as d²).

The covariance computation is 78% of the total cost. This is the bottleneck. At a modern GPU's ~10 TFLOPS, the entire pipeline takes under a microsecond. But scale n to 1 million and d to 5000, and the covariance step alone costs 50 trillion FLOPs — 5 seconds on GPU. Now you see why efficient covariance estimation (randomized methods, streaming PCA) matters at scale.

Center data → Covariance matrix → Eigendecompose → Project onto top-k → Normal equations

You must center before computing covariance (uncentered data gives a biased covariance). You must eigendecompose before projecting (you need the eigenvectors). And you train on the projected data, not the original.

The dominant term is the covariance: O(d²n). Doubling both: (2d)² × 2n = 4d² × 2n = 8d²n. That's 8× the cost. The eigendecomposition is O(d³): (2d)³ = 8d³, also 8×. For large d, eigendecomposition eventually dominates (d³ vs d²n when d > n).

Tiny. But at n=1M and d=50K: X alone is 200 GB. The covariance matrix C [50K, 50K] is 10 GB. This is why large-scale PCA uses randomized algorithms that never materialize the full covariance matrix.