Реализация алгоритма Левенберга-Марквардта для оптимизации нейронных сетей на TensorFlow10.12.2018 08:06

Это tutorial по библиотеке TensorFlow. Рассмотрим её немного глубже, чем в статьях про распознавание рукописных цифр. Это tutorial по методам оптимизации. Совсем без математики здесь не обойтись. Ничего страшного, если вы её совершенно забыли. Вспомним. Не будет никаких формальных доказательств и сложных выводов, только необходимый минимум для интуитивного понимания. Для начала небольшая предыстория о том, чем этот алгоритм может быть полезен при оптимизации нейронной сети.


kavvkhi2wic8-bj1_-6i4zjvkpe.jpeg

Полгода назад друг попросил показать, как на Python сделать нейросеть. Его компания выпускает приборы для геофизических измерений. Несколько различных зондов в процессе бурения измеряют набор сигналов, связаных с параметрами окружающей скважину среды. В некоторых сложных случаях точно вычислить параметры среды по сигналам долго даже на мощном компьютере, а необходимо интерпретировать результаты измерений в полевых условиях. Возникла идея посчитать на кластере несколько сот тысяч случаев, и на них натренировать нейронную сеть. Так как нейросеть работает очень быстро, её можно использовать для определения параметров, согласующихся с измеренными сигналами, прямо в процессе бурения. Детали есть в статье:

Kushnir, D., Velker, N., Bondarenko, A., Dyatlov, G., & Dashevsky, Y. (2018, October 29). Real-Time Simulation of Deep Azimuthal Resistivity Tool in 2D Fault Model Using Neural Networks (Russian). Society of Petroleum Engineers. doi:10.2118/192573-RU

Одним вечером я показал, как keras реализовать простую нейронную сеть, и друг на работе запустил обучение на насчитанных данных. Через пару дней обсудили результат. С моей точки зрения он выглядел перспективно, но друг сказал, что нужны вычисления с точностью прибора. И если средняя квадратичная ошибка (mean squared error) получилась в районе 1, то нужна была 1е-3. На 3 порядка меньше. В тысячу раз.

Эксперименты с архитектурой нейронной сети, нормализацией данных, подходами к оптимизации не дали почти ничего. Через пару недель друг позвонил и сказал, что он установил MatLab и решил проблему методом Левенберга-Марквардта (далее будем называть LM). Оптимизировалось долго (несколько дней), не работало на GPU, но результат получился нужный. Это звучало как вызов.

Беглый поиск готового оптимизатора LM для keras или TensorFlow не дал результата. Наткнулся только на библиотеку pyrenn, но её функционал показался мне бедным. Решил реализовать самостоятельно. На первый взгляд выглядело всё просто, и двух вечеров должно было хватить. Провозился дольше. Проблемы оказалось две:


  1. TensorFlow. Куча статей, но практически все уровня «а давайте напишем hello world распознавалку рукописных цифр».
  2. Математика. Многое забыл, а авторы математических статей совершенно не заботятся о таких как я: сплошные формулы без объяснений, «очевидно!» и тд.

В итоге написал статью для забывших математику и желающих разобраться в TensorFlow чуть глубже, но без хардкора. В статье много текста и мало кода. Обратный вариант, когда мало текста и много кода, есть здесь Jupyter Notebook Levenberg-Marquardt.


Познакомимся с функцией Розенброка

Тренировочные данные будем генерировать функцией Розенброка, которую часто используют в качестве эталонного теста алгоритмов оптимизации:

$f(x, y) = (a-x)^2+b(y-x^2)^2$


crby_vv7ldxutqj5dxfaehtfnlk.jpeg

Чем же она хороша?


  • Красивый график. Называют долиной Розенброка и непереводимо Rosenbrock’s banana function.
  • Глобальный минимум находится внутри длинной, узкой, параболической плоской долины. Найти долину тривиально, а глобальный минимум сложно.
  • Есть многомерный вариант. Придумать хорошую функцию многих переменных не так уж и просто.

С неё и начнём писать код подключив необходимые для дальнейшей работы библиотеки: 

import numpy as np
import tensorflow as tf
import math

def rosenbrock(x, y, a, b):
    return (a - x)**2 + b*(y - x**2)**2


Сформулируем задачу

Раз уж мы говорили об измерительном приборе, давайте дальше использовать аналогию. Наш прибор в вымышленном мире умеет измерять координаты $(x, y)$ и высоту $z$. Физики изучили мир и заявили:»Да это же Розенброк! Зная координаты, можно точно вычислить высоту, вам не надо её измерять.». Иными словами, учёные дали нам модель $z = rosenbrock(x, y, a, b)$, которая зависит от параметров $(a, b)$. Эти параметры хоть и постоянны в вымышленном мире, но неизвестны. Их нужно найти.

