11.3 Градиентный спуск

Хотя градиентный спуск (gradient descent) редко применяют в глубоком обучении напрямую, понимание этого метода — ключ к пониманию алгоритма стохастического градиентного спуска из следующего раздела. Например, из-за слишком большой скорости обучения задача оптимизации может расходиться; это явление уже проявляется в градиентном спуске. Аналогично, предобусловливание (preconditioning) — распространённый приём в градиентном спуске, который переносится и на более совершенные алгоритмы. Начнём с простого одномерного случая.

11.3.1 Одномерный градиентный спуск

Почему алгоритм градиентного спуска может оптимизировать целевую функцию? Одномерный градиентный спуск даёт хорошее представление об этом. Рассмотрим класс непрерывно дифференцируемых вещественных функций $f:\mathbb{R}\to \mathbb{R}$. С помощью разложения в ряд Тейлора получаем

$$f(x+ϵ)=f(x)+ϵf^{\prime }(x)+\mathcal{O}({ϵ}^2).\text{\hspace{1em}(11.3.1)}$$

То есть в первом приближении $f(x+ϵ)$ определяется значением функции $f(x)$ и её первой производной $f^{\prime }(x)$ в точке $x$. Можно предположить, что при движении в направлении антеградиента на малую величину $ϵ$ значение $f$ уменьшится. Для простоты выберем фиксированный шаг $\eta >0$ и положим $ϵ=-\eta f^{\prime }(x)$. Подставив это в разложение Тейлора, получим

$$f(x-\eta f^{\prime }(x))=f(x)-\eta f^{\prime 2}(x)+\mathcal{O}({\eta }^2{f}^{\prime 2}(x)).\text{\hspace{1em}(11.3.2)}$$

Если производная $f^{\prime }(x)\ne 0$ не обращается в нуль, мы можем продолжать разложение, поскольку $\eta f^{\prime 2}(x)>0$. Кроме того, мы всегда можем выбрать $\eta$ достаточно малым, чтобы члены более высокого порядка стали пренебрежимо малы. Поэтому

$$f(x-\eta f^{\prime }(x))⪅f(x).$$

Это означает, что если итеративно обновлять $x$ по правилу

$$x\leftarrow x-\eta f^{\prime }(x),$$

то значение функции $f(x)$, скорее всего, будет уменьшаться. Таким образом, в градиентном спуске мы сначала выбираем начальное значение $x$ и константу $\eta >0$, а затем итеративно обновляем $x$ до выполнения условия остановки — например, когда модуль градиента $|f^{\prime }(x)|$ становится достаточно малым или когда число итераций достигает заданного значения.

Покажем, как реализовать градиентный спуск. Для простоты возьмём целевую функцию $f(x)=x^2$. Хотя мы знаем, что $f(x)$ достигает минимума при $x=0$, мы всё же используем эту простую функцию, чтобы наблюдать за изменением $x$.

%matplotlib inline
from d2l import mxnet as d2l
from mxnet import np, npx
npx.set_np()

%matplotlib inline
from d2l import torch as d2l
import numpy as np
import torch

%matplotlib inline
from d2l import tensorflow as d2l
import numpy as np
import tensorflow as tf

%matplotlib inline
from d2l import paddle as d2l
import warnings
warnings.filterwarnings("ignore")
import numpy as np
import paddle

def f(x):  # Целевая функция
    return x ** 2

def f_grad(x):  # Градиент (производная) целевой функции
    return 2 * x

Далее возьмём $x=10$ в качестве начального значения и предположим, что $\eta =0.2$. Применив градиентный спуск к $x$ на 10 итераций, мы увидим, что значение $x$ в конце концов приблизится к оптимальному решению.

def gd(eta, f_grad):
    x = 10.0
    results = [x]
    for i in range(10):
        x -= eta * f_grad(x)
        results.append(float(x))
    print(f'эпоха 10, x: {x:f}')
    return results

results = gd(0.2, f_grad)

def gd(eta, f_grad):
    x = 10.0
    results = [x]
    for i in range(10):
        x -= eta * f_grad(x)
        results.append(float(x))
    print(f'эпоха 10, x: {float(x):f}')
    return results

results = gd(0.2, f_grad)

Процесс оптимизации $x$ можно изобразить так.

