Как стать автором
Обновить

Комментарии 12

Первое правило оптимизаций — не оптимизировать то, что и так хорошо работает.
Начать нужно с примитивных функций, например с sum.of.squares и лишнего apply.
Вот пример того, что я получаю. Тестирую apply по разным измерениям, затем purrr::map_dbl (самый быстрый) и future::map_dbl (ооооочень медленно).


Reprex
library(tidyverse)
library(microbenchmark)
#> Warning: package 'microbenchmark' was built under R version 4.0.3
library(furrr)
#> Warning: package 'furrr' was built under R version 4.0.3
#> Loading required package: future
#> Warning: package 'future' was built under R version 4.0.3

a <- matrix(rnorm(500000, mean=0, sd=2), 100000, 50)

sum_of_squares <- function(n) {
  n_sq <- n ^ 2
  sum(n_sq)
}

# Verifying summing over different dimensions
abs(
    sum(apply(a, 1, sum_of_squares)) - 
    sum(apply(a, 2, sum_of_squares)))
#> [1] 0

abs(
    sum(map_dbl(1:dim(a)[2], ~sum_of_squares(a[,.x]))) - 
    sum(map_dbl(1:dim(a)[1], ~sum_of_squares(a[.x, ]))))
#> [1] 0

plan(cluster, workers = 4L)
# Verifying cluster setup, should be 4 unique pids
future_map_int(1:20, ~Sys.getpid()) %>% unique
#> [1] 14240  2112  5356 14848

abs(
    sum(future_map_dbl(1:dim(a)[2], ~sum_of_squares(a[,.x]))) -
    sum(future_map_dbl(1:dim(a)[1], ~sum_of_squares(a[.x, ]))))
#> [1] 0

microbenchmark(
    apply_dim_1 = sum(apply(a, 1, sum_of_squares)), 
    apply_dim_2 = sum(apply(a, 2, sum_of_squares)), 
    map_dim_1 = sum(map_dbl(1:dim(a)[1], ~sum_of_squares(a[.x, ]))),
    map_dim_2 = sum(map_dbl(1:dim(a)[2], ~sum_of_squares(a[,.x]))),
    future_dim_1 = sum(future_map_dbl(1:dim(a)[1], ~sum_of_squares(a[.x, ]))),
    future_dim_2 = sum(future_map_dbl(1:dim(a)[2], ~sum_of_squares(a[,.x]))),
    times = 100L)
#> Unit: milliseconds
#>          expr      min       lq      mean   median         uq       max neval
#>   apply_dim_1 404.4141 460.1293 520.03837 510.4109  563.45995  802.6375   100
#>   apply_dim_2 130.5241 150.7894 192.64263 172.8856  234.17170  553.5947   100
#>     map_dim_1 445.9119 476.0938 528.86206 516.6661  549.29145  913.7160   100
#>     map_dim_2  39.1026  53.3201  60.41777  60.3456   68.07625   86.4652   100
#>  future_dim_1 789.9283 877.1860 954.96807 936.8079 1012.73945 1566.2939   100
#>  future_dim_2 601.0447 660.5598 720.77519 692.2518  752.43105 1506.9873   100

Created on 2020-11-09 by the reprex package (v0.3.0)


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


Самое очевидное — нет абсолютно никакого смысла считать вдоль dim = 1, даже если с точки зрения memory layout это верно. map_dim_2 выполняет суммирование всех квадратов вдоль dim = 2 за 60 миллисекунд, вызывая кастомную функцию только 50 раз, что почти в 9 раз быстрее чем вдоль dim = 1, что приводит к вызову той же функции 100_000 раз.


Опять же, если я нигде не допустил ошибки, использование multi-process parallelization в данном случае (на жалких 5 миллионах точек) не оправдано. С другой стороны, если каждый из 50 вызовов суммирующей функции занимал бы, например, 10 секунд, то распараллеливание дало бы ускорение, на вскидку, порядка ~2 раз на моем CPU, и наверное близко к ~10 на CPU автора (в зависимости от оверхеда).

Задача, на которой всё демонстрируется – это не реальный случай, а специально упрощённый пример, придуманный просто для удобства изложения. Обычно на вход приходят не случайные значения с нормальным распределением, а внутри функции сидит какое-то более сложное вычисление, чем просто сумма квадратов. Цель была не в том, чтобы показать, как оптимизировать конкретно этот пример, а чтобы продемонстрировать способ ускорения одинаковых операций над строками за счёт распараллеливания вычислений. При этом хотелось показать всё на максимально понятном примере, чтобы не отвлекаться на объяснения каких-то специфических вычислений.
Предложенная идея с подсчётом по столбцам в данном случае совсем не подходит: так можно поступать, только если мы хотим просуммировать все значения в матрице, но нам нужны отдельные результаты для каждой из строк. Так что никакие способы с перегруппировкой вычислений не походят. Грубо говоря, три велосипеда – это три руля, три рамы и шесть колёс. Но если одному ребёнку дать в подарок три руля, второму – три рамы, а третьему – шесть колёс, это будет для каждого из них совсем не то же, что получить в подарок целый велосипед. Хотя суммы и равны. Допустим, если строки – это показания с набора датчиков на каждый момент времени, то вычисление суммы квадратов показаний с первого датчика за всю историю наблюдений никак не поможет нам получить сумму квадратов показаний всех датчиков в тот или иной момент времени.

