Hi.
I'm looking for an efficient way of writing a function(x,mat1,mat2).
n<-4
m<-4
r<-3
x <- array(sample(1:1000)/10^4,rep.int(n,m))
mat1 <- matrix(sample(1:1000)/10^4,n,n)
mat2 <- matrix(sample(1:1000)/10^4,n,n)
It needs to work for *any* itegers n, r<=m with output (in horrible gory
detail) equivalent to:
ans<-array(0,rep.int(n,(m+1)))
for(i1 in 1:n) for(i2 in 1:n) for(i3 in 1:n) for(i4 in 1:n) for(i5 in 1:n)
ans[i1,i2,i3,i4,i5] <- sum(x[i1,i2,,i3]*mat1[,i4]*mat2[,i5])
Notice how I take a pointwise product of x in the "r^th slot" of x
with the
"1^st slot" of mat1 and mat2, and then sum out. This is the guts of
what I
want to do. It's kind of like an inner product "%*%" but on three
objects
instead of two. Although, the above is what I want to do for my application,
in fact it would be nice to be able to do this for three (or more) arrays in
three (or more) specified slots.
We should also be able do it with outer and apply:
y <- outer(outer(x,mat1),mat2)
len <- length(dim(y))
ans2 <- apply(y,c(1:len)[-c(r,len-3,len-1)],FUN=sum)
But for some reason I can't get this to give the same answer as above even
though dim(ans)=dim(ans2). I think apply is not doing what I think it is
doing. Taking the double outer is bad news efficiency wise anyway.
Any help, ideas, or direction will be greatly appreciated.
Jeremy.
[[alternative HTML version deleted]]
