Quantile Regression

The Pinball Loss and the Population Quantile

Quantile regression rests on a single asymmetric loss function, the check or pinball loss. Where squared loss penalises positive and negative residuals symmetrically and recovers the conditional mean at its minimum, the pinball loss penalises them asymmetrically and recovers the conditional quantile. The next two results make this statement precise.

Numerical verification: sample $Y \sim \mathcal{N}(0, 1)$ and minimize the empirical pinball loss over a fine grid of $q$ for several $\tau$ — the empirical argmin should track $\Phi^{-1}(\tau)$ as $n$ grows.

import numpy as np
from scipy.stats import norm

np.random.seed(0)
y_sample = norm.rvs(size=20_000, random_state=1)
q_grid = np.linspace(-3, 3, 601)

print('Empirical vs theoretical quantile at n = 20,000:')
print(f'{"tau":>6} | {"Phi^-1(tau)":>12} | {"argmin emp loss":>16} | {"diff":>8}')
for tau in [0.10, 0.25, 0.50, 0.75, 0.90]:
    L_grid = np.array([
        np.mean(np.where(y_sample >= q,
                         tau * (y_sample - q),
                         (tau - 1) * (y_sample - q)))
        for q in q_grid
    ])
    q_hat = q_grid[L_grid.argmin()]
    q_true = norm.ppf(tau)
    print(f'{tau:>6.2f} | {q_true:>12.4f} | {q_hat:>16.4f} | {q_hat - q_true:>+8.4f}')

τ0.50Y ∼Normal(0, 1)Exponential(1)Uniform[0, 1]

Linear Quantile Regression

Theorem 1 says the population $\tau$-quantile of $Y$ is recovered by minimising the expected pinball loss. The Koenker–Bassett 1978 quantile regression estimator is the empirical-version twin of that statement: minimise the sample pinball loss over a parametric family of conditional-quantile candidates.

The workhorse fit-and-predict helper used throughout this topic — a thin wrapper around scikit-learn’s QuantileRegressor (which calls HiGHS) on degree-3 polynomial features:

from sklearn.linear_model import QuantileRegressor
from sklearn.preprocessing import PolynomialFeatures

def fit_predict_quantile(x_train, y_train, x_eval, tau, alpha_l2=0.01, degree=3):
    """Quantile regression at level τ on degree-`degree` polynomial features.

    alpha_l2: small L2 penalty for numerical stability (kept tiny so KB78
              asymptotics are not distorted).
    """
    poly = PolynomialFeatures(degree=degree, include_bias=False)
    Phi_train = poly.fit_transform(np.asarray(x_train).reshape(-1, 1))
    Phi_eval = poly.transform(np.asarray(x_eval).reshape(-1, 1))
    qr = QuantileRegressor(quantile=tau, alpha=alpha_l2, solver='highs')
    qr.fit(Phi_train, y_train)
    return qr.predict(Phi_eval)

The KKT condition $L'(q) = F(q) - \tau = 0$ has a clean empirical analog: at the QR optimum, fraction $\tau$ of training residuals are negative.

# At each τ, fraction of training residuals < 0 should ≈ τ.
for tau in [0.10, 0.50, 0.90]:
    fit_at_train = fit_predict_quantile(x_demo, y_demo, x_demo, tau)
    print(f'τ = {tau}:  P(Y < q̂(X)) = {(y_demo < fit_at_train).mean():.3f}')

τ0.50n200degree3multi-τ (K = 11)rearrange (CFV-G)Resample

The LP Reformulation

The pinball loss is piecewise-linear, so the QR optimisation is a linear program. This is more than a curiosity: it explains why sklearn’s QuantileRegressor uses HiGHS (a state-of-the-art LP solver), why QR is exact rather than iterative, and why the “interpolation” structure of the fit is combinatorial — exactly $p$ data points sit on the fitted line in the non-degenerate case.

The trick is to split each residual $r_i = Y_i - X_i^\top \beta$ into its positive and negative parts:

$$r_i \;=\; u_i^+ - u_i^-, \qquad u_i^+ = \max(r_i, 0), \quad u_i^- = \max(-r_i, 0).$$

Then $\rho_\tau(r_i) = \tau\, u_i^+ + (1 - \tau)\, u_i^-$, which is linear in the slack variables. Stacking the $n$ constraints:

Solving this LP directly with scipy.optimize.linprog makes Definition 3 concrete; the result agrees with QuantileRegressor to machine precision since both call HiGHS.

from scipy.optimize import linprog

