我陷入了以下问题：
假设我有三对序列：
{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)
}}}
我不想循环的原因是，可能会有更多的序列对，而不仅仅是这三对。而且我想使代码结构尽可能“高效且灵活”。
如果我们只有两对序列，那么这个问题就很容易，因为它简化为简单的矩阵乘法（用于的(N x T)矩阵(a_i - A_t)和(N x T)用于(b_i - B_t)将第一项与第二项的转置相乘的矩阵）。
但是一旦您拥有多于两对序列，我就不确定是否可以不使用循环而完成，因为的维数output取决于序列对的数量...
------------------------------------------相关问题（2013年11月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一些组观察）。
我知道我可能会要求。但是很长时间以来，我对这个问题一无所知。因此，非常欢迎您提供任何帮助。
实际上，这根本不需要任何矩阵乘法。它只需要获取外部乘积，并且可以利用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具有length 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
本文收集自互联网，转载请注明来源。
如有侵权，请联系[email protected] 删除。
我来说两句