def show_trace(results, f):
    n = max(abs(min(results)), abs(max(results)))
    f_line = d2l.arange(-n, n, 0.01)
    d2l.set_figsize()
    d2l.plot([f_line, results], [[f(x) for x in f_line], [
        f(x) for x in results]], 'x', 'f(x)', fmts=['-', '-o'])

show_trace(results, f)

def show_trace(results, f):
    n = max(abs(min(results)), abs(max(results)))
    f_line = d2l.arange(-n, n, 0.01, dtype='float32')
    d2l.set_figsize()
    d2l.plot([f_line, results], [[f(x) for x in f_line], [
        f(x) for x in results]], 'x', 'f(x)', fmts=['-', '-o'])

show_trace(results, f)

11.3.1.1 Скорость обучения

Скорость обучения (learning rate) определяет, сойдётся ли целевая функция к локальному минимуму и когда это произойдёт. Скорость обучения $\eta$ задаёт разработчик алгоритма. Заметим, что если выбрать слишком маленькую скорость обучения, обновления $x$ будут очень медленными, и понадобится больше итераций. Например, рассмотрим ход той же задачи оптимизации при $\eta =0.05$. Как видно ниже, даже за 10 шагов мы всё ещё далеки от оптимального решения.

show_trace(gd(0.05, f_grad), f)

Напротив, если выбрать слишком большую скорость обучения, $\left|\eta f^{\prime }(x)\right|$ может оказаться слишком большим для линейного приближения первого порядка. То есть член $\mathcal{O}({\eta }^2{f}^{\prime 2}(x))$ в формуле (11.3.2) может стать существенным. В этом случае итерации по $x$ уже не гарантируют уменьшения $f(x)$. Например, при $\eta =1.1$ значение $x$ перескакивает оптимальное решение $x=0$ и постепенно расходится.

show_trace(gd(1.1, f_grad), f)

11.3.1.2 Локальные минимумы

Чтобы продемонстрировать градиентный спуск для невыпуклой функции, рассмотрим функцию $f(x)=x\cdot \cos (cx)$, где $c$ — некоторая константа. Эта функция имеет бесконечно много локальных минимумов. В зависимости от выбранной скорости обучения мы в итоге можем получить лишь одно из многих решений. Следующий пример показывает, как (нереалистично) большая скорость обучения приводит к худшему локальному минимуму.

c = d2l.tensor(0.15 * np.pi)

def f(x):  # Целевая функция
    return x * d2l.cos(c * x)

def f_grad(x):  # Градиент целевой функции
    return d2l.cos(c * x) - c * x * d2l.sin(c * x)

show_trace(gd(2, f_grad), f)

11.3.2 Многомерный градиентный спуск

Теперь, когда одномерный случай понятен, рассмотрим $\mathbf{x}=[x_1,x_2,\dots ,x_d{]}^{\top }$. То есть целевая функция $f:{\mathbb{R}}^d\to \mathbb{R}$ отображает вектор в скаляр. Соответственно, её градиент тоже многомерный: это вектор из $d$ частных производных:

$$\nabla f(\mathbf{x})=[\frac{\partial f(\mathbf{x})}{\partial x_1},\frac{\partial f(\mathbf{x})}{\partial x_2},\dots ,\frac{\partial f(\mathbf{x})}{\partial x_d}{]}^{\top }.$$

Каждый элемент градиента $\partial f(\mathbf{x})/\partial x_i$ показывает скорость изменения $f$ в точке $\mathbf{x}$ по входу $x_i$. Как и в одномерном случае, для многомерных функций можно использовать соответствующее разложение Тейлора. А именно,

$$f(\mathbf{x}+\mathbf{ϵ})=f(\mathbf{x})+{\mathbf{ϵ}}^{\top }\nabla f(\mathbf{x})+\mathcal{O}(\Vert \mathbf{ϵ}{\Vert }^2).\text{\hspace{1em}(11.3.3)}$$

Иными словами, с точностью до членов второго порядка по $\mathbf{ϵ}$ направление наибыстрейшего убывания задаёт антеградиент $-\nabla f(\mathbf{x})$. Выбрав подходящую скорость обучения $\eta >0$, получаем типичный алгоритм градиентного спуска:

