Robust Estimation

Robust estimation fits a model to measurements of which some are not noisy versions of the truth but simply wrong: mismatched features, pixels of a moving object when the camera’s motion is wanted, occluded or specular pixels in an alignment, a wrong loop closure in a map. Least squares gives every such gross error an influence that grows with its size, so a few of them can carry the estimate anywhere. A robust estimator limits that influence. Two families are used in vision: M-estimators, which replace the squared residual with a function that grows more slowly and are solved by iteratively reweighted least squares, and high-breakdown methods such as least median of squares and RANSAC, which search over fits to small random subsets. They can be combined, as in the example below: a high-breakdown search finds the right basin, and an M-estimator refines the answer inside it.

This article covers the framework that the applications share: what an outlier model is, how the influence of a residual is bounded or made to vanish, how the scale of the residuals is estimated, how tuning constants trade efficiency against robustness, what the breakdown point measures, and why the starting point matters. Nonlinear least squares introduces M-estimator costs and reweighting inside Gauss–Newton, L1 and L2 loss covers robust losses for training networks, and the application articles linked below cover their own robust choices.

Purpose

Least squares is the maximum likelihood estimate under independent Gaussian errors of equal variance. Murphy (2022, §11.6) puts the problem with outliers this way: the squared penalty grows quadratically, so points far from the fit pull on it more than points near it, and one remedy is a noise model with heavy tails, which assigns outliers a higher likelihood so that the fit need not move to explain them. Robust estimation in the wider sense keeps a model for the good data, such as a Gaussian, and asks for an estimator whose answer stays close to the right one when a fraction of the data follows some other, unknown process.

Vision produces such data in many places, and the outliers can be structured: not independent random errors but a coherent second population. The pixels of a car moving through a scene obey the car’s motion, not the camera’s; the wrong matches on a repeated facade agree with one another; a wrong loop closure on a second identical floor of a garage, the case Pose Graph Optimization simulates, is a perfectly consistent measurement of the wrong thing. Structured outliers are harder than scattered ones, because they can pull an estimate consistently in one direction or form a rival model of their own.

When to Use It

A fraction of the data is grossly wrong and cannot be removed beforehand: feature correspondences after matching, a flow field with moving objects, a pose graph with doubtful loop closures, photometric residuals with occlusions.
A good starting point exists or can be found: redescending M-estimators are local methods, so they suit refinement after a linear estimate or a random-sampling search.
The fraction of outliers matters for the choice. An M-estimator started from least squares tolerates few outliers, especially clustered ones, while a high-breakdown start such as least median of squares handles up to half the data; RANSAC can go beyond half when the outliers do not agree with one another, as the RANSAC article explains.

When nearly all data are good and the noise is Gaussian, plain least squares is the efficient estimator, and a robust estimator costs a little precision, quantified below. When the “outliers” are a second model of interest, robust estimation of one model is the first step of motion segmentation or multi-model fitting rather than a substitute for it.

How It Works

M-estimators and the influence of a residual

An M-estimator of parameters $\boldsymbol{\theta}$ minimizes

\sum_{i=1}^{m} \rho\!\left(\frac{r_i(\boldsymbol{\theta})}{s}\right),

where $r_i$ are the residuals, $s$ is a scale in the residuals’ units and $\rho$ is a symmetric function that grows more slowly than the square. Nonlinear Least Squares and Optimization use the same name. With $\rho = -\log p$ for a noise density $p$ , the estimate is the maximum likelihood estimate under that density. Setting the gradient to zero gives the estimating equations $\sum_i \psi(r_i/s)\, \partial r_i/\partial\boldsymbol{\theta} = 0$ with $\psi = \rho'$ . The function $\psi$ is the influence of one residual on the estimate, as L1 and L2 Loss describes it: for least squares $\psi(u) = u$ grows without bound, so one gross error has unbounded influence.

Three shapes of $\psi$ matter:

Unbounded: least squares, $\rho(u) = \tfrac{1}{2}u^2$ , $\psi(u) = u$ .
Monotone and bounded: the absolute value, $\psi(u) = \operatorname{sign}(u)$ , and Huber’s function, quadratic for $|u| \le k$ and linear beyond, with $\psi(u) = \operatorname{clip}(u, -k, k)$ . Every outlier keeps a constant pull of size $k$ , however large it is. The cost is convex, so for a linear model it has no spurious minima. The function is named after Huber’s 1964 work on robust estimation of a location parameter (not read here), as L1 and L2 Loss notes; the form is Murphy’s (2022, Eq. 5.36).
Redescending: $\psi$ rises, then falls back to zero, so a very large residual has no influence at all. Tukey’s biweight, which statsmodels’ TukeyBiweight notes is also called the bisquare, is one example:

\rho(u) = \frac{c^2}{6}\left[1 - \left(1 - \frac{u^2}{c^2}\right)^3\right], \quad \psi(u) = u\left(1 - \frac{u^2}{c^2}\right)^2, \quad w(u) = \left(1 - \frac{u^2}{c^2}\right)^2

for $|u| \le c$ , and $\rho = c^2/6$ , $\psi = w = 0$ beyond (the form in statsmodels 0.14.5’s TukeyBiweight). The Cauchy (Lorentzian) and Geman–McClure losses also redescend; L1 and L2 Loss gives them as members of Barron’s general loss. A redescending cost is not convex, and its minimum depends on where the iteration starts.

Iteratively reweighted least squares

