0

Is there any efficient way to calculate 2x2 matrix H without for statement?

    n=10

    a=array(rnorm(n),c(2,1,n))
    b=array(rnorm(n),c(2,1,n))

    H=matrix(0,2,2)
    for(i in 1:n) H=H+a[,,i] %*% t(b[,,i])
user1690124
  • 303
  • 1
  • 2
  • 7
  • You will want to use the function `apply`. [This answer](http://stackoverflow.com/a/7141669/2752888) might help you. – ZNK Dec 17 '13 at 04:20

2 Answers2

3
  H=matrix(0,2,2)
     for(i in 1:n) H=H+a[,,i] %*% t(b[,,i])
 H
#----------
          [,1]       [,2]
[1,] 10.770929 -0.4245556
[2,] -5.613436 -1.7588095


 H2 <-a[ ,1, ] %*% t(b[ ,1, ])
H2
#------------- 
          [,1]       [,2]
[1,] 10.770929 -0.4245556
[2,] -5.613436 -1.7588095

This does depend on the arrays in question having one of their dimensions == 1, and on the fact that "[" will drop length-1 dimensions unless you specify drop=FALSE.

IRTFM
  • 258,963
  • 21
  • 364
  • 487
  • Because there is no information to loose. Look at `length(a)` and `length(a[,1,])` – IRTFM Dec 17 '13 at 05:13
  • Sorry, I was editing my comment when you responded. It would be good to add this to the answer, that each sub-matrix only has column 1, so this gives a reordering of the data which produces correct results (in general). – Matthew Lundberg Dec 17 '13 at 05:14
2

This is the same (up to FAQ 7.31 issues) as what you calculate:

In case the second dimension truly has only 1 level, you can use

tcrossprod( matrix(a,nr=2), matrix(b,nr=2) )

and more generally,

crossprod( matrix( aperm(a, c(3,1,2)), nc=2), matrix( aperm(b, c(3,1,2)), nc=2) )

If you can create 'a' and 'b' ordered so that you do not need the aperm() it will be still faster.

The relative speed of different solutions depends on the dimensions. If the first two are both big and the last one small, a loop like yours (but using crossprod) might be as quick as you can get.

user2087984
  • 1,236
  • 1
  • 7
  • 3