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

[Paul Rathouz]

Thank you for your response, Oliver.  

I think that the Matrix package has its own solve() functions that take advantage of the structure.  If you make a block diagonal matrix, say, with kronecker() and do not impose any special structure on it.  then, calling solve() will take very much longer than when using one of the Matrix-package classes of matrices. This is what puzzles me. -- pr

On Jun 26, 2012, at 12:56 PM, Paul Rathouz wrote:

> 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 
> 

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