13 · Content revision 6

Check gradients before trusting backpropagation

Cross-check selected LLM-training derivatives from actual floating-point probes, with a local-smoothness requirement, scale-aware error, and deterministic tensor coordinates in Rust.

Predict one quadratic derivative

Start with a small scalar gradient check:

q(\theta)=\theta^2,\qquad \theta=3,\qquad h=0.1

The requested step asks for the two probes written as $2.9$ and $3.1$ . An f64 stores the nearby representable values $\theta_- = \operatorname{fl}(3-0.1)$ and $\theta_+ = \operatorname{fl}(3+0.1)$ instead. The implementation therefore derives the actual distances $h_-=3-\theta_-$ and $h_+=\theta_+-3$ rather than assuming that either distance equals the requested $0.1$ . In the usual rounded hand calculation, $q(2.9)=8.41$ and $q(3.1)=9.61$ , and the three-point estimate predicts the numerical derivative $6$ .

That result tests a candidate rather than creating one. The hand derivative of the function $\theta^2$ supplies candidate $6$ , which should pass. Candidate $5.5$ is also a perfectly finite number, but it should return a completed comparison with pass=false. A mathematical mismatch is different from an invalid step or a non-finite objective value.

Define one quadratic prediction and both analytic candidates rust/demos/ch13-gradient-checking/src/lib.rs#quadratic-gradient-prediction

pub fn quadratic(point: f64) -> f64 {
    point * point
}

/// Checks either the expected derivative `6` or a deliberately wrong candidate.
pub fn quadratic_gradient_check(analytic: f64) -> Result<ScalarGradientCheck, GradCheckError> {
    scalar_gradient_check(3.0, analytic, 0.1, STEP_TOLERANCE, quadratic)
}

Use the spacing between the values the computer can represent

The caller supplies a positive requested step $h$ , but floating-point addition and subtraction produce the actual representable probes

\theta_- = \operatorname{fl}(\theta-h),\qquad \theta_+ = \operatorname{fl}(\theta+h).

Their positive distances from the checked point are

h_- = \theta-\theta_-,\qquad h_+ = \theta_+-\theta.

The unequal-spacing three-point formula combines the two one-sided slopes:

f'(\theta)\approx\frac{h_+}{h_-+h_+}\frac{f(\theta)-f(\theta_-)}{h_-}+\frac{h_-}{h_-+h_+}\frac{f(\theta_+)-f(\theta)}{h_+}

Only when the actual distances are equal, $h_-=h_+=\widehat h$ , does the center coefficient cancel and the formula reduce to

f'(\theta)\approx\frac{f(\theta_+)-f(\theta_-)}{2\widehat h}.

The denominator then contains the actual distance $\widehat h$ , not merely the requested $h$ .

The implementation forms the weights without first adding the two possibly huge spacings. With $m=\max(h_-,h_+)$ , it uses $u_-=h_-/m$ and $u_+=h_+/m$ , then left weight $u_+/(u_-+u_+)$ and right weight $u_-/(u_-+u_+)$ . The ratios preserve the unequal-spacing coefficients while avoiding overflow in $h_-+h_+$ .

If $f$ is locally smooth enough—for example, it has a continuous third derivative across the interval from $\theta_-$ to $\theta_+$ —the three-point interpolation has truncation error of order $O(\max(h_-,h_+)^2)$ . Shrinking the requested step helps only until rounded probe locations and subtraction of nearly equal function values dominate.

The identity function $f(x)=x$ is a useful adversary. Choose a finite point and requested step whose two floating-point operations produce unequal nonzero distances. Each one-sided slope in the weighted formula is still $1$ , so the combined estimate must remain $1$ within tolerance. Dividing by the requested $2h$ instead can fail this test because that denominator does not describe the actual probes.

Local smoothness is a caller precondition, not something one passing comparison can prove. At $f(x)=|x|$ and $x=0$ , symmetric probes produce the numerical value $0$ , and candidate $0$ can pass even though the derivative does not exist at the corner. This is a false pass outside the method’s differentiable-point scope; the checker is not a universal detector for corners or other nonsmooth points.

compare_one_sided_slopes compares the recorded slopes with a scale-aware gap. Disagreement can indicate a kink, a step too large for local curvature, or rounding damage; agreement cannot establish differentiability.

The implementation requires finite $\theta$ , positive finite $h$ , a finite requested range, two finite probes that are strictly ordered around $\theta$ , positive finite actual distances, and finite function values at all three points. A collapsed probe is rejected rather than reported as a zero derivative.

Derive actual probe spacing before applying the unequal-spacing estimate rust/crates/llm-from-scratch/src/autograd/gradcheck.rs#central-difference

