Comment accélérer ce type de double boucle?

Question

Je programme une expectation-maximization algorithm avec R. Afin d'accélérer le calcul, je voudrais vectoriser ce goulot d'étranglement. Je sais que N est environ cent fois 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]])) 
} 
} 

Il y a aussi cette liste de matrices:

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

---

J'ai accéléré les choses un peu en utilisant:

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

et:

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

Pour le premier, je suppose que MyF est une fonction que vous avez créée? Si vous pouvez vous assurer qu'il prendra vos matrices et des listes des entrées et la sortie d'une matrice, vous pouvez faire quelque chose comme:

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

Pour la seconde, je pense que vous faites comme liste, il sera plus difficile à vectoriser. Peut-être qu'au lieu d'une liste de matrices, vous pourriez avoir un tableau en trois dimensions? Alors que MyDf.array [i, j, k] a des dimensions N, d, d (ou d, d, N).

Answer 2

Vous pouvez réduire le travail effectué dans la boucle interne si les choses sont symétriques: A[i,j] = A[j,i]

Answer 3

Je déteste même suggérer cela prématurément, mais c'est le genre de chose où la construction d'une extension C dans R pourrait avoir du sens. Pour les matrices avec la taille définie (connue) (que vous avez ici!), Les extensions C ne sont pas que dur pour construire, je promets! Le plus méchant ici passerait probablement dans 'myF'

Ma connaissance de R est assez obsolète, mais pour les boucles (particulièrement celle-là!) C'était brutal.

Peut-être que le chronométrage et la détermination de la partie lente seraient utiles? Est-ce myF? Que faire si vous le changez en une identité?