Robust and Quantile Regression in C++

July 2026 · 50 min read

Every solver so far in this series, OLS, QR/SVD, gradient descent, and the ridge/lasso/elastic-net family from the last post, minimizes some form of squared error. Penalizing the coefficients changes how much a fixed-shape squared-error bowl is allowed to move; it never changes the shape of the bowl itself. That shape is the problem this post addresses: squared error weighs a residual of $10$ a hundred times more heavily than a residual of $1$, so a handful of extreme observations, a data-entry error, a sensor glitch, a genuinely unusual customer, can dominate a fit that every other point would otherwise agree on. This post replaces the loss itself. Huber regression blends squared and absolute error to stay robust to outliers while remaining smooth. Quantile regression (with the pinball, or check, loss) estimates a conditional quantile instead of the conditional mean, and reduces to LAD (least absolute deviations, median regression) at the midpoint. Both are fit by iteratively reweighted least squares (IRLS), a technique that reuses the exact weighted normal equations a diagonal $\Lambda$ was folded into in the ridge post, except the diagonal is recomputed from the residuals at every iteration instead of fixed in advance.

Why Squared Loss Is Fragile

The training objective behind every earlier post is $J(\beta)=\sum_i \tfrac12 r_i^2$ with $r_i = y_i - \hat y_i$. Its derivative with respect to a prediction is $\partial J/\partial \hat y_i = -r_i$: the influence of a single point on the fit grows linearly and without bound as its residual grows. The derivative of absolute error, $\partial |r_i|/\partial \hat y_i = -\operatorname{sign}(r_i)$, behaves oppositely: it is $\pm 1$ regardless of how large $r_i$ gets. A concrete table makes the contrast impossible to miss:

A point with a residual of $100$ pulls a squared-error fit $100\times$ harder than a point with a residual of $1$; the same point pulls an absolute-error fit exactly as hard as any other point with a nonzero residual, no matter how extreme. Robust regression sits between these two extremes: Huber loss caps the influence at a tunable threshold $\delta$ rather than at $1$, and the pinball loss used for quantile regression is a tilted version of the same capped-influence idea.

Huber Loss and Its Gradient

Huber loss is exactly squared error for small residuals and exactly (a shifted) absolute error for large ones, matched at $|r|=\delta$ so the function and its derivative are both continuous:

The derivative, sometimes called the clipped residual, is exactly $r$ near the origin and is capped at $\pm\delta$ everywhere else, which is precisely the bounded-influence property the table above showed absolute error has, but Huber loss keeps the smooth, fully differentiable squared-error behavior near zero, where being wrong by a small amount should still push the fit to improve at the ordinary quadratic rate.

Minimizing $\sum_i L_\delta(y_i - x_i^\top\beta)$ has no closed form once any residual crosses $\delta$, but the derivative has a useful algebraic property: it can always be written as $\dfrac{dL_\delta}{dr} = w(r)\cdot r$ for a weight $w(r) = \min\!\bigl(1,\ \delta/|r|\bigr)$. Both cases confirm it: for $|r|\le\delta$, $w(r)=1$ and $w(r)\cdot r = r$, matching the top branch exactly. For $|r|>\delta$, $w(r)=\delta/|r|$ and $w(r)\cdot r = \delta\cdot r/|r| = \delta\operatorname{sign}(r)$, matching the bottom branch exactly. Setting the gradient of $\sum_i L_\delta(r_i)$ to zero is therefore algebraically identical to setting the gradient of the weighted squared-error objective $\sum_i w(r_i)\, r_i^2$ to zero, as long as the weights $w(r_i)$ are held fixed while differentiating. That equivalence is the entire algorithm: alternate between computing weights from the current residuals and re-solving a weighted least-squares problem with those weights held fixed, and the fixed point of that alternation is a stationary point of the true Huber objective.

IRLS by Hand: Watching an Outlier Lose Its Vote

Six points make the reweighting mechanical. Let $x=[-2.5,-1.5,-0.5,0.5,1.5,2.5]$ with a clean generating line $y=1+2x=[-4,-2,0,2,4,6]$, except the last point is corrupted to $y_6=16$ (a $+10$ outlier at the highest-leverage end of the range). The unweighted OLS fit on this contaminated data is already dragged: $\beta_{\text{OLS}}=(2.6667,\ 3.4286)$, far from the true $(1,2)$. IRLS with $\delta=1.0$ starts from that OLS fit and alternates weighting and refitting:

1. Iteration 1. Residuals from the OLS start: $r=[1.9048,\,0.4762,\,-0.9524,\,-2.3810,\,-3.8095,\,4.7619]$. Weights $w_i=\min(1,\delta/|r_i|)$: $[0.5250,\,1,\,1,\,0.4200,\,0.2625,\,0.2100]$, so the two points with the largest residuals, the corrupted point 6 ($|r_6|=4.76$) and its spatial neighbor point 5 ($|r_5|=3.81$, "masked" by the outlier dragging the fit toward it), are the two most heavily downweighted; the two points OLS already fit well (residuals $\approx 0$) keep full weight $1$. Re-solving the weighted normal equations gives $\beta=(2.2733,\ 3.0310)$.
2. Iteration 2. New residuals $r=[1.3041,\,0.2732,\,-0.7578,\,-1.7887,\,-2.8197,\,6.1494]$; weights $[0.7668,\,1,\,1,\,0.5591,\,0.3546,\,0.1626]$. The outlier's residual actually grew (the fit is un-dragging from it, so its apparent error increases) while its weight keeps falling. $\beta=(1.9206,\ 2.7084)$.
3. Iteration 3. $r=[0.8503,\,0.1420,\,-0.5664,\,-1.2748,\,-1.9832,\,7.3084]$; weights $[1,\,1,\,1,\,0.7844,\,0.5042,\,0.1368]$. Three points now sit fully inside $\pm\delta$ and carry weight $1$; the outlier's weight has fallen further, to $0.1368$. $\beta=(1.6420,\ 2.4895)$.

The fit keeps moving toward the true line while the outlier's vote keeps shrinking. Continuing this alternation (verified by script, not shown step by step) converges at iteration $15$ to exactly $\beta=(1.35,\ 2.30)$, much closer to the true $(1,2)$ than OLS's $(2.6667,\ 3.4286)$, though not exact: with only $6$ points and one of them corrupted by $10$ units, even a robust estimator cannot fully undo the contamination, a point the complete program's larger, more realistic Scenario A revisits.

HuberRegressor Implementation

Every IRLS iteration is one call to a weighted normal-equation solver, which is the unweighted Part-1 fitOLSDirect generalized to accept a per-row weight vector (setting every weight to $1$ recovers Part 1 exactly, which is how IRLS is initialized below):

struct LinFit { double intercept; vector<double> coef; };

LinFit fitWeightedOLS(const Matrix& X, const vector<double>& y, const vector<double>& w) {
    int n = (int)X.size(), d = (int)X[0].size(), p = d + 1;
    Matrix A(p, vector<double>(p, 0.0));
    vector<double> b(p, 0.0);
    for (int i = 0; i < n; ++i) {
        vector<double> row(p);
        row[0] = 1.0;
        for (int j = 0; j < d; ++j) row[j+1] = X[i][j];
        for (int a = 0; a < p; ++a) {
            b[a] += w[i] * row[a] * y[i];
            for (int c = 0; c < p; ++c) A[a][c] += w[i] * row[a] * row[c];
        }
    }
    vector<double> beta = solveLinearSystem(A, b);   // Part 1's partial-pivoting solver
    LinFit out;
    out.intercept = beta[0];
    out.coef.assign(beta.begin()+1, beta.end());
    return out;
}

LinFit fitOLSDirect(const Matrix& X, const vector<double>& y) {
    validateFitInput(X, y, "fitOLSDirect");
    return fitWeightedOLS(X, y, vector<double>(X.size(), 1.0));   // every weight = 1
}

vector<double> residualsOf(const Matrix& X, const vector<double>& y, const LinFit& f) {
    vector<double> r(X.size());
    for (size_t i = 0; i < X.size(); ++i) {
        double pred = f.intercept;
        for (size_t j = 0; j < X[i].size(); ++j) pred += f.coef[j] * X[i][j];
        r[i] = y[i] - pred;
    }
    return r;
}

Both robust regressors below fit in standardized feature coordinates rather than in raw ones. A short train-only scaler, built only from each fit's own training X, centers and scales every non-constant column to zero mean and unit population variance, drops any constant column so the weighted solve cannot go singular on it (its raw coefficient is then fixed at exactly $0$), and records the means and scales needed to map the fit back. The target y is never rescaled, so residuals and $\delta$ stay in raw $y$ units:

