Robust Regression
Huber loss, M-estimators, RANSAC, and how far outliers can go before a fit breaks.
On this page
Beginner: ordinary least squares (2.3.1) fits a line by minimizing the sum of squared residuals — the vertical gaps between the line and each point, squared, then added up. Squaring is what makes a huge miss punished far more than a small one: a residual of 10 contributes 100 to the sum, a residual of 1 contributes only 1. That sounds harmless until one single data point is wildly wrong — a sensor glitch, a typo during data entry, a corrupted transmission — and its one huge squared residual can swamp every other point's contribution combined. OLS, trying to minimize the total, will happily bend the entire fitted line toward that one bad point to shave a little off its enormous squared error, dragging the line away from every well-behaved point in the process. Robust regression is the umbrella term for fitting methods built specifically so that a small number of badly-behaved points cannot do this.
Intermediate: the fix explored in this lesson swaps the loss function. The Huber loss behaves exactly like squared error for small residuals (so it fits clean, well-behaved data just as sensibly as OLS does) but switches to behaving like absolute error — a straight linear penalty instead of a quadratic one — once a residual passes a threshold δ. A residual of 100 under squared loss contributes 10,000; under Huber loss past the threshold it contributes only proportionally to 100, not 100². One badly-placed point simply cannot dominate the sum the way it can under pure squared error. Fitting under this hybrid loss turns out to reduce, as derived below, to repeatedly re-solving an ordinary weighted least-squares problem with weights that quietly shrink for points with large residuals — the same "compute-weights-then-resolve" loop used to fit GLMs (2.3.9), aimed at a different problem.
Advanced: Huber-style reweighting is not the only strategy, and it is not even the most outlier-tolerant one. RANSAC (Random Sample Consensus) throws away the idea of using every point's residual at all. Instead it repeatedly grabs a tiny random subset of points — as few as the model needs to be uniquely determined, just 2 for a line — fits only to that subset, then checks how many of the remaining points happen to agree with that fit. Across enough random trials, a subset drawn entirely from the clean, well-behaved majority will eventually get picked, and its fit will attract far more agreement ("inliers") than any fit contaminated by a bad point. RANSAC never averages over the bad points at all — it just needs to get lucky once, and running enough trials makes that near-certain. That different strategy is why RANSAC tolerates far higher fractions of corrupted data than reweighting-based methods like Huber regression can, at the cost of being a randomized, trial-and-error search rather than a single deterministic optimization.
Worked intuition: picture eleven points that fall almost exactly on a line with slope 1.2, plus one twelfth point that starts near the line and then drifts further and further below it — precisely the scenario animated in the diagram below. Under OLS, that twelfth point's squared residual grows as the square of how far it drifts, so the fit keeps tilting to chase it, and the recovered slope drifts away from 1.2 continuously, without ever settling. Under Huber loss, once that point's residual passes δ its contribution grows only linearly, and the IRLS loop derived below assigns it a weight that shrinks roughly in proportion to 1/|residual| — so however far it drifts, its influence on the fitted slope keeps shrinking rather than growing. The recovered slope barely moves off 1.2 at all. Same eleven-plus-one dataset, two different notions of "how much should one bad point matter," two very different outcomes.
r = y_i - \hat y_i is a single residual and δ > 0 is a threshold you choose. Near zero the Huber loss is exactly \tfrac12 r^2 — quadratic, just like OLS's squared loss. Past δ it switches to a straight line in |r| — the same growth rate as the more outlier-tolerant absolute-error loss. The two pieces are constructed to meet with matching value and matching slope exactly at r = ±δ, so the overall loss is smooth (continuously differentiable) everywhere, with no kink a gradient- based solver would trip over.
Part 1 — the Huber loss's derivative is an "influence function." Differentiate each piece of L_δ(r) with respect to r. For |r| ≤ δ:
and for |r| > δ, differentiating δ|r| (the constant -\tfrac12\delta^2 term drops out) and using d|r|/dr = \operatorname{sign}(r):
Call this derivative ψ(r), the loss's influence function — it measures how much a one-unit change in a point's residual pushes the total loss (and hence the fit) around. For small residuals, ψ(r) = r, growing without limit, exactly like squared loss. Past δ, ψ(r) = δ·sign(r) — a constant magnitude, no matter how much further the residual grows. That capped influence, not the loss value itself, is the entire mechanism behind Huber regression's robustness.
Part 2 — that derivative is exactly a reweighted squared loss. Define a per-point weight:
Now differentiate the weighted squared form \tfrac12\,w(r)\,r^2 treating w(r) as a fixed constant (frozen at its current value, not re-differentiated) — exactly the fiction IRLS relies on at every step:
Plugging in each branch of w(r) reproduces Part 1 exactly:
(the last step uses r/|r| = sign(r)). So minimizing the Huber loss and minimizing a squared loss reweighted by w(r) have identical gradients, point for point — provided the weights are computed from the residuals of a fit you already have, then held fixed while you resolve. That's the entire idea behind Iteratively Reweighted Least Squares (IRLS):
— the exact weighted normal equations, the same machinery from 2.3.1's OLS derivation, now weighted. The loop: (i) compute residuals from the current fit, (ii) turn them into Huber weights, (iii) re-solve the weighted normal equations for a new fit, (iv) repeat until the weights (and thus β̂) stop changing. Points with small residuals keep weight 1 and are fit like ordinary OLS; points whose residual has drifted past δ get progressively down-weighted, shrinking their pull on the next iteration's fit.
Where this is used: notice this is the identical reweight-and-resolve pattern used to fit GLMs via IRLS (2.3.9) — freeze weights, solve a weighted least-squares problem in closed form, recompute weights from the new fit, repeat. Only the reason for the weights differs: GLM-IRLS weights come from the variance function of an assumed exponential-family noise distribution (modeling heteroscedastic noise), while Huber-IRLS weights come from how far each point's residual has drifted (down-weighting points that don't look trustworthy). Same outer loop, different job. This M-estimator machinery underlies robust curve-fitting throughout scientific instrumentation and any regression pipeline where "most points are trustworthy but a few might not be" is a reasonable assumption.
Part 3 — RANSAC: a completely different strategy. RANSAC does not touch every point's residual at all. Its loop is: (i) draw a random minimal sample of n points — the fewest needed to determine the model uniquely, 2 for a line, 3 for a plane; (ii) fit the model to just that sample; (iii) count how many of the remaining points fall within some tolerance of that fit ("inliers"); (iv) after many trials, keep whichever fit had the largest inlier count. A fit built from two clean points is essentially exact, and every other clean point will agree with it — a fit contaminated by even one bad point in its minimal sample is generically wrong and attracts almost no agreement. All RANSAC needs is to get lucky at least once across its trials.
How many trials, k, guarantee that with probability at least p? Let w be the assumed fraction of the data that's genuinely inlying. The chance a single random sample of n points is entirely clean is w^n (independent draws, each needing to land in the inlier fraction); the chance a single trial fails to draw an all-clean sample is 1 - w^n. Running k independent trials, the chance every one of them fails is (1-w^n)^k. Requiring this total-failure probability to be at most 1-p (equivalently: at least one clean sample by the end with probability p) and solving for k:
Every symbol has a direct reading: p is how confident you want to be (e.g. 0.99), w is your assumed inlier fraction, n is the minimal sample size the model needs, and k is how many random trials that confidence costs. Both log(1-p) and log(1-w^n) are negative (since p and w^n are both probabilities less than 1), so their ratio comes out positive, as a trial count must. Notice how sharply k depends on n through the exponent w^n — fitting a homography from 4 point correspondences at 50% inliers needs vastly more trials than fitting a line from 2 points at the same inlier rate, because 0.5^4 = 0.0625 is a much rarer event than 0.5^2 = 0.25.
Where this is used: RANSAC's defining strength is that its logic never assumes the outliers are "close" to correct — a completely wild, arbitrary point that happens to occasionally land near a correct-looking fit purely by chance is just as harmless as one that's nowhere near it, because it either lands in enough random minimal samples to matter (rare, and diluted across many trials) or it doesn't get consensus. That is precisely why RANSAC is the workhorse of geometric computer vision — fitting a line, plane, fundamental matrix, or homography to noisy point correspondences where a sizable fraction can be outright false matches, not just noisy measurements of the right thing.
Ten clean points sit near the dashed grey trend line. On load, the amber point drifts from near the line to far below it while the blue OLS fit and red Huber-IRLS fit are recomputed every frame. The live readout tracks both slopes against the true slope of 1.20. Drag the amber point yourself afterward.
Left: squared loss (OLS), absolute loss, and Huber loss as functions of a single residual r — squared loss keeps growing as r², Huber matches it near zero then switches to a straight line past the dashed δ threshold. Right: each loss's derivative, the influence function ψ(r) derived above — squared loss's influence grows without bound, but Huber's flattens to a constant ±δ, which is the entire mechanism behind the OLS-vs-Huber drift shown in the diagram above. Drag δ to move the threshold.
The from-scratch tabs implement the exact IRLS loop derived above — weighted normal equations solved directly (via numpy.linalg.solve in Python, explicit 2×2 algebra in C++), no robust-regression library called anywhere. Both report Huber's recovered slope landing far closer to the true 1.2 than OLS's. The library tab reproduces the same comparison with sklearn's ready-made HuberRegressor and RANSACRegressor, and also prints RANSAC's inlier mask explicitly flagging the planted outlier as not-an-inlier — the direct, visible output of the consensus-counting step derived above.
#include <cmath>
#include <cstdio>
#include <vector>
#include <random>
// Solve a 2x2 weighted normal-equations system explicitly -- no linear-algebra
// library, since the model here is just an intercept and one slope.
struct Fit { double intercept; double slope; };
Fit weightedFit(const std::vector<double>& x, const std::vector<double>& y,
const std::vector<double>& w) {
double sw = 0, swx = 0, swy = 0, swxx = 0, swxy = 0;
for (size_t i = 0; i < x.size(); ++i) {
sw += w[i];
swx += w[i] * x[i];
swy += w[i] * y[i];
swxx += w[i] * x[i] * x[i];
swxy += w[i] * x[i] * y[i];
}
double denom = sw * swxx - swx * swx;
double slope = (sw * swxy - swx * swy) / denom;
double intercept = (swy - slope * swx) / sw;
return {intercept, slope};
}
Fit huberIrls(const std::vector<double>& x, const std::vector<double>& y,
double delta, int iterations = 25) {
std::vector<double> w(x.size(), 1.0);
Fit fit = weightedFit(x, y, w);
for (int iter = 0; iter < iterations; ++iter) {
for (size_t i = 0; i < x.size(); ++i) {
double r = y[i] - (fit.intercept + fit.slope * x[i]);
double absR = std::fabs(r);
w[i] = (absR <= delta) ? 1.0 : delta / std::max(absR, 1e-9);
}
fit = weightedFit(x, y, w);
}
return fit;
}
int main() {
std::mt19937 rng(0);
std::normal_distribution<double> noise(0.0, 0.3);
const double trueIntercept = 1.0, trueSlope = 1.2;
const int nClean = 40;
std::vector<double> xs, ys;
for (int i = 0; i < nClean; ++i) {
double xi = 10.0 * i / (nClean - 1);
xs.push_back(xi);
ys.push_back(trueIntercept + trueSlope * xi + noise(rng));
}
// One deliberate, wildly wrong point.
xs.push_back(8.5);
ys.push_back(trueIntercept + trueSlope * 8.5 - 16.0);
std::vector<double> onesWeight(xs.size(), 1.0);
Fit ols = weightedFit(xs, ys, onesWeight);
Fit huber = huberIrls(xs, ys, 1.0);
double olsErr = std::fabs(ols.slope - trueSlope);
double huberErr = std::fabs(huber.slope - trueSlope);
std::printf("True slope: %.3f\n", trueSlope);
std::printf("OLS recovered slope: %.3f (error %.3f)\n", ols.slope, olsErr);
std::printf("Huber recovered slope: %.3f (error %.3f)\n", huber.slope, huberErr);
std::printf("Huber's slope is %.1fx closer to the true slope than OLS's.\n",
olsErr / std::max(huberErr, 1e-9));
return 0;
}- Computer vision geometric fitting — RANSAC's original and still primary domain. Estimating a homography or fundamental matrix between two camera views starts from automatically-matched feature points, and a sizable fraction of those matches are simply wrong (two unrelated patches that happened to look similar). RANSAC's minimal- sample-and-count-consensus strategy is built exactly for this: outright false correspondences, not just noisy measurements of correct ones.
- Robotics and LiDAR point-cloud processing. Fitting a ground plane or wall to a 3D point cloud for obstacle detection or SLAM runs into the same problem — some returned points are genuine reflections off the surface being fit, others are noise, multi-path reflections, or hits on unrelated objects. Plane-fitting via RANSAC is a standard preprocessing step for exactly this reason.
- Sensor telemetry with occasional glitches. Industrial and IoT sensors occasionally emit a corrupted reading — a dropped bit during transmission, a momentary power sag, a stuck-at-zero fault — sitting far outside the sensor's true physical range. Fitting a trend line to a stream of such readings with plain least squares lets one glitch swing the whole estimated trend; Huber or RANSAC-based fitting is standard in telemetry pipelines for this reason.
- Financial time-series trend estimation. Daily returns occasionally include rare, extreme events — flash crashes, single-day shocks — that are real data points, not errors, but that a risk or trend model usually shouldn't let dominate an estimate of the "typical" underlying trend. Robust regression is a standard tool for estimating trend or beta coefficients that aren't hostage to a handful of extreme days.
- Survey and manually-logged data. Any dataset built from human data entry — lab notebooks, medical intake forms, manually keyed survey responses — accumulates occasional gross errors: a decimal point typed in the wrong place, a height entered in the wrong units, a stray digit. These errors are rarely small; they're often enormous and easy for a robust fit to shrug off, precisely because their residuals are so large under any reasonable model.
- A/B testing and web analytics contaminated by bot traffic. A small fraction of sessions in an experiment's metrics can come from bots, scrapers, or instrumentation bugs rather than genuine users, and their behavior (session length, click counts) can be wildly unlike a real user's. Robust estimators of the treatment effect are a common defense against a handful of such sessions swinging an otherwise clean experiment's conclusion.
- Picking Huber's
δblindly. Set it too small and genuinely normal points routinely land past the threshold, sliding into the linear-loss region and getting needlessly down-weighted — you lose statistical efficiency on perfectly good data. Set it too large and almost nothing ever crosses the threshold, so the fit behaves like plain OLS and inherits its outlier sensitivity right back. A common default,δ ≈ 1.345σfor an estimated noise scaleσ, is a starting point for tuning, not a substitute for checking it against your own data's residual scale. - Feeding RANSAC a badly wrong assumed inlier fraction
w. The trial-count formula derived above,k = log(1-p)/log(1-w^n), is extremely sensitive towthrough the exponentw^n— overestimate the true inlier fraction and you silently undersample the number of trials needed for your target confidencep, and can walk away with an unlucky, contaminated fit while believing it's reliable at confidencep. - Assuming "robust" means "immune to outliers, full stop, no matter how many there are." It doesn't — every robust method has a breakdown point, the largest fraction of arbitrarily bad data it can tolerate before its output can be made arbitrarily bad too (defined precisely below). Huber-style M-estimators and something like RANSAC or median-based estimators sit at very different points on that scale — treating "robust" as one solved, uniform property rather than a spectrum of trade-offs is the single most common misunderstanding of this whole topic.
Going deeper
The breakdown point of an estimator is defined precisely: it's the smallest fraction of the data that an adversary could corrupt — replacing those points with literally arbitrary values, including values driven to infinity — such that the estimator's output can be forced arbitrarily far from correct. It is a worst-case, adversarial notion, not an average-case one. Concrete reference points make the scale tangible. OLS's breakdown point is exactly 0: move even a single point far enough away and its squared residual can be made to dominate the entire sum without any limit, dragging the fitted line arbitrarily far in the process — there's no floor under how little corruption it takes. Huber-style M-estimators do meaningfully better than 0, since the influence function derived above caps at δ no matter how far a residual drifts — but they still land well short of 50%, because that capped influence only limits how much a large residual can matter; it doesn't protect against a point placed at an extreme position in the input/feature space (a leverage point), which can still swing a fitted slope disproportionately even with a modest residual. Methods built to close that remaining gap — least median of squares, least trimmed squares, and RANSAC's consensus-counting among them — approach or reach the 50% ceiling, at the cost of more computation and, on clean data, worse statistical efficiency than Huber regression.
That 50% figure isn't an engineering shortfall waiting to be improved past — it's a hard mathematical ceiling on any reasonable estimator. Once more than half of a dataset is corrupted, "which half is the real signal and which half is the corruption" stops being a hard estimation problem and becomes genuinely ambiguous — an adversary can construct a second, entirely different "clean majority" out of the corrupted points and there is no data-driven way to tell the two apart. This is the frame worth carrying forward: robust regression is not one solved technique, but a family of methods each trading some statistical efficiency on clean data for tolerance of corruption, positioned at different points along that fundamental 0%-to-50% axis.
A dataset is roughly 35% arbitrarily-corrupted points (say, from a badly malfunctioning sensor mixed in with a working one) and 65% genuinely clean. You fit both a Huber M-estimator and RANSAC. Which is more likely to still recover a good fit, and why?
RANSAC is far more likely to succeed here. Huber-style M-estimators have a breakdown point well short of 50% -- with over a third of the data corrupted, especially if any corrupted points act as leverage points (extreme x-values), the IRLS reweighting derived in this lesson can still be overwhelmed, since its weights only cap the influence of large residuals, not of extreme positions in feature space. RANSAC's strategy is different in kind: it repeatedly fits to tiny random minimal samples and keeps whichever fit gets the most agreement, so as long as its assumed inlier fraction w (~0.65 here) and trial count k (via k = log(1-p)/log(1-w^n)) are set correctly, it only needs to get lucky enough times to draw an all-clean minimal sample -- which, at 65% genuine data and a small n like 2, is still a comfortably common event. RANSAC's practical tolerance for corruption extends much further before it breaks down.
OLS's squared loss punishes large residuals quadratically, so a single bad point can drag the whole fit toward it — its breakdown point is exactly 0. The Huber loss fixes this by behaving quadratically only near zero and linearly beyond a threshold δ; minimizing it reduces to IRLS — the same reweight-and-resolve loop as GLM fitting (2.3.9), here driven by residual size rather than assumed noise variance. RANSAC takes an entirely different strategy — fit to tiny random minimal samples, keep whichever wins by consensus — which is why it tolerates far higher corruption fractions, governed by the trial-count formula k = log(1-p)/log(1-w^n). Every robust method sits somewhere on a single axis running from OLS's 0% breakdown point up toward the hard 50% ceiling beyond which "signal" and "corruption" stop being distinguishable at all — choosing a robust method is choosing a point on that trade-off, not opting into blanket immunity to bad data.