Clear Filters
Clear Filters

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

55 views (last 30 days)
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.

Accepted Answer

John D'Errico
John D'Errico on 15 Nov 2015
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.
  2 Comments
Luke Skywalker
Luke Skywalker on 15 Nov 2015
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
John D'Errico on 15 Nov 2015
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.

Sign in to comment.

More Answers (2)

Marc
Marc on 15 Nov 2015
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....
  3 Comments
John D'Errico
John D'Errico on 15 Nov 2015
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
Marc on 15 Nov 2015
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.

Sign in to comment.


Greg Bishop
Greg Bishop on 15 Feb 2017
Your code doesn't seem to provide the same answer as chol(A) does.

Categories

Find more on Linear Algebra in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!