#### Problem 1: Finding the median as an optimization problem

NLP model for finding the median |
---|

\[\min\>\sum_i \left|\color{darkblue}y_i - \color{darkred}m\right|\] |

- You may want to think a little bit about how this model indeed finds the median (or rather a median with 50% of the data below and 50% above).
- This model looks deceptively simple. Actually, this is a nasty non-differentiable NLP problem. Many general-purpose NLP solvers may have a really difficult time with it. Most of them will not be able to establish (local) optimality. In the box below are some of the frightening messages you may encounter. Standard NLP solvers expect smoothness: functions and gradients should be continuous. Ignoring this may be a bad idea.
- This is a great example of why we always need to reformulate absolute values. Absolute values are wolves in sheep's clothing and they will eat you alive.
- Note that this model works also if the \(\color{darkblue}y_i\) are endogenous instead of just exogenous data.
- Why not just sort the data? This is related to the previous point. Sorting data is easy, sorting decision variables is a different thing altogether and much, much more complicated.

---- 11 PARAMETERydatai1 0.172, i2 0.843, i3 0.550, i4 0.301, i5 0.292 **** SOLVER STATUS 4 Terminated By Solver **** MODEL STATUS 7 Feasible Solution **** OBJECTIVE VALUE 0.9297 LOWER LEVEL UPPER MARGINAL ---- VAR m -INF 0.3011 +INF 1.0000 NOPT ---- VAR z -INF 0.9297 +INF . m median z objective Scary messages from different solvers: Conopt: ** Feasible solution. Convergence too slow. The change in objective has been less than 3.0000E-12 for 20 consecutive iterations MINOS: EXIT - The current point cannot be improved. SNOPT: EXIT - Current point cannot be improved. IPOPT: EXIT: Restoration Failed! Final point is feasible: scaled constraint violation (0) is below tol (1e-08) and unscaled constraint violation (0) is below constr_viol_tol (0.0001). Knitro: EXIT: Primal feasible solution estimate cannot be improved; desired accuracy in dual feasibility could not be achieved.

**reformulate absolute values**in a linear fashion, so we end up with an LP. There are (at least) two ways to do this:

LP model 1 for finding the median |
---|

\[\begin{align}\min&\sum_i\color{darkred}e_i \\ & \color{darkred}e_i\ge\color{darkblue}y_i -\color{darkred}m&& \forall i \\ & \color{darkred}e_i \ge -(\color{darkblue}y_i -\color{darkred}m) && \forall i \\ &\color{darkred}e_i\ge 0\end{align}\] |

**variable splitting**:

LP model 2 for finding the median |
---|

\[\begin{align}\min&\sum_i \left(\color{darkred}e^+_i+\color{darkred}e^-_i\right) \\ & \color{darkred}e^+_i-\color{darkred}e^-_i =\color{darkblue}y_i -\color{darkred}m && \forall i \\ &\color{darkred}e^-_i,\color{darkred}e^+_i\ge 0\end{align}\] |

---- 69 VARIABLEe2.Labsolute value, formulation 2+ - i1 0.129 i2 0.542 i3 0.249 i5 0.009

For every error \(i\) we see that only one of \(\color{darkred}e^+_i,\color{darkred}e^-_i\) will be nonzero. This is what we would expect. Note that there is one error that is zero.

**LP Model 2**is really the building block for Quantile Regression.

#### Problem 2: Finding a quantile as an optimization problem

NLP model for finding the \(\tau\)-th quantile |
---|

\[\min\>\sum_{i|\color{darkblue}y_i\ge\color{darkred}q} \color{darkblue}\tau\left|\color{darkblue}y_i - \color{darkred}q \right| +\sum_{i|\color{darkblue}y_i\lt\color{darkred}q} (1-\color{darkblue}\tau)\left|\color{darkblue}y_i-\color{darkred}q \right| \] |

**variable splitting**.

LP model for finding the \(\tau\)-th quantile |
---|

\[\begin{align}\min&\sum_i \left( \color{darkblue}\tau \cdot \color{darkred}e^+_i +(1-\color{darkblue}\tau) \cdot \color{darkred}e^-_i \right) \\ &\color{darkred}e^+_i -\color{darkred}e^-_i = \color{darkblue}y_i - \color{darkred}q && \forall i \\ &\color{darkred}e^-_i,\color{darkred}e^+_i\ge 0\end{align} \] |

