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.
$\begingroup$
$\endgroup$
8
2 Answers
$\begingroup$
$\endgroup$
4
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.
-
$\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$– jitendraCommented Aug 31, 2011 at 19:37 -
11$\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$– RaxelCommented Mar 6, 2015 at 7:11
-
-
$\begingroup$ @Ga13 Not necessarily (and not on two systems that I tested this on). This might depend on your CPU and which linear algebra routines you are using. $\endgroup$– WolfgangCommented Mar 25, 2020 at 11:26
$\begingroup$
$\endgroup$
5
If you are working with covariance matrix or any positive definite matrix you can use pd.solve is faster.
Following the Wolfgang example:
library(MASS)
library(mnormt)
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))
> system.time(RI1 <- solve(R))
usuário sistema decorrido
13.21 0.03 13.76
> system.time(RI2 <- chol2inv(chol(R)))
usuário sistema decorrido
5.62 0.05 5.80
> system.time(RI3 <- qr.solve(R))
usuário sistema decorrido
20.42 0.09 21.10
> system.time(RI4 <- pd.solve(R))
usuário sistema decorrido
5.53 0.00 5.61
-
-
-
$\begingroup$ This might depend on your CPU and which linear algebra routines you are using. On two systems that I tested this on (both using OpenBLAS),
chol2inv(chol(R)))was slightly faster thanpd.solve(). $\endgroup$– WolfgangCommented Mar 25, 2020 at 11:21 -
1$\begingroup$ Also, it would be better to do this more than once (and average), because single results can be skewed by some other process taking up CPU cycles in the background. Using package
microbenchmarkwould be even better. $\endgroup$– WolfgangCommented Mar 25, 2020 at 11:25 -
3$\begingroup$ pd.solve is just a wrapper for chol2inv(chol(R))), so it can't be faster (look at the code). $\endgroup$– Car LozCommented Oct 6, 2020 at 18:10
solvemethod 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$solve? Of course, absent special structure, you can't escape the theoretical complexity bounds on general matrix inversion. $\endgroup$