Numerical computation, conditioning, and reliable experiments

An equation specifies an operation over an ideal number system. A numerical program stores a finite approximation, performs a sequence of rounded operations, and returns an approximation to that equation’s result. Distinguish input representation, algorithmic arithmetic, and the approximation inherent in the model. Increasing precision can reduce rounding while leaving an incorrect derivative, an ill-conditioned inverse problem, or a biased statistical estimate unchanged. Before debugging decimals, identify the mathematical quantity, its scale, and the error that matters to its use.

A normal binary floating-point number has a sign, a significand with a fixed number p of significant binary digits, and an exponent from a finite range. Within one binade, the interval from 2^e inclusive to 2^(e+1) exclusive, consecutive normal numbers are spaced by 2^(e−p+1). Spacing grows with magnitude; the same absolute increment can change a small number and disappear when added to a large one. The exponent controls range while the significand controls local relative precision. They solve different storage problems.

The familiar IEEE binary16, binary32, and binary64 formats have p equal to 11, 24, and 53 respectively, including the implicit leading digit for normal numbers. Their spacing immediately above one, called eps by NumPy, is 2^(1−p). Under rounding to nearest with ties to even, the unit roundoff u is eps/2. Do not substitute eps for a universal minimum resolvable change: spacing depends on magnitude, and the spacing below a power of two differs from the spacing above it. Inspect the actual dtype after conversions and reductions.

For an elementary operation on normal results, absent overflow and underflow, the standard model is fl(a op b)=(a op b)(1+δ), with |δ|≤u. This is a model of the rounded exact operation on stored inputs. It does not say that the stored a equals the original real measurement, nor that a long algorithm has relative error at most u. When the exact operation gives zero, absolute reasoning is more useful. Near the underflow range, subnormal numbers provide gradually decreasing absolute spacing but can lose relative precision.

Overflow occurs when a result exceeds the finite range and may produce infinity; later infinity minus infinity or infinity divided by infinity produces NaN. Underflow may produce a subnormal value or zero. Some sufficiently tiny probabilities are thus stored as zero even when the model assigns positive mass. NaN is not an ordinary missing number that can safely participate in comparisons: many comparisons with it return false. Check finiteness at defined boundaries, and decide whether an infinite mathematical result is legitimate, as for unsupported predictions in Module 27, or an implementation failure.

Converting an already rounded binary32 array to binary64 preserves its rounded values; it does not restore the lost original bits. Conversely, computing a sensitive reduction in binary64 before storing its final output in binary32 may preserve information that a binary32 reduction discards. State input, accumulation, and output precision separately. Mixed precision is a design choice requiring an error argument for each sensitive stage, rather than a promise that a large dtype name anywhere makes the whole pipeline accurate.

Given an ideal reference q and approximation q_hat, absolute error is |q_hat−q|. Relative error divides this by |q| when q is nonzero. An absolute error of 10^(−8) is tiny relative to a million and large relative to 10^(−12). At a zero target, relative error is undefined; use an absolute scale derived from the application. For vectors choose a norm and report it. A componentwise error can reveal a badly computed small component hidden inside a small overall relative norm error.

Subtraction is not automatically inaccurate. The difference between two nearby stored floating-point values can even be exactly representable. The dangerous issue is that previously introduced errors in their operands can become large relative to a small desired difference. If a and b each carry perturbations, their difference carries the difference of those perturbations. When |a−b| is tiny compared with |a|+|b|, many accurate digits in the inputs may be needed to determine accurate digits of the result. This sensitivity is called cancellation.

Consider sqrt(1+x)−1 for small positive x. Forming 1+x can erase x, and subtracting one then gives zero. Multiplying numerator and denominator by sqrt(1+x)+1 yields x/(sqrt(1+x)+1), the same real function where defined. Its denominator is near two and does not subtract nearby values. This transformation preserves a tiny stored x without first adding it to one. Algebraic equivalence over reals therefore does not imply equal numerical behaviour. Rewriting can remove avoidable algorithmic error without changing the target.