$$\mathbf{x}\leftarrow \mathbf{x}-\eta \nabla f(\mathbf{x}).$$

Как этот алгоритм работает на практике? Рассмотрим целевую функцию $f(\mathbf{x})=x_1^2+2x_2^2$ с двумерным вектором $\mathbf{x}=[x_1,x_2{]}^{\top }$ на входе и скаляром на выходе. Градиент задаётся формулой $\nabla f(\mathbf{x})=[2x_1,4x_2{]}^{\top }$. Мы будем наблюдать траекторию $\mathbf{x}$ при градиентном спуске из начальной точки $[-5,-2]$.

Нам понадобятся две вспомогательные функции: первая — это функция update, которую мы применим к начальному значению 20 раз; вторая — функция, которая отображает траекторию $\mathbf{x}$.

def train_2d(trainer, steps=20, f_grad=None):  #@save
    """Оптимизация 2D-целевой функции с помощью пользовательского тренера"""
    # s1 и s2 — внутренние переменные состояния, которые будут использованы позже
    x1, x2, s1, s2 = -5, -2, 0, 0
    results = [(x1, x2)]
    for i in range(steps):
        if f_grad:
            x1, x2, s1, s2 = trainer(x1, x2, s1, s2, f_grad)
        else:
            x1, x2, s1, s2 = trainer(x1, x2, s1, s2)
        results.append((x1, x2))
    print(f'эпоха {i + 1}, x1: {float(x1):f}, x2: {float(x2):f}')
    return results

def show_trace_2d(f, results):  #@save
    """Отображение траектории 2D-переменных в ходе оптимизации"""
    d2l.set_figsize()
    d2l.plt.plot(*zip(*results), '-o', color='#ff7f0e')
    x1, x2 = d2l.meshgrid(d2l.arange(-5.5, 1.0, 0.1),
                          d2l.arange(-3.0, 1.0, 0.1))
    d2l.plt.contour(x1, x2, f(x1, x2), colors='#1f77b4')
    d2l.plt.xlabel('x1')
    d2l.plt.ylabel('x2')

def show_trace_2d(f, results):  #@save
    """Отображение траектории 2D-переменных в ходе оптимизации"""
    d2l.set_figsize()
    d2l.plt.plot(*zip(*results), '-o', color='#ff7f0e')
    x1, x2 = d2l.meshgrid(d2l.arange(-5.5, 1.0, 0.1),
                          d2l.arange(-3.0, 1.0, 0.1), indexing='ij')
    d2l.plt.contour(x1, x2, f(x1, x2), colors='#1f77b4')
    d2l.plt.xlabel('x1')
    d2l.plt.ylabel('x2')

def train_2d(trainer, steps=20, f_grad=None):  #@save
    """Оптимизация 2D-целевой функции с помощью пользовательского тренера"""
    # s1 и s2 — внутренние переменные состояния, которые будут использованы позже
    x1, x2, s1, s2 = -5, -2, 0, 0
    results = [(x1, x2)]
    for i in range(steps):
        if f_grad:
            x1, x2, s1, s2 = trainer(x1, x2, s1, s2, f_grad)
        else:
            x1, x2, s1, s2 = trainer(x1, x2, s1, s2)
        results.append((x1, x2))
    print(f'эпоха {i + 1}, x1: {float(x1):f}, x2: {float(x2):f}')
    return results

def show_trace_2d(f, results):  #@save
    """Отображение траектории 2D-переменных в ходе оптимизации"""
    d2l.set_figsize()
    d2l.plt.plot(*zip(*results), '-o', color='#ff7f0e')
    x1, x2 = d2l.meshgrid(d2l.arange(-5.5, 1.0, 0.1, dtype='float32'),
                          d2l.arange(-3.0, 1.0, 0.1, dtype='float32'))
    d2l.plt.contour(x1, x2, f(x1, x2), colors='#1f77b4')
    d2l.plt.xlabel('x1')
    d2l.plt.ylabel('x2')

Далее мы наблюдаем траекторию переменной $\mathbf{x}$ при скорости обучения $\eta =0.1$. Видно, что после 20 шагов значение $\mathbf{x}$ приближается к минимуму в точке $[0,0]$. Прогресс идёт довольно уверенно, но медленно.

