least-absolute-deviationslinear programmingoptimizationquantile regressionregression

Here is the least absolute deviation problem under concerned: $ \underset{\textbf{w}}{\arg\min} L(w)=\sum_{i=1}^{n}|y_{i}-\textbf{w}^T\textbf{x}|$. I know it can be rearranged as LP problem in following way:

$\min \sum_{i=1}^{n}u_{i}$

$u_i \geq \textbf{x}^T\textbf{w}- y_{i} \; i = 1,\ldots,n$

$u_i \geq -\left(\textbf{x}^T\textbf{w}-y_{i}\right) \; i = 1,\ldots,n$

But I have no idea to solve it step by step, as I am a newbie to LP. Do you have any idea? Thanks in advance!

EDIT:

Here is the latest stage I have reached at this problem. I am trying to solve the problem following this note:

Step 1: Formulating it to a standard form

$\min Z=\sum_{i=1}^{n}u_{i}$

$ \textbf{x}^T\textbf{w} -u_i+s_1=y_{i} \; i = 1,\ldots,n$

$ \textbf{x}^T\textbf{w} +u_i+s_2=-y_{i} \; i = 1,\ldots,n$

subject to $s_1 \ge 0; s_2\ge 0; u_i \ge 0 \ i=1,…,n$

Step 2: Construct a initial tableau

           |      |    0      |    1   |  0  |   0   |   0    
 basic var | coef |  $p_0$    |  $u_i$ |  W  | $s_1$ | $s_2$ 
      $s_1$| 0    |  $y_i$    |   -1   |  x  |   1   |   0
      $s_2 | 0    |  $-y_i$   |    1   |  x  |   0   |   1
      z    |      |    0      |    -1  |  0  |   0   |   0

Step 3: Choose basic variables

$u_i$ is chosen as input base variable. Here comes a problem. When choosing the output base variable, it is obvious $y_i/-1=-y_i/1=-y_i$. According to the note, if $y_i\ge0$, the problem has unbounded solution.

I am totally lost here. I wonder if there is anything wrong and how should I continue the following steps.

Best Answer

You want an example for solving least absolute deviation by linear programming. I will show you an simple implementation in R. Quantile regression is a generalization of least absolute deviation, which is the case of the quantile 0.5, so I will show a solution for quantile regression. Then you can check the results with the R quantreg package:

    rq_LP  <-  function(x, Y, r=0.5, intercept=TRUE) {
        require("lpSolve")
        if (intercept) X  <-  cbind(1, x) else X <-  cbind(x)
        N   <-  length(Y)
        n  <-  nrow(X)
        stopifnot(n == N)
        p  <-  ncol(X)
        c  <-  c(rep(r, n), rep(1-r, n), rep(0, 2*p))  
                 # cost coefficient vector
        A  <- cbind(diag(n), -diag(n), X, -X)
        res  <-  lp("min", c, A, "=", Y, compute.sens=1)
            ### Desempaquetar los coefs:
        sol <- res$solution
        coef1  <-  sol[(2*n+1):(2*n+2*p)]
        coef <- numeric(length=p)
        for (i in seq(along=coef)) {
             coef[i] <- (if(coef1[i]<=0)-1 else +1) *  
                  max(coef1[i], coef1[i+p])
        }
        return(coef)
        }

Then we use it in a simple example:

    library(robustbase)
    data(starsCYG)
    Y  <- starsCYG[, 2]
    x  <- starsCYG[, 1]
    rq_LP(x, Y)
    [1]  8.1492045 -0.6931818

then you yourself can do the check with quantreg.