Loss Functions and Convex Surrogates

Why gradient-based classification uses a surrogate

When classification accuracy is the objective, the 0-1 loss measures mistakes directly: 1 if the prediction is wrong, 0 if it is right. Accuracy is not the only possible objective (unequal error costs and probability estimates come up later in this article), but it is the natural starting point. The trouble is optimization. The 0-1 loss is piecewise constant, so its gradient is zero wherever it is defined and undefined at the boundary; ordinary gradient-based training gets no signal from it. It can be minimized directly over a small finite set of candidates, or by combinatorial search for some model classes, but exact search can become computationally difficult as the model class grows. A surrogate supplies a graded penalty that makes gradient or subgradient methods useful.

The surrogates below are written in terms of the margin \(m = y \cdot f(x)\), where \(y \in \{-1, +1\}\) is the label and \(f(x)\) is a real-valued score whose sign is the predicted class — not a probability. If \(y = +1\) and \(f(x) = 2\) the margin is \(+2\): correct, with a score two units above the boundary. If \(y = -1\) and \(f(x) = 2\) the margin is \(-2\): a mistake, with a score two units on the wrong side. A correct prediction has positive margin, a mistake has negative margin, and its magnitude measures distance from zero on the score scale. This is not a calibrated confidence level: multiplying all scores by ten preserves every predicted class. The table counts margin zero as an error for either label; an actual binary classifier needs a tie rule that assigns one class at score zero.

import numpy as np

def zero_one(m):   return (m <= 0).astype(float)
def hinge(m):      return np.maximum(0, 1 - m)
def logistic(m):   return np.logaddexp(0, -m) / np.log(2)     # base 2 so it passes through 1
def exponential(m):return np.exp(-m)
def squared(m):    return (1 - m) ** 2

print(f"{'margin':>7} {'0-1':>6} {'hinge':>8} {'logistic':>9} {'exp':>9} {'squared':>9}")
for m in (-3.0, -1.0, 0.0, 1.0, 3.0):
    a = np.array([m])
    print(f"{m:7.1f} {zero_one(a)[0]:6.1f} {hinge(a)[0]:8.3f}"
          f" {logistic(a)[0]:9.3f} {exponential(a)[0]:9.3f} {squared(a)[0]:9.3f}")
#  margin    0-1    hinge  logistic       exp   squared
#    -3.0    1.0    4.000     4.398    20.086    16.000
#    -1.0    1.0    2.000     1.895     2.718     4.000
#     0.0    1.0    1.000     1.000     1.000     1.000
#     1.0    0.0    0.000     0.452     0.368     0.000
#     3.0    0.0    0.000     0.070     0.050     4.000

All four surrogates pass through 1 at margin 0 and upper-bound the 0-1 loss. Being an upper bound is not what makes them useful — the constant function 1 is also a convex upper bound and teaches nothing. What matters is that each is convex in the score and decreasing through the boundary, so pushing it down pushes scores toward the correct sign; with unrestricted scores at each input and both labels having positive conditional probability, their expected-loss minimizers choose the more probable class. At equal probabilities either class is optimal. This property is called classification calibration; it does not guarantee the same boundary for restricted models fitted on finite data. Convexity here is in the score \(f(x)\). Composing the loss with a nonlinear parameterized score need not preserve convexity; neural-network training is generally nonconvex. What differs between the surrogates is their behavior away from the boundary, and each difference is a design decision.

