CODEWITHSUNDEEP
r
user1690124
user1690124

Reputation: 303

Efficient way to calculate array multiplication

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])

Upvotes: 0

Views: 84

Answers (2)

IRTFM
IRTFM

Reputation: 263342

  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.

Upvotes: 3

user2087984
user2087984

Reputation: 1246

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.

Upvotes: 2

Related Questions