Nested Loop für Bootstrapping verschnellern

Wie erweitere ich R um eigene Funktionen oder Pakete? Welches Paket ist passend für meine Fragestellung?

Moderatoren: EDi, jogo

ehi

Nested Loop für Bootstrapping verschnellern

Beitrag von ehi »

Hallo!
Ich bin ziemlich verzweifelt, da ich eine große Rechenoperation durchführen möchte. Bisher habe ich das so gelöst, dass ich einen verschachtelten for-Loop erstellt habe, jedoch scheint das viel zu langsam zu sein. Ich habe jetzt viel gelesen und bemerkt, dass das Problem bekannt ist. Auch wird immer mal gesage, dass "Vektorisierung" die beste Lösung sei - ich weiß jedoch nicht, wie das umzusetzen ist.
Hier ein Beispielcode (der keinen Sinn macht, weil immer das gleiche Ergebnis in die Matrix geschrieben wird):

Code: Alles auswählen

# Bootstrap Funktion
# Mittelwert bootstrappen
bt.fun <- function(x, d){
  dat <- x[d]
  mean(dat)
}

# Matrix in der die Ergebnisse gespeichert werden 
m  <- matrix(rep(0,5*5), ncol = 5)

# Zeitmessung
before = Sys.time()

# verschachtelter For-Loop
for(i in 1:5){
  for(j in 1:5){
    
    set.seed(1)
    # 100 Bootstrap-Samples, Parallel-Computing wird benutzt mit 4 Kernen
    b      <- boot::boot(mtcars$mpg, bt.fun, R = 100, parallel = "snow", ncpus = 4)
    m[i,j] <- b$t0
    
  }
}

# Zeitmessung2
time1 = Sys.time() - before
Was ich gerne machen würde ist, die Ergebnisse des Bootstrappens in eine Matrix zu schreiben. Bei mir ist jedoch die Zahl der Bootstrap Samples bei 10.000 und die matrix viel größer (2295 Zellen). Mit dem beschriebenen Code ist es nicht machbar die ganze Rechenoperation zu beenden..
Ganz herzlichen Dank im Voraus!
jogo
Beiträge: 2085
Registriert: Fr Okt 07, 2016 8:25 am

Re: Nested Loop für Bootstrapping verschnellern

Beitrag von jogo »

Hallo ehi,

bei mir sähe das so aus:

Code: Alles auswählen

# Bootstrap Funktion
# Mittelwert bootstrappen
bt.fun <- function(x, d) mean(x[d])

set.seed(1)
# 100 Bootstrap-Samples, Parallel-Computing wird benutzt mit 4 Kernen
m <- matrix(replicate(n=25, boot::boot(mtcars$mpg, bt.fun, R = 100, parallel = "snow", ncpus = 4)$t0), 5,5)

#######
library("microbenchmark")
microbenchmark(unit="relative", ...)
Bezüglich des Algorithmus habe ich folgende Veränderung vorgenommen:
in Deinem Code setzt Du vor jedem Aufruf von boot() das random seed (IMHO ist das nicht entsprechend der Grundidee des Bootstrapping); in meinem Code setze ich es nur einmal. Die von mir vorgenommene Verkürzung des Codes ist eher loop-hiding - also ist keine wesentliche Veränderung der Laufzeit zu erwarten.
Bitte mache eine Abschätzung der zeitlichen Komplexität des Problems:
dann weißt Du, ob Du bei Deiner großen Matrix nur einen Kaffee trinken kannst, bis der Computer fertig ist, oder ob Du Deine Mittagspause machen kannst oder den Computer abends anschubsen musst, damit Du am Morgen die Ergebnisse hast.

Gruß, Jörg
p.s.: Hier ist mein kompletter Code:

Code: Alles auswählen

library("microbenchmark")

# Bootstrap Funktion
# Mittelwert bootstrappen
bt.fun <- function(x, d){
  dat <- x[d]
  mean(dat)
}

ori = function() {
  # Matrix in der die Ergebnisse gespeichert werden 
  m  <- matrix(rep(0,5*5), ncol = 5)
  
  # verschachtelter For-Loop
  for(i in 1:5){
    for(j in 1:5){
      
      set.seed(1)
      # 100 Bootstrap-Samples, Parallel-Computing wird benutzt mit 4 Kernen
      b      <- boot::boot(mtcars$mpg, bt.fun, R = 100, parallel = "snow", ncpus = 4)
      m[i,j] <- b$t0
      
    }
  }
}