Мы провели ряд экспериментов, которые дали $m$ точек $(x_1, y_1, z_1), (x_2, y_2, z_2),..., (x_m, y_m, z_m)$: 

# (2.5, 2.5) - это реальные параметры, давайте забудем, что они нам известны
data_points = np.array([[x, y, rosenbrock(x, y, 2.5, 2.5)]
                       for x in np.arange(-2, 2.1, 2) 
                       for y in np.arange(-2, 2.1, 2)])
m = data_points.shape[0]

Первый способ оптимизации — попытаться угадать параметры. Воспользуемся библиотекой Numpy: 

x, y = data_points[:, 0], data_points[:, 1]
z = data_points[:, 2]
# а вдруг а=5 и b=5?
a_guess, b_guess = 5., 5.
# Суффикс -hat используется для прогнозов, так проще ориентироваться,
# где реальные данные, а где то, что насчитала модель. В формулах
# соответствующие переменные будут иметь ^ треугольничек сверху -
# шляпку. Отсюда и суффикс hat.
z_hat = rosenbrock(x, y, a_guess, b_guess)

Как понять, насколько мы ошиблись? Посчитать невязки (residuals) — размеры ошибок. $m$ точек дают $m$ невязок — нужен интегральный показатель. Возведём каждую невязку в квадрат и посчитаем среднее:

$MSE(a, b) = \frac{1}{m}\sum_{i=1}^{m}(z_{i} - \widehat{z_{i}})^2$

Такая мера близости называется средняя квадратичная ошибка (mean squared error, далее будем называть mse): 

# r - residuals (невязки)
r = z - z_hat
# mse
loss = np.mean(r**2)
print(loss)

[Out]: 3868.2291666666665

Минимизируя mse, мы решаем задачу о наименьших квадратах (nonlinear squares minimization):


ltqksu80eo7y4h91xedhj8w0djy.png

Видно, что параметры не угадали совершенно.


Сформулируем задачу на TensorFlow

Модель имеет вид $z = rosenbrock(x, y, a, b)$. Приведём её к виду $y = f(x, p)$ (обычно математики пишут $\beta$ вместо $p$, но программисты не используют бету). Теперь модель имеет вид $y = rosenbrock(x, p)$, где $y$ — высота, $x$ — вектор координат из двух элементов (компонент), а $p$ — вектор параметров.

Программисты часто думают о векторах, как об одномерных массивах. Это не совсем корректно. Массив чисел — это средство представления вектора. Представить вектор можно массивом размерностью $N$, двухмерным массивом $1\times N$, и даже массивом $N \times 1$ в тех случаях, когда важен факт, что вектор — это вектор-столбец (например, для умножения матрицы на него):

$ \begin{bmatrix} x_1\\\vdots \\ x_N \end{bmatrix} $

В TensorFlow используется понятие тензор. Тензор, как и массив, может быть одномерным (для представления вектора), двухмерным (для матрицы или вектора-столбца) и любой большей размерности.

# тензоры данных экспериментов ('placeholder' означает, что данные надо 
# передавать каждый раз при запуске вычислений)
x = tf.placeholder(tf.float64, shape=[m, 2])
y = tf.placeholder(tf.float64, shape=[m])
# тензор параметров ('variable' означает, что он хранит состояние)
# инициализируем его нашей догадкой (5, 5)
p = tf.Variable([5., 5.], dtype=tf.float64)
# модель
y_hat = rosenbrock(x[:, 0], x[:, 1], p[0], p[1])
# невязки
r = y - y_hat
# mse (mean squared error)
loss = tf.reduce_mean(r**2)

Код на TensorFlow по форме не отличается от кода на Numpy. По содержанию отличия огромные. Код на Numpy производит вычисление значения mse. Код на TensorFlow не производит никаких вычислений вообще, он формирует граф потока данных (dataflow graph), который может вычислить mse. Очень выносящий мозг момент — работа функции rosenbrock. Мы её используем в обоих случаях. Но когда передаём массивы Numpy, она производит вычисления по формуле и возвращает числа. А когда передаём тензоры TensorFlow, формирует подграф потока данных и возвращает его ребро (edge) в виде тензора. Чудеса полиморфизма, но не стоит ими злоупотреблять:


