21
$\begingroup$

I need to calculate matrix inverse and have been using solve function. While it works well on small matrices, solve tends to be very slow on large matrices. I was wondering if there is any other function or combination of functions (through SVD, QR, LU, or other decomposition functions) that can give me faster results.

$\endgroup$
  • 2
    $\begingroup$ Can you provide more information? What are the approximate dimensions? Does the matrix have any special structure (symmetry, sparsity, etc.)? What is your quantitative definition of "slow"? And "fast"? $\endgroup$ – cardinal Aug 29 '11 at 18:19
  • $\begingroup$ The approximate dimensions are like 2000x2000. The matrix doesn't have any special structure. Well, solve method definitely does my work but I want the algorithm to be faster. So, I am just wondering if there is a more efficient (in time context) function for calculating inverse for such large size matrix. $\endgroup$ – jitendra Aug 29 '11 at 18:27
  • 1
    $\begingroup$ Have you tried any of the other suggestions on the help page for solve? Of course, absent special structure, you can't escape the theoretical complexity bounds on general matrix inversion. $\endgroup$ – cardinal Aug 29 '11 at 18:38
  • 3
    $\begingroup$ @Cardinal The trick is to probe further concerning the actual application, for as you know, in many cases inverting the matrix is unnecessary (and time-consuming and error-prone). $\endgroup$ – whuber Aug 29 '11 at 18:44
  • $\begingroup$ @whuber: This is a very good point. I suppose sometimes I approach these questions a little too directly. $\endgroup$ – cardinal Aug 29 '11 at 18:54
23
$\begingroup$

Have you tried what cardinal suggested and explored some of the alternative methods for computing the inverse? Let's consider a specific example:

library(MASS)

k   <- 2000
rho <- .3

S       <- matrix(rep(rho, k*k), nrow=k)
diag(S) <- 1

dat <- mvrnorm(10000, mu=rep(0,k), Sigma=S) ### be patient!

R <- cor(dat)

system.time(RI1 <- solve(R))
system.time(RI2 <- chol2inv(chol(R)))
system.time(RI3 <- qr.solve(R))

all.equal(RI1, RI2)
all.equal(RI1, RI3)

So, this is an example of a $2000 \times 2000$ correlation matrix for which we want the inverse. On my laptop (Core-i5 2.50Ghz), solve takes 8-9 seconds, chol2inv(chol()) takes a bit over 4 seconds, and qr.solve() takes 17-18 seconds (multiple runs of the code are suggested to get stable results).

So the inverse via the Choleski decomposition is about twice as fast as solve. There may of course be even faster ways of doing that. I just explored some of the most obvious ones here. And as already mentioned in the comments, if the matrix has a special structure, then this probably can be exploited for more speed.

$\endgroup$
  • $\begingroup$ Thanks a lot for this solution. I, at least, know one method that can solve it half the time as compared to solve :-) $\endgroup$ – jitendra Aug 31 '11 at 19:37
  • 8
    $\begingroup$ The Cholesky decomposition is a good choice for covariance/correlation matrices but keep in mind that in general the matrix has to be Hermitian(in case of real matrices that means symmetric), positive definite matrix. That uses half of the memory required for LU decomposition. $\endgroup$ – Raxel Mar 6 '15 at 7:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.