x1o / lr-gd-r Goto Github PK

В курсе статистики доказывается, что оптимальная оценка $\hat{\boldsymbol{\beta}}$ может быть получена аналитически: $$\hat{\boldsymbol{\beta}}=(\mathbf{X}^{\intercal}\mathbf{X})^{-1}\mathbf{X}^{\intercal}\mathbf{y}$$

• Сгенерируем данные для n = 1e6, d = 100
• Попробуем посчитать аналитическим методом. Получилось?
• Сравним время исполнения обоих алгоритмов

Аналитический метод предпочтительней для случаев, когда данные помещаются в память. Ясно, что если $\mathbf{X}$ больше объема оперативной памяти, то посчитать аналитически не получится. Вдобавок, считать matrix inverse вычислительно тоже очень дорого (кстати, сколько?).

• Посчитайте сложность градиентного спуска и сравните со сложностью аналитического метода.

Наконец, градиентный спуск можно легко запараллелить, что будем делать в следующих тикетах.

• Будем изменять количество нод в кластере от 1 до какого-то, подходящего вашей системе; для каждого кластера прогоним алгоритм на достаточно больших данных (можно n = 1e7, d = 10 или даже d = 20); время выполнения сохраним (можно посчитать используя system.time()) и организуем в таблицу.

Прирост в скорости может быть заметен при использовании больше двух нод. На макбуке с 12 ядрами автора задания удалось сократить время выполнения примерно в 1.5 раза на двух нодах по сравнению с одной, причем время исполнения продолжало уменьшаться вплоть до 12 нод. На 4-ядерном ThinkPad прирост производительности не был существенным.

Но вы можете использовать всю ту же технику для оптимизации более сложных моделей на неограниченно большом числе вычислительных узлов. Ура!

• Посчитаем, как будет выглядеть элемент градиента $$\frac{\partial R}{\partial\beta_{i}}$$
• Перепишем его в векторном виде, чтобы считался сразу весь градиент $\nabla R$
• Используя $\mathbf{X}$ из #2, сгенерируем $\mathbf{y}$, однако не будем зашумлять. Попробуем восстановить сгенерированную гиперплоскость.
• Напишем функцию gd_step, реализующую шаг градиентного спуска для заданного вектора параметров b.
• Попробуем для случайного вектора $\boldsymbol{\beta}_{\mathrm{init}}$ прогнать шаг спуска. Считает?
• Попробуем запустить в цикле gd_step (шагов 200). $\alpha$ можно положить равным 1e-2. Полученный вектор похож на тот, который использовался для генерации $\mathbf{y}$? Сколько разница?
• Кстати говоря, не стоит ли попробовать запустить оптимизацию на другом $\boldsymbol{\beta}_\mathrm{init}$ - а ну как попали в локальный минимум?

• Напишем функцию risk, считающую среднее квадрата разницы между f на заданном векторе параметров и данным $\boldsymbol{y}$
• Обернем логику из #4 в функцию gd. Теперь, однако, вдобавок к максимальному числу шагов добавим еще один критерий останова -- сначала надо придумать, какой! Есть по крайней мере два варианта, когда оптимизацию можно прекращать. gd должна возвращать текущее значение риска, оценку $\hat{\boldsymbol{\beta}}$ и количество итераций.
• Напишем функцию plot_risk, рисующую значения риска в зависимости от номера итерации.
• Заоптимизируем и нарисуем plot_risk.
• Попробуем то же самое, только в этот раз на зашумленных данных. Что видно?
• При каких значениях $\alpha$ алгоритм расходится? Исследуем сходимость алгоритма в зависимости от значений $\alpha$.

• Сгенерируем данные с $d = 2$ и соответствующие наблюдения
• Трехмерные данные рисовать сложно. Один вариант - нарисовать "контурную проекцию". Для этого придется генерировать сетку из значений $\mathbf{X}$ и считать f на ней. Можно попробовать нарисовать то что получилось при помощи ggplot2::geom_contour с aes(colour = after_stat(level)) и наложить сверху точки из $\mathbf{X}$