/// Approximates one derivative from the actual representable probe spacings.
///
/// The callback is evaluated exactly three times, in minus, center, plus order.
/// A meaningful gradient check requires a deterministic, side-effect-free
/// objective that is differentiable and sufficiently smooth throughout the
/// closed probe interval. Known or possible kinks must be excluded or examined
/// with [`compare_one_sided_slopes`]; a passing comparison is evidence for the
/// sampled smooth point, not a universal proof of a derivative.
pub fn central_difference(
    point: f64,
    step: f64,
    mut evaluate: impl FnMut(f64) -> f64,
) -> Result<CentralDifference, GradCheckError> {
    let geometry = perturbations(point, step)?;

    let minus_value = evaluate(geometry.minus_point);
    if !minus_value.is_finite() {
        return Err(GradCheckError::NonFiniteEvaluation {
            side: DifferenceSide::Minus,
            value: minus_value,
        });
    }

    let center_value = evaluate(point);
    if !center_value.is_finite() {
        return Err(GradCheckError::NonFiniteEvaluation {
            side: DifferenceSide::Center,
            value: center_value,
        });
    }

    let plus_value = evaluate(geometry.plus_point);
    if !plus_value.is_finite() {
        return Err(GradCheckError::NonFiniteEvaluation {
            side: DifferenceSide::Plus,
            value: plus_value,
        });
    }

    let left_slope = (center_value - minus_value) / geometry.minus_spacing;
    if !left_slope.is_finite() {
        return Err(GradCheckError::NonFiniteOneSidedSlope {
            side: DifferenceSide::Minus,
            value: left_slope,
        });
    }
    let right_slope = (plus_value - center_value) / geometry.plus_spacing;
    if !right_slope.is_finite() {
        return Err(GradCheckError::NonFiniteOneSidedSlope {
            side: DifferenceSide::Plus,
            value: right_slope,
        });
    }

    let derivative = geometry.left_weight * left_slope + geometry.right_weight * right_slope;
    if !derivative.is_finite() {
        return Err(GradCheckError::NonFiniteNumericalGradient { value: derivative });
    }

    Ok(CentralDifference {
        point,
        requested_step: step,
        minus_point: geometry.minus_point,
        plus_point: geometry.plus_point,
        minus_spacing: geometry.minus_spacing,
        plus_spacing: geometry.plus_spacing,
        minus_value,
        center_value,
        plus_value,
        left_slope,
        right_slope,
        left_weight: geometry.left_weight,
        right_weight: geometry.right_weight,
        stencil: geometry.stencil,
        derivative,
    })
}

Name the derivative-check quantities

Symbol	Operational meaning
$f$	The deterministic scalar loss-valued function being probed.
$\theta$	The finite scalar parameter or selected tensor coordinate.
$h$	The positive finite requested step used to form both floating-point probes.
$\theta_-,\theta_+$	The actual representable probes $\operatorname{fl}(\theta-h)$ and $\operatorname{fl}(\theta+h)$ .
$h_-,h_+$	The actual positive distances $\theta-\theta_-$ and $\theta_+-\theta$ .
$f'(\theta)$	The derivative at a locally smooth point approximated by the unequal-spacing three-point formula.

For analytic candidate $a$ and numerical value $n$ , the checker chooses $s=\max(1,|a|,|n|)$ and records $e=|a/s-n/s|$ . The floor at one behaves like absolute error for small gradients; larger gradients are compared relative to their magnitude. The candidate passes exactly when scaled error is no greater than the declared finite nonnegative tolerance.

Normalize finite gradient disagreement by the larger magnitude or one rust/crates/llm-from-scratch/src/autograd/gradcheck.rs#scale-aware-comparison

/// Compares gradients after scaling both by the larger magnitude or one.
pub fn compare_gradients(
    analytic: f64,
    numerical: f64,
    tolerance: f64,
) -> Result<GradientComparison, GradCheckError> {
    validate_tolerance(tolerance)?;
    if !analytic.is_finite() {
        return Err(GradCheckError::NonFiniteAnalyticGradient { value: analytic });
    }
    if !numerical.is_finite() {
        return Err(GradCheckError::NonFiniteNumericalGradient { value: numerical });
    }

    let scale = 1.0_f64.max(analytic.abs()).max(numerical.abs());
    let absolute_error = (analytic - numerical).abs();
    let scaled_error = (analytic / scale - numerical / scale).abs();

    Ok(GradientComparison {
        analytic,
        numerical,
        absolute_error,
        scale,
        scaled_error,
        tolerance,
        passed: scaled_error <= tolerance,
    })
}

/// Compares the recorded left and right secant slopes without claiming proof.
///
/// A failed diagnostic can indicate a nondifferentiable kink, a step that is
/// too large for the local curvature, or floating-point rounding. A passing
/// diagnostic cannot establish differentiability on its own.
pub fn compare_one_sided_slopes(
    difference: &CentralDifference,
    tolerance: f64,
) -> Result<OneSidedSlopeComparison, GradCheckError> {
    validate_tolerance(tolerance)?;
    let left = difference.left_slope;
    if !left.is_finite() {
        return Err(GradCheckError::NonFiniteOneSidedSlope {
            side: DifferenceSide::Minus,
            value: left,
        });
    }
    let right = difference.right_slope;
    if !right.is_finite() {
        return Err(GradCheckError::NonFiniteOneSidedSlope {
            side: DifferenceSide::Plus,
            value: right,
        });
    }
    let scale = 1.0_f64.max(left.abs()).max(right.abs());
    let absolute_gap = (left - right).abs();
    let scaled_gap = (left / scale - right / scale).abs();
    Ok(OneSidedSlopeComparison {
        left,
        right,
        absolute_gap,
        scale,
        scaled_gap,
        tolerance,
        consistent: scaled_gap <= tolerance,
    })
}