For a sum S of n stored inputs, cancellation sensitivity is measured by Σ|x_i|/|S| when S≠0. A large value means small individual perturbations can greatly change the final relative answer. Under the usual rounding model, sequential summation has an absolute error bounded by γ_(n−1)Σ|x_i|, where γ_k=ku/(1−ku) and ku<1. This is a named standard bound, not a universal claim under overflow or arbitrary execution rules. Its division by |S| exposes why relative error can explode for a cancelling sum.

Pairwise reduction reduces the depth of arithmetic dependencies from n to roughly log2 n, giving an improved bound under corresponding assumptions. Compensated schemes retain low-order parts otherwise discarded by additions. Python math.fsum is useful for a high-quality sum of binary floating inputs; it does not infer their original decimal intention. Different reduction orders, including threaded library implementations, can change final low-order bits. Rearranging terms by magnitude is not a universal cure for cancellation, and a different reduction does not repair a badly conditioned data problem.

Numerical comparisons need an explicit reference and scale. A common rule is |a−b|≤atol+rtol|b|, treating b as reference; NumPy allclose uses this asymmetry. Near zero, atol controls acceptance. For quantities with different units, choose tolerances per quantity rather than one dimensionless magic constant. A small rtol with a large atol may accept a completely wrong tiny probability. Equality of integer counts, task IDs, shapes, and exact Boolean contracts should remain exact when their mathematical contract is exact.

Bounds and branch decisions require particular care. An eigenvalue slightly below zero may reflect rounding of a PSD matrix or a genuinely indefinite input. Compare its magnitude with a scale-aware numerical tolerance, report that tolerance, and retain a clear distinction between a mathematical proof and a numerical check. An input outside a model’s domain is not made valid by silently clipping it. Likewise, clipping probabilities before taking logs changes an infinite or very large loss into a different quantity; it should not be passed off as a stable evaluation of the original objective.

Conditioning belongs to a mathematical problem with a declared input/output metric. For a differentiable scalar map f at nonzero x with nonzero f(x), a local relative condition number is |x f’(x)/f(x)|: a small relative input change is amplified by about this factor to first order. For f(x)=x² away from zero it equals two. For f(x)=x−1 near one it becomes |x/(x−1)| and can be enormous. At zeros or nondifferentiable points, another absolute or local definition is needed.

The declared metric matters in practice. A parameter vector mixing metres, seconds and dimensionless coefficients does not have an intrinsically meaningful unweighted Euclidean error. Changing units can change a matrix condition estimate without changing the physical experiment. Choose nondimensional coordinates or an application-motivated weighting before interpreting the number. If predictions are the scientific target, inspect their sensitivity separately from coefficient sensitivity: nearly collinear predictors can allow unstable coefficient trade-offs while their combined fitted predictions remain similar on the observed design. New inputs outside that design can expose the unstable direction again. This is a reason to report the target quantity, rather than one condition number as a verdict on the whole model.

Algorithmic stability concerns how a particular computation responds to arithmetic errors. Forward error compares its output to the exact answer for the intended stored input. Backward error asks whether its output exactly solves a nearby problem and how nearby that input is. A backward-stable algorithm can deliver large forward error for an ill-conditioned problem: it faithfully solves a slightly changed input whose mathematical solution is very different. Accurate arithmetic and insensitive mathematics are separate properties, and reporting one does not establish the other.

For nonsingular A in a linear solve Ax=b, the induced two-norm condition number is κ2(A)=||A||2||A^(−1)||2=σ_max/σ_min. Multiplying the entire matrix and right-hand side by the same nonzero scalar leaves this condition number unchanged. Unequal row or column scalings can change the coordinate metric and conditioning; they must preserve the original model when transformed back. A singular matrix has no ordinary inverse, and this condition number is infinite. Numerical rank uses a declared threshold and is a precision-dependent decision.