cw6unzjtqsl_uvl3auasc-x46zq.jpeg

Благодаря наличию такого графа потока данных, TensorFlow в частности умеет вычислять производные автоматически (техникой reverse mode automatic differentiation).

Минутка математики. Блоки «для тех, кто забыл» буду прятать в спойлер.


Производная (число вошло — число вышло)

Мы построили граф потока данных, давайте запустим вычисление mse: 

# при запуске вычислений нам надо передать данные 
# для подстановки в тензоры типа placeholder (координаты и высота)
feed_dict = {x: data_points[:,0:2], y: data_points[:,2]}
# в сессии запускаются операции и вычисления TensorFlow
session = tf.Session()
# инициализируем глобальные переменные графа
session.run(tf.global_variables_initializer())
# запускаем вычисление (оценку) тензора loss (mse)
current_loss = session.run(loss, feed_dict)
print(current_loss)

[Out]: 3868.2291666666665

Получился тот же результат, что и с Numpy. Значит не ошиблись.


Начнём оптимизировать

К сожалению угадать параметры не получилось. Но зато мы:


  1. Задали критерий оптимальности — минимальное значение mse.
  2. Определили варьирующие параметры: вектор $p$ с компонентами $a$, $b$ функции Розенброка.
  3. Пока не подумали про ограничения, но их нет.

На прошлом шаге мы сконструировали граф потока данных с конечным тензором loss (функция потерь). Цель оптимизации — найти значение вектора параметров $p$, при котором значение функции потерь минимальное. Нам повезло, график этой функции очень простой (вогнутый и без локальных минимумов):


vfdohrpseaocd_vpe6e7sp_6f0u.png

Приступаем к оптимизации. Для начала напишем обобщённый цикл: 

# параметры: целевое значение mse, максимальное число шагов, тензор
# для вычисления mse, операция шага оптимизации и данные для тензоров placeholder
def train(target_loss, max_steps, loss_tensor, train_step_op, inputs):
    step = 0
    current_loss = session.run(loss_tensor, inputs)
    # оптимизируем пока не закончились шаги или не добились нужного результата
    while current_loss > target_loss and step < max_steps:
        step += 1
        # логируем прогресс на 1, 2, 4, 8, 16... шагах
        if math.log(step, 2).is_integer():
            print(f'step: {step}, current loss: {current_loss}')
        # делаем оптимизационный шаг
        session.run(train_step_op, inputs)
        current_loss = session.run(loss_tensor, inputs)
    print(f'ENDED ON STEP: {step}, FINAL LOSS: {current_loss}')


Оптимизируем методом наискорейшего градиентного спуска (SGD)

Действия этого метода можно сравнить с ездой дерзкого горнолыжника, который всегда ставит лыжи по линии падения склона (в самом крутом направлении вниз). При этом учитывает только уклон в точке нахождения. И если уклон сильный, то до следующего перестроения лыжник пролетает большое расстояние. При слабом же уклоне движется маленькими шагами. Может как улететь в дерево (алгоритм расходится), так и застрять в яме (локальном минимуме).


oeetijkwervilr42jsjl9wsgsey.jpeg

Записать можно так (меняем $\boldsymbol{p}$ на $\boldsymbol{p} - ...$):

$\boldsymbol{p} \rightarrow \boldsymbol{p}-\alpha [\nabla_{p} loss(\boldsymbol{p})]$

Жирное $\boldsymbol{p}$ подчёркивает, что это точка фактического нахождения — значение вектора параметров на текущем шаге. На первом шаге это наша догадка (5, 5). В формуле есть два интересных момента: $\alpha$ — скорость обучения (learning rate), $\nabla_{p} loss$ — градиент (gradient) функции потерь по вектору параметров.


qx1hem1cns_udmg0gucq1h9vs9k.jpeg

Градиент (вектор вошёл — число вышло)

# вычисляем градиент функции потерь по вектору параметров
grad = tf.gradients(loss, p)[0]
# скорость обучения
learning_rate = 0.0005
# оптимизатор нам нужен, чтобы воспользоваться его методом apply_gradients -
# обновление вектора параметров на градиент со знаком минус
opt = tf.train.GradientDescentOptimizer(learning_rate=1)
# сдвигаем вектор параметров на значение градиента с учётом скорости обучения
sgd = opt.apply_gradients([(learning_rate*grad, p)])
# запускаем цикл оптимизации, сделаем не больше 40000 шагов
session.run(tf.global_variables_initializer())
train(1e-10, 40000, loss, sgd, feed_dict)
print('PARAMETERS:', session.run(p))

