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