Robust Regression Methods

Ordinary Least Squares (OLS) regression is one of the most elegant and widely-used statistical tools ever developed. Its mathematical properties are beautiful: closed-form solutions, unbiased estimates, and optimal efficiency under Gaussian errors. Yet in the trenches of real-world data analysis, OLS harbors a dangerous vulnerability—a single outlier can completely corrupt your model.

Consider fitting a simple linear regression to predict housing prices. 99 homes in your dataset follow a clean linear relationship, but one data entry error records a $500,000 home as $50,000,000. With OLS, that single point doesn't just influence your model—it dominates it. The resulting regression line will twist dramatically toward the outlier, producing wildly inaccurate predictions for every other observation.

This vulnerability isn't a minor inconvenience. In real-world applications—fraud detection, medical diagnostics, autonomous systems, financial modeling—outliers are not exceptions but certainties. Sensor malfunctions, data entry errors, genuine extreme events, and adversarial inputs all generate observations that violate the assumptions underlying OLS.

Robust regression methods address this fundamental problem. They produce reliable estimates even when a fraction of the data is corrupted or comes from a different distribution entirely. At the heart of modern robust regression lies a deceptively simple idea: replace the squared loss with something less sensitive to large residuals.

The Huber loss function is the canonical solution—a hybrid loss that behaves like squared error for small residuals and like absolute error for large ones. This page provides an exhaustive treatment of Huber loss, from its theoretical foundations to practical implementation considerations.

To understand why Huber loss exists, we must first deeply understand what goes wrong with ordinary squared loss. The OLS objective minimizes the sum of squared residuals:

$$\mathcal{L}{\text{OLS}}(\boldsymbol{\beta}) = \sum{i=1}^{n} (y_i - \mathbf{x}i^\top \boldsymbol{\beta})^2 = \sum{i=1}^{n} r_i^2$$

where $r_i = y_i - \mathbf{x}_i^\top \boldsymbol{\beta}$ is the residual for observation $i$.

The quadratic penalty problem:

When we square the residuals, we impose a penalty that grows quadratically with the size of the error. This seems reasonable at first—larger errors should receive larger penalties. But the quadratic growth rate creates severe problems:

Mathematical analysis of influence:

To understand outlier sensitivity rigorously, we examine the gradient of the squared loss with respect to the prediction:

$$\frac{\partial}{\partial \hat{y}} (y - \hat{y})^2 = -2(y - \hat{y}) = -2r$$

The gradient is linear in the residual. This means:

• Large residuals produce large gradients
• Large gradients drive large parameter updates
• A single outlier can shift the entire model

In optimization terms, outliers create extreme gradients that pull the solution toward them. The model sacrifices fit quality on many observations to reduce error on a few extreme ones.

The statistical perspective:

Squared loss is optimal when errors follow a Gaussian distribution—this is the maximum likelihood estimator under $\epsilon \sim \mathcal{N}(0, \sigma^2)$. But real-world error distributions often have heavy tails—more extreme values than a Gaussian would predict. Common heavy-tailed distributions include:

• Student's t-distribution (financial returns, measurement errors)
• Cauchy distribution (ratio of Gaussians, resonance phenomena)
• Laplace distribution (double exponential, L1 loss MLE)
• Contaminated Gaussians (mixture of normal data and outliers)

Under heavy-tailed errors, OLS is no longer optimal—and can be catastrophically inefficient. The sample mean (equivalent to intercept-only OLS) has infinite variance for Cauchy-distributed data, meaning more data doesn't help improve the estimate.

Before introducing Huber loss, let's examine the natural alternative to squared loss: absolute loss (L1 loss, also called Least Absolute Deviations or LAD).

$$\mathcal{L}{\text{LAD}}(\boldsymbol{\beta}) = \sum{i=1}^{n} |y_i - \mathbf{x}i^\top \boldsymbol{\beta}| = \sum{i=1}^{n} |r_i|$$

