[R] computing time for solve() in Matrix package

[Paul Rathouz]

Hi -- I am wondering why the time to solve (invert) a block diagonal matrix created with bdiag() from the Matrix package does not scale with the number of blocks [all of the same size]. Or, what I am doing wrong. The code / output below creates a series of block diagonal matrices with 4x4 blocks and with 500, 1000, 1500, and 2000 blocks.  Note that I make a copy of the object to have something to write into, which seems to help. The computational time goes up about 10-fold when the number of blocks grows from 500 to 2000.  Is there a different way to encode the block structure? -- pr


Paul Rathouz, PhD
Professor and Chair
Department of Biostatistics & Medical Informatics
University of Wisconsin School of Medicine and Public Health
K6/446 CSC, Box 4675
600 Highland Avenue
Madison, WI 53792-4675
608.263.1706


> library(Matrix)
Loading required package: lattice

Attaching package: ‘Matrix’

The following object(s) are masked from ‘package:base’:

    det

> 
> #### function to construct AR1 correlation matrix
> #### for one individual
> corr.mat = function(alpha,len)  {
+   powers = abs(outer(0:(len-1),0:(len-1),"-"))
+   corr.mat = alpha^powers
+   return(corr.mat)
+ }
> 
> #### Compute AR1 correlation matrix with alpha=.3 and n.obs=4
> R = corr.mat(.3,4)
> 
> #### To optimize time to invert, need an object to write to
> system.time(block.R <- bdiag(rep(list(R),1000)))
   user  system elapsed 
  0.599   0.007   0.610 
> class(block.R)
[1] "dgCMatrix"
attr(,"package")
[1] "Matrix"
> system.time(inv.block <- solve(block.R))
   user  system elapsed 
  1.147   0.266   1.423 
> system.time(inv.block <- solve(block.R))
   user  system elapsed 
  0.694   0.041   0.740 
> 
> #### How does time to solve scale with number of blocks
> block.R <- bdiag(rep(list(R),500))
> inv.block <- block.R
> system.time(inv.block <- solve(block.R))
   user  system elapsed 
  0.177   0.072   0.251 
> 
> block.R <- bdiag(rep(list(R),1000))
> inv.block <- block.R
> system.time(inv.block <- solve(block.R))
   user  system elapsed 
  0.702   0.041   0.745 
> 
> block.R <- bdiag(rep(list(R),1500))
> inv.block <- block.R
> system.time(inv.block <- solve(block.R))
   user  system elapsed 
  2.064   0.745   2.813 
> 
> block.R <- bdiag(rep(list(R),2000))
> inv.block <- block.R
> system.time(inv.block <- solve(block.R))
   user  system elapsed 
  3.182   1.294   4.480