A note on variable splitting: complementarity

This note is about a well-known fact in linear programming: variable splitting and minimizing absolute values lead to complementarity pairs \[r_i^{\bf +} \cdot r_i^{\bf -} = 0\] in the solution. A simple example is implementing LAD (Least Absolute Deviation) regression through Linear Programming.

Let's have a look.

Problem: LAD regression

I created an artificial dataset with a few possible outliers:

----     25 PARAMETER data  

                        x           y

observation1        0.442       6.448
observation2        3.823      10.020
observation3        6.231       7.789
observation4        5.727       9.704
observation5        4.607      12.430
observation6        8.779      11.618
observation7        9.277      14.677
observation8        5.934      28.955
observation9        2.727       6.060
observation10       8.390      15.843
observation11       7.665      11.440
observation12       5.798       8.444
observation13       9.096       7.256
observation14       7.773      11.288
observation15       4.459       7.596
observation16       1.593       6.750
observation17       3.951       7.491
observation18       2.772      11.990
observation19       1.752      11.374
observation20       4.911      11.558
observation21       5.613      12.916
observation22       6.168      11.604
observation23       7.551      11.726
observation24       5.057       6.750
observation25       2.547       6.266
observation26       8.157      11.052
observation27       2.398       8.139
observation28       5.162      10.734
observation29       8.787      12.720
observation30       2.949      11.679

The outliers will heavily impact standard least-squares regression, but using absolute deviations reduces its impact.

The LAD-regression problem can be formulated and solved as a simple LP model, using a variable-splitting technique. We split the residual \({\color{darkred}r}_i\) into its positive and negative part: \[{\color{darkred}r}_i = {\color{darkred}r}^+_i - {\color{darkred}r}^-_i\]

The trick is that we can write \[\left|{\color{darkred}r}_i\right| = {\color{darkred}r}^+_i + {\color{darkred}r}^-_i\] This approach only works if at most one of the pair \(({\color{darkred}r}^+_i, {\color{darkred}r}^-_i)\) is nonzero. This is automatically the case for our LAD model as we will discuss below.

The complete LP model is: 

LAD Model

\[\begin{align}\min\>&{\color{darkred}z}=\sum_i \left({\color{darkred}r}^+_i + {\color{darkred}r}^-_i \right) \\ & {\color{darkblue}y}_i ={\color{darkred}\alpha} +{\color{darkred}\beta}{\color{darkblue}x}_i +{\color{darkred}r}^+_i -{\color{darkred}r}^-_i & \forall i \\ &{\color{darkred}\alpha}, {\color{darkred}\beta}\>{\bf free} \\& {\color{darkred}r}^+_i,  {\color{darkred}r}^-_i \ge 0 \end{align}\]

As usual for statistical models, \({\color{darkblue}x}\) and \({\color{darkblue}y}\) are data (exogenous).

Solution

If we take a look at the solution, we see:

----     63 VARIABLE a.L                   =        6.119  intercept
            VARIABLE b.L                   =        0.743  slope

----     63 VARIABLE r.L  residuals

                        +           —

observation2        1.062
observation3                    2.956
observation4                    0.667
observation5        2.890
observation6                    1.020
observation7        1.669
observation8       18.430
observation9                    2.084
observation10       3.495
observation11                   0.370
observation12                   1.980
observation13                   5.617
observation14                   0.603
observation15                   1.834
observation16                   0.552
observation17                   1.562
observation18       3.813
observation19       3.954
observation20       1.792
observation21       2.629
observation22       0.905
observation24                   3.124
observation25                   1.745
observation26                   1.124
observation27       0.239
observation28       0.782
observation29       0.076
observation30       3.370


----     63 VARIABLE z.L                   =       70.345  sum of absolute deviations

Notice that observations 1 and 23 are missing. The display is skipping zero elements. The residuals are zero for these observations.

Also note the typical "complementarity pattern". Only one of the plus or minus components of each residual can be nonzero.

Why does this work?

There are two reasons why we see this complementarity property:

1. The objective pushes \({\color{darkred}r}^+_i + {\color{darkred}r}^-_i\) down. If both were positive, we could subtract the smaller from both and improve the objective. I.e., the solution would not be optimal.
2. This second reason is about how the simplex method works. We can't have both components of the residuals being basic. If we allow that, the basic matrix \(B\) would be singular (not invertible), as the coefficients of \(A\) for these two columns  \({\color{darkred}r}^+_i\) and \({\color{darkred}r}^-_i\) are identical except for the sign. Since only one of these two columns can be basic, the variable corresponding to the other must be zero (i.e., non-basic). A caveat is using an interior point method without crossover or a first-order method. Of course, they can still rely on the objective. Also, I assume here that the objective row is not part of the matrix \(A\).

I typically don't see reason 2 being mentioned. Maybe a singular basis matrix is a little bit too technical.

 

To verify that reason number 2 is indeed working, try solving the problem with a different objective. E.g. just \(\min\>{\color{darkred}z}=0\). Then check the complementarity in the solution.

 

Notes:

1. Using the alternative bounding approach for LAD regression \[\begin{align}\min & \sum_i  {\color{darkred}t}_i \\ & -{\color{darkred}t}_i \le {\color{darkred}r}_i \le {\color{darkred}t}_i \\ &  {\color{darkblue}y}_i ={\color{darkred}\alpha} +{\color{darkred}\beta}{\color{darkblue}x}_i +{\color{darkred}r}_i \\ & {\color{darkred}t}_i \ge 0 \end{align}\] does not have this issue or property. 
2. Of course, minimizing absolute values happens in applications outside of LAD regression. 
3. If we were to maximize the absolute values, we would no longer have a convex problem. In that case, we would need heavy machinery such as binary variables or SOS1 sets. 

References

1. Least absolute deviations, https://en.wikipedia.org/wiki/Least_absolute_deviations

Appendix: GAMS Implementation 

$ontext         LAD (Least Absolute Deviation) Regression          LP model using variable splitting    $offtext   option lp=highs;   option seed=1234;   *--------------------------------------------------------------------- * data *---------------------------------------------------------------------   set i 'observations' /observation1*observation30/;   parameter data(i,*);   data(i,'x') = uniform(0,10); data(i,'y') = 10+0.25*data(i,'x')+normal(0,3); * for a few observations, generate more wild y data data(i,'y')$(uniform(0,1)<0.1) = 10+0.25*data(i,'x')+normal(0,10); display data;   *--------------------------------------------------------------------- * derived data *---------------------------------------------------------------------   parameters     x(i) 'independent variable'     y(i) 'dependent variable' ; x(i) = data(i,'x'); y(i) = data(i,'y');   *--------------------------------------------------------------------- * LAD regression LP model *---------------------------------------------------------------------   set pm 'plus/minus' /'+','-'/;   variables     a 'intercept'     b 'slope'     z 'sum of absolute deviations' ; positive variables     r(i,pm) 'residuals' ;   equations     objective 'objective function'     fit(i)    'model to be fitted' ;   objective.. z =e= sum((i,pm), r(i,pm)); fit(i)..    y(i) =e= a + b * x(i) + r(i,'+') - r(i,'-');   model lad /all/; solve lad using lp minimizing z; display a.l, b.l, r.l, z.l;