MATLAB: Cholesky Decomposition Column-Wise Algorithm Implementation

cholesky decompositionMATLAB

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.

Best Answer

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.
Related Question