MATLAB: Is the built-in Cholesky function so much faster than the own implementation

cholesky

Hi,
I am currently investigating the efficiency of matrix-inversion-methods and came across the Cholesky decomposition. I then implemented the cholesky in Matlab and compared it to the built-in chol()-function.
What you can see in the graph below is a Benchmark of my implemented Cholesky decompositions and the chol()-function: – "Cholesky" is the regular Cholesky Decomposition – "Incremental Cholesky" is a method where an old Cholesky decomp of a Matrix A is used to calculate the decomposition of an incremented Matrix B with one extra row and column
Both functions are also implemented in Mex-C which improved performance a little bit.
Here is the code of my Cholesky implementation:
function L = cholMatlab(M)
n = length( M );
L = zeros( n, n );
for i=1:n
L(i, i) = sqrt(M(i, i) - L(i, :)*L(i, :)');
for j=(i + 1):n
L(j, i) = (M(j, i) - L(i,:)*L(j ,:)')/L(i, i);
end
end
end
Can you tell me why the chol()-function performs so much better?
Best regards
Edit: Should have mentioned this: The calculations are based on random full symmetric positive definite matrices.

Best Answer

And, why are you surprised? chol will have been carefully written using tools like the blas for maximum speed, by professionals who have had many opportunities to learn every trick. As far as your MATLAB code goes, it is basic multiply looped MATLAB code, which is rarely that efficient. Even compiling it will not gain much.