Hold A fixed and perturb b by δb. Then δx=A^(−1)δb, so ||δx||≤||A^(−1)||||δb||. Since ||b||≤||A||||x||, division gives relative solution error at most κ(A) times relative right-hand-side perturbation. This is a worst-direction bound, not an equality for every perturbation. Direction matters: an input change aligned with a small singular mode is amplified most. A high condition number predicts a possible sensitivity; it does not prove every observed solve is inaccurate.

If A changes as well, write (A+ΔA)(x+δx)=b+δb. Rearranging gives δx=A^(−1)(δb−ΔA x−ΔA δx). Apply norms and move the last term left. If κ(A) ε_A<1, with ε_A=||ΔA||/||A|| and ε_b=||δb||/||b||, this yields a relative bound κ(A)(ε_A+ε_b)/(1−κ(A)ε_A). The denominator states when the perturbative argument is useful. Near singularity, even tiny changes can invalidate a reassuring-looking bound.

The residual r=b−A x_hat is available without knowing the true x. With fixed A, x_hat exactly solves the perturbed right-hand side b−r, so ||r||/||b|| is a right-hand-side backward-error measure when b≠0. Forward error obeys ||x_hat−x||/||x||≤κ(A)||r||/||b||. A small relative residual therefore needs a condition estimate before it implies a small solution error. Evaluate the residual with sufficient precision; computing it in the same low precision as the solve can hide discrepancies.

For a square nonsingular system, solve using a suitable factorisation rather than explicitly forming an inverse and multiplying it by b. The inverse requires extra work and storage and adds another stage of rounding. A dense general solve typically uses pivoted elimination; a symmetric positive-definite system admits Cholesky. Structure must be established, not inferred from a variable name. A reported factorisation failure can mean an invalid mathematical assumption, a near-singular rounded input, or a threshold decision that should be examined.

Least squares minimises ||Xw−y||² with n rows and d columns. For full column rank, its normal equations are XᵀXw=Xᵀy. If X=UΣVᵀ, the nonzero singular values of X are σ_i and the eigenvalues of XᵀX are σ_i². Therefore κ2(XᵀX)=κ2(X)² exactly for the real full-rank matrix. Forming the Gram matrix can lose a small mode before the subsequent solve begins. Its computed condition number need not equal the exact square once rounding has changed that matrix.

With a reduced QR decomposition X=QR and orthonormal columns of Q, split y into its projection QQᵀy and its orthogonal remainder. Pythagoras shows that minimising the residual is equivalent to minimising ||Rw−Qᵀy||². For nonsingular R, solve Rw=Qᵀy by triangular substitution. This avoids explicitly squaring the singular values through a Gram matrix. Backward stability is a property of suitable implementations and assumptions, rather than of every hand-written routine that merely calls itself QR.

SVD reveals rank and directions of sensitivity directly. In full rank the least-squares coefficients are VΣ^(−1)Uᵀy; small σ_i multiply projected noise by large reciprocals. At deficient rank, a pseudoinverse returns a minimum-norm solution under a chosen rank threshold. Truncating a small singular value changes the target solution, which may be desirable as regularisation but must be stated. NumPy’s pinned lstsq call uses rcond=None explicitly and returns rank and singular values; its residual array can be empty, so compute a residual directly for a consistent diagnostic.

For large systems, iterative methods can avoid a dense factorisation. Gradient iteration for a quadratic with SPD Hessian H has error update e_(t+1)=(I−ηH)e_t. Each eigenmode multiplies by 1−ηh_i; a fixed rate contracts all modes exactly when 0<η<2/h_max. If h_min is very small, the safe rate can still produce slow progress in that mode. Conjugate gradients exploits SPD structure and in exact arithmetic terminates within the dimension; finite arithmetic, stopping tolerances, and preconditioning qualify that statement.

Preconditioning transforms a system so iterative progress is better balanced, while preserving how the original solution is recovered. For SPD problems, symmetric transformations maintain the structure needed by the solver. Feature standardisation in learning similarly changes coordinates and can reduce curvature imbalance, but penalties expressed in the new coordinates represent a different geometry unless transformed consistently. Measure condition and convergence in the declared coordinates. Do not confuse a changed estimator, such as ridge, with a purely algebraic solver repair.