The L1 loss penalizes residuals linearly rather than quadratically:

Properties of L1 loss:

The L1 loss has remarkable robustness properties:

1. Bounded influence: The gradient of L1 loss is $\frac{\partial |r|}{\partial \hat{y}} = -\text{sign}(r)$, which has magnitude 1 regardless of residual size. Outliers cannot exert disproportionate influence.

2. Median regression: The L1 estimator finds the conditional median rather than the conditional mean. The median is far more robust to outliers than the mean.

3. Breakdown point: L1 regression can tolerate nearly 50% outliers before breaking down (exact breakdown point depends on design matrix geometry).

But L1 has drawbacks:

1. Non-differentiability: The absolute value function has a kink at zero: $|r|$ is not differentiable at $r=0$. This complicates optimization—standard gradient descent doesn't apply directly.

2. Efficiency loss: When errors truly are Gaussian, L1 is less efficient than L2. You need approximately 50% more data to achieve the same precision with L1 vs L2 under Gaussian errors.

3. Non-unique solutions: L1 solutions may not be unique when the data has certain configurations—a technical nuisance for interpretation.

Peter Huber introduced his eponymous loss function in 1964 as a way to achieve the "best of both worlds"—efficiency of L2 near zero while maintaining robustness of L1 for large residuals.

Definition:

The Huber loss function $\rho_\delta(r)$ is defined as:

$$\rho_\delta(r) = \begin{cases} \frac{1}{2}r^2 & \text{if } |r| \leq \delta \ \delta |r| - \frac{1}{2}\delta^2 & \text{if } |r| > \delta \end{cases}$$

where $\delta > 0$ is the threshold parameter (also called the tuning constant or transition point).

Intuition:

• For small residuals ($|r| \leq \delta$): Use quadratic loss, just like OLS
• For large residuals ($|r| > \delta$): Switch to linear loss, like LAD

The transition point $\delta$ determines what counts as "small" vs "large".

The derivative (influence function):

The derivative of Huber loss determines how much each observation influences the parameter estimates:

$$\psi_\delta(r) = \frac{\partial \rho_\delta(r)}{\partial r} = \begin{cases} r & \text{if } |r| \leq \delta \ \delta \cdot \text{sign}(r) & \text{if } |r| > \delta \end{cases}$$

This is called the Huber psi-function or influence function. Key observations:

1. Continuous at transition: At $r = \delta$: from below, $\psi = \delta$; from above, $\psi = \delta \cdot 1 = \delta$. The function is continuous.

2. Bounded influence: For $|r| > \delta$, the influence is capped at $\pm\delta$. No matter how extreme an outlier, its maximum influence is bounded.

3. Smooth near zero: Unlike L1 loss, Huber loss is differentiable everywhere, including at zero. This enables standard gradient-based optimization.

The second derivative:

$$\psi'\delta(r) = \frac{\partial^2 \rho\delta(r)}{\partial r^2} = \begin{cases} 1 & \text{if } |r| < \delta \ 0 & \text{if } |r| > \delta \end{cases}$$

The second derivative is discontinuous at $\pm\delta$, so Huber loss is $C^1$ but not $C^2$ (once but not twice continuously differentiable). This is sufficient for most optimization algorithms.

The threshold parameter $\delta$ is the critical tuning knob of Huber regression. It controls the tradeoff between efficiency (using more of the smooth L2 region) and robustness (treating more residuals as potential outliers).

Extreme cases:

• $\delta \to \infty$: Huber loss approaches L2 loss (OLS behavior)
• $\delta \to 0$: Huber loss approaches L1 loss (LAD behavior)

Practical guidelines:

The standard choice in robust statistics is $\delta = 1.345\sigma$, where $\sigma$ is the standard deviation of the errors. This value has a special property: 95% asymptotic efficiency under Gaussian errors.

The chicken-and-egg problem:

