Как ускорить применение функции в слишком много циклов - PullRequest
Новая
       23

Как ускорить применение функции в слишком много циклов

2 голосов
/ 03 ноября 2019

Во-первых, я хотел бы извиниться. Я изучаю R самостоятельно, поэтому я не смог упростить свою проблему и решил просто написать короткую версию моих реальных переменных здесь. Я пытаюсь реализовать вариант классификатора максимального правдоподобия в R. Итак, у меня есть несколько переменных для каждого класса, записанных в векторах и списках (каждая позиция относится к одному классу), и я хочу применить функцию к строкамматрица, которая содержит данные, которые я хочу классифицировать. Проблема в том, что мне нужны результаты этой функции, разделенные по классам. Пока я делаю это: 

cc<-vector(length=2)

mm<-list(length=2)

ii<-list(length=2)

temp1<-matrix(nrow=16,ncol=6)
temp1<-as.data.frame(temp1)
temp1[]<-c(256,235,194,235,215,173,215,215,194,215,215,215,194,173,152,215,
           430,388,388,388,388,430,430,430,388,346,346,388,388,388,346,388,
           283,317,283,283,248,283,283,283,214,214,248,283,214,283,214,248,
           3701,3450,3576,3826,3534,3450,3868,4035,3450,3493,3450,3701,3534,3242,3032,3116,
           1646,1589,1589,1646,1646,1589,1646,1732,1560,1475,1589,1589,1675,1532,1503,1418,
           474,556,556,515,556,556,597,637,556,515,515,515,515,515,434,434)


temp2<- matrix(nrow=11,ncol=6)
temp2<-as.data.frame(temp2)
temp2[]<-c(422,463,462,483,546,525,483,566,546,483,546,
           770,812,770,812,854,854,812,939,939,854,981,
           1038,1175,1004,1141,1209,1209,1038,1311,1311,1175,1311,
           2359,2359,2275,2359,2359,2359,2359,2401,2359,2401,2401,
           2445,2531,2417,2588,2759,2617,2388,2674,2730,2645,2731,
           1413,1413,1373,1495,1618,1535,1413,1535,1659,1535,1618)


cc[1]<-det(cov(temp1))
cc[2]<-det(cov(temp2))

mm[[1]]<-as.numeric(sapply(temp1,"mean"))
mm[[2]]<-as.numeric(sapply(temp2,"mean"))



ii[[1]]<-solve(cov(temp1))
ii[[2]]<-solve(cov(temp2))




data<-matrix(nrow=10,ncol=6)
data<-as.data.frame(data)
data[]<-c(181,203,224,203,203,224,181,181,161,161,
          338,338,338,338,296,296,338,381,338,296,
          208,242,208,208,208,208,208,242,208,173,
          3164,2954,2660,2787,2744,2787,2534,3457,2870,2912,
          1476,1505,1391,1332,1304,1391,1132,1591,1448,1304,
          474,474,474,515,392,432,432,556,515,474)



for (k in 1:2){
  Pxi<-apply(data,1,function(x)1/(2*pi^(6/2)*cc[k]^(1/2))*exp(-1/2*t(as.numeric(x-mm[[k]]))%*%ii[[k]]%*%(as.numeric(x-mm[[k]]))))

  if (k==1) {rule<-Pxi} else {rule<-cbind(rule,Pxi)}  
}

Итак, я понял: 

rule
              rule           Pxi
 [1,] 4.316396e-13  0.000000e+00
 [2,] 6.835553e-15 7.970888e-284
 [3,] 8.674921e-21 2.687251e-145
 [4,] 5.923777e-19 8.020048e-189
 [5,] 5.627127e-16 8.064007e-184
 [6,] 2.495667e-17 5.738550e-209
 [7,] 6.311390e-22  8.913098e-97
 [8,] 1.413893e-12  0.000000e+00
 [9,] 5.521715e-15 1.619401e-221
[10,] 5.212091e-17 5.810407e-254