/// Runs a central difference and compares it with one analytic candidate.
///
/// The objective has the same smoothness, determinism, and side-effect
/// requirements as [`central_difference`].
pub fn scalar_gradient_check(
    point: f64,
    analytic: f64,
    step: f64,
    tolerance: f64,
    evaluate: impl FnMut(f64) -> f64,
) -> Result<ScalarGradientCheck, GradCheckError> {
    validate_tolerance(tolerance)?;
    let difference = central_difference(point, step, evaluate)?;
    let comparison = compare_gradients(analytic, difference.derivative, tolerance)?;
    Ok(ScalarGradientCheck {
        difference,
        comparison,
    })
}

From next-word backpropagation to checked Transformer training

Bengio et al.’s neural language model maximizes next-word log-likelihood with an explicit backward/update phase over output, hidden, and learned word-feature parameters. Those propagated derivatives make repeated training updates practical, but the implemented derivative path is not an independent check of itself.

The training procedure is detailed by Bengio et al., A Neural Probabilistic Language Model. Bengio et al. maximize next-word log-likelihood and publish a backward/update phase that propagates gradients through output units, hidden weights, and learned word-feature vectors.

The Transformer carries gradient-based training into repeated attention and feed-forward layers, using Adam for 100,000 base-model or 300,000 big-model steps. Baydin et al. distinguish numerical differentiation from reverse-mode automatic differentiation: finite differences estimate one local derivative from repeated evaluations, while reverse mode efficiently produces a scalar objective’s gradient over many parameters. Here, a materially separate numerical route can reveal a local mistake in a candidate derivative.

The later model evidence comes from Vaswani et al., Attention Is All You Need. Vaswani et al. train Transformer base models for 100,000 steps and big models for 300,000 steps, using Adam with an explicit learning-rate schedule.

For the numerical-method distinction, see Baydin et al., Automatic Differentiation in Machine Learning: a Survey. Baydin et al. describe centered finite differences, the truncation-versus-round-off step-size trade-off, poor scaling for full numerical gradients, and reverse mode’s efficiency for a scalar objective with many parameters.

This chapter uses three-point finite differences only as a slow sampled numerical cross-check for locally smooth objectives, including selected Chapter 12 indexed-mean-NLL derivatives, before Chapter 14 builds reverse mode. A passing sample is evidence for the chosen coordinates, probes, step, tolerance, and fixture—not proof of the complete gradient or local differentiability. The check does not train or run the decoder, and its policies are course-local.

The same passing sample also does not prove the formula’s derivation, either complete implementation, or behavior at nearby points and different inputs.

Neither model paper claims to use this checker. Finite differences are not a later LLM invention, a replacement training algorithm, or decoder runtime work. Their role here is narrower: give future backward code a materially separate, sampled forward-loss comparison before its values can update language-model parameters.

Implement a sampled numerical cross-check

Typed errors separate a bad configuration from a failed finite comparison. The checker names the rejected side and wraps tensor failures with the affected coordinate. It validates the complete selected coordinate set before the first objective call, so one late non-finite candidate cannot leave a partial result.

Name invalid steps, probes, values, shapes, samples, views, and coordinates rust/crates/llm-from-scratch/src/autograd/gradcheck.rs#gradcheck-errors

/// A rejected numerical derivative, comparison, sample request, or tensor read.
#[derive(Clone, Debug, PartialEq)]
pub enum GradCheckError {
    InvalidStep {
        step: f64,
    },
    NonFinitePoint {
        point: f64,
    },
    PerturbationNotFinite {
        side: DifferenceSide,
        point: f64,
        step: f64,
    },
    PerturbationUnchanged {
        side: DifferenceSide,
        point: f64,
        step: f64,
    },
    InvalidActualSpacing {
        side: DifferenceSide,
        point: f64,
        probe: f64,
        spacing: f64,
    },
    NonFiniteStencilWeight {
        side: DifferenceSide,
        value: f64,
    },
    NonFiniteEvaluation {
        side: DifferenceSide,
        value: f64,
    },
    NonFiniteOneSidedSlope {
        side: DifferenceSide,
        value: f64,
    },
    NonFiniteNumericalGradient {
        value: f64,
    },
    InvalidTolerance {
        tolerance: f64,
    },
    NonFiniteAnalyticGradient {
        value: f64,
    },
    GradientShapeMismatch {
        parameters: Vec<usize>,
        analytic: Vec<usize>,
    },
    ShapeOverflow,
    EmptyTensor,
    ZeroSamples,
    SampleAllocationFailed {
        samples: usize,
        rank: usize,
    },
    View(TensorViewError),
    AtCoordinate {
        coordinate: Vec<usize>,
        source: Box<GradCheckError>,
    },
}

