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
Hello,
My guess would be that solve() does not take advantage of the
special structure of the matrix and that you may want a sparse matrix
representation.
Take care
Oliver
On Tue, Jun 26, 2012 at 1:56 PM, Paul Rathouz <rathouz at
biostat.wisc.edu> 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
>
> ______________________________________________
> R-help at 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.
--
Oliver Ruebenacker, Bioinformatics and Network Analysis Consultant
President and Founder of Knowomics
(http://www.knowomics.com/wiki/Oliver_Ruebenacker)
Consultant at Predictive Medicine
(http://predmed.com/people/oliverruebenacker.html)
SBPAX: Turning Bio Knowledge into Math Models (http://www.sbpax.org)
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