---- 15 PARAMETERydatai1 25.457, i2 85.894, i3 59.534, i4 37.102, i5 36.299, i6 30.165, i7 41.485, i8 87.064 i9 16.040, i10 55.019, i11 99.831, i12 62.086, i13 99.202, i14 78.603, i15 21.762, i16 67.575 i17 24.357, i18 32.507, i19 70.204, i20 49.182, i21 42.373, i22 41.630, i23 21.834, i24 23.509 i25 63.020 ---- 15 SETtquantile levels0 , 0.25, 0.5 , 0.75, 1 ---- 52 PARAMETERquantilesSolution0 16.040, 0.25 30.165, 0.5 42.373, 0.75 67.575, 1 99.831

*type=1*:

#### Problem 3: Quantile regression

NLP model for quantile regression |
---|

\[\begin{align}\min&\sum_{i|\color{darkred}e_i\ge 0} \color{darkblue}\tau|\color{darkred}e_i| +\sum_{i|\color{darkred}e_i\lt 0} (1-\color{darkblue}\tau) |\color{darkred}e_i| \\ & \color{darkred}e = \color{darkblue}y-\color{darkblue}X\cdot\color{darkred}\beta\end{align}\] |

LP model for quantile regression |
---|

\[\begin{align}\min&\sum_i \left( \color{darkblue}\tau \cdot\color{darkred}e^+_i + (1-\color{darkblue}\tau) \cdot \color{darkred}e^-_i \right)\\ & \color{darkred}e^+_i - \color{darkred}e^-_i = \color{darkblue}y_i-\sum_j \color{darkblue}X_{i,j}\cdot\color{darkred}\beta_j && \forall i\\ &\color{darkred}e^+_i,\color{darkred}e^-_i\ge 0\end{align}\] |

---- 96 PARAMETERestimatessuppins age white female totchr intercept 0.25 453.444 16.083 338.083 16.056 782.472 -1412.889 0.5 687.222 35.111 632.889 -260.556 1332.833 -2252.556 0.75 708.409 87.364 801.682 -554.591 2855.318 -4512.045

**totchr**. Indeed the coefficient for this variable changes a lot for different quantile levels. The R code from [3] produces:

>results <- rq(Y ~ X, data=mydata, tau=c(0.25, 0.5, 0.75))>summary(results)Call: rq(formula = Y ~ X, tau = c(0.25, 0.5, 0.75), data = mydata) tau: [1] 0.25 Coefficients: Value Std. Error t value Pr(>|t|) (Intercept) -1412.88889 433.20179 -3.26150 0.00112 Xsuppins 453.44444 75.05348 6.04162 0.00000 Xtotchr 782.47222 37.55769 20.83388 0.00000 Xage 16.08333 6.19162 2.59760 0.00943 Xfemale 16.05556 72.20278 0.22237 0.82404 Xwhite 338.08333 71.51522 4.72743 0.00000 Call: rq(formula = Y ~ X, tau = c(0.25, 0.5, 0.75), data = mydata) tau: [1] 0.5 Coefficients: Value Std. Error t value Pr(>|t|) (Intercept) -2252.55556 846.23023 -2.66187 0.00781 Xsuppins 687.22222 137.29264 5.00553 0.00000 Xtotchr 1332.83333 74.77913 17.82360 0.00000 Xage 35.11111 11.29450 3.10869 0.00190 Xfemale -260.55556 150.46285 -1.73169 0.08343 Xwhite 632.88889 243.05734 2.60387 0.00926 Call: rq(formula = Y ~ X, tau = c(0.25, 0.5, 0.75), data = mydata) tau: [1] 0.75 Coefficients: Value Std. Error t value Pr(>|t|) (Intercept) -4512.04545 2350.56284 -1.91956 0.05501 Xsuppins 708.40909 375.76929 1.88522 0.05950 Xtotchr 2855.31818 196.12587 14.55860 0.00000 Xage 87.36364 30.98410 2.81963 0.00484 Xfemale -554.59091 378.71501 -1.46440 0.14319 Xwhite 801.68182 370.96108 2.16109 0.03077

**dual**version of this model [1]. It looks like:

Dual LP model for quantile regression |
---|

\[\begin{align}\max&\sum_i \color{darkblue}y_i \cdot \color{darkred}d_i\\ & \sum_i \color{darkblue}X_{i,j} \cdot \color{darkred}d_i = 0 \perp \color{darkred} \beta_j && \forall j\\ & \color{darkred}d_i \in [\color{darkblue}\tau-1,\color{darkblue}\tau] \end{align}\] |

