Outer/tensor product in R

Question

Given p vectors x1,x2,...,xp each of dimension d, what's the best way to compute their tensor/outer/Kruskal product (the p-array X with entries X[i1,i2,..ip] = x1[i1]x2[i2]...xp[ip])? Looping is trivial, but stupid. Using repeated calls to outer works OK, but doesn't seem like the optimal solution (and will get slower as p increases, obviously). Is there a better way?

Edit:

My current best is

array(apply(expand.grid(x1, x2, x3), 1, prod), dim=rep(d, 3))

which at least "feels better"...

Edit 2: In response to @Dwin, here's a complete example

d=3
x1 = 1:d
x2 = 1:d+3
x3 = 1:d+6
array(apply(expand.grid(x1, x2, x3), 1, prod), dim=rep(d, 3))

, , 1

     [,1] [,2] [,3]
[1,]   28   35   42
[2,]   56   70   84
[3,]   84  105  126

, , 2

     [,1] [,2] [,3]
[1,]   32   40   48
[2,]   64   80   96
[3,]   96  120  144

, , 3

     [,1] [,2] [,3]
[1,]   36   45   54
[2,]   72   90  108
[3,]  108  135  162

Answer 1

It will be hard to beat the performance of outer. This ends up doing a matrix multiplication which is done by the BLAS library. Calling outer repeatedly doesn't matter either, since the last call will dominate both speed and memory wise. For example, for vectors of length 100, the last call is at least 100x slower than the previous one...

Your best bet to get the best performance here is to get the best BLAS library for R. The default one isn't very good. On Linux, you can fairly easily configure R to use ATLAS BLAS. On Windows it is harder, but possible. See R for Windows FAQ.

# multiple outer
mouter <- function(x1, ...) { 
    r <- x1
    for(vi in list(...)) r <- outer(r, vi)
    r
}

# Your example
d=3
x1 = 1:d
x2 = 1:d+3
x3 = 1:d+6 
mouter(x1,x2,x3)

# Performance test
x <- runif(1e2)
system.time(mouter(x,x,x))   # 0 secs (less than 10 ms)
system.time(mouter(x,x,x,x)) # 0.5 secs / 0.35 secs (better BLAS)

I replaced my Windows Rblas.dll with the DYNAMIC_ARCH version of GOTO BLAS at this place which improved the time from 0.5 to 0.35 secs as seen above.

Answer 2

You can use tensor package.

And also %o% function

A <- matrix(1:6, 2, 3)
D <- A %o% A

Answer 3

I find myself wondering if the kronecker product is what you want. I cannot tell from your problem description exactly what is desired, but the elements from this on a small set of arguments are the same (albeit in a different arrangement as those in produced by Chalasani solution you were criticizing as slow:

kronecker( outer(LETTERS[1:2], c(3, 4, 5),FUN=paste), letters[6:8] ,FUN=paste)
     [,1]    [,2]    [,3]   
[1,] "A 3 f" "A 4 f" "A 5 f"
[2,] "A 3 g" "A 4 g" "A 5 g"
[3,] "A 3 h" "A 4 h" "A 5 h"
[4,] "B 3 f" "B 4 f" "B 5 f"
[5,] "B 3 g" "B 4 g" "B 5 g"
[6,] "B 3 h" "B 4 h" "B 5 h"

If you want products, then substitute either prod or "*". In any case offering a sample set of vectors and the desired output is a best practice in posing questions.