Thanks!
On Mon, Feb 21, 2011 at 11:40 AM, <rex.dw...@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-boun...@r-project.org [mailto:r-help-boun...@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@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.
>
>
______________________________________________
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.
