2

I want to inverse a square symmetric positive definite matrix. I know there are two functions solve() and chol2inv() in R but their results is different. I need to know why this happen?

Thank you.

Richard Chambers
  • 15,607
  • 3
  • 71
  • 101
Mohammad
  • 497
  • 1
  • 7
  • 18
  • Please make your situation reproducible, i.e. provide us with the data and the code needed to mimic your situation. See http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example for more tips on how to do this. – Paul Hiemstra Mar 11 '13 at 10:22
  • how different are the biggest and smallest eigen values? i.e. for a matrix A, `range(eigen(A)$values)` the bigger the range, the harder it is to invert "accurately". – Sean Mar 11 '13 at 10:42

2 Answers2

6

Here are several ways to compute matrix inverse, including solve() and chol2inv():

> A <- matrix(c(2, -1, 0, -1, 2, -1, 0, -1, 2), 3)

> solve(A)
     [,1] [,2] [,3]
[1,] 0.75  0.5 0.25
[2,] 0.50  1.0 0.50
[3,] 0.25  0.5 0.75

> chol2inv(chol(A))
     [,1] [,2] [,3]
[1,] 0.75  0.5 0.25
[2,] 0.50  1.0 0.50
[3,] 0.25  0.5 0.75

> library(MASS)
> ginv(A)
     [,1] [,2] [,3]
[1,] 0.75  0.5 0.25
[2,] 0.50  1.0 0.50
[3,] 0.25  0.5 0.75
NPE
  • 464,258
  • 100
  • 912
  • 987
5

For solve you need to give your original matrix, but for chol2inv you use precomputed cholesky decomposition:

set.seed(1)
a<-crossprod(matrix(rnorm(9),3,3))
a_chol<-chol(a)
solve(a)
            [,1]        [,2]       [,3]
[1,]  1.34638151 -0.02957435  0.8010735
[2,] -0.02957435  0.32780020 -0.1786295
[3,]  0.80107345 -0.17862950  1.4533671
chol2inv(a_chol)
            [,1]        [,2]       [,3]
[1,]  1.34638151 -0.02957435  0.8010735
[2,] -0.02957435  0.32780020 -0.1786295
[3,]  0.80107345 -0.17862950  1.4533671
Jouni Helske
  • 6,337
  • 28
  • 52