struct Standardizer {
    vector<double> mu, sigma;   // per original feature column
    vector<int> keep;           // original indices of the non-constant columns
    int d = 0;
};

Standardizer standardizeFeatures(const Matrix& X, Matrix& Z) {
    int n = (int)X.size(), d = (int)X[0].size();
    Standardizer s; s.d = d;
    s.mu.assign(d, 0.0);
    s.sigma.assign(d, 1.0);
    for (int i = 0; i < n; ++i) for (int j = 0; j < d; ++j) s.mu[j] += X[i][j];
    for (int j = 0; j < d; ++j) s.mu[j] /= n;
    vector<double> maxDev(d, 0.0);
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < d; ++j) maxDev[j] = max(maxDev[j], fabs(X[i][j] - s.mu[j]));
    for (int j = 0; j < d; ++j) {
        if (maxDev[j] == 0.0) { s.sigma[j] = 1.0; continue; }   // constant column: dropped, raw coef 0
        double acc = 0.0;
        for (int i = 0; i < n; ++i) { double z = (X[i][j] - s.mu[j]) / maxDev[j]; acc += z * z; }
        s.sigma[j] = maxDev[j] * sqrt(acc / n);                 // population standard deviation
        s.keep.push_back(j);
    }
    int k = (int)s.keep.size();
    Z.assign(n, vector<double>(k));
    for (int i = 0; i < n; ++i)
        for (int c = 0; c < k; ++c) { int j = s.keep[c]; Z[i][c] = (X[i][j] - s.mu[j]) / s.sigma[j]; }
    return s;
}

// Standardized-space fit -> raw units: coef_j = betaTilde_c / sigma_j for kept
// columns, exactly 0 for constant columns; the intercept absorbs the centering.
LinFit toRawFit(const LinFit& z, const Standardizer& s) {
    LinFit raw; raw.coef.assign(s.d, 0.0);
    raw.intercept = z.intercept;
    for (size_t c = 0; c < s.keep.size(); ++c) {
        int j = s.keep[c];
        raw.coef[j] = z.coef[c] / s.sigma[j];
        raw.intercept -= raw.coef[j] * s.mu[j];
    }
    return raw;
}