Costs also depend on structure. Dense QR for n≥d typically requires order nd² work and nd storage for the design; a triangular solve costs order d². A dense square factorisation costs order d³, while reusing its factors for another right-hand side is cheaper. Sparse matrix-vector products can cost proportional to stored nonzeros, but sparse factorisation can introduce fill. Report the operation model and memory shape before comparing a theoretical count with wall time. Hardware and library choices affect constants and parallelism.

The sigmoid σ(z)=1/(1+exp(−z)) can overflow internally for a large negative z. Evaluate its positive branch directly, and for negative z use exp(z)/(1+exp(z)). Both branches have exponent arguments at most zero. Finite arithmetic can still round extremely small values to zero and extremely large probabilities to one. Returning a plausible endpoint does not mean its logarithm remains a faithful loss calculation. Compute log probabilities from logits when that is the quantity required by the objective.

For finite logits z_j, let m=max_j z_j. Softmax equals exp(z_j−m)/Σ exp(z_k−m), since the common exp(m) factor cancels. Each exponential is at most one and at least one equals one. Log-sum-exp is m+log Σ exp(z_k−m). The class-c negative log-probability can be evaluated more accurately as log Σ exp(z_k−m)−(z_c−m), avoiding subtraction of a large common offset after adding it back. This derivation assumes a nonempty set of finite logits; all −infinity or NaN requires an explicit separate contract.

For Bernoulli y∈{0,1}, negative log-likelihood is softplus(z)−yz. A stable form is max(z,0)−yz+log1p(exp(−|z|)). The function log1p accurately evaluates log(1+t) for small t without first erasing t in 1+t. Its derivative is σ(z)−y. For a mean loss with X of shape n by d and w of length d, the gradient is Xᵀ(σ(Xw)−y)/n, plus any declared penalty gradient. Stable arithmetic must preserve the same objective and mean normalisation used in the derivation.

Gradient checking compares an analytic derivative with an independent numerical approximation of the same smooth scalar objective. A central difference in coordinate j is [F(w+he_j)−F(w−he_j)]/(2h). Taylor expansion with a bounded third derivative gives truncation error order h². Rounding in the two function evaluations is divided by h and typically contributes order u/h times a relevant function scale. Decreasing h indefinitely therefore makes the check worse. A rough balance suggests h of order u^(1/3) in suitably scaled coordinates, not a universal fixed value.

Use a sweep of h, both absolute and scale-aware relative errors, and several coordinates or random directions. A directional check compares gᵀv with a central difference along v and can be cheap for high-dimensional models. Scale v and h explicitly. Verify input shapes, objective reduction, regularisation, and intercept treatment before blaming rounding. A missing factor n produces a systematic discrepancy across a sensible h range; increasing precision will not fix it. Two implementations sharing the same wrong formula can agree, so independence of the check matters.

Coordinate steps also carry units. Adding the same h to a coefficient near a million and another near a billionth probes very different relative changes, and one perturbation may not alter its stored value at all. A practical check can start with h_j=h_base max(1,|w_j|) in declared nondimensional coordinates, then sweep h_base and inspect each result. The maximum is a heuristic scale choice, not a theorem ensuring accuracy for every objective. Report the actual perturbations and a useful error plateau. A single pass/fail threshold hides whether a disagreement tracks truncation, rounding, or a persistent algebraic mistake.

At a kink, such as absolute value at zero, the symmetric finite difference can be zero while no ordinary derivative exists. At piecewise branches, check away from boundaries unless a subgradient contract is specifically intended. Freeze stochastic sampling, dropout masks, and data order during a smooth derivative test; otherwise different function evaluations represent different random objectives. A check on a small controlled example establishes local implementation evidence, rather than a proof of every branch, every input, or optimiser convergence.

