Squared, Absolute, and Quantile Loss Compared by Their Gradients

Every regression loss is a function of the residual $r = y - \hat y$ that defines two things: the quantity the fit converges to, and the influence each observation has on the fitted values. The first is a property of where the loss is minimised. The second is a property of the loss’s derivative, because for any gradient-based fit the update applied per observation is proportional to $-\rho’(r)$. Comparing losses therefore means comparing their derivatives, and the derivative of the squared loss and the derivative of the absolute loss differ in a way that changes the behaviour of a fit on the same data by an order of magnitude.

The derivative is the influence

For the squared loss $\rho(r) = r^2$, the derivative is $2r$: unbounded in the residual. An observation sitting a thousand units from the current fit contributes a thousand times the gradient of one sitting one unit away, so any single point with an extreme response can pull the entire solution. For the absolute loss $\rho(r) = |r|$, the derivative is $\operatorname{sign}(r)$: constant in magnitude. Every observation contributes the same influence regardless of how far off it is, so an outlier’s effect is bounded by that of any other point. This is the definition of an influence function that is bounded, and it is why the absolute loss is robust and the squared loss is not.

The two losses also converge to different functionals. The squared loss is minimised at the conditional mean, because the derivative condition $\sum (y_i - \hat y) = 0$ says the residual sum is zero. The absolute loss is minimised at the conditional median, because the derivative condition says the count of positive residuals balances the count of negative ones. For a symmetric conditional distribution the two coincide; for a skewed one they do not, and a model fitted with the squared loss on a skewed response predicts something no observation resembles.

Absolute error has a genuine obstacle. The subgradient at $r = 0$ is the interval $[-1, 1]$, so the objective is non-differentiable wherever a prediction is exact, and a plain gradient method has to choose a convention there. Taking the sign function to be zero at the origin means exact fits stop updating, which is what a subgradient method does and why it converges to a point near the optimum rather than exactly at it. The alternative is to change the algorithm rather than patch the derivative, since the absolute-loss regression problem is a linear program: with residuals split into positive and negative parts, the objective and all constraints are linear and a simplex or interior-point solver returns an exact optimum. The linear-programming formulation costs $2n$ auxiliary nonnegative variables and $n$ equality constraints, which is cheaper than it sounds for moderate $n$ and considerably more robust than iterating on a subgradient.

Huber: quadratic near zero, linear in the tails

The Huber loss joins the two:

$$\rho_\delta(r) = \begin{cases} \tfrac12 r^2 & |r| \le \delta \ \delta(|r| - \tfrac12\delta) & |r| > \delta \end{cases}$$

with derivative $r$ inside the band and $\delta \operatorname{sign}(r)$ outside. The parameter $\delta$ is a scale threshold, and the fit it defines depends on it: as $\delta \to \infty$ the loss becomes the squared loss, and as $\delta \to 0$ it becomes a scaled absolute loss. Since $\delta$ is in the units of the response, choosing it requires an estimate of the residual scale, and the estimate must not itself be computed from a squared-loss fit, which the outliers would corrupt. The usual approach is to iterate: fit, estimate scale from the median absolute deviation of the residuals, reset $\delta$ at a fixed multiple of that estimate, and refit. Median absolute deviation is consistent for the standard deviation of a normal sample after multiplication by roughly 1.4826, a factor that matters because it is the only way to compare a robust scale estimate with a conventional one.

Huber’s influence function is bounded, so it is robust, but it is not resistant in the sense of having a high breakdown point: a large enough group of identical outliers can still shift the fit arbitrarily, because each contributes the full $\delta$. The theoretical maximum breakdown point of any regression estimator with a bounded influence is $1/p$ in the number of parameters, which is why high-breakdown methods are a separate line of work rather than a loss choice.

import numpy as np

def huber(r, delta):
    a = np.abs(r)
    loss = np.where(a <= delta, 0.5 * a**2, delta * (a - 0.5 * delta))
    grad = np.clip(r, -delta, delta)      # equal to sign(r)*delta outside the band
    return loss, grad

The clipped residual is the practical definition, and it also shows that the loss’s second derivative is 1 inside the band and 0 outside: the curvature vanishes in the tails, so a residual that is currently an outlier exerts bounded force but no restoring curvature, and a Newton step built from this Hessian will not shrink the outlier’s effect.

Pinball loss targets a quantile, not a centre

Quantile regression generalises the median. For a target level $\tau \in (0,1)$, the pinball loss is