# Bootstrap Funktion
# Mittelwert bootstrappen
my.bt.fun <- function(x, d) mean(x[d])

my <- function() {
  set.seed(1)
  # 100 Bootstrap-Samples, Parallel-Computing wird benutzt mit 4 Kernen
  m <- matrix(replicate(n=25, boot::boot(mtcars$mpg, my.bt.fun, R = 100, parallel = "snow", ncpus = 4)$t0), 5,5)
}

microbenchmark(times=1, unit="relative", ori=ori(), my=my())
Unter den Bedingungen in meiner Umgebung ist Dein Code 3% langsamer als meiner. Dies bestätigt meine Vermutung bzgl. des loop-hidings
ehi

Re: Nested Loop für Bootstrapping verschnellern

Beitrag von ehi »

Vielen vielen Dank schon einmal!
Ich habe beide Versionen getestet und sie brauchen beide fast gleich lange.
Was ich beobachtet habe ist, dass die Zahl der Bootstrap-Samples kaum einen Unterschied macht. Folgende Zeiten habe ich da gemessen (Anzahl der Bootstrap Samples):
10: 28s
100: 28s
1000: 34s

Was jedoch einen riesen Unterschied macht ist wenn das selbe für eine größere Matrix gemacht wird. Für mehrere Werte von i und j habe ich auch mal die Zeiten gemessen:
3: 11s
5: 28s
7: 57s
--> dementsprechend bei >2000 Zellen wird es kritisch

Um ehrlich zu sein habe ich mir schon Kaffee gekocht und viele Male über meinen Code geschlafen - er läuft schon seit mehreren Wochen..
Gibt es nicht eine Möglichkeit wie man das schneller bekommt? Ich habe gelesen, dass For-Loops eine sehr schlechte Performance haben sollen.
Viele Grüße und nochmals vielen Dank!
jogo
Beiträge: 2085
Registriert: Fr Okt 07, 2016 8:25 am

Re: Nested Loop für Bootstrapping verschnellern

Beitrag von jogo »

ehi hat geschrieben: Do Mai 18, 2017 8:22 pm Vielen vielen Dank schon einmal!
Ich habe beide Versionen getestet und sie brauchen beide fast gleich lange.
Was ich beobachtet habe ist, dass die Zahl der Bootstrap-Samples kaum einen Unterschied macht. Folgende Zeiten habe ich da gemessen (Anzahl der Bootstrap Samples):
10: 28s
100: 28s
1000: 34s

Was jedoch einen riesen Unterschied macht ist wenn das selbe für eine größere Matrix gemacht wird. Für mehrere Werte von i und j habe ich auch mal die Zeiten gemessen:
3: 11s
5: 28s
7: 57s
--> dementsprechend bei >2000 Zellen wird es kritisch
als p.s. habe ich eine Messung in meine vorherige Nachricht geschrieben.
Es scheint mit, dass boot::boot() sich parallelisieren lässt - die for-Schleife jedoch erstmal nicht.
Um ehrlich zu sein habe ich mir schon Kaffee gekocht und viele Male über meinen Code geschlafen - er läuft schon seit mehreren Wochen..
Gibt es nicht eine Möglichkeit wie man das schneller bekommt? Ich habe gelesen, dass For-Loops eine sehr schlechte Performance haben sollen.
Das ist so pauschal nicht richtig. Es hängt immer davon ab, ob die Operationen innerhalb der for-Schleife sich vektorisieren lassen oder nicht.
Wenn sich die Operationen nicht parallelisieren lassen, dann bringt die Umwandlung der for-Schleife in andere Konstrukte keinen wesentlichen Geschwindigkeitsgewinn. In Deinem Fall steckt der Hauptaufwand in dem Aufruf von boot::boot() (mit der entsprechenden Anzahl von Replikationen R=...). Wesentlich schöner wäre es, man könne boot::boot() anweisen, gleich alle Werte der Matrix auszurechenen (also in Deinem Minimalbeispiel 25 Werte). Dann kann boot::boot() die Parallelisierung nutzen.

Gruß, Jörg
ehi

Re: Nested Loop für Bootstrapping verschnellern

