# Cholesky decomposition

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?

• Welcome to Code Review! I changed the title so that it describes what the code does per site goals: "State what your code does in your title, not your main concerns about it.". Please check that I haven't misrepresented your code, and correct it if I have. Nov 18, 2021 at 9:01
• Doesn't MATLAB have a built-in function for the Cholesky decomposition? Nov 18, 2021 at 9:23
• @MartinR how can I access the source code of the "Chol" function in MatLab? Nov 18, 2021 at 10:38
• What @MartinR is asking is, is there a reason why you wrote this Cholesky function rather than just calling chol? Nov 18, 2021 at 12:03
• @MartinR This is part of my academic project, and I just want to do it efficiently on MatLab by myself Nov 18, 2021 at 13:19

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