def solve_qr_linprog(X, y, tau):
    """Solve QR via the LP reformulation using scipy.optimize.linprog (HiGHS).

    X: (n, p) design matrix (include constant column for intercept).
    y: (n,) responses.
    tau: quantile level in (0, 1).
    Returns: (beta, residuals, slack_pos, slack_neg).
    """
    n, p = X.shape
    # z = [β (p), u⁺ (n), u⁻ (n)],  total p + 2n variables.
    c = np.concatenate([np.zeros(p), tau * np.ones(n), (1 - tau) * np.ones(n)])
    # Equality constraint: X β + u⁺ - u⁻ = y    (slacks split residual sign)
    A_eq = np.hstack([X, np.eye(n), -np.eye(n)])
    b_eq = y
    # β unbounded; slacks ≥ 0
    bounds = [(None, None)] * p + [(0, None)] * (2 * n)
    res = linprog(c, A_eq=A_eq, b_eq=b_eq, bounds=bounds, method='highs')
    z = res.x
    beta = z[:p]
    return beta, y - X @ beta, z[p:p + n], z[p + n:]

Equivariance Under Monotone Transformations

Conditional quantiles have a property the conditional mean lacks: they commute with monotone transformations of $Y$. If $Y$ is income and we want the median of $\log(Y)$, the answer is $\log$ of the median of $Y$ — no Jensen correction required. This makes quantile regression unusually robust to the scale on which the response is modelled.

Asymptotic Theory: Koenker–Bassett 1978, Knight 1998

We just saw that QR returns the conditional $\tau$-quantile in expectation. The next question is: at what rate does the estimator converge, and what is the limit distribution? The answer, established by Koenker–Bassett 1978 with the clean modern proof due to Knight 1998, is the same root-$n$ / Gaussian shape we know from OLS, with two structural differences. The “noise variance” is replaced by $\tau(1 - \tau)$ — the Bernoulli variance of the indicator $\mathbb{1}\{Y < \hat q(X)\}$ — and the design Gram matrix is sandwich-replaced by a density-weighted Gram matrix evaluated at the conditional quantile.

The bootstrap is the most direct way to see Theorem 3 — resample $(X, Y)$ pairs $B$ times, refit QR, and read off the empirical distribution of any coefficient.

def fit_qr_coefs(x_train, y_train, tau, degree=3):
    """Return all polynomial coefficients (intercept + slopes) of QR at level τ."""
    poly = PolynomialFeatures(degree=degree, include_bias=False)
    Phi = poly.fit_transform(np.asarray(x_train).reshape(-1, 1))
    qr = QuantileRegressor(quantile=tau, alpha=0.0, solver='highs')
    qr.fit(Phi, y_train)
    return np.concatenate([[qr.intercept_], qr.coef_])  # [a₀, a₁, a₂, a₃]


def bootstrap_qr_coefs(x, y, tau, B, rng, degree=3):
    """B (X, Y)-bootstrap resamples; returns (B, p+1) coefficient matrix."""
    n = len(x)
    out = np.empty((B, degree + 1))
    for b in range(B):
        idx = rng.integers(0, n, n)
        out[b] = fit_qr_coefs(x[idx], y[idx], tau, degree=degree)
    return out

At fixed $\tau$, the empirical std of the slope coefficient $\hat\beta_1(\tau)$ across bootstrap replicates should scale as $1/\sqrt{n}$ (so $\hat\sigma \cdot \sqrt{n}$ stabilises). Across $\tau$, the constant inflates as $\tau \to 0$ or $\tau \to 1$ — Remark 8.

τ0.50n200Run B = 200 bootstraps(sliders auto-run B = 50 for fast feedback)

Multi-Quantile Estimation, Crossing, and Rearrangement

A natural use of QR is to fit multiple quantile levels simultaneously and read the resulting bands as a non-parametric estimate of the conditional distribution of $Y$ given $X$. Three or five evenly-spaced quantiles already give a usable picture; for high-resolution density estimation, fit at a fine grid in $(0, 1)$.

But there is a problem. The population conditional quantiles are weakly increasing in $\tau$ by definition: $\xi_{\tau_1}(Y \mid X = x) \leq \xi_{\tau_2}(Y \mid X = x)$ whenever $\tau_1 \leq \tau_2$. The marginal QR estimates do not enforce this. Each $\tau$-fit is a separate optimisation, and at finite $n$ in finite-dimensional function classes, the resulting curves can — and frequently do — cross: the $\tau = 0.7$ fit can dip below the $\tau = 0.6$ fit at some $x$.

Two helpers — fit a $K$-quantile bundle independently and then sort along the $\tau$-axis at every evaluation point:

