R中的外/张量产品

Question

给定每个维度的p向量,计算其张量/外部/ Kruskal产品的最佳方法是什么(带有条目的-array X ？循环是微不足道的,但是很愚蠢.使用重复调用可以正常工作,但似乎不是最佳解决方案(显然,随着p增加,速度会变慢).有更好的方法吗？x1,x2,...,xpdpX[i1,i2,..ip] = x1[i1]x2[i2]...xp[ip])outer

编辑:

我目前最好的是

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

Answer 1

它将难以击败其表现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)

我Rblas.dll在这个地方使用DYNAMIC_ARCH版本的GOTO BLAS 替换了我的Windows ,如上所示,它将时间从0.5秒提高到0.35秒.