Cholesky Decomposition Column-Wise Algorithm Implementation
43 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
J
le 25 Sep 2019
Réponse apportée : Imane AITSITAHAR
le 8 Avr 2022
Hello I am trying to implement the following algorithm for Cholesky Decomposition Column-Wise Method:
for j=1:n
for i=1:j-1
end
end
My attempt so far at implementing the above:
A=[4 -1 1;
-1 4.25 2.75;
1 2.75 16;];
% Check R matches with col(A);
count = 0;
[n,n] = size(A);
R=zeros(n,n)
for j=1:n
for i=1:j-1
sum1 = 0
for k=1:i-1
sum1 = sum1 + R(k,i)*R(k,j);
end
R(i,j)=(A(i,j)-sum1)/R(i,i);
end
sum2 = 0;
for k=1:j-1
sum2 = sum2 + R(k,j)*R(k,j);
end
R(j,j)=sqrt(A(j,j)-sum2);
end
Q=transpose(R);
S=Q*R;
EDIT: I have modified the code and it runs properly, many thanks to the helpful feedback I received.
3 commentaires
David Goodmanson
le 25 Sep 2019
I takes more than a positive determinant for a symmetric matrix to be positive definite. It also has to have all positive eigenvalues. However,
A = [4 -1 1;
-1 4.25 2.75;
1 2.75 16;];
eig(A)
ans =
2.5946
4.9978
16.6577
so it qualifies.
Réponse acceptée
David Goodmanson
le 28 Sep 2019
Modifié(e) : David Goodmanson
le 28 Sep 2019
Hi J,
I added to the code in your last comment by including the obvious missing 'for' statements, etc.
After that, it's pretty close, with two adjustments needed, one in the code, one in the (supposed) algorithm that you are replicating. First, take a look at
sum1 = 0
for k=1:i-1
sum1 = R(k,i)*R(k,j)
end
The problem here is that sum1 is just set to whatever the last value of R(k,i)*R(k,j) is, and there is no sum of terms. You need to keep a running sum, so replace the sum1 line with
sum1 = sum1 + R(k,i)*R(k,j)
and the same applies to sum2. The resulting code works, almost. The problem is that the algorithm you cited in your original posting is incorrect. It's missing an all-important square root. The expression should be
R(j,j) = sqrt(A(j,j)-sum2)
After that it works. This makes sense since R'*R = A and (for example) when A is a 1x1 scalar, then you are solving R^2 = A for a scalar R. So in general there has to be a square root in there somewhere.
0 commentaires
Plus de réponses (2)
Steven Lord
le 25 Sep 2019
The algorithm you've been given performs a summation twice, once inside both loops and once inside just the outermost loop. Your code does not include the sum function and does not include loops over k.
As a first pass, I recommend writing your code as closely to the algorithm given in your homework / class notes / textbook. [If you're trying to compute the Cholesky decomposition and it's not part of school work, I strongly recommend simply calling chol instead of building your own.] Once you have that working, then you could start modifying it to reduce the number of loops, vectorize some operations, etc.
3 commentaires
Imane AITSITAHAR
le 8 Avr 2022
A=[4 -1 1;
-1 4.25 2.75;
1 2.75 16;];
% Check R matches with col(A);
count = 0;
[n,n] = size(A);
R=zeros(n,n)
for j=1:n
for i=1:j-1
sum1 = 0
for k=1:i-1
sum1 = sum1 + R(k,i)*R(k,j);
end
R(i,j)=(A(i,j)-sum1)/R(i,i);
end
sum2 = 0;
for k=1:j-1
sum2 = sum2 + R(k,j)*R(k,j);
end
R(j,j)=sqrt(A(j,j)-sum2);
end
Q=transpose(R);
S=Q*R;
0 commentaires
Voir également
Catégories
En savoir plus sur Sparse Matrices dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!