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

$\endgroup$
  • $\begingroup$ This could be done in a few lines using CVX in Matlab. $\endgroup$ – littleO Nov 9 '13 at 22:49
  • $\begingroup$ 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}\ }$ ?. $\endgroup$ – Felix Marin Nov 9 '13 at 22:52
3
$\begingroup$

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

$\endgroup$
  • $\begingroup$ Thank you for your answer. It solves my problem. $\endgroup$ – ReV Nov 12 '11 at 10:29
  • $\begingroup$ 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? $\endgroup$ – ReV Nov 16 '11 at 12:08
  • $\begingroup$ @Thomas: You sure? Did you compare the results of linprog() with setting things up explicitly with, say fminbnd()? $\endgroup$ – J. M. is a poor mathematician Nov 16 '11 at 12:18
  • $\begingroup$ 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. $\endgroup$ – ReV Nov 16 '11 at 13:46
4
$\begingroup$

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

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – ReV Nov 16 '11 at 16:03
  • 1
    $\begingroup$ I formatted the code line and edited out broken links; but if you happen to have a substitute for them, that would be nice. $\endgroup$ – user147263 Aug 19 '15 at 13:30
  • $\begingroup$ any chance your n's and m's are reversed? (probably my mistake though) $\endgroup$ – Eran Medan Nov 5 '15 at 0:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.