Beitrag von ehi »

Ah okay - das verstehe ich.

Ich habe mal mit Rprof ein Profiling gemacht, das genau diese Theorie bestätigt und zwar kommt dabei heraus:

Bild

Aber wie kann man das schaffen? Also dass boot:boot gleich alle Werte aus der Matrix berechnet?
Dateianhänge
Bild1.png
Zuletzt geändert von ehi am Do Mai 18, 2017 9:09 pm, insgesamt 1-mal geändert.
jogo
Beiträge: 2085
Registriert: Fr Okt 07, 2016 8:25 am

Re: Nested Loop für Bootstrapping verschnellern

Beitrag von jogo »

Da habe ich leider keine Ahnung.

Da Du für eine 7*7-Matrix ca. 1 Minute benötigst, hast Du bei einer 2000*2000-Matrix

Code: Alles auswählen

> (2000/7)^2/60/24
[1] 56.68934
Tage Zeit. Wenn Du also mit der bisherigen Methode weiter arbeiten möchtest, brauchst Du wesentlich mehr Rechenleistung:
mehr CPUs und oder mehr GHz
.... oder sehr viel Geduld.

Gruß, Jörg
Benutzeravatar
EDi
Beiträge: 1599
Registriert: Sa Okt 08, 2016 3:39 pm

Re: Nested Loop für Bootstrapping verschnellern

Beitrag von EDi »

Wieso t0 rausschreiben (ist doch immer gleich :? ).

Ich verstehe noch nicht ganz was gewollt ist... So ergibt das keinen Sinn für mich.

Hier mal meine Interpretation und einpaar vorschläge (1) boot() kann man auch durch sample() ersetzen (2) mean ist langsam (3) Rcpp ist auch schnell (macht aber bei mean keinen unterschied). danach kann man immernoch versuchen zu parallelisieren.

Code: Alles auswählen

library("microbenchmark")
library(boot)

# originial ---------------------------------------------------------------


# Bootstrap Funktion
# Mittelwert bootstrappen
bt.fun <- function(x, d){
  dat <- x[d]
  mean(dat)
}

ori = function() {
  # Matrix in der die Ergebnisse gespeichert werden 
  m  <- matrix(rep(0,5*5), ncol = 5)
  
  # verschachtelter For-Loop
  for(i in 1:5){
    for(j in 1:5){
      
      set.seed(1)
      # 100 Bootstrap-Samples, Parallel-Computing wird benutzt mit 4 Kernen
      b      <- boot::boot(mtcars$mpg, bt.fun, R = 100, parallel = "snow", ncpus = 4)
      m[i,j] <- b$t0
      
    }
  }
}


# jogo --------------------------------------------------------------------

# Bootstrap Funktion
# Mittelwert bootstrappen
my.bt.fun <- function(x, d) mean(x[d])

my <- function() {
  set.seed(1)
  # 100 Bootstrap-Samples, Parallel-Computing wird benutzt mit 4 Kernen
  m <- matrix(replicate(n=25, 
                        boot::boot(mtcars$mpg, my.bt.fun, R = 100, 
                                   parallel = "snow", ncpus = 4)$t0), 5,5)
}

microbenchmark(times=1, unit="relative", ori=ori(), my=my())


# edi ---------------------------------------------------------------------

# generate bootstrap samples
gen_samples <- function() sample(mtcars$mpg, replace = T)
# calc means of 100 bootstrap samples
boot_means <- function() {
  replicate(100, mean(gen_samples()))
  }
boot_means()
# repeat
myboot <- function() {
  replicate(25, mean(boot_means()))
}
set.seed(123)
myboot()


# edi2 --------------------------------------------------------------------
gen_samples <- function() sample(mtcars$mpg, replace = T)
# calc means of 100 bootstrap samples
boot_means2 <- function() {
  replicate(100, {x <- gen_samples(); sum(x) / length(x)})
}
# repeat
myboot2 <- function() {
  replicate(25, {x <- boot_means2(); sum(x) / length(x)})
}
set.seed(123)
myboot2()

# edi3 --------------------------------------------------------------------

library(Rcpp)
cppFunction('double meanC(NumericVector x) {
  int n = x.size();
  double total = 0;
  
  for(int i = 0; i < n; ++i) {
    total += x[i];
  }
  return Mean(x);
}')

