How to Make Linear Regression Survive Outliers

A simple model with a serious weakness

A straight line can look surprisingly convincing—until a few bad measurements pull it somewhere it should never have gone.

Linear regression is often one of the first predictive models practitioners learn—and one of the first they set aside when more sophisticated machine-learning methods become available. Yet linear models remain valuable when coefficients need a physical interpretation, predictions must run on a resource-constrained device, computational latency matters, or a simple benchmark is needed before introducing a higher-capacity model. They are also useful as local approximations: even a complex nonlinear relationship may behave approximately linearly over a sufficiently small region.

Its simplicity, however, comes with an important weakness:

Ordinary Least Squares treats every observation as trustworthy.

In real data, that assumption is easy to violate. A faulty sensor, communication error, calibration problem, or biased measurement can produce observations far from the relationship we actually want to estimate. Because Ordinary Least Squares (OLS) squares every residual, a few such observations can have a disproportionate effect. Robust estimators try to prevent those observations from dominating the fit.

Figures 1 and 2 show how quickly the picture can change. With clean observations, OLS follows the nominal relationship closely. After 30% of the responses are replaced by outliers, the same estimator is pulled sharply away from it. A representative robust fit, however, stays much closer to the relationship supported by the nominal observations.

Figure 1. With clean observations, the OLS fit remains close to the nominal linear relationship.

Figure 2. After 30% of the responses are replaced by outliers, the OLS fit is pulled away from the nominal relationship, while a representative robust fit remains close to it.

That leads to the practical question I want to explore in this article:

How do different robust estimators behave when we do not know the outlier statistics in advance—and when the contamination becomes progressively harder?

Knowing that a dataset contains outliers is only part of the problem. In practice, we rarely know their percentage, bias, variance, distribution, or structure beforehand. A method that works well for one convenient outlier model may behave very differently under another. The estimators here are therefore tested across several deliberately different forms of contamination.

This article compares OLS as the non-robust baseline (Legendre, 1805) with five robust estimators: Huber regression (Huber, 1964), Random Sample Consensus (RANSAC; Fischler and Bolles, 1981), Graduated Non-Convexity with the Geman–McClure loss (GNC-GM), Graduated Non-Convexity with the Truncated Least-Squares loss (GNC-TLS; Yang et al., 2020), and Adaptive Selective Outlier Rejecting (ASOR; Chughtai et al., 2024).

Huber regression and RANSAC are classical robust-estimation methods, whereas GNC-GM, GNC-TLS, and ASOR represent more recent approaches based on non-convex continuation and adaptive residual weighting. Their central algorithmic steps are implemented directly so that the weighting, rejection, sampling, stopping, and continuation mechanisms remain visible.

The estimators are evaluated using prediction error and runtime to capture both statistical accuracy and computational efficiency.

Disclosure. The author developed ASOR in the original study cited here. To ensure a transparent comparison, all estimators are evaluated on the same Monte Carlo realizations using fixed and documented settings.

Why squared loss can be dominated by a few observations

To see why a few bad measurements can have so much influence, consider the scalar linear model:

yi=β0+β1xi+ϵiy_i = \beta_0 + \beta_1 x_i + \epsilon_i

Here, β0β₀ is the intercept, β1β_1 is the slope, and the nominal measurement error follows a Gaussian distribution:

ϵi∼N(0,σin2)\epsilon_i \sim \mathcal{N}\left(0,\sigma_{\mathrm{in}}^2\right)

For NN independent observations, maximizing the likelihood with respect to ββ is equivalent to minimizing the OLS objective:

β^=arg min⁡β∥y−Xβ∥22\hat{\boldsymbol{\beta}} = \underset{\boldsymbol{\beta}}{\operatorname{arg\,min}} \left\|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\right\|_2^2

where

X=[11⋯1x1x2⋯xN]⊤, y=[y1y2⋯yN]⊤, β=[β0β1]⊤% 1. Transpose of the Design Matrix (X^T) \mathbf{X} = \begin{bmatrix} 1 & 1 & \cdots & 1 \\ x_1 & x_2 & \cdots & x_N \end{bmatrix}^{\top},\ % 2. Transpose of the Target Vector (y^T) \mathbf{y} = \begin{bmatrix} y_1 & y_2 & \cdots & y_N \end{bmatrix}^{\top},\ % 3. Transpose of the Parameter Vector (\beta^T) \boldsymbol{\beta} = \begin{bmatrix} \beta_0 & \beta_1 \end{bmatrix}^{\top}

When X\boldsymbol{X} has full column rank, the familiar closed-form expression is:

β^=(X⊤X)−1X⊤y\hat{\boldsymbol{\beta}} = \left(\mathbf{X}^{\top}\mathbf{X}\right)^{-1} \mathbf{X}^{\top}\mathbf{y}

A numerical least-squares solver is preferable to explicitly forming the inverse:

The same squared-loss objective that makes OLS simple and efficient also creates its main weakness. The contribution of an observation grows quadratically with its residual magnitude:

12=1, 102=100, 1002=100001² = 1,\ 10² = 100,\ 100² = 10000

