Perché la complessità temporale di questo ciclo non è lineare?

Question

Perché la complessità temporale di questo ciclo non è lineare e perché è così lento? Il ciclo prende ~38s for N=50k, e ~570s for N=200k. C'è un modo più veloce per farlo? Rprof() sembra indicare che scrivere in memoria è molto lento.

df <- data.frame(replicate(5, runif(200000))) 
df[,1:3] <- round(df[,1:3]) 

Rprof(line.profiling = TRUE); timer <- proc.time() 
x <- df; N <- nrow(df); i <- 1 
ind <- df[1:(N-1),1:3] == df[2:N,1:3]; 
rind <- which(apply(ind,1,all)) 
N <- length(rind) 
while(i <= N) 
{ 
    x$X4[rind[i]+1] <- x$X4[rind[i]+1] + x$X4[rind[i]] 
    x$X5[rind[i]+1] <- x$X4[rind[i]+1] * x$X3[rind[i]+1] 
    x$X5[rind[i]+1] <- trunc(x$X5[rind[i]+1]*10^8)/10^8 
    x$X1[rind[i]] <- NA 
    i <- i + 1 
};x <- na.omit(x) 
proc.time() - timer; Rprof(NULL) 
summaryRprof(lines = "show") 

Lo scopo di questo algoritmo è di scorrere il frame di dati e combinare righe adiacenti che corrispondono a determinati elementi. Cioè, rimuove una delle righe e aggiunge alcuni dei valori di quella riga all'altra riga. Il frame dati risultante dovrebbe avere n righe in meno, dove n è il numero di righe adiacenti corrispondenti nel frame dati originale. Ogni volta che viene combinata una coppia di righe, l'indice del frame di dati di origine e il nuovo frame di dati vengono fuori sincrono di 1, poiché una riga viene rimossa/omessa dal nuovo frame, quindi i tiene traccia della posizione sui dati di origine frame e q tiene traccia della posizione sul nuovo frame di dati.

Il codice sopra è aggiornato grazie al commento di @ joran. Le prestazioni sono sostanzialmente migliorate a ~5.5s for N=50k e ~88s for N=200k. Tuttavia, la complessità del tempo è ancora non lineare, cosa che non riesco a capire. Ho bisogno di eseguire questo a N = 1 milione o più, quindi la sua non ancora grande velocità.

Answer 1

Solo l'aggiornamento X4 colonna dipende dai valori precedenti, quindi il ciclo può essere principalmente 'vettorializzato' (con un po 'di ottimizzazione, evitando aggiunta 1 a rind in ogni iterazione) come

rind1 <- rind + 1L 
for (i in seq_len(N)) 
    x$X4[rind1[i]] <- x$X4[rind1[i]] + x$X4[rind[i]] 

x$X5[rind1] <- x$X4[rind1] * x$X3[rind1] 
x$X5[rind1] <- trunc(x$X5[rind1] * 10^8)/10^8 
x$X1[rind] <- NA 
na.omit(x) 

X4 è un valore numerico e l'aggiornamento può essere resa più efficiente aggiornando come un vettore piuttosto che una colonna di un data.frame

X4 <- x$X4 
for (i in seq_len(N)) 
    X4[rind1[i]] <- X4[rind1[i]] + X4[rind[i]] 
x$X4 <- X4 

Per confronto, abbiamo

f0 <- function(nrow) { 
    set.seed(123) 
    df <- data.frame(replicate(5, runif(nrow))) 
    df[,1:3] <- round(df[,1:3]) 
    x <- df; N <- nrow(df); i <- 1 
    ind <- df[1:(N-1),1:3] == df[2:N,1:3]; 
    rind <- which(apply(ind,1,all)) 
    N <- length(rind) 

    while(i <= N) 
    { 
     x$X4[rind[i]+1] <- x$X4[rind[i]+1] + x$X4[rind[i]] 
     x$X5[rind[i]+1] <- x$X4[rind[i]+1] * x$X3[rind[i]+1] 
     x$X5[rind[i]+1] <- trunc(x$X5[rind[i]+1]*10^8)/10^8 
     x$X1[rind[i]] <- NA 
     i <- i + 1 
    } 
    na.omit(x) 
} 

f1a <- function(nrow) { 
    set.seed(123) 
    df <- data.frame(replicate(5, runif(nrow))) 
    df[,1:3] <- round(df[,1:3]) 
    x <- df; N <- nrow(df) 
    ind <- df[1:(N-1),1:3] == df[2:N,1:3]; 
    rind <- which(apply(ind,1,all)) 

    rind1 <- rind + 1L 
    for (i in seq_along(rind)) 
     x$X4[rind1[i]] <- x$X4[rind1[i]] + x$X4[rind[i]] 

    x$X5[rind1] <- x$X4[rind1] * x$X3[rind1] 
    x$X5[rind1] <- trunc(x$X5[rind1] * 10^8)/10^8 
    x$X1[rind] <- NA 
    na.omit(x) 
} 

