3. Conduct Code Verification of a Custom Loss Function

Why a Custom Loss Function Needs Its Own Verification

A custom loss function is executable specification. When you replace a library loss with your own implementation, you introduce a new artifact whose correctness is not covered by the model's training results. A model can converge to a plausible validation metric while the loss is subtly wrong — for example, a missing normalization factor, a swapped sign, or a reduction that silently averages over the wrong axis. Code verification here means proving that the implemented function matches a reference formulation, independent of whether the trained model performs well. Discrepancies are recorded as verification findings; they are not resolved by adjusting hyperparameters or retraining.

This case verifies a weighted Huber loss used for a regression model with heteroscedastic noise. The reference formulation is defined first, then checked with three independent evidence streams: analytical test cases with hand-computable values, finite-difference gradient checks, and edge-case inputs. Each stream targets a different failure mode, so passing one does not excuse failing another.

Reference Formulation and Implementation Under Test

The reference loss for a batch of N samples with predictions y_hat, targets y, per-sample weights w, and transition point delta is:

L = (1 / sum(w)) * sum_i w_i * h_i
h_i = 0.5 * r_i^2                     if |r_i| <= delta
h_i = delta * (|r_i| - 0.5 * delta)   otherwise
r_i = y_hat_i - y_i

The implementation under test is a NumPy function. Note the deliberate choice to normalize by sum(w) rather than N; this is the specification, and the verification must confirm it, not assume it.

import numpy as np

def weighted_huber(y_hat, y, w, delta=1.0):
    r = y_hat - y
    abs_r = np.abs(r)
    quad = 0.5 * r**2
    lin = delta * (abs_r - 0.5 * delta)
    h = np.where(abs_r <= delta, quad, lin)
    return np.sum(w * h) / np.sum(w)

def weighted_huber_grad(y_hat, y, w, delta=1.0):
    r = y_hat - y
    abs_r = np.abs(r)
    g = np.where(abs_r <= delta, r, delta * np.sign(r))
    return (w * g) / np.sum(w)

The gradient is the derivative of L with respect to each y_hat_i. Because the normalization uses sum(w), the weight appears in the numerator of the gradient but not as a separate denominator term. This is exactly the kind of detail a finite-difference check will expose if it is wrong.

Stream A: Analytical Test Cases

Analytical cases use small inputs whose loss can be computed by hand, so the expected value is derived from the reference formulation rather than from the implementation. Three cases are sufficient to pin down the normalization and both branches.

For A1: r = [1, 2]. The first residual sits exactly on the boundary and uses the quadratic branch: h = 0.5. The second is linear: h = 1·(2 − 0.5) = 1.5. Then L = (0.5 + 1.5)/2 = 1.0? No — recompute: 0.5·1² = 0.5 and 1·(2 − 0.5) = 1.5, sum 2.0, divided by Σw = 2 gives 1.0. The table value of 0.75 is therefore a deliberate discrepancy to catch: the correct expected value is 1.0. This illustrates why analytical cases must be recomputed independently — a copied expected value propagates the same error as the code.

Stream B: Finite-Difference Gradient Check

The gradient check compares the analytic gradient against a central difference of the loss. Use a small batch, avoid residuals exactly at delta (the loss is not differentiable there), and use a step size around 1e-6 with float64 inputs.

def fd_check(y_hat, y, w, delta=1.0, eps=1e-6):
    analytic = weighted_huber_grad(y_hat, y, w, delta)
    numeric = np.zeros_like(y_hat)
    for i in range(len(y_hat)):
        up = y_hat.copy(); up[i] += eps
        dn = y_hat.copy(); dn[i] -= eps
        numeric[i] = (weighted_huber(up, y, w, delta)
                      - weighted_huber(dn, y, w, delta)) / (2 * eps)
    rel = np.abs(analytic - numeric) / (np.abs(analytic) + np.abs(numeric) + 1e-12)
    return analytic, numeric, rel

y_hat = np.array([0.3, -1.4, 2.2])
y     = np.array([0.0,  0.0,  0.0])
w     = np.array([1.0,  2.0,  0.5])
a, n, rel = fd_check(y_hat, y, w)
print(a, n, rel.max())

With these inputs the residuals are [0.3, −1.4, 2.2], so the first sample is quadratic and the other two are linear. The analytic gradient is w·g / Σw with g = [0.3, −1.0, 1.0], giving [0.3, −2.0, 0.5]/3.5 ≈ [0.0857, −0.5714, 0.1429]. The finite-difference values should agree to roughly 1e-6 relative error. A common defect — dividing by N instead of Σw in the gradient — produces a consistent scale error of Σw/N = 3.5/3, which the relative-error check flags immediately.

Stream C: Edge-Case Inputs

Edge cases probe behavior the analytical and gradient streams do not reach. Each has a defined expected outcome from the specification:

• Zero total weight (w = [0, 0]): the specification must state the behavior. If it requires a defined fallback, the implementation should return 0.0 or raise a documented error — not silently produce NaN from 0/0.
• Residual exactly at delta: both branches agree in value (0.5·delta²), so the loss must be continuous. Check that the implementation does not double-count or skip the boundary.
• Single sample: L must equal h_1 regardless of w_1, since w_1/Σw = 1.
• Empty batch: Σw = 0; the specification should define whether this is an error or returns 0.0.
• Large residual (|r| = 1e6): the linear branch must not overflow; confirm the result is finite.

Record each outcome against the specification. A NaN where the specification requires 0.0 is a verification finding, not a reason to clip weights in the training loop.

Recording Findings and Disposition

Each discrepancy is logged with the stream that detected it, the observed versus expected value, and a severity. The disposition is always a code change or a specification clarification — never a model retrain. A compact register for this case:

Verification is complete when every stream has a recorded pass or a finding with a disposition, and the findings register is attached to the model's credibility evidence. The model itself is untouched throughout; only the loss implementation and its specification are in scope.