А, ну я значит упустил момент про абсолютную наобходимость суммирования по строкам. Если величины настолько разные, почему тогда матрица а не data.frame? Есть какие-то ограничения? Работать с прямоугольными данными проще, если они — data.frame.


В любом случае, я лишь хотел отметить, что просто так добавить future — не панацея.
Например, в вашем случае у вас прирост примерно в 4 раза, хотя логических ядер у вас 12. Т.е., предполагая что availableWorkers() у вас тоже 12, можно еще что-то оптимизировать прежде чем распараллеливать. Возможно, меньшее количество процессов даст больший прирост производительности на данном конкретрном демо-примере.


Отдельный вопрос, полезный для изучения — как future_apply* распределяет работу?
Например, функция может отгружать новые задачи по мере завершения предыдущих, по одной. Альтернативно — сразу разделить весь датасет и назначить каждому процессу свой блок, после чего ожидать завершения всех участвующих процессов. В первом случае большое количество легковесных операций (как в данном примере) затормозит работу, и, возможно, потребует дополнительных оптимизаций.

Основное преимущество таблиц данных (data.frame) — возможность хранить в разных столбцах значения различающихся типов (числовые значения, категориальные переменные, текст и т.п.). По всей видимости, эта возможность не бесплатная, поэтому если она не нужна (например, если мы обрабатываем показания с набора датчиков и это заведомо будут только числа), то использование матрицы может дать некоторый прирост производительности. Опять же, я не занимался специальной оптимизацией конкретно этого примера, он исключительно демонстрационный, но даже в нём переход с матрицы на таблицу данных даёт замедление расчёта большее, чем можно отыграть распараллеливанием.
Уменьшение числа используемых ядер (по умолчанию используются все ядра) выглядит не самой хорошей идеей. Непонятно, каким образом это может повысить производительность. В моём случае наилучшие результаты были с параметрами по умолчанию. А насчёт возможностей тонкой настройки параметров параллельного запуска согласен — там вполне может быть ещё какой-то запас пространства для оптимизаций. Единственное, это уже скорее всего будет специфично конкретно для данной аппаратной конфигурации и решаемой задачи. Но попробовать можно.

Ну вот с таким детальным разъяснением становится понятно, зачем демо-пример именно в таком виде. Если вы применили такой подход на реальных данных, можете поделиться результатом? Хотя бы порядок ускорения на N ядер (как у вас тут, 3.8 раз на 12 логических ядер)?

Степень ускорения в первую очередь зависит от количества строк входных данных и «тяжести» обрабатывающей процедуры, куда в качестве составляющей в некотором смысле входит и количество элементов строки. Замедление происходит либо при слишком малом числе строк (допустим, если будет не 500000 строк, а 5000 – прирост будет уже менее чем двукратный), либо при слишком простой процедуре, требующей мало вычислений. Со случаями, где сильно влияла бы структура данных или порядок значений, я не сталкивался, так что условно можно считать, что значительного влияния нет и случайные значения в этом плане мало отличимы от реальных.
Но в любом случае, есть какой-то предел, выше которого подняться в конкретном случае не получается. При работе с числами ускорение на 12 ядрах у меня в лучшем случае обычно приближается к 4. Судя по тому, что при преобразовании матрицы в таблицу данных ускорение стремится уже не к 4, а к 8 (при существенно общей потере производительности), 4 – это не какой-то абсолютный предел, и на данных других типов ускорение может быть и выше.
И точно, даже в данном маленьком примере матрица отработала быстрее «быстрых» таблиц:

На моём ПК в данном примере самая быстрая реализация оказалась в пакете matrixStats,
хотя есть ошибка на пределе машинной точности (от -9.095e-13 до 9.095e-13).

library(matrixStats)
start_time <- Sys.time()
a2 <- a^2
bRowSums <- rowSums2(a2)
timediff <- difftime(Sys.time(), start_time)
cat("Расчёт занял: ", timediff, units(timediff))
# Расчёт занял:  0.01400399 secs
identical(b, bRowSums)
# FALSE
summary(b - bRowSums)
#       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
# -9.095e-13 -2.842e-14  0.000e+00  8.440e-17  2.842e-14  9.095e-13


В data.table расчёт дольше, хотя при работе с десятками миллионов записей широких таблиц data.table меня много раз выручал.