A numerical report should identify the equation, stored inputs, dtype at each sensitive stage, solver or update rule, and diagnostic tolerances. Include conditioning or singular values when interpreting coefficients and distinguish exact identities from measured floating outputs. A rounded display with six digits is insufficient evidence for a claimed twelve-digit error bound. Retain a machine-readable result or reproducible script, and label values that are specific to the tested platform. The mathematical explanation should account for a failure before the repaired result is presented.

Lower precision can reduce memory and increase throughput, but its value depends on an operation’s arithmetic intensity, bandwidth, supported hardware instructions, and accumulation precision. Array storage uses n×d×bytes_per_element before accounting for intermediates and factorisation workspaces. A conversion can temporarily keep two arrays alive. Computing a dense d by d Gram matrix may be a memory bottleneck even when storing the original data is feasible. Numerical quality and runtime should be compared under the same problem, stopping criteria, and budget.

Approximation error has several sources with different remedies. Discretisation or finite iteration error can decrease with refinement or more steps; rounding may increase as that extra work accumulates. Monte Carlo error decreases according to the sampling law from Module 24, whereas deterministic floating error depends on arithmetic and conditioning. Statistical estimation uncertainty comes from the observed data. Reporting all of these as a single “precision” number hides their meanings. State which error estimate is measured, bounded, or merely diagnosed, and which assumptions justify it.

Seeds identify a pseudorandom experiment only together with the generator, draw order, packages, and input generation. A fixed seed does not control every threaded numerical reduction or make samples independent. Different library builds can choose different algorithms and produce small final-bit changes; different generator interfaces can produce different datasets from the same integer seed. Use tolerance-based regression checks for floating quantities, exact checks for contracts, and structural mathematical checks such as probability sums, rank assumptions, and residuals. Record why each tolerance is reasonable.

Record package and Python versions, shapes, split rules, random generator, solver thresholds, and whether measured time includes data generation or factorisation. Separate compilation or warm-up from repeated timing when relevant, and report a distribution of measured times rather than an unsupported speed claim from a single run. For teaching labs, deterministic small CPU examples make the derivation inspectable. They are not benchmark evidence for all devices, all precisions, or the fastest deployment design.

Lab A · rounding, reduction, and rationalisation

Use the pinned NumPy preparation. Predict each format’s loss of a unit at its large-spacing boundary. Compare sequential sum, fsum and a Decimal reference for the actual binary inputs. Explain the rationalised expression; do not call display rounding a mathematical repair.

Download lab1_rounding_summation_and_cancellation.py

"""CPU rounding experiments; Decimal references the actual binary inputs."""
import math
from decimal import Decimal, localcontext
import numpy as np

print('NumPy', np.__version__, '; rounding to nearest on the tested CPU')
for dtype in (np.float16, np.float32, np.float64):
    info = np.finfo(dtype)
    large = dtype(2 ** (info.nmant + 1))
    terms = [large, dtype(1), -large]
    total = dtype(0)
    for term in terms:
        total = dtype(total + term)
    print(f'{dtype.__name__}: eps={info.eps:.9e}; spacing(large)={np.spacing(large):g}; sequential={total:g}; exact=1')
    assert total == 0
terms = [1e16, 1., -1e16]
with localcontext() as ctx:
    ctx.prec = 80
    exact = sum(map(Decimal.from_float, terms))
print('binary64 sequential:', sum(terms), '; fsum:', math.fsum(terms), '; binary-input Decimal reference:', exact)
assert math.fsum(terms) == float(exact) == 1
x = 1e-16
naive = math.sqrt(1+x)-1
stable = x/(math.sqrt(1+x)+1)
print(f'sqrt(1+x)-1 at x={x:g}: naive={naive:.9e}; rationalised={stable:.9e}')
assert naive == 0 and stable > 0
print('0.1+0.2:', repr(.1+.2), '; difference from float(0.3):', (.1+.2)-.3)
print('Subnormal != normal minimum:', np.finfo(np.float64).smallest_subnormal < np.finfo(np.float64).tiny)