def f_2d(x1, x2):  # Целевая функция
    return x1 ** 2 + 2 * x2 ** 2

def f_2d_grad(x1, x2):  # Градиент целевой функции
    return (2 * x1, 4 * x2)

def gd_2d(x1, x2, s1, s2, f_grad):
    g1, g2 = f_grad(x1, x2)
    return (x1 - eta * g1, x2 - eta * g2, 0, 0)

eta = 0.1
show_trace_2d(f_2d, train_2d(gd_2d, f_grad=f_2d_grad))

11.3.3 Адаптивные методы

Как мы видели в разделе 11.3.1.1, выбрать «правильную» скорость обучения $\eta$ непросто. Если её выбрать слишком маленькой, прогресса почти не будет; если слишком большой — решение начнёт колебаться и может даже расходиться. А что, если бы мы могли определять $\eta$ автоматически или вообще не выбирать скорость обучения? Выход может дать метод второго порядка, который учитывает не только значение и градиент целевой функции, но и её кривизну. Хотя из-за вычислительных затрат такие методы напрямую неприменимы в глубоком обучении, они дают полезную интуицию для построения продвинутых алгоритмов оптимизации, которые могут воспроизводить многие желаемые свойства описанных ниже алгоритмов.

11.3.3.1 Метод Ньютона

Вспомним разложение Тейлора для функции $f:{\mathbb{R}}^d\to \mathbb{R}$. Фактически его можно записать так:

$$f(\mathbf{x}+\mathbf{ϵ})=f(\mathbf{x})+{\mathbf{ϵ}}^{\top }\nabla f(\mathbf{x})+\frac{1}{2}{\mathbf{ϵ}}^{\top }{\nabla }^{2}f(\mathbf{x})\mathbf{ϵ}+\mathcal{O}(\Vert \mathbf{ϵ}{\Vert }^3).\text{\hspace{1em}(11.3.4)}$$

Чтобы не громоздить обозначения, определим $\mathbf{H}\stackrel{\mathrm{def}}{=}{\nabla }^{2}f(\mathbf{x})$ как гессиан $f$; это матрица размера $d\times d$. Для малых $d$ и простых задач $\mathbf{H}$ вычислить легко. Однако для глубоких нейронных сетей $\mathbf{H}$ может быть очень большой: хранение $\mathcal{O}(d^2)$ элементов обходится дорого, а вычисление через обратное распространение может сделать ситуацию ещё хуже. Тем не менее пока отвлечёмся от этих соображений и посмотрим, какой алгоритм из этого получится.

В конце концов, минимум $f$ удовлетворяет условию $\nabla f=0$. Следуя правилам математического анализа из раздела 2.4, возьмём производную от $\mathbf{ϵ}$ по выражению (11.3.4) и отбросим несущественные члены высшего порядка. В результате получаем

$$\nabla f(\mathbf{x})+\mathbf{H}\mathbf{ϵ}=0\text{ и, следовательно, }\mathbf{ϵ}=-{\mathbf{H}}^{-1}\nabla f(\mathbf{x}).$$

То есть в рамках задачи оптимизации нам нужно обратить гессиан $\mathbf{H}$.

Для простого примера $f(x)=\frac{1}{2}{x}^2$ имеем $\nabla f(x)=x$ и $\mathbf{H}=1$. Поэтому для любого $x$ получаем $ϵ=-x$. Иными словами, для идеальной сходимости достаточно одного шага без каких-либо настроек. Здесь нам повезло: разложение Тейлора точное, так как $f(x+ϵ)=\frac{1}{2}{x}^2+ϵx+\frac{1}{2}{ϵ}^2$.

Рассмотрим другие задачи. Пусть задана выпуклая функция $c(x)=\cosh (cx)$ (гиперболический косинус), где $c$ — некоторая константа. Видно, что после нескольких итераций мы получаем глобальный минимум в точке $x=0$.

c = d2l.tensor(0.5)

def f(x):  # Целевая функция
    return d2l.cosh(c * x)

def f_grad(x):  # Градиент целевой функции
    return c * d2l.sinh(c * x)

def f_hess(x):  # Гессиан целевой функции
    return c**2 * d2l.cosh(c * x)