The unpenalized Huber and quantile estimands are invariant to this rescaling: centering and scaling the columns is an invertible linear reparametrization of $\beta$, so the achievable predictions, and therefore the raw coefficients that toRawFit recovers, are exactly what a raw-space fit would produce. What standardization changes is the conditioning of the weighted normal equations $Z^\top W Z$ that IRLS re-solves every iteration. The weights swing across orders of magnitude as the loss reweights rows (Huber pushes outlier weights toward zero; the quantile weight $c/\lvert r\rvert$ becomes very large near zero residual until the implementation's explicit $\varepsilon$ floor caps it), and on an un-centered, differently-scaled design that makes $X^\top W X$ badly conditioned, so each Gaussian-elimination solve sheds precision. Standardized columns share one scale and are decorrelated from the intercept column, holding the condition number down so every re-solve stays accurate. Scaling never moves the optimum; it makes the repeated linear algebra that finds it reliable, which is exactly why the penalized estimators of the previous post standardized too.

struct HuberRegressor {
    bool fitted = false;
    int nFeatures = 0;
    double intercept = 0.0;
    vector<double> coef;
    double delta = 1.0;
    int iterationsRun = 0;
    bool converged = false;
    vector<double> weights;   // final IRLS weights, one per training row

    void fit(const Matrix& X, const vector<double>& y, double deltaIn,
             int maxIter = 200, double tol = 1e-9) {
        fitted = false;
        validateFitInput(X, y, "HuberRegressor");
        if (!isfinite(deltaIn) || deltaIn <= 0.0)
            throw invalid_argument("HuberRegressor::fit: delta must be finite and positive");
        delta = deltaIn;
        int n = (int)X.size(), d = (int)X[0].size();

        Matrix Z;
        Standardizer scaler = standardizeFeatures(X, Z);   // means/scales from THIS training X only
        LinFit cur = fitWeightedOLS(Z, y, vector<double>(n, 1.0));   // unit-weight OLS start
        vector<double> w(n, 1.0);
        int iter = 0;
        bool didConverge = false;
        for (; iter < maxIter; ++iter) {
            vector<double> r = residualsOf(Z, y, cur);   // residuals stay in raw y units
            for (int i = 0; i < n; ++i) {
                double a = fabs(r[i]);
                w[i] = (a <= delta) ? 1.0 : delta / a;
            }
            LinFit next = fitWeightedOLS(Z, y, w);
            double maxDelta = fabs(next.intercept - cur.intercept);
            for (size_t j = 0; j < next.coef.size(); ++j) maxDelta = max(maxDelta, fabs(next.coef[j] - cur.coef[j]));
            cur = next;
            if (maxDelta < tol) { ++iter; didConverge = true; break; }
        }
        iterationsRun = iter;
        converged = didConverge;
        LinFit raw = toRawFit(cur, scaler);   // convert coefficients/intercept back to raw units
        intercept = raw.intercept; coef = raw.coef; weights = w;
        nFeatures = d; fitted = true;
    }
    // predict(x) / predict(X): identical shape to every other regressor in the series.
};

Quantile Regression and the Pinball Loss

Ordinary regression estimates $E[y\mid x]$, the conditional mean. Quantile regression at level $\tau\in(0,1)$ estimates the conditional $\tau$-quantile instead, the value below which a fraction $\tau$ of the conditional distribution of $y$ falls. The loss that makes this estimator fall out of a minimization, the pinball (or check) loss, penalizes a positive residual by $\tau$ per unit and a negative residual by $(1-\tau)$ per unit:

The second form (half the absolute value plus a tilt term) makes the asymmetry easy to read off: at $\tau=0.5$ the tilt term vanishes and $\rho_{0.5}(r)=\tfrac12|r|$, exactly (half of) absolute error, so LAD is quantile regression at $\tau=0.5$ and needs no separate implementation. For $\tau\ne 0.5$, under-predicting ($r=y-\hat y>0$) and over-predicting ($r<0$) are penalized at different rates, $\tau$ versus $1-\tau$, which is exactly what pulls the fit toward the $\tau$-quantile rather than the mean: a large $\tau$ (say $0.9$) penalizes under-prediction nine times more than over-prediction, so the fit rises until only about $10\%$ of points remain above the line.

Differentiating $\sum_i \rho_\tau(r_i)$ with $r_i = y_i - x_i^\top\beta$ gives the subgradient $\psi_\tau(r) = \tau - \mathbb{1}[r<0]$ (equivalently $\tau$ when $r\ge 0$ and $\tau - 1$ when $r<0$), the standard result for the check-function derivative. The same weight trick as Huber applies: writing $\psi_\tau(r) = c(r)\cdot r / |r|\cdot|r|$… more directly, define $c(r) = \tau$ if $r\ge 0$ else $1-\tau$, and observe $\psi_\tau(r) = c(r)\cdot\operatorname{sign}(r)$. Setting $w(r) = c(r)/\max(|r|,\varepsilon)$ (with a small floor $\varepsilon$ guarding only against an exactly-zero residual) gives $w(r)\cdot r = c(r)\cdot r/|r| = c(r)\cdot\operatorname{sign}(r) = \psi_\tau(r)$ for every $r\ne 0$, so once again the true objective's stationarity condition coincides with the stationarity condition of a weighted squared-error problem, $\sum_i w(r_i)\,r_i^2$, and the same IRLS alternation applies. At $\tau=0.5$, $c(r)=0.5$ for every point regardless of sign, and this reduces to the classical IRLS scheme for LAD.

A single check confirms the asymmetry is doing what it should. For $\tau=0.75$: a point with $r=+2$ (the fit under-predicted) gets $c(r)=0.75$ and weight $w=0.75/2=0.375$; a point with $r=-2$ (over-predicted, same magnitude) gets $c(r)=0.25$ and weight $w=0.25/2=0.125$, one-third as large. Under-prediction is punished three times harder than an equally large over-prediction at $\tau=0.75$, exactly $\tau/(1-\tau)=0.75/0.25=3$, which is what pulls the fitted line up until it sits above about $75\%$ of the points.

IRLS by Hand: The 90th-Percentile Line Wants the Outlier

Reusing the identical six-point dataset from the Huber walkthrough makes a sharp point: robustness and quantile-seeking are different goals, and the same "outlier" is treated oppositely by each. Fitting $\tau=0.9$ from the same OLS point $\beta=(2.6667,\ 3.4286)$ proceeds as follows:

1. Iteration 1. Residuals $r=[1.9048,\,0.4762,\,-0.9524,\,-2.3810,\,-3.8095,\,4.7619]$. Since $\tau=0.9$, every non-negative residual gets $c=0.9$ and every negative residual gets $c=0.1$: $c=[0.9,\,0.9,\,0.1,\,0.1,\,0.1,\,0.9]$. Dividing by $|r_i|$ gives weights $[0.4725,\,1.8900,\,0.1050,\,0.0420,\,0.0263,\,0.1890]$: point $2$, which sits just above the line, gets the largest weight of all, while every point below the line is suppressed. Re-solving gives $\beta=(4.0855,\ 3.8431)$: the line jumps up, not down.
2. Iteration 2. $r=[1.5224,\,-0.3208,\,-2.1639,\,-4.0070,\,-5.8502,\,2.3067]$; $c=[0.9,\,0.1,\,0.1,\,0.1,\,0.1,\,0.9]$; weights $[0.5912,\,0.3118,\,0.0462,\,0.0250,\,0.0171,\,0.3902]$. The line keeps rising toward the points that remain above it. $\beta=(5.2064,\ 3.9985)$.

Continuing this alternation (verified by script) drives $\beta$ toward $(6,\ 4)$, the exact minimizer of the pinball objective on this data, a reference optimum fixed by the geometry rather than by the solver: with two parameters the pinball loss is piecewise linear in $\beta$, so its minimizer sits at a vertex of that surface, which generically means the fitted line interpolates exactly two training points. Here that vertex is the line through point $1$ ($x=-2.5$: $6+4(-2.5)=-4=y_1$) and the corrupted point $6$ ($x=2.5$: $6+4(2.5)=16=y_6$). The $\varepsilon$-floored IRLS only approaches this vertex; its coefficient-change stopping rule reports that the iterate stopped moving, not a certificate that the linear program's optimum was reached exactly (the implementation section makes this approximation precise), so $(6,\ 4)$ is best read as the limiting reference value the iteration tracks rather than a guaranteed-exact return. Huber's IRLS pushed the outlier's weight toward zero because its job is to resist being pulled by extreme points; quantile IRLS at $\tau=0.9$ gives the same point elevated influence, because being asked for the line below which $90\%$ of points should fall is, definitionally, a request to pay the most attention to the points near the top.

QuantileRegressor Implementation

struct QuantileRegressor {
    bool fitted = false;
    int nFeatures = 0;
    double intercept = 0.0;
    vector<double> coef;
    double tau = 0.5;
    int iterationsRun = 0;
    bool converged = false;
    vector<double> weights;

    void fit(const Matrix& X, const vector<double>& y, double tauIn,
             int maxIter = 1000, double tol = 1e-9) {
        fitted = false;
        validateFitInput(X, y, "QuantileRegressor");
        if (!isfinite(tauIn) || tauIn <= 0.0 || tauIn >= 1.0)
            throw invalid_argument("QuantileRegressor::fit: tau must lie strictly between 0 and 1");
        tau = tauIn;
        int n = (int)X.size(), d = (int)X[0].size();

        double yScale = 1.0;
        for (double v : y) yScale = max(yScale, fabs(v));
        const double epsFloor = yScale * 1e-9;   // floor on the raw residual only

        Matrix Z;
        Standardizer scaler = standardizeFeatures(X, Z);
        LinFit cur = fitWeightedOLS(Z, y, vector<double>(n, 1.0));
        vector<double> w(n, 1.0);
        int iter = 0;
        bool didConverge = false;
        for (; iter < maxIter; ++iter) {
            vector<double> r = residualsOf(Z, y, cur);
            for (int i = 0; i < n; ++i) {
                double c = (r[i] >= 0.0) ? tau : (1.0 - tau);
                w[i] = c / max(fabs(r[i]), epsFloor);
            }
            LinFit next = fitWeightedOLS(Z, y, w);
            double maxDelta = fabs(next.intercept - cur.intercept);
            for (size_t j = 0; j < next.coef.size(); ++j) maxDelta = max(maxDelta, fabs(next.coef[j] - cur.coef[j]));
            cur = next;
            if (maxDelta < tol) { ++iter; didConverge = true; break; }
        }
        iterationsRun = iter;
        converged = didConverge;
        LinFit raw = toRawFit(cur, scaler);
        intercept = raw.intercept; coef = raw.coef; weights = w;
        nFeatures = d; fitted = true;
    }
    // predict(x) / predict(X): identical shape to every other regressor in the series.
};

The only structural difference from HuberRegressor::fit is the weight formula: a fixed threshold $\delta$ and a symmetric cap versus a per-row asymmetric factor $c(r)$ and a division floored at $\varepsilon$. Everything else, the train-only standardization, the OLS initialization, the weighted-refit loop, the max-coefficient-change convergence test, and the raw-unit conversion, is identical in shape to the Huber code the previous section already validated.

That $\varepsilon$ floor is not a cosmetic guard; it is what makes this routine an IRLS approximation to pinball-loss minimization rather than an exact solver for it. Minimizing $\sum_i\rho_\tau(r_i)$ is a linear program whose optimum is characterized by an LP optimality (KKT) certificate: at the solution the fitted hyperplane rests against a subset of points whose subgradients balance and $d{+}1$ residuals are exactly zero. The subgradient weight $c(r)/\lvert r\rvert$ diverges as those residuals approach zero, so flooring it at $\varepsilon$ replaces each true infinite weight with a large finite one; the fixed point of the floored iteration minimizes a smoothed surrogate that only approaches the true pinball optimum as $\varepsilon\to 0$. The stopping rule compounds this: convergence here means the coefficient vector moved by less than tol between iterations, a Cauchy-style numerical criterion, not a verified LP/KKT optimality certificate. The returned coefficients therefore track the pinball optimum closely but are not guaranteed correct to full precision, and the fit is sensitive to $\varepsilon$, to the OLS initialization, and to tol, since tightening or loosening any of them shifts both the iterate the loop stops at and the number of iterations it takes. The per-$\tau$ pinball loss printed by the complete program is included precisely so a lower value from a reference solver can expose the gap; for production work where the exact conditional quantile and a genuine optimality certificate matter, a dedicated linear-program or interior-point quantile solver is the right tool, with this IRLS scheme serving as a transparent, from-scratch approximation.

The Delta Continuum: Huber Between OLS and LAD

Huber's $\delta$ is not just a tuning knob, it interpolates between two estimators this post has now fully implemented. As $\delta\to\infty$, the condition $|r_i|\le\delta$ holds for every finite residual, every weight is exactly $1$, and IRLS converges in a single iteration to precisely the unweighted OLS fit. As $\delta\to 0$, almost every nonzero residual exceeds $\delta$, and the weight $\delta/|r_i|$ becomes proportional to $1/|r_i|$ for every point, exactly the LAD weighting up to an overall constant that cancels in the stationarity condition. On the contaminated dataset used in the complete program's Scenario A, sweeping $\delta$ over five decades makes this continuum concrete (all six rows share the same training data; only $\delta$ changes):

At $\delta=50$ and $\delta=5000$ the fit matches the OLS coefficients ($2.8391$, $[2.0571,-1.5649]$) printed later in the complete program's output to every displayed digit, confirming the $\delta\to\infty$ limit lands on the unweighted estimator (the two paths differ only in the last floating-point bits, since the robust fit reaches the OLS coefficients through the standardized solve rather than the raw one), and each converges in a single iteration because every weight is $1$ on the very first pass. Choosing $\delta$ in practice is a statement about how large a residual is still plausibly ordinary noise rather than contamination: a common rule of thumb is a small multiple (around $1.345$) of a robust estimate of the noise scale, such as the median absolute deviation of the residuals, though the complete program below simply fixes $\delta=1.0$ for a noise scale of $\sigma=0.5$.

The complete program turns that rule of thumb into code rather than leaving it as advice. It estimates the noise scale from the OLS residuals with the median absolute deviation, $\hat\sigma=\operatorname{median}_i\lvert r_i-\operatorname{median}(r)\rvert/0.67449$, where the $0.67449$ divisor rescales the MAD so it estimates the standard deviation of Gaussian noise. Because the MAD discards everything past the median, it stays near the true $\sigma=0.5$ even with $8\%$ of the residuals corrupted: the program reports $\hat\sigma=0.5895$ and the implied $\delta=1.345\hat\sigma=0.7929$, and the Huber fit at that data-driven $\delta$ lands at $(2.9717,\ [1.9545,-1.5419])$, essentially the fixed-$\delta$ result. The MAD carries one degenerate case worth guarding: when more than half the residuals coincide the MAD is exactly zero, so the scale routine falls back to the mean absolute deviation when that value is positive. A residual vector with no dispersion at all is rejected explicitly because the data contains no positive noise scale from which to infer $\delta$. The program checks the non-degenerate zero-MAD path on $[0,0,0,0,2,-2]$, whose fallback scale is $0.6667$. This is also where the estimator's contract is worth stating outright: although HuberRegressor standardizes its features internally (train-only) and converts the fit back to raw units, it never rescales $y$, so each $\lvert r_i\rvert$ is compared directly against $\delta$ in raw $y$ units and $\delta$, $\hat\sigma$, and the residuals all share the scale of $y$. The internal feature standardization is the same train-only device the penalized estimators in the previous post rely on, used here purely to condition the repeatedly re-solved weighted normal equations rather than to make a penalty scale-fair.

Validation and Errors

• Row count mismatch, ragged rows, empty input, or non-finite entries throw std::invalid_argument, exactly as in every earlier post.
• HuberRegressor::fit rejects a non-finite or non-positive delta.
• madScale rejects an empty or fully zero-dispersion residual vector; neither case contains enough information to infer a positive robust noise scale.
• QuantileRegressor::fit rejects a tau outside the open interval $(0,1)$; $\tau=0$ or $\tau=1$ would ask for an infinite quantile of a continuous distribution and is refused rather than silently clamped.
• Non-convergence within maxIter is not an exception, following the same convention as gradient descent and elastic net: iterationsRun and converged are exposed for the caller to inspect.
• A non-finite intercept or coefficient at the end of fitting, which would indicate a genuine solver failure, throws std::runtime_error.
• Calling predict before fit throws std::logic_error; a feature-count mismatch or non-finite feature throws std::invalid_argument.

Complete Program and Output

The program below runs two scenarios. Scenario A trains on a two-feature dataset ($n=150$) with $8\%$ of rows corrupted by a large additive outlier (magnitude $8$-$15$, random sign) and evaluates OLS, Huber ($\delta=1.0$), and LAD on a separate, clean $60$-row test set from the same true model, so the comparison measures recovery of the true relationship rather than fit to contaminated data. It also sweeps $\delta$ across six decades and reports a MAD-based noise scale with the data-driven $\delta$ it implies. Scenario B fits a single-feature heteroskedastic dataset ($n=300$, noise standard deviation growing linearly with $x$) at five quantile levels, then checks the empirical coverage of the resulting conditional quantile bands on a separate deterministic $300$-row heteroskedastic holdout set and counts how many holdout rows show a quantile-order crossing.

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <limits>
#include <numeric>
#include <random>
#include <stdexcept>
#include <string>
#include <vector>
using namespace std;
using Matrix = vector<vector<double>>;

// ============================================================
// Shared validation and metrics (same conventions as Parts 1-4).
// ============================================================
static void validateFitInput(const Matrix& X, const vector<double>& y, const string& who) {
    if (X.size() != y.size())
        throw invalid_argument(who + "::fit: X has " + to_string(X.size()) +
                               " rows but y has " + to_string(y.size()) + " entries");
    int n = (int)X.size();
    if (n == 0) throw invalid_argument(who + "::fit: empty training set");
    int d = (int)X[0].size();
    if (d == 0) throw invalid_argument(who + "::fit: zero feature columns");
    for (int i = 0; i < n; ++i) {
        if ((int)X[i].size() != d)
            throw invalid_argument(who + "::fit: row " + to_string(i) +
                                   " has " + to_string(X[i].size()) +
                                   " columns, expected " + to_string(d));
        for (int j = 0; j < d; ++j)
            if (!isfinite(X[i][j]))
                throw invalid_argument(who + "::fit: non-finite feature in row " + to_string(i));
        if (!isfinite(y[i]))
            throw invalid_argument(who + "::fit: non-finite target at row " + to_string(i));
    }
    if (n < d + 1)
        throw invalid_argument(who + "::fit: need at least " + to_string(d + 1) +
                               " rows for " + to_string(d) + " features plus intercept, got " +
                               to_string(n));
}

static void validateMetricInputs(const vector<double>& yTrue, const vector<double>& yPred, const string& name) {
    if (yTrue.size() != yPred.size()) throw invalid_argument(name + ": size mismatch");
    if (yTrue.empty()) throw invalid_argument(name + ": vectors must be non-empty");
    for (size_t i = 0; i < yTrue.size(); ++i)
        if (!isfinite(yTrue[i]) || !isfinite(yPred[i]))
            throw invalid_argument(name + ": non-finite value at index " + to_string(i));
}

static double mse(const vector<double>& yTrue, const vector<double>& yPred) {
    validateMetricInputs(yTrue, yPred, "mse");
    double s = 0.0;
    for (size_t i = 0; i < yTrue.size(); ++i) { double e = yTrue[i]-yPred[i]; s += e*e; }
    return s / yTrue.size();
}
static double rmse(const vector<double>& yTrue, const vector<double>& yPred) { return sqrt(mse(yTrue, yPred)); }
static double r2(const vector<double>& yTrue, const vector<double>& yPred) {
    validateMetricInputs(yTrue, yPred, "r2");
    double meanY = accumulate(yTrue.begin(), yTrue.end(), 0.0) / yTrue.size();
    double ssRes = 0.0, ssTot = 0.0;
    for (size_t i = 0; i < yTrue.size(); ++i) {
        ssRes += (yTrue[i]-yPred[i])*(yTrue[i]-yPred[i]);
        ssTot += (yTrue[i]-meanY)*(yTrue[i]-meanY);
    }
    if (ssTot == 0.0) throw invalid_argument("r2: target values are constant, so R^2 is undefined");
    return 1.0 - ssRes / ssTot;
}
static double mae(const vector<double>& yTrue, const vector<double>& yPred) {
    validateMetricInputs(yTrue, yPred, "mae");
    double s = 0.0;
    for (size_t i = 0; i < yTrue.size(); ++i) s += fabs(yTrue[i]-yPred[i]);
    return s / yTrue.size();
}

// ============================================================
// Partial-pivoting Gaussian solver (identical to Part 1).
// ============================================================
static vector<double> solveLinearSystem(Matrix A, vector<double> b) {
    int n = (int)A.size();
    if (n == 0) throw invalid_argument("solveLinearSystem: A must be non-empty");
    if (b.size() != A.size()) throw invalid_argument("solveLinearSystem: A and b dimensions do not match");

    double matrixScale = 0.0;
    for (int i = 0; i < n; ++i) {
        if (A[i].size() != A.size())
            throw invalid_argument("solveLinearSystem: A must be square; row " +
                                   to_string(i) + " has " + to_string(A[i].size()) +
                                   " entries, expected " + to_string(n));
        for (double v : A[i]) {
            if (!isfinite(v)) throw invalid_argument("solveLinearSystem: A contains a non-finite value");
            matrixScale = max(matrixScale, fabs(v));
        }
        if (!isfinite(b[i])) throw invalid_argument("solveLinearSystem: b contains a non-finite value");
    }
    const double tol = numeric_limits<double>::epsilon() * max(1.0, matrixScale) * n * 100.0;

    for (int k = 0; k < n; ++k) {
        int piv = k;
        double best = fabs(A[k][k]);
        for (int i = k + 1; i < n; ++i) if (fabs(A[i][k]) > best) { best = fabs(A[i][k]); piv = i; }
        if (best <= tol)
            throw runtime_error("solveLinearSystem: matrix is singular or near-singular "
                                 "(no usable pivot in column " + to_string(k) + "); "
                                 "check for duplicated/collinear features");
        if (piv != k) { swap(A[k], A[piv]); swap(b[k], b[piv]); }
        for (int i = k + 1; i < n; ++i) {
            double factor = A[i][k] / A[k][k];
            if (factor == 0.0) continue;
            for (int j = k; j < n; ++j) A[i][j] -= factor * A[k][j];
            b[i] -= factor * b[k];
        }
    }
    vector<double> beta(n, 0.0);
    for (int k = n - 1; k >= 0; --k) {
        double s = b[k];
        for (int j = k + 1; j < n; ++j) s -= A[k][j] * beta[j];
        beta[k] = s / A[k][k];
    }
    return beta;
}

// ============================================================
// Weighted normal equations: (X_b^T W X_b) beta = X_b^T W y.
// With w[i] == 1 for all i this is exactly the Part-1 OLS solve;
// every IRLS iteration below is one call to this function.
// ============================================================
struct LinFit { double intercept; vector<double> coef; };

static LinFit fitWeightedOLS(const Matrix& X, const vector<double>& y, const vector<double>& w) {
    int n = (int)X.size(), d = (int)X[0].size(), p = d + 1;
    Matrix A(p, vector<double>(p, 0.0));
    vector<double> b(p, 0.0);
    for (int i = 0; i < n; ++i) {
        vector<double> row(p);
        row[0] = 1.0;
        for (int j = 0; j < d; ++j) row[j+1] = X[i][j];
        for (int a = 0; a < p; ++a) {
            b[a] += w[i] * row[a] * y[i];
            for (int c = 0; c < p; ++c) A[a][c] += w[i] * row[a] * row[c];
        }
    }
    vector<double> beta = solveLinearSystem(A, b);
    LinFit out;
    out.intercept = beta[0];
    out.coef.assign(beta.begin()+1, beta.end());
    return out;
}

static LinFit fitOLSDirect(const Matrix& X, const vector<double>& y) {
    validateFitInput(X, y, "fitOLSDirect");
    return fitWeightedOLS(X, y, vector<double>(X.size(), 1.0));
}

static vector<double> residualsOf(const Matrix& X, const vector<double>& y, const LinFit& f) {
    vector<double> r(X.size());
    for (size_t i = 0; i < X.size(); ++i) {
        double pred = f.intercept;
        for (size_t j = 0; j < X[i].size(); ++j) pred += f.coef[j] * X[i][j];
        r[i] = y[i] - pred;
    }
    return r;
}

static double predictOne(const LinFit& f, const vector<double>& x) {
    double out = f.intercept;
    for (size_t j = 0; j < x.size(); ++j) out += f.coef[j] * x[j];
    return out;
}
static vector<double> predictMany(const LinFit& f, const Matrix& X) {
    vector<double> out; out.reserve(X.size());
    for (auto& row : X) out.push_back(predictOne(f, row));
    return out;
}

// ============================================================
// Train-only feature standardization for IRLS. Feature means and
// scales are computed from THIS fit's training X only. Every
// non-constant column is mapped to zero mean / unit population
// variance; a constant column (maximum deviation exactly 0) is
// dropped from the standardized design so the weighted solve never
// goes singular on it, and its raw coefficient is fixed at exactly 0.
// y is never rescaled, so residuals (and delta) stay in raw y units.
// The unpenalized estimand is invariant to this rescaling; the point
// is conditioning: the repeatedly re-solved weighted normal equations
// Z^T W Z are far better conditioned than X^T W X.
// ============================================================
struct Standardizer {
    vector<double> mu, sigma;   // per original feature column
    vector<int> keep;           // original indices of the non-constant columns
    int d = 0;                  // original feature count
};

static Standardizer standardizeFeatures(const Matrix& X, Matrix& Z) {
    int n = (int)X.size(), d = (int)X[0].size();
    Standardizer s; s.d = d;
    s.mu.assign(d, 0.0);
    s.sigma.assign(d, 1.0);
    for (int i = 0; i < n; ++i) for (int j = 0; j < d; ++j) s.mu[j] += X[i][j];
    for (int j = 0; j < d; ++j) s.mu[j] /= n;
    vector<double> maxDev(d, 0.0);
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < d; ++j) maxDev[j] = max(maxDev[j], fabs(X[i][j] - s.mu[j]));
    for (int j = 0; j < d; ++j) {
        if (maxDev[j] == 0.0) { s.sigma[j] = 1.0; continue; }  // constant column: dropped, raw coef 0
        double acc = 0.0;
        for (int i = 0; i < n; ++i) { double z = (X[i][j] - s.mu[j]) / maxDev[j]; acc += z * z; }
        s.sigma[j] = maxDev[j] * sqrt(acc / n);                // population standard deviation
        s.keep.push_back(j);
    }
    int k = (int)s.keep.size();
    Z.assign(n, vector<double>(k));
    for (int i = 0; i < n; ++i)
        for (int c = 0; c < k; ++c) { int j = s.keep[c]; Z[i][c] = (X[i][j] - s.mu[j]) / s.sigma[j]; }
    return s;
}

