Why is the built-in Cholesky function so much faster than my own implementation?

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.

Answer 1

[John D'Errico]

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.

[Luke Skywalker]

I was indeed surprised, as at least the incremental Cholesky, where just one additional row needs to be calculated, could have potentially significantly faster.

Do you se a way to increase speed on my calculations? Is it possible to deploy the blas-tool?

Best

[John D'Errico]

You simply won't come near to the built-in chol using MATLAB code. Accept that as a fact.

In order to do better, you might try writing C code directly, rather than compiling MATLAB code. Then you could insert calls to the blas. Since I cannot even spell C, I can't/won't offer advice in this respect.

Answer 2

[Marc]

I am just guessing and based on the help information for chol() is that it uses sparseness as you increase the size of the matrix.

Several years ago, I ran into a fsolve() problem where an in house fortran based solver wasn't flinching as I increased the size of the matrix but fsolve() after n = 100,000 equations really got bogged down. One of TMW application engineers looked at my problem, tweaked things to take advantage of sparseness and it flew within milliseconds that I could not tell the difference between ML and the in house program.

Some of these guys are ninjas with this stuff....

[Luke Skywalker]

I didn't mention it above, but the calculations were based on random full symmetric positive-definite matrices. So I guess that at least in this case sparseness can not be used?

[John D'Errico]

sparsity information is ONLY used when you specify your matrix as sparse when you create the matrix. MATLAB does NOT decide on the fly that your matrix can benefit because it appears sparse. You need to specify that. As well, if you convert your matrix to sparse form and it is not truly sparse (i.e., a sparse matrix should have MANY zeros) then code will often take longer to run. So, just blithely using sparse form on all of your matrices will probably result in slower code overall.

[Marc]

You guys may be right. I have no idea. I did not look into what you wrote but in the help file it says that chol() uses only the diagonal and upper triangle, assuming the lower is a complex conjugate transpose of the upper.

Does this mean it is doing "half" the computations that you are doing?

In general, I never disagree with John. His answer above wasn't showing when I typed mine in but I would always defer to him over me.

Answer 3

[Greg Bishop]

Your code doesn't seem to provide the same answer as chol(A) does.

[John D'Errico]

This is not an answer. Please don't answer a question, just to add a comment. Use a comment for that.