impl fmt::Display for GradCheckError {
    fn fmt(&self, formatter: &mut fmt::Formatter<'_>) -> fmt::Result {
        match self {
            Self::InvalidStep { step } => write!(
                formatter,
                "central-difference requested step {step:?} must be positive and finite"
            ),
            Self::NonFinitePoint { point } => {
                write!(formatter, "gradient-check point {point:?} must be finite")
            }
            Self::PerturbationNotFinite { side, point, step } => write!(
                formatter,
                "{side} perturbation from point {point:?} by step {step:?} is not finite"
            ),
            Self::PerturbationUnchanged { side, point, step } => write!(
                formatter,
                "{side} perturbation from point {point:?} by step {step:?} rounds back to the point"
            ),
            Self::InvalidActualSpacing {
                side,
                point,
                probe,
                spacing,
            } => write!(
                formatter,
                "{side} actual spacing from point {point:?} to probe {probe:?} must be positive and finite, got {spacing:?}"
            ),
            Self::NonFiniteStencilWeight { side, value } => {
                write!(formatter, "{side} stencil weight is not finite: {value:?}")
            }
            Self::NonFiniteEvaluation { side, value } => {
                write!(formatter, "{side} function evaluation returned {value:?}")
            }
            Self::NonFiniteOneSidedSlope { side, value } => {
                write!(formatter, "{side} one-sided slope is not finite: {value:?}")
            }
            Self::NonFiniteNumericalGradient { value } => {
                write!(
                    formatter,
                    "central difference produced non-finite gradient {value:?}"
                )
            }
            Self::InvalidTolerance { tolerance } => write!(
                formatter,
                "gradient-check tolerance {tolerance:?} must be finite and nonnegative"
            ),
            Self::NonFiniteAnalyticGradient { value } => {
                write!(formatter, "analytic gradient {value:?} must be finite")
            }
            Self::GradientShapeMismatch {
                parameters,
                analytic,
            } => write!(
                formatter,
                "parameter shape {parameters:?} does not match analytic-gradient shape {analytic:?}"
            ),
            Self::ShapeOverflow => formatter.write_str("sampled tensor shape overflows usize"),
            Self::EmptyTensor => {
                formatter.write_str("sampled tensor gradient check needs at least one value")
            }
            Self::ZeroSamples => {
                formatter.write_str("sampled tensor gradient check needs at least one sample")
            }
            Self::SampleAllocationFailed { samples, rank } => write!(
                formatter,
                "cannot allocate {samples} sampled coordinates of rank {rank}"
            ),
            Self::View(error) => error.fmt(formatter),
            Self::AtCoordinate { coordinate, source } => {
                write!(
                    formatter,
                    "gradient check at coordinate {coordinate:?} failed: {source}"
                )
            }
        }
    }
}

impl Error for GradCheckError {
    fn source(&self) -> Option<&(dyn Error + 'static)> {
        match self {
            Self::View(error) => Some(error),
            Self::AtCoordinate { source, .. } => Some(source.as_ref()),
            _ => None,
        }
    }
}

impl From<TensorViewError> for GradCheckError {
    fn from(error: TensorViewError) -> Self {
        Self::View(error)
    }
}

Let $N$ be the number of elements in the nonempty tensor, let $R$ be the maximum coordinate count requested by the caller, and let $S=\min(R,N)$ be the number of coordinates the sampler actually selects. The Rust argument max_samples supplies $R$ . Zero requests and empty tensors are rejected, so $S\geq1$ .

When $S>1$ , the sampler evaluates the flat offset $\left\lfloor k(N-1)/(S-1)\right\rfloor$ for every integer $k\in\{0,1,\ldots,S-1\}$ . For shape [2,3], $N=6$ ; with $R=4$ , the actual count is $S=4$ , and $k=0,1,2,3$ gives [0,1,3,5]. Those offsets map to [[0,0],[0,1],[1,0],[1,2]]. No random generator or hidden state is involved. When $S=1$ , the formula’s denominator would be zero, so a separate branch selects $\lfloor N/2\rfloor$ ; for even $N$ , this is the larger of the two central flat offsets. The general sampler uses an overflow-safe intermediate and the same suffix-product shape semantics as the cumulative tensor storage.

Choose unique ordered tensor coordinates without randomness rust/crates/llm-from-scratch/src/autograd/gradcheck.rs#sample-tensor-coordinates

/// Selects unique, ordered coordinates without randomness or hidden state.
pub fn sample_tensor_coordinates(
    shape: &[usize],
    max_samples: usize,
) -> Result<Vec<SampledCoordinate>, GradCheckError> {
    if max_samples == 0 {
        return Err(GradCheckError::ZeroSamples);
    }
    let elements = element_count(shape)?;
    if elements == 0 {
        return Err(GradCheckError::EmptyTensor);
    }
    let samples = max_samples.min(elements);
    let mut coordinates = Vec::new();
    coordinates
        .try_reserve_exact(samples)
        .map_err(|_| GradCheckError::SampleAllocationFailed {
            samples,
            rank: shape.len(),
        })?;

    for sample in 0..samples {
        let flat_index = if samples == 1 {
            elements / 2
        } else {
            let numerator = (sample as u128) * ((elements - 1) as u128);
            (numerator / ((samples - 1) as u128)) as usize
        };
        coordinates.push(SampledCoordinate {
            flat_index,
            coordinate: coordinate_from_offset(shape, flat_index)?,
        });
    }
    Ok(coordinates)
}

For each selected coordinate, the tensor checker evaluates the same objective at the actual minus probe $\theta_-$ , the unperturbed value $\theta$ , and the actual plus probe $\theta_+$ , in that order. It restores the exact source bits immediately after each temporary probe, before the result is inspected and before every ordinary error return. A panic inside the learner-supplied objective is explicitly outside this small reference guarantee.