def fit_multiple_quantiles(x_train, y_train, x_eval, taus, alpha_l2=0.01, degree=3):
    """Fit QR at each τ in `taus` on degree-`degree` polynomial features.

    Returns a (K, |x_eval|) array Q where Q[k, j] = q̂_{taus[k]}(x_eval[j]).
    """
    poly = PolynomialFeatures(degree=degree, include_bias=False)
    Phi_train = poly.fit_transform(np.asarray(x_train).reshape(-1, 1))
    Phi_eval = poly.transform(np.asarray(x_eval).reshape(-1, 1))
    out = np.empty((len(taus), len(x_eval)))
    for k, tau in enumerate(taus):
        qr = QuantileRegressor(quantile=tau, alpha=alpha_l2, solver='highs')
        qr.fit(Phi_train, y_train)
        out[k] = qr.predict(Phi_eval)
    return out


def rearrange_quantile_estimates(Q):
    """CFV-G 2010 rearrangement: sort along the τ-axis at each evaluation point."""
    return np.sort(Q, axis=0)

τ0.50n200degree3multi-τ (K = 11)rearrange (CFV-G)Resample

Penalized Quantile Regression

Two settings push QR away from its plain Koenker–Bassett 1978 form. First, when $p > n$ or $p \sim n$ the unregularised QR LP is under-determined and badly behaved. Second, even at moderate $p$ with $n$ large enough for asymptotics, an $L^1$ penalty drives variable selection that simple shrinkage cannot match. Both settings use the same trick: add a penalty term to the pinball-loss objective.

The $L^2$ penalty (already familiar from sklearn’s alpha parameter):

$$\min_\beta \; \sum_i \rho_\tau(Y_i - X_i^\top \beta) \;+\; \lambda\, \lVert \beta \rVert_2^2.$$

This has a closed-form Hessian (modulo the non-smoothness of $\rho_\tau$) and can be solved either by augmenting the LP with linear-quadratic structure or via the smoothed-check-loss accelerated-gradient solver we use in the in-browser viz components.

The $L^1$ penalty (the focus of this section, “QR-lasso” or “L1-QR”):

$$\min_\beta \; \sum_i \rho_\tau(Y_i - X_i^\top \beta) \;+\; \lambda\, \lVert \beta \rVert_1.$$

This stays a linear program, since $\lVert \beta \rVert_1 = \mathbf{1}^\top \beta^+ + \mathbf{1}^\top \beta^-$ splits cleanly with $\beta = \beta^+ - \beta^-$. Augment the QR LP with these auxiliary variables and solve for the full path of $\lambda$ values.

The L1-QR LP splits $\beta = \beta^+ - \beta^-$ with $\beta^+, \beta^- \geq 0$, then adds the slack pair $u^+, u^-$ from §3. Four nonnegative variable groups give a standard-form LP that scipy / HiGHS solves directly:

def solve_lasso_qr_linprog(X, y, tau, lam, no_penalty_mask=None):
    """L1-penalised QR via the standard LP reformulation.

    Variables z = [β⁺ (p), β⁻ (p), u⁺ (n), u⁻ (n)].
    Objective:   λ·(1ᵀβ⁺ + 1ᵀβ⁻) + τ·1ᵀu⁺ + (1−τ)·1ᵀu⁻,
    with the λ coefficients zeroed out at indices in no_penalty_mask
    (typically the intercept column).
    Constraints: X(β⁺ − β⁻) + u⁺ − u⁻ = y, all variables ≥ 0.
    """
    n, p = X.shape
    if no_penalty_mask is None:
        no_penalty_mask = np.zeros(p, dtype=bool)
    pen_costs = np.where(no_penalty_mask, 0.0, lam)
    c = np.concatenate([
        pen_costs,                  # cost on β⁺
        pen_costs,                  # cost on β⁻
        tau * np.ones(n),           # cost on u⁺
        (1 - tau) * np.ones(n),     # cost on u⁻
    ])
    A_eq = np.hstack([X, -X, np.eye(n), -np.eye(n)])
    bounds = [(0, None)] * (2 * p + 2 * n)
    res = linprog(c, A_eq=A_eq, b_eq=y, bounds=bounds, method='highs')
    z = res.x
    return z[:p] - z[p:2 * p]      # β = β⁺ − β⁻

Quantile Regression as the Base Learner of Conformalized QR

We close the loop with the use that motivated this topic in the T4 track: quantile regression as the base learner inside Conformalized Quantile Regression (Topic 1, §6). CQR uses two QR fits — at levels $\alpha/2$ and $1 - \alpha/2$ — to produce a heteroscedastic prediction band, and then applies conformal calibration on top to guarantee finite-sample marginal coverage at level $1 - \alpha$. Conformal’s coverage theorem holds for any base learner (Conformal Prediction, Theorem 1); the role of QR specifically is to give the band the right shape.