NumPy 1.26.4 ; rounding to nearest on the tested CPU
float16: eps=9.765625000e-04; spacing(large)=2; sequential=0; exact=1
float32: eps=1.192092896e-07; spacing(large)=2; sequential=0; exact=1
float64: eps=2.220446049e-16; spacing(large)=2; sequential=0; exact=1
binary64 sequential: 0.0 ; fsum: 1.0 ; binary-input Decimal reference: 1
sqrt(1+x)-1 at x=1e-16: naive=0.000000000e+00; rationalised=5.000000000e-17
0.1+0.2: 0.30000000000000004 ; difference from float(0.3): 5.551115123125783e-17
Subnormal != normal minimum: True

Lab B · sensitivity, residuals, and solvers

Calculate the small right-hand-side perturbation and its amplified solution change. Compare QR, lstsq and normal equations on exactly generated least-squares data. Numerical coefficient errors can vary by library; the full-rank conditioning identity and small-residual counterexample are the mathematical lessons.

Download lab2_solves_conditioning_and_residuals.py

"""Small deterministic systems: residual, sensitivity and solver choice."""
import numpy as np

epsilon = 1e-8
A = np.array([[1., 1.], [1., 1.+epsilon]])
truth = np.ones(2)
b = A @ truth
delta = np.array([0., 1e-8])
x = np.linalg.solve(A, b+delta)
relative_input = np.linalg.norm(delta)/np.linalg.norm(b)
relative_error = np.linalg.norm(x-truth)/np.linalg.norm(truth)
print(f'cond2(A)={np.linalg.cond(A):.6e}; relative RHS perturbation={relative_input:.6e}; relative solution change={relative_error:.6e}')
wrong = np.array([0., 2.])
print(f'wrong x relative residual={np.linalg.norm(A@wrong-b)/np.linalg.norm(b):.6e}; relative error={np.linalg.norm(wrong-truth)/np.linalg.norm(truth):.6e}')
assert relative_input < 1e-8 and relative_error > .9
t = np.linspace(-1, 1, 30)
X = np.column_stack([np.ones_like(t), t, t+1e-7*t*t])
beta = np.array([1., 2., -1.])
y = X @ beta
Q, R = np.linalg.qr(X, mode='reduced')
qr = np.linalg.solve(R, Q.T@y)
svd, _, rank, singular = np.linalg.lstsq(X, y, rcond=None)
print(f'X rank={rank}; cond2(X)={singular[0]/singular[-1]:.6e}; computed cond2(Gram)={np.linalg.cond(X.T@X):.6e}')
for label, estimate in [('QR', qr), ('lstsq', svd)]:
    print(f'{label}: coefficient error={np.linalg.norm(estimate-beta):.6e}; residual={np.linalg.norm(X@estimate-y):.6e}')
    assert np.linalg.norm(estimate-beta) < 1e-6
try:
    normal = np.linalg.solve(X.T@X, X.T@y)
    print(f'normal equations: coefficient error={np.linalg.norm(normal-beta):.6e}; residual={np.linalg.norm(X@normal-y):.6e}')
except np.linalg.LinAlgError:
    print('normal equations: rounded Gram reported singular')
print('Exact full-column-rank identity: cond2(Gram)=cond2(X)^2; computed Gram may have lost its smallest mode.')
H = np.diag([1., 100.])
for rate in (.01, .03):
    w = np.ones(2)
    for _ in range(20):
        w -= rate * (H@w)
    print(f'quadratic eta={rate:.2f}; w20={w}; stable={0 < rate < 2/100}')

cond2(A)=4.000000e+08; relative RHS perturbation=3.535534e-09; relative solution change=1.000000e+00
wrong x relative residual=3.535534e-09; relative error=1.000000e+00
X rank=3; cond2(X)=4.444812e+07; computed cond2(Gram)=2.444580e+15
QR: coefficient error=3.437643e-09; residual=1.464485e-15
lstsq: coefficient error=2.904832e-09; residual=6.466036e-15
normal equations: coefficient error=7.856742e-02; residual=9.681657e-09
Exact full-column-rank identity: cond2(Gram)=cond2(X)^2; computed Gram may have lost its smallest mode.
quadratic eta=0.01; w20=[0.81790694 0.        ]; stable=True
quadratic eta=0.03; w20=[5.43794343e-01 1.04857600e+06]; stable=False