f4a <- function(nrow) { 
    set.seed(123) 
    df <- data.frame(replicate(5, runif(nrow))) 
    df[,1:3] <- round(df[,1:3]) 
    x <- df; N <- nrow(df) 
    ind <- df[1:(N-1),1:3] == df[2:N,1:3]; 
    rind <- which(apply(ind,1,all)) 

    rind1 <- rind + 1L 
    X4 <- x$X4 
    for (i in seq_along(rind)) 
     X4[rind1[i]] <- X4[rind1[i]] + X4[rind[i]] 
    x$X4 <- X4 

    x$X1[rind] <- NA 
    x$X5[rind1] <- X4[rind1] * x$X3[rind1] 
    x$X5[rind1] <- trunc(x$X5[rind1] * 10^8)/10^8 

    na.omit(x) 
} 

I risultati sono gli stessi

> identical(f0(1000), f1a(1000)) 
[1] TRUE 
> identical(f0(1000), f4a(1000)) 
[1] TRUE 

L'aumento di velocità è sostanziale (usando library(microbenchmark))

> microbenchmark(f0(10000), f1a(10000), f4a(10000), times=10) 
Unit: milliseconds 
     expr  min  lq  mean median  uq  max neval 
    f0(10000) 346.35906 354.37637 361.15188 363.71627 366.74944 373.88275 10 
f1a(10000) 124.71766 126.43532 127.99166 127.39257 129.51927 133.01573 10 
f4a(10000) 41.70401 42.48141 42.90487 43.00584 43.32059 43.83757 10 

La ragione per la differenza può essere visto quando R è stato compilato con profiling memoria abilitato -

> tracemem(x) 
[1] "<0x39d93a8>" 
> tracemem(x$X4) 
[1] "<0x6586e40>" 
> x$X4[1] <- 1 
tracemem[0x39d93a8 -> 0x39d9410]: 
tracemem[0x6586e40 -> 0x670d870]: 
tracemem[0x39d9410 -> 0x39d9478]: 
tracemem[0x39d9478 -> 0x39d94e0]: $<-.data.frame $<- 
tracemem[0x39d94e0 -> 0x39d9548]: $<-.data.frame $<- 
> 

Ogni riga indica una copia di memoria, quindi l'aggiornamento di una cella in un frame di dati comporta 5 copie della struttura esterna o del vettore stesso. Al contrario, un vettore può essere aggiornato senza alcuna copia.

> tracemem(X4) 
[1] "<0xdd44460>" 
> X4[1] = 1 
tracemem[0xdd44460 -> 0x9d26c10]: 
> X4[1] = 2 
> 

(La prima assegnazione è costoso perché rappresenta la duplicazione della colonna data.frame; successivi aggiornamenti sono X4, solo X4 riferisce al vettore essendo aggiornato, e il vettore non deve essere duplicata) .

attuazione Il data.frame sembra scalare non lineare

> microbenchmark(f1a(100), f1a(1000), f1a(10000), f1a(100000), times=10) 
Unit: milliseconds 
     expr   min   lq  mean  median   uq 
    f1a(100) 2.372266 2.479458 2.551568 2.524818 2.640244 
    f1a(1000) 10.831288 11.100009 11.210483 11.194863 11.432533 
f1a(10000) 130.011104 138.686445 139.556787 141.138329 141.522686 
f1a(1e+05) 4092.439956 4117.818817 4145.809235 4143.634663 4172.282888 
     max neval 
    2.727221 10 
    11.581644 10 
    147.993499 10 
4216.129732 10 

La ragione è evidente nella seconda riga dell'output tracemem sopra - aggiornare una riga innesca una copia dell'intera colonna.Quindi l'algoritmo scala come il numero di righe da aggiornare per il numero del numero di righe in una colonna, approssimativamente quadratico.

f4a() sembra scalare linearmente

> microbenchmark(f4a(100), f4a(1000), f4a(10000), f4a(100000), f4a(1e6), times=10) 
Unit: milliseconds 
     expr   min   lq  mean  median   uq 
    f4a(100) 1.741458 1.756095 1.827886 1.773887 1.929943 
    f4a(1000) 5.286016 5.517491 5.558091 5.569514 5.671840 
f4a(10000) 42.906895 43.025385 43.880020 43.928631 44.633684 
f4a(1e+05) 467.698285 478.919843 539.696364 552.896109 576.707913 
f4a(1e+06) 5385.029968 5521.645185 5614.960871 5573.475270 5794.307470 
     max neval 
    2.003700 10 
    5.764022 10 
    44.983002 10 
    644.927832 10 
5823.868167 10 

Si potrebbe cercare di essere intelligente sulla vettorizzazione del ciclo, ma è ora necessario?