gen_samples <- function() sample(mtcars$mpg, replace = T)
# calc means of 100 bootstrap samples
boot_means3 <- function() {
  replicate(100, meanC(gen_samples()))
}
# repeat
myboot3 <- function() {
  replicate(25, meanC(boot_means3()))
}
set.seed(123)
myboot3()

microbenchmark(times=1, unit="relative", ori=ori(), my=my(), 
               edi = myboot(), edi2 = myboot2(), edi3 = myboot3())

Code: Alles auswählen

Unit: relative
 expr        min         lq       mean     median         uq        max neval
  ori 633.454173 633.454173 633.454173 633.454173 633.454173 633.454173     1
   my 634.870799 634.870799 634.870799 634.870799 634.870799 634.870799     1
  edi   1.332128   1.332128   1.332128   1.332128   1.332128   1.332128     1
 edi2   1.000000   1.000000   1.000000   1.000000   1.000000   1.000000     1
 edi3   1.083489   1.083489   1.083489   1.083489   1.083489   1.083489     1
Bitte immer ein reproduzierbares Minimalbeispiel angeben. Meinungen gehören mir und geben nicht die meines Brötchengebers wieder.

Dieser Beitrag ist lizensiert unter einer CC BY 4.0 Lizenz
Bild.
jogo
Beiträge: 2085
Registriert: Fr Okt 07, 2016 8:25 am

Re: Nested Loop für Bootstrapping verschnellern

Beitrag von jogo »

EDi hat geschrieben: Do Mai 18, 2017 11:56 pm Wieso t0 rausschreiben (ist doch immer gleich :? ).

Ich verstehe noch nicht ganz was gewollt ist... So ergibt das keinen Sinn für mich.
deshalb hielt ich auch meine Änderung am Algorithmus (set.seed() nach außen) für angebracht.
Hier mal meine Interpretation und einpaar vorschläge (1) boot() kann man auch durch sample() ersetzen (2) mean ist langsam (3) Rcpp ist auch schnell (macht aber bei mean keinen unterschied). danach kann man immernoch versuchen zu parallelisieren.
Auch ich mache sonst immer Bootstrapping mit eigener Programmierung entsprechend Deiner Variante (1)
... deshalb habe ich auch keine Ahnung, was boot() so alles anstellt.
Die Geschwindigkeitszuwächse schon bei Variante (1) sind beachtlich (statt in 56 Tagen ist die gleiche Arbeit jetzt in 2,5 Stunden erledigt). Die Ergebnisse des Profiling könnten schon den passenden Hinweis enthalten:
boot() betreibt (wenn ich das richtig verstehe) einen riesigen Aufwand, um das Parallelisieren zu organisieren und zu verwalten (zumindest riesig im Vergleich zu der dann anfallenden winzigen Aufgabe) - das mag in anderen Situatione völlig anders sein. Ob die Ergebnisse mit boot() besser, wenn man von der Parallelisierung keinen Gebrauch macht? Kann man boot() dazu bringen, keine zusätzlichen statistischen Kennzahlen zu berechnen, wenn man nur den Wert $t0 haben möchte?

Gruß, Jörg
p.s.: Danke für die Demonstration von Rcpp: sowas habe ich bisher auch noch nicht gemacht!

Nachtrag:

Code: Alles auswählen

library("microbenchmark")
library(boot)

# originial ---------------------------------------------------------------


# Bootstrap Funktion
# Mittelwert bootstrappen
bt.fun <- function(x, d){
  dat <- x[d]
  mean(dat)
}

ori = function() {
  # Matrix in der die Ergebnisse gespeichert werden 
  m  <- matrix(rep(0,5*5), ncol = 5)
  
  # verschachtelter For-Loop
  for(i in 1:5){
    for(j in 1:5){
      
      set.seed(1)
      # 100 Bootstrap-Samples, Parallel-Computing wird benutzt mit 4 Kernen
      b      <- boot::boot(mtcars$mpg, bt.fun, R = 100, parallel = "snow", ncpus = 4)
      m[i,j] <- b$t0
      
    }
  }
}


# jogo --------------------------------------------------------------------

# Bootstrap Funktion
# Mittelwert bootstrappen
my.bt.fun <- function(x, d) mean(x[d])

