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@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.
