[R] Building an array from matrix blocks

[rex.dwyer at syngenta.com]

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.