• Hinge is exactly zero past margin 1. Examples with margin strictly above 1 contribute zero to the loss gradient. In a standard SVM, the support vectors determine the fitted decision function; they lie on or inside the margin, including misclassified examples.
• Logistic is never zero. Each example contributes a nonzero derivative at any finite margin, though contributions can cancel across examples. On strictly separable data, scaling a separating linear score drives the unregularized loss toward zero while its coefficients grow without bound. There is no finite minimizer; this is a statement about the objective, not a guarantee that an optimizer will run indefinitely.
• Exponential grows fastest on mistakes: 20.1 at margin \(-3\) against hinge’s 4.0. AdaBoost uses it, and persistently negative margins can receive very large weight. This helps explain AdaBoost’s sensitivity to mislabeled examples; boosting methods using other losses need separate assessment.
• Squared is the only one that increases again past margin 1: 0 at margin 1 and 4.0 at margin 3. It penalizes being confidently right. That does not make it a wrong classification loss in general — its population minimizer is \(\mathbb{E}[Y \mid x] = 2P(Y = 1 \mid x) – 1\), whose sign is the optimal decision under equal costs, and the Brier score in the last section is squared loss on probabilities. But with a restricted score model, such as a linear score, the penalty on large correct margins can move the boundary at the expense of points near it. Exercise 1 shows this.

What each loss estimates

A loss does not just rank models — it determines which summary of the possible outcomes at a given input the optimal prediction targets. Such a summary, for example a mean or a quantile, is called a functional of the conditional distribution. At the population level, with unrestricted predictions, squared error targets the conditional mean when the conditional second moment is finite. Absolute error targets a conditional median, and pinball loss at level \(0 < q < 1\) targets a conditional \(q\)-quantile; finite conditional absolute expectation suffices for these latter losses. Medians and quantiles need not be unique. A fitted model with a restricted family, regularization, and a finite sample approximates these targets rather than reproducing them exactly.

The simplest check is the case with no input \(x\): a single constant prediction \(c\) for a sample. The pinball loss charges \(q \cdot (z – c)\) when the outcome exceeds the prediction and \((1 – q) \cdot (c – z)\) when it falls short, so at \(q = 0.9\) an under-prediction costs nine times as much as an over-prediction of the same size, so a minimizing prediction is a 90th quantile of the empirical distribution.

import numpy as np
from scipy.optimize import minimize_scalar

rng = np.random.default_rng(1)
z = rng.lognormal(0, 1, size=200_000)

for q in (0.1, 0.5, 0.9):
    pinball = lambda c: np.mean(np.maximum(q * (z - c), (q - 1) * (z - c)))
    best = minimize_scalar(pinball, bounds=(0, 20), method="bounded").x
    print(f"q={q}  minimizer {best:.4f}   empirical quantile {np.quantile(z, q):.4f}")
print(f"squared-loss minimizer (the mean) {np.mean(z):.4f}")
# q=0.1  minimizer 0.2782   empirical quantile 0.2782
# q=0.5  minimizer 0.9979   empirical quantile 0.9979
# q=0.9  minimizer 3.5769   empirical quantile 3.5769
# squared-loss minimizer (the mean) 1.6438

The numerical minimizer and NumPy’s quantile estimate agree to four decimals in this large sample. This is not a general exact equality: sample pinball loss can have an interval of minimizers, and NumPy’s default quantile uses linear interpolation. For example, on the sample 0 and 10 at \(q=0.9\), pinball loss is minimized at 10, while the default np.quantile returns 9. The mean in our log-normal sample is 1.64, well above the median of 1.00 and below the 90th percentile of 3.58. On a skewed target the mean, the median, and the 90th percentile are three different answers, and squared-error training targets the mean even when a median or an upper quantile would better serve the decision.

This is also a direct route to prediction intervals: fit models at \(q = 0.05\) and \(q = 0.95\), and the two estimated conditional quantiles form a nominal 90% interval, with a \(q = 0.5\) model added if a central prediction is wanted alongside it. No distributional family has to be chosen. The nominal level is not a guarantee, though: the true 5th and 95th quantiles bracket 90% of outcomes under mild continuity conditions, but the estimated ones carry sampling and model error, and separately fitted quantiles can cross. Coverage should be checked on held-out data.

Robustness: what an outlier costs

Squared error grows quadratically with the residual, so a single far-off point can dominate the entire objective. Absolute error grows linearly. Writing the residual as outcome minus prediction, the Huber loss is quadratic near zero and linear beyond a threshold \(\delta\), combining the smooth optimization of one with a bounded derivative in the residual, which limits the effect of increasing a residual while its feature values stay fixed.

