least absolute deviation when we have data set

I have this data

x = (1:10)';
y = 10 - 2*x + randn(10,1);
y(10) = 0;

how can I use least absolute value regression?

Answer 1

[Bjorn Gustavsson]

You can do it rather straight-forwardly with fminsearch (or other similar tools on the file exchange: fminsearchbnd, minimize etc):

M = [ones(size(x)),x]; % Matrix for linear LSQ-regression, we could do centering and scaling etc...
p0 = M\y;           % straight least-square fit - to get ourselves a sensible start-guess (hopefully)
errfcn = @(p,y,M) sum(abs(y-M*p)); % L1 error-function 
p1 = fminsearch(@(p) errfcn(p,y,M),p0); % L1-optimization
subplot(2,1,1)
plot(x,y,'.-')
hold on
plot(x,M*p0)
plot(x,M*p1)
subplot(2,1,2)
plot(x,y*0,'.-')
hold on
plot(x,y-M*p0,'.-')
plot(x,y-M*p1,'.-')
% For my test I got L1-error-function-value for the least-square-fit p0:
% errfcn(p0,y,M)
% ans =
%        22.058
% and for the L1-optimal parameters:
% >> errfcn(p1,y,M)
% ans =
%       20.067

This would generalize to more interesting problems too. Also have a look at Huber-norms, for an error-norm kind of intermediate between L1 and L2.

HTH

[NA]

Thank you. For error-norm kind of intermediate between L1 and L2, I add this line

plot(x,(M*p0+M*p1)/2)

is it correct?

[Bjorn Gustavsson]

When I use Huber-norms I typically do just as above but with a modification of the error-function to use thepseudo-Huber-loss-function instead of the absolute deviation. That is possible to use independently of the stats-toolbox. It will require a bit more work in terms of calculating the statistics of the fit and such-like.

[NA]

How can I calculate residual of least absolute deviation in fminsearch function?

[Bjorn Gustavsson]

If you take a look at the error-function, errfnc, above:

errfcn = @(p,y,M) sum(abs(y-M*p));

You see that it will calculate the sum of the absolute value of the residuals (y-M*p), in order to calculate the

residuals you will simply have to calculate them after fminsearch finished:

resL1 = y-M*p1;

[NA]

Thank you for your answer.

I used this code for pseudo-Huber-loss-function

errfcn = @(p,y,M) ((abs(y-M*p)).^2).*(sqrt(1+((y-M*p)/abs(y-M*p)).^2)-1); % pseudo-Huber-loss-function 

But I do not think this is correct. Would it be possible to tell me how I can do it for Huber.

[Bjorn Gustavsson]

My Huber-like-norm function looks something like this:

function errs = normHuberlike(y,ymod,sigma)
errs = sigma.^2.*(sqrt(1+(y-ymod).^2./sigma.^2)-1)
end

So that defines how the errors in errs vary with the scaled deviation. Then you'll have to use that in your error-function, something like this:

M = [ones(size(x)),x]; % Matrix for linear LSQ-regression, we could do centering and scaling etc...
p0 = M\y;           % straight least-square fit - to get ourselves a sensible start-guess (hopefully)
sigma = 1; % just an arbitrary setting for where the transition between quadrtic and linear behaviour should be
errfcn = @(p,y,M,sigma) sum(normHuberlike(y,M*p,sigma)); % L1 error-function 
p1 = fminsearch(@(p) errfcn(p,y,M,sigma),p0); % L1-optimization

One thing you should "always" do is to plot your norm-function to get a sense of how they weigh residuals of different magniture:

dy = linspace(-4,4,401);
plot(dy,[abs(dy);dy.^2/2;normHuberlike(dy,0*dy,1)])

[NA]

Thank you for taking the time to answer my question. I used your code and compare with 'robustfit'.

Why I could not get same regression?

[Bjorn Gustavsson]

Because they use different algorithms, and from the robust-fit documentation you can look up the weighting used for its different settings of wfun and tune. Do the regressions differ by much? How do they vary if you vary the different tuning-parameters? When using robust fitting you should always check the residuals and their relative contributions to the total error-function.

Answer 2

[Bruno Luong]

% Test data
x = (1:10)';
y = 10 - 2*x + randn(10,1);
y(10) = 0;

order = 1; %  polynomial order
M = x(:).^(0:order);
m = size(M,2);
n = length(x);
Aeq = [M, speye(n,n), -speye(n,n)];
beq = y(:);
c = [zeros(1,m) ones(1,2*n)]';
%
LB = [-inf(1,m) zeros(1,2*n)]';
% no upper bounds at all.
UB = [];

sol =  linprog(c, [], [], Aeq, beq, LB, UB);
Pest = sol(m:-1:1); % here is the polynomial

% Check
clf(figure(1));
plot(x, y, 'or', x, polyval(Pest,x), 'b');

[Terry nichols]

I have 1 comment, 1 question, and 1 interestig result

(actually I have many qustions, but I'll just keep it to one for now)

1) Nice piece of work! I tried it for my problem and it worked seamlessly. The L1 norm just blew through all of the outliers. Although I am still trying to understand how you calcualted the values that need to get passed to linprog.m (my long list of questions).

2) I am totally new to the ways of linear programming, so I am wondering how come you have no inequality constraints? I am guessing you are saying the solution must adhere to the objective function, precisely?

3) I tested your code for a quadratic using both the L2 (polyfit, polyval) and L1 norm provided in your code and got the following.

[Bruno Luong]

"2) I am totally new to the ways of linear programming, so I am wondering how come you have no inequality constraints? I am guessing you are saying the solution must adhere to the objective function, precisely? "

Because I don't need it. I formulate the problem as

M*P - u + v = y

where u and v a extra variables, they meant to be positive

v =( M*P - y) = u

so

argmin (u + v)  is sum(abs( M*P - y)) is L1 norm of the fit.

I could formulate with inequality but they are equivalent. There is no unique way to formulate LP, as long as it does what we want.

And as comment; all LP can be showed to be equivalent to a "canonical form" where all the inequalities are replaced by only linear equalities + positive bounds

argmin f'*x
    A*x = b
    x >= 0

[Terry nichols]

Thanks much for all of your help!