6
$\begingroup$

I would like to solve a linear system $Ax = b$ under the $L_1$ norm constraint $\min(|Ax-b|)$. All that I can find about $L_1$ minimization is a way to minimize $|x|_1$ subject to $Ax=b$.

I wanted to use linear programming in matlab to solve this problem. This lead me to solving a new system with linprog in matlab:

$$My = k$$

So I did some transformations, knowing that Linprog can solve the following problem:

$$\min f(x) \ \ \textrm{s.t.} \ \ Ax \leq b$$

I want to minimize this problem:

$$\min ||My - k||_1$$

for $y$

To convert this to the linprog form, I introduced a new vector I want to minimize, $z$:

$$\min \sum z$$

s.t.

$$My - k \leq z$$ $$-My + k \leq z$$

However, the constraints include both the variable $y$ and the constants $k$. We can formulate this in the required form by creating a bigger system:

$Ax \leq b$, where the matrix $A$ is:

$$\begin{pmatrix} M & -I \\ -M & I \end{pmatrix} $$

the vector $x$ is equal to $y$ and $z$ stacked: $$\begin{pmatrix} y \\ z \end{pmatrix} $$

and the constraints are $k$ and $-k$ stacked: $$\begin{pmatrix} k \\ -k \end{pmatrix} $$

The problem is that this doesn't work. (Matlab can't find a good solution). I would like to know if there is another way to solve this problem? If you have heard about a matlab code that already does it?

Thanks a lot

  • 0
    This could be done in a few lines using CVX in Matlab.2013-11-09
  • 0
    What is your norm $\large\left\vert\left\vert\cdots\right\vert\right\vert$. Is it $\large\left\vert\left\vert{\bf x}\right\vert\right\vert \equiv \sqrt{\vphantom{\Large A}{\bf x}^{\sf T}{\bf x}\ }$ ?.2013-11-09

2 Answers 2

3

You forgot about lower bound constraints $z \ge 0$. After I added this constraint to LP linprog managed to solve your problem.

  • 0
    Thank you for your answer. It solves my problem.2011-11-12
  • 0
    I am going back to this question again. Also your answer helped me solve the problem I described, it seems that the way I derived the new system is not good. In fact, it always solve the system under the $L_2$ norm. What I am looking for is really minimizing $|Ax-b|$. Do you have any idea how to minimize the system this way?2011-11-16
  • 0
    @Thomas: You sure? Did you compare the results of `linprog()` with setting things up explicitly with, say `fminbnd()`?2011-11-16
  • 0
    I did some comparisons after my previous post. I think it finally works... It is just, that for my problem, the solution of both minimisation algorithms ($L_1$ and $L_2$) are almost the same. But it seems correct (even if I didn't expect that) Thank you a lot.2011-11-16
4

If you have min $\|Ax-b\|_1$ where $A$ is $n\times m$ for $n>m$ and $b$ is $n\times 1$ (so the sum of the absolute values of all the components of the result vector $y=Ax-b$ where y is $n\times 1$) then you can solve this as an LP with linprog with a command like this:

linprog([zeros(m,1);ones(n,1)],[+A,-eye(n);-A,-eye(n)],[b;-b])

which is just a translation into Matlab of the usual trick described on Wikipedia.

Here is a polynomial fitting example that compares least absolute deviations to linear least squares.

clc
x = [0:1/100:1]';

%random polynomial of degree 4
poly = [7/960;-1/8;227/320;-23/15;1]/0.0583;
A = [x.^0,x.^1,x.^2,x.^3,x.^4];
b = A*poly;

% bigger k means less outliers. k = 2 has too many outliers
% for the least absolute deviations solution to lie on the initial
% polynomial
k = 3;
b(k*[1:101/k]) = x(k*[1:101/k]);

% linear least squares
poly = A\b;
y = A*poly;

% least absolute deviations
[n,m] = size(A);
poly = linprog([zeros(m,1);ones(n,1)],[+A,-eye(n);-A,-eye(n)],[b;-b]); poly = poly(1:m);
z = A*poly;

% blue lines: original data
% green lines: least squares solution
% magenta crosses: least absolute deviations solution
% red circles: outliers
plot(x,b,'b-',x,y,'g-',x,z,'m+',x(k*[1:101/k]),x(k*[1:101/k]),'ro')

% check least absolute deviations
check = fminsearch(@(x) norm(A*x-b,1), poly);
[[poly; norm(A*poly-b,1)], [check; norm(A*check-b,1)]]

This is what the output looks like enter image description here

  • 0
    Thanks for this answer. Actually, it shows me that I have to use this matrix: $$\begin{pmatrix} M & -I \\ -M & -I \end{pmatrix} $$ instead of $$\begin{pmatrix} M & -I \\ -M & I \end{pmatrix} $$ Your code is looking exactly like mine after this transformation. The result of my problem is slightly different but still of the same order of magnitude. Thanks for the help.2011-11-16
  • 1
    I formatted the code line and edited out broken links; but if you happen to have a substitute for them, that would be nice.2015-08-19
  • 0
    any chance your n's and m's are reversed? (probably my mistake though)2015-11-05