Una versione accordata parte di elaborazione di dati della funzione di indicizzazione utilizza negativo (ad esempio, -nrow(df)) per rimuovere le righe dal frame di dati, invece di rowSums()apply() e unname() modo che le operazioni di sottoinsieme non portano nomi intorno inutilizzati:

g0 <- function(df) { 
    ind <- df[-nrow(df), 1:3] == df[-1, 1:3] 
    rind <- unname(which(rowSums(ind) == ncol(ind))) 
    rind1 <- rind + 1L 

    X4 <- df$X4 
    for (i in seq_along(rind)) 
     X4[rind1[i]] <- X4[rind1[i]] + X4[rind[i]] 

    df$X4 <- X4 
    df$X1[rind] <- NA 
    df$X5[rind1] <- trunc(df$X4[rind1] * df$X3[rind1] * 10^8)/10^8 

    na.omit(df) 
} 

Rispetto alla soluzione suggerita da data.table @Khashaa

g1 <- function(df) { 
    x <- setDT(df)[, r:=rleid(X1, X2, X3),] 
    x <- x[, .(X1=X1[.N], X2=X2[.N], X3=X3[.N], X4=sum(X4), X5=X5[.N]), by=r] 
    x <- x[, X5:= trunc(X3 * X4 * 10^8)/10^8] 
    x 
} 

la versione base R esibisce favorevolmente con tempi

> n_row <- 200000 
> set.seed(123) 
> df <- data.frame(replicate(5, runif(n_row))) 
> df[,1:3] <- round(df[,1:3]) 
> system.time(g0res <- g0(df)) 
    user system elapsed 
    0.247 0.000 0.247 
> system.time(g1res <- g1(df)) 
    user system elapsed 
    0.551 0.000 0.551 

(La versione di pre-tuning in f4a richiede circa 760 ms, quindi più di due volte più lentamente).

I risultati dell'attuazione data.table non sono corretti

> head(g0res) 
    X1 X2 X3  X4  X5 
1 0 1 1 0.4708851 0.8631978 
2 1 1 0 0.8977670 0.8311355 
3 0 1 0 0.7615472 0.6002179 
4 1 1 1 0.6478515 0.5616587 
5 1 0 0 0.5329256 0.5805195 
6 0 1 1 0.8526255 0.4913130 
> head(g1res) 
    r X1 X2 X3  X4  X5 
1: 1 0 1 1 0.4708851 0.4708851 
2: 2 1 1 0 0.8977670 0.0000000 
3: 3 0 1 0 0.7615472 0.0000000 
4: 4 1 1 1 0.6478515 0.6478515 
5: 5 1 0 0 0.5329256 0.0000000 
6: 6 0 1 1 0.8526255 0.8526255 

e io non sono abbastanza di una procedura guidata data.table (appena un utente data.table) per sapere qual è la formulazione corretta è.

Compilazione (benefici esclusivamente dal ciclo for?) Aumenta la velocità di circa il 20%

> g0c <- compiler::cmpfun(g0) 
> microbenchmark(g0(df), g0c(df), times=10) 
Unit: milliseconds 
    expr  min  lq  mean median  uq  max neval 
    g0(df) 250.0750 262.941 276.1549 276.8848 281.1966 321.3778 10 
    g0c(df) 214.3132 219.940 228.0784 230.2098 235.4579 242.6636 10 

Answer 2

seguito è solo una riscrittura della risposta di @ Martin Morgan, utilizzando il veloce sottoinsiemi di data.table. È circa 3 volte più veloce rispetto all'approccio data.frame.

library(data.table) 
library(matrixStats) # for efficient rowAlls function 

g01 <- function(df) { 
    setDT(df) 
    ind <- df[-nrow(df), 1:3, with=FALSE] == df[-1, 1:3, with=FALSE] 
    rind <- which(rowAlls(ind)) + 1L 

    X4 <- df$X4 
    for (i in seq_along(rind)) 
    X4[rind[i]] <- X4[rind[i]] + X4[rind[i] - 1L] 

    df$X4 <- X4 
    df$X5[rind] <- trunc(df$X4[rind] * df$X3[rind] * 10^8)/10^8 
    df[-rind + 1L,] 
} 

g01c <- compiler::cmpfun(g01) 

n_row <- 1e6 
set.seed(123) 
df <- data.frame(replicate(5, runif(n_row))) 
df[,1:3] <- round(df[,1:3]) 
# data.frame 
system.time(g0(df)) 
# user system elapsed 
# 1.14 0.00 1.14 
system.time(g0c(df)) 
# user system elapsed 
# 0.82 0.03 0.86 

# data.table 
system.time(g01(df)) 
# user system elapsed 
# 0.40 0.02 0.43 
system.time(g01c(df)) 
# user system elapsed 
# 0.12 0.03 0.16