9

What is the best way of calculating the diagonal of the inverse of a symmetric dense matrix (2000 * 2000)? Currently I calculate the inverse first using solve(x) and then extract the diagonal (diag(y)). Even though it works but I'm wondering whether there is a better way to do it so the code runs faster. I tried chol2inv() but it didn't work since my matrix is not positive-definite.

Update: For anyone who may be interested, I was able to speed up the matrix inversion by using an optimized math library Intel MKL. It takes 3 seconds to inverse a 2000 * 2000 matrix on my machine. Intel MKL is available with Microsoft R Open.

Katherine
  • 111
  • 4
  • In the short google search I did I found quite some links to algorithms and other implementations related to this problem. Do you explicitly need each element of the diagonal or do you intend to use it in another calculation? If you want to use it in further calculation, that knowledge would be crucial for speeding up execution. – Vandenman Jul 18 '17 at 20:58
  • I use the diagonal in other calculations. Specifically, I divide each element in another vector by each element of the diagonal of the inverse matrix and sum them up. – Katherine Jul 18 '17 at 22:36
  • How do you know that each element of the diagonal of the inverse is non zero? Does the matrix have some property that guarantees this? – dmuir Jul 24 '17 at 11:25

2 Answers2

1

If you use LU decomposition then it is possible to get a 30-40% speed up over sum(diag(solve(x)) since it is not necessary to construct the full inverse to get its diagonal.

library( Matrix )

N   <- 2e3; 
mat <- matrix( runif( N*N), ncol = N ); 

t       <- Sys.time(); 
decomp  <- expand(lu( mat)); 
diag_lu <- colSums((solve(decomp$L) %*% t(decomp$P)) * t(solve(decomp$U)));
t_lu    <- Sys.time()-t;

t           <- Sys.time(); 
diag_simple <- diag(solve( mat));  
t_simple    <- Sys.time()-t;  

cat( all.equal( diag_simple, diag_lu), t_lu, t_simple )

Note the matrices generated by expand(lu( mat)) are triangular or permutation matrices so the subsequent matrix function use faster algorithms exploiting this.

Rob
  • 169
  • 7
0

If your matrix has no nice properties like being symmetric, diagonal, or positive-definite, your only choice sadly is to do sum(diag(solve(x)))

How long does that take to run on your matrix?

Alex Braksator
  • 308
  • 1
  • 11
  • It's symmetric but unfortunately not positive-definite. Chol2inv would be a lot faster. It took about 15 seconds to run. The runtime is not terrible but it adds up quickly since I have an iterative process so the matrix is inversed every time in the process. – Katherine Jul 18 '17 at 22:12