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
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
