Least Absolute Value based regression

Hi ,

I want to use linear regression based on least absolute value deviation to find the coefficients of my model with the help of measured data and 3 independent variables. The number of measurements I have are much greater than 3( number of cofficients). I have been trying to implement LAV method for this, but with no success. I need urgent help. Can someone please guide me in this. If someone already has the code for LAV, I would be grateful if you could share it with me.

Thank you!

Answer 1

[Matt J]

With only 3 unknowns, fminsearch should work fairly well

fminsearch(@(x) norm(A*x-b,1), x0 )

[abc]

That is what I expected but it is not working. In my case A is a say 1500*3 matrix (assuming 1500 measured data points), x is a 3*1 matrix of the coefficients and b is a 1500*1 vector. When I give the initial values in x0, fminsearch simply returns to me the initial values. Is there a way around it?

[Matt J]

I don't see why that would be happening. How about you attach your A,b data in a .mat file so we can try ourselves.

[abc]

The excel files are the A and b matrices. A MATLAB file cannot be uploaded here. There is one file with both A and b combined. Thanks a lot for the help!!

[abc]

I apologize for the mistake but the file Matrix A. csv has both the A and b vectors. I am attaching the file with the A matrix alone if it helps. The file A and b.csv has both the matrices.

[Matt J]

Please put A and b in a .mat file and attach it. If you .zip the mat file, you can upload it.

[Matt J]

And, as of today, it is now possible to post .mat files unzipped, see

http://www.mathworks.com/matlabcentral/answers/88363#comment_179826

[abc]

Hi Matt, Thank you for informing about the update. I have uploaded the A and b matrices in the attached .mat file. Thank you once again for your help.

[Matt J]

Shruti,

The solution I get with fminsearch, and which I believe is correct, is

x =

     -0.0014
     -0.0383
     -0.7048

A more rigorous way to solve the problem is using linprog as below.

    [m,n]=size(A);
    f=[zeros(1,n), ones(1,m)];
    Ac=[A,-speye(m);-A,-speye(m)];
    bc=[b;-b];

    lb(1:3)=-inf; lb(4:m+n)=0;

    xz=linprog(f,Ac,bc,[],[],lb,[]);

    x=xz(1:3);

In general, LINPROG is a more reliable solver than fminsearch, but I get the same answer this way as with fminsearch, and the computation using fminsearch is faster.

[Puneet]

Hello Matt,

I am implementing the similar algorithm for my problem I have matrix A and b in formulation Ax = b,

min Ax-b, such that x>0

What I understood from your algorithm is the rest of the variables except x(basic) are slack variables(4 to m+n) in the formulation of (when you break equality as:)

A*x <= b

-A*x<= -b

So I replaced the formulation to Ac = [A,speye(m);-A,speye(m)] and in my problem the variables had a lower bound of 0. This modification returned me better estimate of variables.

What is your opinion on this modification? If I am correct. Please correct me as I am just a beginner in optimization.

Thanks in advance.

[Matt J]

Hi Puneet,

No, it's not clear to me that the auxiliary variables can be interpreted as slack. Slack variables are used to convert inequalities to equalities, whereas my formulation contains only inequalities.

As far as I can see, the only change you've made is to flip the sign of speye(m), which is equivalent to changing the sign of the auxiliary variables in the solution. It would be equivalent to redefining the objective weight vector as

 f=[zeros(1,n), -ones(1,m)];

But this, in turn, is the same as flipping the sign of f. So, ultimately, it looks like you are maximizing the same objective function that I was minimizing. Not sure how maximizing instead of minimizing led to better results...

[Puneet]

Hello Matt,

Thanks for your reply. I have not changed the sign of f in your formulation. Can you tell me the reference as how you came up with your formulation with the auxiliary variables? I am not sure on the use of auxiliary variables. I just flipped the sign of speye(m) assuming it to be slack variables. Since I tested on a small system and it gave me very good results but before I apply this to very large formulations I need to be sure if I am fundamentally correct in formulation.

Thanks in advance.

[Matt J]

The problem

   min_x  norm(A*x-b,1)

is equivalent to the constrained problem,

   min_xz sum(z)

   subject to

    -z <= A*x-b <= z

