We can exploit the structure of a real, positive definite, symmetric matrix by using the Cholesky decomposition to compute the inverse. The standard MATLAB inv function uses LU decomposition which requires twice as many operations as the Cholesky decomposition and is less accurate. Rarely does one need to compute the inverse of a matrix (e.g. when solving a linear system, we should use \), but when it is needed (e.g. least squares or Kalman Filtering applications), the matrix is positive definite and symmetric. Thus, we can speed up and improve the accuracy of the computation by using the Cholesky decomposition. Because the Cholesky decomposition takes half as many operations, the algorithm will decrease run-time by about 50% for large matrices. Accuracy is also improved with some examples showing orders of magnitude improvement (see Felix Govaers's comment).
Eric Blake (2021). Fast and Accurate Symmetric Positive Definite Matrix Inverse Using Cholesky Decomposition (https://www.mathworks.com/matlabcentral/fileexchange/34511-fast-and-accurate-symmetric-positive-definite-matrix-inverse-using-cholesky-decomposition), MATLAB Central File Exchange. Retrieved .
Short addition: The method I mentioned below is only faster for small matrices ~< 400x400. For larger matrices, the invChol_mex function is faster.
Guys, you should consider using
opts.POSDEF = true;
opts.SYM = true;
I = eye(dim_of_A);
Ainv = linsolve(A, I, opts);
Works another 20% faster than this one for me!
I have to fast-invert this matrix Z= (R+jwL) where R is a real diagonal matrix and L is a real, symmetric, positive definite matrix. Is there any way to use this kind of inversion for the Z matrix?
Thank you in advance!
I wonder if even more speed could be achieved by using intel mkl?
works great thanks a lot! Is there any way to also get the decomposition itself from it? (this is needed in maximum likelihood applications to calculate the log determinant)
Thanks for the code.
I would like to know what is the reference of this code? I mean the technical details of the mentioned fast and accurate computation. Unfortunately it is not mentioned in the code. In case someone is interested to cite it, what reference (paper/text book) should be cited?
I am curious to know the complexity analysis of this inversion method.
Thanks a lot for the code, very much appreciated. It should be noted that the Cholesky inversion is also of higher precision. This made a huge difference in my application!
err_chol_inv = A_cholinv * A - eye(matdim);
err_inv = A_inv * A - eye(matdim);
Max diff chol inv: 5.704351e-09
Max diff inv: 1.980871e-04
This can also be achieved by using Matlab code only (without the speed up advantage):
L = inv(chol(A));
A_inv = L * L';
Thanks for the updated code, Eric!
Works like a charm.
On a side note, and I think it is worth mentioning in the description that the speed-up is noticeable when the matrices are large (in the current example they are small). In my case, I have to solve a big least-squares problem - and invChol, literally, saves the day.
When I run the test.m file i receive an illegal value warning in the dpotrf procedure. Any Help would be greatly appreciated
Excellent code, simple yet efficient, thank you!
One comment. When I compile this in Linux (x86) on MATLAB 2009, I'm getting these warnings:
/opt/matlab/bin/mex invChol_mex.c -lmwlapack
invChol_mex.c: In function ‘mexFunction’:
invChol_mex.c:75: warning: implicit declaration of function ‘dpotrf_’
invChol_mex.c:87: warning: implicit declaration of function ‘dpotri_’
invChol_mex.c:111: warning: implicit declaration of function ‘spotrf_’
invChol_mex.c:123: warning: implicit declaration of function ‘spotri_’
Not critical, of course, but still.
Find the treasures in MATLAB Central and discover how the community can help you!Start Hunting!