Cholesky decomposition

Question

This code is Cholesky decomposition in numerical linear algebra.

function LLT = Cholesky(A)

    if is_pos_def(A) == 1
        disp('Cholesky Decomposition is used only for positive definite matrixes')
    else
        n = rank(A);
        L = zeros(rank(A));
        L(1,1) = sqrt(A(1,1));
        for i=2:n
            L(i,1)=A(i,1)/L(1,1);
        end
    
        for i = 2:n
        
            L(i,i)= sqrt(A(i,i)- sum(power(L(i,1:i-1),2)));

            for j = i+1:n

                L(j,i) = (A(i,j)-(dot(L(j,1:i-1), L(i,1:i-1))))/L(i,i);
            end
        
        end
    end
    LLT = L;
end

I'd prefer to use vectorized code instead of the second for loop.

Is it possible to replace the inner (j) loop?

Answer 1

First I'll point out some general code improvements, then we'll get into vectorization.

There is a bug in your code, in case the first condition is true:

function LLT = Cholesky(A)
    if is_pos_def(A) == 1
        disp('Cholesky Decomposition is used only for positive definite matrixes')
    else
        % ...computations...
    end
    LLT = L; % <--- L is not defined!
end

Since this is supposed to be an error condition, use error to print your error message and exit the function abnormally. You then also don't need else:

function LLT = Cholesky(A)
    if is_pos_def(A) == 1
        error('Cholesky Decomposition is used only for positive definite matrixes')
    end
    % ...computations...
    LLT = L;
end

Note we don't have the definition of is_pos_def(), so we can't discuss it. However, I'd expect it to return true if the matrix is positive definite, in which case you do want to run your computations, so I'd expect a test

    if ~is_pos_def(A)
        error('Cholesky Decomposition is used only for positive definite matrixes')
    end

This reads more meaningfully: "if not is positive definite..."

The spacing around operators is inconsistent, as is the use of empty lines. In particular, I see L(i,1)=A(i,1)..., L(i,i)= sqrt(... and L(j,i) = (A(i,j)... (no spaces around the assignment operator, a space after but not before, and a space on either side). Pick one style and stick to it. I personally prefer spaces around assignment operator, as it helps me read the code.

The assignment at the end, LLT = L, is necessary for the function to produce an output value. But I find it really ugly because the only purpose of this assignment is to change the name of a variable. I would suggest one of two alternatives:

1. declare your function as function L = Cholesky(A), then the variable L is the output variable; or
2. use LLT directly instead of L everywhere in your code. Both solutions yield the same result: the ugly assignment is no longer needed.

The code

n = rank(A);
L = zeros(rank(A));

is awkward as well. Why compute the rank of A twice? The second line should read L = zeros(n);.

The loop

for i=2:n
    L(i,1)=A(i,1)/L(1,1);
end

can be trivially vectorized. Here you apply a division with the same number to a series of values, which can be written as

L(2:n,1) = A(2:n,1)/L(1,1);

Vectorizing the loop above makes the code more compact and easier to read (IMO), which is the main reason to vectorize code. It will also be a bit faster, but the difference has shrunk a lot in recent years. MATLAB first included a JIT (Just In Time) compiler in 2006 I think, and it has been steadily improving over the years. In MATLAB R2015b, they introduced a completely new JIT, and since then the difference between trivial loops and vectorized code is no longer important, and a more complex loop is oftentimes faster than the vectorized version (usually the case when large intermediate arrays are necessary to vectorize).

Finally, vectorizing the main inner loop,

for j = i+1:n
    L(j,i) = (A(i,j)-(dot(L(j,1:i-1), L(i,1:i-1))))/L(i,i);
end

is a bit more complex. Let's do this step by step. dot(L(j,1:i-1), L(i,1:i-1)) is the same as L(j,1:i-1) * L(i,1:i-1).', except the latter is a matrix product and so will generalize to be applied to many vectors at once (i.e. a range of j values). If we fill in j = i+1:n in that matrix multiplication expression, the result is a column vector. A(i,j) for that range of j values is a row vector, we will have to transpose it to make it match the other result. Finally, we arrive at:

L(i+1:n,i) = (A(i,i+1:n).' - (L(i+1:n,1:i-1) * L(i,1:i-1).')) / L(i,i);

I do think that this expression is harder to read than the loop version, and I might prefer to keep the loop version at least in a comment to aid understanding. The vectorized code here likely is a bit faster, at least for smaller arrays.

The final function is:

function L = Cholesky(A)
n = rank(A);
L = zeros(n);
L(1,1) = sqrt(A(1,1));
L(2:n,1) = A(2:n,1)/L(1,1);
for i = 2:n
   L(i,i) = sqrt(A(i,i)- sum(power(L(i,1:i-1),2)));
   L(i+1:n,i) = (A(i,i+1:n).' - (L(i+1:n,1:i-1) * L(i,1:i-1).')) / L(i,i);
end
end