Here z is a vector of auxiliary variables. The string of inequalities -z <= A*x-b <= z when re-arranged leads to

   Ac*[x;z]<=bc

where Ac and bc are defined as earlier and sum(z) can be implemented as dot(f,[x;z]) using the f I also defined earlier.

[NA]

Hello Matt,

I have a regression model like:

zi = Ai1*x1+Ai2*x2+ei i=1,2,3,4,5

The observation are given as table below:

i = [1;2;3;4;5];
zi = [-3.01;3.52;-5.49;4.03;5.01];
Ai1 = [1;0.5;-1.5;0;1];
Ai2 = [1.5;-0.5;0.25;-1;-0.5];

AA = [i,zi,Ai1,Ai2];
T = array2table(AA);
T.Properties.VariableNames(1:4) = {'i','z_i','A_i1','A_i2'};

I used this code for LAV estimation

A = Ai1+Ai2;
b = zi;

len_x = size(A,2);
c = [zeros(len_x,1);ones(size(b))];
F = [A -eye(size(b,1)); -A -eye(size(b,1))];
g = [b; -b];

z = linprog(c,F,g);
xlav = z(1:2);

The result is not correct as xlav and residual should be

xlav = [3.005;-4.010]
residual = [0; 0.0125; 0.2; 0.02; 1]

[Amin Nassaj]

The problem is the way you have defined matrix A. It must be:

A=[Ai1 Ai2];

Then it works!

Answer 2

[Matt J]

@NA In the years after this thread, I implemented minL1lin, which you can download from

https://www.mathworks.com/matlabcentral/fileexchange/52795-constrained-minimum-l1-norm-solutions-of-linear-equations?s_tid=srchtitle

A quick test shows that it gives your expected answer:

zi = [-3.01;3.52;-5.49;4.03;5.01];
Ai1 = [1;0.5;-1.5;0;1];
Ai2 = [1.5;-0.5;0.25;-1;-0.5];

[xlav,~,res]=minL1lin([Ai1,Ai2],zi)

Optimal solution found.

xlav = 2×1

3.0050 -4.0100

res = 5×1

0.0000 -0.0125 -0.0200 -0.0200 0.0000

[NA]

Thank you for your time. Is this a 'Simplex Based Algorithm'?

[Matt J]

It uses whichever one of linprog's algorithms you select in the options input. The default, I believe, is 'dual-simplex'.

[NA]

I have a regression model like:

Z=Ax+e

[z1;z2;z3] = [A11 A12 A13 A14; A21 A22 A23 A24; A31 A32 A33 A34]*[x1;x2;x3;x4]

where, Z is m*1 vector, A is m*n matrix, and x is n*1 vector.

If we have covariance matrix R (m*m)

In this context,

[xlav,~,res]=minL1lin([A],Z)

But how we use covariance matrix in this function?

[Matt J]

What is the intended computation, and how would the covariance matrix be involved?

[NA]

What is the intended computation,

For state estimation.

how would the covariance matrix be involved?

for weighted least square the covariance matrix involved as:

F = A' * inv(R) * (z - z_est);
J = A' * inv(R) * A;
dx = (J \ F);

[Matt J]

Yes, but what would be the covariance-weighted objective function that you are considering in the case of L1 minimization? Would it be

norm(inv(R)^0.5 * (A*x-y),1)

[NA]

I think, but how can I use this objective function?

[Matt J]

Well, it's just a direct application of minL1lin with

C = inv(R)^0.5 * A;
d = inv(R)^0.5 * y;

[NA]

I checked the formulation, the weighted least absolute value minimizes this cost function:

r is the residual vector anf Wj is reciprocal of the covariance of entry j.

If I have W, is this true?

C = W * A;
d = W * y;

[Matt J]

If W is dagonal, then yes, it's true.

[NA]

Thank you. If we change one of the entry of matrix A, such that

A_d = [A11 A12 A13 A14; A21 4*A22 A23 A24; A31 A32 A33 A34]

When we calculate residuals, how minL1lin removes the entry to find the regression. I mean, what is the threshold level that minL1lin use for regression?

[Matt J]

I'm not really following you. What leads you to believe that elements of A are being removed?