Acelere a operação do loop em R

193

Eu tenho um grande problema de desempenho em R. Eu escrevi uma função que itera sobre um data.frameobjeto. Ele simplesmente adiciona uma nova coluna a data.framee acumula algo. (operação simples). O data.frametem aproximadamente 850K linhas. Meu PC ainda está funcionando (cerca de 10 horas agora) e não faço ideia do tempo de execução.

dayloop2 <- function(temp){
    for (i in 1:nrow(temp)){    
        temp[i,10] <- i
        if (i > 1) {             
            if ((temp[i,6] == temp[i-1,6]) & (temp[i,3] == temp[i-1,3])) { 
                temp[i,10] <- temp[i,9] + temp[i-1,10]                    
            } else {
                temp[i,10] <- temp[i,9]                                    
            }
        } else {
            temp[i,10] <- temp[i,9]
        }
    }
    names(temp)[names(temp) == "V10"] <- "Kumm."
    return(temp)
}

Alguma idéia de como acelerar essa operação?

Kay
fonte

Respostas:

433

O maior problema e raiz da ineficácia é a indexação de data.frame, quero dizer todas essas linhas em que você usa temp[,].
Tente evitar isso o máximo possível. Peguei sua função, mudei de indexação e aqui version_A

dayloop2_A <- function(temp){
    res <- numeric(nrow(temp))
    for (i in 1:nrow(temp)){    
        res[i] <- i
        if (i > 1) {             
            if ((temp[i,6] == temp[i-1,6]) & (temp[i,3] == temp[i-1,3])) { 
                res[i] <- temp[i,9] + res[i-1]                   
            } else {
                res[i] <- temp[i,9]                                    
            }
        } else {
            res[i] <- temp[i,9]
        }
    }
    temp$`Kumm.` <- res
    return(temp)
}

Como você pode ver, crio um vetor resque reúne resultados. No final, eu o adiciono data.framee não preciso mexer em nomes. Então, como é melhor?

Eu corro cada função data.framecom nrowde 1.000 a 10.000 por 1.000 e meço o tempo comsystem.time

X <- as.data.frame(matrix(sample(1:10, n*9, TRUE), n, 9))
system.time(dayloop2(X))

O resultado é

desempenho

Você pode ver que sua versão depende exponencialmente de nrow(X). A versão modificada possui relação linear e o lmmodelo simples prevê que, para 850.000 linhas, a computação leva 6 minutos e 10 segundos.

Poder da vetorização

Como Shane e Calimo declaram em suas respostas, a vetorização é a chave para um melhor desempenho. No seu código, você pode sair do loop:

  • condicionamento
  • inicialização dos resultados (que são temp[i,9])

Isso leva a esse código

dayloop2_B <- function(temp){
    cond <- c(FALSE, (temp[-nrow(temp),6] == temp[-1,6]) & (temp[-nrow(temp),3] == temp[-1,3]))
    res <- temp[,9]
    for (i in 1:nrow(temp)) {
        if (cond[i]) res[i] <- temp[i,9] + res[i-1]
    }
    temp$`Kumm.` <- res
    return(temp)
}

Compare o resultado dessas funções, desta vez nrowde 10.000 a 100.000 por 10.000.

desempenho

Ajustando o sintonizado

Outro ajuste é a alteração de uma indexação de loop temp[i,9]para res[i](que são exatamente iguais na i-ésima iteração de loop). Novamente, é uma diferença entre indexar um vetor e indexar a data.frame.
Segunda coisa: quando você olha no loop, pode ver que não há necessidade de repetir tudo i, mas apenas os que se encaixam na condição.
Aqui vamos nos

dayloop2_D <- function(temp){
    cond <- c(FALSE, (temp[-nrow(temp),6] == temp[-1,6]) & (temp[-nrow(temp),3] == temp[-1,3]))
    res <- temp[,9]
    for (i in (1:nrow(temp))[cond]) {
        res[i] <- res[i] + res[i-1]
    }
    temp$`Kumm.` <- res
    return(temp)
}

O desempenho que você obtém altamente depende de uma estrutura de dados. Precisamente - em porcentagem de TRUEvalores na condição. Para meus dados simulados, leva tempo de computação para 850.000 linhas abaixo de um segundo.

desempenho

