Un forte limite nell'utilizzo di R (in particolare per chi è già abituato a sviluppare in altri linguaggi di programmazione) è la lentezza dell'esecuzione di algoritmi basati sui classici "cicli di for". Eseguire un blocco di codice ripetutamente, infatti, è un'operazione concettualmente semplice che, quindi, ci permette di gestire facilmente la ripetizione di un insieme di "operazioni" (eventualmente effettuando il controllo di alcune condizioni con dei sempici "if"). Putroppo in R è fortemente sconsigliato questo approccio e per questo motivo si parla di vectorizing computing: trasformare il ciclo di for in operazioni tra vettori o matrici. Quindi la ripetizione di un blocco di codice controllato da un contatore deve essere trasformato, ad esempio, in un prodotto di matrici. Questo metodo di sviluppo ha come fine quello di sfruttare al meglio il motore di calcolo di R, che è appunto ottimizzato per il calcolo matriciale. Ovviamente questo approccio può portare ad altro tipo di problemi, ossia alla gestione dei dati mediante matrici "troppo grandi" e difficilmente gestibili. L'esperienza comunque permette di arrivare alla soluzione megliore, e non credo si possa stilare un elenco di situazioni in cui preferire un approccio ad un altro: ogni soluzione va valutata singolarmente "sul campo" (...bella novità...).
Riporto di seguito alcuni esempi di for in R con le relative "vettorizzazioni".
Creo un vettore p composto da 500.000 valori estratti da una normale standardizzata e, da questo, ne derivo il vettore formato dagli stessi 500.000 valori rapportati alla loro somma. Forse la vettorizzazione di questo esempio è veramente banale, ma serve a rendersi conto dei tempi di esecuzione.
Il relativo ciclo di for è semplicissimo:
p<-rnorm(500000)
s<-sum(p)
for (i in 1:500000) p[i]<-p[i]/s
Per la valutazione dei tempi di esecuzione basta osservare il risultato di system.time(),che esegue il codice in parentesi e riporta il tempo in secondi. Nel nostro caso ci limitiamo ad osservare il tempo elapsed (in blu è roportato l'output di R):
system.time(for (i in 1:500000) p[i]<-p[i]/s)
user system elapsed
3.2 0.0 3.2
La cosiddetta vettorizzazione è in tal caso banale (p/s), ma il risultato è notevole:
system.time(p<-p/s)
user system elapsed
0.0 0.02 0.02
In termini assoulti, un risparmio di 3.18 secondi può sembrare inutile, ma se si ragiona in termini relativi si osserva un risparmio del 99% del tempo di esecuzione, quindi molto importante per calcoli più complessi!
Prendiamo in esame un esempio in cui la vettorizzazione è meno immediata.
Consideriamo un insieme di 1.000.000 di dati estratti da una normale standardizzata. Pensate ad esempio ad un tipica applicazione in ambito industriale, in cui si controlla che il prodotto rispetti determinate caratteristiche (peso, dimensioni, ecc. ...) e in cui ci si aspetti che gli "errori" nella produzione siano distribuiti normalmente. In elaborazioni effettuate in tale ambito, quindi, si dispone di un vettore di questo tipo. I calcoli che svolgo di seguito, comunque, sono solo a titolo di esempio e non rappresentano alcun indicatore o controllo tipico della statistica industriale.
Immaginiamo di avere i dati in forma matriciale: 1000 casi controllati (in riga) per 1000 caratteristiche (in colonna). Da questo dataset vogliamo ricavare una matrice di valori dicotomici, con 1 o 0 a seconda che il valore sia maggiore o uguale a zero (1), oppure negativo (0), rispettivamente.
Anche in questo caso il codice da usare per un ciclo di for è immediato:
N<-1000
mtx<-matrix(rnorm(N*N),N,N)
mtx2<-matrix(0,N,N)
system.time(for (i in 1:N) for (j in 1:N) if(mtx[i,j]>=0) mtx2[i,j]<-1)
user system elapsed
6.20 0.00 6.52
La vettorizzazione si può effettuare mediante il comando ifelse(), che analizza ogni singolo elemento della matrice e esegue un'operazione a seconda che il test sia positivo o negativo (da notare che è differente da if() else):
system.time(mtx2<-ifelse(mtx>=0,1,0))
user system elapsed
0.90 0.05 0.98
Anche qui la riduzione del tempo di elaborazione è molto significativa: 85%.
Riporto infine un esempio di vectorizing computing in cui si fa uso di matrici triangolari, molto utili in questo tipo di programmazione in R.
Consideriamo sempre la matrice mtx di 1000 casi per 1000 caratteristiche. Per ogni singolo caso (vettore riga) vogliamo calcolare il rapporto tra l'(i)-esimo elemento e la somma dei successivi 1000-i elementi (quindi la somma dei casi da (i+1) a 1000). La soluzione con un for necessita lo scorrimento dell'intera matrice in questo semplice modo:
N<-1000
mtx<-matrix(rnorm(N*N),N,N)
mtx2<-matrix(0,N,N) system.time(for (i in 1:N) for (j in 1:(N-1)) mtx2[i,j]<-mtx[i,j]/sum(mtx[i,(j+1):N]))
user system elapsed
37.99 0.17 39.93
Per la trasformazione in codice vettorizzato faccio uso del comando lower.tri(), che crea una matrice triangolare inferiore con valori TRUE e FALSE. Dal relativo help di R: "Returns a matrix of logicals the same size of a given matrix with entries TRUE in the lower or upper triangle".
Questo il codice:
N<-1000
mtx<-matrix(rnorm(N*N),N,N)
mtx2<-matrix(0,N,N)
system.time(mtx2<-mtx/(mtx%*%(1*lower.tri(matrix(0,N,N)))))
user system elapsed
2.55 0.00 2.67
Ho praticamente sostituito lo scorrimento dei dati mediante for, con un prodotto matriciale in cui vi è un'opportuna matrice triangolare inferiore. La riduzione del tempo di elaborazione è del 93%.
Riporto di seguito alcuni esempi di for in R con le relative "vettorizzazioni".
Creo un vettore p composto da 500.000 valori estratti da una normale standardizzata e, da questo, ne derivo il vettore formato dagli stessi 500.000 valori rapportati alla loro somma. Forse la vettorizzazione di questo esempio è veramente banale, ma serve a rendersi conto dei tempi di esecuzione.
Il relativo ciclo di for è semplicissimo:
p<-rnorm(500000)
s<-sum(p)
for (i in 1:500000) p[i]<-p[i]/s
Per la valutazione dei tempi di esecuzione basta osservare il risultato di system.time(),che esegue il codice in parentesi e riporta il tempo in secondi. Nel nostro caso ci limitiamo ad osservare il tempo elapsed (in blu è roportato l'output di R):
system.time(for (i in 1:500000) p[i]<-p[i]/s)
user system elapsed
3.2 0.0 3.2
La cosiddetta vettorizzazione è in tal caso banale (p/s), ma il risultato è notevole:
system.time(p<-p/s)
user system elapsed
0.0 0.02 0.02
In termini assoulti, un risparmio di 3.18 secondi può sembrare inutile, ma se si ragiona in termini relativi si osserva un risparmio del 99% del tempo di esecuzione, quindi molto importante per calcoli più complessi!
Prendiamo in esame un esempio in cui la vettorizzazione è meno immediata.
Consideriamo un insieme di 1.000.000 di dati estratti da una normale standardizzata. Pensate ad esempio ad un tipica applicazione in ambito industriale, in cui si controlla che il prodotto rispetti determinate caratteristiche (peso, dimensioni, ecc. ...) e in cui ci si aspetti che gli "errori" nella produzione siano distribuiti normalmente. In elaborazioni effettuate in tale ambito, quindi, si dispone di un vettore di questo tipo. I calcoli che svolgo di seguito, comunque, sono solo a titolo di esempio e non rappresentano alcun indicatore o controllo tipico della statistica industriale.
Immaginiamo di avere i dati in forma matriciale: 1000 casi controllati (in riga) per 1000 caratteristiche (in colonna). Da questo dataset vogliamo ricavare una matrice di valori dicotomici, con 1 o 0 a seconda che il valore sia maggiore o uguale a zero (1), oppure negativo (0), rispettivamente.
Anche in questo caso il codice da usare per un ciclo di for è immediato:
N<-1000
mtx<-matrix(rnorm(N*N),N,N)
mtx2<-matrix(0,N,N)
user system elapsed
6.20 0.00 6.52
La vettorizzazione si può effettuare mediante il comando ifelse(), che analizza ogni singolo elemento della matrice e esegue un'operazione a seconda che il test sia positivo o negativo (da notare che è differente da if() else):
system.time(mtx2<-ifelse(mtx>=0,1,0))
user system elapsed
0.90 0.05 0.98
Anche qui la riduzione del tempo di elaborazione è molto significativa: 85%.
Riporto infine un esempio di vectorizing computing in cui si fa uso di matrici triangolari, molto utili in questo tipo di programmazione in R.
Consideriamo sempre la matrice mtx di 1000 casi per 1000 caratteristiche. Per ogni singolo caso (vettore riga) vogliamo calcolare il rapporto tra l'(i)-esimo elemento e la somma dei successivi 1000-i elementi (quindi la somma dei casi da (i+1) a 1000). La soluzione con un for necessita lo scorrimento dell'intera matrice in questo semplice modo:
N<-1000
mtx<-matrix(rnorm(N*N),N,N)
mtx2<-matrix(0,N,N) system.time(for (i in 1:N) for (j in 1:(N-1)) mtx2[i,j]<-mtx[i,j]/sum(mtx[i,(j+1):N]))
user system elapsed
37.99 0.17 39.93
Per la trasformazione in codice vettorizzato faccio uso del comando lower.tri(), che crea una matrice triangolare inferiore con valori TRUE e FALSE. Dal relativo help di R: "Returns a matrix of logicals the same size of a given matrix with entries TRUE in the lower or upper triangle".
Questo il codice:
N<-1000
mtx<-matrix(rnorm(N*N),N,N)
mtx2<-matrix(0,N,N)
system.time(mtx2<-mtx/(mtx%*%(1*lower.tri(matrix(0,N,N)))))
user system elapsed
2.55 0.00 2.67
Ho praticamente sostituito lo scorrimento dei dati mediante for, con un prodotto matriciale in cui vi è un'opportuna matrice triangolare inferiore. La riduzione del tempo di elaborazione è del 93%.
