Is There an $ {L}_{1} $ Norm Equivalent to Ordinary Least Squares?

Question

The ordinary least squares (OLS) method is very useful. It gives you the solution to the problem

$$ \arg \min_{x} {\left\| A x - b \right\|}_{2}^{2} $$

Now, if the problem is the same, but the $1$-norm is used instead:

$$ \arg \min_{x} {\left\| A x - b \right\|}_{1} $$

Is there a known (approximate or not) solution to that problem? Any time efficient algorithm to get this optimum? I've read of the Theil-Sen estimator which should do the trick in dimension $2$, and some multidimentionnal extension of it, but the algorithm computation time increases hugely with dimension, I don't think I'll get any solution before a year if I use that.

Answer 1

The problem is given by:

$$ \arg \min_{x} {\left\| A x - b \right\|}_{1} $$

If we define $ r = A x - b $ we can rewrite the problem as:

$$\begin{align*} \arg \min_{x, t} \quad & \boldsymbol{1}^{T} t \\ \text{subject to} \quad &-t \preceq A x - b \preceq t \end{align*}$$

This can be easily solved by any Linear Programming solver.

You may see some more related approaches in my answer to Mathematics Q1639716 - How Can $ {L}_{1} $ Norm Minimization with Linear Equality Constraints (Basis Pursuit / Sparse Representation) Be Formulated as Linear Programming?

Answer 2

Not equivalent, but one approach is iteratively re-weighted least squares, where you solve a sequence of ordinary least squares problems

$$x_{i}= \min_x\{\|W_i(Ax-b)\|_2\}$$

And re-weighting means to calculate $W_{i+1}$ based on $x_i$

Basically a loop

1. Calculate $W_k$ based on previous solution
2. Calculate $x_k$
3. If $\|x_k-x_{k-1}\|>\epsilon$ loop back to 1.

Answer 3

The answer seems to be: The problem is not particularly hard to solve, but there are no standard solutions. Though there are many ways to do this, I feel like @Royi's anwer is most "correct", since it reduces the problem to a simple linear programming problem, as opposed to the more general IRLS-solution which can be used to solve much harder problems. However, @royi's answer is short and can be a bit confusing, so I wanted to elaborate a bit.

Firstly, why is $$ \arg \min_x \|Ax-b\|_1 $$ equivalent to $$ \begin{align*} \arg \min_{x, t} \quad & \boldsymbol{1}^{T} t \\ \text{subject to} \quad &-t \preceq A x - b \preceq t \end{align*} $$ ?

Because $\|Ax-b\|_1$ is the sum of the absolute values of the vector $r=Ax-b$, and, at the optimal solution, $t_{opt}=|r|$ (where i use $|\cdot|$ to denote the element-wise absolute value). It requires a little bit of thought, but you can probably convince yourself that $-t \preceq r \preceq t$ is equivalent to $|r| \preceq t$. So the elements of $t$ must be equal to or greater than the absolute value of the corresponding elements of $r$, but, since we are minimizing the sum of the elements of $t$, they must take the smallest value possible, therefore they are not greater than but equal to the elements of $|r|$ at the solution $t_{opt}$.

Secondly, how do we convert this into a standard linear programming problem, where we are optimizing for $x$?

According to the problem statement, the cost function $\mathbf{1}^Tt$ should be optimized for both $x$ and $t$, but the function only contains $t$. For a practical programming approach, we need to stack $x$ and $t$ into a single vector $q=[t^T\quad x^T]^T$ which will be the optimization target. We can then rewrite the objective function as $$ \arg\min_q \left[\begin{array}{c}\mathbf{1}\\ \mathbf{0}\end{array}\right]^Tq. $$ We also need to rewrite the constraints in terms of $q$. Treating the upper and lower bounds separately, we have $-t\preceq Ax-b \implies -t -Ax\preceq -b\implies [-I\quad -A]q\preceq -b$. Equivalently, we have $Ax-b\preceq t\implies Ax-t\preceq t\implies [-I\quad A]q \preceq b$.These can be combined to a single inequality $$ \left[\begin{array}{cc}-I & -A \\ -I & A\end{array}\right]q\preceq \left[\begin{array}{c}-b \\ b \end{array}\right] $$ This form of the objective function and constraints in terms of $q$ can be plugged into a linear programming solver, and the optimal $x$ can be extracted as the final elements of the optimal $q$.

The following python code does this using scipy.optimize.linprog:

import scipy.optimize as opt
def linear_minimum_mean_absolute_error(A, b):
    n = len(b) 
    m = A.shape[1]
    c = np.zeros(n + m)
    c[0:n] = 1.0
B = np.concatenate((np.concatenate((-np.eye(n), -A), 1), 
                    np.concatenate((-np.eye(n), A), 1)), 0)
w = np.concatenate((-b, b))

res = opt.linprog(c, A_ub=B, b_ub=w, bounds = None)

return res.x[n:]