import numpy as np
from scipy.optimize import minimize_scalar

rng = np.random.default_rng(0)
x = rng.normal(size=200)
y = 2 * x + 0.5 * rng.normal(size=200)

x_out, y_out = x.copy(), y.copy()
x_out[:4], y_out[:4] = 6.0, -10.0              # 2% high-leverage outliers

def fit(loss, X, Y):
    return minimize_scalar(lambda b: np.mean(loss(Y - b * X)),
                           bounds=(-5, 5), method="bounded").x

squared  = lambda r: r ** 2
absolute = lambda r: np.abs(r)
huber    = lambda r, d=1.0: np.where(np.abs(r) <= d, 0.5 * r ** 2,
                                     d * (np.abs(r) - 0.5 * d))

print(f"{'loss':12s} {'clean':>9} {'contaminated':>13} {'shift':>8}")
for name, L in (("squared", squared), ("absolute", absolute), ("huber(1.0)", huber)):
    c, o = fit(L, x, y), fit(L, x_out, y_out)
    print(f"{name:12s} {c:9.4f} {o:13.4f} {o - c:8.4f}")
# loss             clean  contaminated    shift
# squared         1.9647        0.3736  -1.5911
# absolute        1.9297        1.8530  -0.0767
# huber(1.0)      1.9508        1.8158  -0.1350

Four contaminated points out of two hundred move the squared-error slope from 1.96 to 0.37, a shift of 1.59. Absolute error moves by 0.08 and Huber by 0.14. The contamination is 2% of the data and changes the squared-error answer by 81%.

Leverage, meaning an observation’s unusual position in the input space, also matters, and a robust loss does not remove it. The gradient of this fit with respect to the slope is proportional to \(-x \cdot \psi(y – bx)\), where \(\psi\) is the derivative of the loss. Huber bounds \(\psi\), the residual part; absolute loss has bounded slopes and a bounded choice of subgradient at zero. Both limit this residual factor, but not \(x\). These outliers were placed at \(x = 6\), far from the bulk of the data, and that factor multiplies whatever influence the residual is allowed. In this model without an intercept, the same vertical displacement at exactly \(x = 0\) would contribute nothing to the slope gradient, although swapping out the original four observations would still change the fit slightly. This is why residual plots should be read against fitted values and against each feature, not only as a histogram, and why bounded-influence claims about robust losses refer to the residual direction only.

The trade is that \(\delta\) is another hyperparameter, and that large residuals may be genuine observations rather than errors. Choose the loss according to the quantity you want to predict and how strongly large errors should influence the fit. Squared loss targets the mean, absolute loss a median, and Huber loss defines a different target through its threshold; these targets can differ, especially when the outcome distribution is asymmetric.

Asymmetric costs

Fraud screening involves unequal error costs. Missing a fraudulent transaction and blocking a legitimate one are not the same mistake. A symmetric training loss does not have to assert that they are, because estimating probabilities and making decisions are separate steps: a symmetric proper loss can estimate \(p = P(Y = 1 \mid x)\), and the cost asymmetry can be applied when the probability is turned into a decision.

Class weights change the training objective; a decision threshold changes how an existing score becomes a label. With a proper probability loss, weighting changes the population probability that the loss targets. The fitted values still depend on sample size, regularization, and the model family. In an unrestricted population calculation, suitable class weights and a threshold change can yield the same decisions, but finite fitted models need not agree.

For nonnegative costs, with a false-positive cost \(C_{\mathrm{FP}}\), a false-negative cost \(C_{\mathrm{FN}}\), and zero cost for correct decisions, predicting 1 has expected cost \((1-p)C_{\mathrm{FP}}\): you pay only when the case is actually negative. Predicting 0 has expected cost \(pC_{\mathrm{FN}}\). Comparing these costs gives \((1-p)C_{\mathrm{FP}} < pC_{\mathrm{FN}}\), or, for positive total error cost,

