3

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
  • 16,643
  • 4
  • 81
  • 106
Mohammad
  • 571
  • 1
  • 8
  • 24
  • 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

7

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
  • 486,780
  • 108
  • 951
  • 1,012
6

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,427
  • 29
  • 52