Thus, a normalized residual of 100100 contributes as much to the OLS objective as 1000010000 normalized residuals of 11. A small number of severe outliers can therefore pull the fitted model away from the relationship supported by most observations, as demonstrated in Figure 2.

Practical note on missing values. Rows containing a missing predictor or response can be excluded before fitting:

This treatment is appropriate when missing values are limited and non-systematic. Outliers, however, are different. Missing observations can usually be identified before fitting, whereas outliers must be inferred from residuals that depend on the unknown regression model. Since the fitted model is itself influenced by the outliers, model estimation and outlier identification must be performed jointly.

Six estimators, one shared idea

Table 1. Robustness mechanisms and principal limitations of the six estimators.

| Method

|

Robustness mechanism

|

Limitations

| | --- | --- | --- | |

OLS

|

Assigns equal weight, wi=1w_i=1, to every observation.

|

Residual influence is unbounded, so a small number of severe outliers can substantially shift the fitted model.

| |

Huber

|

Smoothly reduces the influence of large residuals.

|

Severe outliers retain nonzero influence, and performance depends on the selected thresholdδ\delta.

| |

RANSAC

|

Fits random minimal subsets, selects the largest consensus set, and refits using its observations.

|

The method is randomized, requires an inlier threshold, and becomes more expensive as the inlier fraction decreases.

| |

GNC-GM

|

Uses continuation toward a non-convex soft-weighting loss.

|

Weights remain nonzero, so severe outliers may retain influence. The weighting also depends on the nominal-noise scale.

| |

GNC-TLS

|

Uses continuation toward truncated least squares and eventual hard rejection.

|

The result depends on the inlier threshold, and valid observations with unusually large residuals may receive zero weight.

| |

ASOR

|

Uses adaptive posterior probabilities to assign soft observation weights.

|

Its convergence effort can vary across datasets, and its behavior depends on the assumed or estimated nominal-noise scale.

|

A common way to achieve this joint treatment is to control the influence of each observation through a residual-dependent weight. Most methods in this comparison therefore repeatedly solve a weighted least-squares problem:

β^=arg min⁡β∑i=1Nwi ri2(β)\hat{\boldsymbol{\beta}} = \underset{\boldsymbol{\beta}}{\operatorname{arg\,min}} \sum_{i=1}^{N} w_i\,r_i^2(\boldsymbol{\beta})

The sum runs over all NNobservations, and the normalized residual is:

ri(β)=yi−β0−β1xiσinr_i(\boldsymbol{\beta}) = \frac{y_i - \beta_0-{\beta_1}{x}_i } {\sigma_{\mathrm{in}}}

Here, wiw_i controls the influence of the iith observation. The methods differ mainly in how these weights are determined, or whether weighting is replaced by a selected consensus set. Their robustness mechanisms and main limitations are summarized in Table 1.

A shared weighted least-squares engine

Most estimators in this comparison repeatedly solve the same weighted least-squares problem. To keep their method-specific weighting, sampling, and continuation mechanisms visible, they use the following shared numerical solver:

The complete implementations, reproducible notebook, generated figures, and software requirements are available in the public GitHub repository. The repository contains the complete weight updates, sampling rules, stopping criteria, and continuation schedules, while the focused snippets presented here emphasize the distinguishing operation of each estimator.

OLS: use every observation equally

OLS assigns wi=1w_i=1 to every observation, so it fits one line to the entire dataset without distinguishing between nominal measurements and outliers.

How it works. Construct the design matrix, solve one least-squares problem, and use all observations at full weight. OLS requires no iterative stopping rule.

OLS is fast, interpretable, and statistically efficient when the Gaussian model is appropriate. Its limitation is unbounded residual influence: a small number of severe observations can move the fitted line substantially. For straight-line regression, its cost is approximately O(N)\mathcal{O}(N).

Huber regression: reduce influence smoothly

OLS fails because every residual receives its full quadratic penalty. The simplest response is not necessarily to reject suspicious observations completely, but to reduce how strongly large residuals can influence the fit. Huber regression does exactly that: it is quadratic for small residuals and linear for large ones.

The Huber loss is:

