I was wondering what is your recommended way to compute the inverse of a matrix?
The ways I found seem not satisfactory. For example,
> c=rbind(c(1, -1/4), c(-1/4, 1)) > c [,1] [,2] [1,] 1.00 -0.25 [2,] -0.25 1.00 > inv(c) Error: could not find function "inv" > solve(c) [,1] [,2] [1,] 1.0666667 0.2666667 [2,] 0.2666667 1.0666667 > solve(c)*c [,1] [,2] [1,] 1.06666667 -0.06666667 [2,] -0.06666667 1.06666667 > qr.solve(c)*c [,1] [,2] [1,] 1.06666667 -0.06666667 [2,] -0.06666667 1.06666667 Thanks!
15 Answers
solve(c) does give the correct inverse. The issue with your code is that you are using the wrong operator for matrix multiplication. You should use solve(c) %*% c to invoke matrix multiplication in R.
R performs element by element multiplication when you invoke solve(c) * c.
You can use the function ginv() (Moore-Penrose generalized inverse) in the MASS package
3Note that if you care about speed and do not need to worry about singularities, solve() should be preferred to ginv() because it is much faster, as you can check:
require(MASS) mat <- matrix(rnorm(1e6),nrow=1e3,ncol=1e3) t0 <- proc.time() inv0 <- ginv(mat) proc.time() - t0 t1 <- proc.time() inv1 <- solve(mat) proc.time() - t1 1solve(matrix) = Inverse of the matrix, it does the job perfectly.
1Use solve(matrix) if the matrix is larger than 1820x1820. Using inv() from matlib or ginv() from MASS takes longer or will not solve at all because of RAM limits.