Wenn R-Code zu langsam läuft, sind es oft nur ganz bestimmte Stellen, die optimiert werden müssen. Nicht immer ist sofort klar, welche Codezeilen das sind. Daher ist es sinnvoll zu wissen, wie man solche Flaschenhälse (oder „Bremsklötze“) effizient und elegant finden kann. Ein hilfreiches Werkzeug dafür ist das sogenannte Profiling: Das automatisierte Erstellen eines Profils, das auf einen Blick zeigt, wo der R-Code „hängen bleibt“.
Profiling schützt vor vorzeitigem Optimieren
Ein zwar altes, aber weiterhin gültiges Zitat nach Donald Knuth beschreibt eine Falle beim Programmieren:
Premature optimization is the root of all evil (or at least most of it) in programming.
Vorzeitiges Optimieren ist die Quelle allen Übels (zumindest eines großen Teils) beim Programmieren.
Donald Knuth, Computer Programming as an Art (1974)
Es ist daher sinnvoll, zunächst funktionierenden, leicht zu lesenden Code zu schreiben, anstatt von Anfang an jede Zeile optimieren zu wollen. Tut man letzteres, stehen die Chancen gut, sich zu verheddern.
Anwendungsbeispiel: Anpassungslinien in Streudiagrammen mit unterschiedlichen Parameter-Einstellungen
Im Anwendungsbeispiel geht es darum, in eine Reihe von Streudiagrammen Anpassungslinien einzuzeichnen. Dazu möchte ich verschiedene Methoden vergleichen, bei denen ich von unterschiedlichen Laufzeiten ausgehe. Eine lineare Anpassungsgerade sollte schnell eingetragen sein: Es gibt nur einen Steigungsparameter, der über den gesamten Wertebereich konstant bleibt. Anders sieht es bei nichtlinearen Anpassungslinien aus, die je nach Datenlage unterschiedliche Kurvenverläufe modellieren können.
Simulierte Daten: 100 Prädiktoren, die schwach mit y korrelieren
Zunächst simuliere ich Zufallsdaten:
n <- 2e6
set.seed(2020)
ygroup <- data.frame(y = rnorm(mean = 100, sd = 20, n = n),
group = sample(LETTERS[1:16], size = n, replace = TRUE))
simulation <- data.frame(replicate(100, rnorm(mean = 50, sd = 10, n = n) + 0.1 * ygroup$y))
simulation <- cbind(ygroup, simulation)
rm(ygroup)
Die Fallzahl soll n = 2 Mio. (2e6) betragen. In ygroup wird eine normalverteilte Zufallsvariable y erzeugt (Mittelwert = 100, Standardabweichung = 20) sowie eine Gruppenvariable group mit 16 Ausprägungen, einfach mit Buchstaben von A bis P bezeichnet.
Dazu kommen in simulation 100 Zufallsvariablen, die von der replicate()-Funktion mit X1 bis X100 benannt werden. Durch die Formel (0.1 * ygroup$y) ist eine schwache Korrelation mit der abhängigen Variable y sichergestellt.
Mit cbind() werden die beiden Teile zu einem Datensatz zusammengefügt.
Anschließend bestimme ich die X-Variable, die am stärksten mit y korreliert, und nenne sie X:
index <- which.max(cor(simulation$y, simulation[, -(1:2)])) + 2 names(simulation[index]) cor(simulation$y, simulation[index]) simulation$X <- simulation[, index]
Streudiagramme mit Anpassungslinien
Nach dieser Vorbereitung kann ich verschiedene Streudiagramme für die Untergruppen erstellen und dabei jeweils Anpassungslinien einzeichnen, etwa so:
ggplot(simulation, aes(x = X, y = y)) +
geom_jitter(alpha = 0.6, size = 1.2) +
geom_smooth(method = "loess") +
facet_wrap(~ group, nrow = 4) +
labs(title = "loess, span = default = 0.75")
loess ist eine Variante von lowess (locally weighted scatterplot smoother), einer lokal gewichteten Glättung. Dabei wird schrittweise eine Reihe von linearen Regressionsgeraden berechnet, wobei jeweils die Punkte in der näheren Umgebung höher gewichtet werden, sodass nichtlineare Verläufe entstehen. Der Parameter span gibt das Ausmaß der Glättung an.
Um Rechenaufwand bzw. Laufzeiten zu vergleichen, variiere ich den span-Parameter, teste in geom_smooth() auch die auto-Methode, die auf gam zurückgreift (gam = generalized additive model, ein leistungsfähiger Algorithmus für nichtlineare Zusammenhänge, siehe Machine Learning-Algorithmen verstehen: Interaktionseffekte), und nehme ein simples lineares Modell auf (method = „lm“).
Profiling: Wie wirken sich verschiedene Glättungsmethoden auf die Code-Laufzeit aus?
Nun möchte ich durch Profiling herausfinden, wie sich die verschiedenen Methoden für Anpassungslinien auf die Code-Laufzeit auswirken. Meine Hypothese ist, dass das lineare Modell deutlich schneller verarbeitet wird als die flexibleren nichtlinearen Modelle.
Um das Profiling zu aktivieren, habe ich zwei Möglichkeiten. Erstens kann ich den gesamten Codeblock in einen profvis-Aufruf einbetten:
profvis({
# Datensimulation
n <- 2e4
set.seed(2020)
ygroup <- ...
...
ggplot(simulation, aes(x = X, y = y)) +
geom_jitter(alpha = 0.6, size = 1.2) +
geom_smooth(method = "loess") +
facet_wrap(~ group, nrow = 4) +
labs(title = "loess, span = default = 0.75")
...
}, interval = 0.005)
Bei den drei Punkten […] habe ich Code-Abschnitte weggelassen, um vor allem das Wesentliche zu zeigen: Die profvis-Funktion. Das Intervall für die Zeitmessung kann angegeben werden – es wird empfohlen, nicht unter 5 Millisekunden zu gehen, da darunter die Genauigkeit leidet.
Zweitens unterstützt RStudio Profiling bequem über die Menüführung. Ich muss lediglich die Codezeilen markieren (üblicherweise ein nicht zu kurzer Codeblock) und über Profile – Profile Selected Lines aufrufen.
Profiling: Wo ist der Flaschenhals? Ein überraschendes Ergebnis
profvis berichtet sowohl Laufzeiten als auch Speicherbedarf. Es gibt mehrere Ansichten. Im Bild ist der sog. flamegraph dargestellt, der einen schnellen visuellen Überblick gibt. Hier sehen wir anhand der Zeitspalte (Time, rechts) auf einen Blick, wo unser R-Code stockt: Bei der Datensimulation! Der Balken ist skaliert auf den prozentualen Anteil der Laufzeit, der auf die jeweilige Stelle im Code entfällt. Die ggplot-Aufrufe fallen gegenüber der Datensimulation kaum ins Gewicht – es spielt also für die Laufzeit des gesamten Skripts praktisch keine Rolle, für welche Glättungsvariante mit welchen Parameter-Einstellungen ich mich entscheide – die Zeit geht ganz woanders verloren!
[Ich habe die Fallzahlen in n variiert, da die Laufzeiten bei unterschiedlicher Hardware zum Teil deutlich voneinander abweichen.]
Es gibt noch weitere Ansichten: Neben dem flamegraph gibt es einen Reiter Data, der die Ergebnisse des Profilings numerisch darstellt.
Darunter findet sich eine Zeitachse, die den Call Stack (deutsch etwa Aufrufstapel, Stapelspeicher) abbildet und interaktiv per mouse-over detaillierte Informationen zeigt.
Im unserem Beispiel führt das in interne Funktionsaufrufe (replicate etwa ruft sapply und lapply auf). Ich denke es hängt vom Anwendungsfall ab, ob diese Details weiterhelfen. Meine Hypothese wurde bereits mit der groben Übersicht im flamegraph widerlegt.
Profiling: Fazit
Das Profiling mit profvis erweist sich als leistungsfähiges und anwenderfreundliches Werkzeug, das mit sehr wenig Aufwand zeigt, an welcher Stelle sich eine Code-Optimierung lohnen könnte. So muss man keine Zeit damit vergeuden, an Code-Zeilen zu feilen, die für die Gesamtlaufzeit kaum relevant sind. Zudem ist das Profiling einfacher und schneller umzusetzen, als verschiedene Code-Blöcke / Funktionsaufrufe jeweils manuell mit separaten Laufzeitmessungen zu untersuchen, wie es spezialisierte Pakete wie bench oder microbenchmark tun. (Inzwischen bevorzuge ich das neuere bench-Paket gegenüber dem älteren microbenchmark – siehe den Beitrag Doubletten ausschließen in R: unique() und wie man es schneller macht).
Viel Erfolg mit Euren R-Projekten!
Literatur zur effizienten R-Programmierung:
Hadley Wickham: Advanced R, Second Edition
Roger Peng: R Programming for Data Science
Hadleys Buch gibt es kostenlos online zu lesen, ebenso wie das von Roger Peng.