За последнее десятилетие 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)
![dotplot](https://habrastorage.org/getpro/habr/post_images/32b/793/9c0/32b7939c000e77d392d7b1c23a55f521.png)
![boxplot](https://habrastorage.org/getpro/habr/post_images/26c/9b2/d5e/26c9b2d5ef8b61ed2847be5d3b76d4e2.png)
![violin plot](https://habrastorage.org/getpro/habr/post_images/cfb/3b7/742/cfb3b7742f8d8e84f2973a7d2d96bac8.png)
Итоги
Итак, по результатам замеров можно сформулировать следующие выводы.
- Опорное (наивное) решение работает медленно, поскольку двойной цикл
for
производитвысокоуровневых итераций, что не соответствует векторизированной направленности языка R;
- Прямой перенос цикла в функцию
outer
с принудительной векторизацией черезVectorize
не решает проблему, потому что решение семантически идентично наивному; - Рекурсивное решение быстрее на порядок, потому что высокоуровневых итераций становится
а операции с подматрицами в R достаточно быстрые;,
- Использование
outer
в сочетании с естественной векторизацией функцииpmin
на два порядка быстрее за счёт низкоуровневой реализации (скомпилированный C); - Прямое обращение к
Rcpp
на три порядка быстрее (без учёта времени на компиляцию С++), поскольку пакетRcpp
дополнительно оптимизирован при помощи современных возможностей языкаC++.
А теперь самое важное, что я хотел бы донести этой статьёй. Если знать, как устроен язык R и какую философию он исповедует; если помнить про естественную векторизацию и мощный набор функций базового R; наконец, если вооружиться современными средствами замера времени исполнения, то даже самая тяжелая задача может оказаться вполне подъёмной, и вам не придётся проводить долгие минуты и даже часы в ожидании завершения работы вашего кода.
Небольшой disclaimer: пакет benchr
задуман и реализован Артёмом Клевцовым (@artemklevtsov) с добавлениями и улучшениями от меня (Антон Антонов, @tonytonov) и Филиппа Управителева (@konhis). Мы будем рады баг репортам и предложениям по функциональности.
Задачи, подобные этой, легли в основу курса «Основы программирования на R» на платформе Stepik.org. Курс открыт без дедлайнов (его можно проходить в свободном темпе) и включен в программу профессиональной переподготовки «Анализ данных» (СПбАУ РАН, Институт биоинформатики).
Если вы находитесь в Санкт-Петербурге, то будем рады видеть вас на одной из регулярных встреч 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!
Комментарии (11)
nikitazvonarev
30.11.2016 16:34+2Ребята, вот же красивое решение безо всяких циклов. Строит ганкелеву матрицу размера LxK, состоящую из номеров побочных диагоналей, в три строчки безо всяких циклов, Rcpp, Vectorize и прочего:
K <- N - L + 1 i <- rep(1:L, K) j <- rep(1:K, each = L) matrix(i + j - 1, nrow = L)
Естественно, если последнюю строчку перед взятием matrix поменять на что-то типа:
matrix(min(1:N, N:1)[i + j - 1], nrow = L)
то получим решение вашей задачи.tonytonov
30.11.2016 16:46-1Спасибо, хороший вариант. Задача, как я упоминал, отчасти учебная — на ней удобно показывать спектр различных приёмов. А вариантов её решить миллион, и всегда найдётся простор для ещё большей оптимизации. Так, например, заменить первый
rep
наrep.int
— будет ещё быстрее.nikitazvonarev
30.11.2016 19:00+2В любом случае, спасибо за статью! Для новичка это must-read, потому что типичный императивный «сишно-паскальный» подход приводит к крайне медленным и страшным программам, и нужно с самого начала учиться «думать» на R :)
Cekory
30.11.2016 23:57+1Быстрый и безболезненный способ повысить производительность циклов — использовать пакет
compiler
. Он присутствует в R из коробки, с версии R 3.4 обещают включать компиляцию по умолчанию. На моей машине строчкаlibrary(compiler);enableJIT(3)
перед всеми расчетами сразу ускоряетbuild_naive
раза в четыре — после включения компиляции она вычисляется даже чуточку быстрее скомпилированного жеbuild_recursive
.
pro100olga
01.12.2016 15:04Большое спасибо за статью! Я работаю в R, а в университете учила С++, и эти циклы из меня буквально лезут, не могу справиться :)
jzha
01.12.2016 16:26+1Спасибо за статью!
Еще немного оптимизации — можно воспользоваться тем, что матрица является персимметричной и работать с побочными диагоналями. В таком случае сокращается число итераций в цикле и не надо использовать min
IntegerMatrix build_cpp_PerSym(const int n) { IntegerMatrix x(n, n); for (int i=0; i<n; i++) for (int j=0; j<=i; j++) x(j, i-j) = x(n-1-j, n-1-(i-j)) = i+1; return x; }
У меня microbenchmark показывает, что этот метод немного быстрее, чем
build_cpp
.
В R также можно использовать из коробки замечательные cpp библиотеки для работы с линейной алгеброй — Armadillo и Eigen. В Armadillo требуемая матрица генерируется, можно сказать, в одну строчку
umat flipToeplitz(const int n) { return fliplr(toeplitz(linspace<uvec>(n, 1, n))); }
Здесь столбцы матрицы Тёплица от вектора n:1 записываются в обратном порядке.
Этим методом матрицы не очень большого размера (~ 100) получается генерировать примерно за то же время, что и методомbuild_cpp
.tonytonov
01.12.2016 17:49-1Спасибо за идеи — абсолютно согласен. У меня есть опыт работы с обеими библиотеками, и он весьма положительный, так что присоединяюсь к рекомендациям.
Для полноты картины прямые ссылки на пакеты, реализующие соответствующие интерфейсы: RcppArmadillo, RcppEigen.
ikbrain
Статья про оптимизацию кода на R. Вывод статьи: пишите на C++. Ставлю класс!
Шучу. А если серьезно, так оно и есть, и в силу удобства и мощи Rcpp крайне рекомендуется для всех пользователей R. В сложных проектах часто даже «векторные» эрные функции переписывал на C++, и работали они потом на порядок быстрее.