Check deterministic tensor coordinates and restore every temporary perturbation rust/crates/llm-from-scratch/src/autograd/gradcheck.rs#sampled-tensor-gradient-check

/// Checks deterministic tensor coordinates while restoring every probed value.
///
/// The parameter tensor is restored on every ordinary `Ok` or `Err` path. A
/// panic inside `objective` is deliberately outside this reference guarantee.
/// The objective must be deterministic, side-effect-free, and differentiable
/// with sufficient smoothness across every sampled probe interval; known or
/// possible kinks are outside this comparison's contract.
pub fn sampled_tensor_gradient_check(
    parameters: &mut Tensor,
    analytic: &TensorView<'_>,
    step: f64,
    tolerance: f64,
    max_samples: usize,
    mut objective: impl FnMut(&Tensor) -> f64,
) -> Result<TensorGradientCheck, GradCheckError> {
    validate_step(step)?;
    validate_tolerance(tolerance)?;
    if parameters.shape() != analytic.shape() {
        return Err(GradCheckError::GradientShapeMismatch {
            parameters: parameters.shape().to_vec(),
            analytic: analytic.shape().to_vec(),
        });
    }

    let samples = sample_tensor_coordinates(parameters.shape(), max_samples)?;
    let mut candidates = Vec::new();
    candidates.try_reserve_exact(samples.len()).map_err(|_| {
        GradCheckError::SampleAllocationFailed {
            samples: samples.len(),
            rank: parameters.rank(),
        }
    })?;

    // Validate every selected coordinate before the first objective call.
    for sample in &samples {
        let point = parameters.as_slice()[sample.flat_index];
        perturbations(point, step).map_err(|error| at_coordinate(&sample.coordinate, error))?;
        let analytic_value = *analytic
            .get(&sample.coordinate)
            .map_err(GradCheckError::View)?;
        if !analytic_value.is_finite() {
            return Err(at_coordinate(
                &sample.coordinate,
                GradCheckError::NonFiniteAnalyticGradient {
                    value: analytic_value,
                },
            ));
        }
        candidates.push((sample.clone(), point, analytic_value));
    }

    let mut checks = Vec::new();
    checks.try_reserve_exact(candidates.len()).map_err(|_| {
        GradCheckError::SampleAllocationFailed {
            samples: candidates.len(),
            rank: parameters.rank(),
        }
    })?;

    for (sample, point, analytic_value) in candidates {
        let difference = central_difference(point, step, |probe| {
            parameters.as_mut_slice()[sample.flat_index] = probe;
            let value = objective(parameters);
            parameters.as_mut_slice()[sample.flat_index] = point;
            value
        })
        .map_err(|error| at_coordinate(&sample.coordinate, error))?;
        let comparison = compare_gradients(analytic_value, difference.derivative, tolerance)
            .map_err(|error| at_coordinate(&sample.coordinate, error))?;
        checks.push(CoordinateGradientCheck {
            flat_index: sample.flat_index,
            coordinate: sample.coordinate,
            difference,
            comparison,
        });
    }

    Ok(TensorGradientCheck {
        shape: parameters.shape().to_vec(),
        requested_samples: max_samples,
        passed: checks.iter().all(|check| check.comparison.passed),
        checks,
    })
}

The numerical objective uses Chapter 12 indexed mean NLL. Its shape [2,3] logits are [0,1,-1,2,0,-2], its targets are [0,2], and its forward mean loss is 2.775268796472. The separate analytic routine starts from each raw logit row. For row $r$ , it computes

m_r=\max_j \ell_{r,j},\qquad p_{r,c}=\frac{\exp(\ell_{r,c}-m_r)}{\sum_j\exp(\ell_{r,j}-m_r)}.

For vocabulary class $c$ and target $t_r$ , the manually derived candidate is

\frac{\partial\mathcal L}{\partial \ell_{r,c}} =\frac{p_{r,c}-\mathbf 1[c=t_r]}{2}.

At each target-logit position, subtract one from the locally computed probability, then divide every row gradient by the two target rows. This analytic routine does not call the production softmax or indexed_mean_nll; the perturbed numerical objective does call indexed_mean_nll. Separating those material algorithm paths can expose a mistake shared by neither one, but it is not total independence: both sides still share the raw logits and targets, f64 elementary arithmetic, Tensor storage, and the same row-major index convention.

Compute the analytic NLL candidate locally from raw logits rust/demos/ch13-gradient-checking/src/lib.rs#hand-derived-nll-gradient

