# 🔍 Gradient Checking with a Finite-Difference Harness

By Better Wise

> Sign in to view the full content of this lesson.

## A reference answer for the gradient

Last lesson you kept the normal-equation solution as the arbiter: if gradient descent misbehaved, you had an exact $w^\star$ to compare against. This lesson builds the same safety net for the gradient itself. **You do not get to run an optimizer until the derivative passes the harness.**

By the end of this lesson you will be able to:

- Derive the central-difference estimate of a gradient and explain why it beats the forward difference.
- Write one reusable `check_grad` function that takes a loss, an analytic gradient and a parameter array.
- Pick `eps` and a relative tolerance from the dtype, not from vibes.
- Read a failure pattern (all coordinates vs one coordinate) and name the bug class it implies.

> **Note:** 
> Everything below assumes `float64`. Cast your arrays before you check, or the harness will report roundoff noise as a gradient bug.

## The finite-difference harness, built once

The definition of a partial derivative is a limit you can approximate. Central difference, one coordinate at a time:

$$
\frac{\partial f}{\partial x_i} \approx \frac{f(x + \varepsilon e_i) - f(x - \varepsilon e_i)}{2\varepsilon}
$$

Truncation error is $O(\varepsilon^2)$ instead of the $O(\varepsilon)$ of the forward difference, so you can use a comfortably large $\varepsilon$ and still get 8-9 correct digits.

```python
import numpy as np

def numeric_grad(f, x, eps=1e-5):
    """Central-difference gradient of scalar f at x. Restores x on exit."""
    x = np.asarray(x, dtype=np.float64)
    g = np.zeros_like(x)
    it = np.nditer(x, flags=["multi_index"])
    for _ in it:
        i = it.multi_index
        orig = x[i]
        x[i] = orig + eps; fp = f(x)
        x[i] = orig - eps; fm = f(x)
        x[i] = orig                      # <-- never skip this
        g[i] = (fp - fm) / (2 * eps)
    return g

def check_grad(f, grad_fn, x, eps=1e-5, tol=1e-6):
    x = np.array(x, dtype=np.float64)     # work on a copy
    num = numeric_grad(f, x, eps)
    ana = np.asarray(grad_fn(x), dtype=np.float64)
    assert ana.shape == num.shape, (ana.shape, num.shape)
    rel = np.abs(num - ana) / np.maximum(1e-8, np.abs(num) + np.abs(ana))
    idx = np.unravel_index(rel.argmax(), rel.shape)
    return bool(rel.max() < tol), float(rel.max()), idx
```

The denominator `|num| + |ana|` is what makes this a *relative* test. Absolute differences are useless here: a 1e-3 error is fatal on a gradient entry of 1e-2 and irrelevant on one of 1e4.

Now the worked example, ridge regression, the same model as lesson 1:

```python
rng = np.random.default_rng(0)
n, d, lam = 50, 3, 0.1
X, y = rng.normal(size=(n, d)), rng.normal(size=n)

def loss(w):
    r = X @ w - y
    return 0.5 * (r @ r) / n + 0.5 * lam * (w @ w)

def grad(w):
    return X.T @ (X @ w - y) / n + lam * w

print(check_grad(loss, grad, rng.normal(size=d)))
# (True, 1.7e-09, (1,))
```

Worst relative error 1.7e-09: the gradient is correct to more digits than you will ever need. Now inject the classic bug, dropping the `1/n` on the data term only:

```python
def grad_bug(w):
    return X.T @ (X @ w - y) + lam * w
```

This reports roughly `0.96`, and here is the diagnostic value: **every coordinate fails by the same factor**, because the analytic gradient is uniformly inflated by `n`. A uniform failure is a scaling bug. A single failing coordinate is an indexing bug in that slice.

The central difference estimate is _____.
Its truncation error is _____, which is why eps around 1e-5 is safe in float64.
Gradient comparison uses the _____ error |num - ana| / max(1e-8, |num| + |ana|).
A worst relative error below _____ is a pass for a float64 central-difference check.
When every coordinate fails by the same factor, suspect a _____ bug.

| dtype | `eps` | pass if worst rel err |
| --- | --- | --- |
| float64 | 1e-5 | < 1e-6 |
| float32 | 1e-3 | < 1e-2 |

## Exercises: break it on purpose, then fix it

**1. Predict, then run.** Before executing `grad_bug`, work out the expected relative error in terms of `n` and `lam`. Does it depend on which coordinate you look at?