Lab C · stable loss and a derivative sweep

Repair raw exponential overflow and compare three central-difference steps on the identical mean objective. The smallest step fails through rounding. Diagnose the missing mean factor separately and explain the nondifferentiable absolute-value example.

Download lab3_stable_losses_and_gradient_checks.py

"""Repair overflow and audit a smooth central-difference gradient check."""
import numpy as np

def sigmoid(z):
    z = np.asarray(z, dtype=np.float64)
    out = np.empty_like(z)
    positive = z >= 0
    out[positive] = 1/(1+np.exp(-z[positive]))
    e = np.exp(z[~positive])
    out[~positive] = e/(1+e)
    return out

logits = np.array([1000., 1001., 1002.])
with np.errstate(over='ignore', invalid='ignore'):
    broken = np.exp(logits)/np.exp(logits).sum()
shifted = logits-logits.max()
softmax = np.exp(shifted)/np.exp(shifted).sum()
loss = np.log(np.exp(shifted).sum())-shifted[0]
print('raw softmax finite:', bool(np.isfinite(broken).all()))
print('shifted softmax:', softmax, '; class-0 NLL:', f'{loss:.9f}')
assert np.isfinite(softmax).all() and np.isclose(softmax.sum(), 1)
extreme = np.array([-1000., 1000.])
labels = np.array([1., 0.])
losses = np.maximum(extreme, 0)-labels*extreme+np.log1p(np.exp(-np.abs(extreme)))
print('stable sigmoid:', sigmoid(extreme), '; stable Bernoulli losses:', losses)
assert np.all(losses == 1000)
X = np.array([[1., .2], [1., -1.], [1., 2.]])
y = np.array([1., 0., 1.])
w = np.array([.3, -.4])
def objective(v):
    z = X@v
    return np.mean(np.logaddexp(0, z)-y*z)
analytic = X.T@(sigmoid(X@w)-y)/len(y)
for step in (1e-2, 1e-5, 1e-16):
    numeric = np.array([(objective(w+step*np.eye(2)[j])-objective(w-step*np.eye(2)[j]))/(2*step) for j in range(2)])
    error = np.linalg.norm(numeric-analytic)
    print(f'h={step:.0e}; numerical={numeric}; error={error:.6e}')
    if step == 1e-5:
        assert error < 1e-9
print('missing mean factor discrepancy:', np.linalg.norm(3*analytic-analytic))
print('abs at zero: central difference=0 while derivative does not exist; a symmetric slope is not a differentiability certificate.')

raw softmax finite: False
shifted softmax: [0.09003057 0.24472847 0.66524096] ; class-0 NLL: 2.407605964
stable sigmoid: [0. 1.] ; stable Bernoulli losses: [1000. 1000.]
h=1e-02; numerical=[-0.13316435 -0.66738056]; error=2.981183e-06
h=1e-05; numerical=[-0.13316411 -0.66738353]; error=1.168331e-11
h=1e-16; numerical=[0. 0.]; error=6.805391e-01
missing mean factor discrepancy: 1.3610781824396445
abs at zero: central difference=0 while derivative does not exist; a symmetric slope is not a differentiability certificate.

Check the pinned official NumPy 1.26 finfo, lstsq, and cond references for precision, rank thresholds and norm selection. The SciPy 1.14 logsumexp reference is an optional library comparison; the labs need NumPy alone. The conditioning and derivative arguments above are derived here.

Exit task: produce a numerical audit containing a spacing example, a cancellation repair, a residual/sensitivity counterexample, a solver justification, a stable probability calculation and a smooth derivative sweep.

Ready to move on: explain what the computation establishes and what uncertainty remains. Generalisation next connects mathematical risk statements to finite data and kernel models. Return to the course overview.