To set $\delta = 1.345\sigma$, we need to know $\sigma$—but we're trying to estimate the model to get residuals to estimate $\sigma$. This circular dependency is resolved through iterative reweighting:

1. Initial estimate: Run OLS or LAD to get initial residuals
2. Robust scale: Estimate $\sigma$ using a robust scale estimator (e.g., MAD)
3. Update δ: Set $\delta = 1.345 \cdot \hat{\sigma}_{\text{robust}}$
4. Re-estimate: Fit Huber regression with updated $\delta$
5. Iterate: Repeat 2-4 until convergence

Robust scale estimation:

The Median Absolute Deviation (MAD) provides a robust estimate of scale:

$$\text{MAD} = \text{median}(|r_i - \text{median}(r_i)|)$$

Under Gaussian errors: $\sigma \approx 1.4826 \cdot \text{MAD}$

The factor 1.4826 makes MAD a consistent estimator of $\sigma$ for Gaussian data.

Unlike OLS, Huber regression does not have a closed-form solution. The non-differentiable transition at $\pm\delta$ (in the second derivative) prevents simple matrix algebra. However, the objective is convex, ensuring that any local minimum is the global minimum.

Iteratively Reweighted Least Squares (IRLS):

The most elegant approach is IRLS, which solves Huber regression through a sequence of weighted least squares problems.

Key insight: The Huber estimating equation can be written as:

$$\sum_{i=1}^{n} \psi_\delta(r_i) \mathbf{x}_i = \mathbf{0}$$

This is equivalent to a weighted least squares problem where the weights are:

$$w_i = \frac{\psi_\delta(r_i)}{r_i} = \begin{cases} 1 & \text{if } |r_i| \leq \delta \ \frac{\delta}{|r_i|} & \text{if } |r_i| > \delta \end{cases}$$

Observations with small residuals receive full weight (1), while outliers receive down-weighted influence proportional to $1/|r_i|$.

Huber estimators possess well-characterized statistical properties that justify their use in practice. Understanding these properties helps practitioners know when and how to apply Huber regression appropriately.

Consistency:

Under mild regularity conditions, the Huber estimator $\hat{\boldsymbol{\beta}}_H$ is consistent:

$$\hat{\boldsymbol{\beta}}_H \xrightarrow{p} \boldsymbol{\beta}_0 \text{ as } n \to \infty$$

This holds even when the error distribution has heavier tails than Gaussian, provided the tails are symmetric and the contamination fraction is below the breakdown point.

Asymptotic normality:

The Huber estimator is asymptotically normal:

$$\sqrt{n}(\hat{\boldsymbol{\beta}}_H - \boldsymbol{\beta}_0) \xrightarrow{d} \mathcal{N}(\mathbf{0}, V_H)$$

The asymptotic covariance matrix involves the Fisher consistency factor:

$$V_H = \frac{\mathbb{E}[\psi_\delta(\epsilon)^2]}{(\mathbb{E}[\psi'_\delta(\epsilon)])^2} (\mathbf{X}^\top \mathbf{X})^{-1}$$

Efficiency:

The asymptotic relative efficiency (ARE) of Huber vs OLS under Gaussian errors is:

$$\text{ARE}{\text{Huber/OLS}} = \frac{(\mathbb{E}[\psi'\delta(\epsilon)])^2}{\mathbb{E}[\psi_\delta(\epsilon)^2]}$$

For $\delta = 1.345\sigma$, this yields ARE ≈ 0.95 (95% efficiency).

We have established a deep understanding of Huber loss as the foundational loss function for robust regression. Let's consolidate the key insights:

What's next:

Huber loss is just the beginning of robust regression. In the next page, we'll explore M-estimators—the general framework of which Huber is a special case. M-estimators allow us to design custom loss functions for different robustness-efficiency tradeoffs, including completely bounded-influence estimators like Tukey's bisquare.