• Данная формула описывает модель линейной регрессии в скалярном виде: $$f(\mathbf{x}_i\mid\boldsymbol{\beta})=\beta_{0}+\sum_{j=1}^{d}\beta_{j}x_{ij}.$$ Давайте напишем то же, но используя скалярное произведение, чтобы считался сразу вектор $\mathbf{y}$. Соглашение такое: в матрице данных $\mathbf{X}$ по строкам стоят $\mathbf{x}_i^{\intercal}$; соответственно, номер столбца соответствует номеру признака.
• Реализуем функцию make_X, генерирующую $\mathbf{X}$.
• Реализуем функцию f, реализующую модель линейной регрессии.
• Сгенерируем данные большей размерности: пусть $n = 1000$, $d = 5$, $x_{\min} = -5, x_{\max} = 5$, $\beta_{\min} = -5$, $\beta_{\max} = 5$. Выберем случайные значения для вектора параметров $\boldsymbol{\beta}$, посчитаем наблюдения и зашумим.
• Исследуем распределения полученных данных. Для визуализации хорошо подходит тип графиков "pairs plot". Он позволяет смотреть, как соотносятся признаки друг с другом и с наблюдениями. Рекомендация: GGally::ggpairs.

Как можно было заметить по pairs, $\mathbf{y}$ распределен вовсе не равномерно. А как? Фактически, результирующее распределение есть сумма равномерно распределенных случайных величин, однако их сумма не есть равномерное распределение. Вывести результат сложно даже со знанием теории вероятностей. Для ознакомления можно почитать вот эту статью https://arxiv.org/pdf/math/0411298v1.pdf

См. также #3.

• Как и что будем параллелить? Напишите формально. Опишите идею алгоритма.
• Модифицируем gd, gd_step под параллельное исполнение.
• Попробуем прогнать достаточно большой пример, который использовали на последовательной реализации. Ну как? В зависимости от конфигурации, прироста в скорости может не быть совсем.

• Начнем с простого и сгенерируем данные для "парной" линейной регрессии (с одним признаком). Пусть количество наблюдений n = 10, параметры модели $\beta_0 = 2, \beta_1 = 4$), значения признака выбираем с одинаковой вероятностью из отрезка $[-5, 5]$, дисперсия шума $\sigma^2 = 2$.
• Нарисуем полученные данные.

Чтобы при разных выполнениях кода, генерирующего случайные значения, получать одинаковые результаты, удобно выполнить set.seed(1) в самом начале.

В R можно делать параллельные вычисления различными способами. Нам нужен такой, который бы позволял запускать вычислительные ноды, копировать на них произвольные данные и держать из запущенными между итерациями градиентного спуска. Подойдет пакет parallel, входящий в базовую поставку R, см. ?parallel::`parallel-package` .

Рабочие ноды могут быть где угодно, но для нашего случая будем запускать ноды на рабочей машине в виде R процессов.

• Определим количество доступных ядер: detectCores(). Условно, будем использовать по одной ноде на одно ядро. Одно ядро оставим для остальной системы.
• Создадим кластер: cl <- makePSOCKcluster(n_workers). Тут удобно открыть свой менджер процессов и убедиться, что запустилось нужно количество R сессий.
• На каждую ноду нужно скопировать те объекты, которыми будем пользоваться, при помощи clusterExport. Удобно складывать все объекты во временное окружение, например так:

temp_env <- new.env()
temp_env$x <- 1:10
temp_env$g <- function(x) sum(x)
clusterExport(cl, c('g', 'x'), envir = temp_env)

• Можно также копировать данные на ноды выборочно, к примеру вот так можно отправить переменную x только на первую ноду:

temp_env$x <- 1:5
clusterExport(cl[1], 'x', envir = temp_env)

• Чтобы выполнить функцию g(x) на всех кластерах, используем clusterEvalQ: clusterEvalQ(cl, g(x))

Можно использовать и clusterCall, и другие функции, но это сложно, потому что приходится особенным образом формировать выполняемую функцию (чтобы на каждую ноду не копировалось всё её окружение).

• Потушим кластер: stopCluster(cl)

Другие пакеты для параллельного программирования:

• multidplyr для распараллеливания вычислений на дата-фрейме
• furrr для переноса map-reduce логики в параллельный мир на stateless кластерах