За последнее десятилетие R прошёл большой путь: от нишевого (как правило, академического) инструмента до мейнстримной «большой десятки» самых популярных языков программирования. Такой интерес вызван многими причинами, среди которых и принадлежность к open source, и деятельное коммьюнити, и активно растущий сегмент применения методов machine learning / data mining в разнообразных бизнес-задачах. Приятно видеть, когда один из твоих любимых языков уверенно завоёвывает новые позиции, и когда даже далёкие от профессиональной разработки пользователи начинают интересоваться R. Но здесь есть, однако, одна большая проблема: язык R написан, по выражению создателей, «статистиками для статистиков». Поэтому он высокоуровневый, гибкий, с огромным и надёжным функционалом «из коробки», не говоря уже о тысячах дополнительных пакетов; и одновременно с этим – весьма фривольный по части семантики, соединящий в себе принципиально разные парадигмы и обладающий широкими возможностями для написания чудовищно неэффективного кода. А если учесть тот факт, что для многих в последнее время это первый «настоящий» язык программирования, то ситуация становится ещё более угрожающей.
Впрочем, улучшение качества и скорости исполнения кода на R – универсальная задача, с которой постоянно сталкиваются не только новички, но и опытные пользователи. Именно поэтому очень важно время от времени возвращаться к самым простым структурам и основным принципам работы R, чтобы помнить об источниках неэффективностей и способах их устранения. В этой статье я расскажу об одном кейсе, который сначала встретился мне в ходе работы над исследованием финансовых стратегий, а потом был адаптирован мной в качестве учебного для онлайн-курса.
Как не увязнуть в матрице
Опуская все нерелевантные подробности, мы можем сформулировать нашу тестовую задачу следующим образом. Для произвольного натурального необходимо быстро генерировать матрицы размером
вот такого вида (приведён пример для
):
1 2 3 4 5 2 3 4 5 4 3 4 5 4 3 4 5 4 3 2 5 4 3 2 1
Это частный случай так называемых ганкелевых матриц. Здесь левый верхний и правый нижний углы всегда содержат единичку, вся побочная диагональ (от левого нижнего до правого верхнего угла) состоит из числа ,
Естественно ожидать, что параметр может варьироваться в некотором разумном диапазоне, скажем, от единиц до сотен. Поэтому решение задачи мы будем офомлять в виде отдельной функции, а
будет выступать в ней аргументом.
Простейший подсчёт показывает, что нужная нам матрица определяется поэлементно при помощи условия .
build_naive <- function(n) { m <- matrix(0, n, n) for (i in 1:n) for (j in 1:n) m[i, j] <- min(i + j - 1, 2*n - i - j + 1) m }
Даже несмотря на то, что мы заранее инициализируем матрицу нужного размера (отсутствие первой строки – преаллокации памяти – частая ошибка новичков), такое решение работает весьма медленно. Для знатоков языка такой результат вполне предсказуем: поскольку R – интерпретируемый язык, каждая итерация цикла влечёт множество накладных расходов на исполнение. Примеры с неоптимальным использованием циклов встречаются, пожалуй, в каждом учебнике по R. И вывод всегда один: по возможности использовать так называемые векторизированные функции. Они обычно реализованы на низком уровне (на C) и предельно оптимизированы.
Так, в нашем случае, когда элемент матрицы есть функция от двух индексов и
,
outer и pmin:
build_outer_pmin <- function(n) { outer(1:n, 1:n, function(i, j) pmin(i + j - 1, 2*n - i - j + 1)) }
Функция pmin отличается от min тем, что она оперирует векторами и на входе, и на выходе. Сравните: min(3:0, 1:4) ищет общий минимум (ноль в данном случае), а pmin(3:0, 1:4) вернёт вектор из попарных минимумов: (1, 2, 1, 0). Понятие векторизированной функции непросто сформулировать формально, но общая идея, надеюсь, понятна. Именно такую функцию ожидает в качестве аргумента функция outer. Именно поэтому передать в неё min, как в наивном решении, не получится – попробуйте это сделать и убедитесь, что возникнет ошибка. К слову, если мы всё-таки хотим использовать сочетание outer и min, мы можем выполнить принудительную векторизацию функции min, обернув её в функцию Vectorize:
build_outer_min <- function(n) { outer(1:n, 1:n, Vectorize(function(i, j) min(i + j - 1, 2*n - i - j + 1))) }
Надо помнить при этом, что принудительная векторизация такого рода всегда будет работать медленнее, чем естественная. Однако в некоторых ситуациях без Vectorize не обойтись – например, если определение элемента матрицы содержит условный оператор if.
Ещё один способ избавиться неэффективности наивного решения – каким-то образом убрать вложенность циклов – сводится к тому, чтобы итерировать не по единичным элементам матрицы, а по векторам или подматрицам. Вариантов может быть много, я покажу только один из них, рекурсивный. Мы будем строить матрицу «поэтажно», от центра к краям:
build_recursive <- function(n) { if (n == 0) return(0) m <- matrix(0, n, n) m[-1, -n] <- build4(n - 1) + 1 m[1, ] <- 1:n m[, n] <- n:1 m }
В этом случае эффективность достигается за счёт того, что количество высокоуровневых итераций не ,
(роль цикла на себя берёт стек вызовов).
Наконец, в том случае, когда оптимизация на уровне базовых функций R не приводит к желаемому результату, всегда есть ещё один вариант: реализовать медленную часть на максимально низком уровне. То есть написать кусочек кода на чистом C или С++, скомпилировать его и вызывать из R. Такой подход популярен среди разработчиков пакетов в силу своей чрезвычайной эффективности и простоты, а также потому, что для работы с Rcpp.
library(Rcpp) cppFunction("NumericMatrix build_cpp(const int n) { NumericMatrix x(n, n); for (int i=0; i<n; i++) for (int j=0; j<n; j++) x(i, j) = std::min(i + j + 1, 2*n - i - j - 1); return x;}")
Когда мы подобным образом вызываем cppFunction, произойдёт компиляция кода на build_cpp в глобальном окружении, которую можно вызывать как любую другую определённую нами функцию. Из очевидных минусов такого подхода можно отметить необходимость компиляции (это тоже занимает время), не говоря уже о том, что необходимо владеть языком Rcpp не является «серебряной пулей», и написать медленный код на C/
Пакет microbenchmark
Когда мы рассуждали о быстродействии решений, мы опирались на априорные знания о том, как устроен язык. Теория теорией, но любые рассуждения о производительности должны сопровождаться экспериментами – измерением времени исполнения на конкретных примерах. Более того, в подобных случаях всегда лучше всё проверять и перепроверять самому, потому что любые такие замеры не являются абсолютной истиной: они зависят от железа и состояния процессора, от операционной системы и её нагрузки в момент измерения.
Есть несколько типичных способов измерения скорости исполнения кода на R. Один из наиболее распространённых – это функция system.time. Когда используют эту функцию, обычно делают что-то вроде system.time(replicate(expr, 100)). Это неплохой вариант, но есть несколько пакетов, которые занимаются тем же самым, но в гораздо более удобном и настраиваемом исполнении. Фактически стандартом является пакет под названием microbenchmark, его мы и будем использовать. Всю работу будет выполнять одноимённая функция, которой можно передать произвольное количество выражений (см. unevaluated expression) для измерения. По умолчанию каждое выражение будет выполнено сто раз, и мы увидим собранную статистику по времени исполнения.
library(microbenchmark) measure <- function(n) { microbenchmark(build_naive(n), build_outer_min(n), build_outer_pmin(n), build_recursive(n), build_cpp(n)) } m <- measure(100) m Unit: microseconds expr min lq mean median uq max neval build_naive(n) 20603.310 21825.7975 24380.0332 22884.3255 24215.368 66753.031 100 build_outer_min(n) 22624.873 25430.4985 28164.4496 26662.8955 29028.744 84972.926 100 build_outer_pmin(n) 159.755 187.8325 295.3822 204.1985 225.069 3710.103 100 build_recursive(n) 1741.992 2822.2910 5990.9897 3241.1980 3918.055 55059.974 100 build_cpp(n) 21.321 23.4230 33.6600 34.2335 38.588 91.289 100
В полученном выводе результаты замера представлены в виде таблицы с набором описательных статистик. Кроме того, полученный объект можно рисовать: plot(m) или ggplot2::autoplot(m).
Пакет benchr
В настоящее время в разработке находится ещё один пакет, предназначенный для замера времени исполнения выражений – пакет benchr. Функционал и синтаксис практически идентичен пакету microbenchmark, но с некоторыми ключевыми отличиями. Кратко эти отличия можно описать так:
- Мы используем единый кроссплатформенный механизм измерения, основанный на таймере из
С++11; - Пользователю предоставляется больше опций для конфигурирования процесса измерения;
- Вывод текстовых и графических результатов становится более подробным и наглядным.
Пакет benhcr находится в CRAN (для версии R старше 3.3.2), а development-версия доступна на gitlab, там же есть инструкции по установке. Использование benchr в нашем примере будет выглядеть так (всё, что нужно – заменить функцию microbenchmark на функцию benchmark):
install.packages("benchr") library(benchr) measure2 <- function(n) { benchmark(build_naive(n), build_outer_min(n), build_outer_pmin(n), build_recursive(n), build_cpp(n)) } m2 <- measure2(100) m2 Benchmark summary: Time units : microseconds expr n.eval min lw.qu median mean up.qu max total relative build(n) 100 21800 24300 25600.0 27400.0 28100.0 79300.0 2740000 906.00 build2(n) 100 24400 28800 30100.0 31600.0 33700.0 44000.0 3160000 1070.00 build3(n) 100 157 189 198.0 738.0 217.0 43400.0 73800 7.02 build4(n) 100 1820 1980 3380.0 4910.0 4050.0 50000.0 491000 120.00 build5(n) 100 21 27 28.2 29.4 30.5 50.1 2940 1.00
Помимо удобного текстового вывода (столбец relative позволяет быстро оценить преимущество самого быстрого решения), мы можем воспользоваться различными визуализациями. Кстати, если у вас установлен пакет ggplot2, он будет использован для отрисовки графиков. Если вы по каким-то причинам не хотите использовать ggplot2, такое поведение может быть изменено при помощи опций пакета benchr. Сейчас поддерживаются три типичных визуализации результатов замера: dotplot, boxplot и violin plot.
plot(m2) boxplot(m2) boxplot(m2, violin = TRUE) 
Итоги
Итак, по результатам замеров можно сформулировать следующие выводы.
- Опорное (наивное) решение работает медленно, поскольку двойной цикл
forпроизводитвысокоуровневых итераций, что не соответствует векторизированной направленности языка R;
- Прямой перенос цикла в функцию
outerс принудительной векторизацией черезVectorizeне решает проблему, потому что решение семантически идентично наивному; - Рекурсивное решение быстрее на порядок, потому что высокоуровневых итераций становится
,
- Использование
outerв сочетании с естественной векторизацией функцииpminна два порядка быстрее за счёт низкоуровневой реализации (скомпилированный C); - Прямое обращение к
Rcppна три порядка быстрее (без учёта времени на компиляцию С++), поскольку пакетRcppдополнительно оптимизирован при помощи современных возможностей языкаC++.
А теперь самое важное, что я хотел бы донести этой статьёй. Если знать, как устроен язык R и какую философию он исповедует; если помнить про естественную векторизацию и мощный набор функций базового R; наконец, если вооружиться современными средствами замера времени исполнения, то даже самая тяжелая задача может оказаться вполне подъёмной, и вам не придётся проводить долгие минуты и даже часы в ожидании завершения работы вашего кода.
Небольшой disclaimer: пакет benchr задуман и реализован Артёмом Клевцовым (@artemklevtsov) с добавлениями и улучшениями от меня (Антон Антонов, @tonytonov) и Филиппа Управителева (@konhis). Мы будем рады баг репортам и предложениям по функциональности.
Задачи, подобные этой, легли в основу курса «Основы программирования на R» на платформе stepik. Курс открыт без дедлайнов (его можно проходить в свободном темпе) и включен в программу профессиональной переподготовки «Анализ данных» (СПбАУ РАН, Институт биоинформатики).
Если вы находитесь в Санкт-Петербурге, то будем рады видеть вас на одной из регулярных встреч SPb R user group (VK, meetup).
Материалы и дополнительное чтение по теме:
- Norman Matloff, The Art of R Programming: A Tour of Statistical Software Design
- Patrick Burns, The R Inferno
- Hadley Wickham, Advanced R
- Dirk Eddelbuettel, Seamless R and
C++ Integration with Rcpp - А. Б. Шипунов, Е. М. Балдин и др. Наглядная статистика. Используем R!
ссылка на оригинал статьи https://habrahabr.ru/post/316516/
Добавить комментарий