This is the division of labour:

• QR base fits provide the band SHAPE — wide where conditional variance is large, narrow where it is small. A symmetric residual-conformal interval does not have this property; it produces constant-width bands.
• Conformal calibration provides the band WIDTH — adjust the QR fits by an additive constant so that the empirical miscoverage rate on a held-out calibration set matches the target $\alpha$.

Conformal Prediction’s Theorem 1 (split-conformal validity) gives marginal coverage $1 - \alpha$ regardless of how badly QR is misspecified — even constant fits. What QR contributes is conditional approximate validity: when the linear class spans the true conditional quantile function, the conformal correction is a small constant and the resulting band tracks the true conditional coverage rate uniformly. When QR is misspecified, the conformal correction absorbs the misspecification globally; the resulting band still has marginal coverage $1 - \alpha$ but conditional coverage will be uneven. The full treatment of prediction-interval procedures — comparing CQR with locally adaptive variants (CQR-r, CQR-m), conditional-coverage methods, and base learners beyond linear QR — lives in Prediction Intervals (coming soon).

Recall (callback to Conformal Prediction §6, Definition 4) the CQR prediction set at level $1 - \alpha$ based on QR fits $\hat q_{\alpha/2}$ and $\hat q_{1 - \alpha/2}$:

$$\hat C_\alpha(x) \;=\; \bigl[\,\hat q_{\alpha/2}(x) - \hat Q_\alpha,\quad \hat q_{1 - \alpha/2}(x) + \hat Q_\alpha\,\bigr],$$

where $\hat Q_\alpha$ is the $(1 - \alpha)$ empirical quantile (with the standard $1/n_{\text{cal}}$ correction) of the calibration nonconformity scores

$$E_i \;=\; \max\bigl\{\,\hat q_{\alpha/2}(X_i) - Y_i,\ \ Y_i - \hat q_{1 - \alpha/2}(X_i)\,\bigr\}.$$

The two QR fits define the band shape; $\hat Q_\alpha$ is a single scalar that inflates or deflates the band uniformly to hit the target coverage.

Connections and Further Reading

Within formalML, this topic feeds into:

• Conformal Prediction (T4 #1). The two-sided QR fits at $\alpha/2$ and $1 - \alpha/2$ are the base learner inside Conformalized Quantile Regression. Theorem 1 of conformal-prediction guarantees marginal coverage of the CQR band regardless of QR’s specification; the role of QR is to give the band the right shape — heteroscedastic where the data are heteroscedastic. See §8 above.

• Prediction Intervals (coming soon). A unifying treatment of frequentist prediction intervals — fixed-width residual-based, conformal, and quantile-regression-based — and their coverage guarantees under various assumptions (i.i.d., exchangeable, group-conditional). Quantile regression is one of three “spokes” feeding into that umbrella topic.

• Extreme Value Theory (coming soon). Theorem 3’s variance formula tells us QR’s asymptotic variance inflates as $\tau \to 0$ or $\tau \to 1$, because the density-weighted matrix $D(\tau)$ shrinks in the tail. Beyond moderate quantile levels — typically $\tau$ outside roughly $[0.05, 0.95]$ for a Gaussian-tailed $Y$ — direct QR is not the right framework. EVT replaces it with the generalised Pareto / generalised extreme-value families and a peaks-over-threshold estimator that targets the tail directly.

Cross-site prerequisites:

• formalStatistics: Linear Regression — the OLS analog. The KB78 estimator is the pinball-loss replacement of squared loss; the LP reformulation is the QR replacement of the normal equations; Theorem 3’s asymptotic Gaussian is the QR replacement of OLS’s $\hat\beta_{OLS} \sim \mathcal{N}(\beta, \sigma^2 (X^\top X)^{-1})$.

• formalStatistics: Order Statistics and Quantiles — the no-covariate special case. With a single intercept feature, KB78 reduces to the empirical $\tau$-quantile of $\{Y_i\}$; Theorem 3 reduces to the classical Bahadur–Ghosh asymptotics.

• formalStatistics: Empirical Processes — the toolkit behind Theorem 3’s proof sketch (Knight’s identity, the empirical-process limit of the rescaled objective, the argmin lemma) and behind Theorem 4’s restricted-eigenvalue / oracle-inequality argument.

Internal prerequisites:

• Convex Analysis — the pinball loss is convex but non-smooth; the LP reformulation exploits its piecewise-linear structure.

• Gradient Descent — the smoothed-check-loss accelerated-gradient solver is what powers the in-browser visualisation widgets, where running an LP solver in the user’s browser is impractical.