I compared various methods to compute the inverse of a symmetric matrix:
- solve (from the package LAPCK)
- solve (but using a higher machine precision)
- qr.solve (said to be faster)
- ginv (MASS package, implementation of the Moore-Penrose algo)
- chol2inv (using the Cholesky decomposition)
The inverse-matrix was compared through their eigenvalues:
R
library(MASS)
## Create the matrix
m = replicate(10, runif(n=10))
m[lower.tri(m)] = t(m)[lower.tri(m)]
## Inverse the matrix
inv1 = solve(m)
inv2 = solve(m, tol = .Machine$double.eps)
inv3 = qr.solve(m)
inv4 = ginv(m)
inv5 = chol2inv(m)
## Eigenvalues of the inverse
em1=eigen(inv1)
em2=eigen(inv2)
em3=eigen(inv3)
em4=eigen(inv4)
em5=eigen(inv5)
## Plot the abs of the eigenvalues (may be complex)
myPch=c( 20, 15, 17, 25, 3 )
plot(abs(em1$values),pch=myPch[1],cex=1.5)
points(abs(em2$values),pch=myPch[2], cex=1.5)
points(abs(em3$values),pch=myPch[3], cex=1.5)
points(abs(em4$values),pch=myPch[4], cex=1.5)
points(abs(em5$values),pch=myPch[5], cex=1.5)
legend( "topright", c("solve","solve-double","solve-fast","Moore-Penrose","Cholesky"), pch=myPch )
As you can see the inverse given by the Cholesky method is clearly different from the other.
According to this post, if the matrix is symmetric (in our case yes), the Cholesky method is to be preferred: Matrix inversion or Cholesky?
but solve() being the "official-wellspread" R method to invert method, I may rather misunderstand something...
Any good tip?
Thanks in advance,