library(data.table)
dt <- as.data.table(a)
start_time <- Sys.time()
dt2 <- dt^2
colNames <- names(dt)
dt2[, SUM := rowSums(.SD), .SDcols = colNames]
timediff <- difftime(Sys.time(), start_time)
cat("Расчёт занял: ", timediff, units(timediff))
# Расчёт занял:  0.05101109 secs
identical(dt2$SUM, b)
# TRUE
Согласен, что пример очень неудачный, но благодарю за примеры в пакете, о котором не знал. Сам привык выполнять параллельные расчёты в R в doParallel, который на данном примере вообще не успевал ядра нагрузить, но расчёт через %dopar% выполнился аж за 24 секунды, а при монопотоке как в публикации в 4 секунды:
library(tictoc)
tic()
for(i in 1:dim(a)[1]) {
  b[i] <- sum.of.squares(a[i,])
}
toc()
# 4.16 sec elapsed

library(doParallel)
cl <- makeCluster(16)
registerDoParallel(cl)
tic()
parB <- foreach(i=seq(dim(a)[1]), .inorder = FALSE, .combine = rbind) %dopar% {
                    sum.of.squares(a[i,])
                  }
toc()
# 23.62 sec elapsed


При этом в более пригодном формате хранения больших таблиц (в data.table) пример считается в разы быстрее, если нам действительно нужно что-то построчно сделать:


library(data.table)
library(matrixStats)
dt <- as.data.table(a)
dt2 <- dt^2
colNames <- names(dt)
tic()
dt2[, SUM := rowSums(.SD), .SDcols = colNames]
toc()
# 0.03 sec elapsed
identical(dt2$SUM, b)
# TRUE


Реализация в пакете future.apply у меня в разы медленнее оказалась на 3950х:

library("future.apply")
plan(multisession)
tic()
bFuture <- future_apply(a, 1, function(x) sum.of.squares(x))
toc()
# 2.61 sec elapsed
Спасибо, интересно взглянуть на вопрос с другой стороны. Трюк с поэлементной обработкой значений пригоден только конкретно для этого учебного примера, и в большинстве реальных случаев будет неприменим. Попробовал переложить исходные данные в data.table, но само по себе это это дало только небольшое замедление. Возможно выигрыш будет на более крупных массивах, при случае проверю.
Вообще, когда доходит до вычислительных задач, то обилие пакетов в R из преимущества становится скорее проблемой. Почти никто не пользуется чистым R, все ставят пакеты. Пакетов тысячи, их функции пересекаются и дублируют друг друга, как эти функции соотносятся в плане производительности — во многих случаях просто неизвестно, некоторые вещи никто никогда не сравнивал. У каждого из пользователей — своя уникальная экосистема из десятков привычных ему пакетов, то есть фактически вместо единого языка существует огромное число частично пересекающихся между собой диалектов. Можно стать экспертом по одному диалекту, можно владеть несколькими, но практически невозможно охватить их все. Из-за этого теряется какая-либо возможность сделать путеводитель по правильным техникам решения тех или иных задач в R, а изучение возможностей языка начинает носить характер натурных экспериментов и какой-то алхимии. Иногда это даже забавно, но когда дело касается вычислений — весёлого мало. Хочется идти по правильному пути, писать сразу более-менее начисто — а официального правильного пути нет, есть тысячи пересекающихся тропинок и надо самому сравнивать их между собой. Это становится реальной проблемой, слишком уж много времени уходит на сравнения различных способов. Многие советуют переписывать тяжёлые части на каком-то более производительном языке и вызывать их из R, возможно в некотором смысле это самое правильное решение, хотя и достаточно громоздкое. Пока что для меня это открытая проблема.
У вас выбран не самый удачный демонстрационный пример.
Такой вариант в лоб, который дает тот же результат,
rowSums(a*a)
вычисляется на моей машине в среднем за 0.016 секунд.
Среднее время расчета вашего исходного примера с циклом у меня занимает 4.43 секунды.
Я взял максимально упрощённый пример, чтобы не нужно было устраивать лирическое отступление и рассказывать, что в нём вычисляется. При этом хотелось обойтись без псевдокода, а привести какое-то вычисление, которое действительно можно запустить и выполнить. Поскольку это только демонстрационный пример, цели максимально затюнить конкретно это вычисление не ставилось, я даже не придумывал ему защиту от решения другими способами. Допустим, завтра нам понадобится посчитать для каждой строки не сумму квадратов всех значений, а сумму квадратов разниц между текущим значением и средним значением по ряду. Послезавтра — медианное значение дельт между текущим значением и всеми остальными. Тут трюки с поэлементной обработкой уже не подойдут, а предложенный в статье способ довольно универсален, ему не очень важно, какое вычисление работает внутри
Да-да ваша цель ясна, полностью поддерживаю вас. Просто пример вызывал некоторое недоумение. Он как раз демонстрирует сильную сторону базового R.

На мой взгляд, пример из виньетки future.apply тоже прост и показателен — средствами базового R здесь future.apply не превзойти.
x <- replicate(1e5, runif(100), simplify = F)
y <- lapply(x, quantile, probs = 1:3/4)
yy <- future_lapply(x, quantile, probs = 1:3/4)

all.equal(y, yy)
Зарегистрируйтесь на Хабре , чтобы оставить комментарий

Публикации

Истории