$$\rho_\tau(r) = \begin{cases} \tau r & r \ge 0 \ (\tau - 1) r & r < 0 \end{cases}$$

with subgradient $-\tau$ for positive residuals and $1-\tau$ for negative ones. The asymmetry is the entire mechanism. At $\tau = 0.9$, a positive residual costs 0.9 per unit and a negative one costs 0.1 per unit, so the fit is pushed upward until roughly 90 percent of the observations lie below the prediction. Setting $\tau = 0.5$ recovers the absolute loss up to a factor of one half.

The loss is a proper scoring rule for the quantile functional, which is why the fit converges to the conditional $\tau$-quantile rather than to some compromise between the quantiles. Fitting several levels and reading them off as a band produces a prediction interval without any distributional assumption, at the cost that the separately fitted quantiles can cross, since nothing in the separate objectives prevents the $\tau = 0.6$ fit from landing below the $\tau = 0.4$ fit on some inputs, and monotonicity has to be imposed afterwards or by a joint formulation.

import numpy as np

def pinball(r, tau):
    return np.where(r >= 0, tau * r, (tau - 1.0) * r)

Both the absolute loss and the pinball loss are evaluated with the residual in the direction of the prediction error. Sign conventions are the single most common source of a quantile fit that appears to target $1-\tau$ instead of $\tau$, and the mistake is invisible on symmetric data.

Where the choice of loss shows up

Loss selection changes more than the point estimate. With the squared loss, the fitted model is the conditional mean and the residual variance is a consistent estimator of the conditional variance; with the absolute loss, the residuals have no such interpretation and their scale is related to the conditional spread only through the spacing of the quantiles. Downstream calculations that assume homoscedastic normal errors, including confidence intervals built from residual standard errors and prediction intervals of the form $\hat y \pm 1.96 s$, are invalid after a switch, and the invalidation is silent.

The functional also decides what a comparison between models measures. Two models evaluated on mean squared error are ranked by how well they predict the mean, and a model that predicts the median accurately can score badly on that metric even when it represents the data better, because the mean is dominated by the tail the median ignores. The relationship between the loss used to fit and the loss used to evaluate is not symmetric, and metrics that decompose into the same functional are the only ones that a fit under that loss will optimise directly. This is the connection to proper scoring rules for probability forecasts: the squared loss is proper for the mean, and no loss is proper for a functional it does not target.

Log-transforming the response changes the loss’s effective target without changing the code. Squared error on $\log y$ estimates the conditional geometric mean of $y$, and exponentiating the prediction gives a biased estimate of the conditional arithmetic mean. The bias is approximately half the conditional variance of the log, and correcting for it requires an estimate of that variance. The modelling consequence is documented in link functions and Poisson regression, where the log link performs the same transformation with a different implied error structure.

Outlier handling interacts with the fitting algorithm as much as with the loss. A squared-loss fit computed by an iterative method that updates every coefficient on each pass will let one extreme observation move the iterate from the first step onward, and the solution it converges to has already been displaced. The absolute and pinball losses bound that displacement per step but slow convergence when many residuals are small, since the gradient magnitude does not shrink with the residual and the last few digits take as much work as the first.

A useful diagnostic applies regardless of which loss is fitted. The residual distribution under a robust fit separates two situations that look alike in aggregate error. If the robust residuals are symmetric and unimodal except for a small group of extreme values, those extremes are candidates for investigation as data errors, since a loss that bounds their influence has already prevented them from distorting the rest of the fit. If the entire distribution is heavy-tailed and asymmetric, the loss is describing the wrong functional, and a family of quantile fits will show it: a set of levels whose fitted curves diverge from each other as the covariate grows indicates a conditional distribution whose spread is not constant, which no single point estimate represents.

The connection to probability forecasts is worth stating, because it is what makes the loss choice a statement about the target. A model fitted by squared error estimates the mean, and the mean is the value that minimises squared error in expectation; no other loss makes a mean estimate optimal. So a forecast of the mean is only usable by a consumer whose own loss is squared error. A consumer with an asymmetric loss, one that pays more for under-provisioning than over-provisioning, as in inventory or staffing, needs a quantile, and fitting a quantile requires the pinball loss rather than a transformation of the mean estimate. For forecasts that are probabilities rather than point estimates, the same reasoning takes the form of propriety, set out in proper scoring rules, where the requirement is that the loss be minimised in expectation by the honest forecast.