Пример ускорения расчётов в R путём многопоточности

Введение

Как следует из Википедии:

R — язык программирования для статистической обработки данных и работы с графикой, а также свободная программная среда вычислений с открытым исходным кодом в рамках проекта GNU.

Данный язык, в настоящее время, нашёл широкое применение во многих практических и чисто научных областях. Однако, исторически сложилось, что скорость многих ресурсоёмких вычислений оставляет желать лучшего. Тема параллельных вычислений в R на habrahabr уже поднималась. В этой статье я попытаюсь показать применение подобного подхода на конкретном примере с использованием пакета для многопоточных вычислений — parallel.

Постановка задачи

Потребность в ускорении вычислений возникла при изучении связи генетических маркеров с фенотипическими признаками человека с помощью пакета для построения нейронных сетей — neuralnet. Однако предложенное решение может быть использовано и для других задач.

Материалы и методы

В работе использовался R 2.15.2, neuralnet — 1.32, parallel — 2.15.2, rbenchmark — 1.0.0. Компьютер под управлением Windows 7, процессором Intell Core i5-2550K CPU @ 3.40GHz (4 ядра) и 8Гб оперативной памяти. Встроенный набор данных data(infert, package=«datasets»).

Результаты и обсуждение

Для начала посмотрим нагрузку на процессор при обычном вычислении нейронной сети. Пример взят из описания пакета neuralnet с небольшой модификацией: threshold=0.0001 — один из критериев остановки, и rep=20 — 20 повторов вычисления, для того, чтобы расчёты были более «сложные».

nn<-function() {   print("Обычное выполнение")   neuralnet(case~parity+induced+spontaneous, infert,             err.fct="ce", linear.output=FALSE, likelihood=TRUE,rep=20, threshold=0.0001) } nn() 

Из Диспетчера задач Windows видно — при выполнении используется только одно из ядер процессора и загрузка ЦП составляет 25%.
Теперь измерим время выполнения этого кода. Для этого я использовал пакет rbenchmark. В общем виде код для измерения будет выглядеть так:

within(benchmark(test.name=test.function(), # название теста и тестируемая функция                  replications=c(3), # число повторов в тесте (+1 повтор на "разогрев")                  columns=c('test', 'replications', 'elapsed'), # стандартные колонки в таблице вывода результатов                  order=c('elapsed', 'test')), # параметры сортировки вывода { average = elapsed/replications }) # нестандартная колонка среднего времени выполнения одного повтора теста 

В результате получается:

                test replications elapsed     average 1 Обычное_выполнение            3   47.83 15.94333333 

Теперь сравним это время с выполнением 20 вычислений не через опцию rep=20, а через, например, sapply.

nn.s<-function() {   print("Последовательное выполнение")   nets<-sapply(1:20, function(X)     neuralnet(case~parity+induced+spontaneous, infert,               err.fct="ce", linear.output=FALSE, likelihood=TRUE,rep=1, threshold=0.0001)) } 

Получается:

                         test replications elapsed average 1          Обычное_выполнение            3   46.05   15.35 2 Последовательное_выполнение            3   47.52   15.84 

Время выполнения примерно одинаковое. Видимо никакой оптимизации выполнения повторов внутри функции neuralnet нет. Чтобы задействовать остальные ядра используем пакет parallel. Он входит в состав R c версии 2.14.0 и позволяет выполнять код в нескольких независимых потоках. Каждый из потоков будет выполнять функцию neuralnet.

nn.p<-function() {   print("Параллельное выполнение")   cl <- makeCluster(getOption("cl.cores", 4)) # создание кластера из четырёх ядер процессора   clusterExport(cl,"infert") # передача данных внутрь кластера   clusterEvalQ(cl,library(neuralnet)) # загрузка пакета neuralnet в кластере   parSapply(cl, 1:20, function(X) # параллельная версия sapply     neuralnet(case~parity+induced+spontaneous, infert,               err.fct="ce", linear.output=FALSE, likelihood=TRUE,rep=1, threshold=0.0001)   )   stopCluster(cl) } 

Сравним время выполнения:

                         test replications elapsed      average 3     Параллельное_выполнение            3   17.38  5.793333333 2 Последовательное_выполнение            3   45.88 15.293333333 1          Обычное_выполнение            3   46.61 15.536666667 

Таким образом, параллельный вариант выполняется более чем в 2,5 раза быстрее. При этом загрузка процессора составляет 100%.

Для удобства можно создать свою обёртку для функции neuralnet.

pneuralnet <- function(formula, data, rep=1, ..., cl) {   clusterExport(cl,"data")   clusterEvalQ(cl,library(neuralnet))   nets <- parLapply(cl, 1:rep, function(X)     neuralnet(formula, data, rep=1, ...)   )   # сортировка списка сетей в зависимости от reached.threshold   nets <- nets[order(sapply(1:rep,function(i){nets[[i]]$result.matrix["reached.threshold", ]}))]   return(nets) }  cl <- makeCluster(getOption("cl.cores", 4)) nets <- pneuralnet(case~parity+induced+spontaneous, infert,                    err.fct="ce", linear.output=FALSE, likelihood=TRUE,rep=4,                    threshold=0.0001, cl=cl) stopCluster(cl) 

library(parallel) library(neuralnet) library(rbenchmark)  data(infert, package="datasets")  nn<-function() {   print("Обычное выполнение")     neuralnet(case~parity+induced+spontaneous, infert,              err.fct="ce", linear.output=FALSE, likelihood=TRUE,rep=20, threshold=0.0001) }  nn.s<-function() {   print("Последовательное выполнение")   nets<-sapply(1:20, function(X)       neuralnet(case~parity+induced+spontaneous, infert,                err.fct="ce", linear.output=FALSE, likelihood=TRUE,rep=1, threshold=0.0001)) }  nn.p<-function() {   print("Параллельное выполнение")   cl <- makeCluster(getOption("cl.cores", 4))   clusterExport(cl,"infert")   clusterEvalQ(cl,library(neuralnet))   parSapply(cl, 1:20, function(X)             neuralnet(case~parity+induced+spontaneous, infert,                err.fct="ce", linear.output=FALSE, likelihood=TRUE,rep=1, threshold=0.0001)                )   stopCluster(cl)   }  pneuralnet <- function(formula, data, rep=1, ..., cl) {     clusterExport(cl,"data")   clusterEvalQ(cl,library(neuralnet))   nets <- parLapply(cl, 1:rep, function(X)             neuralnet(formula, data, rep=1, ...)             )   # сортировка списка сетей в зависимости от reached.threshold   nets <- nets[order(sapply(1:rep,function(i){nets[[i]]$result.matrix["reached.threshold", ]}))]   return(nets) }  within(benchmark(Обычное_выполнение=nn(),                  Последовательное_выполнение=nn.s(),                  Параллельное_выполнение=nn.p(),                  replications=c(3),                  columns=c('Тест', 'replications', 'elapsed'),                  order=c('elapsed', 'test')), { average = elapsed/replications })  cl <- makeCluster(getOption("cl.cores", 4)) nets <- pneuralnet(case~parity+induced+spontaneous, infert,                     err.fct="ce", linear.output=FALSE, likelihood=TRUE,rep=4,                     threshold=0.0001, cl=cl) stopCluster(cl)