r - R中的外部/张量积

给定每个维度p的x1,x2,...,xp向量d，计算其张量/外部/ Kruskal乘积的最佳方法是什么（具有条目p的X[i1,i2,..ip] = x1[i1]x2[i2]...xp[ip])数组X？循环很简单，但是很愚蠢。重复调用outer可以，但是似乎不是最佳解决方案（显然，随着p的增加，它会变慢），是否有更好的方法？

编辑：

我目前最好的是

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

至少“感觉更好” ...

编辑2：回应@Dwin，这是一个完整的示例

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

最佳答案

很难击败outer的表现。最终完成了由BLAS库完成的矩阵乘法。重复调用outer也不重要，因为最后一次调用将同时控制速度和内存。例如，对于长度为100的向量，最后一次调用比上一次调用至少慢100倍...

在这里获得最佳性能的最佳选择是获得R的最佳BLAS库。默认库不是很好。在Linux上，您可以相当容易地将R配置为使用ATLAS BLAS。在Windows上很难，但是可以。请参见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)

我在this place用GOTO BLAS的DYNAMIC_ARCH版本替换了Windows Rblas.dll，将时间从0.5秒缩短到0.35秒，如上所示。

关于r - R中的外部/张量积，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/8764954/