Eu quero que você possa ir mais longe, vejo pelo menos duas coisas que podem ser feitas:

  • escreva um Ccódigo para fazer cumsum condicional
  • se você sabe que na sequência máxima de dados não é grande, pode alterar o loop para vetorizado enquanto, algo como

    while (any(cond)) {
        indx <- c(FALSE, cond[-1] & !cond[-n])
        res[indx] <- res[indx] + res[which(indx)-1]
        cond[indx] <- FALSE
    }

O código usado para simulações e figuras está disponível no GitHub .

Marek
fonte
2
Como não consigo encontrar uma maneira de perguntar a Marek em particular, como esses gráficos foram gerados?
carbontwelve
@carbontwelve Você está perguntando sobre dados ou parcelas? As parcelas foram feitas com embalagem de treliça. Se eu tiver tempo, coloco o código em algum lugar da Web e aviso.
Marek
@carbontwelve Opa, eu estava errado :) Estes são gráficos padrão (da base R).
Marek
@ Gregor Infelizmente não. É cumulativo, portanto você não pode vetorizá-lo. Exemplo simples: res = c(1,2,3,4)e condé tudo TRUE, então resultado final deve ser: 1, 3(causa 1+2), 6(segunda causa é agora 3, e em terceiro lugar é 3também), 10( 6+4). Fazendo somatório simples que você tem 1, 3, 5, 7.
Marek
Ah, eu deveria ter pensado nisso com mais cuidado. Obrigado por me mostrar o erro.
Gregor Thomas
132

Estratégias gerais para acelerar o código R

Primeiro, descubra onde realmente está a parte lenta. Não é necessário otimizar o código que não está sendo executado lentamente. Para pequenas quantidades de código, simplesmente pensar nele pode funcionar. Se isso falhar, o RProf e ferramentas de perfil semelhantes podem ser úteis.

Depois de descobrir o gargalo, pense em algoritmos mais eficientes para fazer o que deseja. Os cálculos devem ser executados apenas uma vez, se possível, portanto:

O uso de funções mais eficientes pode produzir ganhos de velocidade moderados ou grandes. Por exemplo, paste0produz um pequeno ganho de eficiência, mas .colSums()seus familiares produzem ganhos um pouco mais pronunciados. meané particularmente lento .

Então você pode evitar alguns problemas particularmente comuns :

  • cbind irá atrasá-lo muito rapidamente.
  • Inicialize suas estruturas de dados e preencha-as, em vez de expandi-las a cada vez .
  • Mesmo com a pré-alocação, você pode mudar para uma abordagem de passagem por referência em vez de uma abordagem de passagem por valor, mas pode não valer a pena.
  • Dê uma olhada no R Inferno para evitar mais armadilhas.

Tente uma melhor vetorização , o que geralmente pode ajudar, mas nem sempre. Nesse sentido, comandos inerentemente vetorizados como ifelse, diffe similares fornecerão mais melhorias do que a applyfamília de comandos (que fornecem pouco ou nenhum aumento de velocidade em um loop bem escrito).

Você também pode tentar fornecer mais informações para funções R . Por exemplo, use em vapplyvez desapply e especifique colClassesao ler dados baseados em texto . Os ganhos de velocidade serão variáveis, dependendo de quanto você adivinhar.

Em seguida, considere os pacotes otimizados : O data.tablepacote pode produzir enormes ganhos de velocidade sempre que possível, na manipulação de dados e na leitura de grandes quantidades de dados ( fread).

Em seguida, tente obter ganhos de velocidade por meios mais eficientes de chamar R :

  • Compile seu script R. Ou use os pacotes Rae jitem conjunto para compilação just-in-time (Dirk tem um exemplo nesta apresentação ).
  • Verifique se você está usando um BLAS otimizado. Eles fornecem ganhos de velocidade gerais. Honestamente, é uma pena que o R não use automaticamente a biblioteca mais eficiente na instalação. Espero que o Revolution R contribua com o trabalho que eles fizeram aqui de volta à comunidade em geral.
  • Radford Neal fez várias otimizações, algumas das quais foram adotadas no R Core e muitas outras que foram implementadas no pqR .

E, por último, se todas as opções acima ainda não o levarem tão rápido quanto necessário, você pode precisar mudar para um idioma mais rápido para o snippet de código lento . A combinação de Rcppe inlineaqui facilita muito a substituição da parte mais lenta do algoritmo pelo código C ++. Aqui, por exemplo, é a minha primeira tentativa de fazê-lo , e surpreende até as soluções R altamente otimizadas.