Ну, как вы можете себе представить, data на самом деле намного большечем в моем примере, и этот последний цикл занимает очень много времени, когда k слишком велик. Любые предложения о том, как сделать это быстрее?

Ответы [ 2 ]

1 голос
/ 05 ноября 2019

Должно быть быстрее, если вы работаете в матрицах. Вот предложение заменить for цикл 

data <- as.matrix(data)
const <- 2*pi^(6/2)
do.call(cbind, lapply(1L:2L, function(k) {
    m <- sweep(data, 2L, mm[[k]])
    #1/(const*cc[k]^(1/2))* exp(-1/2 * diag(m %*% ii[[k]] %*% t(m)))
    1/(const*cc[k]^(1/2))* exp(-1/2 * rowSums((m %*% ii[[k]]) * m))
}))

Использование rowSums (вместо оригинального diag(m %*% ii[[k]] %*% t(m)) было из вычислить только диагонали умножения матриц в R

вывод: 

              [,1]          [,2]
 [1,] 4.316396e-13  0.000000e+00
 [2,] 6.835553e-15 7.970888e-284
 [3,] 8.674921e-21 2.687251e-145
 [4,] 5.923777e-19 8.020048e-189
 [5,] 5.627127e-16 8.064007e-184
 [6,] 2.495667e-17 5.738550e-209
 [7,] 6.311390e-22  8.913098e-97
 [8,] 1.413893e-12  0.000000e+00
 [9,] 5.521715e-15 1.619401e-221
[10,] 5.212091e-17 5.810407e-254
1 голос
/ 03 ноября 2019

Использование cbind() в цикле очень дорого. Вместо этого вы должны назначить результаты промежуточного цикла списку, а затем do.call(cbind, rule) в конце:

Что касается медленного оператора apply(), для каждой строки необходимо выполнить много операций. data перебрал. Вместо этого лучше попытаться выполнить матричные операции (или функцию) одновременно.

Используется функция mahalanobis() для упрощения вызова exp(). Оказывается, что функция использует тот же точный подход, что и @ chinsoon12. 

1 / (2*pi^(6/2)*det(cov(temp1))^(1/2))*exp(-1 / 2 * mahalanobis(data, colMeans(temp1), cov(temp1)))

mahalanobis
#function (x, center, cov, inverted = FALSE, ...) 
#{
#    x <- if (is.vector(x)) 
#        matrix(x, ncol = length(x))
#    else as.matrix(x)
#    if (!isFALSE(center)) 
#        x <- sweep(x, 2L, center)
#    if (!inverted) 
#        cov <- solve(cov, ...)
#    setNames(rowSums(x %*% cov * x), rownames(x))
#}
#<bytecode: 0x000000000c217d80>
#<environment: namespace:stats>

Я бы подошел к этому, сначала сделав list() ваших временных data.frames, а затем используя lapply() для их циклического просмотра: 

tmps <- list(temp1, temp2)

do.call(cbind,
        lapply(tmps,
               function(tmp) {
                 n = length(tmp)
                 cov_tmp <- cov(tmp)
                 1 / (2*pi^(n/2)*det(cov_tmp)^(1/2))*exp(-1 / 2 * mahalanobis(data, colMeans(tmp), cov_tmp))
                 }
               )
)

              [,1]          [,2]
 [1,] 4.316396e-13  0.000000e+00
 [2,] 6.835553e-15 7.970888e-284
 [3,] 8.674921e-21 2.687251e-145
 [4,] 5.923777e-19 8.020048e-189
 [5,] 5.627127e-16 8.064007e-184
 [6,] 2.495667e-17 5.738550e-209
 [7,] 6.311390e-22  8.913098e-97
 [8,] 1.413893e-12  0.000000e+00
 [9,] 5.521715e-15 1.619401e-221
[10,] 5.212091e-17 5.810407e-254

Ссылка: http://sar.kangwon.ac.kr/etc/rs_note/rsnote/cp11/cp11-7.htm

...