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

[Oliver Ruebenacker]

     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)