*Answer:* For small `lam` the data term dominates, so `ana ≈ n · num` and rel err `≈ (n-1)/(n+1) ≈ 0.96` at `n = 50`. It is the same for every coordinate, which is exactly why the failure *pattern* is a diagnostic signal.

**2. Sweep `eps`.** Run `check_grad` with `eps` in `[1e-2, 1e-4, 1e-6, 1e-8, 1e-10]` on the correct `grad` and plot worst relative error against `eps`.

*Answer:* You get a V shape. Large `eps` is dominated by truncation error growing like `eps^2`; tiny `eps` is dominated by roundoff growing like `eps_mach / eps`. The minimum sits near `eps_mach^(1/3) ≈ 6e-6`, which is why 1e-5 is the default almost everywhere.

**3. Project thread: check, then train.** Once `check_grad` passes, run gradient descent on `loss` and compare the final objective against the closed-form ridge solution $(X^\top X + n\lambda I)^{-1} X^\top y$. If GD converges to a different value, where is the bug?

*Answer:* Gradient checking only validates the derivative. If the gradient passes but training lands somewhere else, the fault is in the update rule (`w -= lr * grad`, not `+`), the learning rate, or the loss you are actually minimizing (for example, forgetting `lam` in `loss` but keeping it in `grad`, which the harness would have caught).

For a model with several parameter arrays, flatten them into one vector with a small `pack`/`unpack` pair and reuse `check_grad` unchanged. You will need that in Module 3.

## Where gradient checks lie to you

- **Checking at a kink.** ReLU at exactly 0, L1 penalties, `log` of a value near 0, or any clipping boundary. The one-sided derivatives disagree and the central difference averages them. Pick a random point strictly inside the domain.
- **Forgetting to restore the perturbed entry.** The next coordinate's evaluation inherits the leftover perturbation and everything after it is garbage silently.
- **Comparing sums instead of entries.** `abs(num.sum() - ana.sum())` can hide two large opposite-sign errors perfectly. Check elementwise, then report the max.
- **Mismatched regularization.** If `lam` appears in `grad` but not in `loss`, the harness will catch it; if `lam = 0`, it will not test it at all. Use a nonzero `lam`.
- **Checking in float32.** Roundoff is around 1e-7 relative, so the difference quotient drowns. Either cast to float64 or use the float32 row of the table above.

## Rules of thumb before you trust a gradient

1. **Smallest possible input.** Check on a 3-4 sample, 2-3 feature toy problem first. A bug that shows up at `n = 4` is one you can read in the debugger; the same bug at `n = 10000` is one you will bisect for an hour.
2. **Random points, not zeros.** All-zero parameters switch off whole branches (any multiplicative path) and can make a broken gradient look correct. Sample from a normal and run 2-3 seeds.
3. **1e-6 is pass, 1e-4 is suspicious, 1e-2 is broken.** ReLU networks land around 1e-7 to 1e-5 because of kinks; anything above 1e-4 means keep looking.
4. **Check the regularizer separately.** Set `lam = 0`, confirm the data term, then set `X = 0` and confirm only the penalty remains. Two clean checks beat one muddy one.
5. **Re-run the harness after every refactor.** It costs a second and it is the cheapest test in the project.
6. **Timing tells you it is working.** A full check costs two forward passes per parameter. Slow but honest; never run it inside a training loop.

## What to carry forward

You now have one function that converts "I think the gradient is right" into a pass/fail with a number attached. Every model in this course gets run through it before it gets trained: ridge here, cross-entropy next, then PCA and k-means objectives, then the MLP.

Keep the V-shaped `eps` sweep in mind too. It is the first place to look any time a check that used to pass starts drifting.

Next up: **Logistic Regression to Softmax: Cross-Entropy Done Right**, where the loss gets a `log` and the gradient gets a `sigmoid`, which is precisely the combination that produces both silent numerical overflow and gradients that pass the harness while still being wrong.

Front: Central-difference formula for one coordinate of the gradient
Back: (f(x + eps e_i) - f(x - eps e_i)) / (2 eps)

Front: Why use central difference over forward difference?
Back: Truncation error is O(eps^2) instead of O(eps), so a larger eps still gives many correct digits.

Front: Default eps for a float64 gradient check
Back: About 1e-5, close to the cube root of machine epsilon (6e-6).

Front: Relative error formula for comparing two gradients
Back: |num - ana| / max(1e-8, |num| + |ana|), computed elementwise.

Front: Passing tolerance for a float64 central-difference check
Back: Worst relative error below 1e-6; above 1e-2 means a genuine bug.