Writing $\psi(u) = w(u)\,u$ turns the estimating equations into those of a weighted least-squares problem with weights $w_i = w(r_i/s)$ . Iteratively reweighted least squares (IRLS) alternates: compute the residuals and weights at the current estimate, solve the weighted least-squares problem (or take a weighted Gauss–Newton step for a nonlinear model), and repeat. For Huber this is the weight $\min(1, k/|u|)$ that Nonlinear Least Squares derives; for Tukey, the weight above, which is exactly zero beyond $c$ . Nonlinear Least Squares also gives Triggs and colleagues’ remark that reweighting can converge more slowly than a step that also models how the weights change. Libraries differ here:

SciPy 1.16.3’s least_squares evaluates its losses on $z = (f/C)^2$ , with $C$ = f_scale and a total cost $\tfrac{1}{2}C^2\sum\rho(z)$ ; its 'huber' loss, $\rho(z) = z$ for $z \le 1$ and $2\sqrt{z} - 1$ beyond, read on its own is twice the textbook form, but with these factors it is exactly the textbook Huber cost with $k = C$ , as the example checks numerically. Its robust losses rescale the Jacobian by $\sqrt{\rho' + 2\rho'' f^2}$ (scale_for_robust_loss_function in _lsq/common.py), a correction beyond plain IRLS weights, and method='lm' accepts only the plain squared loss.
Ceres Solver 2.2.0 applies its LossFunction to the squared norm $s = \lVert f_i\rVert^2$ of a whole residual block, with cost $\tfrac{1}{2}\rho(s)$ and a scale $a$ entering as $a^2\rho(s/a^2)$ . Its TukeyLoss, $\rho(s) = \tfrac{1}{3}(1 - (1-s)^3)$ for $s \le 1$ , is then the biweight above applied to the unscaled residual with $c = a$ , in the residual’s units (our derivation), acting on the vector norm rather than on each component. Its documentation describes the Triggs correction of the residual and Jacobian, and notes, for LossFunctionWrapper, that starting with a large scale and then reducing it can converge better than using a small scale from the start.
statsmodels’ RLM.fit is a textbook IRLS: by default an ordinary least-squares start, the MAD scale re-estimated after every weighted fit, at most 50 iterations.

Scale

The weights depend on $r/s$ , so an M-estimate of regression needs a scale that is itself robust. The standard deviation of the residuals is not: one gross error inflates it without bound. The median absolute deviation (MAD), divided by $\Phi^{-1}(3/4) \approx 0.6745$ so that it estimates $\sigma$ for Gaussian data, is one such scale; statsmodels’ mad uses this constant, and its RLM, whose default scale estimate it is, computes the MAD of the residuals about zero. The scale may be fixed from a robust start, or re-estimated at each iteration (statsmodels’ default). A third option estimates $\sigma$ jointly with the parameters: scikit-learn’s HuberRegressor minimizes a Huber criterion in which the coefficients, intercept and scale “are parameters to be optimized” (its docstring), citing Huber and Ronchetti’s book on concomitant scale estimates and a 2006 report by Owen (neither read here). Each choice gives a different estimate on the same data, as the example’s cross-check prints.

Efficiency and the tuning constants

Robustness costs precision when the data are in fact Gaussian. For an M-estimator with known scale, a standard argument for estimating equations (our derivation) gives the estimate’s asymptotic variance as $E[\psi^2]/(E[\psi'])^2$ times that of least squares, so the asymptotic efficiency at the Gaussian, relative to least squares, is

\mathrm{eff} = \frac{\big(E[\psi'(Z)]\big)^2}{E[\psi(Z)^2]}, \qquad Z \sim \mathcal{N}(0, 1).

For Huber’s $\psi$ this has a closed form (our derivation), $\mathrm{eff} = P^2/\big(P - 2k\varphi(k) + 2k^2(1 - \Phi(k))\big)$ with $P = 2\Phi(k) - 1$ . The tuning constants $k = 1.345$ (Huber) and $c = 4.685$ (Tukey) are statsmodels’ defaults, in units of the scale; the example computes both efficiencies by quadrature and prints 0.950000 and 0.949997, so each constant gives 95% efficiency at the Gaussian. scikit-learn’s HuberRegressor defaults to epsilon=1.35, on its own jointly estimated scale. A smaller constant is less efficient (our derivation) and, as scikit-learn’s docstring says of epsilon, more robust to outliers. Murphy (2022, §11.6.3) notes that the Huber threshold is usually set by hand or by cross-validation.

Breakdown point

The breakdown point of an estimator, in the sense used here, is the smallest fraction of the data that, replaced by arbitrary values, can carry the estimate arbitrarily far. As Median Filter states, the median’s is 50% and the mean’s is zero; by the same argument, one arbitrarily large residual moves a least-squares fit arbitrarily far. For regression, a bounded $\psi$ is not enough in our reading: clustered outliers all pull in the same direction, and outliers at extreme positions in the design (leverage points) pull the fit through the slope parameters, so a Huber estimate can be carried far from the right answer. The example shows Huber estimates more than a pixel from the background motion once a clustered object covers 0.203 of the field at the image edge, and 0.4 px from it at 0.297 in the centre. Least median of squares (LMedS) minimizes the median of the squared residuals instead of their sum; it fits the model to many small random subsets and keeps the one with the smallest median residual, as Motion Segmentation describes. RANSAC notes that it fails once more than half the data are outliers, because the median residual is then determined by the outliers. The example also shows LMedS’s other side: its estimate, fitted to three points, is far less precise than an M-estimate.

Starting points

A redescending M-estimator judges each residual at the current estimate. Started from least squares that the outliers have already moved, it may count outliers as inliers and inliers as outliers, and converge to the wrong fit; started inside the right basin, it ignores the outliers completely. Motion Segmentation reports Szeliski’s caution that with too many outliers M-estimators may not reach the global optimum, and Nonlinear Least Squares and Pose Graph Optimization make the same point for their costs. Hence the standard pipelines: a linear or random-sampling estimate first, then robust refinement, and graduated non-convexity, which Pose Graph Optimization describes, when no good start exists.

How to Apply It

Example: dominant affine motion with a moving object

Estimating the camera’s motion from a flow field in which an object moves on its own is the dominant-motion problem of Motion Segmentation; here it serves to compare estimators. The script builds a $128 \times 128$ flow field. The background moves by a known affine motion, $u = 0.5 + \tilde{x} - 0.5\tilde{y}$ and $v = -0.5 + 0.5\tilde{x} + \tilde{y}$ pixels, in normalized coordinates $\tilde{x} = (x - 63.5)/63.5$ and $\tilde{y} = (y - 63.5)/63.5$ , with pixel centres at integer column $x$ and row $y$ ( $y$ down). A full-height band, the object, translates by $(-3, +3)$ pixels, so every object pixel lies at least 2 pixels, 20 noise standard deviations, from the background model in each component. Both components carry Gaussian noise of $\sigma = 0.1$ px. The band is $2\,\mathrm{round}(64p)$ columns wide for $p = 0, 0.1, \ldots, 0.6$ , which gives object fractions of 0.094 to 0.594, and it sits either at the left edge of the image or in its centre, since the position decides how much leverage the outliers have. Each flow component is a separate linear regression on $[1, \tilde{x}, \tilde{y}]$ with 16,384 rows.

Six estimators run on each field: ordinary least squares; Huber IRLS with $k = 1.345$ ; least median of squares from 500 random three-pixel fits; and Tukey biweight IRLS with $c = 4.685$ started from each of the other three. Every IRLS fit re-estimates the MAD scale about zero at each iteration and stops when no parameter changes by more than $10^{-10}$ . The measure is the RMS error of the fitted motion over the background pixels, against the noise-free background motion, over seeds 0–9; flags are residuals beyond $k s$ or $c s$ at the final fit. Before the sweep, the script checks the generated data, checks that every estimator is exact on noise-free data, compares its Huber fit with SciPy’s least_squares and with scikit-learn’s HuberRegressor on one case, and computes the asymptotic efficiencies and an empirical efficiency over seeds 0–19 without an object. It uses NumPy 2.3.5, SciPy 1.16.3 and scikit-learn 1.7.2 and runs in about 30 seconds on a laptop.

import os
for var in ("OPENBLAS_NUM_THREADS", "OMP_NUM_THREADS", "MKL_NUM_THREADS"):
    os.environ[var] = "1"                      # one BLAS thread, set before NumPy loads
import warnings
import numpy as np
from scipy import integrate, stats
from scipy.optimize import least_squares
from sklearn.linear_model import HuberRegressor

# A 128 x 128 flow field: the background moves by a known affine motion, a band-shaped
# object by its own translation. Pixel centres at integer (x, y), y down; the motion is
# written in normalized coordinates xt, yt in [-1, 1].
H = W = 128
N = H * W
yy, xx = np.mgrid[0:H, 0:W].astype(float)
X = np.column_stack([np.ones(N), ((xx - 63.5) / 63.5).ravel(), ((yy - 63.5) / 63.5).ravel()])
B_BG = np.array([[0.5, 1.0, -0.5],     # u = 0.5 + 1.0 xt - 0.5 yt
                 [-0.5, 0.5, 1.0]])    # v = -0.5 + 0.5 xt + 1.0 yt
B_OBJ = np.array([[-3.0, 0.0, 0.0],    # the object translates by (-3, +3) px
                  [3.0, 0.0, 0.0]])
BG, OBJ = X @ B_BG.T, X @ B_OBJ.T      # (N, 2) noise-free motions
SIGMA, K, C = 0.1, 1.345, 4.685         # noise (px); Huber and Tukey tuning constants
PS = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6]
MAD_C = stats.norm.ppf(0.75)            # MAD / 0.6745 estimates a Gaussian sigma

def width(p):                           # band width in columns, even so it can be centred
    return 2 * int(round(64 * p))

def object_mask(p, place):
    m = np.zeros((H, W), bool)
    c0 = 0 if place == "edge" else 64 - width(p) // 2
    m[:, c0:c0 + width(p)] = True
    return m.ravel()

def make_flow(p, place, seed, noise=True):
    m = object_mask(p, place)
    f = np.where(m[:, None], OBJ, BG)
    if noise:                           # the same noise for every p and placement of a seed
        f = f + SIGMA * np.random.default_rng(seed).standard_normal((2, N)).T
    return f, m

def bg_rms(B, m):                       # RMS error of the fitted motion on background pixels
    d = X[~m] @ B.T - BG[~m]
    return np.sqrt(np.mean(np.sum(d ** 2, axis=1)))

# ---- estimators; each flow component u, v is its own regression y = X beta -------------
def ols(y):
    return np.linalg.lstsq(X, y, rcond=None)[0]

def mad_scale(r):                       # MAD about zero; floor only matters for exact data
    return max(np.median(np.abs(r)) / MAD_C, 1e-12)

def huber_w(z):
    a = np.abs(z)
    return np.where(a <= K, 1.0, K / np.maximum(a, K))

def tukey_w(z):
    return np.where(np.abs(z) <= C, (1 - (z / C) ** 2) ** 2, 0.0)

CAPPED = []                             # IRLS fits that reach the iteration cap

def irls(y, beta, weight, tuning, maxit=500):
    for _ in range(maxit):
        r = y - X @ beta
        sw = np.sqrt(weight(r / mad_scale(r)))       # scale re-estimated every iteration
        new = np.linalg.lstsq(X * sw[:, None], y * sw, rcond=None)[0]
        step, beta = np.max(np.abs(new - beta)), new
        if step < 1e-10:
            break
    else:
        CAPPED.append(1)
    r = y - X @ beta
    s = mad_scale(r)
    return beta, s, np.abs(r) > tuning * s           # parameters, scale, flagged residuals

def lmeds(y, samples):                  # least median of squares over random 3-pixel fits
    A = X[samples]
    keep = np.abs(np.linalg.det(A)) > 1e-6
    betas = np.linalg.solve(A[keep], y[samples[keep]][..., None])[..., 0]
    med = np.median((y[:, None] - X @ betas.T) ** 2, axis=0)
    return betas[np.argmin(med)]

NAMES = ["OLS", "Huber", "LMedS", "Tukey<-OLS", "Tukey<-Huber", "Tukey<-LMedS"]

def fit_all(f, seed):
    rng = np.random.default_rng(10_000 + seed)
    samples = np.array([rng.choice(N, 3, replace=False) for _ in range(500)])
    B = {n: np.zeros((2, 3)) for n in NAMES}
    flags = {n: np.zeros((N, 2), bool) for n in NAMES[3:] + ["Huber"]}
    scale = {}
    for j in range(2):
        y = f[:, j]
        B["OLS"][j] = ols(y)
        B["Huber"][j], scale[j], flags["Huber"][:, j] = irls(y, B["OLS"][j], huber_w, K)
        B["LMedS"][j] = lmeds(y, samples)
        for start in ("OLS", "Huber", "LMedS"):
            n = "Tukey<-" + start
            B[n][j], _, flags[n][:, j] = irls(y, B[start][j], tukey_w, C)
    return B, flags, scale

# ---- checks -----------------------------------------------------------------------------
ok = lambda c: "held" if c else "FAILED"
counts = [(int(object_mask(p, pl).sum()), 128 * width(p)) for p in PS for pl in ("edge", "centre")]
print(f"D1 object pixels = 128 x band width, every p and placement: {ok(all(a == b for a, b in counts))}")
print("   object fractions:", " ".join(f"{128 * width(p) / N:.4f}" for p in PS))
gap = np.abs(OBJ - BG).min(axis=0)
print(f"D2 min |object - background| motion: u {gap[0]:.6f}, v {gap[1]:.6f} px: "
      f"{ok(np.all(np.abs(gap - 2.0) < 1e-12))}")
sd = SIGMA * np.random.default_rng(0).standard_normal((2, N)).std()
print(f"D3 seed-0 noise SD {sd:.5f} px: {ok(abs(sd - SIGMA) < 0.002)}")

f, m = make_flow(0.0, "edge", 0, noise=False)
B, _, _ = fit_all(f, 0)
worst = max(bg_rms(B[n], m) for n in NAMES)
print(f"C1a noise-free, no object: largest background RMS over all estimators {worst:.1e} px: "
      f"{ok(worst < 1e-9)}")
f, m = make_flow(0.3, "edge", 0, noise=False)
B, _, _ = fit_all(f, 0)
e1, e2 = bg_rms(B["LMedS"], m), bg_rms(B["Tukey<-LMedS"], m)
print(f"C1b noise-free, edge band 0.297: LMedS {e1:.1e}, Tukey<-LMedS {e2:.1e} px: "
      f"{ok(max(e1, e2) < 1e-9)}")

f, m = make_flow(0.3, "edge", 0)                     # the cross-check case
B, _, scale = fit_all(f, 0)
for j, comp in enumerate("uv"):
    y = f[:, j]
    res = least_squares(lambda b: X @ b - y, B["OLS"][j], jac=lambda b: X, loss="huber",
                        f_scale=K * scale[j], method="trf", ftol=1e-14, xtol=1e-14, gtol=1e-14)
    r = np.abs(X @ res.x - y)
    kk = K * scale[j]
    textbook = np.sum(np.where(r <= kk, 0.5 * r ** 2, kk * r - 0.5 * kk ** 2))
    d_sp = np.max(np.abs(res.x - B["Huber"][j]))
    rel = abs(res.cost - textbook) / textbook
    print(f"C2 {comp}: our Huber IRLS vs SciPy least_squares(loss='huber', f_scale={kk:.4f}): "
          f"max param diff {d_sp:.1e}, SciPy cost vs textbook Huber sum {rel:.1e}: "
          f"{ok(d_sp < 1e-6 and rel < 1e-10)}")
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        hr = HuberRegressor(epsilon=K, alpha=0.0, fit_intercept=False,
                            max_iter=10_000, tol=1e-12).fit(X, y)
    print(f"C3 {comp}: HuberRegressor(alpha=0) param diff {np.max(np.abs(hr.coef_ - B['Huber'][j])):.1e}, "
          f"its scale {hr.scale_:.4f} vs MAD scale {scale[j]:.4f} px, {hr.n_iter_} L-BFGS iterations")

phi = stats.norm.pdf
def efficiency(psi, dpsi, cuts):        # (E psi')^2 / E psi^2 under N(0, 1), our derivation
    q = lambda g: sum(integrate.quad(lambda z: g(z) * phi(z), a, b, epsabs=1e-13)[0]
                      for a, b in zip(cuts[:-1], cuts[1:]))
    return q(dpsi) ** 2 / q(lambda z: psi(z) ** 2)
e_hub = efficiency(lambda z: np.clip(z, -K, K), lambda z: float(abs(z) <= K), [-12, -K, K, 12])
u = lambda z: (z / C) ** 2
e_tuk = efficiency(lambda z: z * (1 - u(z)) ** 2, lambda z: (1 - u(z)) * (1 - 5 * u(z)), [-C, C])
P = 2 * stats.norm.cdf(K) - 1
e_closed = P ** 2 / (P - 2 * K * phi(K) + 2 * K ** 2 * stats.norm.sf(K))
print(f"C4 Huber efficiency: quadrature {e_hub:.6f}, closed form {e_closed:.6f}: "
      f"{ok(abs(e_hub - e_closed) < 1e-10)}")
print(f"Asymptotic Gaussian efficiency: Huber k=1.345 {e_hub:.6f}, Tukey c=4.685 {e_tuk:.6f}")

rms0 = {n: [] for n in NAMES}
for seed in range(20):
    f, m = make_flow(0.0, "edge", seed)
    B, _, _ = fit_all(f, seed)
    for n in NAMES:
        rms0[n].append(bg_rms(B[n], m))
ms = {n: np.sum(np.square(v)) for n, v in rms0.items()}
print("Empirical efficiency, no object, seeds 0-19 (sum of OLS MSE / sum of MSE):",
      ", ".join(f"{n} {ms['OLS'] / ms[n]:.4f}" for n in NAMES[1:]))

# ---- the sweep over object fraction ---------------------------------------------------------
for place in ("edge", "centre"):
    rms = {n: np.zeros((len(PS), 10)) for n in NAMES}
    fobj = {n: np.zeros((len(PS), 10)) for n in NAMES[3:] + ["Huber"]}
    fbg = {n: np.zeros((len(PS), 10)) for n in NAMES[3:] + ["Huber"]}
    for i, p in enumerate(PS):
        for seed in range(10):
            f, m = make_flow(p, place, seed)
            B, flags, _ = fit_all(f, seed)
            for n in NAMES:
                rms[n][i, seed] = bg_rms(B[n], m)
            for n in fobj:
                fobj[n][i, seed] = flags[n][m].mean() if m.any() else np.nan
                fbg[n][i, seed] = flags[n][~m].mean()
    row = lambda label, vals, fmt: print(f"{label:24s}" + "".join(format(v, fmt) for v in vals))
    print(f"\nBackground RMS error of the fitted motion (px), {place} band, seeds 0-9")
    row("object fraction", [128 * width(p) / N for p in PS], "8.3f")
    for n in NAMES:
        for stat, fn in (("mean", np.mean), ("min", np.min), ("max", np.max)):
            row(f"{n if stat == 'mean' else '':16s}{stat}", fn(rms[n], axis=1), "8.4f")
    row("seeds Huber < OLS", np.sum(rms["Huber"] < rms["OLS"], axis=1), "8d")
    row("seeds LMedS > T<-LMedS", np.sum(rms["LMedS"] > rms["Tukey<-LMedS"], axis=1), "8d")
    print("flagged residuals, mean over seeds (object, background)")
    for n in fobj:
        row(f"{n:16s}obj", np.mean(fobj[n], axis=1), "8.4f")
        row(f"{'':16s}bg", np.mean(fbg[n], axis=1), "8.5f")
    row("Tukey<-LMedS    obj min", np.min(fobj["Tukey<-LMedS"], axis=1), "8.4f")
    row("                bg max", np.max(fbg["Tukey<-LMedS"], axis=1), "8.5f")
print(f"\nC5 IRLS fits that reached the 500-iteration cap: {len(CAPPED)}: {ok(len(CAPPED) == 0)}")

Output:

D1 object pixels = 128 x band width, every p and placement: held
   object fractions: 0.0000 0.0938 0.2031 0.2969 0.4062 0.5000 0.5938
D2 min |object - background| motion: u 2.000000, v 2.000000 px: held
D3 seed-0 noise SD 0.09978 px: held
C1a noise-free, no object: largest background RMS over all estimators 9.3e-14 px: held
C1b noise-free, edge band 0.297: LMedS 8.0e-17, Tukey<-LMedS 1.8e-14 px: held
C2 u: our Huber IRLS vs SciPy least_squares(loss='huber', f_scale=1.2696): max param diff 9.9e-12, SciPy cost vs textbook Huber sum 1.7e-16: held
C3 u: HuberRegressor(alpha=0) param diff 2.2e-02, its scale 0.6333 vs MAD scale 0.9440 px, 10 L-BFGS iterations
C2 v: our Huber IRLS vs SciPy least_squares(loss='huber', f_scale=1.7971): max param diff 3.2e-12, SciPy cost vs textbook Huber sum 0.0e+00: held
C3 v: HuberRegressor(alpha=0) param diff 6.9e-02, its scale 0.8791 vs MAD scale 1.3361 px, 8 L-BFGS iterations
C4 Huber efficiency: quadrature 0.950000, closed form 0.950000: held
Asymptotic Gaussian efficiency: Huber k=1.345 0.950000, Tukey c=4.685 0.949997
Empirical efficiency, no object, seeds 0-19 (sum of OLS MSE / sum of MSE): Huber 0.9035, LMedS 0.0072, Tukey<-OLS 0.8902, Tukey<-Huber 0.8902, Tukey<-LMedS 0.8902

Background RMS error of the fitted motion (px), edge band, seeds 0-9
object fraction            0.000   0.094   0.203   0.297   0.406   0.500   0.594
OLS             mean      0.0018  0.7142  1.1768  1.3329  1.3268  1.2462  1.2242
                min       0.0011  0.7135  1.1762  1.3324  1.3262  1.2456  1.2234
                max       0.0024  0.7152  1.1776  1.3335  1.3272  1.2468  1.2252
Huber           mean      0.0019  0.0587  1.1179  1.3395  1.3379  1.2409  1.1953
                min       0.0012  0.0578  1.1166  1.3388  1.3371  1.2402  1.1938
                max       0.0024  0.0598  1.1204  1.3401  1.3388  1.2414  1.1966
LMedS           mean      0.0203  0.0270  0.0299  0.0345  0.0379  1.3832  5.2655
                min       0.0120  0.0203  0.0225  0.0217  0.0177  1.2171  5.2174
                max       0.0312  0.0335  0.0411  0.0507  0.0578  2.3329  5.3210
Tukey<-OLS      mean      0.0019  0.0019  1.1255  1.3437  1.3370  1.2421  1.1950
                min       0.0012  0.0012  1.1246  1.3430  1.3364  1.2415  1.1941
                max       0.0024  0.0024  1.1267  1.3443  1.3378  1.2425  1.1959
Tukey<-Huber    mean      0.0019  0.0019  1.1255  1.3437  1.3370  1.2421  1.1950
                min       0.0012  0.0012  1.1246  1.3430  1.3364  1.2415  1.1941
                max       0.0024  0.0024  1.1267  1.3443  1.3378  1.2425  1.1959
Tukey<-LMedS    mean      0.0019  0.0019  0.0019  0.0019  0.0022  1.2421  5.2465
                min       0.0012  0.0012  0.0012  0.0012  0.0016  1.2415  5.2427
                max       0.0024  0.0024  0.0026  0.0026  0.0029  1.2425  5.2512
seeds Huber < OLS              0      10      10       0       0      10      10
seeds LMedS > T<-LMedS        10      10      10      10      10       8       8
flagged residuals, mean over seeds (object, background)
Tukey<-OLS      obj          nan  1.0000  0.0000  0.0000  0.0000  0.0000  0.0000
                bg       0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
Tukey<-Huber    obj          nan  1.0000  0.0000  0.0000  0.0000  0.0000  0.0000
                bg       0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
Tukey<-LMedS    obj          nan  1.0000  1.0000  1.0000  1.0000  0.0000  0.0000
                bg       0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 1.00000
Huber           obj          nan  1.0000  0.4679  0.1262  0.1519  0.1821  0.1798
                bg       0.17780 0.12976 0.03308 0.10332 0.16521 0.13940 0.09010
Tukey<-LMedS    obj min      nan  1.0000  1.0000  1.0000  1.0000  0.0000  0.0000
                bg max   0.00003 0.00000 0.00000 0.00000 0.00000 0.00000 1.00000

Background RMS error of the fitted motion (px), centre band, seeds 0-9
object fraction            0.000   0.094   0.203   0.297   0.406   0.500   0.594
OLS             mean      0.0018  0.4678  1.0139  1.4820  2.0287  2.4983  2.9701
                min       0.0011  0.4667  1.0128  1.4809  2.0275  2.4971  2.9689
                max       0.0024  0.4692  1.0153  1.4834  2.0301  2.4997  2.9715
Huber           mean      0.0019  0.0254  0.0778  0.4187  2.0287  2.4983  2.9701
                min       0.0012  0.0242  0.0761  0.4109  2.0275  2.4971  2.9689
                max       0.0024  0.0271  0.0792  0.4267  2.0301  2.4997  2.9715
LMedS           mean      0.0203  0.0238  0.0248  0.0269  0.0284  5.4759  5.0763
                min       0.0120  0.0121  0.0136  0.0154  0.0166  5.2376  5.0475
                max       0.0312  0.0352  0.0354  0.0357  0.0397  5.7363  5.1058
Tukey<-OLS      mean      0.0019  0.0020  0.0020  0.0023  1.9667  2.4963  3.0238
                min       0.0012  0.0016  0.0014  0.0018  1.9657  2.4952  3.0227
                max       0.0024  0.0024  0.0029  0.0034  1.9682  2.4978  3.0252
Tukey<-Huber    mean      0.0019  0.0020  0.0020  0.0023  1.9667  2.4963  3.0238
                min       0.0012  0.0016  0.0014  0.0018  1.9657  2.4952  3.0227
                max       0.0024  0.0024  0.0029  0.0034  1.9682  2.4978  3.0252
Tukey<-LMedS    mean      0.0019  0.0020  0.0020  0.0023  0.0026  2.4963  5.0742
                min       0.0012  0.0016  0.0014  0.0018  0.0017  2.4952  5.0727
                max       0.0024  0.0024  0.0029  0.0034  0.0036  2.4978  5.0757
seeds Huber < OLS              0      10      10      10       0       0       0
seeds LMedS > T<-LMedS        10      10      10      10      10      10       6
flagged residuals, mean over seeds (object, background)
Tukey<-OLS      obj          nan  1.0000  1.0000  1.0000  0.0000  0.0000  0.0000
                bg       0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
Tukey<-Huber    obj          nan  1.0000  1.0000  1.0000  0.0000  0.0000  0.0000
                bg       0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
Tukey<-LMedS    obj          nan  1.0000  1.0000  1.0000  1.0000  0.0000  0.0000
                bg       0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 1.00000
Huber           obj          nan  1.0000  1.0000  1.0000  0.0000  0.0000  0.0000
                bg       0.17780 0.12995 0.07021 0.00002 0.00000 0.00000 0.00000
Tukey<-LMedS    obj min      nan  1.0000  1.0000  1.0000  1.0000  0.0000  0.0000
                bg max   0.00003 0.00000 0.00000 0.00000 0.00000 0.00000 1.00000

C5 IRLS fits that reached the 500-iteration cap: 0: held

Checks. The data checks D1–D3, the exact-data checks C1a and C1b, the SciPy cross-check C2 for both components, the closed-form check C4 and the iteration-cap check C5 held. C2 confirms the convention above: started from least squares with f_scale equal to $k$ times our final MAD scale, SciPy’s minimizer agrees with our IRLS to $10^{-11}$ , and its reported cost equals the textbook Huber sum. HuberRegressor with alpha=0 (C3, printed without a bar) differs from our fit by 0.022 and 0.069 in the parameters, and its scale, which it estimates jointly with the coefficients, is 0.6333 and 0.8791 px against our MAD scales of 0.9440 and 1.3361 px: the same Huber function on a different scale gives a different estimate. Its default alpha=1e-4 adds a ridge penalty $\alpha\lVert\mathbf{w}\rVert^2$ to the coefficients (the source’s docstring), so alpha=0 is needed for a pure M-estimate.

Efficiency. Both asymptotic efficiencies lie within 0.005 of 0.95, so that hypothesis held; it is an illustration rather than a test, because a planning script had computed the same two values before the pilot. With no object, the empirical efficiency over 20 seeds, the summed squared error of least squares divided by that of the estimator, is 0.9035 for Huber and 0.8902 for each Tukey fit, within 0.1 of the asymptotic values, so that hypothesis held too. Least median of squares reaches 0.0072: its fit rests on three pixels, which is why it is used as a start rather than an answer. Huber flags 0.17780 of the background residuals at $p = 0$ , against $2(1 - \Phi(1.345)) = 0.17863$ for Gaussian residuals beyond 1.345 standard deviations (our derivation; held, within 0.01): Huber’s “outliers” include many good measurements, which still contribute with reduced weight.

Least squares. From the smallest object fraction on, every seed’s error exceeds 0.1 px at both placements (held), from 0.7135 px at the edge and 0.4667 px in the centre at a fraction of 0.094. The hypothesis that the mean error increases strictly with the fraction held for the centre band (0.4678 rising to 2.9701 px) and failed for the edge band, where the means are 0.7142, 1.1768, 1.3329, 1.3268, 1.2462 and 1.2242 px: they peak at 0.297 and then fall.

Huber. Huber’s error exceeds 0.01 px for every seed at every nonzero fraction at both placements (held: it is biased from the first outliers on, 0.0587 px at the edge and 0.0254 px in the centre at 0.094). The hypothesis that Huber stays below least squares for every seed failed in five conditions and held in the other seven. At the edge the counts of seeds with Huber below least squares are 10, 10, 0, 0, 10 and 10 over the six fractions, so it failed at 0.297 and 0.406, where Huber’s means are 1.3395 and 1.3379 px against least squares’ 1.3329 and 1.3268. In the centre the counts are 10, 10, 10, 0, 0 and 0, so it failed at 0.406, 0.500 and 0.594, where Huber’s means equal least squares’ (2.0287, 2.4983 and 2.9701 px). The hypothesis that Huber’s mean error increases strictly with the fraction failed at the edge (0.0587, 1.1179, 1.3395, 1.3379, 1.2409, 1.1953 px) and held in the centre (0.0254, 0.0778, 0.4187, 2.0287, 2.4983, 2.9701 px). The flag rows show where Huber stops separating the object: at the edge it flags all object residuals at 0.094 but 0.4679 of them at 0.203 and 0.1262 at 0.297; in the centre all of them up to 0.297 and none from 0.406 on, together with no background residuals. In the cross-check case (edge, 0.297, seed 0), the printed MAD scales at the Huber fit are 0.9440 and 1.3361 px, about nine and thirteen times the noise; these are printed values, and the example does not separate their part in the result. A bounded influence did not keep the estimate near the background motion once the clustered object was a large fraction of the field.

Starting points and breakdown. Tukey’s biweight started from least median of squares stays below 0.01 px for every seed at every fraction up to 0.406 at both placements, with a largest error of 0.0036 px, and every seed exceeds 0.1 px at 0.594 (5.2427 px at the edge and 5.0727 px in the centre at the best seed), where it fits the object: held. At 0.5, printed but not scored, its means are 1.2421 px at the edge and 2.4963 px in the centre, and least median of squares alone gives 1.3832 and 5.4759 px. The same biweight started from least squares fails at a smaller fraction, 0.203 at the edge and 0.406 in the centre, against 0.5 from least median of squares at both placements (held). Tukey started from Huber, which has no hypothesis, matched Tukey started from least squares in every printed row. The flags show the two outcomes cleanly: where the start lies in the right basin, the biweight flags every object residual and no background residual (from least median of squares, every object residual for every seed and at most 0.00003 of background residuals up to 0.406: held); where it does not, it either flags nothing, with an error of the same order as that of least squares, or, started from least median of squares at 0.594, flags every background residual and no object residual. At every fraction up to 0.406 and for every seed, least median of squares alone has a larger error than its Tukey refinement (held at both placements): from 0.0120 to 0.0578 px at the edge, against 0.0029 px or less after refinement.

In short, the high-breakdown start and the efficient refinement do different jobs, and neither does both here: least median of squares finds the background motion up to an object fraction of 0.406 here but locates it to a few hundredths of a pixel, and the biweight locates it to a few thousandths of a pixel only when started inside its basin.

Process. The script was run as planned. Its set-up, measures, checks and hypotheses, with every bar, were written before any script ran; one planning script then computed the reachable values (the error of least squares on the known background pixels, the error of a fit that follows the object, the two efficiencies, the MAD constant and the Gaussian flag fraction). Of the bars, only the numerical value 0.17863 of the Huber flag fraction was first written after that script printed it; the hypothesis named it as $2(1 - \Phi(1.345))$ before. The 95% efficiency bar came from the batch plan and was met by the planning script, as stated above. The pilot and the first run used the same script and seeds, and their outputs are identical. Twelve points of the set-up depart from the batch plan’s description of this example, each fixed before the pilot: (1) the object is a full-height band $2\,\mathrm{round}(64p)$ columns wide, so the fractions are 0.094–0.594 rather than exactly 0.1–0.6; (2) it is run at two placements, edge and centre; (3) the motions were fixed from the construction, the object as a pure translation 2 px or more from the background model; (4) u and v are separate regressions, with least median of squares sharing its samples between them and choosing per component; (5) seeds 0–9 for the sweep and 0–19 without an object; (6) least median of squares uses 500 samples and skips degenerate ones; (7) the scale is the MAD about zero, re-estimated each iteration and floored at $10^{-12}$ px for exact data, with a $10^{-10}$ stopping rule and a 500-iteration cap; (8) “near the noise level” became an error below 0.01 px and “fails” an error above 0.1 px, with the fraction 0.5 printed but not scored; (9) the empirical efficiency is the ratio of summed squared errors over seeds, printed for least median of squares without a bar; (10) the cross-check case and solver settings were fixed, with SciPy’s cost also checked against the textbook Huber sum; (11) background flag fractions, seed counts, data checks, exact-data checks, the closed-form check and the iteration-cap count were added; (12) the hypotheses on least median of squares alone and on the flag fractions, and the bar on the empirical efficiency, were added.

Variants

Random sampling and consensus. RANSAC maximizes the number of inliers within a threshold over models fitted to minimal samples; least median of squares minimizes the median residual over the same kind of samples. RANSAC needs a threshold; least median of squares needs none but stops at half the data. Both are usually followed by a refinement on the inliers, as in the homography pipeline.
Heavy-tailed likelihoods. Murphy (2022, §11.6) treats robust regression as maximum likelihood with a Laplace likelihood, whose estimate can be computed by linear programming (§11.6.1), or a Student-t likelihood (§11.6.2). Direct Methods describes Kerl, Sturm and Cremers’ (“Robust Odometry Estimation for RGB-D Cameras”, 2013) t-distribution weights for photometric residuals.
Robust kernels in graph optimization. Pose Graph Optimization covers Huber, Cauchy, dynamic covariance scaling, switchable constraints and graduated non-convexity for wrong loop closures, and Dellaert and Kaess’s report on which kernels the common back-ends offer and recommend; for bundle adjustment, Least Squares notes that Triggs and colleagues’ survey covers the robust costs often used in place of the square.
Robust penalties in dense estimation. Black and Anandan made both terms of variational optical flow robust, as Brightness Constancy and Horn–Schunck describe; robust Lucas–Kanade variants reweight occluded or specular pixels.
Robust losses in training. L1 and L2 Loss covers Huber, smooth L1, Charbonnier and Barron’s general loss as training losses, and Loss Function shows bounded influence on a burst of frames with impulse noise.
Mode seeking. Mean shift computes a robust location M-estimate, as Comaniciu and Meer show.

Best Practices

Start robustly, then refine efficiently. Use a high-breakdown method (random sampling, least median of squares) or a linear estimate on cleaned data for the start, and an M-estimator with an efficient tuning constant for the final fit.
Estimate the scale robustly, with the MAD (divided by 0.6745) or a scale from the high-breakdown start, and state whether it is fixed, re-estimated or estimated jointly.
State tuning constants in units of the scale (1.345 and 4.685 for 95% Gaussian efficiency) or in residual units, and say which; a constant in pixels is wrong for a residual in millimetres.
Look at the weights. The fraction of zero-weight residuals and where they lie (the object, the border, one image) tell whether the fit separated the outliers.
Check the result against the noise level, the cheap check Nonlinear Least Squares suggests: in our reading, a robust fit whose inlier residuals are much larger than the expected noise has probably converged to the wrong basin.

Pitfalls

Treating a bounded influence as a high breakdown point. In the example, Huber’s error exceeds that of least squares at some object fractions.
Starting a redescending estimator from least squares when the outliers are many or clustered: the biweight then keeps the wrong fit and flags nothing.
Mixing conventions. SciPy and Ceres apply $\rho$ to a squared, scaled residual with a factor $\tfrac{1}{2}$ outside; textbooks and statsmodels apply it to the residual divided by the scale. Ceres applies it to the norm of a residual block, statsmodels and scikit-learn to each scalar residual.
Library defaults. HuberRegressor adds a ridge penalty by default (alpha=1e-4) and estimates its own scale; SciPy’s method='lm' refuses robust losses.
Reading a flag as a verdict. Huber flags about 18% of perfectly Gaussian residuals at $k = 1.345$ , by construction.

Evidence and Impact

In our reading, robust estimation is what lets geometric vision use real matches: the standard pipelines of feature matching, homography and camera-pose estimation run random sampling to find the inliers and robust or plain least squares to refine them, and dense alignment, optical flow and SLAM back-ends weight their residuals robustly. The evidence for particular choices is application-specific and is reported in the articles that own them: Kerl, Sturm and Cremers’ comparison of weights for RGB-D odometry in Direct Methods, the loop-closure experiments of switchable constraints and dynamic covariance scaling in Pose Graph Optimization, Barron’s adaptive loss in L1 and L2 Loss. The example here shows the general pattern on one controlled problem: the 95% efficiency of the usual tuning constants costs little on clean data, a bounded influence does not survive a large clustered object, and a redescending estimator is as good as its starting point.