2

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,

Community
  • 1
  • 1
Xavier Prudent
  • 1,570
  • 3
  • 25
  • 54
  • 1
    Your example isn't reproducible because you didn't set a random seed. However, check if the matrix is positive definite. Also, general advice is to avoid inverting a matrix. – Roland Nov 03 '14 at 16:31
  • 3
    Also, it should of course be `inv5 = chol2inv(chol(m))`. – Roland Nov 03 '14 at 16:34
  • 1
    @Roland Nitpicking; general advice is to avoid *explicitly* inverting a matrix. Using a decomposition is often a way to avoid doing the explicit inversion of the original matrix. – Gavin Simpson Nov 03 '14 at 16:34
  • 1
    @Roland & if one does `chol2inv(chol(m))` as you suggest, for the random example I just ran, R states: `Error in chol.default(m) : the leading minor of order 3 is not positive definite`, thus highlighting the problem the OP faces. +1 – Gavin Simpson Nov 03 '14 at 16:38
  • A poor-man way to improve my code so that the matrix be positive definite and symmetric: library(caper); t = rtree(n=10); m = vcv( t ) – Xavier Prudent Nov 04 '14 at 09:08

1 Answers1

5

You need to pass the Cholesky decomposition to chol2inv:

inv5 = chol2inv(chol(m))

If m is positive definite (which it probably isn't for your not-reproducible input) this should give the same result as the other methods.

Roland
  • 127,288
  • 10
  • 191
  • 288