Numerical Methods

An 11% error from a single subtraction! Both operands were stored with ~15.9 digits of precision, but the subtraction wiped out all but ~1 digit. The 0.110223... comes from 2^-53 × 2⁵⁰ = 2^-3 ≈ 0.125, a rounding artifact of binary representation.

The answer is dominated by the max term (30). Without the trick, you'd compute exp(30) ≈ 10¹³ which is fine, but exp(1000) would overflow. The trick ensures the largest exponent is always 0, so no overflow is possible regardless of input magnitude. This is the backbone of stable cross-entropy loss.

javascript
function logSumExp(arr) {
  const maxVal = Math.max(...arr);
  const sumExp = arr.reduce((s, x) =>
    s + Math.exp(x - maxVal), 0);
  return maxVal + Math.log(sumExp);
}

Line 8 uses the textbook formula Var = E[X²] - E[X]². This is mathematically correct but numerically catastrophic. When the mean is large (e.g., arr = [1000000.0, 1000000.1]), sumSq/n ≈ 10¹² and (sum/n)² ≈ 10¹². Subtracting two huge nearly-equal numbers → catastrophic cancellation → the result can go negative.

The fix: use the two-pass algorithm. First compute the mean, then compute ∑(x_i - mean)² / n. The subtraction happens at each element (small numbers), not at the aggregate (huge numbers).

We started with 7 digits of precision in each operand and lost ~6 of them, leaving ~1 reliable digit. This formula is a quick mental check: if a/|a-b| ≈ 10^k, you lose k digits. In FP32 (7 digits), losing 6 leaves 1 — barely usable. In FP64 (15.9 digits), losing 6 leaves ~10 — usually fine.

A condition number of 40,000 means a perturbation of 10^-7 in b (FP32 machine epsilon) gets amplified to 10^-7 × 4×10⁴ = 4×10^-3 relative error in x. In FP32, you lose ~4.6 of your 7 digits. This matrix is ill-conditioned because its rows are nearly parallel.

A 0.005% change in b produced a 100% change in x. The solution completely flipped! This is the condition number at work: κ ≈ 40,000 means small perturbations get amplified by up to 40,000×.

The identity matrix is perfectly conditioned. Solving Ix = b gives x = b exactly (no amplification). κ = 1 is the best possible. Orthogonal matrices (Q^TQ = I) also have κ = 1. This is why QR decomposition is preferred over direct inversion — it factors A into a product involving orthogonal matrices.

javascript
function conditionNumber(a, b, c, d) {
  const normA = Math.max(
    Math.abs(a) + Math.abs(b),
    Math.abs(c) + Math.abs(d)
  );
  const det = a * d - b * c;
  const normAinv = Math.max(
    Math.abs(d/det) + Math.abs(b/det),
    Math.abs(c/det) + Math.abs(a/det)
  );
  return normA * normAinv;
}

Jacobi formula: x₁^(k+1) = (9 - x₂^(k))/4,   x₂^(k+1) = (11 - x₁^(k))/3

The true solution is x = [1, 3]. After just 2 iterations we have [1.333, 2.917] — already close. Convergence is guaranteed here because A is strictly diagonally dominant (|4| > |1| and |3| > |1|).

After 3 iterations: [1.021, 3.222]. The error in x₁ went from 1.333 to 0.021 — it's converging nicely. True solution: [1, 3].

The spectral radius tells you the convergence rate. Each iteration reduces the error by a factor of ~0.289. After k iterations, error ≈ 0.289^k × initial error. To get 6 digits of accuracy: 0.289^k < 10^-6 → k > 11.1. So about 12 iterations suffice.

javascript
function jacobiStep(A, b, x) {
  const n = b.length;
  const xNew = new Array(n);
  for (let i = 0; i < n; i++) {
    let sum = 0;
    for (let j = 0; j < n; j++) {
      if (j !== i) sum += A[i][j] * x[j];
    }
    xNew[i] = (b[i] - sum) / A[i][i];
  }
  return xNew;
}

