- Подписка на печатную версию:
- Подписка на электронную версию:
- Подшивки старых номеров журнала (печатные версии)
LXF130:R
Материал из Linuxformat.
- R Высокопроизводительные вычисления на современном персональном компьютере
Содержание |
R: Ускоряем на примере
| R |
|---|
| R на примере |
|---|
|
- Завершая миникурс под кодовым названием «R на примере», Сергей Петров и Евгений Балдин рассмотрят высокопроизводительные вычисления.
Компьютеры ускоряются, диски становятся просторнее, но темпы роста объёма данных и усложнение их обработки всё равно вынуждают ускорять вычисления.
По традиции, сложные вычисления выполняют кластеры. Многие наверняка слышали про список TOP500, но речь пойдет не о нём: кластерные технологии потихоньку проникают и на обычные компьютеры, можно сказать, в дома. Почти у каждого процессора сейчас несколько ядер. Отчего же этим не воспользоваться? Но сперва прикинем, как оценить выигрыш от подобных приёмов.
Анализ эффективности
Результаты профилирования функции lm()...
...и lm.fit(). Как говорится, почувствуйте разницу!
Измерить скорость вычислений в R и, соответственно, оценить эффективность написанного кода можно несколькими способами:
- Применить system.time() для простых измерений.
- Профилировать выполнение кода с помощью Rprof().
- Профилировать потребление памяти с помощью Rprofmem().
Для визуализации накопленных с помощью Rprof() данных можно применять пакеты profr и proftools.
Приложим средства профилирования к простой задаче. Пусть в ходе моделирования многократно вычисляются параметры линейной регрессии, через функцию lm(). Её недостаток – много лишних действий: она просто слишком универсальна. Для простой линейной регрессии достаточно уточнённого вызова lm.fit().
Насколько же эффективнее прямой вызов? Попробуем обработать набор макроэкономических показателей longley с помощью обоих функций. Построим зависимость между числом работающих и остальными показателями. Для оценки проделаем по 1000 вычислений за раз.
# Загружаем данные > data(longley) # Записываем профиль в файл lm.out > Rprof(“lm.out”) # Выполняем lm() 1000 раз > invisible(replicate(1000,lm(Employed ~ .1, data=longley))) # Отключаем профилирование > Rprof(NULL) # Готовим данные для lm.fit() > longleydm < data.matrix(data.frame(longley)) # Записываем профиль в файл lm.fit.out > Rprof(“lm.fit.out”) # Выполняем lm.fit() 1000 раз > invisible(replicate(1000,lm. fit(longleydm[,7],longleydm[,7]))) # Отключаем профилирование > Rprof(NULL)
Записанные в файлах профилирования данные можно проанализировать с помощью встроенной команды summaryRprof(). Например, для того, чтобы выяснить время работы программы, достаточно обратиться к переменной sampling.time:
> summaryRprof(“lm.out”)$sampling.time [1] 6.42 > summaryRprof(“lm.fit.out”)$sampling.time [1] 0.44
Данные профилирования можно отобразить и графически (см. рисунок), с помощью пакета profr:
# Устанавливаем пакет (если нужно) > install.packages(“profr”) > library(“profr”) > plot(parse_rprof(“lm.out”),main=”Profile of lm()”) > plot(parse_rprof(“lm.fit.out”),main=”Profile of lm.fit()”)
Из рисунков видно, сколько времени проводит программа в каждой из функций. Неудивительно, что lm.fit() работает в четырнадцать (!) раз быстрее.
Для представления зависимостей вызовов в виде графа можно воспользоваться пакетом proftools. Для этого необходимо установить системный пакет graphviz-dev:
=> aptitude install graphvizdev
а также пакет Rgraphviz и сам proftools:
> source(“http://bioconductor.org/biocLite.R”) > biocLite(“Rgraphviz”) > install.packages(“proftools”)}
Для отображения графа нужно выполнить такие команды:
> library(“Rgraphviz”) > library(“proftools”) > lmfitprod < readProfileData(“lm.fit.out”) > plotProfileCallGraph(lmfitprod)
Для отладки потребления памяти нужно использовать специально модифицированную версию R, которая собрана с опцией --enable-memory-profiling. При этом для анализа применяется команда Rprofmem, и все происходит аналогично использованию Rprof.
Перерасход дисковой памяти отследит функция tracemem: она срабатывает при каждом копировании какого-либо объекта.
Ключ к ускорению
Достичь максимальных скоростей поможет знание встроенных функций R. Обычно такие функции создаются, когда выполняемые ими задачи начинают частить в обработке данных. Подобный подход существенно сокращает усилия по вводу кода программы, а также ускоряет выполнение кода, за счёт реализованных во встроенных функциях оптимальных алгоритмов.
Как правило, встроенные функции реализованы на более низкоуровневом языке, предоставляющем углубленный доступ к возможностям оборудования. Поясним это на примере: напишем программу сложения чисел от 1 до указанного, в обычном цикле, и сравним с программой, использующей встроенные функции R:
> mysum <- function(N) { a <- 0;
+ for (i in 1:N) a <- a + i;
+ return(a) }
> system.time(mysum(10000000))
user system elapsed
7.456 0.024 7.482
> system.time(sum(as.numeric(seq(1,10000000))))
user system elapsed
0.052 0.060 0.112
Разница очевидна: второй вариант подсчёта суммы работает более чем в 70 раз быстрее! То, что сейчас было проделано, носит название «векторизация». На такой способ расчётов трудно переключиться после императивного стиля программирования, но именно векторизация обуславливает необычайно высокую продуктивность работы в R. Заметьте: не только уменьшилось время выполнения программы, но и сократился её размер!
Возьмем простую тестовую задачу: «Найти распределение детерминанта матрицы 2×2, в которую занесены независимо и случайно изменяющиеся значения из диапазона 0, 1, 2, … , 9». Она эквивалентна задаче о нахождении всех сочетаний ab-cd, где a, b, c, d – это цифры.
Наивное императивное решение «с циклами» выглядит привычно, и благодаря синтаксису R достаточно «компактно»:
Граф вызовов, полученный с помощью пакета proftools.
> dd.for.0 <- function()
+ {
+ val <- NULL
+ for (a in 0:9) for (b in 0:9) for (d in 0:9) for (e in 0:9)
+ val <- c(val, a*b - d*e)
+ table(val)
+ }
> system.time(dd.for.0())
user system elapsed
0.196 0.000 0.195
От запуска к запуску время слегка изменяется, так что распределение отличается от нормального, и лучше всего находить медиану нескольких попыток:
median(replicate(20, system.time(dd.for.0())[“elapsed”])) [1] 0.177
Попробуем добиться от этого наивного варианта большего: например, выделим память под расчёты сразу в начале программы.
> dd.for.1 <- function()
+ {
+ val <- double(10000) # преаллоцируем val
+ nval <- 0
+ for (a in 0:9) for (b in 0:9) for (d in 0:9) for (e in 0:9)
+ val[nval <- nval + 1] <- a*b - d*e
+ table(val)
+ }
> median(replicate(20, system.time(dd.for.1())[“elapsed”]))
[1] 0.059
Эффект от проделанных изменений очевиден: код ускорился более чем в три раза. Поскольку наши данные – целые числа, давайте посмотрим, что даст нам использование встроенной функции tabulate():
> dd.for.3 <- function()
+ {
+ val <- double(10000)
+ nval <- 0
+ for (a in 0:9) for (b in 0:9) for (d in 0:9) for (e in 0:9)
+ val[nval <- nval + 1] <- a*b - d*e
+ tabulate(val)
+ }
> median(replicate(20, system.time(dd.for.3())[“elapsed”]))
[1] 0.057
Время выполнения сократилось еще чуть-чуть, но это всё полумеры. Давайте сделаем решительный шаг и вспомним, что прародителем R был и язык APL. Запишем решение задачи как операцию над массивами:
> dd.fast <- function()
+ {
+ val <- outer(0:9, 0:9, “*”)
+ val <- outer(val, val, “-”)
+ tabulate(val)
+ }
> median(replicate(20, system.time(dd.fast())[“elapsed”]))
[1] 0.001
Лучшее из наших решений, использующих циклы, обойдено более чем в 50 раз, а «традиционное» – почти в 200!
Что же делать, если избавиться от цикла нельзя? Один из вариантов – использовать среду R, в которую встроена возможность компиляции кода. Сборка R с функцией JIT (Just-in-Time Compilation, компиляция во время выполнения) позволяет рассчитывать на ускорение кода, содержащего циклы, раза в полтора.
Операции с матрицами являются в R такими производительными, потому что они опираются на процедуры библиотеки BLAS (Basic Linear Algebra Subprogram). R может быть скомпилирова на с различными вариантами реализации BLAS: это и свободная библиотека Atlas (пакет atlas3‑base), и платная Goto, и библиотеки от двух основных производителей процессоров: Intel (Intel Math Kernel Library, http://software.intel.com/ru-ru/intel-mkl/) и AMD (AMD Core Math Library, http://www.amd.com/acml). Они не только имеют более производительный код, но и автоматически задействуют в вычислениях все имеющиеся ядра процессора персонального компьютера. Дополнительный прирост производительности можно получить, настроив Atlas под конкретно используемый в расчётах процессор.
Помимо библиотеки BLAS, можно «попытать счастья» с экспериментальным пакетом pnmath0 (http://www.stat.uiowa.edu/~luke/R/experimental/) от Люка Тьерни [Luke Tierney]. Он распараллеливает выполнение векторных функций R, используя потоки Pthreads. Эта возможность пока не встроена в R – её придётся добавить самостоятельно. Учтите, что параллельные вычисления активируются только при достаточной длине векторов аргументов.
Если на компьютере установлена видеокарта, которая поддерживает вычисления на своем GPU (пока сюда входят только NVIDIA CUDA и CUBLAS), то при установке пакета gputools появляется возможность выполнять иерархический кластерный анализ, классификацию с обучением (по алгоритму 'SVM) и расчёт коэффициентов корреляции с очень высокой скоростью.
Несмотря на использование высокопроизводительных векторных операций и компиляции в режиме Just-in-Time, бывают моменты, когда на счету каждый такт процессора. В этом случае есть два способа оперативно встроить в расчет, выполненный в среде R, низкоуровневый код императивного языка программирования:
- inline — для простой вставки небольшого фрагмента кода;
- Rcpp – для облегчённого процесса интеграции сложного кода на C++.
Пакет inline предоставляет функцию cfunction(), умеющую автоматически встраивать код, написанный на Fortran, C, C++. Для выполнения следующего простого примера на Fortran, естественно необходимо установить и загрузить сам inline:
# Не забудьте про отсту пы! Fortran есть Fortran. > code < “ + do i = 1, n(1) + x(i) = x(i)**3 + end do” > cubefn < cfunction(signature(n=”integer”, x=”numeric”), + code, convention=”.Fortran”) > x < as.numeric (1:10) > n < as.integer(10) > cubefn(n,x)$x [1] 1 8 27 64 125 216 343 512 729 1000
В параллель
Среда R, работающая в 64-битном окружении, практически не имеет ограничений на объём обрабатываемых данных. Современным ответом в области вычисления с гигантскими массивами данных является пакет iterators от REvolution Computing. С его возможностью поэлементно обработать структуру, не вмещающуюся в память, крайне удачно сочетается второй пакет, foreach. Он вводит возможность циклически обработать созданный итератор и вернуть суммарный результат. Отсутствие побочных эффектов позволяет выполнить оптимизируемую операцию параллельно.
Отсутствие побочных эффектов – именно то, что позволяет не заботиться, где выполняется та или иная часть кода. Компьютеры не только стали мощнее и у них теперь больше ядер. Компьютеров, прежде всего, стало много больше. Под рукой практически у каждого из нас имеются 5–10 машин, и их многочисленные процессоры, если приглядеться внимательно, не загружены и на 10%. Вся эта мощь доступна пользователю при работе с R.
Кластерные вычисления настолько естественны для векторных операций R, что существует несколько способов их реализации в этой среде:
- Rmpi – реализация Message Passing Interface (MPI), который является стандартом в области паралельных вычислений;
- NWS – написанная на Python альтернативная реализация MPI;
- snow – высокоуровневая надстройка над MPI, PVM, NWS и сокетами (sockets);
- papply – параллелизация функции apply через MPI;
- multicore – параллельные вычисления на многоядерных машинах.
Поскольку наша цель – это быстро и «глобально» для нашей же пользы задействовать мощь окружающих нас и фактически простаивающих компьютеров, сразу же воспользуемся высокоуровневым средством, а именно, пакетом snow. А чтобы рассказ не был пустым теоретизированием, попробуем решить реальную задачу, имеющую практическую ценность.
Задача
Число телефонных пар, проходящих рядом и передающих сигнал без помех, ограничено. Безусловно, влияет и длина кабеля, и его ёмкость, и диаметр жил. Однако в номограммах оценки характеристик телефонной пары не предусмотрено свыше 25 ADSL-пар в одном кабеле.
Нормативные документы российских провайдеров ADSL определяют емкость одиночного кабеля в 18 %. Скорее всего, для кабелей небольшой ёмкости всё же верен предел на 18 жил, занятых одновременно работающими ADSL-модемами, в одном кабеле. А для кабелей ёмкостью более 100 пар верен процентный лимит.
Пусть в городе Н-ске ёмкость ADSL-сети Н-ск-телеком превысила 10 000 абонентов. Попробуем оценить, какие проблемы встретит сеть при своем развитии.
Модель
Город Н-ск насчитывает 300 000 населения, в средней семье 3 человека, что определяет следующие исходные условия:
# количество абонентов услуги > n.abonentov <- 10000 # количество домов > n.domov <- 1000 # количество квартир в доме > n.kvartir <- 100 # критическое число абонентов для ADSL в одном доме > n.kritic <- 18
Мы получаем модель города в виде вектора, каждый элемент которого – квартира, помеченная номером дома:
> gorod <- rep(1:n.domov, each = n.kvartir)
Поскольку вычисления ресурсоёмкие, сразу загрузим библиотеку для параллелизации основных функций пакета (предполагается, что эта экспериментальная библиотека уже установлена):
> library(“pnmath0”)
Делаем выборку случайных n.abonentov в векторе gorod.
> vyborka <- as.factor(gorod)[sample(1:length(gorod), + n.abonentov, replace= FALSE)]
Подсчитываем, сколько абонентов попало в нее в каждом доме, и строим гистограмму.
> hist(as.numeric(tapply(rep(1,n.abonentov), vyborka, sum)))
Вычислим, сколько домов испытывают трудности:
> length(as.numeric(tapply(rep(1,n.abonentov), + vyborka, sum))[as.numeric(tapply(rep(1,n.abonentov), + vyborka, sum))>n.kritic]) [1] 4
а также оценим число «проблемных» абонентов:
> sum(as.numeric(tapply(rep(1,n.abonentov), + vyborka, sum))[as.numeric(tapply(rep(1,n.abonentov), + vyborka, sum))>n.kritic]) [1] 81
Это только одна реализация. Построим бутстреп-процедуру (LXF128) и оценим долю домов, где будет свыше n.kritic абонентов. Тут нужно не менее 10 000 вычислительных экспериментов:
> nn <- 10000
> bstr.dom <- numeric(nn)
> bstr.abonent <- numeric(nn)
> for (n in 1:nn) {
+ vyborka <- as.factor(gorod)[sample(1:length(gorod), n.abonentov, replace= FALSE)]
+ rr <- as.numeric(tapply(rep(1,n.abonentov),vyborka , sum))
+ bstr.abonent[n] <- sum(rr[rr>n.kritic])
+ bstr.dom[n] <- length(rr[rr>n.kritic])
+ }
# Проблемные абоненты
> hist(bstr.abonent)
# Проблемные дома
> hist(bstr.dom)
Проведём эксперимент с подсчётом: сколько процентов абонентов и какое количество домов окажется с плохим качеством услуги при росте абонентской базы от 10 000 до 20 000 абонентов.
Оформим функцию:
> my.boot.adsl < function (n.abonentov) {
+ nn < 10000
+ bstr.dom < numeric(nn)
+ bstr.abonent < numeric(nn)
+ for (n in 1:nn) {
+ vyborka < as.factor(gorod)[sample(1:length(gorod),
+ n.abonentov, replace= FALSE)]
+ rr < as.numeric(tapply(rep(1,n.abonentov),vyborka , sum))
+ bstr.abonent[n] < sum(rr[rr>n.kritic])
+ bstr.dom[n] < length(rr[rr>n.kritic])
+ }
+ return(abonent=bstr.abonent, dom=bstr.dom)
+}
Рассчитаем оценку числа проблемных абонентов в диапазоне размера абонентской базы от 10 000 до 20 000:
> xx < c(NA)
> for (n in 1:10) {
+ xx[n] < my.boot.adsl(10000+(n*1000))
+}
Отнормируем на размер абонентской базы:
> xx.norm < xx > xx.norm1 < xx1/11000 > for (n in 2:10) { > xx.normn < xxn/(10000+(n*1000)) >}
Наконец, отобразим в виде боксплота:
> boxplot(xx.norm,names=seq(11000,20000,by = 1000))
Данный модельный город предполагал наличие только 100-квартирных домов. Повторим вычисления на данных о реальном количестве квартир в городе Н-ске.
Реальность
Ввиду ресурсоёмкости вычислений напишем параллельную версию программы. Для этого загрузим библиотеку snow:
> library(snow)
Создадим кластер из двух узлов:
> cl <- makeCluster(c(“localhost”,”localhost”), type = “SOCK”)
Критическое число абонентов для ADSL в одном доме –
> n.kritic <- 18 > clusterExport(cl, “n.kritic”)
Получаем модель города в виде вектора, каждый элемент которого – квартира, помеченная номером дома. Обработаем реальные дома города Гродно (Беларусь). Загрузим список максимальных номеров квартир в доме:
> Nsk <- na.omit(read.table(“data_Nsk.txt”)) > Nsk.sort <- sort(t(Nsk))
Файл data_Nsk.txt – просто колонка чисел, которую можно скачать, например, отсюда: http://www.inp.nsk.su/~baldin/data_Nsk.txt.
Создадим модель города:
> gorod <- unlist(mapply(rep, + 1:length(Nsk.sort),Nsk.sort)) > clusterExport(cl, “gorod”)
Отведём место под результаты оценки числа проблемных домов и абонентов, для уменьшения накладных расходов при выделении памяти.
> nn <- 10000 > bstr.dom <- numeric(nn) > clusterExport(cl,”bstr.dom”) > bstr.abonent <- numeric(nn) > clusterExport(cl,”bstr.abonent”)
Функция расчёта числа проблемных абонентов при случайном распределении выглядит так:
> my.boot.func < function (n, n.abonentov) {
+ vyborka < as.factor(gorod)[sample(1:length(gorod),
+ n.abonentov, replace= FALSE)]
+ rr < as.numeric(tapply(rep(1,n.abonentov),vyborka , sum))
+ bstr.abonent[n] < sum(na.omit(rr[rr>n.kritic]))
+}
Рассчитаем оценку числа проблемных абонентов в диапазоне абонентской базы от 3000 до 32 000 человек
> xx < cbind(sapply(1000+((1:10)*3000), + function (n.abonentov) parSapply(cl, + c(1:nn), my.boot.func, n.abonentov )))
Вновь отнормируем на размер абонентской базы:
> xx.norm < xx
> xx.norm[,1] < xx[,1]/4000
> for (n in 2:10) {
+ xx.norm[,n] < xx[,n]/(1000+(n*3000))
+}
И, наконец, отобразим в виде боксплота, после чего остановим кластер.
> boxplot(xx.norm, names=seq(4000,31000,by = 3000)) > stopCluster(cl)
А в чём выигрыш?
По приведённым выше вычислениям тестирование выполнялось на процессоре семейства Intel Core Duo, модель T2050 с частотой 1,60 ГГц. При использовании распараллеливания время вычисления составляло 6470 секунд, а без оного – 12754 секунд. Иными словами, два ядра примерно в два раза лучше, чем одно. Что и Требовалось Доказать.