\[p > \frac{C_{\mathrm{FP}}}{C_{\mathrm{FP}} + C_{\mathrm{FN}}}\]

The formula uses the conditional probability for the case being decided. If a missed fraud costs 9 and a false alarm costs 1, the threshold is \(1/(1+9)=0.1\). At \(p=0.2\), flagging costs \(0.8\times1=0.8\) on average, while not flagging costs \(0.2\times9=1.8\), so flagging is cheaper. At the threshold the two actions tie. With estimated probabilities, this is an estimated-cost rule: good average calibration does not guarantee accurate probabilities for every case or subgroup.

For a fixed, known cost ratio and useful probability estimates, threshold adjustment is a simple starting point: it is one number, it can be changed without retraining, and it keeps the probability estimates available for other consumers of the model. Class weights can change which cases the fitted model prioritizes, including when rare positives have little influence on the original fit. Their value should be assessed against the application’s evaluation metric or cost. A model that outputs only a decision score, not a probability, can still have its threshold moved; what it cannot do is take the cost formula above directly, because the score is not a probability.

Choosing one

Using the convention that a score is a loss to minimize, a proper scoring rule is one whose expected value is minimized by reporting the true probability; a strictly proper one is minimized only there, so the true probability is the unique target. Log loss and Brier score are strictly proper. Binary log loss uses labels 0 and 1 and a reported probability \(r\): \(-y\log r-(1-y)\log(1-r)\). It is the earlier logistic margin loss after converting a score to \(r=1/(1+e^{-f(x)})\) and recoding the labels, apart from the log-base scale factor. Multiplying an unregularized loss by a positive constant leaves its minimizers unchanged; with a fixed additive regularization penalty, it changes the relative penalty strength. Hinge loss is not a scoring rule for probabilities, and thresholded accuracy is not strictly proper: it cannot separate the true probability from any other report on the same side of the threshold.

Calibration means that among cases predicted at 0.7, about 70% turn out positive. Training on a strictly proper loss makes the true probability the target, but a finite sample, a misspecified model, or regularization can still leave the fitted probabilities miscalibrated, so calibration is something to measure on held-out data, and recalibration can be fitted on a separate calibration set and evaluated on another held-out set. The proper loss specifies the probability target; the calibration check assesses how closely the fitted reports achieve it.

Exercises

1. Squared error and far-away correct points. Fit a linear score to \(\pm 1\) labels by squared error, logistic loss, and hinge loss, with and without a few far-away correctly labeled points, and compare how each decision boundary moves.

You should get: the squared-error boundary shifts toward the far points, which were already classified correctly; the logistic and hinge boundaries show no change at the reported precision.

import numpy as np
from scipy.optimize import minimize

rng = np.random.default_rng(0)
x_base = np.concatenate([rng.normal(-1, 0.7, 40), rng.normal(1, 0.7, 40)])
y_base = np.concatenate([-np.ones(40), np.ones(40)])           # overlapping classes
x_far = np.concatenate([x_base, [8.0] * 5])
y_far = np.concatenate([y_base, np.ones(5)])                   # the 5 are correct

def boundary(loss, x, y):
    f = lambda p: np.mean(loss(y * (p[0] * x + p[1])))
    w, b = minimize(f, [1.0, 0.0], method="Nelder-Mead").x
    return -b / w

sq = lambda m: (1 - m) ** 2
lg = lambda m: np.logaddexp(0, -m)          # stable log(1 + exp(-m))
hi = lambda m: np.maximum(0, 1 - m)

print(f"{'loss':9s} {'without':>9} {'with 5 far':>11} {'shift':>8}")
for name, L in (("squared", sq), ("logistic", lg), ("hinge", hi)):
    b0, b1 = boundary(L, x_base, y_base), boundary(L, x_far, y_far)
    print(f"{name:9s} {b0:+9.4f} {b1:+11.4f} {b1 - b0:+8.4f}")
