Why Quantile Regression Is a Linear Program
Quantile regression estimates how a chosen conditional quantile of an outcome varies with regressors. Its computational advantage is that, for a fixed quantile level, minimizing the check-loss objective can be written as a linear program (LP). The regression coefficients need not be restricted in sign, and the residuals need not be known in advance to be positive or negative: the LP handles both possibilities with nonnegative variables.
Let \(y_i\) be the outcome for observation \(i\), let \(x_i\) be a \(k\)-vector of regressors (including a constant if the model has an intercept), and let \(X\) be the \(n\times k\) matrix whose \(i\)th row is \(x_i^\top\). For a quantile level \(\tau\) with \(0<\tau<1\), the sample quantile-regression estimate minimizes the sum of check losses of the residuals \(r_i=y_i-x_i^\top\beta\). This is the same check loss used in the earlier material on Check Loss and the Quantile Score Equation.
A positive residual means the fitted value is below the observed outcome and receives weight \(\tau\). A negative residual means the fitted value is above the outcome and receives weight \(1-\tau\). This asymmetry is what makes the fitted line target a quantile other than the median.
Turning Residuals Into LP Variables
For each observation, write its residual as the difference of two nonnegative variables: \(r_i=u_i-v_i\), where \(u_i\geq 0\) represents the positive part and \(v_i\geq 0\) represents the magnitude of the negative part. The residual equation becomes \(y_i=x_i^\top\beta+u_i-v_i\). The loss is then linear in \(u_i\) and \(v_i\).
Here \(\beta\) is unrestricted in sign; \(u\) and \(v\) are \(n\)-vectors, and \(y\) is the outcome vector. If a solver accepts unrestricted variables, enter \(\beta\) directly. If it only accepts nonnegative variables, write each coefficient as \(\beta_j=\beta_j^+-\beta_j^-\), with \(\beta_j^+,\beta_j^-\geq0\). The fitted coefficient is their difference.
Why does this LP have the same minimum as quantile regression? For any fixed residual \(r_i\), the constraint \(u_i-v_i=r_i\) permits \(u_i=\max(r_i,0)\) and \(v_i=\max(-r_i,0)\). These values give exactly \(\rho_\tau(r_i)\). If both \(u_i\) and \(v_i\) were positive, reducing both by their smaller value would preserve their difference and lower the objective, since both objective weights are positive. Thus an optimal decomposition does not need both to be positive, and its objective equals the check loss. Conversely, every regression coefficient vector supplies feasible residual variables by this positive-part and negative-part decomposition. The two optimization problems therefore have the same attainable minimum.
This equivalence concerns computation of the sample estimator. It requires finite observed outcomes and regressors, but it does not require normally distributed errors or a particular error variance. Interpreting the fit as a conditional quantile typically uses the model specification \(Q_\tau(Y\mid X=x)=x^\top\beta_\tau\) for the covariate values of interest. The LP itself does not establish that this specification is correct.
Worked Example: A Binary Regressor at Two Quantile Levels
Consider six hypothetical observations. Let \(x_i=0\) for the first group and \(x_i=1\) for the second, with a model containing an intercept and this binary regressor.
| Group | \(x_i\) | Observed \(y_i\) |
|---|---|---|
| 0 | 0 | 2 |
| 0 | 0 | 4 |
| 0 | 0 | 9 |
| 1 | 1 | 5 |
| 1 | 1 | 8 |
| 1 | 1 | 12 |
The design matrix has rows \((1,0)\) for group 0 and \((1,1)\) for group 1. The LP has \(\beta=(\beta_0,\beta_1)^\top\), six \(u_i\)'s, and six \(v_i\)'s, with constraints \(y_i=\beta_0+x_i\beta_1+u_i-v_i\). The fitted values in the two groups are \(\beta_0\) and \(\beta_0+\beta_1\). These can vary independently: any pair of group fitted values \(q_0,q_1\) corresponds to \(\beta_0=q_0\) and \(\beta_1=q_1-q_0\). The objective therefore separates into one check-loss minimization for each group.
At \(\tau=0.25\), group 0 is fit at 2. Its residuals are \(0,2,7\), so its loss is \(0.25(0+2+7)=2.25\). Group 1 is fit at 5, giving residuals \(0,3,7\) and loss \(0.25(0+3+7)=2.50\). The total objective is \(2.25+2.50=4.75\). The coefficients are \(\hat\beta_0=2\) and \(\hat\beta_1=5-2=3\). In LP variables, \(u=(0,2,7,0,3,7)^\top\), \(v=0\), and the objective is \(0.25\sum_i u_i=4.75\).
At \(\tau=0.75\), group 0 is fit at 9, with residuals \(-7,-5,0\); group 1 is fit at 12, with residuals \(-7,-4,0\). The total loss is \(0.25(7+5+7+4)=5.75\). Thus \(\hat\beta_0=9\) and \(\hat\beta_1=12-9=3\). Here \(u=0\), \(v=(7,5,0,7,4,0)^\top\), and the LP objective is \(0.25\sum_i v_i=5.75\).
Within each group, these fitted values minimize the check loss: with three ordered outcomes, moving a fitted value from below the lowest observation toward it lowers the loss, and moving it above the appropriate group quantile raises it. The result is a fitted conditional quantile of 2 versus 5 at \(\tau=0.25\), and 9 versus 12 at \(\tau=0.75\), for groups 0 and 1 respectively. The coefficient 3 is the difference between the two group quantiles at each level, not a causal effect.
Reading and Checking an LP Solution
A solver returns coefficient values, residual variables, and an objective value. Check the solution by confirming that each constraint \(y_i=x_i^\top\hat\beta+u_i-v_i\) holds to numerical tolerance, all \(u_i,v_i\) are nonnegative to numerical tolerance, and the reported objective equals \(\tau\sum u_i+(1-\tau)\sum v_i\). Small discrepancies can result from the solver's stopping tolerance.
The coefficient \(\hat\beta_{j,\tau}\) is a change in the fitted conditional \(\tau\)-quantile per one-unit change in regressor \(j\), holding the other regressors fixed, under the stated linear conditional-quantile specification. Its units are outcome units per regressor unit. It is not generally a change in the conditional mean. A full-column-rank design matrix rules out exact linear dependence among regressors, but by itself does not guarantee that the LP has a unique coefficient solution.
Worked Example: Intercept-Only Quantiles From an LP
Suppose four hypothetical outcomes are \(1,2,4,8\), and the only regressor is an intercept. The LP has one unrestricted coefficient \(b\) and constraints \(y_i=b+u_i-v_i\), with \(u_i,v_i\geq0\).
For \(\tau=0.4\), the solution is \(b=2\). The residuals are \(-1,0,2,6\), so \(u=(0,0,2,6)\), \(v=(1,0,0,0)\), and the objective is \(0.4(8)+0.6(1)=3.8\). To check optimality, for \(b\) between 1 and 2, the objective decreases as \(b\) increases: its slope is \(0.6(1)-0.4(3)=-0.6\). Between 2 and 4, its slope is \(0.6(2)-0.4(2)=0.4\), so it increases. Therefore \(b=2\) minimizes the objective.
For \(\tau=0.7\), the solution is \(b=4\). Its residuals are \(-3,-2,0,4\), and its objective is \(0.3(3+2)+0.7(4)=4.3\). Between 2 and 4 the slope is \(0.3(2)-0.7(2)=-0.8\); between 4 and 8 it is \(0.3(3)-0.7(1)=0.2\). The objective decreases up to 4 and increases beyond 4. In this intercept-only model, the fitted conditional quantile is constant across covariate values; the different fitted values at the two levels show how the target changes with \(\tau\).
A Design-Matrix Optimality Check
The LP constraints can sometimes provide a concise certificate that a proposed fit is optimal. The next example shows the idea without relying on a solver. Suppose \(X\) has three rows \((1,0),(1,1),(1,3)\), and \(y=(1,2,5)^\top\). For any coefficient vector, the residuals satisfy
because the same weighted combination of each column of \(X\) is zero.
Worked Example: Verifying Two Fits With the Residual Constraint
At \(\tau=0.25\), take residuals \((0.5,0,0)\). They satisfy \(2(0.5)-3(0)+0=1\), and correspond to \(\beta_0=0.5,\ \beta_1=1.5\): the fitted values are \(0.5,2,5\). The LP objective is \(0.25(0.5)=0.125\).
For any feasible residual vector, \(\rho_{0.25}(r)\geq0.125r\) for every real \(r\). For \(r\geq0\), this follows from \(\rho_{0.25}(r)=0.25r\); for \(r<0\), the loss is \(-0.75r\), which is also at least \(0.125r\). Therefore every feasible objective satisfies
so use the design-matrix relation with its coefficients: for every residual, \(\rho_{0.25}(r)\geq0.125(2r)\) when the contribution is positive through \(r_1\), and more generally the needed componentwise bound is \(\rho_{0.25}(r_i)\geq0.125c_i r_i\) for \(c=(2,-3,1)\). For \(c_i=2\), this holds for both signs; for \(c_i=-3\), the right side is nonpositive when \(r_i\geq0\), and when \(r_i<0\) it is \(0.375|r_i|\leq0.75|r_i|\); for \(c_i=1\), it also holds for both signs. Summing gives an objective of at least \(0.125(2r_1-3r_2+r_3)=0.125\). The proposed fit attains this lower bound and is therefore optimal.
At \(\tau=0.75\), take residuals \((0,-1/3,0)\). They satisfy \(2(0)-3(-1/3)+0=1\), and correspond to \(\beta_0=1,\ \beta_1=4/3\), with fitted values \(1,7/3,5\). The objective is \(0.25(1/3)=1/12\). For each residual, \(\rho_{0.75}(r_i)\geq(1/12)c_i r_i\) with \(c=(2,-3,1)\): for the negative second residual this is equality, and for each other sign the right side is no larger than the loss. Summing yields an objective at least \((1/12)(2r_1-3r_2+r_3)=1/12\). This fit is also optimal. The example illustrates how the design matrix restricts feasible residuals; the LP is the general computational method.
Practical Pitfalls
- Giving the residual variables the wrong weights. A positive residual is weighted by \(\tau\), while the magnitude of a negative residual is weighted by \(1-\tau\). Reversing them targets the opposite quantile level.
- Forcing regression coefficients to be nonnegative. Coefficients are generally unrestricted. If the solver requires nonnegative inputs, represent each coefficient as the difference of two nonnegative variables.
- Reading a coefficient as a mean effect or a causal effect. It describes a fitted conditional quantile under the model. A causal interpretation requires additional design and identification conditions; it does not follow from solving the LP.
- Assuming the solver's output is automatically correct. Check the residual constraints, nonnegativity, and objective value. A numerical optimizer can report a solution within tolerance, not necessarily exact arithmetic.
Check Your Understanding
Use the LP formulation and interpretation rules to answer the following questions.
- Write the objective weights on \(u_i\) and \(v_i\) for quantile level \(\tau=0.8\). Which variable represents an overprediction?
- If a solver requires nonnegative variables, how can an unrestricted slope coefficient be represented?
- In the binary-regressor example, what are the fitted group quantiles at \(\tau=0.75\), and what is the slope coefficient in outcome units per unit of the binary regressor?
- Why does the LP computation itself not require a normal-error assumption?
- What three checks can you apply to an LP solver's returned coefficients and residual variables?