3

Я программирую expectation-maximization algorithm с 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) 
} 

ответ

4

Для первого я предполагаю, что MyF - это функция, которую вы создали? Если вы можете убедиться, что он будет принимать ваши матрицы и списки в качестве входов и выходов матрицы, вы могли бы сделать что-то вроде:

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

Для второго, я думаю, потому что вы делаете это, как список будет более сложным для векторизации. Может быть, вместо списка матриц у вас может быть 3-мерный массив? Так что MyDf.array [i, j, k] имеет размеры N, d, d (или d, d, N).

2

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

+0

Спасибо за совет. К сожалению, здесь нет симметрии. – Wok 2010-11-22 18:31:38

3

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

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

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

+0

Спасибо за совет. `myF` - это вызов функции плотности (сокращенный pdf в Википедии) [Многомерное нормальное распределение] (http://en.wikipedia.org/wiki/Multivariate_normal_distribution). Это быстрый лайнер. Это действительно цикл по N, который занимает большую часть времени. N равно 500, в то время как k равно 4. – Wok 2010-11-22 19:01:11

+0

Я также предложил бы собрать весь пример и отправить в R-список. Я больше не делаю R, но E-M - это хорошо понятая область! (например, с простой функцией e-m, такой как mvnorm!) Я предполагаю, что это Solved (tm). – 2010-11-22 19:05:04

 Смежные вопросы

  • Нет связанных вопросов^_^