# loss        without  with 5 far    shift
# squared     +0.0843     +0.3505  +0.2662
# logistic    -0.0636     -0.0636  +0.0000
# hinge       -0.0758     -0.0758  +0.0000

2. Threshold or class weight. On an imbalanced problem, compare shifting the decision threshold against reweighting the classes, matched to flag the same number of observations. Report the recall, precision, Brier score, and mean predicted probability of each.

You should get: two routes to the same alert volume and, in this run, the same recall and precision, but probability losses that target different populations.

Recall is the fraction of actual positives flagged; precision is the fraction of flagged cases that are positive. The Brier score averages \( (r-y)^2 \) for probability reports \(r\) and labels \(y\in\{0,1\}\), so lower is better. Moving a threshold changes the first two metrics but leaves the underlying probabilities and their Brier score unchanged.

import numpy as np
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import precision_score, recall_score, brier_score_loss

rng = np.random.default_rng(0)
n = 20_000
X = rng.normal(size=(n, 4))
logit = -3.0 + 1.5 * X[:, 0] + 1.0 * X[:, 1]
y = (rng.random(n) < 1 / (1 + np.exp(-logit))).astype(int)
tr, te = slice(0, 15_000), slice(15_000, None)

plain = LogisticRegression().fit(X[tr], y[tr])
weighted = LogisticRegression(class_weight="balanced").fit(X[tr], y[tr])

p_plain = plain.predict_proba(X[te])[:, 1]
p_weight = weighted.predict_proba(X[te])[:, 1]

thr = np.quantile(p_plain, 1 - (p_weight > 0.5).mean())   # match the alert volume
pred_thr, pred_w = p_plain > thr, p_weight > 0.5
for name, pred, prob in (("threshold", pred_thr, p_plain),
                         ("class weight", pred_w, p_weight)):
    print(f"{name:13s} recall {recall_score(y[te], pred):.4f}"
          f"  precision {precision_score(y[te], pred):.4f}"
          f"  brier {brier_score_loss(y[te], prob):.4f}"
          f"  mean prob {prob.mean():.4f}")
print(f"actual positive rate {y[te].mean():.4f}")
print(f"alerts {pred_thr.sum()} vs {pred_w.sum()}, flagged differently {(pred_thr != pred_w).sum()}")
# threshold     recall 0.8093  precision 0.3230  brier 0.0733  mean prob 0.1184
# class weight  recall 0.8093  precision 0.3230  brier 0.1486  mean prob 0.3435
# actual positive rate 0.1164
# alerts 1458 vs 1458, flagged differently 2

3. Proper scoring rules. Show that the expected log loss and Brier score are minimized only by reporting the true probability, and that thresholded accuracy cannot single it out.

You should get: two scores with a unique minimum at the true probability, and one that is flat across every report on the same side of 0.5.

import numpy as np

true_p = 0.7
outcomes = np.array([1, 0])
weights = np.array([true_p, 1 - true_p])

for report in (0.5, 0.6, 0.7, 0.8, 1.0):
    r = np.clip(report, 1e-12, 1 - 1e-12)
    log_loss = -(weights * np.log(np.where(outcomes == 1, r, 1 - r))).sum()
    brier = (weights * (r - outcomes) ** 2).sum()
    accuracy = (weights * (outcomes == (r > 0.5))).sum()
    print(f"report {report:.1f}  log loss {log_loss:.4f}"
          f"  brier {brier:.4f}  accuracy {accuracy:.4f}")
# report 0.5  log loss 0.6931  brier 0.2500  accuracy 0.3000
# report 0.6  log loss 0.6325  brier 0.2200  accuracy 0.7000
# report 0.7  log loss 0.6109  brier 0.2100  accuracy 0.7000
# report 0.8  log loss 0.6390  brier 0.2200  accuracy 0.7000
# report 1.0  log loss 8.2893  brier 0.3000  accuracy 0.7000