Cholesky Decomposition Column-Wise Algorithm Implementation

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.

[Fabio Freschi]

When j = 1, you have to evaluate R(1:j-1,j) but the indexing of rows is 1:0 that is

  1×0 empty double row vector  

Note also that your A matrix is not SPD, so Cholesky cannot be applied

[J]

Re: Fabio Freschi

When j = 1, you have to evaluate R(1:j-1,j) but the indexing of rows is 1:0 that is

1×0 empty double row vector

Note also that your A matrix is not SPD, so Cholesky cannot be applied

I don't understand why is A not Positive Semi Definite, it has positive deteriminant and symmetric. Any help with how to directly implement the last line of the original algorithm seems to be the answer I am looking for, any help with this would be appreciated.

[David Goodmanson]

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.

Answer 1

[David Goodmanson]

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.

Answer 2

[Steven Lord]

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.

[Steven Lord]

Your modified code is overwriting R(i, j) every iteration through the two loops over k. Before you enter the first for loop over k, make a temporary value to store the sum. Inside the for loop over k, add R(k, i)*R(k, j) to your temporary value. After the for loop is finished, then update R(i, j) using the sum stored in the temporary value.

Use that same technique for the second loop as well.

[J]

Re: Steven Lord

Thanks you for your insights on my problem, I have again modified the two lines dealing with the sums in the code:

        sum1 = 0
        for k=1:i-1
        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 = (R(k,j)*R(k,j));
    end
    R(j,j)=A(j,j)-sum2

The resultant matrix it produces is now the expected size but is unfortunately still different from directly invoking the chol(A) command. Am I implementing what you're suggesting correctly or no? Further help would be appreciated, thanks.

[J]

Thank you for your feedback and help, I will update the original post to include a working code.

Answer 3

[Imane AITSITAHAR]

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;