Overview of M-Estimation

Robust statistics provides tools for dealing with data that contains outliers or deviations from standard model assumptions. Among the various approaches developed—such as L-estimators, R-estimators, and M-estimators—M-estimators (Maximum Likelihood-type estimators) are among the most widely utilized in regression analysis. Introduced by Huber, they generalize the maximum likelihood estimation principle to reduce the influence of outliers.

The core objective of an M-estimator is to minimize a cost function defined as:

S(β) = min ∑ i=1n ρ(ri)

Where β represents the parameters to be estimated, ri is the residual (error) for the i-th observation, and ρ is a robust loss function. To ensure robustness, ρ must satisfy specific properties: it must be continuous, non-negative, symmetric (even function), pass through the origin, and be monotonically increasing for positive values. Common choices for ρ include Huber, Cauchy, Tukey, and Welsch functions.

Solving M-Estimators via IRLS

Directly minimizing the sum of the robust loss function is often computationally difficult. Instead, the standard approach is the Iteratively Reweighted Least Squares (IRLS) algorithm. This method transforms the robust optimization problem into a sequence of weighted least squares problems.

To derive the IRLS formulation, we take the partial derivative of the objective function with respect to the parameters β and set it to zero:

∂S / ∂β = ∑ i=1n ψ(ri) · ∂ri / ∂β = 0

Here, ψ(r) is the derivative of ρ(r), often referred to as the influence function. We define a weight function w(r) as:

w(r) = ψ(r) / r

Substituting this into the derivative equation yields:

∑ i=1n w(ri) ri ∂ri / ∂β = 0

This equation is mathematically equivalent to the normal equations of a weighted least squares problem where the weights are determined by the residuals from the previous iteration. Therefore, the solution can be found iteratively:

β(k+1) = (XT W(k) X)-1 XT W(k) y

The algorithm proceeds as follows:

1. Initialization: Compute an initial parameter estimate β(0), typically using standard Ordinary Least Squares (OLS).
2. Calculate Weights: Using the current parameters, compute residuals ri and determine the corresponding weights w(ri) using the chosen robust function.
3. Update Parameters: Solve the weighted least squares problem to obtain new parameters β(k+1).
4. Convergence Check: If the change in parameters is below a threshold ε or the maximum iterations are reached, stop. Otherwise, return to step 2.

The Huber Loss Function

The Huber loss function is a popular choice because it behaves quadratically for small residuals (like OLS) and linear for large residuals, reducing the impact of outliers.

Given a threshold parameter k:

• ρ(r) (Loss):
  If |r| ≤ k: r2 / 2
  If |r| > k: k|r| - k2 / 2
• ψ(r) (Influence):
  If r < -k: -k
  If |r| ≤ k: r
  If r > k: k
• w(r) (Weight):
  If r < -k: -k / r
  If |r| ≤ k: 1
  If r > k: k / r

C++ Implementation with Eigen

The folllowing C++ example demonstrates IRLS regression using the Huber loss function. We use the Eigen library for matrix operations. The data represents a linear relationship y = 2x + 1 with significant outliers added to test robustness.

#include <iostream>
#include <vector>
#include <cmath>
#include <Eigen/Dense>

int main() {
    // Configuration: Huber threshold and max iterations
    const double huber_threshold = 1.5;
    const int max_iterations = 20;
    const double tolerance = 1e-3;
    const int sample_count = 7;

    // Data Points: y approx 2x + 1, with outliers
    // Indices 5 and 6 (index 0-based) are outliers relative to the ground truth
    std::vector<double> x_vals = {1.0, 2.1, 2.9, 5.01, 8.093, 6.0, 3.0};
    std::vector<double> y_vals = {3.02, 4.97, 7.1, 10.88, 17.06, 2.0, 17.6};

    // Construct the Design Matrix X (n x 2) and Observation Vector Y (n x 1)
    // Column 0: x values, Column 1: bias (1.0)
    Eigen::MatrixXd design_matrix(sample_count, 2);
    Eigen::VectorXd observations(sample_count);

    for (int i = 0; i < sample_count; ++i) {
        design_matrix(i, 0) = x_vals[i];
        design_matrix(i, 1) = 1.0;
        observations(i) = y_vals[i];
    }

    // Step 1: Initial guess using Ordinary Least Squares (OLS)
    // params = (X^T * X)^-1 * X^T * y
    Eigen::VectorXd params = (design_matrix.transpose() * design_matrix).inverse() 
                             * design_matrix.transpose() * observations;

    std::cout << "Initial OLS Result:\n" << params << "\n" << std::endl;

    // IRLS Iteration Loop
    int iter = 0;
    bool converged = false;

    while (iter < max_iterations && !converged) {
        // Calculate residuals: r = X * params - y
        Eigen::VectorXd residuals = design_matrix * params - observations;

        // Calculate Weights Matrix W (diagonal matrix)
        Eigen::DiagonalMatrix<double, Eigen::Dynamic> weight_matrix(sample_count);

        for (int i = 0; i < sample_count; ++i) {
            double r = residuals(i);
            
            if (std::abs(r) <= huber_threshold) {
                // Small residuals: weight is 1.0 (quadratic zone)
                weight_matrix.diagonal()(i) = 1.0;
            } else {
                // Large residuals: weight decreases linearly (k/|r|)
                weight_matrix.diagonal()(i) = huber_threshold / std::abs(r);
            }
        }

        // Step 3: Solve Weighted Least Squares
        // params_new = (X^T * W * X)^-1 * X^T * W * y
        Eigen::VectorXd params_new = (design_matrix.transpose() * weight_matrix * design_matrix).inverse()
                                    * design_matrix.transpose() * weight_matrix * observations;

        // Step 4: Check for convergence
        double diff = (params_new - params).norm();
        if (diff < tolerance) {
            converged = true;
        }
        
        params = params_new;
        iter++;
    }

    std::cout << "Final IRLS Result (Robust):\n" << params << std::endl;
    std::cout << "Iterations: " << iter << std::endl;

    return 0;
}

Analysis of Results

Running the algorithm typically yields an initial OLS solution heavily skewed by the outliers (e.g., slope near 1.09, intercept near 4.56). After applying the Huber-based IRLS process, the parameters converge closer to the ground truth (slope near 1.84, intercept near 1.59). This demonstrates the robustness of M-estimation compared to standard least squares, as the iterative weighting mechanism effectively down-weights the influence of the anomalous data points.