// Convert a standardized-space fit (intercept plus one coefficient per KEPT
// column) back to raw units: coef_j = betaTilde_c / sigma_j for kept columns,
// exactly 0 for constant columns, with the intercept absorbing the centering.
// Predictions on raw X are identical to the standardized-space predictions.
static LinFit toRawFit(const LinFit& z, const Standardizer& s) {
    LinFit raw; raw.coef.assign(s.d, 0.0);
    raw.intercept = z.intercept;
    for (size_t c = 0; c < s.keep.size(); ++c) {
        int j = s.keep[c];
        raw.coef[j] = z.coef[c] / s.sigma[j];
        raw.intercept -= raw.coef[j] * s.mu[j];
    }
    return raw;
}

// ============================================================
// HuberRegressor: IRLS with weight w_i = 1 if |r_i| <= delta,
// else delta/|r_i|. Features are standardized (train-only) before the
// loop and the fit is mapped back to raw units afterward; each iteration
// re-solves the standardized weighted normal equations above. The
// sequence of iterates is a majorize-minimize scheme whose fixed point
// is the exact Huber M-estimator.
// ============================================================
struct HuberRegressor {
    bool fitted = false;
    int nFeatures = 0;
    double intercept = 0.0;
    vector<double> coef;
    double delta = 1.0;
    int iterationsRun = 0;
    bool converged = false;
    vector<double> weights;  // final IRLS weights, one per training row

