Низкоуровневая оптимизация и измерение производительности кода на R

За последнее десятилетие 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

Это частный случай так называемых ганкелевых матриц. Здесь левый верхний и правый нижний углы всегда содержат единичку, вся побочная диагональ (от левого нижнего до правого верхнего угла) состоит из числа , а все остальные значения расположены ступеньками. В итоге имеем нечто похожее на горный хребет или скат крыши.

Естественно ожидать, что параметр может варьироваться в некотором разумном диапазоне, скажем, от единиц до сотен. Поэтому решение задачи мы будем офомлять в виде отдельной функции, а будет выступать в ней аргументом.

Простейший подсчёт показывает, что нужная нам матрица определяется поэлементно при помощи условия . Если бы мы писали на одном из классических императивных языков (таких, как C или C++), то решение на этом было бы закончено, поскольку оно сводится к элементарному двойному циклу. На R это «наивное решение» выглядит так:

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. Такой подход популярен среди разработчиков пакетов в силу своей чрезвычайной эффективности и простоты, а также потому, что для работы с C++ существует замечательный пакет 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, произойдёт компиляция кода на C++, и по её завершении мы получим функцию build_cpp в глобальном окружении, которую можно вызывать как любую другую определённую нами функцию. Из очевидных минусов такого подхода можно отметить необходимость компиляции (это тоже занимает время), не говоря уже о том, что необходимо владеть языком C++. Разумеется, использование Rcpp не является «серебряной пулей», и написать медленный код на C/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, но с некоторыми ключевыми отличиями. Кратко эти отличия можно описать так:

1. Мы используем единый кроссплатформенный механизм измерения, основанный на таймере из С++11;
2. Пользователю предоставляется больше опций для конфигурирования процесса измерения;
3. Вывод текстовых и графических результатов становится более подробным и наглядным.

Пакет 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)

Итоги

Итак, по результатам замеров можно сформулировать следующие выводы.

1. Опорное (наивное) решение работает медленно, поскольку двойной цикл for производит высокоуровневых итераций, что не соответствует векторизированной направленности языка R;
2. Прямой перенос цикла в функцию outer с принудительной векторизацией через Vectorize не решает проблему, потому что решение семантически идентично наивному;
3. Рекурсивное решение быстрее на порядок, потому что высокоуровневых итераций становится , а операции с подматрицами в R достаточно быстрые;
4. Использование outer в сочетании с естественной векторизацией функции pmin на два порядка быстрее за счёт низкоуровневой реализации (скомпилированный C);
5. Прямое обращение к Rcpp на три порядка быстрее (без учёта времени на компиляцию С++), поскольку пакет Rcpp дополнительно оптимизирован при помощи современных возможностей языка C++.

А теперь самое важное, что я хотел бы донести этой статьёй. Если знать, как устроен язык R и какую философию он исповедует; если помнить про естественную векторизацию и мощный набор функций базового R; наконец, если вооружиться современными средствами замера времени исполнения, то даже самая тяжелая задача может оказаться вполне подъёмной, и вам не придётся проводить долгие минуты и даже часы в ожидании завершения работы вашего кода.

Небольшой disclaimer: пакет benchr задуман и реализован Артёмом Клевцовым (@artemklevtsov) с добавлениями и улучшениями от меня (Антон Антонов, @tonytonov) и Филиппа Управителева (@konhis). Мы будем рады баг репортам и предложениям по функциональности.

Задачи, подобные этой, легли в основу курса «Основы программирования на R» на платформе stepik. Курс открыт без дедлайнов (его можно проходить в свободном темпе) и включен в программу профессиональной переподготовки «Анализ данных» (СПбАУ РАН, Институт биоинформатики).

Если вы находитесь в Санкт-Петербурге, то будем рады видеть вас на одной из регулярных встреч SPb R user group (VK, meetup).

Материалы и дополнительное чтение по теме:

1. Norman Matloff, The Art of R Programming: A Tour of Statistical Software Design
2. Patrick Burns, The R Inferno
3. Hadley Wickham, Advanced R
4. Dirk Eddelbuettel, Seamless R and C++ Integration with Rcpp
5. А. Б. Шипунов, Е. М. Балдин и др. Наглядная статистика. Используем R!

Автор: tonytonov

Источник