[R] Building an array from matrix blocks

[Eduardo de Oliveira Horta]

Thanks!

On Mon, Feb 21, 2011 at 11:40 AM,  <rex.dwyer at syngenta.com> wrote:
> Well, you can lose B by just adding to X in the first for-loop, can't you?
> For (...) X <- X + A[...]
>
> But if you want elegance, you could try:
>
> X = Reduce("+",lapply(1:(p+1), function(i) A[i:(n-p-1+i),i:(n-p-1+i)]))
>
> I imagine someone can be even more eleganter than this.
>
> rad
>
> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf Of Eduardo de Oliveira Horta
> Sent: Saturday, February 19, 2011 9:49 AM
> To: r-help
> Subject: [R] Building an array from matrix blocks
>
> Hello,
>
> I've googled for a while and couldn't find anything on this topic: say
> I have a matrix A and want to build matrices B1, B2,... using blocks
> from A (or equivalently an array B with B[,,i] being a block from A),
> and that I must sum the B[,,i]'s.
>
> I've come up with this rather non-elegant code:
>
>> n = 6
>> p = 3
>>
>> A <- matrix(1:(n^2), n, n, byrow=TRUE)
>>
>> B <- array(0, c(n-p, n-p, p+1))
>> for (i in 1:(p+1)) B[,,i] <- A[i:(n-p-1+i), i:(n-p-1+i)]
>>
>> X <- matrix(0, n-p, n-p)
>> for (i in 1:(p+1)) X <- X + B[,,i]
>> A
>     [,1] [,2] [,3] [,4] [,5] [,6]
> [1,]    1    2    3    4    5    6
> [2,]    7    8    9   10   11   12
> [3,]   13   14   15   16   17   18
> [4,]   19   20   21   22   23   24
> [5,]   25   26   27   28   29   30
> [6,]   31   32   33   34   35   36
>> B
> , , 1
>
>     [,1] [,2] [,3]
> [1,]    1    2    3
> [2,]    7    8    9
> [3,]   13   14   15
>
> , , 2
>
>     [,1] [,2] [,3]
> [1,]    8    9   10
> [2,]   14   15   16
> [3,]   20   21   22
>
> , , 3
>
>     [,1] [,2] [,3]
> [1,]   15   16   17
> [2,]   21   22   23
> [3,]   27   28   29
>
> , , 4
>
>     [,1] [,2] [,3]
> [1,]   22   23   24
> [2,]   28   29   30
> [3,]   34   35   36
>
>> X
>     [,1] [,2] [,3]
> [1,]   46   50   54
> [2,]   70   74   78
> [3,]   94   98  102
>
> Note that the blocks B[,,i] are obtained by sweeping the diagonal of
> A. I wonder if there is a better and faster way to achieve this using
> block matrix operations for instance. Actually what matters most for
> me is getting to the matrix X, so if it is possible to do this without
> having to construct the array B it would be ok as well...
>
> Interesting observation:
>
>> system.time(for (j in 1:10000) {X <- matrix(0, n-p, n-p); for (i in 1:(p+1)) X <- X + B[,,i]})
>   user  system elapsed
>   0.27    0.00    0.26
>> system.time(for (j in 1:10000) {X <- apply(B,c(1,2),sum)})
>   user  system elapsed
>   1.82    0.02    1.86
>
> Thanks in advance, and best regards,
>
> Eduardo Horta
>
>> sessionInfo()
> R version 2.11.1 (2010-05-31)
> x86_64-pc-mingw32
>
> locale:
> [1] LC_COLLATE=Portuguese_Brazil.1252  LC_CTYPE=Portuguese_Brazil.1252
> [3] LC_MONETARY=Portuguese_Brazil.1252 LC_NUMERIC=C
> [5] LC_TIME=Portuguese_Brazil.1252
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
> [1] Revobase_4.2.0   RevoScaleR_1.1-1 lattice_0.19-13
>
> loaded via a namespace (and not attached):
> [1] grid_2.11.1       pkgXMLBuilder_1.0 revoIpe_1.0       tools_2.11.1
> [5] XML_3.1-0
>
> ______________________________________________
> 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.
>
>
>
>
> message may contain confidential information. If you are not the designated recipient, please notify the sender immediately, and delete the original and any copies. Any use of the message by you is prohibited.
>
>