    void fit(const Matrix& X, const vector<double>& y, double deltaIn,
             int maxIter = 200, double tol = 1e-9) {
        fitted = false;
        validateFitInput(X, y, "HuberRegressor");
        if (!isfinite(deltaIn) || deltaIn <= 0.0)
            throw invalid_argument("HuberRegressor::fit: delta must be finite and positive");
        if (maxIter <= 0) throw invalid_argument("HuberRegressor::fit: maxIter must be positive");
        if (!isfinite(tol) || tol <= 0.0) throw invalid_argument("HuberRegressor::fit: tol must be finite and positive");
        delta = deltaIn;
        int n = (int)X.size(), d = (int)X[0].size();

        Matrix Z;
        Standardizer scaler = standardizeFeatures(X, Z);  // means/scales from THIS training X only

        // IRLS runs entirely in standardized feature coordinates; residuals
        // and delta stay in raw y units because y is never rescaled.
        LinFit cur = fitWeightedOLS(Z, y, vector<double>(n, 1.0));  // unit-weight OLS start
        vector<double> w(n, 1.0);
        int iter = 0;
        bool didConverge = false;
        for (; iter < maxIter; ++iter) {
            vector<double> r = residualsOf(Z, y, cur);
            for (int i = 0; i < n; ++i) {
                double a = fabs(r[i]);
                w[i] = (a <= delta) ? 1.0 : delta / a;
            }
            LinFit next = fitWeightedOLS(Z, y, w);
            double maxDelta = fabs(next.intercept - cur.intercept);
            for (size_t j = 0; j < next.coef.size(); ++j) maxDelta = max(maxDelta, fabs(next.coef[j] - cur.coef[j]));
            cur = next;
            if (maxDelta < tol) { ++iter; didConverge = true; break; }
        }
        iterationsRun = iter;
        converged = didConverge;
        LinFit raw = toRawFit(cur, scaler);  // convert coefficients/intercept back to raw units
        if (!isfinite(raw.intercept)) throw runtime_error("HuberRegressor::fit: intercept became non-finite");
        for (double c : raw.coef) if (!isfinite(c)) throw runtime_error("HuberRegressor::fit: coefficient became non-finite");
        intercept = raw.intercept; coef = raw.coef; weights = w;
        nFeatures = d; fitted = true;
    }
    double predict(const vector<double>& x) const {
        if (!fitted) throw logic_error("HuberRegressor::predict: model is not fitted");
        if ((int)x.size() != nFeatures)
            throw invalid_argument("HuberRegressor::predict: expected " + to_string(nFeatures) +
                                   " features, got " + to_string(x.size()));
        double out = intercept;
        for (int j = 0; j < nFeatures; ++j) {
            if (!isfinite(x[j])) throw invalid_argument("HuberRegressor::predict: non-finite feature at index " + to_string(j));
            out += coef[j] * x[j];
        }
        if (!isfinite(out)) throw runtime_error("HuberRegressor::predict: prediction became non-finite");
        return out;
    }
    vector<double> predict(const Matrix& X) const {
        vector<double> out; out.reserve(X.size());
        for (auto& row : X) out.push_back(predict(row));
        return out;
    }
};

