arrays - 无循环的动态矩阵乘法

Question

我遇到了以下问题:

---

假设我有三对序列:

{a_1, ..., a_N}, {A_1, ..., A_T}, {b_1, ..., b_N}, {B_1, ..., B_T}, {c_1, . .., c_N}, {C_1, ..., C_T}.

我的目标是执行以下操作(不循环!):

for (i in 1:N) {
  for (j in 1:N) {
    for (k in 1:N) {
      ret[i,j,k] <- \sum_{t=1}^T (a_i - A_t) * (b_j - B_t) * (c_k - C_t)
}}}

我不想循环的原因是可能有比这三个序列更多的序列对。我希望尽可能“高效和灵活”地构建代码。

如果我们只有两对序列，那么问题就很简单了，因为它简化为简单的矩阵乘法((a_i - A_t)< 的(N x T) 矩阵 和 (N x T) 矩阵用于 (b_i - B_t)，其中第一个乘以第二个的转置)。

但是一旦你有两对以上的序列，我不确定它是否可以在没有循环的情况下完成，因为 output 的维数取决于序列对的数量......

---

---------------------------------------- 相关问题(十一月2013 年 8 号)----------------------------------------

感谢 @-mrip，我成功实现了第一部分。但是如果我想要以下内容，代码将如何更改:

for (i in 1:N) {
  for (j in 1:N) {
    for (k in 1:N) {
      ret[i,j,k] <- \sum_{t=1}^T foo(a_i, a_i - A_t) * foo(b_j, b_j - B_t) * foo(c_k, c_k - C_t)
}}}

其中 foo(a, a-A) 是一些常见的双变量函数。是否有“通用”解决方案，或者您是否需要有关 foo(a, a-A) 结构的更多信息？

我尝试使用简单的解决方案并简单地实现循环。当然，这既不灵活(因为我必须提前限制可能的对/维度数量)也不快速(因为 a 和 A 都可能很大 -尽管 a 可能只是一个标量，而 A 可能是一些观察集)。

我知道我可能要求很高。但是很长一段时间以来我完全陷入了这个问题...因此非常欢迎任何帮助。

Answer 1

这实际上根本不需要任何矩阵乘法。它只需要取外积，利用R矩阵和多维数组的列优先布局，可以高效简洁地拼凑出最终结果。通过将循环中的乘积扩展为单独的被加数，可以提高计算效率。这导致以下实现。

 vecout<-function(...)as.vector(outer(...))

 f2<-function(a,b,c,A,B,C){
 N<-length(a)
 t<-length(A)

 ab<-vecout(a,b)
 ret<-array(vecout(ab,c),c(N,N,N))*t

 ret<-ret - ab * sum(C)
 ret<-ret - vecout(a,rep(c,each = N)) * sum(B)
 ret<-ret - rep(vecout(b,c) * sum(A),each=N)
 ret<-ret + a * sum(B*C)
 ret<-ret + rep(b * sum(A*C),each=N)
 ret<-ret + rep(c * sum(A*B),each=N^2)
 ret<-ret - sum(A*B*C)
 ret
 }

我们可以比较运行时间并检查正确性，如下所示。这是天真的实现:

f1<-function(a,b,c,A,B,C){
 N<-length(a)
 ret<-array(0,c(N,N,N))
 for(i in 1:N)
  for(j in 1:N)
   for(k in 1:N)
    ret[i,j,k]<-sum((a[i]-A)*(b[j]-B)*(c[k]-C))
 ret
 }

优化后的版本运行速度大大加快并产生相同的结果，但数值精度错误:

> a<-rnorm(100)
> b<-rnorm(100)
> c<-rnorm(100)
> A<-rnorm(200)
> B<-rnorm(200)
> C<-rnorm(200)
> system.time(r1<-f1(a,b,c,A,B,C))
   user  system elapsed 
  9.006   1.125  10.204 
> system.time(r2<-f2(a,b,c,A,B,C))
   user  system elapsed 
  0.203   0.033   0.242 
> max(abs(r1-r2))
[1] 1.364242e-12

如果每个序列都超过三个，同样的想法也行得通。编写通用解决方案应该不会太难，事实上，用比硬编码的 3 序列解决方案更少的代码行编写通用解决方案是可能的，尽管需要一点思考才能得到所有的索引操作都恰到好处。

关于编辑:好的，无法抗拒。这是一个通用的解决方案，具有任意数量的对，作为两个矩阵的列传递:

f3<-function(a,A){
  subsets<-as.matrix(expand.grid(rep(list(c(F,T)),ncol(a))))
  ret<-array(0,rep(nrow(a),ncol(a)))

  for(i in 1:nrow(subsets)){
    sub<-as.logical(subsets[i,])
    temp<-Reduce(outer,as.list(data.frame(a[,sub,drop=F])),init=1)
    temp<-temp*sum(apply(A[,!sub,drop=F],1,prod))
    temp<-aperm(array(temp,dim(ret)),order(c(which(sub),which(!sub))))
    ret<-ret+temp*(-1)^(sum(!sub))
  } 
  ret
}

> system.time(r3<-f3(cbind(a,b,c),cbind(A,B,C)))
   user  system elapsed 
  0.258   0.056   0.303 
> max(abs(r3-r1))
[1] 9.094947e-13

--------------------再次编辑(2013 年 11 月 8 日)-------------------- --

当数组 A、B 和 C 很大时，上面给出的答案是最有效的。如果A、B和C远小于a、b，和 c，那么这是一个替代的、更简洁的解决方案:

f4<-function(a,A){
  ret<-array(0,rep(nrow(a),ncol(a)))
  for(i in 1:nrow(A)){
    temp<- Reduce(outer,as.list(data.frame(a-rep(A[i,],each=nrow(a)))),init=1)
    ret<-ret + as.vector(temp )
  }
  ret
}

在上面的设置中，a、b、c 的长度分别为 100 和 A、B、C 的长度为 200，这比其他解决方案慢:

> system.time(r3<-f3(cbind(a,b,c),cbind(A,B,C)))
   user  system elapsed 
  0.704   0.092   0.256 
> system.time(r4<-f4(cbind(a,b,c),cbind(A,B,C)))
   user  system elapsed 
 65.824  19.060   3.553 
> max(abs(r3-r4))
[1] 2.728484e-12

但是，如果 A、B 和 C 的长度为 1，那么它会快得多:

> A<-rnorm(1)
> B<-rnorm(1)
> C<-rnorm(1)
> system.time(r3<-f3(cbind(a,b,c),cbind(A,B,C)))
   user  system elapsed 
  0.796   0.172   0.222 
> system.time(r4<-f4(cbind(a,b,c),cbind(A,B,C)))
   user  system elapsed 
  0.180   0.012   0.017 
> max(abs(r3-r4))
[1] 7.105427e-15