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 <rath...@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@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) ______________________________________________ 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.