// ============================================================
// QuantileRegressor: IRLS approximation to pinball-loss minimization,
// weight w_i = c_i / max(|r_i|, eps), c_i = tau if r_i >= 0 else (1 - tau).
// tau = 0.5 is LAD (median regression). Features are standardized
// (train-only) and the fit is mapped back to raw units, exactly as in
// HuberRegressor. eps is a small data-scaled floor on the raw residual;
// it makes the scheme an educational numerical approximation, not an
// exact linear-program solver (see the discussion below).
// ============================================================
struct QuantileRegressor {
    bool fitted = false;
    int nFeatures = 0;
    double intercept = 0.0;
    vector<double> coef;
    double tau = 0.5;
    int iterationsRun = 0;
    bool converged = false;
    vector<double> weights;

    void fit(const Matrix& X, const vector<double>& y, double tauIn,
             int maxIter = 1000, double tol = 1e-9) {
        fitted = false;
        validateFitInput(X, y, "QuantileRegressor");
        if (!isfinite(tauIn) || tauIn <= 0.0 || tauIn >= 1.0)
            throw invalid_argument("QuantileRegressor::fit: tau must lie strictly between 0 and 1");
        if (maxIter <= 0) throw invalid_argument("QuantileRegressor::fit: maxIter must be positive");
        if (!isfinite(tol) || tol <= 0.0) throw invalid_argument("QuantileRegressor::fit: tol must be finite and positive");
        tau = tauIn;
        int n = (int)X.size(), d = (int)X[0].size();

        double yScale = 1.0;
        for (double v : y) yScale = max(yScale, fabs(v));
        const double epsFloor = yScale * 1e-9;  // floor on the raw residual only

        Matrix Z;
        Standardizer scaler = standardizeFeatures(X, Z);  // means/scales from THIS training X only

        LinFit cur = fitWeightedOLS(Z, y, vector<double>(n, 1.0));
        vector<double> w(n, 1.0);
        int iter = 0;
        bool didConverge = false;
        for (; iter < maxIter; ++iter) {
            vector<double> r = residualsOf(Z, y, cur);
            for (int i = 0; i < n; ++i) {
                double c = (r[i] >= 0.0) ? tau : (1.0 - tau);
                w[i] = c / max(fabs(r[i]), epsFloor);
            }
            LinFit next = fitWeightedOLS(Z, y, w);
            double maxDelta = fabs(next.intercept - cur.intercept);
            for (size_t j = 0; j < next.coef.size(); ++j) maxDelta = max(maxDelta, fabs(next.coef[j] - cur.coef[j]));
            cur = next;
            if (maxDelta < tol) { ++iter; didConverge = true; break; }
        }
        iterationsRun = iter;
        converged = didConverge;
        LinFit raw = toRawFit(cur, scaler);
        if (!isfinite(raw.intercept)) throw runtime_error("QuantileRegressor::fit: intercept became non-finite");
        for (double c : raw.coef) if (!isfinite(c)) throw runtime_error("QuantileRegressor::fit: coefficient became non-finite");
        intercept = raw.intercept; coef = raw.coef; weights = w;
        nFeatures = d; fitted = true;
    }
    double predict(const vector<double>& x) const {
        if (!fitted) throw logic_error("QuantileRegressor::predict: model is not fitted");
        if ((int)x.size() != nFeatures)
            throw invalid_argument("QuantileRegressor::predict: expected " + to_string(nFeatures) +
                                   " features, got " + to_string(x.size()));
        double out = intercept;
        for (int j = 0; j < nFeatures; ++j) {
            if (!isfinite(x[j])) throw invalid_argument("QuantileRegressor::predict: non-finite feature at index " + to_string(j));
            out += coef[j] * x[j];
        }
        if (!isfinite(out)) throw runtime_error("QuantileRegressor::predict: prediction became non-finite");
        return out;
    }
    vector<double> predict(const Matrix& X) const {
        vector<double> out; out.reserve(X.size());
        for (auto& row : X) out.push_back(predict(row));
        return out;
    }
};

static double pinballLoss(const vector<double>& yTrue, const vector<double>& yPred, double tau) {
    double s = 0.0;
    for (size_t i = 0; i < yTrue.size(); ++i) {
        double r = yTrue[i] - yPred[i];
        s += (r >= 0.0) ? tau * r : (tau - 1.0) * r;
    }
    return s / yTrue.size();
}

// ============================================================
// MAD-based robust scale of a residual vector, with an explicit zero-MAD
// fallback. sigma_hat = median(|r_i - median(r)|) / 0.67449; the divisor
// rescales the MAD to estimate the standard deviation of Gaussian noise.
// ============================================================
static double medianOf(vector<double> v) {
    size_t n = v.size();
    sort(v.begin(), v.end());
    return (n % 2) ? v[n/2] : 0.5 * (v[n/2 - 1] + v[n/2]);
}

static double madScale(const vector<double>& r) {
    if (r.empty())
        throw invalid_argument("madScale: residual vector must be non-empty");
    double med = medianOf(r);
    vector<double> ad(r.size());
    for (size_t i = 0; i < r.size(); ++i) ad[i] = fabs(r[i] - med);
    double mad = medianOf(ad);
    if (mad > 0.0) return mad / 0.67449;             // usual case
    double meanAd = 0.0; for (double a : ad) meanAd += a; meanAd /= (double)ad.size();
    if (meanAd > 0.0) return meanAd;                 // zero-MAD fallback: mean absolute deviation
    throw invalid_argument(
        "madScale: residuals have zero dispersion; choose delta explicitly");
}

