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

This might answer your question -- http://stackoverflow.com/questions/6192848/how-to-generalize-outer-to-n-dimensions/6193836#6193836 — Prasad Chalasani, Jan 06 '12 at 22:54
The accepted answer is quite slow; it seems to be slower than repeated calls to outer (unsurprising given how general it is). But I think maybe the second can be adapted... — MMM, Jan 07 '12 at 00:44
I would be very surprised if the `expand.grid` solution were faster than the `outer` solution. — IRTFM, Jan 07 '12 at 02:55
@DWin you'd be right to be surprised; I wrote too soon. Repeated applications of outer() have it at the moment - it's about twice as fast in some limited testing. Oh well. — MMM, Jan 07 '12 at 05:27
Generally matrix operations (kronecker and outer) will outperform dataframe operations (expand.grid), often by a substantial margin. — IRTFM, Jan 07 '12 at 08:22
Indeed. I thought maybe calling outer(outer(outer... would catch up to me, but it hasn't (not for reasonably sized problems, anyway). I thought perhaps there was a single function, or a tensor library calling compiled code, but apparently not. — MMM, Jan 07 '12 at 15:28

Tommy · Accepted Answer · 2012-01-08T09:49:40.450

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.

Makes sense. I was mostly wondering if there was something in base I was missing, but I guess it really is kind of a specialized thing. Thanks. — MMM, Jan 08 '12 at 22:49
Rather than creating `mouter` you can use `Reduce`: `Reduce("%o%",list(x1,x2,x3))`. Don't think performance would change much though. — James, Jan 09 '12 at 13:06
`/ref:` Estimates of the speedup from an optimised BLAS can be an order of magnitude higher: http://blog.felixriedel.com/2012/10/speed-up-r-by-using-a-different-blas-implementation/ — isomorphismes, Nov 09 '14 at 18:58

MYaseen208 · Answer 2 · 2012-01-06T22:13:58.160

1

You can use tensor package.

And also %o% function

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

edited Jan 06 '12 at 22:13

answered Jan 06 '12 at 22:08

MYaseen208

22,666
37
165
309

Well, `%o%` is exactly `outer(x,y,'*')` so that won't help the speed problem. As always, :-) I have to ask,"what is the problem you're trying to solve?" It's possible there's a different way through your data. – Carl Witthoft Jan 06 '12 at 22:26
I know how to take the outer product of a two matrices (or two vectors); that's not what I need. Looking at the tensor package it isn't immediately clear to me how it would help. – MMM Jan 06 '12 at 22:54
@Carl Witthoft The resultant array is (almost) the joint pmf of a multivariate discrete variable. I need this approximate pmf, or a way to return all its entries greater than some value, and also the sum of the entries. – MMM Jan 06 '12 at 23:05

score 1 · Answer 3 · answered Jan 07 '12 at 01:41

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.

Sorry, I thought it was clear from the question that the output should an array, not a matrix. I've added an example to clarify. All my inputs are vectors, so kronecker and outer are equivalent. Replacing the kronecker function in your answer with another call to outer was the solution I originally had, which does return an array. — MMM, Jan 07 '12 at 02:42
Yeah, it was the reason I mentioned the different arrangement. Your use case was not described so I thought it might satisfy. — IRTFM, Jan 07 '12 at 02:57

Outer/tensor product in R

3 Answers3