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

За останнє десятиліття R пройшов великий шлях: від нішевого (як правило, академічного) інструменту до мейнстрімної «великої десятки» найпопулярніших мов програмування. Такий інтерес викликаний багатьма причинами, серед яких і приналежність до open source, і діяльну коммьюніті, і активно зростаючий сегмент застосування методів machine learning / data mining у різноманітних бізнес-завдання. Приємно бачити, коли один з твоїх улюблених мов впевнено завойовує нові позиції, і коли навіть далекі від професійної розробки користувачі починають цікавитися R. Але тут є, проте, одна велика проблема: мова R написаний, за висловом творців, «статистиками для статистиків». Тому він високорівневий, гнучкий, з величезним і надійним функціоналом «з коробки», не кажучи вже про тисячі додаткових пакетів; і одночасно з цим – досить фривольний по частині семантики, соединящий в собі принципово різні парадигми і володіє широкими можливостями для написання жахливо неефективного коду. А якщо врахувати той факт, що для багатьох останнім часом це перший «справжній» мова програмування, то ситуація стає ще більш загрозливою.
Втім, поліпшення якості і швидкості виконання коду на R – універсальна завдання, з якою постійно стикаються не тільки початківці, але й досвідчені користувачі. Саме тому дуже важливо час від часу повертатися до найбільш простим структурам і основним принципам роботи R, щоб пам'ятати про джерела неэффективностей і способи їх усунення. У цій статті я розповім про одному кейсі, який спочатку зустрівся мені в ході роботи над дослідженням фінансових стратегій, а потім був адаптований мною в якості навчального для онлайн-курсу.
Як не загрузнути в матриці
Опускаючи всі релевантні подробиці, ми можемо сформулювати нашу тестову задачу наступним чином. Для довільного натурального nнеобхідно швидко генерувати матриці розміром n \times nось такого виду (наведено приклад n=5):
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

Це приватний випадок так званих ганкелевых матриць. Тут верхній лівий і правий нижній кути завжди містять одиницю, вся побічна діагональ (від лівого нижнього до верхнього правого кута) складається з числа n, а всі інші значення розташовані сходами. У підсумку маємо щось схоже на гірський хребет або скат даху.
Природно очікувати, що параметр nможе змінюватись в деякому розумному діапазоні, скажімо, від одиниць до сотень. Тому рішення задачі ми будемо офомлять у вигляді окремої функції, а nбуде виступати в ній аргументом.
Простий підрахунок показує, що потрібна нам матриця визначається поелементно за допомогою умови M_{i, j} = \min(i + j - 1, 2n - i - j + 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 – інтерпретується в microsoft мова, кожна ітерація циклу тягне безліч накладних витрат на виконання. Приклади з неоптимальним використанням циклів зустрічаються, мабуть, у кожному підручнику по R. І висновок завжди один: по можливості використовувати так звані векторизированные функції. Вони зазвичай реалізовані на низькому рівні © і гранично оптимізовані.
Так, в нашому випадку, коли елемент матриці є функція від двох індексів ij, добре підійде поєднання функцій
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
}

В цьому випадку ефективність досягається за рахунок того, що кількість високорівневих ітерацій не n^2, як в наївному рішенні, а тільки n(роль циклу на себе бере стек викликів).
Нарешті, в тому випадку, коли оптимізація на рівні базових функцій 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)

dotplot
boxplot
violin plot
Підсумки
Отже, за результатами дослідження можна сформулювати наступні висновки.
  1. Опорне (наївна) рішення працює повільно, оскільки подвійний цикл
    for
    виробляє n^2високорівневих ітерацій, що не відповідає векторизированной спрямованості мови R;
  2. Прямий перенос циклу в функцію
    outer
    з примусовою векторизацией через
    Vectorize
    не вирішує проблему, тому що рішення семантично ідентично наївному;
  3. Рекурсивне рішення швидше на порядок, тому що високорівневих ітерацій стає n, а операції з подматрицами в 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!
Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.