/// Applies `(normalized probabilities - one_hot(target)) / batch_size` through a separate local path.
///
/// This frozen two-row analytic path deliberately implements its own row traversal,
/// maximum shift, exponential sum, normalization, target subtraction, and
/// batch scaling without calling the production probability or indexed-NLL
/// helpers. Both paths still share the input values and targets, IEEE `f64`
/// arithmetic and its primitive exponential, `Tensor` storage, and row-major
/// index conventions.
pub fn hand_derived_nll_gradient(logits: &Tensor) -> Result<Tensor, Box<dyn Error>> {
    if logits.shape() != LOGIT_SHAPE {
        return Err(format!(
            "the local Chapter 13 analytic path requires shape {LOGIT_SHAPE:?}, got {:?}",
            logits.shape()
        )
        .into());
    }

    let columns = LOGIT_SHAPE[1];
    let mut values = vec![0.0; LOGIT_VALUES.len()];
    for (row, &target) in TARGETS.iter().enumerate() {
        let start = row * columns;
        let row_logits = &logits.as_slice()[start..start + columns];
        let maximum = row_logits.iter().copied().fold(f64::NEG_INFINITY, f64::max);
        if !maximum.is_finite() {
            return Err(format!("oracle row {row} has no finite maximum").into());
        }

        let mut normalizer = 0.0;
        for (column, &logit) in row_logits.iter().enumerate() {
            if !logit.is_finite() {
                return Err(format!("oracle logit at [{row}, {column}] is not finite").into());
            }
            let weight = (logit - maximum).exp();
            values[start + column] = weight;
            normalizer += weight;
        }
        if !normalizer.is_finite() || normalizer <= 0.0 {
            return Err(
                format!("oracle row {row} has invalid normalization {normalizer:?}").into(),
            );
        }

        for column in 0..columns {
            values[start + column] /= normalizer;
        }
        values[start + target] -= 1.0;
    }
    for gradient in &mut values {
        *gradient /= TARGETS.len() as f64;
    }
    Ok(Tensor::from_vec(LOGIT_SHAPE.to_vec(), values)?)
}

Compare four locally derived candidates with perturbed indexed-NLL evaluations rust/demos/ch13-gradient-checking/src/lib.rs#sampled-nll-gradient-check

/// Checks four deterministic vocabulary-logit coordinates against indexed NLL.
pub fn tiny_nll_gradient_example() -> Result<TinyNllGradientExample, Box<dyn Error>> {
    let mut logits = logits()?;
    let analytic = hand_derived_nll_gradient(&logits)?;
    let loss = indexed_mean_nll(&logits.view(), 1, &TARGETS)?;
    let original_bits = logits
        .as_slice()
        .iter()
        .map(|value| value.to_bits())
        .collect::<Vec<_>>();
    let check = sampled_tensor_gradient_check(
        &mut logits,
        &analytic.view(),
        TENSOR_STEP,
        TENSOR_TOLERANCE,
        TENSOR_SAMPLES,
        |candidate| {
            indexed_mean_nll(&candidate.view(), 1, &TARGETS)
                .expect("the frozen finite logits and targets remain valid")
        },
    )?;
    let restored_exactly = logits
        .as_slice()
        .iter()
        .map(|value| value.to_bits())
        .eq(original_bits);

    Ok(TinyNllGradientExample {
        logits,
        analytic,
        loss,
        check,
        restored_exactly,
    })
}

The demo prepares both scalar candidates, the six-step scan, all four NLL coordinates, exact restoration, and one collapsed-step error:

Prepare the complete deterministic learner evidence before printing rust/demos/ch13-gradient-checking/src/main.rs#learner-gradient-check-output

    let correct = quadratic_gradient_check(6.0)?;
    let wrong = quadratic_gradient_check(5.5)?;
    let rounded = rounded_identity_gradient_check()?;
    let kink = absolute_kink_diagnostic()?;
    let scan = cubic_step_scan()?;
    let nll = tiny_nll_gradient_example()?;
    let collapsed = central_difference(1.0, 1.0e-20, |point| point * point).unwrap_err();

./course run cargo run --quiet --locked -p ch13-gradient-checking

quadratic: theta=3.000000000000 requested_h=1.00000000000000006e-1 actual_h_minus=1.00000000000000089e-1 actual_h_plus=1.00000000000000089e-1 f_minus=8.410000000000 f_center=9.000000000000 f_plus=9.610000000000 left_slope=5.900000000000 right_slope=6.100000000000 left_weight=5.00000000000000000e-1 right_weight=5.00000000000000000e-1 numerical=6.000000000000
wrong candidate: analytic=5.500000000000 scaled_error=8.333333333333e-2 tolerance=1.000000000000e-6 pass=false
nll logits: shape=[2, 3] values=[0.0, 1.0, -1.0, 2.0, 0.0, -2.0] targets=[0, 2] loss=2.775268796472
sampled coordinates: [[0, 0], [0, 1], [1, 0], [1, 2]]
tensor restored exactly: true

Watch the step size help, then hurt

For $g(\theta)=\theta^3-2\theta$ at $\theta=1.5$ , the analytic derivative is $4.75$ . The Rust scan starts with requested $h=1$ , improves through $10^{-1}$ and $10^{-3}$ , and reaches its lowest displayed scaled error at requested $h=10^{-5}$ : the numerical value is $4.750000000100$ , the scaled error is $2.103583973678\times10^{-11}$ , and the record passes. Rounding then makes requested $h=10^{-13}$ produce $4.751111111111$ with scaled error $2.338634237605\times10^{-4}$ , so it fails. At requested $h=10^{-15}$ , the numerical value is $4.800000000000$ and scaled error is $1.041666666667\times10^{-2}$ , which also fails. The lesson does not promote $10^{-5}$ as a universal optimum; it is the trustworthy middle for this worked example.