[Out]: step: 1, current loss: 3868.2291666666665
       step: 2, current loss: 1381.5379689135807
       [...]
       ENDED ON STEP: 582, FINAL LOSS: 9.698531012270816e-11
       PARAMETERS: [2.50000205 2.49999959]

Потребовалось 582 шага:


b2j3jbxzfk6_ll66vpr83vklwjs.jpeg

Движение в направлении антиградиента


Оптимизируем методом Adam

Не будем дальше углубляться в градиентные методы, но вариаций есть масса. Почитать про них можно в статье Методы оптимизации нейронных сетей. В TensorFlow многие оптимизаторы уже реализованы. Например, Adam: 

# не будем сами вычислять и применять градиенты, 
# а сразу сконструируем шаг оптимизации
adm = tf.train.AdamOptimizer(15).minimize(loss)
# запускаем цикл оптимизации, сделаем не больше 40000 шагов
session.run(tf.global_variables_initializer())
train(1e-10, 40000, loss, adm, feed_dict)
print('PARAMETERS:', session.run(p))

[Out]: step: 1, current loss: 3868.2291666666665
       step: 2, current loss: 34205.72916492336
       [...]
       ENDED ON STEP: 317, FINAL LOSS: 2.424142714263483e-12
       PARAMETERS: [2.49999969 2.50000008]

Справились за 317 шагов. Гораздо быстрее.


Оптимизируем методом Ньютона

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


p_iukfzoa0idjm_cbycfmlye4mw.jpeg

Фактически, и методы градиентного спуска, и методы второго порядка пытаются угадать (аппроксимировать) функцию по текущей точке. Градиентные методы ориентируются только на уклон (slope) графика функции в точке — первую производную. Методы же второго порядка кроме уклона учитывают и кривизну (curvature) — вторую производную: «если кривизна сохранится, то где будет минимум?». Вычисляем и направляемся туда:


xbsnnb1o5cui4x10wspkrgjinaw.jpeg

Построить такую аппроксимацию и вычислить предполагаемую точку минимума можно с помощью ряда Тейлора (Taylor series). Для одномерного случая аппроксимация полиномом второго порядка в точке $a$ выглядит так:

