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.