[R] Building an array from matrix blocks

[Eduardo de Oliveira Horta]

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