jogo <- function() {
  set.seed(1)
  # 100 Bootstrap-Samples, Parallel-Computing wird benutzt mit 4 Kernen
  m <- matrix(replicate(n=25, boot::boot(mtcars$mpg, my.bt.fun, R = 100, parallel = "snow", ncpus = 4)$t0), 5,5)
}

jogo2 <- function() {
  set.seed(1)
  # 100 Bootstrap-Samples, no Parallel
  m <- matrix(replicate(n=25, boot::boot(mtcars$mpg, my.bt.fun, R = 100, parallel="no")$t0), 5,5)
}

jogo3 <- function() {
  set.seed(1)
  # 100 Bootstrap-Samples, no Parallel
  m <- matrix(replicate(n=25, boot::boot(mtcars$mpg, my.bt.fun, R = 100, parallel="no", ncpus = 4)$t0), 5,5)
}


# edi ---------------------------------------------------------------------

# generate bootstrap samples
gen_samples <- function() sample(mtcars$mpg, replace = T)
# calc means of 100 bootstrap samples
boot_means <- function() {
  replicate(100, mean(gen_samples()))
}
boot_means()
# repeat
myboot <- function() {
  replicate(25, mean(boot_means()))
}
set.seed(123)
myboot()


# edi2 --------------------------------------------------------------------
gen_samples <- function() sample(mtcars$mpg, replace = T)
# calc means of 100 bootstrap samples
boot_means2 <- function() {
  replicate(100, {x <- gen_samples(); sum(x) / length(x)})
}
# repeat
myboot2 <- function() {
  replicate(25, {x <- boot_means2(); sum(x) / length(x)})
}
set.seed(123)
myboot2()

microbenchmark(times=1, unit="relative", jogo=jogo(), jogo2=jogo2(), jogo3=jogo3(), edi = myboot(), edi2 = myboot2())
und die Ergebnisse:

Code: Alles auswählen

Unit: relative
  expr        min         lq       mean     median         uq        max neval
  jogo 613.684526 613.684526 613.684526 613.684526 613.684526 613.684526     1
 jogo2   1.000000   1.000000   1.000000   1.000000   1.000000   1.000000     1
 jogo3   1.132280   1.132280   1.132280   1.132280   1.132280   1.132280     1
   edi   1.562846   1.562846   1.562846   1.562846   1.562846   1.562846     1
  edi2   1.526096   1.526096   1.526096   1.526096   1.526096   1.526096     1
ehi

Re: Nested Loop für Bootstrapping verschnellern

Beitrag von ehi »

Aah superIdee! Ganz herzlichen Dank EDi! Das hört sich sehr vielversprechend an und ich denke, dass sich das sehr gut auf mein Problem anwenden lässt!
t0 habe ich nur raus geschrieben, um den Beispiel-Code kurz zu halten - Sinn machts natürlich keinen.

Ich habe übrigens noch eine sehr merkwürdige Sache beobachtet - vllt interessiert das noch jemanden: Und zwar dauert mein langsamer verschachtelter Loop sehr sehr viel länger, wenn ich parallelisiere. Hier Anzahl der Kerne und Dauer für den unten stehenden Code auf meinem Rechner:
1: 0.4s
4: 30.84s

Kann sich das jemand erklären?

Code: Alles auswählen

bt.fun <- function(x, d){
  dat <- x[d]
  mean(dat)
}

ori = function(cores) {
  m  <- matrix(rep(0,5*5), ncol = 5)
  
  for(i in 1:5){
    for(j in 1:5){
      
      set.seed(1)
      b      <- boot::boot(mtcars$mpg, bt.fun, R = 1000, parallel = "snow", ncpus = cores)
      m[i,j] <- b$t0
      
    }
  }
}

t4 <- system.time(ori(cores = 4))
t1 <- system.time(ori(cores = 1))
jogo
Beiträge: 2085
Registriert: Fr Okt 07, 2016 8:25 am

Re: Nested Loop für Bootstrapping verschnellern

Beitrag von jogo »

Hallo ehi,

meine Gedanken dazu habe ich bereits in der vorherigen Nachricht geschrieben. Es geht sogar soweit, dass ..., parallel="no", ncpus = 4) 13% langsamer ist als ..., parallel="no") Ob boot() so lange überlegen muss, was diese Parameterkombination bedeuten soll?

Gruß, Jörg
Antworten