The figure places requested and actual probe spacing, the rounded identity check, the known nondifferentiable corner at $|x|$ when $x=0$ , separate analytic and numerical NLL paths with their shared assumptions, every error regime, verdict, selected coordinate, and rejected request in one view. Its scan records rounded estimates and scaled errors only for readability; the complete executable output retains their exact values together with the probe points and function values. Every displayed $h$ is the requested step; the derivative uses the two actual representable distances. The scan reveals the useful smooth middle region between truncation and rounding error, while the token-loss samples show the deliberately limited scope of a successful coordinate check.

Run one cubic derivative check across all six step sizes rust/demos/ch13-gradient-checking/src/lib.rs#step-size-scan

/// Runs the same central difference from truncation-dominated to rounded probes.
pub fn cubic_step_scan() -> Result<Vec<StepScanRecord>, GradCheckError> {
    STEP_SCAN
        .iter()
        .zip(STEP_PHASES)
        .map(|(&step, phase)| {
            scalar_gradient_check(CUBIC_POINT, CUBIC_ANALYTIC, step, STEP_TOLERANCE, cubic)
                .map(|check| StepScanRecord { phase, check })
        })
        .collect()
}

Compare actual probe spacing before trusting a sampled check

Compare requested and actual probe spacing, a rounded identity check, a known nondifferentiable corner, separate NLL paths and their shared assumptions, selected token-loss coordinates, exact restoration, and invalid inputs.

Quadratic numerical derivative: 6.000000000000
Cubic scan point: $\theta=1.500000000000$
Indexed mean NLL: 2.775268796472

Check actual spacing around theta equals three

Actual minus probe: $\theta_-=2.9$; $q(\theta_-)=8.41$
Actual left spacing: 1.00000000000000089e-1
Left one-sided slope: 5.900000000000
Center: $\theta=3,\;q(\theta)=9$
Requested step: 1.00000000000000006e-1
Probe-spacing stencil: symmetric; $w_-$ 5.00000000000000000e-1 $w_+$ 5.00000000000000000e-1
Numerical gradient: 6.000000000000
Actual plus probe: $\theta_+=3.1$; $q(\theta_+)=9.61$
Actual right spacing: 1.00000000000000089e-1
Right one-sided slope: 6.100000000000

Scan six step sizes around theta equals one point five

$h=1$ over tolerance

Numerical gradient

$5.75$

Scaled error

$1.74\times10^{-1}$

Actual left spacing / Actual right spacing

1.00000000000000000e0 1.00000000000000000e0

Error regime

truncation
$h=10^{-1}$ over tolerance

Numerical gradient

$4.76$

Scaled error

$2.10\times10^{-3}$

Actual left spacing / Actual right spacing

1.00000000000000089e-1 1.00000000000000089e-1

Error regime

truncation
$h=10^{-3}$ within tolerance

Numerical gradient

$4.750001$

Scaled error

$2.11\times10^{-7}$

Actual left spacing / Actual right spacing

9.99999999999889866e-4 9.99999999999889866e-4

Error regime

converging
$h=10^{-5}$ within tolerance

Numerical gradient

$4.7500000001$

Scaled error

$2.10\times10^{-11}$

Actual left spacing / Actual right spacing

1.00000000000655120e-5 1.00000000000655120e-5

Error regime

trusted range
$h=10^{-13}$ over tolerance

Numerical gradient

$4.751111$

Scaled error

$2.34\times10^{-4}$

Actual left spacing / Actual right spacing

9.99200722162640886e-14 9.99200722162640886e-14

Error regime

rounding
$h=10^{-15}$ over tolerance

Numerical gradient

$4.8$

Scaled error

$1.04\times10^{-2}$

Actual left spacing / Actual right spacing

1.11022302462515654e-15 1.11022302462515654e-15

Error regime

rounding

Compare two finite gradient candidates

Numerical gradient $6$ Tolerance $10^{-6}$

Analytic candidate $6$ within tolerance

Scaled error

$0$
Analytic candidate $5.5$ over tolerance

Scaled error

$8.33\times10^{-2}$

Cross-check four selected coordinates of the token loss

Each record is evidence only for its selected coordinate, objective, probes, step, tolerance, and fixture.

Original bits restored yes, exactly

$0\to[0,0]$ ✓

Analytic candidate

-0.377635764473

Numerical gradient

-0.377635764481

Scaled error

$8.75\times10^{-12}$
$1\to[0,1]$ ✓

Analytic candidate

0.332620477887

Numerical gradient

0.332620477894

Scaled error

$6.43\times10^{-12}$
$3\to[1,0]$ ✓

Analytic candidate

0.433406666099

Numerical gradient

0.433406666087

Scaled error

$1.21\times10^{-11}$
$5\to[1,2]$ ✓

Analytic candidate

-0.492061880012

Numerical gradient

-0.492061879994

Scaled error

$1.75\times10^{-11}$

Read the method boundaries before trusting agreement

Unequal rounded probes require actual spacing; a known kink can expose a false pass, and sampled agreement is not proof.

$f(x)=x,\;x=1$ within tolerance

Requested step: 1.33226762955018780e-16
Actual left spacing: 1.11022302462515654e-16
Actual right spacing: 2.22044604925031308e-16
Probe-spacing stencil: unequal
Numerical gradient: 1.000000000000

$f(x)=|x|,\;x=0$

Left one-sided slope: -1.000000000000
Right one-sided slope: 1.000000000000
Numerical gradient: 0.000000000000
Tolerance: 1.000000000000e-12

