I am currently using a library (gnu gsl) which takes great pains in both API design and documentation to guide people away from using matrix inverses explicitly and towards using the LU decomposition. I have code which is going to calculate Mahalanobis distance, requiring the inverse of covariance matrix. Is there any reason to use a LU decomposition in calculating Mahalanobis distance? At present I'm a little bit skeptical, as I will invert a (say 100x100) covariance matrix once, then calculate *many* M distances. While I don't know exactly how to formulate Mahalanobis distance using the LU decomposition, or even know if it's possible (googling on both yields nothing that obviously answers my question), I would presume that any rewrite would require additional computation and therefore time. Would there be a significant improvement in numerical accuracy? At present my default inclination is to find the actual inverse and then calculate distance by the standard formula.
Also, I have to solve some systems of linear equations. But the gsl libray includes library routines to do this given a LU decomposition, so unless corrected I'm assuming it's a no-brainer to solve them using a LU decomposition.