Quantile Regression: Check Loss, KKT Conditions, and Scalable Optimization
Fix a quantile level $\tau \in (0, 1)$. The quantile regression estimator $\hat{\beta}$ minimizes the check-loss objective: $$\hat{\beta} = \arg\min_{\beta} \sum_{i=1}^{n} \rho_\tau(y_i - x_i^\top \beta), \quad \rho_\tau(u) = u(\tau - \mathbf{1}\{u < 0\})$$
(i) Show that $\hat{\beta}$ defined above is indeed the sample $\tau$-quantile regression estimator. Derive the first-order (subgradient) optimality conditions.
(ii) Reformulate the minimization as a linear program (LP). Define the auxiliary variables and write the LP explicitly.
(iii) For large-scale problems where the LP is too slow, describe coordinate descent with a Huberized approximation of the check loss. Explain why the Huberized approximation is needed and how to choose the smoothing parameter.
Hints
- The check loss $\rho_\tau(u) = \tau u^+ + (1-\tau)u^-$ is convex and piecewise linear. Its subgradient at $u=0$ is the interval $[-(1-\tau), \tau]$. Setting the subgradient of the total loss to zero gives you the optimality conditions.
- For the LP formulation, split each residual $r_i = y_i - x_i^\top \beta$ into its positive part $u_i = r_i^+$ and negative part $v_i = r_i^-$. The check loss becomes $\tau u_i + (1-\tau) v_i$ with equality constraint $u_i - v_i = r_i$ and $u_i, v_i \geq 0$.
- For coordinate descent, the check loss is non-smooth at zero, which causes slow convergence. Replace it with a Huberized version that adds a small quadratic region of width $h$ around zero. This makes the gradient Lipschitz continuous with constant $1/h$, enabling fast coordinate updates.
Worked Solution
How to Think About It: Quantile regression generalizes median regression ($\tau = 0.5$) to arbitrary quantile levels. The check loss $\rho_\tau(u)$ is an asymmetric absolute-value function: it penalizes positive residuals with weight $\tau$ and negative residuals with weight $(1-\tau)$. When $\tau = 0.5$, this reduces to $|u|/2$ (least absolute deviations). The function is convex but not differentiable at $u=0$, which makes subgradient calculus the right tool. The LP reformulation is the classic trick for turning piecewise-linear convex problems into linear programs by splitting residuals into positive and negative parts.
Key Insight: The check loss is piecewise linear, so its subgradient conditions have a natural interpretation: at the solution, the fraction of observations with negative residuals (fitted values above $y_i$) must equal exactly $\tau$, weighted by the design.
Part (i): Subgradient optimality conditions.
The check loss can be written: $$\rho_\tau(u) = \tau u^+ + (1-\tau)u^-, \quad u^+ = \max(u,0), \quad u^- = \max(-u,0)$$
So the objective is $Q(\beta) = \tau \sum_i r_i^+ + (1-\tau) \sum_i r_i^-$ where $r_i = y_i - x_i^\top\beta$.
The subgradient of $\rho_\tau$ with respect to $u$ is: $$\partial_u \rho_\tau(u) = \begin{cases} \tau & \text{if } u > 0 \\ -(1-\tau) & \text{if } u < 0 \\ [-(1-\tau),\, \tau] & \text{if } u = 0 \end{cases}$$
By the chain rule, $\partial Q / \partial \beta = -\sum_i x_i \partial_u \rho_\tau(r_i)$. Setting $0 \in \partial Q$, the KKT (subgradient) conditions are: $$\sum_{i: r_i > 0} \tau x_i - \sum_{i: r_i < 0} (1-\tau) x_i + \sum_{i: r_i = 0} s_i x_i = 0$$
where $s_i \in [-(1-\tau), \tau]$ for zero-residual observations. Rearranging: $$\sum_{i=1}^n x_i \psi_\tau(r_i) = 0, \quad \psi_\tau(u) = \tau - \mathbf{1}\{u < 0\}$$
(with the convention that $\psi_\tau(0)$ is any value in $[-(1-\tau), \tau]$). This says the weighted score function is zero at the solution -- a natural generalization of the OLS normal equations.
Part (ii): LP reformulation.
Introduce non-negative slack variables $u_i, v_i \geq 0$ for the positive and negative parts of each residual: $$r_i = y_i - x_i^\top\beta = u_i - v_i, \quad u_i, v_i \geq 0, \quad u_i v_i = 0$$
The constraint $u_i v_i = 0$ is automatically satisfied at optimality (since the objective penalizes both, it never pays to have both positive). The LP is: $$\min_{\beta,\, u,\, v \geq 0} \quad \tau \mathbf{1}^\top u + (1-\tau) \mathbf{1}^\top v$$ $$\text{subject to:} \quad X\beta + u - v = y, \quad u \geq 0, \quad v \geq 0$$
where $X \in \mathbb{R}^{n \times p}$, $y \in \mathbb{R}^n$. This is a standard LP with $p + 2n$ variables and $n$ equality constraints. It can be solved by the simplex method or interior-point methods.
Practical note: The LP formulation handles the non-differentiability cleanly, but for $n$ in the millions and $p$ in the thousands, the LP can be slow. This motivates Part (iii).
Part (iii): Coordinate descent with Huberized check loss.
Why Huberization? Coordinate descent requires computing one-dimensional updates of the form $\arg\min_t f(\beta + t e_j)$ where $e_j$ is the $j$-th basis vector. For the check loss, this is a piecewise-linear 1D problem solvable in closed form -- but it is non-smooth, and the zero-subgradient update can cycle or converge slowly near the solution (especially with correlated predictors).
The Huberized check loss replaces the kink at $u=0$ with a smooth quadratic transition: $$\rho_{\tau,h}(u) = \begin{cases} \tau u - h/2 & \text{if } u > h \\ u^2 / (2h) + (\tau - 1/2) u + h/8 & \text{if } |u| \leq h \\ -(1-\tau) u - h/2 & \text{if } u < -h \end{cases}$$
for a smoothing parameter $h > 0$. As $h \to 0$, $\rho_{\tau,h} \to \rho_\tau$. The key properties: $\rho_{\tau,h}$ is convex, continuously differentiable, and its gradient is Lipschitz with constant $1/h$.
Coordinate descent update: With a smooth objective, the 1D update for coordinate $j$ is: $$\hat{t} = \arg\min_t \sum_i \rho_{\tau,h}(r_i - t x_{ij})$$
This has a closed-form solution (or near-closed-form via a simple Newton step) because the smoothed loss is locally quadratic. Cycle through all coordinates until convergence.
Choosing $h$: A standard data-driven choice is $h \propto n^{-1/3}$ (the minimax-optimal smoothing bandwidth for quantile problems). In practice, start with $h = \hat{\sigma} / \sqrt{n}$ where $\hat{\sigma}$ is the residual standard deviation from an initial OLS fit, then decrease $h$ across iterations (continuation method) to tighten the approximation.
Answer: - (i) $\hat{\beta}$ satisfies $\sum_i x_i(\tau - \mathbf{1}\{y_i < x_i^\top \hat{\beta}\}) = 0$, the quantile regression score equation. - (ii) LP: minimize $\tau \mathbf{1}^\top u + (1-\tau)\mathbf{1}^\top v$ subject to $X\beta + u - v = y$, $u, v \geq 0$. - (iii) Smooth the check loss near zero using a quadratic Huber kernel with bandwidth $h \to 0$; run coordinate descent on the smoothed objective; decrease $h$ via a continuation schedule for convergence to the exact quantile solution.
Intuition
Quantile regression is the natural generalization of median regression: instead of penalizing squared deviations (which targets the mean), you penalize asymmetric absolute deviations, and the degree of asymmetry ($\tau$ vs $1-\tau$) controls which quantile you are targeting. This gives you a way to model the entire conditional distribution of $y$ given $x$, not just its conditional mean -- which is enormously useful in risk management (value-at-risk is a quantile), economics (distributional treatment effects), and forecasting (prediction intervals).
The LP connection is fundamental: because the loss is piecewise linear, the solution always lies at a vertex of the feasible polytope (a basic feasible solution), which means at most $p$ of the $n$ residuals are exactly zero at the optimum. This is the quantile regression analog of OLS having a unique solution: generically, $p$ observations are 'pinned' to the fitted line. The Huberized approximation for large-scale problems trades a tiny bias for dramatically faster convergence -- a classic bias-variance tradeoff in optimization.