def newton(eta=1):
    x = 10.0
    results = [x]
    for i in range(10):
        x -= eta * f_grad(x) / f_hess(x)
        results.append(float(x))
    print('эпоха 10, x:', x)
    return results

show_trace(newton(), f)

Теперь рассмотрим невыпуклую функцию, например $f(x)=x\cos (cx)$, где $c$ — некоторая константа. Обратите внимание: в методе Ньютона в конце концов мы делим на гессиан. Это значит, что если вторая производная отрицательна, значение $f$ может начать расти. Это фатальный недостаток алгоритма! Посмотрим, что происходит на практике.

c = d2l.tensor(0.15 * np.pi)

def f(x):  # Целевая функция
    return x * d2l.cos(c * x)

def f_grad(x):  # Градиент целевой функции
    return d2l.cos(c * x) - c * x * d2l.sin(c * x)

def f_hess(x):  # Гессиан целевой функции
    return - 2 * c * d2l.sin(c * x) - x * c**2 * d2l.cos(c * x)

show_trace(newton(), f)

Результат получился катастрофическим. Как это исправить? Один из способов — заменить гессиан его абсолютным значением; другая стратегия — снова ввести скорость обучения. Это может показаться противоречащим исходной идее, но это не совсем так: знание второго порядка позволяет быть осторожным там, где кривизна велика, и использовать большую скорость обучения там, где целевая функция более пологая. Посмотрим, как это работает при несколько меньшей скорости обучения, например $\eta =0.5$. Как видим, получился довольно эффективный алгоритм.

show_trace(newton(0.5), f)

11.3.3.2 Анализ сходимости

Здесь мы на примере некоторых выпуклых функций $f$ проанализируем скорость сходимости метода Ньютона. Эти функции трижды дифференцируемы, а вторая производная не равна нулю, то есть $f^{\prime \prime }>0$. Поскольку доказательство для многомерного случая является прямым обобщением доказательства для приведённого ниже одномерного случая и не добавляет понимания задачи, мы опускаем многомерный случай.

Обозначим через $x^{(k)}$ значение $x$ на $k$-й итерации, а через $e^{(k)}\stackrel{\mathrm{def}}{=}{x}^{(k)}-x^{\ast }$ — расстояние до оптимальности на $k$-й итерации. Разложив в ряд Тейлора, получаем, что условие $f^{\prime }(x^{\ast })=0$ можно записать в виде

$$0=f^{\prime }(x^{(k)}-e^{(k)})=f^{\prime }(x^{(k)})-e^{(k)}{f}^{\prime \prime }(x^{(k)})+\frac{1}{2}(e^{(k)}{)}^2{f}^{\prime \prime \prime }({\xi }^{(k)}),$$

где ${\xi }^{(k)}\in [x^{(k)}-e^{(k)},x^{(k)}]$. Разделив это разложение на $f^{\prime \prime }(x^{(k)})$, получаем

$$e^{(k)}-\frac{f^{\prime }(x^{(k)})}{f^{\prime \prime }(x^{(k)})}=\frac{1}{2}(e^{(k)}{)}^2\frac{f^{\prime \prime \prime }({\xi }^{(k)})}{f^{\prime \prime }(x^{(k)})}.$$

Вспомним предыдущее уравнение обновления $x^{(k+1)}=x^{(k)}-f^{\prime }(x^{(k)})/f^{\prime \prime }(x^{(k)})$. Подставив его и взяв модули обеих частей, получаем

$$\left|e^{(k+1)}\right|=\frac{1}{2}(e^{(k)}{)}^2\frac{\left|f^{\prime \prime \prime }({\xi }^{(k)})\right|}{f^{\prime \prime }(x^{(k)})}.$$

Следовательно, всякий раз, когда мы находимся в ограниченной области, где $\left|f^{\prime \prime \prime }({\xi }^{(k)})\right|/(2f^{\prime \prime }(x^{(k)}))\le c$, ошибка убывает квадратично:

$$\left|e^{(k+1)}\right|\le c(e^{(k)}{)}^2.$$

