如何加速这种双重循环的速度？（涉及IT技术）

Question

我正在使用R编程一个期望最大化算法。为了加速计算，我希望对这个瓶颈进行向量化处理。我知道N大约是k的一百倍。

MyLoglik = 0
for (i in c(1:N))
{
 for (j in c(1:k))
 {
  MyLoglik = MyLoglik + MyTau[i,j]*log(MyP[j]*MyF(MyD[i,], MyMu[j,], MyS[[j]]))
 }
}

还有这个矩阵列表：

MyDf.list <- vector("list", k)
for(i in 1:k)
{
 MyDf.list[[i]] <- matrix(0,d,d)
 for (j in c(1:N))
 {
  MyDf.list[[i]] = MyDf.list[[i]] + MyTau[j,i]*as.numeric((MyD[j,]-MyMu[i,])) %*% t(as.numeric(MyD[j,]-MyMu[i,]))  
 }
 MyDf.list[[i]] = MyDf.list[[i]] / MyM[i]
}

---

我稍微加快了速度，使用了：

MyLoglik = 0
for (j in c(1:k))
{
 MyR= apply(MyD, 1, function(x) log(MyP[j]*MyF(x, MyMu[j,], MyS[[j]])))
 MyLoglik = MyLoglik + sum(MyTau[,j]*MyR)
}

并且：

d = dim(MyD)[2]
MyDf.list <- vector("list", k)
for(i in 1:k)
{
 MyDf.list[[i]] <- matrix(0,d,d)
 MyR= apply(MyD, 1, function(x) as.numeric((x-MyMu[i,])) %*% t(as.numeric(x-MyMu[i,])))
 MyDf.list[[i]] = matrix(rowSums(t(MyTau[,i]*t(MyR))) / MyM[i],d,d)
}

- Wok

Answer 1

首先，我假设MyF是您编写的一个函数？如果您能确保它可以将矩阵和列表作为输入，并输出一个矩阵，您可以这样做：

MyLoglik = sum(MyTau%*%log(MyP)) + sum(MyTau*log(MyF(MyD, MyMu, MyS)))

对于第二个问题，我认为因为你将其作为列表处理，所以向量化会更加困难。也许你可以使用一个三维数组来代替矩阵的列表？这样MyDf.array[i,j,k]就有了N、d、d（或d、d、N）这些维度。

- Wine

---

第一个建议不错！至于第二个，我认为列表是像Matlab中获取数组的唯一方法。 - Wok

请检查?array - 它可以处理多个维度。 - Wine

Answer 2

我不愿意过早地提出这个建议，但对于与已知大小的矩阵（您在这里就是这种情况！）构建R中的C扩展可能是有意义的。 我保证，构建C扩展并不是那么难！ 在这里最困难的部分可能是传递“myF”。

我的R知识已经过时了，但是循环（特别是像这样的循环！）曾经非常棘手。

也许计时和找出哪一部分是慢的会有所帮助？ 是myF吗？ 如果将其更改为identity会怎样？

- Gregg Lind

---

感谢您的建议。myF 是对多元正态分布密度函数（在维基百科上缩写为pdf）的调用。这是一个快速的一行代码。实际上，循环 N 才是最耗时间的部分。其中 N 为 500，而 k 为 4。 - Wok

我建议将整个示例粘贴到pastebin并提交到R-list。虽然我不再经常使用R，但E-M是一个众所周知的领域！（特别是像mvnorm这样简单的e-m函数！）我猜这个问题已经解决了（商标）。 - Gregg Lind

Answer 3

如果事物是对称的，您可以减少内部循环中执行的工作量：A[i,j] = A[j,i]

- duffymo

---

谢谢建议。不幸的是，这里没有对称性。 - Wok