---- 131 PARAMETERestimates2from dualsuppins age white female totchr intercept 0.25 453.444 16.083 338.083 16.056 782.472 -1412.889 0.5 687.222 35.111 527.556 -260.556 1332.833 -2147.222 0.75 708.409 87.364 801.682 -554.591 2855.318 -4512.045

#### Conclusion

#### References

- Koenker, R., and Bassett, G. W. (1978). “Regression Quantiles.”
*Econometrica*46:33–50. *Quantile Regression*, https://en.wikipedia.org/wiki/Quantile_regression- Ani Katchova, Quantile Regression, https://sites.google.com/site/econometricsacademy/econometrics-models/quantile-regression, this is a very good introduction with examples
- Linear Programming and L1 Regression, https://yetanothermathprogrammingconsultant.blogspot.com/2017/11/lp-and-lad-regression.html
*Median*, https://en.wikipedia.org/wiki/Median.*Quantile*, https://en.wikipedia.org/wiki/Quantile. Note: in case you thought you knew what quantiles are, R supports 9 different types of them.

#### Appendix A: GAMS code for non-smooth NLP model to calculate the median

*abs()*inside an NLP model, we need to declare the model as a

**DNLP**. This is a warning that trouble may be ahead.

$ontext This
model will have difficulties due to non-differentiabilityof the
abs() function.References:http://yetanothermathprogrammingconsultant.blogspot.com/2021/06/median-quantiles-and-quantile.htmlerwin@amsterdamoptimization.com$offtext set i /i1*i5/;parameter y(i) 'data';y(i) = uniform(0,1); display y;variablem 'median'z 'objective'; equation sumabs;sumabs.. z =e= sum(i, abs(y(i)-m));* initial valuem.l = 0.5; model median /sumabs/;option dnlp=conopt;solve median minimizing z
using dnlp;display m.l; |

#### Appendix B: GAMS code for finding the median as an LP

$ontext Use two
formulations for the absolute valuesReferences:http://yetanothermathprogrammingconsultant.blogspot.com/2021/06/median-quantiles-and-quantile.htmlerwin@amsterdamoptimization.com$offtext seti /i1*i5/ ; parameter y(i) 'data';y(i) = uniform(0,1); display y;*-------------------------------------------------------------* formulation 1: bounding*-------------------------------------------------------------variablesm 'median'z 'objective'; positive variablese(i) 'absolute value, formulation
1'; equationsobj1 'objective, formulation 1'bound1(i) bound2(i) ; obj1.. z =e= sum(i, e(i));bound1(i).. e(i) =g= m-y(i); bound2(i).. e(i) =g= -(m-y(i)); model medianLP1 /obj1,bound1,bound2/;solve medianLP1 minimizing z
using lp;display m.l;*-------------------------------------------------------------* formulation 2: variable splitting*-------------------------------------------------------------setpm /'+','-'/ ; positive variablese2(i,pm) 'absolute value, formulation
2'; equationsobj2 'objective, formulation 2'split(i) 'variable splitting'; obj2.. z =e= sum((i,pm), e2(i,pm));split(i).. e2(i,'+')-e2(i,'-') =e= m-y(i); model medianLP2 /obj2,split/;solve medianLP2 minimizing z
using lp;display m.l;display e2.l; |

#### Appendix C: GAMS code for finding quantiles.

$ontext References:http://yetanothermathprogrammingconsultant.blogspot.com/2021/06/median-quantiles-and-quantile.htmlerwin@amsterdamoptimization.com$offtext seti 'observations' /i1*i25/t 'quantile levels' /'0','0.25','0.5','0.75','1'/; parameter y(i) 'data';y(i) = uniform(10,100); display y,t;*-------------------------------------------------------------* variable splitting LP Model*-------------------------------------------------------------setpm /'+','-'/ ; scalar tau 'quantile level';positive variables e(i,pm) 'absolute value';variablez 'objective variable'q 'quantile'; equationsobj 'objective'split(i) 'variable splitting'; obj.. z =e= sum(i, tau*e(i,'+') + (1-tau)*e(i,'-'));split(i).. e(i,'+')-e(i,'-') =e= y(i)-q; model quantileLP /obj,split/;*-------------------------------------------------------------* solve loop*-------------------------------------------------------------parameter quantiles(t)
'Solution';loop(t,tau = t.val; solve quantileLP minimizing z using lp;quantiles(t) = q.l; ); display quantiles; |