The analytic NLL path is local and does not call production probability or indexed-NLL helpers; both paths still share input values and targets, IEEE f64 arithmetic and its primitive exponential, Tensor storage, and row-major index conventions.

Analytic candidate: local-row-max-exp-sum-normalize-target-gradient
objective-path: indexed-mean-nll
shared-primitives: f64-exp,frozen-inputs-and-targets

Reject unsafe numerical requests

Unsafe requests stop before a derivative is accepted.

zero step

$h=0$
step leaves the point unchanged

minus probe $\theta=1$ $h=10^{-20}$
non-finite objective value

minus probe $f=\mathrm{NaN}$
shape mismatch

$[2]\ne[1,2]$

Predict before running Rust

Calculate $q(2.9)$ and $q(3.1)$ , then explain why the implementation must derive $h_-$ and $h_+$ instead of dividing by the requested $2h$ .
For $f(x)=x$ and unequal nonzero actual distances, predict both one-sided slopes and the weighted three-point result.
Predict the pass states of candidates $6$ and $5.5$ at tolerance $10^{-6}$ .
Select the lowest-error result in this six-step scan without calling its requested step universally optimal.
Explain why requested $h=10^{-13}$ and $h=10^{-15}$ both fail after the larger trusted step passes.
Compute the wrong candidate’s scale and scaled error before reading the answer.
For a tensor with $N=6$ elements, max_samples is $R=4$ . Compute $S=\min(R,N)$ , then evaluate $\left\lfloor k(N-1)/(S-1)\right\rfloor$ for every $k\in\{0,1,2,3\}$ and map the four flat offsets to shape [2,3] coordinates.
Predict zero-step and collapsed-probe behavior plus tensor restoration on ordinary errors.
For $f(x)=|x|$ at $x=0$ , explain why numerical value $0$ can agree with candidate $0$ without establishing a derivative.
Misconception check: what exactly do four passing NLL coordinates establish, and what do they leave shared or untested?

Check the predictions

$q(2.9)=8.41$ and $q(3.1)=9.61$ in the rounded hand calculation. Actual floating-point subtraction and addition can place the two probes at unequal distances, so those distances—not the requested $h$ —must determine the formula.
Both quotients are $1$ : $(f(x)-f(x_-))/h_-=1$ and $(f(x_+)-f(x))/h_+=1$ . Their weights sum to one, so the result remains $1$ within floating-point tolerance even when $h_-\ne h_+$ .
Candidate $6$ passes. Candidate $5.5$ returns a finite comparison with pass=false.
Requested $h=10^{-5}$ has the lowest displayed scaled error in this fixed scan: $2.103583973678\times10^{-11}$ for numerical value $4.750000000100$ . Another function, point, floating type, or pair of actual probes can have a different best region.
At requested $h=10^{-13}$ , the numerical value $4.751111111111$ has scaled error $2.338634237605\times10^{-4}$ and fails. At $h=10^{-15}$ , the actual spacing is only $1.11022302462515654\times10^{-15}$ on each side; the one-sided slopes become $4.4$ and $5.2$ , so their weighted value $4.800000000000$ has scaled error $1.041666666667\times10^{-2}$ and also fails. Rounded nearby values no longer retain enough low-bit information for a trustworthy slope.
In exact arithmetic, $\max(1,5.5,6)=6$ and $|5.5/6-6/6|=1/12$ ; the printed f64 result is 8.333333333333e-2, still far above $10^{-6}$ .
Here $S=\min(4,6)=4$ . For $k=0,1,2,3$ , $\lfloor k(6-1)/(4-1)\rfloor=\lfloor 5k/3\rfloor$ gives flats [0,1,3,5], hence coordinates [[0,0],[0,1],[1,0],[1,2]].
Zero $h$ is invalid; $h=10^{-20}$ at point $1$ rounds the minus probe back to the point. On every ordinary error, a temporary tensor value is restored before return.
The symmetric secant has slope $0$ , but the left slope is $-1$ and the right slope is $1$ . The derivative at the corner does not exist, so agreement with $0$ is a false pass outside the local-smoothness precondition.
The results support only the four selected coordinates for this objective, fixture, actual probes, requested step, and tolerance. They neither prove the full gradient nor train or run the decoder. The analytic and numerical routes also share inputs, f64, Tensor, and index conventions, so those shared assumptions remain untested.

Run the example after predicting:

./course run cargo run --quiet --locked -p ch13-gradient-checking

Prepare reverse-mode differentiation

The cumulative project can now compare selected hand-derived vocabulary-logit derivatives from a locally implemented analytic path with perturbed evaluations of the production indexed mean NLL. The analytic path does not call the production softmax or indexed_mean_nll, but both sides still share the raw inputs, f64 arithmetic, Tensor storage, and index conventions. Agreement is therefore useful sampled evidence, not proof of either complete implementation. Chapter 14 will apply the same boundary to reverse-mode scalar derivatives.

The next chapter will store scalar computation nodes, apply local derivative rules in reverse topological order, and accumulate contributions when one value is reused. Its analytic results will not be trusted merely because the backward pass runs: the three-point checker from this chapter remains a materially separate sampled numerical comparison at locally smooth points.