// ============================================================
// main
// ============================================================
int main() {
    // === Scenario A: outlier contamination, evaluated on a CLEAN held-out test set ===
    printf("=== Scenario A: 8%% contaminated training data, clean test set ===\n");
    Matrix Xtr, Xte; vector<double> ytr, yte; vector<int> outlierIdx;
    {
        mt19937 rng(11);
        uniform_real_distribution<double> ux(-3.0, 3.0);
        normal_distribution<double> noise(0.0, 0.5);
        uniform_real_distribution<double> outMag(8.0, 15.0);
        uniform_int_distribution<int> outSign(0, 1);
        const int N = 150;
        for (int i = 0; i < N; ++i) {
            double x1 = ux(rng), x2 = ux(rng);
            Xtr.push_back({x1, x2});
            ytr.push_back(3.0 + 2.0*x1 - 1.5*x2 + noise(rng));
        }
        int nOut = (int)(0.08 * N);
        // Deterministically pick which rows become outliers via a shuffled index list.
        vector<int> order(N); iota(order.begin(), order.end(), 0);
        shuffle(order.begin(), order.end(), rng);
        for (int k = 0; k < nOut; ++k) {
            int i = order[k];
            double mag = outMag(rng);
            double sign = outSign(rng) ? 1.0 : -1.0;
            ytr[i] += sign * mag;
            outlierIdx.push_back(i);
        }
        sort(outlierIdx.begin(), outlierIdx.end());

        const int NT = 60;
        for (int i = 0; i < NT; ++i) {
            double x1 = ux(rng), x2 = ux(rng);
            Xte.push_back({x1, x2});
            yte.push_back(3.0 + 2.0*x1 - 1.5*x2 + noise(rng));  // clean, no contamination
        }
    }
    printf("contaminated rows (%zu of %zu): ", outlierIdx.size(), Xtr.size());
    for (size_t i = 0; i < outlierIdx.size(); ++i) printf("%s%d", i ? "," : "", outlierIdx[i]);
    printf("\n");

    LinFit olsA = fitOLSDirect(Xtr, ytr);
    printf("[OLS]   intercept=%.4f  coef=[%.4f, %.4f]\n", olsA.intercept, olsA.coef[0], olsA.coef[1]);
    printf("[OLS]   test RMSE=%.4f  test MAE=%.4f  test R^2=%.4f\n",
           rmse(yte, predictMany(olsA, Xte)), mae(yte, predictMany(olsA, Xte)), r2(yte, predictMany(olsA, Xte)));

    HuberRegressor huber; huber.fit(Xtr, ytr, 1.0);
    printf("[Huber] intercept=%.4f  coef=[%.4f, %.4f]  (delta=1.0, iters=%d, converged=%s)\n",
           huber.intercept, huber.coef[0], huber.coef[1], huber.iterationsRun, huber.converged ? "true" : "false");
    printf("[Huber] test RMSE=%.4f  test MAE=%.4f  test R^2=%.4f\n",
           rmse(yte, huber.predict(Xte)), mae(yte, huber.predict(Xte)), r2(yte, huber.predict(Xte)));
    double meanWOutlier = 0.0, meanWClean = 0.0;
    {
        int oi = 0; double sumOut = 0.0, sumClean = 0.0; int nOut = 0, nClean = 0;
        for (int i = 0; i < (int)huber.weights.size(); ++i) {
            bool isOut = (oi < (int)outlierIdx.size() && outlierIdx[oi] == i);
            if (isOut) { sumOut += huber.weights[i]; ++nOut; ++oi; }
            else { sumClean += huber.weights[i]; ++nClean; }
        }
        meanWOutlier = sumOut / nOut; meanWClean = sumClean / nClean;
    }
    printf("[Huber] mean IRLS weight: contaminated rows=%.4f, clean rows=%.4f\n", meanWOutlier, meanWClean);

    QuantileRegressor lad; lad.fit(Xtr, ytr, 0.5);
    printf("[LAD]   intercept=%.4f  coef=[%.4f, %.4f]  (iters=%d, converged=%s)\n",
           lad.intercept, lad.coef[0], lad.coef[1], lad.iterationsRun, lad.converged ? "true" : "false");
    printf("[LAD]   test RMSE=%.4f  test MAE=%.4f  test R^2=%.4f\n",
           rmse(yte, lad.predict(Xte)), mae(yte, lad.predict(Xte)), r2(yte, lad.predict(Xte)));

    // Delta continuum: Huber -> LAD as delta shrinks, Huber -> OLS as delta grows.
    printf("\n-- Huber delta continuum on the same contaminated training data --\n");
    for (double d0 : {0.05, 0.2, 1.0, 5.0, 50.0, 5000.0}) {
        HuberRegressor h; h.fit(Xtr, ytr, d0, 500);
        printf("  delta=%-8g intercept=%.4f  coef=[%.4f, %.4f]  iters=%d\n", d0, h.intercept, h.coef[0], h.coef[1], h.iterationsRun);
    }

    // MAD-based noise scale of the OLS residuals and the rule-of-thumb Huber delta
    // it implies. HuberRegressor standardizes its features internally but never
    // rescales y, so delta, sigma_hat, and the residuals all share y's raw scale.
    {
        vector<double> rOls = residualsOf(Xtr, ytr, olsA);
        double sigmaHat = madScale(rOls);
        double deltaRot = 1.345 * sigmaHat;
        HuberRegressor hMad; hMad.fit(Xtr, ytr, deltaRot);
        printf("\n-- MAD-based scale and the Huber delta it suggests --\n");
        printf("sigma_hat (MAD/0.6745 of OLS residuals) = %.4f; rule-of-thumb delta = 1.345*sigma_hat = %.4f\n", sigmaHat, deltaRot);
        printf("[Huber, MAD delta] intercept=%.4f  coef=[%.4f, %.4f]  (delta=%.4f, iters=%d)\n",
               hMad.intercept, hMad.coef[0], hMad.coef[1], deltaRot, hMad.iterationsRun);
        vector<double> degenerate = {0.0, 0.0, 0.0, 0.0, 2.0, -2.0};
        printf("zero-MAD guard: madScale([0,0,0,0,2,-2]) = %.4f (MAD=0, falls back to mean absolute deviation)\n", madScale(degenerate));
    }

    // === Scenario B: heteroskedastic noise, quantile fan / conditional quantile bands ===
    printf("\n=== Scenario B: heteroskedastic noise, quantile fan (n=300) ===\n");
    Matrix Xb; vector<double> yb;
    {
        mt19937 rng(202);
        uniform_real_distribution<double> ux(0.0, 10.0);
        normal_distribution<double> stdNoise(0.0, 1.0);
        const int N = 300;
        for (int i = 0; i < N; ++i) {
            double x = ux(rng);
            double sigma = 0.3 + 0.35 * x;
            double yy = 2.0 + 1.0 * x + sigma * stdNoise(rng);
            Xb.push_back({x});
            yb.push_back(yy);
        }
    }
    vector<double> taus = {0.1, 0.25, 0.5, 0.75, 0.9};
    vector<QuantileRegressor> fits;
    for (double t : taus) {
        QuantileRegressor q; q.fit(Xb, yb, t);
        fits.push_back(q);
        printf("[tau=%.2f] intercept=%.4f  slope=%.4f  iters=%d  converged=%s  pinball=%.4f\n",
               t, q.intercept, q.coef[0], q.iterationsRun, q.converged ? "true" : "false",
               pinballLoss(yb, q.predict(Xb), t));
    }
    // Coverage and crossings are measured on a separate deterministic heteroskedastic
    // holdout set (same true model, disjoint RNG stream, seed 909), never on the
    // training rows the quantile lines were fit to, so the numbers report genuine
    // out-of-sample calibration of the conditional quantile bands, not in-sample fit.
    Matrix Xh; vector<double> yh;
    {
        mt19937 rng(909);
        uniform_real_distribution<double> ux(0.0, 10.0);
        normal_distribution<double> stdNoise(0.0, 1.0);
        const int M = 300;
        for (int i = 0; i < M; ++i) {
            double x = ux(rng);
            double sigma = 0.3 + 0.35 * x;
            double yy = 2.0 + 1.0 * x + sigma * stdNoise(rng);
            Xh.push_back({x});
            yh.push_back(yy);
        }
    }
    auto predAt = [&](int fitIdx, double x)->double { vector<double> v{x}; return fits[fitIdx].predict(v); };
    int lo10 = 0, hi90 = 4, lo25 = 1, hi75 = 3;
    int inside80 = 0, inside50 = 0, crossingRows = 0;
    for (size_t i = 0; i < Xh.size(); ++i) {
        double x = Xh[i][0];
        double lo = predAt(lo10, x), hi = predAt(hi90, x);
        if (yh[i] >= lo && yh[i] <= hi) ++inside80;
        double lo2 = predAt(lo25, x), hi2 = predAt(hi75, x);
        if (yh[i] >= lo2 && yh[i] <= hi2) ++inside50;
        // Crossing convention: a holdout row counts once if the five fitted predictions
        // at that row's x, read in increasing tau, are not non-decreasing, i.e. some
        // adjacent higher-tau line predicts strictly below its lower-tau neighbour.
        bool crossed = false;
        for (size_t k = 0; k + 1 < fits.size(); ++k)
            if (predAt(static_cast<int>(k + 1), x) < predAt(static_cast<int>(k), x)) { crossed = true; break; }
        if (crossed) ++crossingRows;
    }
    printf("holdout empirical coverage (n=%d): [tau=0.10,0.90] band=%.1f%% (target 80%%), [tau=0.25,0.75] band=%.1f%% (target 50%%)\n",
           static_cast<int>(Xh.size()), 100.0 * inside80 / Xh.size(), 100.0 * inside50 / Xh.size());
    printf("holdout quantile-order crossings: %d of %d rows have >=1 adjacent-tau crossing\n",
           crossingRows, static_cast<int>(Xh.size()));
    printf("conditional quantile band width [tau=0.10,0.90] at x=0: %.3f   at x=10: %.3f\n",
           predAt(hi90, 0.0) - predAt(lo10, 0.0), predAt(hi90, 10.0) - predAt(lo10, 10.0));

    // Error-handling demonstrations.
    printf("\n-- Error handling --\n");
    try { HuberRegressor bad; bad.fit(Xb, yb, -1.0); }
    catch (const invalid_argument& e) { printf("caught: %s\n", e.what()); }
    try { QuantileRegressor bad; bad.fit(Xb, yb, 1.0); }
    catch (const invalid_argument& e) { printf("caught: %s\n", e.what()); }
    try { QuantileRegressor bad; bad.fit(Xb, yb, 0.0); }
    catch (const invalid_argument& e) { printf("caught: %s\n", e.what()); }
    try { HuberRegressor m; m.fit(Xb, yb, 1.0); m.predict(vector<double>{1.0, 2.0}); }
    catch (const invalid_argument& e) { printf("caught: %s\n", e.what()); }
    try { QuantileRegressor m; m.predict(vector<double>{1.0}); }
    catch (const logic_error& e) { printf("caught: %s\n", e.what()); }

    return 0;
}