#### Appendix D: GAMS code for quantile regression

---- 34 PARAMETER data all data totexp ltotexp suppins age white female totchr 93193020 3.000 1.099 1.000 69.000 1.000 72072017 6.000 1.792 1.000 65.000 1.000 1.000 25296013 9.000 2.197 85.000 1.000 1.000 23628011 14.000 2.639 76.000 1.000 1.000 95041014 18.000 2.890 71.000 1.000 1.000 1.000 25090018 20.000 2.996 81.000 1.000 1.000 96478017 24.000 3.178 74.000 1.000 96815023 25.000 3.219 1.000 83.000 1.000 1.000 93118011 29.000 3.367 80.000 1.000 1.000 24669013 30.000 3.401 73.000 1.000 1.000 3.000 93317018 30.000 3.401 1.000 78.000 1.000 97173010 31.000 3.434 1.000 69.000 1.000 . . .

$ontext Data
from:https://sites.google.com/site/econometricsacademy/econometrics-models/quantile-regressionReferences:http://yetanothermathprogrammingconsultant.blogspot.com/2021/06/median-quantiles-and-quantile.htmlerwin@amsterdamoptimization.com$offtext *-------------------------------------------------------------* data from csv file*-------------------------------------------------------------$set csv quantile_health.csv $set gdx quantile_health.gdx setsi 'observations'j0 'column names in csv file'; parameter data(i,*) 'all data';$call csv2gdx %csv% output=%gdx% id=data useHeader=T index=1 values=(2..8) $gdxin %gdx% $loaddc i=Dim1 j0=Dim2 data display i,j0,data;*-------------------------------------------------------------* setup of regression data y,X*-------------------------------------------------------------setj 'independent variables' /intercept,suppins,totchr,age,female,white/; parametery(i) 'dependent variable'X(i,j) 'independent variables'; y(i) = data(i,'totexp'); X(i,j) = data(i,j); X(i,'intercept') = 1; *-------------------------------------------------------------* quantile regression LP model (primal)*-------------------------------------------------------------setpm /'+','-'/ ; scalar tau 'quantile';positive variables e(i,pm) 'absolute value';variablesz 'objective variable'beta(j) 'estimates'; equationsobj 'objective'split(i) 'variable splitting'; obj.. z =e= sum(i, tau*e(i,'+') + (1-tau)*e(i,'-'));split(i).. e(i,'+')-e(i,'-') =e= y(i)- sum(j,X(i,j)*beta(j));model quantileLP /obj,split/;*-------------------------------------------------------------* solve loop*-------------------------------------------------------------set q 'quantile levels' /'0.25','0.5','0.75'/;parameter
estimates(q,j);loop(q,tau = q.val; solve quantileLP minimizing z using lp;estimates(q,j) = beta.l(j); ); display estimates;*-------------------------------------------------------------* quantile regression LP model (dual)*-------------------------------------------------------------variable d(i) 'variables of dual problem';equationsobj2 'objective of dual'eqdual(j) 'constraint of dual problem'; obj2.. z =e= sum(i, y(i)*d(i));eqdual(j).. sum(i, x(i,j)*d(i)) =e= 0;* still missing are the bounds on d.* we populate them in the solve loop belowmodel quantileDual /obj2,eqdual/;*-------------------------------------------------------------* solve loop*-------------------------------------------------------------parameter
estimates2(q,j) 'from
dual';loop(q,tau = q.val; d.lo(i) = tau-1; d.up(i) = tau; solve quantileDual maximizing z using lp;estimates2(q,j) = eqdual.m(j); ); display estimates2; |

I did something very similar to "LP model 1 for finding the median" in my bachelor's thesis. I don't think I realized that you could also formulate it by splitting variables.

ReplyDeleteIn most cases it does not really matter which formulation you choose. Here we use the variable splitting approach so we can put different weights on the positive and negative deviation in the subsequent models.

DeleteThe NLP model for finding the median is not differentiable, that is of course correct. But it is convex in m and it has a subgradient defined at every point. A subgradient solver will reliably find the global minimum for that function.

ReplyDeleteThanks. Indeed, the standard solvers I have access to are all gradient based. I could not find off-the-shelf subgradient solvers (I assume users develop them for their problem at hand).

Delete