С другой стороны, специалисты по оптимизации называют сходимость «линейной» при условии $\left|e^{(k+1)}\right|\le \alpha \left|e^{(k)}\right|$, а при этом говорят о «постоянной» скорости сходимости. Заметим, что скорость сходимости в целом оценить нельзя, но как только мы оказываемся вблизи минимума, сходимость становится очень быстрой. Кроме того, такой анализ требует, чтобы $f$ хорошо вела себя при высоких порядках производных, то есть чтобы $f$ не обладала какими-либо «экстремальными» свойствами в том, как быстро меняется её значение.

11.3.3.3 Предобусловливание

Вычислять и хранить полный гессиан очень дорого, и один из способов смягчить эту проблему — «предобусловливание». Оно позволяет обойтись без вычисления всего гессиана и учитывает только его «диагональные» элементы; алгоритм обновляется по формуле:

$$\mathbf{x}\leftarrow \mathbf{x}-\eta \mathrm{diag}(\mathbf{H}{)}^{-1}\nabla f(\mathbf{x}).$$

Это не так точно, как полный метод Ньютона, но всё же гораздо лучше, чем обходиться без него. Почему предобусловливание работает? Представьте, что одна переменная измеряет высоту в миллиметрах, а другая — в километрах. Если оба естественных масштаба выражены в метрах, наша параметризация сильно рассогласована. К счастью, предобусловливание устраняет это несоответствие. Эффективное предобусловливание градиентного спуска равносильно выбору разной скорости обучения для каждой координаты вектора $\mathbf{x}$. Мы увидим в одном из следующих разделов, что предобусловливание стимулировало ряд новшеств в алгоритмах стохастического градиентного спуска.

11.3.3.4 Градиентный спуск и линейный поиск

Одна из ключевых проблем градиентного спуска в том, что можно перескочить цель или продвигаться слишком медленно. Простое решение этой проблемы — комбинировать линейный поиск с градиентным спуском. А именно, мы используем направление, заданное $\nabla f(\mathbf{x})$, и затем проводим двоичный поиск, чтобы определить такую скорость обучения $\eta$, при которой $f(\mathbf{x}-\eta \nabla f(\mathbf{x}))$ принимает минимальное значение.

С точки зрения анализа и доказательств этот алгоритм сходится быстро (см. Boyd and Vandenberghe, 2004). Однако для глубокого обучения он малопригоден, поскольку каждый шаг линейного поиска требует вычисления целевой функции на всём наборе данных, и такой способ реализации слишком дорог.

Резюме

• Величина скорости обучения важна: слишком большая скорость приводит к расходимости модели, слишком малая — к отсутствию прогресса.
• Градиентный спуск может застрять в локальном минимуме и не найти глобальный минимум.
• В моделях высокой размерности настроить скорость обучения сложно.
• Предобусловливание помогает выравнивать масштабы.
• В выпуклых задачах метод Ньютона, как только начинает работать корректно, сходится гораздо быстрее.
• Для невыпуклых задач не следует применять метод Ньютона без дополнительных настроек.

Упражнения

1. Проведите эксперименты с градиентным спуском при разных скоростях обучения и разных целевых функциях.
2. Реализуйте линейный поиск на отрезке $[a,b]$ для минимизации выпуклой функции.
   1. Нужны ли производные для двоичного поиска, то есть для выбора между $[a,(a+b)/2]$ и $[(a+b)/2,b]$?
   2. Насколько быстро сходится алгоритм?
   3. Реализуйте этот алгоритм и примените его для нахождения минимума $\log (\exp (x)+\exp (-2x-3))$.
3. Постройте целевую функцию, определённую на ${\mathbb{R}}^2$, для которой градиентный спуск будет очень медленным. Подсказка: масштабы разных координат должны различаться.
4. Реализуйте облегчённую версию метода Ньютона с предобусловливанием.
   1. Используйте диагональный гессиан в качестве предобусловливателя.
   2. Используйте его абсолютное значение вместо фактического (возможно, со знаком).
   3. Примените это к задаче выше.
5. Примените описанный алгоритм к нескольким целевым функциям (выпуклым и невыпуклым). Что будет, если повернуть координаты на $45$ градусов?

MXNet: Обсуждение

PyTorch: Обсуждение

TensorFlow: Обсуждение

PaddlePaddle: Обсуждение

Литература

• Boyd, S., & Vandenberghe, L. (2004). Convex Optimization. Cambridge University Press.