Как ускорить этот вид двойного цикла?

Question

Я программирую алгоритм максимизации ожидания с R. Чтобы ускорить вычисления, я хотел бы векторизовать это узкое место. Я знаю, что N примерно в сто раз больше.

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)
}

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).

Answer 2

Я ненавижу даже предлагать это преждевременно, но это такая вещь, где создание C-расширения в R может иметь смысл. Для матриц с определенным (известным) размером (который у вас есть здесь!) C-расширения не сложнее построить, обещаю! Самый неприятный бит здесь, вероятно, будет передаваться в «myF»

Мои R-знания довольно устарели, но для циклов (особенно таких, как этот!) Раньше были жестокими.

Может быть, поможет время и выяснить, какая часть медленная? Это myF? Что делать, если вы измените его на личность?

Answer 3

Вы можете сократить работу, выполненную во внутреннем цикле, если вещи симметричны: A[i,j] = A[j,i]