Se você ainda tiver problemas depois de tudo isso, precisará de mais poder de computação. Examine a paralelização ( http://cran.r-project.org/web/views/HighPerformanceComputing.html ) ou mesmo as soluções baseadas em GPU ( gpu-tools).

Links para outras orientações

Ari B. Friedman
fonte
35

Se você estiver usando forloops, provavelmente está codificando R como se fosse C ou Java ou qualquer outra coisa. O código R adequadamente vetorizado é extremamente rápido.

Tomemos, por exemplo, esses dois bits simples de código para gerar uma lista de 10.000 números inteiros em sequência:

O primeiro exemplo de código é como alguém codificaria um loop usando um paradigma de codificação tradicional. Demora 28 segundos para concluir

system.time({
    a <- NULL
    for(i in 1:1e5)a[i] <- i
})
   user  system elapsed 
  28.36    0.07   28.61 

Você pode obter uma melhoria quase 100 vezes com a ação simples de pré-alocar memória:

system.time({
    a <- rep(1, 1e5)
    for(i in 1:1e5)a[i] <- i
})

   user  system elapsed 
   0.30    0.00    0.29 

Mas, usando a operação do vetor R base, usando o operador dois pontos, :essa operação é praticamente instantânea:

system.time(a <- 1:1e5)

   user  system elapsed 
      0       0       0 
Andrie
fonte
+1, embora eu considere seu segundo exemplo não convincente, pois a[i]não muda. Mas system.time({a <- NULL; for(i in 1:1e5){a[i] <- 2*i} }); system.time({a <- 1:1e5; for(i in 1:1e5){a[i] <- 2*i} }); system.time({a <- NULL; a <- 2*(1:1e5)})tem um resultado semelhante.
Henry
@ Henry, comentário justo, mas como você aponta, os resultados são os mesmos. Modifiquei o exemplo para inicializar um para rep(1, 1e5)- os tempos são idênticos.
28711 Andrie
17

Isso pode ser feito muito mais rapidamente, ignorando os loops usando índices ou ifelse()instruções aninhadas .

idx <- 1:nrow(temp)
temp[,10] <- idx
idx1 <- c(FALSE, (temp[-nrow(temp),6] == temp[-1,6]) & (temp[-nrow(temp),3] == temp[-1,3]))
temp[idx1,10] <- temp[idx1,9] + temp[which(idx1)-1,10] 
temp[!idx1,10] <- temp[!idx1,9]    
temp[1,10] <- temp[1,9]
names(temp)[names(temp) == "V10"] <- "Kumm."
Shane
fonte
Obrigado pela resposta. Eu tento entender suas declarações. A linha 4: "temp [idx1,10] <- temp [idx1,9] + temp [que (idx1) -1,10]" causou um erro porque o comprimento do objeto mais longo não é um múltiplo do comprimento do objeto mais curto. "temp [idx1,9] = num [1: 11496]" e "temp [which (idx1) -1,10] = int [1: 11494]", portanto, faltam 2 linhas.
Kay
Se você fornecer uma amostra de dados (use dput () com algumas linhas), eu a corrigirei. Por causa do qual () - 1 bit, os índices são desiguais. Mas você deve ver como funciona a partir daqui: não há necessidade de loop ou aplicação; basta usar funções vetorizadas.
Shane
1
Uau! Acabei de alterar um bloco de função aninhado if..else e mapply, para uma função aninhada ifelse e obtive uma aceleração de 200x!
2626 James
Seu conselho geral é correto, mas no código você perdeu o fato, esse ivalor -th depende de i-1-th, para que não possam ser definidos da maneira que você faz (usando which()-1).
Marek
8

Não gosto de reescrever código ... Também é claro que ifelse e lapply são melhores opções, mas às vezes é difícil ajustá-lo.

Freqüentemente eu uso data.frames como se usasse listas como df$var[i]

Aqui está um exemplo inventado:

nrow=function(x){ ##required as I use nrow at times.
  if(class(x)=='list') {
    length(x[[names(x)[1]]])
  }else{
    base::nrow(x)
  }
}

system.time({
  d=data.frame(seq=1:10000,r=rnorm(10000))
  d$foo=d$r
  d$seq=1:5
  mark=NA
  for(i in 1:nrow(d)){
    if(d$seq[i]==1) mark=d$r[i]
    d$foo[i]=mark
  }
})

system.time({
  d=data.frame(seq=1:10000,r=rnorm(10000))
  d$foo=d$r
  d$seq=1:5
  d=as.list(d) #become a list
  mark=NA
  for(i in 1:nrow(d)){
    if(d$seq[i]==1) mark=d$r[i]
    d$foo[i]=mark
  }
  d=as.data.frame(d) #revert back to data.frame
})

data.frame versão:

   user  system elapsed 
   0.53    0.00    0.53

versão da lista:

   user  system elapsed 
   0.04    0.00    0.03 

17 vezes mais rápido para usar uma lista de vetores que um data.frame.

Quaisquer comentários sobre por que internamente os data.frames são tão lentos nesse sentido? Alguém poderia pensar que eles operam como listas ...

Para um código ainda mais rápido, faça isso em class(d)='list'vez de d=as.list(d)eclass(d)='data.frame'

system.time({
  d=data.frame(seq=1:10000,r=rnorm(10000))
  d$foo=d$r
  d$seq=1:5
  class(d)='list'
  mark=NA
  for(i in 1:nrow(d)){
    if(d$seq[i]==1) mark=d$r[i]
    d$foo[i]=mark
  }
  class(d)='data.frame'
})
head(d)
Chris
fonte
1
Provavelmente é graças à sobrecarga de [<-.data.frame, que é chamada de alguma forma quando você faz d$foo[i] = marke pode acabar fazendo uma nova cópia do vetor de possivelmente todo o data.frame em cada <-modificação. Seria uma pergunta interessante sobre SO.
Frank
2
@Frank (i) deve garantir que o objeto modificado ainda seja um data.frame válido e (ii) a afaik faça pelo menos uma cópia, possivelmente mais de uma. Sabe-se que a sub-atribuição de quadro de dados é lenta e se você olhar para o longo código-fonte, não é realmente surpreendente.
Roland
@Frank, @Roland: A df$var[i]notação passa pela mesma [<-.data.framefunção? Eu notei que é bastante longo. Caso contrário, que função ele usa?
Chris
@ Chris Creio que d$foo[i]=marké traduzido aproximadamente d <- `$<-`(d, 'foo', `[<-`(d$foo, i, mark)), mas com algum uso de variáveis ​​temporárias.
Tim Goodman
7

Como Ari mencionou no final de sua resposta, os pacotes Rcppe inlinetornam incrivelmente fácil acelerar as coisas. Como exemplo, tente este inlinecódigo (aviso: não testado):

body <- 'Rcpp::NumericMatrix nm(temp);
         int nrtemp = Rccp::as<int>(nrt);
         for (int i = 0; i < nrtemp; ++i) {
             temp(i, 9) = i
             if (i > 1) {
                 if ((temp(i, 5) == temp(i - 1, 5) && temp(i, 2) == temp(i - 1, 2) {
                     temp(i, 9) = temp(i, 8) + temp(i - 1, 9)
                 } else {
                     temp(i, 9) = temp(i, 8)
                 }
             } else {
                 temp(i, 9) = temp(i, 8)
             }
         return Rcpp::wrap(nm);
        '

settings <- getPlugin("Rcpp")
# settings$env$PKG_CXXFLAGS <- paste("-I", getwd(), sep="") if you want to inc files in wd
dayloop <- cxxfunction(signature(nrt="numeric", temp="numeric"), body-body,
    plugin="Rcpp", settings=settings, cppargs="-I/usr/include")

dayloop2 <- function(temp) {
    # extract a numeric matrix from temp, put it in tmp
    nc <- ncol(temp)
    nm <- dayloop(nc, temp)
    names(temp)[names(temp) == "V10"] <- "Kumm."
    return(temp)
}

Existe um procedimento semelhante para #includeing coisas, onde você apenas passa um parâmetro

inc <- '#include <header.h>

para cxxfunção, como include=inc. O que é realmente legal nisso é que ele faz todos os links e compilações para você, portanto a criação de protótipos é muito rápida.

Isenção de responsabilidade: não tenho certeza absoluta de que a classe de tmp deva ser numérica e não matricial numérica ou qualquer outra coisa. Mas tenho quase certeza.

Edit: se você ainda precisar de mais velocidade depois disso, o OpenMP é um recurso de paralelização adequado C++. Eu não tentei usá-lo inline, mas deve funcionar. A idéia seria, no caso de nnúcleos, fazer com que a iteração de loop kfosse realizada k % n. A introdução adequada é encontrado em Matloff é A Arte de R programação , disponível aqui , no capítulo 16, Recorrer a C .

jclancy
fonte
3

As respostas aqui são ótimas. Um aspecto menor não abordado é que a pergunta declara " Meu PC ainda está funcionando (cerca de 10 horas agora) e não tenho idéia do tempo de execução ". Eu sempre coloco o código a seguir em loops durante o desenvolvimento para ter uma ideia de como as mudanças parecem afetar a velocidade e também para monitorar quanto tempo levará para ser concluído.

dayloop2 <- function(temp){
  for (i in 1:nrow(temp)){
    cat(round(i/nrow(temp)*100,2),"%    \r") # prints the percentage complete in realtime.
    # do stuff
  }
  return(blah)
}

Funciona com lapply também.

dayloop2 <- function(temp){
  temp <- lapply(1:nrow(temp), function(i) {
    cat(round(i/nrow(temp)*100,2),"%    \r")
    #do stuff
  })
  return(temp)
}

Se a função no loop for bastante rápida, mas o número de loops for grande, considere apenas imprimir de vez em quando, pois a impressão no próprio console tem uma sobrecarga. por exemplo

dayloop2 <- function(temp){
  for (i in 1:nrow(temp)){
    if(i %% 100 == 0) cat(round(i/nrow(temp)*100,2),"%    \r") # prints every 100 times through the loop
    # do stuff
  }
  return(temp)
}
novato
fonte
Uma opção semelhante, imprima a fração i / n. Eu sempre tenho algo assim, cat(sprintf("\nNow running... %40s, %s/%s \n", nm[i], i, n))já que geralmente estou repetindo coisas nomeadas (com nomes nm).
24518 Frank
2

No R, é possível acelerar o processamento do loop usando as applyfunções da família (no seu caso, provavelmente seria replicate). Dê uma olhada no plyrpacote que fornece barras de progresso.

Outra opção é evitar completamente os loops e substituí-los por aritmética vetorizada. Não sei exatamente o que você está fazendo, mas você provavelmente pode aplicar sua função a todas as linhas de uma vez:

temp[1:nrow(temp), 10] <- temp[1:nrow(temp), 9] + temp[0:(nrow(temp)-1), 10]

Isso será muito mais rápido e, em seguida, você poderá filtrar as linhas com sua condição:

cond.i <- (temp[i, 6] == temp[i-1, 6]) & (temp[i, 3] == temp[i-1, 3])
temp[cond.i, 10] <- temp[cond.i, 9]

A aritmética vetorizada requer mais tempo e reflexão sobre o problema, mas às vezes você pode salvar várias ordens de magnitude no tempo de execução.

Calimo
fonte
14
você sabe que as funções vetoriais serão mais rápidas que os loops ou os apply (), mas não é verdade que o apply () seja mais rápido que os loops. Em muitos casos, apply () está simplesmente abstraindo o loop do usuário, mas ainda fazendo loop. Veja esta pergunta anterior: stackoverflow.com/questions/2275896/…
JD Long
0

Processar com data.tableé uma opção viável:

n <- 1000000
df <- as.data.frame(matrix(sample(1:10, n*9, TRUE), n, 9))
colnames(df) <- paste("col", 1:9, sep = "")

library(data.table)

dayloop2.dt <- function(df) {
  dt <- data.table(df)
  dt[, Kumm. := {
    res <- .I;
    ifelse (res > 1,             
      ifelse ((col6 == shift(col6, fill = 0)) & (col3 == shift(col3, fill = 0)) , 
        res <- col9 + shift(res)                   
      , # else
        res <- col9                                 
      )
     , # else
      res <- col9
    )
  }
  ,]
  res <- data.frame(dt)
  return (res)
}

res <- dayloop2.dt(df)

m <- microbenchmark(dayloop2.dt(df), times = 100)
#Unit: milliseconds
#       expr      min        lq     mean   median       uq      max neval
#dayloop2.dt(df) 436.4467 441.02076 578.7126 503.9874 575.9534 966.1042    10

Se você ignorar os possíveis ganhos com a filtragem de condições, é muito rápido. Obviamente, se você pode fazer o cálculo no subconjunto de dados, isso ajuda.

Bulat
fonte
2
Por que você está repetindo a sugestão de usar data.table? Já foi feito várias vezes nas respostas anteriores.
IRTFM