Производительность кода
Время выполнения кода
При написании сложных функций и нетривиальном процессинге данных необходимо контролировать время выполнения кода и корректировать те операции, которые занимают слишком много времени или памяти. Для анализа времени выполнения кода в R можно использовать как базовую функцию system.time(), так и функции более подходящих пакетов bench или microbenchmark.
Функция system.time() является обёрткой над двумя вызовами proc.time() и возвращает, сколько времени заняло выполнение выражения. user time (пользователь) — сколько времени потребовалось для выполнения всего выражения (в том числе парсинг выражения и запуск процесса), system time (система) — сколько потребовалось времени процессора для выполнения процесса операционной системой.
set.seed(123)
system.time(for(i in 1:1000) log(runif(1000)))## user system elapsed
## 0.043 0.001 0.042Несмотря на простоту, system.time() не лучший инструмент для измерений: как минимум, могут быть параллельные процессы, более приоритетные для системы. При единичном измерении легко ошибиться, и даже аргумент gcFirst = TRUE, вызывающий сборщик мусора при выполнении, неконтролируемо полезен. В качестве альтернативы, лучше пользоваться функциями microbenchmark::microbenchmark() и bench::mark(), которые выполняют измерение заданное количество раз и в выводе возвращают опиcательные статистики по времени выполнения выражения.
Функция microbenchmark() одноимённого пакета измеряет время выполнения выражения в заданном количестве испытаний и возвращает квартили распределения времени выполнения. Обычно бенчмарки используют для сравнения нескольких разных выражений (функций). Неоптимальный цикл в функции f1() задан умышленно:
library(microbenchmark)
f1 <- function(n) {
res <- vector(mode = 'numeric')
for (i in 1:n) res <- c(res, log(runif(1)))
res
}
f2 <- function(n)
sapply(seq_len(n), function(x) log(runif(1)))
microbenchmark(f1(1000), f2(1000), times = 100)## Unit: milliseconds
## expr min lq mean median uq max neval cld
## f1(1000) 2.921679 3.086750 3.841297 3.188077 3.348235 8.404629 100 b
## f2(1000) 2.193340 2.313839 2.809881 2.382432 2.526168 7.831445 100 aФункция mark() пакета bench возвращает медианное время выполнения выражения в ходе итераций, а также сопутствующую информацию. При этом у пакета есть одно не очень явное ограничение, что-то вроде “защиты от дурака”: нельзя сравнивать блоки кода, которые возвращают разный (например, по классу объектов) результат. mark() учитывает и позволяет выводить результаты измерений с разными параметрами: и относительные оценки (коэффициенты относительно самого быстрого выражения), и надо ли учитывать работу сборщика мусора, и сколько места в памяти потребовалось для измерения.
## # A tibble: 1 × 6
## expression min median `itr/sec` mem_alloc `gc/sec`
## <bch:expr> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl>
## 1 f1(100) 211µs 316µs 3316. 292KB 33.5Расширение функции mark(), функция press(), позволяет использовать mark() в указанном диапазоне параметров. То есть вместо цикла с перебором количества итераций можно просто воспользоваться press():
## Running with:
## rep## 1 1## 2 10## 3 100## # A tibble: 3 × 7
## expression rep min median `itr/sec` mem_alloc `gc/sec`
## <bch:expr> <dbl> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl>
## 1 f1(100) 1 204µs 219µs 4273. 292KB 30.8
## 2 f1(100) 10 203µs 234µs 3570. 292KB 25.9
## 3 f1(100) 100 202µs 217µs 4441. 292KB 31.1Профилирование кода
Оценка производительности кода обычно лишь первый этап — бывает так, что функция или блок кода очень большие и включают в себя множество разнообразных операций. Соответственно, определить, какая из этих операций требует много времени для выполнения, без дополнительного инструментария достаточно сложно.
В базовом R есть функция Rprof(), с помощью которой можно записать результаты профилирования в файл, чтобы потом с помощью summaryRprof() его прочитать и проинтерпретировать. Тем не менее, ощутимо удобнее использовать пакет profvis, который также использует Rprof()/summaryRprof(), но результаты представляет в удобном интерактивном виде.
Посмотрим, какие операции требуют больше всего времени в функции подсчёта, сколько раз встречается та или иная цифра в заданном векторе. Для усложнения возьмем длинный вектор с чисткой прочих символов, кроме цифр.
# определяем функцию
numbers <- function(x) {
x <- as.character(x)
x <- gsub('[^0-9]', '', x)
x <- strsplit(x, '')
x <- unlist(x)
digit_vector <- factor(x, levels = 0:9)
res <- table(digit_vector)
as.numeric(res)
}
# пример работы функции
numbers(c('097-8', '97', '2-25-0'))## [1] 2 0 2 0 0 1 0 2 1 2Для того чтобы отследить время выполнения каждого выражения функции, передадим в profvis() не саму функцию numbers(), а её тело (при работе с большими скриптами можно просто указывать source(your_script.R)).
library(profvis)
vec <- sample(10000, 1e5, replace = TRUE)
profvis(
expr = {
x <- as.character(vec)
x <- gsub('[^0-9]', '', x)
x <- strsplit(x, '')
x <- unlist(x)
digit_vector <- factor(x, levels = 0:9)
res <- table(digit_vector)
as.numeric(res)
},
height = 500
)Rprof() и, следовательно, profvis() не напрямую измеряют время выполнение кода, а с определённым интервалом записывают состояние интерпретатора (для profvis() это 10 мс по умолчанию), какое выражение обрабатывается в этом интервале. Как следствие, для профилирования очень быстрых выражений потребуется определённым образом замедлить выполнение, например, увеличить длину вектора, как в примере выше. Без такого замедления можно получить ошибку Error in parse_rprof(prof_output, expr_source) : No parsing data available. Maybe your function was too fast?.
Оптимизация кода
R нередко называется медленным языком. В этом есть некоторая доля истины, однако в немалой части случаев неторопливость вызвана низким качеством кода и неэффективными конструкциями. К наиболее частым ошибкам или причинам невысокой скорости работы кода можно отнести игнорирование векторизации и ряда других особенностей языка, неоптимальное использование циклов, игнорирование специализированных пакетов и так далее. В сообществе пользователей R нередко даже встречается фраза “пишите на R, а не на %другой язык% средствами R”.
Оптимизация циклов
Очень часто циклы в R пишут неоптимально, вплоть до того, что установка “никогда не пользуйтесь циклами” известна даже начинающим пользователям. Меж тем, это не совсем верно: как правило, плохой цикл начинается с создания пустого объекта и потом последовательного добавления в цикле к нему новых элементов.
Например, у нас есть строковый вектор, каждый элемент которого мы хотим разделить по какому-то символу с помощью strsplit() и извлечь вторую часть каждого элемента.
# создаём строковый вектор
N <- 1e4
set.seed(12345)
uids <- paste0('uid_', sample(N), sample(letters, N, TRUE))
uids[1:3]## [1] "uid_8243b" "uid_720t" "uid_8922g"
# для примера из первого элемента извлекаем второе значение
uids_splitted <- strsplit(uids, '_')
uids_splitted[[1]][2]## [1] "8243b"Неоптимальный цикл будет выглядеть следующим образом:
uids_second <- c()
for (i in seq_len(N))
uids_second <- c(uids_second, uids_splitted[[i]][2])
uids_second[1:3]## [1] "8243b" "720t" "8922g"Такой цикл плох тем, что на каждой итерации создается новый объект uids_second, в который записывается старый объект uids_second и присоединённое к нему на этой итерации новое значение. Старый объект удаляется сборщиком мусора. При больших датасетах эта операция становится очень дорогой, особенно когда дело касается добавления новой строки в большую таблицу.
uids_second <- c()
for (i in seq_len(3)) {
uids_second <- c(uids_second, uids_splitted[[i]][2])
print(tracemem(uids_second))
}## [1] "<0x55b853185d58>"
## [1] "<0x55b8566547d8>"
## [1] "<0x55b8570d1208>"Корректным в данном случае было бы создание объекта и изменение его во время итераций. В англоязычной литературе по R этот ход обычно называют memory preallocation, что не очень точно, так как в данном случае не происходит реального выделения памяти, а всего лишь избегается копирование и удаление объектов на каждой итерации.
Целевой вектор мы создаём с помощью vector() и изменяем его по индексам. Для примера возьмём только первые три элемента, чтобы не удлинять цикл, и посмотрим, как меняется адрес объекта uids_second в памяти.
uids_second <- vector('character', 3)
for (i in seq_len(3)) {
uids_second[i] <- uids_splitted[[i]][2]
print(tracemem(uids_second))
}## [1] "<0x55b856672a78>"
## [1] "<0x55b856672a78>"
## [1] "<0x55b856672a78>"
uids_second## [1] "8243b" "720t" "8922g"Неявные циклы, к которым относится семейство функций *pply (lapply() в частности), работают практически аналогично циклам с предварительным выделением памяти, только с рядом дополнительных оптимизаций, в первую очередь за счёт вызова С-функции и преаллокации на уровне C.
Сравним скорость трёх вариантов цикла, различие на порядки:
loops_bench <- microbenchmark(
'for + c()' = {
uids_second <- c()
for (i in seq_len(N))
uids_second <- c(uids_second, uids_splitted[[i]][2])
},
'for + vector()' = {
uids_second <- vector('character', N)
for (i in seq_len(N))
uids_second[i] <- uids_splitted[[i]][2]
},
'lapply()' = unlist(lapply(uids_splitted, '[[', 2)),
times = 100
)
print(loops_bench, signif = 3)## Unit: milliseconds
## expr min lq mean median uq max neval cld
## for + c() 324.00 334.00 356.00 340.00 357.00 496.00 100 b
## for + vector() 5.53 5.79 6.11 5.94 6.17 9.74 100 a
## lapply() 3.44 3.61 3.89 3.73 3.89 8.06 100 aВекторизация выражений
Ещё один способ ускорить выполнение выражений — не игнорировать особенности языка. В частности, в R очень активно используется векторизация, когда функция применяется ко всем элементам вектора.
В частности, задачу по извлечению второй части составного строкового идентификатора эффективнее всего решать с помощью векторизованной функции gsub(), которая находит и заменяет записанный регулярными выражениями паттерн на тот или иной набор символов (или же её аналогов из других пакетов, stringi, в частности).
Сравним самый быстрый вариант цикла lapply() и решение задачи с помощью gsub(). Для корректности внесём в lapply() ещё и разбивку строковых элементов вектора на части (собственно, эта функция также векторизована). В предыдущем сравнении циклов эта операция была вынесена и не учитывалась в бенчмарках. Как мы видим, решение задачи с помощью более релевантной языковой конструкции оказывается эффективнее, чем перебор элементов даже в оптимальном цикле.
loops_bench <- microbenchmark(
'lapply()' = unlist(lapply(strsplit(uids, '_'), '[[', 2)),
'gsub()' = gsub('^.*_', '', uids) ,
times = 100
)
print(loops_bench, signif = 3)## Unit: milliseconds
## expr min lq mean median uq max neval cld
## lapply() 12.20 13.00 13.60 13.30 13.60 23.7 100 b
## gsub() 4.66 5.19 5.48 5.37 5.67 11.1 100 aПараллелизация кода
Помимо просто написания хорошего кода с правильными конструкциями, можно распараллеливать выполнение кода — выполнять одинаковые операции в доступных параллельных потоках в зависимости от количества ядер процессора или количества процессоров. При этом некоторые пакеты, в частности, data.table, по умолчанию используют половину ядер (и есть отдельная опция для этого).
Для параллелизации прочих частей кода можно воспользоваться как базовой функцией mclapply() (на Windows использует только одно ядро, адекватная замена — parLapply()), которая обращается к пакету parallel, так и спецпакетами foreach и doParallel, которые служат надстройкой над parallel — первый пакет необходим, собственно, для параллелизации, второй служит бэкендом для него.
## Loading required package: iterators## Loading required package: parallel
# задаём количество проб и готовим датасет
trials <- 100
iris_df <- iris[iris$Species != 'setosa', ]
# пишем функцию оценки lm-коффициентов на семпле
lm_bootstrap <- function(df) {
iris_sample <- df[sample(100, replace = TRUE), ]
lm_fit <- glm(Species ~ Sepal.Length, iris_sample, family = binomial)
coefficients(lm_fit)
}
# указываем, какое количество ядер использовать в параллелизации
cores <- detectCores()
cl <- makeCluster(cores)
registerDoParallel(cl)
getDoParWorkers()## [1] 8
# оцениваем каждое выражение
parallel_bench <- microbenchmark(
'lapply' = lapply(seq_len(trials), function(x) lm_bootstrap(iris_df)),
'mclapply()' = mclapply(seq_len(trials), function(x) lm_bootstrap(iris_df)),
'foreach()' = foreach(icount(trials)) %dopar% lm_bootstrap(iris_df),
times = 10
)
# останавливаем параллелизацию по ядрам
stopCluster(cl)
# смотрим результат
print(parallel_bench, signif = 3)## Unit: milliseconds
## expr min lq mean median uq max neval cld
## lapply 152 169 182 183 199 202 10 b
## mclapply() 121 122 127 125 129 150 10 a
## foreach() 113 121 124 124 127 138 10 aКак мы видим, mclapply(), работает сопоставимо с foreach(). В целом, вполне может получаться в некоторых случаях так, что foreach() будет чуть-чуть медленнее. Это объяснимое поведение, так как при работе foreach() появляются дополнительные расходы на запуск и взаимодействие с кластером. Как следствие, foreach() может быть не очень эффективен для параллелизации некоторых (обычно быстрых) операций.
Снижение накладных расходов
Большинство функций пакетов, помимо основных операций, содержат и дополнительные: проверки типов входных данных, обработку ошибок, проверку на наличие и применение дополнительных аргументов и так далее. Это, конечно же, повышает надёжность кода, однако может увеличивать время выполнения.
Достаточно легко роль таких сопутствующих операций можно увидеть на примере вычисления арифметического среднего. Функция mean() — обобщённый метод (в данном случае функция расчёта среднего) для большого набора классов, для каждого из которых метод может быть реализован по-своему. Соответственно, функция mean() будет сначала определять, с каким классом объекта имеет дело, и потом только вызывать соответствующий метод. Функция mean.default() - это метод по умолчанию, когда конкретная имплементация не найдена или у объекта вообще отсутствует атрибут класса. Соответственно, mean.default() должна быть немного быстрее mean(). К тому же обе функции имеют свои допаргументы (trim, na.rm или ...), поэтому будут медленнее, чем прямолинейное вычисление среднего по вектору как отношение суммы элементов и длины вектора.
vec <- rnorm(1e6)
mean_bench <- microbenchmark(
'mean()' = mean(vec),
'mean.default()' = mean.default(vec),
'sum() / len()' = sum(vec) / length(vec),
times = 100
)
print(mean_bench, signif = 3)## Unit: microseconds
## expr min lq mean median uq max neval cld
## mean() 1850 1980 2100 2080 2150 3000 100 b
## mean.default() 1840 1990 2080 2060 2150 2660 100 b
## sum() / len() 933 976 1040 1010 1060 1550 100 aВставки на С++
Ещё одним эффективным методом оптимизации кода может быть использование вставок на других языках программирования, в частности, на C++ (несмотря на наличие интерфейсов к другим языкам, например, Java или Python). Для вставок C++ используется пакет Rcpp или его аналоги.
Сравним функцию вычисления среднего c помощью функций R и с помощью вставки на C++.
# include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
double mean_Rcpp(NumericVector x) {
int n = x.size();
double total = 0;
for(int i = 0; i < n; ++i) {
total += x[i];
}
return total / n;
}Как мы видим, функция получилась быстрая, но всё равно чуть-чуть медленнее, чем базовая реализация sum() / len() на R.
library(Rcpp)
mean_bench <- microbenchmark(
'mean()' = mean(vec),
'mean.default()' = mean.default(vec),
'sum() / len()' = sum(vec) / length(vec),
'mean_Rcpp()' = mean_Rcpp(vec),
times = 100
)
print(mean_bench, signif = 3)## Unit: microseconds
## expr min lq mean median uq max neval cld
## mean() 1840 1930 2000 2000 2060 2240 100 b
## mean.default() 1830 1940 2000 2000 2070 2200 100 b
## sum() / len() 920 956 996 990 1020 1240 100 a
## mean_Rcpp() 913 935 1030 1030 1080 2380 100 a