$f(x) \approx f(a) + \frac{f'(a)(x-a)}{1!} + \frac{f''(a)(x-a)^2}{2!}$

Минимум достигается при $x = a - \frac{f'(a)}{f''(a)}$. Многомерный случай выглядит более серьёзно:


dvpd31r1vt2cbts_a5oamiiwgoc.jpeg

Матрица Гессе (вектор вошёл — число вышло)

Матрица Гессе (Hessian matrix) — это квадратная матрица, составленная из вторых производных:

$\boldsymbol{H}y_{x} = \begin{bmatrix} \frac{\partial^2y}{\partial x_1^2} & \frac{\partial^2y}{\partial x_1 \partial x_2} & \cdots & \frac{\partial^2y}{\partial x_1 \partial x_N} \\ \frac{\partial^2y}{\partial x_2 \partial x_1} & \frac{\partial^2y}{\partial x_2^2} & \cdots & \frac{\partial^2y}{\partial x_2 \partial x_N} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial^2y}{\partial x_N \partial x_1} & \frac{\partial^2y}{\partial x_N \partial x_2} & \cdots & \frac{\partial^2y}{\partial x_N^2} \end{bmatrix}$

Аппроксимация полиномом второго порядка для функции от вектора через градиент и матрицу Гессе в точке $a$ выглядит так:

$f(x) \approx f(a) + (x-a)^\intercal [\nabla_{x}f(a)] + \frac{1}{2!}(x-a)^\intercal [\boldsymbol{H}f_{x}(a)](x-a)$

Минимум достигается при $x =a - [\boldsymbol{H}f_{x}(a)]^{-1}[\nabla_{x}f(a)]$. Форма практически совпадает с одномерным случаем: заменили первую производную на градиент, вторую на матрицу Гессе и сделали поправку на работу с векторами. Делить вектор на матрицу нельзя, поэтому используется умножение на обратную (inverse) матрицу. Т означает транспонирование (transpose). В формуле подразумевается, что по умолчанию вектор — это столбец. Транспонирование превращает вектор-столбец в вектор-строку. При реализации на TensorFlow надо это учитывать, но в обратную сторону: по умолчанию вектор — это строка (одномерный тензор). На всякий случай: транспонирование — это не поворот на 90 градусов, это превращение строк в столбцы в том же порядке.

Итак, шаг метода Ньютона имеет следующий вид:

$\boldsymbol{p} \rightarrow \boldsymbol{p}-[\boldsymbol{H}loss_{p}(\boldsymbol{p})]^{-1}[\nabla_{p}loss(\boldsymbol{p})]$

В TensorFlow есть всё для реализации этого метода: 

# матрица Гессе для функции потерь по параметрам 
hess = tf.hessians(loss, p)[0]
# градиент превращаем в вектор-столбец
grad_col = tf.expand_dims(grad, -1)
# вычислим, насколько надо изменить вектор параметров
dp = tf.matmul(tf.linalg.inv(hess), grad_col)
# конвертируем вектор-столбец в вектор-строку
dp = tf.squeeze(dp)
# сдвигаем p на dp со знаком минус
newton = opt.apply_gradients([(dp, p)])
# запускаем цикл оптимизации, сделаем не больше 40000 шагов
session.run(tf.global_variables_initializer())
train(1e-10, 40000, loss, newton, feed_dict)
print('PARAMETERS:', session.run(p))

[Out]: step: 1, current loss: 3868.2291666666665
       step: 2, current loss: 105.04357496954218
       step: 4, current loss: 9.96663526704236
       ENDED ON STEP: 6, FINAL LOSS: 5.882202372519996e-20
       PARAMETERS: [2.5 2.5]

Хватило 6 шагов:


to4-eorewt5uoh_sdhokzqx0kuq.jpeg

Оптимизируем алгоритмом Гаусса-Ньютона

В методе Ньютона есть один недостаток — матрица Гессе. Благодаря TensorFlow мы можем посчитать её в одну строку кода. Согласно wiki, первое упоминание о своём методе Иоганн Карл Фридрих Гаусс сделал в 1809 году. Вычисление матрицы Гессе для нескольких параметров для метода наименьших квадратов тогда могло занять кучу времени. Сейчас можно считать, что алгоритм Гаусса-Ньютона использует аппроксимацию матрицы Гессе через матрицу Якоби, чтобы упростить вычисления. Но с точки зрения истории это не так: Людвиг Отто Гессе (который разработал матрицу, названную в честь него) родился 1811 году — через 2 года после первого упоминания алгоритма. А Карлу Густаву Якоби было 5 лет.

Алгоритм Гаусса-Ньютона не работает с функцией потерь. Он работает с функцией невязок $r(p)$. Эта функция принимает на вход вектор параметров $p$ и возвращает вектор невязок (residuals). В нашем случае, вектор $p$ состоит из 2 компонент (параметры $a$ и $b$ функции Розенброка), а вектор невязок из $m$ компонент (по числу экспериментов). Получается векторная функция векторного аргумента. Её производная:


Матрица Якоби (вектор вошёл — вектор вышел)

Чтобы не путаться в обилии символов, будем считать, что $\boldsymbol{J}_{r}$ — матрица Якоби функции невязок в текущей точке $\boldsymbol{p}$. Тогда алгоритм Гаусса-Ньютона можно записать так:

$\boldsymbol{p} \rightarrow \boldsymbol{p}-[\boldsymbol{J}_{r}^\intercal\boldsymbol{J}_{r}]^{-1}\boldsymbol{J}_{r}^\intercal r(\boldsymbol{p})$

Запись по форме полностью совпадает с записью метода Ньютона. Только вместо матрицы Гессе используется $\boldsymbol{J}_{r}^\intercal\boldsymbol{J}_{r}$, а вместо градиента $\boldsymbol{J}_{r}^\intercal r(\boldsymbol{p})$. Далее посмотрим, почему можно использовать такую аппроксимацию. Пока же приступим к реализации на TensorFlow: 

# к сожалению, в TensorFlow нет реализации вычисления матрицы Якоби,
# но мы помним, что эта матрица состоит из градиентов компонентов 
# функции невязок. В итоге, матрица вычисляется так:
# 1) разбиваем функцию невязок на отдельные компоненты tf.unstack(r)
# 2) для каждого компонента вычисляем градиент tf.gradients(r_i, p)
# 3) объединяем все градиенты в одну матрицу tf.stack
# такой способ вычисления не очень эффективный, более того мы очень 
# сильно увеличиваем размер графа потока данных
j = tf.stack([tf.gradients(r_i, p)[0]
              for r_i in tf.unstack(r)])
jT = tf.transpose(j)
# вектор невязок переводим в вектор-столбец
r_col = tf.expand_dims(r, -1)
# аппроксимация матрицы Гессе и градиента
hess_approx = tf.matmul(jT, j)
grad_approx = tf.matmul(jT, r_col)
# вычислим, насколько надо изменить вектор параметров
dp = tf.matmul(tf.linalg.inv(hess_approx), grad_approx)
# конвертируем вектор-столбец в вектор-строку
dp = tf.squeeze(dp)
# сдвигаем p на dp со знаком минус
ng = opt.apply_gradients([(dp, p)])
# запускаем цикл оптимизации, сделаем не больше 40000 шагов
session.run(tf.global_variables_initializer())
train(1e-10, 40000, loss, ng, feed_dict)

[Out]: step: 1, current loss: 3868.2291666666665
       step: 2, current loss: 14.653025157673625
       step: 4, current loss: 4.3918079172783016e-07
       ENDED ON STEP: 4, FINAL LOSS: 3.374364957618591e-17
       PARAMETERS: [2.5 2.5]

Хватило 4х шагов. Меньше чем для метода Ньютона.


iqns5lhcfqamjgy_f6isppnaxo4.jpeg

Как видно из кода, функция потерь не используется при оптимизации, только для критерия остановки и логирования. Как же оптимизационный алгоритм узнаёт, какую функцию минимизировать? Ответ удивителен: никак! Гаусс-Ньютон минимизирует только mean squared error.


Закрепим математическую часть статьи

Мы повторили всю необходимую нам математику. Давайте немного закрепим, чтобы дальше сконцентрироваться только на программировании и TensorFlow. Возможно вам понадобится карандаш, чтобы проследить последовательность математических действий.

Есть модель $y=f(x,p)$, где $x$ — вектор, $p$ — вектор параметров размерности $n$, а $y$ — скаляр. Из экспериментов получили $m$ точек $(x_{1}, y_{1}), ..., (x_{m}, y_{m})$ (data pairs). Векторная функция невязок зависит только от вектора параметров: $r(p) = (r_{1}(p), ... r_{m}(p))$, где $r_{k}(p) = y_{k} - \widehat{y_{k}} = y_{k} - f(x_{k}, p)$. Может возникнуть вопрос, почему функция невязок зависит только от $p$, ведь в формуле есть и $x_{k}, y_{k}$? Дело в том, что $x_{k}, y_{k}$ уже известны из экспериментов и зафиксированы, и только вектор параметров является переменным.

Надо найти такое значение вектора параметров $p$, чтобы сумма квадратов невязок (sum of squared error — sse или residual sum-of-squares — rss) была минимальная. Вместо привычной нам mse мы будем использовать sse, чтобы не таскать везде деление на $m$. Эти две функции достигают своего минимума на одном и том же векторе параметров. Наша функция потерь имеет вид:

$loss(p) = r_{1}^2(p) + \cdots + r_{m}^2(p) = \sum_{k=1}^{m} r_{k}^2(p)$

Далее для простоты нигде не будем писать $p$ в скобочках $(p)$.

Чтобы воспользоваться методом Ньютона, нам нужен градиент этой функции и матрица Гессе. Градиент — это вектор частных производных. Производная суммы — это сумма производных, а производная квадрата функции $r^2$ равна $2r \frac{\partial r}{\partial p}$. Таким образом получаем:

$\nabla_{p}loss = (\sum_{k=1}^{m}2r_{k}\frac{\partial r_{k}}{\partial p_{1}}, \cdots, \sum_{k=1}^{m}2r_{k}\frac{\partial r_{k}}{\partial p_{n}}) $

Матрица Гессе состоит из всех комбинаций вторых производных. Ячейка матрицы Гессе имеет такой вид:

$[\boldsymbol{H}loss_{p}]_{ij} = \frac{\partial^2 loss}{\partial p_{i} \partial p_{j}} = \sum_{k=1}^{m}(2\frac{\partial r_{k}}{\partial p_{i}}\frac{\partial r_{k}}{\partial p_{j}} + 2r_{k}\frac{\partial^2 r_{k}}{\partial p_{i}\partial p_{j}})$

Матрица эта диагонально симметричная. При её вычислении мы воспользовались теми же правилами, что и при вычислении градиента, плюс школьной формулой ${(uv)}'={u}'v+u{v}'$.
Отлично! Мы можем воспользоваться методом Ньютона.

Заметим, что в формуле для матрицы Гессе есть часть, которая уменьшается по мере приближения к минимуму, — слагаемое $2r_{k}\frac{\partial^2 r_{k}}{\partial p_{i}\partial p_{j}}$. Чем ближе к минимуму, тем меньше невязка, тем меньше $r_{k}$, а следовательно и всё произведение. И именно этот компонент доставляет больше всего вычислительных неприятностей — он содержит настоящую вторую производную. А что будет, если этот компонент откинуть? Мы получим алгоритм Гаусса-Ньютона.

Для него нужна матрица Якоби:

$\boldsymbol{J}_{r} = \begin{pmatrix} \frac{\partial r_{1}}{\partial p_{1}} & \cdots & \frac{\partial r_{1}}{\partial p_{n}}\\ \vdots & \ddots & \vdots \\ \frac{\partial r_{m}}{\partial p_{1}} & \cdots & \frac{\partial p_{m}}{\partial p_{n}} \end{pmatrix}$

Согласитесь, она выглядит гораздо проще, чем градиент и матрица Гессе. Заметим, что:

$2\boldsymbol{J}_{r}^\intercal\boldsymbol{J}_{r} \approx \boldsymbol{H}loss_{p}$

Заметить это «глазами» не так уж и просто. Лучше взять карандаш и начать производить матричное умножение (строка умножается на столбец). Тогда будет лучше видно, что результат — это выражение для матрицы Гессе без компонента $2r_{k}\frac{\partial^2 r_{k}}{\partial p_{i}\partial p_{j}}$, который мы и хотели откинуть.
С градиентом ещё проще (надо умножить матрицу на вектор):

$2\boldsymbol{J}_{r}^\intercal r = \nabla_{p}loss$

Таким образом, у нас сошлось утверждение, что алгоритм Гаусса-Ньютона — это метод Ньютона с аппроксимацией матрицы Гессе через матрицу Якоби, который работает только с функцией потерь mse.


Переформулируем и усложним задачу

Вымышленный прибор измерял координаты и высоту. Мы решили упростить железо, измерять только координаты, а высоту вычислять. Провели $m$ экспериментов $(x_{1}, y_{1}), ..., (x_{m}, y_{m})$, физики дали модель $y = rosenbrock(x, p)$. Разными методами оптимизации мы успешно нашли значение вектора параметров $p$, при котором модель работает достаточно точно.

А теперь представьте, что учёные отказались давать модель:»Мы совершенно не представляем физику вашего вымышленного мира. Давайте как-нибудь сами! ». В ситуациях, когда недостаточно знаний о природе вещей, или эта природа слишком сложная, можно попробовать решить задачу техникой машинного обучения (supervised learning). Например, с помощью нейронной сети. Такое решение имеет свою цену: данные (training set) — их надо много; сложность интерпретации — не всегда понятно почему модель (prediction model) работает или не работает; проблемы экстраполяции — модель стабильно работает только на данных из того же распределения, что и обучающая выборка.

В статье будем использовать многослойный перцептрон (multi-layer perceptron neural network или mlp). Есть аспекты, которые не будем учитывать в статье, но которые необходимо учитывать в реальной работе:


  • Начальные значения (starting values) весов. Будем инициализировать веса алгоритмом Xavier’a, но наверняка есть и более удачные варианты.
  • Переобучение (overfitting). Тема статьи — методы оптимизации. А их задача искать минимум, несмотря ни на что. А задача разработчика — не допустить переобучения.
  • Масштабирование входных данных (scaling of the input). Не будем делать, но это может значительно улучшить результат.

В прошлый раз хватило 9 экспериментов. В этот раз проведём 500: 

# эксперименты будут случайными
def get_random_rosenbrock_data_points(m):
    result = np.zeros((m, 3))
    result[:, 0] = np.random.uniform(-2, 2, m)
    result[:, 1] = np.random.uniform(-2, 2, m)
    result[:, 2] = rosenbrock(result[:, 0], result[:, 1], 2.5, 2.5)
    return result

m = 500
data_points = get_random_rosenbrock_data_points(m)
# overfitting не тема статьи, но совсем без валидации нельзя
validation_data_points = get_random_rosenbrock_data_points(m)

Из экспериментов получили 500 точек. Наша задача — ничего не зная о функции Розенброка обучить модель (learner), которая будет прогнозировать высоту (outcome measurement) по координатам (features) для новых ещё неизвестных данных.


Реализуем однослойный перцептрон

Часто структуру перцептрона рисуют в виде кружочков со стрелочками (network diagram). Мне больше нравится такой вариант визуализации от MatLab:


npntxk7brhnpkk7p5v1p4bjur-a.png

Вектор из двух компонент подаётся на вход (input). Он умножается на матрицу весов $W$ (weights) размера 2×10, к результату прибавляется вектор смещения $b$ (bias) размера 10, и результат передаётся нелинейной функции активации (activation). Так получается первый (единственный) скрытый уровень (hidden layer) из 10 нейронов. Во втором блоке к нему применяется другая матрица весов, другой вектор смещения, и уже без нелинейной функции получается скалярный результат (output).

Самая удобная запись, на мой взгляд, такая (будем использовать $tanh$):

$\begin{matrix} h_{1} = tanh(xW_{1} + b_{1})\\ \widehat{y} = h_{1}W_{2} + b_{2} \end{matrix} $

Или подробнее:

$h_1 = tanh(\begin{bmatrix} x_1 & x_2 \end{bmatrix}\begin{bmatrix} w^{(1)}_{1,1} & \cdots& w^{(1)}_{1,10} \\ w^{(1)}_{2,1} &\cdots& w^{(1)}_{2,10} \end{bmatrix} + \begin{bmatrix} b^{(1)}_1 & \cdots & b^{(1)}_{10} \end{bmatrix}) \\ \widehat{y} = \begin{bmatrix}h^{(1)}_1 & \cdots & h^{(1)}_{10}\end{bmatrix} \begin{bmatrix} w^{(2)}_{1,1} \\ \vdots \\ w^{(2)}_{1,10} \\ \end{bmatrix} + b_2 $

Здесь используются матричные операции. Размерность $W_{1}$ однозначно определяется размерностью входного вектора и желаемой «шириной» скрытого уровня $h_{1}$, а из неё однозначно определяется размерность вектора-столбца $W_{2}$. Получили 41 параметр. Регулируя количество скрытых слоёв и их ширину, можно изменять количество параметров.

Если входной вектор заменить на матрицу $m \times 2$, то за один проход можно вычислить прогнозы для всех точек. Результатом буде вектор-столбец $\widehat{y}$ из $m$ компонентов: 

# хотим скрытый слой из 10 "нейронов"
n_hidden = 10
# начальную догадку будем конструировать алгоритмом Xavier'a
initializer = tf.contrib.layers.xavier_initializer()

# тензоры данных экспериментов
x = tf.placeholder(tf.float64, shape=[m, 2])
y = tf.placeholder(tf.float64, shape=[m, 1])

# тензоры весов и смещения для вычисления скрытого слоя
W1 = tf.Variable(initializer([2, n_hidden], dtype=tf.float64))
b1 = tf.Variable(initializer([1, n_hidden], dtype=tf.float64))
# вычисление скрытого слоя, используем tanh для активации
h1 = tf.nn.tanh(tf.matmul(x, W1) + b1)
# тензоры весов и смещения для вычисления прогноза
W2 = tf.Variable(initializer([n_hidden, 1], dtype=tf.float64))
b2 = tf.Variable(initializer([1], dtype=tf.float64))
# вычисление прогноза
y_hat = tf.matmul(h1, W2) + b2
# невязки
r = y - y_hat
# также используем mse в качестве функции потерь
loss = tf.reduce_mean(tf.square(r))

# данные для подстановки в тензоры placeholder
feed_dict = {x: data_points[:,0:2], 
             y: data_points[:,2:3]}
validation_feed_dict = {x: validation_data_points[:,0:2], 
                        y: validation_data_points[:,2:3]}


Обучим нейросеть методом Adam

Обучение нейронной сети оптимизатором Adam ничем не отличается от поиска параметров для модели $rosenbrock$. Добавим подсчёт mse для валидационной выборки в конце: 

# сконструируем шаг оптимизации
adm = tf.train.AdamOptimizer(1e-2).minimize(loss)
session.run(tf.global_variables_initializer())
# запускаем цикл оптимизации, сделаем не больше 40000 шагов
train(1e-10, 40000, loss, adm, feed_dict)
print('VALIDATION LOSS: '+str(session.run(loss, validation_feed_dict)))

[Out]: step: 1, current loss: 671.4242576535694
       [...]
       ENDED ON STEP: 40000, FINAL LOSS: 0.22862158574440725
       VALIDATION LOSS: 0.29000289644978866

У вас могут получиться совершенно другие числа. Они во многом зависят от удачи: случайная инициализация весов, тренировочных данных, данных для валидации.


Вычислим матрицу Якоби для нейронной сети

Для модели $rosenbrock$ матрица Якоби вычислялась в 2 строки кода. Теперь ситуация сильно усложнилась:


  • Данные. Раньше было 9 эксперимен

    © Habrahabr.ru