ρδ(r)={12r2,∣r∣≤δ,δ∣r∣−12δ2,∣r∣>δ.\rho_{\delta}(r) = \begin{cases} \dfrac{1}{2}r^2, & |r| \leq \delta, \\[6pt] \delta |r| - \dfrac{1}{2}\delta^2, & |r| > \delta. \end{cases}

Its Iteratively Reweighted Least-Squares update is:

wi={1,∣ri∣≤δ,δ∣ri∣,∣ri∣>δ.w_i = \begin{cases} 1, & |r_i| \leq \delta, \\[6pt] \dfrac{\delta}{|r_i|}, & |r_i| > \delta. \end{cases}

How it works. Starting from OLS, Huber regression computes normalized residuals, assigns unit weight below the threshold, reduces the weights above it, and resolves the weighted least-squares problem. The experiments use δ=1.35δ = 1.35 and stop when the normalized change in the regression coefficients is at most 10−510⁻⁵; no fixed iteration cap is imposed.

Huber is a smooth and comparatively inexpensive improvement over OLS. It never assigns exactly zero weight, so severe or systematically biased outliers can continue to influence the estimate. Its performance also depends on the threshold. If IHI_H iterations are required, the straight-line cost is O(IHN)\mathcal{O}(I_H N).

RANSAC: adaptively search for a consensus

Huber still allows every observation to influence the estimate, even when some receive much smaller weights. RANSAC takes a more aggressive view: instead of softening every large residual, it searches directly for a subset of observations that agrees with one model.

For a line, two observations with distinct predictor values define one model hypothesis. RANSAC repeatedly samples two observations, evaluates all residuals, and retains the model with the largest inlier consensus. Whenever a larger consensus is found, the estimated inlier fraction is updated and the required number of trials is recomputed.

How it works. Randomly select an unseen pair of observations, fit a candidate line, compute the normalized residuals, and form a consensus set using a residual threshold. Whenever a larger consensus is found, update the estimated inlier fraction and recompute the number of trials required to achieve confidence p. This adaptive trial control can terminate the search early when a strong consensus is identified. Finally, refit the model using every observation in the winning consensus set.

The experiments use a normalized residual threshold of 33 and confidence p=0.999p=0.999. The initial trial limit is the number of unique two-point subsets (N2)\binom{N}{2}

The approximate number of required hypotheses is:

K≈log⁡(1−p)log⁡(1−ws)K \approx \frac{\log(1-p)} {\log\left(1-w^s\right)}

Here, pp is the desired confidence, ww is the estimated inlier fraction, and s=2s=2 for straight-line regression. Since ww is initially unknown, it is updated adaptively as w^=Nbest/N\hat{w} = {N_{\mathrm{best}}}/{N} where NbestN_{best} is the size of the largest consensus found so far. The updated trial requirement is therefore:

Kadaptive=⌈log⁡(1−p)log⁡(1−w^ 2)⌉K_{\mathrm{adaptive}} = \left\lceil \frac{\log(1-p)} {\log\left(1-\hat{w}^{\,2}\right)} \right\rceil

The ceil operation rounds upward to the nearest integer.

RANSAC is effective when the nominal observations form a distinct and sufficiently large consensus. Adaptive trial control avoids unnecessary hypotheses when a strong consensus is identified early. However, the method remains randomized and threshold-dependent, and its cost increases as the inlier fraction decreases. A large coherent outlier cluster can also become the winning consensus. For KK evaluated hypotheses, the approximate cost is O(KN)\mathcal{O}(K N).

GNC-GM: introduce non-convexity gradually

RANSAC approaches robustness through random sampling and consensus. Graduated Non-Convexity (GNC) takes a different route: instead of searching over subsets, it gradually transforms an easier optimization problem into a more strongly robust, non-convex one.

GNC avoids optimizing a strongly non-convex robust loss in one step. GNC-GM begins with a smoother surrogate and gradually reduces the continuation parameter μ\mu. Its weights are:

wi=(cˉ 2μri2+cˉ 2μ)2w_i = \left( \frac{\bar{c}^{\,2}\mu} {r_i^2+\bar{c}^{\,2}\mu} \right)^2

Following the interpretation of cˉ\bar{c} as an inlier-error bound, this benchmark chooses cˉ2\bar{c}^2 as the 99th99\text{th}percentile of a chi-squared distribution with one degree of freedom: cˉ2≈6.635\bar{c}^2 \approx 6.635. This corresponds to 99%99\% coverage under the assumed nominal Gaussian noise model. The 0.990.99 coverage level is a benchmark setting rather than a value prescribed by the original GNC formulation.

How it works. Initialize with OLS, select a large μ\mu, update the Geman–McClure weights, and solve weighted least squares. After each update, divide μ\mu by 1.41.4. The continuation rule stops the procedure when μ<1\mu<1, with no separate numerical convergence threshold or fixed iteration cap.

The soft weights make GNC-GM deterministic, comparatively stable, and inexpensive, but severely biased observations may retain enough influence to shift the solution. If IGMI_{GM} iterations are required, the straight-line cost is O(IGMN)\mathcal{O}(I_{GM} N).

GNC-TLS: continue toward hard rejection

GNC-GM reduces the influence of large residuals but keeps their weights nonzero. GNC-TLS pushes the same continuation idea further by gradually moving toward hard rejection through the Truncated Least-Squares objective.

For a given μ\mu, its weights are:

wi={1,ri2≤μμ+1cˉ 2,0,ri2≥μ+1μcˉ 2,cˉ 2μ(μ+1)ri2−μ,otherwise.w_i = \begin{cases} 1, & r_i^2 \leq \dfrac{\mu}{\mu+1}\bar{c}^{\,2}, \\[8pt] 0, & r_i^2 \geq \dfrac{\mu+1}{\mu}\bar{c}^{\,2}, \\[8pt] \sqrt{\dfrac{\bar{c}^{\,2}\mu(\mu+1)}{r_i^2}}-\mu, & \text{otherwise}. \end{cases}

How it works. Initialize with OLS and a small continuation parameter, compute the piecewise weights, resolve weighted least squares, and multiply μ\mu by 1.41.4 after each update. The experiments again choose cˉ2\bar{c}^2as the 99th99\text{th} percentile of a chi-squared distribution with one degree of freedom. They stop when the normalized change in the weighted objective, ∑iwiri2\sum_{i} w_{i}r_{i}^{2} , is at most 10−510⁻⁵; no fixed iteration cap is imposed.

This aggressive rejection is useful under strong biased contamination, because sufficiently large residuals receive zero weight. It can be computationally expensive and depends on a suitable inlier threshold. A coherent false structure can still attract the estimate. If ITLSI_{TLS} iterations are required, the straight-line cost is O(ITLSN)\mathcal{O}(I_{TLS} N).

ASOR: update probabilistic soft weights

GNC-GM and GNC-TLS obtain robustness through continuation and residual-dependent weights. ASOR approaches the same problem probabilistically. Rather than immediately deciding whether an observation is nominal or corrupted, it estimates how strongly each explanation is supported by the data and uses that evidence to determine the observation’s influence on the regression model.

In scalar regression, the weight update is:

wi=ωi+(1−ωi)αβiw_i = \omega_i + (1-\omega_i)\frac{\alpha}{\beta_i}

Here, ωi\omega_i is the posterior nominal-component probability, and:

βi=b+12ri2,α=a0+12. \beta_i = b+\frac{1}{2}r_i^2, \qquad \alpha = a_0+\frac{1}{2}.

The posterior nominal-component probability and its scaling factor are computed as:

ωi=[1+ζba0βi−αexp⁡(ri22)]−1, ζ=(1θ−1)Γ(α)Γ(a0). \omega_i = \left[ 1 + \zeta b^{a_0}\beta_i^{-\alpha} \exp\left(\frac{r_i^2}{2}\right) \right]^{-1}, \ \zeta = \left( \frac{1}{\theta}-1 \right) \frac{\Gamma(\alpha)}{\Gamma(a_0)}.

In these expressions, θ\theta is the prior nominal-component probability and Γ\Gamma denotes the gamma function. The parameter bb controls the outlier-component scale and is updated jointly with the regression coefficients. Larger residuals reduce ωi\omega_i and therefore assign greater probability to the outlier explanation.

How it works. Initialize all weights to one, estimate the regression coefficients, compute the normalized squared residuals, update the posterior nominal-component probabilities, update the outlier-scale parameter bb, and form new probabilistic weights.

The experiments use a0=0.5a_0=0.5, b0=104b_0=10^4, θ=0.5\theta=0.5, prior_a=104\texttt{prior\_a} = 10^4, prior_b=103\texttt{prior\_b} = 10^3. The procedure stops when the normalized change in the weighted objective, ∑iwiri2\sum_{i} w_{i}r_{i}^{2} , is at most 10−510⁻⁵; no fixed iteration cap is imposed. These settings are held fixed across all experiments.

ASOR adaptively balances nominal and outlier explanations without forcing an immediate hard decision. Its convergence effort can vary across datasets, and its behavior depends on the assumed or estimated nominal-noise scale. A coherent alternative structure can also attract the estimate. If IAI_{\mathrm{A}} iterations are required, the cost is O(IAN)O(I_{\mathrm{A}}N).

How I stress-tested the estimators

A robust estimator can look impressive under one convenient outlier model and fail badly under another. Rather than relying on a single contaminated dataset, I deliberately vary the amount, distribution, bias, and structure of the corruption.

The nominal relationship throughout the experiments is:

yi=20+3xi+ϵin,i, ϵin,i∼N(0,σin2).\begin{aligned} y_i &= 20+3x_i+\epsilon_{\mathrm{in},i},\ \epsilon_{\mathrm{in},i} \sim \mathcal{N}\left(0,\sigma_{\mathrm{in}}^2\right). \end{aligned}

The predictor values xix_i are evenly spaced over [−100,100][−100, 100]. Two broad contamination families are used.

Independent replacement outliers

For selected observations, the nominal error is replaced by an outlier error. Gaussian outliers follow:

ϵout,i∼N(μout,σout2)\epsilon_{\mathrm{out},i} \sim \mathcal{N}\left( \mu_{\mathrm{out}}, \sigma_{\mathrm{out}}^2 \right)

Uniform outliers follow:

ϵout,i∼U(aout,bout)\epsilon_{\mathrm{out},i} \sim \textit{U}\left( a_{\mathrm{out}}, b_{\mathrm{out}} \right)

The bounds are selected to match the desired mean and variance:

μout=aout+bout2, σout2=(bout−aout)212. \mu_{\mathrm{out}} = \frac{a_{\mathrm{out}}+b_{\mathrm{out}}}{2}, \ \sigma_{\mathrm{out}}^2 = \frac{\left(b_{\mathrm{out}}-a_{\mathrm{out}}\right)^2}{12}.

The outlier statistics are expressed relative to the nominal noise:

μout=κμσin, σout2=κσσin2.\begin{aligned} \mu_{\mathrm{out}} &= \kappa_{\mu}\sigma_{\mathrm{in}}, \ \sigma_{\mathrm{out}}^2 = \kappa_{\sigma}\sigma_{\mathrm{in}}^2. \end{aligned}

Here, κμ \kappa_{\mu} is the mean-shift multiplier and κσ \kappa_{\sigma} is the variance multiplier. A compact version of the generator is:

A coherent competing line

Randomly scattered outliers are only one kind of failure. A more difficult case appears when the corrupted observations agree with one another and form a plausible alternative relationship. To test that situation, corrupted observations also follow:

yiout=20+15xi+ϵout,i, ϵout,i∼N(0,σin2).y_i^{\mathrm{out}} = 20+15x_i+\epsilon_{\mathrm{out},i},\ \epsilon_{\mathrm{out},i} \sim \mathcal{N}\left(0,\sigma_{\mathrm{in}}^2\right).

The nominal and competing relationships have the same intercept and noise variance, but the competing slope is five times larger. This case tests whether an estimator can recover the nominal relationship in the presence of a coherent alternative structure.

Benchmark design

The experiment families and their main configurations are summarized in Table 2.

Table 2. Experiment families and configurations used in the benchmark.

| Experiment

|

Configuration

| | --- | --- | |

Shared setup

|

Nominal model yi=20+3xi+ϵin,i\begin{aligned} y_i &= 20+3x_i+\epsilon_{\mathrm{in},i} \end{aligned} with  ϵin,i∼N(0,σin2)\ \epsilon_{\mathrm{in},i} \sim \mathcal{N}\left(0,\sigma_{\mathrm{in}}^2\right) and xix_i evenly spaced over [−100,100][−100, 100]. Each condition uses 3030 Monte Carlo realizations.

| |

Robustness sweep

|

N=100N=100, σin2=1000\sigma_{in}^2=1000, and outlier percentages of 0%0\%, 10%10\%, 30%30\%, 50%50\%, 70%70\%, and 90%90\%.

| |

Zero-mean Gaussian

|

Gaussian outliers with variance multiplier κσ=10κ_σ = 10 and mean-shift multiplier κμ=0κ_\mu = 0 .

| |

Biased Gaussian

|

Gaussian outliers with κσ=10κ_σ = 10 and κμ=3κ_\mu = 3 .

| |

Biased uniform

|

Uniform outliers matched to the mean and variance of the biased Gaussian case, with κσ=10κ_σ = 10 and κμ=3κ_\mu = 3 .

| |

Competing line

|

Outliers follow yiout=20+15xi+ϵout,iy_i^{\mathrm{out}} = 20+15x_i+\epsilon_{\mathrm{out},i} with ϵout,i∼N(0,σin2) \epsilon_{\mathrm{out},i} \sim \mathcal{N}\left(0,\sigma_{\mathrm{in}}^2\right) .

| |

Sample-size scaling

|

50%50\% biased Gaussian contamination, σin2=1000\sigma_{\mathrm{in}}^2=1000, κσ=10κ_σ = 10, κμ=3κ_\mu = 3 and N=50N = 50, 100100 , 200200, 500500, 10001000, 20002000, 50005000.

| |

Noise scaling

|

N=500N = 500, 50%50\%biased Gaussian contamination, κσ=10κ_σ = 10, κμ=3κ_\mu = 3 and σin2=10,100,1000,10000\sigma_{\mathrm{in}}^2=10, 100, 1000, 10000.

| | | |

All methods receive the same dataset within each Monte Carlo realization to ensure a fair comparison. All random experiments use deterministic seeds derived from base seed 50005000. Their central algorithmic steps are implemented directly and evaluated using fixed and documented settings.

The goal is not to declare a universal winner from one dataset. I want to see which conclusions survive when the contamination mechanism, sample size, and nominal-noise scale change.

Evaluation metrics

I care about two things: does the method recover the right relationship, and how much computation does that robustness cost? Prediction error measures the first, while runtime captures the computational overhead introduced by the robustness mechanism.

This overhead matters when large volumes of high-rate data must be processed or the model is updated repeatedly on a resource-constrained device. For example, short-term SNR prediction for GPS jamming detection may require timely processing of continuously arriving measurements under limited latency, computing, and energy budgets, where excessive computation can reduce battery life. The relative computational costs observed in this simple regression problem therefore provide an early indication of how well each estimator may scale in practical online applications.

Each fitted model is evaluated against the true noiseless relationship on a grid of J=500J=500 equally spaced points over[−100,100][−100, 100]:

RMSE⁡pred=1J∑j=1J[f^(xj)−f(xj)]2\operatorname{RMSE}_{\mathrm{pred}} = \sqrt{ \frac{1}{J} \sum_{j=1}^{J} \left[ \hat{f}(x_j)-f(x_j) \right]^2 }

The sum runs over the J evaluation points. To compare experiments across different nominal-noise scales, the prediction error is normalized as:

NRMSE⁡=RMSE⁡predσin\operatorname{NRMSE} = \frac{\operatorname{RMSE}_{\mathrm{pred}}} {\sigma_{\mathrm{in}}}

Runtime is recorded from one execution of each method in every Monte Carlo realization and summarized using the median of 3030 measurements in each exact condition, because occasional slow executions can distort the mean.

The experiments were performed on an HP ProBook 455 G10 equipped with an AMD Ryzen 7 7730U processor—8 cores, 16 logical processors, 2.0 GHz—and 32 GB of RAM, running Microsoft Windows 11 Pro, Build 22631.

Absolute execution times depend on processor utilization, power mode, operating-system scheduling, Python and library versions, and the BLAS implementation. The relative ordering and scaling trends are therefore more transferable than the exact millisecond values.

Controlled assumption. The benchmark supplies the true nominal noise scale σᵢₙ to residual-normalized methods. In practice, this scale may be calibrated offline from representative clean measurements. Errors in estimating the nominal noise scale alter the normalized residuals and can therefore change threshold-based decisions and probabilistic weights. The reported comparison assumes a known nominal scale to isolate the behavior of the estimators from errors in noise-scale estimation.

What the experiments reveal

Once the contamination mechanism changes, the relative behavior of the estimators changes with it. Five patterns stand out across the experiments.

Finding 1: Biased outliers cause greater systematic distortion than zero-mean outliers

Figure 3 shows zero-mean Gaussian replacement outliers. Positive and negative errors partly cancel, so the estimated relationship is not consistently pushed in one direction. The robust estimators remain closely matched through moderate contamination, and OLS remains more accurate than it does under biased contamination.

Figure 3. Normalized prediction RMSE versus outlier percentage for zero-mean Gaussian replacement outliers. Each box contains 30 Monte Carlo realizations, and the vertical axis is logarithmic.

Figure 4. Normalized prediction RMSE versus outlier percentage for Gaussian replacement outliers with mean shift 3σᵢₙ. Each box contains 30 Monte Carlo realizations, and the vertical axis is logarithmic.

The pattern changes in Figure 4, where the outlier mean is shifted by 3σin3\sigma_{in}. ASOR and GNC-TLS exhibit the lowest errors. At 10% contamination, ASOR and GNC-TLS are nearly tied, with median normalized RMSE values of 0.108 and 0.109. At 30%, 50%, and 70%, GNC-TLS gives median errors of 0.215, 0.327, and 0.751, while the corresponding ASOR values are 0.244, 0.438, and 0.974.

The important distinction is directional consistency. Zero-mean outliers inflate variability, but biased outliers repeatedly pull the fit in the same direction. OLS and Huber are affected most strongly; GNC-TLS and ASOR remain more accurate through moderate contamination. At 90%, the corrupted observations dominate, and the median errors of all six estimators become similar.

Figure 5. Normalized prediction RMSE versus outlier percentage for biased uniform replacement outliers with variance ratio 10 and mean shift 3σᵢₙ. Each box contains 30 Monte Carlo realizations, and the vertical axis is logarithmic.

But is that pattern specific to Gaussian outliers? To check, Figure 5 repeats the experiment with a biased uniform distribution whose mean and variance are matched to the biased Gaussian case.

The broad ranking remains similar to the biased Gaussian case. At low outlier percentages, all robust estimators remain relatively accurate because the nominal observations still dominate the fit. As contamination increases, the differences become more pronounced: GNC-TLS and ASOR maintain the lowest prediction errors through moderate and high outlier levels, while OLS and Huber deteriorate more rapidly.

Finding 2: A coherent competing line is a model-identification problem

So far, the corrupted observations have been independent. The next experiment is harder: what happens when the outliers themselves form a coherent alternative model? Figure 6 evaluates outliers that follow a second line with five times the nominal slope.

Figure 6. Normalized prediction RMSE when corrupted observations follow a competing line with five times the nominal slope. Both relationships have noise variance σᵢₙ² = 1,000, and the vertical axis is logarithmic.

At 10% and 30%, the robust estimators recover the nominal relationship. At 30%, median normalized RMSE values are 0.157 for RANSAC, 0.152 for ASOR, 0.130 for GNC-GM, and 0.154 for GNC-TLS. Huber is less effective, while OLS is pulled strongly toward the competing line.

The problem becomes fundamentally ambiguous at 50% contamination, where the nominal and competing structures contain the same number of observations. The wide boxes indicate that different Monte Carlo realizations may lead robust estimators toward either of the two coherent relationships.

At 70% and 90%, the competing line is the dominant structure. The robust estimators generally select it, producing normalized RMSE values near 22 relative to the nominal model. This is an identifiability limitation, not merely a numerical failure. Without labels, physical constraints, temporal information, or a multi-model formulation, the data alone do not reveal which coherent relationship is the intended one.

Finding 3: Robustness requires additional computation

Accuracy is only one side of the story. Robustness requires additional computation, and the methods pay very different prices for it. Figure 7 shows runtime under biased Gaussian contamination.

Figure 7. Execution time versus outlier percentage under biased Gaussian contamination. The vertical axis is logarithmic.

OLS remains fastest because it requires one least-squares solve. The robust estimators perform repeated weighted solves or evaluate multiple RANSAC hypotheses. GNC-GM generally provides the lowest and most stable iterative cost, while GNC-TLS is generally the most expensive because its continuation schedule requires more updates. ASOR occupies an intermediate accuracy–runtime region.

Absolute sub-millisecond timings depend on the processor, power mode, operating-system scheduling, Python version, and BLAS implementation. The relative ordering and scaling trends are more transferable than the exact millisecond values.

Table 3 provides a representative comparison across all four contamination models at 50% contamination.

Table 3. Median execution time in milliseconds at 50% contamination for N = 100 and 30 Monte Carlo realizations.

| Method

|

Zero-mean Gaussian

|

Biased Gaussian

|

Biased uniform

|

Competing line

| | --- | --- | --- | --- | --- | |

OLS

|

0.085

|

0.094

|

0.084

|

0.086

| |

Huber

|

0.877

|

1.272

|

1.183

|

9.593

| |

RANSAC

|

1.054

|

1.347

|

1.302

|

2.534

| |

ASOR

|

1.229

|

2.018

|

1.489

|

1.773

| |

GNC-GM

|

0.939

|

1.083

|

0.876

|

1.477

| |

GNC-TLS

|

2.451

|

2.732

|

2.251

|

2.891

|

OLS remains substantially faster than the robust estimators. Among the iterative methods, GNC-GM generally has the lowest and most stable runtime. Huber, RANSAC, and ASOR have intermediate computational costs, whereas GNC-TLS is typically the most expensive because its continuation schedule requires multiple weighted least-squares solves.

The competing-line case produces an unusually large Huber runtime at 50% contamination, indicating slower convergence for this particular configuration. This isolated result should not be interpreted as a general runtime property of Huber regression.

Finding 4: More observations reduce variability, not systematic bias

A natural question is whether simply collecting more data makes the contamination problem disappear. Figures 8 and 9 fix biased Gaussian contamination at 50% and vary the sample size from 50 to 5000.

Figure 8. Normalized prediction RMSE versus sample size under 50% biased Gaussian contamination. Each box contains 30 Monte Carlo realizations, and the vertical axis is logarithmic.

Figure 9. Execution time versus sample size under 50% biased Gaussian contamination. Each box contains 30 timing measurements, and the vertical axis is logarithmic.

OLS remains biased as N increases: its median normalized RMSE is 1.533 at N=50N=50 and 1.500 at N=5000N=5000. More observations do not remove bias when the same contamination mechanism persists in the larger dataset.

GNC-TLS has the lowest median error at every tested sample size. Its median decreases from 0.411 at N=50N=50 to 0.254 at N=200N=200 and then stabilizes near 0.24–0.27. ASOR decreases from 0.582 at N=50N=50 to 0.403 at N=5000N=5000, while RANSAC decreases from 0.753 to 0.485. The most visible benefit of increasing N is narrower variability among the robust estimates.

Runtime increases with sample size. From N=50N=50 to N=5000N=5000, the median execution time increases from 0.084 to 0.182 ms for OLS, 1.220 to 2.936 ms for Huber, 1.228 to 1.652 ms for RANSAC, 1.745 to 4.081 ms for ASOR, 0.946 to 2.733 ms for GNC-GM, and 2.136 to 9.603 ms for GNC-TLS.

RANSAC grows more slowly because its adaptive stopping rule evaluates a similar number of hypotheses across the tested sample sizes, although each hypothesis becomes more expensive as N increases.

Finding 5: Normalization preserves the ranking across noise scales

Finally, I change the absolute noise scale while preserving the relative contamination strength. This checks whether the observed ranking is tied to one particular measurement scale. The final experiment varies:σin2∈{10,100,1000,10000}\sigma_{\mathrm{in}}^2 \in \{10, 100, 1000, 10000\}. The outlier mean and variance are scaled relative to the nominal noise. Table 4 shows the median normalized prediction error.

Table 4. Median normalized prediction RMSE as the nominal-noise variance changes. Boldface marks the two lowest values in each column.

| Method

|

σin2=10\sigma_{\mathrm{in}}^2=10

|

σin2=100\sigma_{\mathrm{in}}^2=100

|

σin2=1000\sigma_{\mathrm{in}}^2=1000

|

σin2=10000\sigma_{\mathrm{in}}^2=10000

| | --- | --- | --- | --- | --- | |

OLS

|

1.523

|

1.499

|

1.490

|

1.521

| |

Huber

|

0.891

|

0.840

|

0.836

|

0.877

| |

RANSAC

|

0.529

|

0.596

|

0.594

|

0.490

| |

ASOR

|

0.410

|

0.409

|

0.380

|

0.427

| |

GNC-GM

|

0.571

|

0.540

|

0.512

|

0.560

| |

GNC-TLS

|

0.291

|

0.287

|

0.265

|

0.300

|

The normalized results remain broadly stable as σin2\sigma_{\mathrm{in}}^2changes from 10 to 10000. Across the four variance levels, the median normalized RMSE lies between 1.490 and 1.523 for OLS, 0.836 and 0.891 for Huber, 0.490 and 0.596 for RANSAC, 0.380 and 0.427 for ASOR, 0.512 and 0.571 for GNC-GM, and 0.265 and 0.300 for GNC-TLS.

The leading methods remain consistent across noise scales: GNC-TLS has the lowest median normalized error, followed by ASOR, whereas OLS has the largest error. The small variations between columns are consistent with finite Monte Carlo variation and show no systematic dependence on the absolute measurement scale.

Which estimator should you start with?

There is no universal winner. The right starting point depends on what you know about the data, how aggressively you are willing to reject observations, and how much computation you can afford. Domain knowledge, residual diagnostics, computational constraints, and sensitivity analysis should all inform the choice. Practical starting points are summarized in Table 5.

Table 5. Practical starting points for estimator selection.

| Observed need or evidence

|

Reasonable starting point

| | --- | --- | |

Data appear clean and speed is critical.

|

OLS. Inspect the residuals and influential observations before trusting the result.

| |

Mild contamination is plausible and a smooth fit is preferred.

|

Huber. Compare its result with OLS and at least one stronger robust estimator.

| |

An acceptable residual tolerance can be specified.

|

RANSAC. It is suitable when valid observations are expected to lie within a known tolerance of the underlying relationship.

| |

A nominal-noise scale can be specified and deterministic soft weighting is preferred.

|

GNC-GM. It provides a moderate-cost compromise without hard rejection, but its weighting depends on the specified or estimated nominal-noise scale.

| |

Aggressive rejection is acceptable and thresholds can be validated.

|

GNC-TLS. It achieved the strongest overall accuracy here, but requires greater computation and threshold-sensitivity checks.

| |

Probabilistic adaptive weighting is desired.

|

ASOR. It provides an accuracy–runtime compromise when the nominal-noise scale can be specified or robustly estimated.

| |

Two coherent structures may be present.

|

Use a multi-model approach. Consider mixture regression, multi-model fitting, labels, temporal continuity, or physical constraints rather than relying on a single robust line.

|

The table is a starting point, not a decision rule. In practice, I would fit more than one robust estimator, compare the resulting coefficients and residual patterns, and test sensitivity to reasonable scale and threshold choices.

Agreement across methods increases confidence in the recovered relationship; strong disagreement is itself useful evidence that the data may contain multiple structures or that the assumed nominal-noise scale needs to be reconsidered.

What this benchmark does not establish

There are also clear limits to what these experiments tell us. This is a controlled scalar-regression benchmark, not a universal leaderboard for robust estimation. The benchmark studies a scalar linear model with one-dimensional residuals, known nominal noise, fixed algorithmic settings, and synthetic contamination.

Different conclusions may emerge for high-dimensional regression, leverage points in the predictor space, heteroscedastic noise, nonlinear models, correlated errors, or real datasets with unknown ground truth.

The runtime results are implementation- and machine-dependent. The experiments also do not solve model identification when multiple coherent structures are present. Robust residual weighting can suppress isolated corruption, but it cannot determine the intended model without additional information once an alternative structure dominates.

What should you take away from this?

The most useful lesson from these experiments is broader than the ranking of the six methods: the percentage of outliers alone does not determine how difficult a regression problem is. Zero-mean outliers can partly cancel, while biased observations repeatedly pull the estimate in the same direction. Even Gaussian and uniform outliers with matched first two moments can produce different errors. Distribution, bias, and structure all matter.

More data help with variability, but not necessarily with bias. When the same contamination mechanism persists, increasing the sample size does not make OLS converge back to the nominal relationship.

In these experiments, GNC-TLS achieved the strongest overall accuracy when aggressive rejection was beneficial, ASOR provided a favorable accuracy–runtime compromise, and GNC-GM offered relatively stable computational behavior. These rankings are useful, but they depend on the assumptions and contamination models used here; they should not be treated as a universal leaderboard.

The competing-line experiment gives the strongest warning. Once corrupted observations form a coherent alternative relationship, robust regression is no longer simply an outlier-rejection problem. It becomes a model-identification problem, and residual weighting alone cannot tell us which coherent structure is the one we intended to recover.

Since we rarely know exactly how outliers will appear in practice, robust estimators are best judged across several plausible contamination patterns, with both accuracy and computational cost in view.

Next in the series. The robust regression series continues with polynomial regression, extending the comparison to nonlinear relationships and higher model complexity.

References

1. Legendre, A. M. (1805). Nouvelles méthodes pour la détermination des orbites des comètes. F. Didot.

2. Huber, P. J. (1964). Robust Estimation of a Location Parameter. The Annals of Mathematical Statistics.

3. Fischler, M. A., and Bolles, R. C. (1981). Random Sample Consensus: A Paradigm for Model Fitting with Applications to Image Analysis and Automated Cartography. Communications of the ACM.

4. Yang, H., Antonante, P., Tzoumas, V., and Carlone, L. (2020). Graduated Non-Convexity for Robust Spatial Perception: From Non-Minimal Solvers to Global Outlier Rejection. IEEE Robotics and Automation Letters.

5. Chughtai, A. H., Tahir, M., and Uppal, M. (2024). Bayesian Heuristics for Robust Spatial Perception. IEEE Transactions on Instrumentation and Measurement.

Get in touch 👋

For more of my work, explore my GitHub, or connect with me on LinkedIn. I welcome questions, new ideas, and opportunities to collaborate in Data science, AI and Statistical Signal Processing. If you enjoy my articles, sharing them with others helps these conversations reach more people.

···

All data-driven figures were generated via simulations by the author. The featured image was created with the assistance of AI.