Front: Why not use absolute differences to compare gradients?
Back: Gradient entries vary in scale, so absolute error hides faults in small entries and exaggerates large ones.

Front: Float32 settings for the harness
Back: eps around 1e-3 with tolerance around 1e-2, or cast to float64 first.

Front: Cost of a full finite-difference check
Back: Two loss evaluations per parameter, so O(P) forward passes.

Front: Where must you never run a gradient check?
Back: At kinks or domain boundaries such as ReLU at 0, L1 at 0, or log near 0.

Front: Every coordinate fails by the same factor. What kind of bug?
Back: A scale or normalization bug, such as a missing 1/n on the whole gradient.

Front: A single coordinate fails. What kind of bug?
Back: An indexing or slicing error affecting that entry only.

Front: Why must the harness restore the perturbed entry?
Back: Later coordinates reuse the same array, so leftover perturbations silently corrupt the rest of the check.

Front: What does a passing gradient check not prove?
Back: That the update rule, learning rate, or loss formulation is correct; it validates only the derivative.

---

## Contents

Part of [🧮 Mathematics for Machine Learning 4: Build Core ML Algorithms from Scratch](https://everwise.study/courses/mathematics-for-machine-learning-4-build-core-ml-algorithms-from-scratch.md)

- **📈 Linear and Cross-Entropy Models from Scratch**
  - [📉 Linear Regression with Normal Equations and Gradient Descent](https://everwise.study/courses/mathematics-for-machine-learning-4-build-core-ml-algorithms-from-scratch/linear-regression-normal-equations-and-gradient-descent.md)
  - [🔍 Gradient Checking with a Finite-Difference Harness](https://everwise.study/courses/mathematics-for-machine-learning-4-build-core-ml-algorithms-from-scratch/gradient-checking-finite-difference-harness.md) (sign-in required) — this lesson
  - [🎚️ Logistic Regression to Softmax: Cross-Entropy Done Right](https://everwise.study/courses/mathematics-for-machine-learning-4-build-core-ml-algorithms-from-scratch/logistic-regression-to-softmax-cross-entropy.md)
  - [✅ Checkpoint: Linear and Cross-Entropy Models](https://everwise.study/courses/mathematics-for-machine-learning-4-build-core-ml-algorithms-from-scratch/checkpoint-linear-and-cross-entropy-models.md)
- **🔵 Unsupervised Learning from Scratch**
  - [🧭 PCA via SVD: Explained Variance and Component Selection](https://everwise.study/courses/mathematics-for-machine-learning-4-build-core-ml-algorithms-from-scratch/pca-via-svd-explained-variance.md)
  - [🌀 k-means: Lloyd's Algorithm and k-means++ Initialization](https://everwise.study/courses/mathematics-for-machine-learning-4-build-core-ml-algorithms-from-scratch/k-means-lloyd-and-k-means-plus-plus.md)
  - [🔔 Gaussian Mixture Models with Expectation Maximization](https://everwise.study/courses/mathematics-for-machine-learning-4-build-core-ml-algorithms-from-scratch/gaussian-mixture-models-expectation-maximization.md)
  - [✅ Checkpoint: PCA, k-means and GMM](https://everwise.study/courses/mathematics-for-machine-learning-4-build-core-ml-algorithms-from-scratch/checkpoint-pca-kmeans-and-gmm.md)
- **🧠 Backpropagation and Benchmark Capstone**
  - [🔗 Backpropagation for a One-Hidden-Layer MLP](https://everwise.study/courses/mathematics-for-machine-learning-4-build-core-ml-algorithms-from-scratch/backpropagation-for-one-hidden-layer-mlp.md)
  - [📊 Cross-Validation and a Fair Baseline Comparison Protocol](https://everwise.study/courses/mathematics-for-machine-learning-4-build-core-ml-algorithms-from-scratch/cross-validation-and-baseline-comparison-protocol.md)
  - [🏁 Capstone: NumPy Implementations vs Library Baselines](https://everwise.study/courses/mathematics-for-machine-learning-4-build-core-ml-algorithms-from-scratch/capstone-numpy-vs-library-baselines.md) (premium content)
  - [✅ Checkpoint: Backpropagation and Model Comparison](https://everwise.study/courses/mathematics-for-machine-learning-4-build-core-ml-algorithms-from-scratch/checkpoint-backpropagation-and-model-comparison.md) (premium content)
- [Course cheatsheet](https://everwise.study/courses/mathematics-for-machine-learning-4-build-core-ml-algorithms-from-scratch/cheatsheet.md) (premium content)