The true value is √2 = 1.41421356237... By x₃ we have 6 correct digits, and x₄ gives 12+. Watch the quadratic convergence: error went 0.414 → 0.086 → 0.0025 → 0.000002 → 0.000000000004. Each step roughly squares the error.

The constant C ≈ 0.35 = |f''(root)|/(2|f'(root)|) = |2|/(2×2√2) = 1/(2√2) ≈ 0.354. This is the convergence constant for Newton's method on x²-2.

javascript
function newton(f, fPrime, x0, nIter) {
  let x = x0;
  for (let i = 0; i < nIter; i++) {
    x = x - f(x) / fPrime(x);
  }
  return x;
}

Line 6 divides by fp(x) without checking if the derivative is zero. If f'(x) = 0 at any iteration (a horizontal tangent), the step becomes ±Infinity, and all subsequent iterates are NaN. The fix: check if (Math.abs(fp(x)) < 1e-12) return x; before dividing — or throw an error indicating the method has stalled at a stationary point.

Equivalently: x₅ = 0.8⁵ × 5 = 0.32768 × 5 = 1.6384. After 5 steps with α=0.1, the loss went from f(5)=25 down to f(1.6384)=2.684 — about a 10× reduction. To reach f < 0.01, you need 0.8^2k × 25 < 0.01, giving k ≥ 35 steps.

For general f with Lipschitz-continuous gradient with constant L, the stability bound is α < 2/L. For f(x) = x², L = 2, so α < 2/2 = 1. This is why learning rate schedules start small — to ensure stability.

javascript
function gradientDescent(grad, x0, lr, nSteps) {
  let x = x0;
  for (let i = 0; i < nSteps; i++) {
    x = x - lr * grad(x);
  }
  return x;
}

31 steps to get 3 digits of accuracy. With α=0.3, the factor is 0.4, and you'd need only ln(0.001)/ln(0.4) = 7.5 → 8 steps. Larger learning rate = faster convergence (until you hit the stability limit).

Error = |0.34375 - 0.33333| = 0.01042. With n=4, the error is O(h²) = O(0.0625). The ratio error/h² ≈ 0.01042/0.0625 = 0.167 ≈ 1/6, which matches the error formula -(b-a)h²f''(c)/12 = -1×0.0625×2/12 = -0.01042. Theory confirmed!

Simpson's rule gives the exact answer for x²! This is because Simpson's rule exactly integrates polynomials up to degree 3 (its error involves the 4th derivative, and f''''(x²) = 0). This is a general pattern: higher-order quadrature rules are exact for higher-degree polynomials.

javascript
function trapezoidal(f, a, b, n) {
  const h = (b - a) / n;
  let sum = f(a) + f(b);
  for (let i = 1; i < n; i++) {
    sum += 2 * f(a + i * h);
  }
  return (h / 2) * sum;
}

In 1D, trapezoidal needs ~289 points while Monte Carlo needs a trillion. Deterministic quadrature wins overwhelmingly in low dimensions. But in d dimensions, trapezoidal needs n^d points (curse of dimensionality), while Monte Carlo still needs just 1/ε², regardless of d. This crossover at d ≈ 4-8 is why stochastic methods dominate in ML.

The error is 0.0601 ≈ 3×2×0.01 = 0.06 = f''(2)×h/2 (the leading error term for forward difference is f''(c)h/2). This is O(h) — first-order accurate.

Central difference error: 0.0001. Forward difference error: 0.0601. That's 601× more accurate for the same h! The central difference is O(h²), so its error is ~h²×f'''(2)/6 = 0.0001×6/6 = 0.0001. This is why gradient checking always uses central differences.

(We got 601 empirically due to higher-order terms.) The ratio scales as O(1/h). Smaller h → bigger advantage for central differences. At h = 10^-5 (typical for gradient checking), central differences are ~10⁵× more accurate than forward differences. One extra function evaluation buys five orders of magnitude.

javascript
function gradCheck(f, x, analyticGrad, h) {
  const numGrad = (f(x + h) - f(x - h)) / (2 * h);
  const denom = Math.max(
    Math.abs(numGrad), Math.abs(analyticGrad), 1e-8
  );
  return Math.abs(numGrad - analyticGrad) / denom;
}

Error = |2.721 - 2.718| = 0.003. The second-derivative formula has error O(h²), so with h=0.1, error ≈ h² f''''(1)/12 = 0.01×e/12 ≈ 0.0023. Close to our observed error. This formula is used in Hessian-vector products for second-order optimization.

Using CSR, we only touch 4 nonzero entries instead of all 16. For each row i, we iterate from row_ptr[i] to row_ptr[i+1], multiply values[j] by x[col_idx[j]], and accumulate. Row 3: j goes from 2 to 4, giving values[2]×x[1] + values[3]×x[3] = 4×2 + 1×4 = 12.

Hmm, 0.20 with 90% sparsity. At 99% sparsity (10,000 nonzeros), the ratio drops to ~0.02 (50× savings). At 99.9% sparsity (typical for large graphs), it's ~0.002 (500× savings). The breakeven point (CSR uses same memory as dense) is around 50% sparsity — below that, just use dense.

javascript
function spMV(values, colIdx, rowPtr, x) {
  const nRows = rowPtr.length - 1;
  const y = new Array(nRows).fill(0);
  for (let i = 0; i < nRows; i++) {
    for (let j = rowPtr[i]; j < rowPtr[i + 1]; j++) {
      y[i] += values[j] * x[colIdx[j]];
    }
  }
  return y;
}

Even exp(11.1) ≈ 66,686 barely overflows FP16. Anything above ~11.09 in the exponent overflows. Real logits of 89-91 would produce Inf in FP16, giving NaN after softmax. This is why the stable softmax trick (subtract max) is not optional — it's a correctness requirement in low precision.

All intermediate values are in [0, 1] — perfectly safe for FP16. The max exponent is 0, so exp can never overflow. This is a concrete example of why the stable trick works: it shifts the dynamic range of the computation to a safe region.

This gradient (0.024) is fine for FP16. But imagine a 1000-class classifier where the true class gets softmax = 0.999 and most other classes get ~10^-6. Those gradients ARE below the FP16 subnormal threshold. They become zero — the model can't learn to adjust those logits. This is gradient underflow, and it's why loss scaling multiplies the loss by 1024 or more before backprop.

In practice, loss scales of 1024 or higher are common (AMP default starts at 65536 and dynamically adjusts). The optimizer divides by S after backprop to recover the true gradient. Dynamic loss scaling: if you see Inf/NaN in gradients, halve S; if no overflow for N steps, double S. This is how PyTorch's GradScaler works.

1. Reproduce with FP32 — Confirm the bug is precision-related, not a logic error. If FP32 also NaNs, the problem is elsewhere.

2. Locate first NaN layer — Insert hooks to check each layer's output for Inf/NaN. The first layer that produces NaN is the culprit.

3. Check value ranges — Log the max absolute activation/gradient at the NaN layer. If max > 65504 (FP16) or gradients < 6e-5, you've found the mechanism.

4. Apply fix — Clamp pre-softmax logits, enable loss scaling, switch to BF16, or add gradient clipping — depending on whether the issue is overflow or underflow.

5. Verify fix converges — Train for the full schedule and confirm loss matches FP32 baseline within acceptable tolerance.

After summing ~128 values (reaching sum ≈ 0.128), adding 0.001 rounds away to zero in BF16. The true sum should be 10.0, but BF16 gives ~0.128 — a 92% error! This is why gradient accumulators MUST be FP32, even when weights and activations are BF16. Every mixed-precision training recipe uses FP32 accumulators. This is the single most important numerical insight for ML training.