0

Let us say I have a very large correlation matrix of this form:

t1.rep1 = rnorm(n=100,mean=10,sd=)
t2.rep1 = t1.rep1 + rnorm(n=100,mean=3,sd=2)
t3.rep1 = t1.rep1 + rnorm(n=100,mean=2,sd=2)
t1.rep2 = rnorm(n=100,mean=2,sd=1)
t2.rep2 = t1.rep2 + rnorm(n=100,mean=0.5,sd=0.5)
t3.rep2 = t1.rep2 + rnorm(n=100,mean=0.7,sd=0.9)
t1.rep3 = rnorm(n=100,mean=2,sd=1)
t2.rep3 = t1.rep3 + rnorm(n=100,mean=0.5,sd=0.5)
t3.rep3 = t1.rep3 + rnorm(n=100,mean=0.7,sd=0.9)
Sigma = matrix(
  c(cov(t1.rep1, t1.rep1), 0, 0, cov(t1.rep1, t2.rep1), 0, 0, cov(t1.rep1, t3.rep1), 0, 0,
  0, cov(t1.rep2, t1.rep2), 0, 0, cov(t1.rep2, t2.rep2), 0, 0, cov(t1.rep2, t3.rep2), 0,
  0, 0, cov(t1.rep3, t1.rep3), 0, 0, cov(t1.rep3, t2.rep3), 0, 0, cov(t1.rep3, t3.rep3),
  cov(t2.rep1, t1.rep1), 0, 0, cov(t2.rep1, t2.rep1), 0, 0, cov(t2.rep1, t3.rep1), 0, 0,
  0, cov(t2.rep2, t1.rep2), 0, 0, cov(t2.rep2, t2.rep2), 0, 0, cov(t2.rep2, t3.rep2), 0,
  0, 0, cov(t2.rep3, t1.rep3), 0, 0, cov(t2.rep3, t2.rep3), 0, 0, cov(t2.rep3, t3.rep3),
  cov(t3.rep1, t1.rep1), 0, 0, cov(t3.rep1, t2.rep1), 0, 0, cov(t3.rep1, t3.rep1), 0, 0,
  0, cov(t3.rep2, t1.rep2), 0, 0, cov(t3.rep2, t2.rep2), 0, 0, cov(t3.rep2, t3.rep2), 0,
  0, 0, cov(t3.rep3, t1.rep3), 0, 0, cov(t3.rep3, t2.rep3), 0, 0, cov(t3.rep3, t3.rep3)),
  nrow = 9, ncol = 9)

My correlation matrix is Sigma.

And I want to compute its inverse, i.e., 

Sigma.inv = solve(Sigma)

In reality my sigma is much much bigger and taking its inverse it's taking a long time.

Is there anyway to use the sparsity and structure of the matrix to compute the inverse of Sigma in a much faster/efficient way?

I need to compute this inverse iteratively so a fast way of computing the inverse would help my algorithm speed considerably.

Dnaiel
  • 7,622
  • 23
  • 67
  • 126

1 Answers1

1

The Sigma you have provided is actually block-diagonal:

x = c(1,4,7,2,5,8,3,6,9)
Sigma[x,x]

          [,1]     [,2]      [,3]     [,4]     [,5]     [,6] ...
[1,] 0.9494388 1.130673 0.9825316 0.000000 0.000000 0.000000 ...
[2,] 1.1306727 4.983144 1.2112634 0.000000 0.000000 0.000000 ...
[3,] 0.9825316 1.211263 5.0771423 0.000000 0.000000 0.000000 ...
[4,] 0.0000000 0.000000 0.0000000 1.211892 1.223293 1.328587 ...
[5,] 0.0000000 0.000000 0.0000000 1.223293 1.469146 1.242400 ...
[6,] 0.0000000 0.000000 0.0000000 1.328587 1.242400 2.377406 ...
...

And block-diagonal matrices can be inverted much faster block-wise. Just replace each block with its inverse.

Andrey Shabalin
  • 4,389
  • 1
  • 19
  • 18
  • @Shabalin, thanks so much, looks like a great suggestion. Any trick to do the inverse of each block fast in R? Just do an apply function to each block? or do a for loop? for loops tend to be quite slow in R... – Dnaiel Feb 25 '14 at 18:09
  • Maybe you can store the matrix in blocks to begin with. – Andrey Shabalin Feb 25 '14 at 23:54
  • Also, take a look at http://stats.stackexchange.com/questions/14951/efficient-calculation-of-matrix-inverse-in-r – Andrey Shabalin Feb 25 '14 at 23:55
  • @Shabalin Thanks a lot, all pointers into the right direction. However, the block by block inversion is still quite slow. I tried cholesky and it was not very good. Is there a way to speed the inversion of a matrix further, is multithreading a good take? I'll post this as a separate question. – Dnaiel Feb 26 '14 at 17:40