Compiled with g++ -O2 -std=c++17 -Wall -Wextra -pedantic, this prints:

Reading the Two Scenarios

Scenario A shows the payoff at realistic scale, not just in the six-point hand example. With $12$ of $150$ training rows corrupted, OLS's coefficients ($2.8391$, $[2.0571,-1.5649]$) are visibly off from the true $(3.0,\ 2.0,\ -1.5)$, and its clean-test-set RMSE is $0.5664$. Huber ($\delta=1.0$) recovers $(2.9701,\ [1.9538,-1.5428])$, much closer to the truth, with test RMSE $0.5167$, an $8.8\%$ reduction, and its IRLS weights make the mechanism visible directly: the $12$ contaminated rows average weight $0.0923$ while the $138$ clean rows average weight $0.9898$, a clean, almost binary separation between "trust this row" and "barely trust this row." LAD lands even closer to the true coefficients ($2.9627$, $[1.9819,-1.5205]$) and posts a lower test MAE ($0.4054$ vs. Huber's $0.4105$), yet its test RMSE is very slightly higher ($0.5173$ vs. $0.5167$). That split is not noise: LAD directly minimizes absolute error, so it should be expected to win on MAE and on the coefficients themselves, while Huber blends absolute and squared error, so it keeps a slight edge on the squared-error-based RMSE. Each estimator looks best under the metric closest to its own training objective.

Scenario B turns the same machinery into a set of calibrated conditional quantile bands rather than a single robust point estimate. Fitting $\tau\in\{0.1,0.25,0.5,0.75,0.9\}$ on the heteroskedastic training set produces a monotonically increasing slope, from $0.4844$ at $\tau=0.1$ to $1.5265$ at $\tau=0.9$, a direct, quantitative fingerprint of heteroskedasticity: if the noise had constant variance, every quantile's slope would agree with the mean regression's slope, not spread out like this. Coverage is then measured out-of-sample, on a separate deterministic $300$-row heteroskedastic holdout set the lines were never fit to: the band between the $\tau=0.1$ and $\tau=0.9$ fitted lines contains $79.0\%$ of the holdout points, and the band between $\tau=0.25$ and $\tau=0.75$ contains $47.7\%$, both within a couple of points of their theoretical targets of $80\%$ and $50\%$, which is the strongest evidence available that quantile regression is doing exactly what it claims: producing conditional quantile bands with the coverage their quantile levels promise, not merely bands that happen to look reasonable. The same holdout also quantifies the one structural weakness of fitting each $\tau$ independently: $13$ of the $300$ holdout rows fall in the narrow left-edge strip where the fitted lines cross out of quantile order, counted as any holdout row whose five predictions, read in increasing $\tau$, are not non-decreasing. The conditional quantile band width also does the right thing structurally: it grows from $0.241$ at $x=0$ to $10.662$ at $x=10$, tracking the widening true noise scale ($\sigma(x)=0.3+0.35x$ grows about $12.7\times$ over that range) far better than a single constant-width band from an OLS-plus-Gaussian-assumption approach ever could.

When to Use Each

Complexity and Memory

Unlike coordinate descent's per-sweep cost of $O(nd)$, each IRLS iteration here re-solves a full $(d{+}1)\times(d{+}1)$ system, $O(nd^2)$ for assembly plus $O(d^3)$ for the solve, because it is rebuilding the weighted normal equations from scratch every time rather than updating one coordinate. This is the right trade for the modest $d$ typical of a regression problem: IRLS needs far fewer iterations ($8$-$95$ in the scenarios above) than coordinate descent needed sweeps in the previous post's collinear example, because each IRLS iteration solves its weighted subproblem exactly rather than taking one coordinate step.

Practical Diagnostics

• The final IRLS weights, inspected directly. HuberRegressor::weights and QuantileRegressor::weights are ordinary training-row-indexed vectors; sorting them and looking at the smallest handful is a direct, model-based way to flag candidate outliers or high-leverage points, no separate diagnostic machinery required.
• The delta continuum. Refitting HuberRegressor across a wide range of $\delta$ and watching the coefficients move (as the table above does) reveals how sensitive a given dataset is to the choice: a coefficient path that is flat across most of the range and only moves near the extremes means $\delta$ barely matters; one that moves continuously means the contamination level and $\delta$ are close, and the choice deserves more care.
• Quantile monotonicity. Fitted slopes (and predictions at any fixed $x$) should be non-decreasing in $\tau$; a "crossing" where a higher $\tau$'s prediction falls below a lower $\tau$'s at some $x$ signals either too few points, too extreme a $\tau$, or a linear model that is a poor fit to genuinely nonlinear conditional quantiles.
• Empirical coverage as a cross-check. The coverage computation in Scenario B, the fraction of held-out points falling inside a $[\tau_{\text{lo}},\tau_{\text{hi}}]$ conditional quantile band, is measured on a separate deterministic holdout set rather than the training rows, a cheap and direct way to confirm a quantile-based band is calibrated out-of-sample rather than merely plausible-looking; a coverage far from $\tau_{\text{hi}}-\tau_{\text{lo}}$ indicates too few data points, a poorly chosen $\tau$ pair, or a linear model missing real curvature in the conditional quantiles.

Limitations and Edge Cases

• Convexity characterizes the optimum, not every iterative path. Both the Huber and pinball objectives are convex, so any certified stationary optimum is global. The Huber reweighting used here is a standard majorize-minimize iteration and converges under the full-rank, positive-weight conditions enforced by the implementation. The quantile routine is an $\varepsilon$-regularized IRLS approximation, so convexity of the target pinball objective does not by itself turn its coefficient-change stopping rule into an exact optimality certificate. Non-convex robust losses used elsewhere lose even the global-optimum property.
• Quantile IRLS is an approximation to the pinball linear program, not an exact solver for it. A training point that lands numerically on the fitted line uses the finite cap $c/\varepsilon$ instead of the unbounded $c/|r|$ weight implied by the IRLS representation as $r\to0$; the true pinball subdifferential at $r=0$ remains the finite interval $[\tau-1,\tau]$. The loop also stops on a coefficient-change tolerance rather than an LP/KKT optimality certificate, so the returned fit approaches but does not certify the exact conditional quantile and shifts with $\varepsilon$, the initialization, and tol. This is immaterial for the well-separated continuous targets used here, but a production quantile fit should defer to a dedicated linear-program or interior-point solver; the per-$\tau$ pinball loss printed above is the diagnostic for checking how close this approximation lands.
• Huber's $\delta$ and quantile regression's $\tau$ are both scale- and problem-dependent choices. Neither this post nor the complete program automates choosing them; a real workflow would select $\delta$ from a robust noise-scale estimate and select which $\tau$ levels to report from the decision the prediction interval is meant to support.
• None of this addresses leverage in the design matrix itself. IRLS reweights based on the size of the residual; a high-leverage point with an unusual $x$ but a small residual (because the fit bent toward it) is not downweighted by anything in this post. Leverage diagnostics are a separate concern, covered in the next post.
• Linear quantile fits can cross, and this dataset shows it. Because each $\tau$ is fit independently, nothing enforces $\hat y_{\tau_1}(x) \le \hat y_{\tau_2}(x)$ for $\tau_1<\tau_2$ at every $x$. Checking all ten pairs of the five fitted lines in Scenario B against the actual $x\in[0,10]$ range: the median ($\tau=0.5$) line dips below both the $\tau=0.1$ and $\tau=0.25$ lines for $x\lesssim 0.36$ (crossing $\tau=0.1$ at $x\approx0.160$ and $\tau=0.25$ at $x\approx0.358$), before the expected ordering reasserts itself for every $x\gtrsim0.36$; every other pair only crosses at a negative $x$, outside the training range entirely. The program makes this concrete by counting, on the same deterministic holdout set used for coverage, how many rows land in that strip: $13$ of the $300$ holdout rows have predictions that are not non-decreasing in $\tau$. The affected region is a narrow strip near the left edge of the data, where the two smaller-quantile lines happen to be close together and least well separated, but it is a real violation, not a hypothetical one. Independent per-$\tau$ optimization supplies no cross-quantile ordering constraint, so monotonicity must be checked or imposed separately.

Next Step

Every post in this series has now covered a different way to estimate $\beta$: closed-form, numerically robust, iterative, penalized, and loss-robust. What has not been covered systematically is how to check any of these fits after the fact, residual patterns, leverage, influence, and which of $R^2$, adjusted $R^2$, or something else is the honest number to report. The next post closes the series with exactly that shared diagnostic layer, applicable to every regressor built across all five previous parts.