How to perform mathematical operation of arrays inside loop?

2 views (last 30 days)
I know it is a very basic question, but I have not been able to solve it. I have an array X(28*1). Using values of X for the first time, I would calculate E=B-A*X where E=(352*1), A(352*28) and B=(352*1). Then if the sum of E is larger than 10^-40, The follwoing operations are done:
C=diag(E*E');
W=diag(C)^-1;
X=(A'*W*A)^-1*A'*W*B
Then the updated value of X will be used for E=B-A*X operation. I understand whole thing is required to be inside a loop. But I can not find the exact process to write the code. If anyone has any idea, It would be much appriciated. I have attached the X for 1st iteration, A, and B as .mat file.
  6 Comments
Bruno Luong
Bruno Luong on 19 Jul 2023
If you want to do robust linear fit with outliers, I recommend you do L1 fitting
minimize norm(A*X - B, 1)
rather than weighted 2-norm (more robust than l2, but still not very robust).
This requires linear programing linprog of optimization toolbox.
Matt J
Matt J on 19 Jul 2023
Edited: Matt J on 19 Jul 2023
I recommend you do L1 fitting...This requires linear programing linprog of optimization toolbox.
But, I'd like to add that there would be no need to implement it from scratch. You could instead use minL1lin() from this FEX download,

Sign in to comment.

Answers (3)

Matt J
Matt J on 18 Jul 2023
Edited: Matt J on 18 Jul 2023
for i=1:N
E=B-A*X;
if norm(E,1)<1e-40, continue; end
C=E.^2;
X=lscov(A,B,1./C);
end
  3 Comments
Matt J
Matt J on 20 Jul 2023
Edited: Matt J on 20 Jul 2023
@ANANTA BIJOY BHADRA I notice that the second column of your A matrix is zero, which iof course makes A ill-conditioned.
A=load('A').A_primary_A;
[l,h]=bounds(A(:,2))
l = 0
h = 0
cond(A)
ans = 2.5955e+16
cond(A(:, [1,3:end])) %skip second column
ans = 1.4142
Bruno Luong
Bruno Luong on 20 Jul 2023
Edited: Bruno Luong on 20 Jul 2023
Good catch. But I beleive the solution of "\" or lscov, both based on QR decomposition would not be affected by 0 column and remain accurate, despite the warning. The solution corresponds to 0-column is simply nil.
A=rand(10,3);
b=rand(10,1);
As=A; As(:,end+1)=0;
w = 0.1+rand(size(b));
format long
A\b
ans = 3×1
0.305726431573425 0.655213115202102 -0.126807101744095
As\b
Warning: Rank deficient, rank = 3, tol = 4.365185e-15.
ans = 4×1
0.305726431573425 0.655213115202102 -0.126807101744095 0
lscov(A,b,w)
ans = 3×1
0.365558371139105 0.511882753808717 -0.133776410641447
lscov(As,b,w)
Warning: A is rank deficient to within machine precision.
ans = 4×1
0.365558371139105 0.511882753808717 -0.133776410641447 0

Sign in to comment.


Mrutyunjaya Hiremath
Mrutyunjaya Hiremath on 18 Jul 2023
E = B-A*X;
while(sum(E(:))> 10^-40)
C=diag(E*E');
W=diag(C)^-1;
X=(A'*W*A)^-1*A'*W*B;
end

Voss
Voss on 18 Jul 2023
E = B-A*X;
while sum(E,'all') > 1e-40
C = diag(E*E');
W = diag(C)^-1;
X = (A'*W*A)^-1*A'*W*B;
E = B-A*X;
end

Categories

Find more on Programming in Help Center and File Exchange

Products


Release

R2023a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!