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!
采纳的回答
Matt J
2013-11-10
With only 3 unknowns, fminsearch should work fairly well
fminsearch(@(x) norm(A*x-b,1), x0 )
14 个评论
abc
2013-11-10
abc
2013-11-12
abc
2013-11-12
Matt J
2013-11-14
And, as of today, it is now possible to post .mat files unzipped, see
abc
2013-11-17
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
2014-4-6
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.
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...
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.
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.
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
2023-1-26
The problem is the way you have defined matrix A. It must be:
A=[Ai1 Ai2];
Then it works!
更多回答（1 个）
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)
12 个评论
It uses whichever one of linprog's algorithms you select in the options input. The default, I believe, is 'dual-simplex'.
NA
2021-11-8
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?
NA
2021-11-8
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);
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
2021-11-8
I think, but how can I use this objective function?
Matt J
2021-11-8
Well, it's just a direct application of minL1lin with
C = inv(R)^0.5 * A;
d = inv(R)^0.5 * y;
NA
2021-11-9
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
2021-11-9
If W is dagonal, then yes, it's true.
NA
2021-11-10
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?
类别
在 帮助中心 和 File Exchange 中查找有关 Fit Postprocessing 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
