Преобразование Лежандра со всех сторон

Учась в институте, я встречал преобразование Лежандра в совершенно неожиданных местах: от методов оптимизации до физики. Было очевидно, что существует некоторая фундаментальная структура, сшивающая все воедино. В этой статье пробую нащупать эту нить: показываю, как преобразование Лежандра проявляется в различных областях, стараюсь идти от конкретных примеров и делать акцент на связях — чтобы выстроить цельную картину. Приветствую любые комментарии и указания на ошибки/неточности

Геометрия

Для начала рассмотрим функцию $f(x) = x^{2}$ . В точке $x_{0}$ касательная: $h(x) = f\left( x_{0} \right) + f^{'}\left( x_{0} \right)\left( x - x_{0} \right)$

$x^{2}$ – типичный пример гладкой выпуклой функции. Формализуем интуицию вокруг выпуклости:

Во-первых, выпуклая функция – это та, чей надграфик выпуклый, т.е. содержит любую хорду. Двигаясь между точками $\left( x,f(x) \right)$ и $\left( y,f(y) \right)$ , мы всегда оказываемся над графиком. Это называют неравенством Йенсена:

$f \left( tx + (1 - t)y \right) \leq t \cdot f(x) + (1 - t) \cdot f(y)\ \forall x,y\ \forall t \in \lbrack 0,1\rbrack$

Во-вторых, для гладких функций выпуклость можно определить так: график функции лежит над касательными:

$f(y) \geq f(x) + f^{'}(x)(y - x)\ \ \forall x,y$

В-третьих, геометрически это означает, что наклон f растет быстрее, чем наклон касательной: $f^{''}(x) \geq 0\ \ \forall x$

Остановимся на (2). Заметим, что функцию можно описать двумя способами:

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

Оказывается, для выпуклых функций оба описания эквивалентны. Пока фиксируем, а чем это полезно — ниже

Алгебра

Рассмотрим некоторую хитрую конструкцию. Возьмем функцию и сопоставим ей функцию $f^{*}(p)$ по такому правилу (пока примем это как черный ящик, постепенно поймем глубинные смыслы):

$f^{*}(p) = \sup_{x}\left( px - f(x) \right)$

Посмотрим на примерах:

$f(x) \equiv 0\ \ \Rightarrow \ \ f^{*}(p) = \sup_{x}(px) = \left\{ \begin{matrix} 0,\ \ p = 0 \\ \infty,\ \ p \neq 0 \end{matrix} \right.$

Получилась некоторая «вывернутая дельта-функция»: во всех точках $\infty$ , в нуле 0. Попробуем сдвинуть исходную функцию f по вертикали:

$f(x) \equiv b\ \ \Rightarrow \ \ f^{*}(p) = \sup_{x}(px - b) = \left\{ \begin{matrix} - b,\ \ p = 0 \\ \infty,\ \ p \neq 0 \end{matrix} \right.$

Наша конструкция $f^{*}(p)$ просела на некоторый сдвиг, но все еще центрирована в . Усложним пример:

$f(x) = x\ \ \Rightarrow \ \ f^{*}(p) = \sup_{x}\left( (p - 1)x \right) = \left\{ \begin{matrix} 0,\ \ p = 1 \\ \infty,\ \ p \neq 1 \end{matrix} \right.$

Получилась та же «вывернутая дельта-функция», но теперь с центром в 1. Видно, что этот пример отличается от самого первого в единственном месте: множитель перед x в первом супремуме – просто p, а в текущем – (p-1). Эта единица возникла именно как коэффициент перед x, то есть наклон функции f. Обобщим и проверим интуицию:

$f(x) = kx + b\ \ \Rightarrow \ \ f^{*}(p) = \sup_{x}\left( (p - k)x - b \right) = \left\{ \begin{matrix} - b,\ \ p = k \\ \infty,\ \ p \neq k \end{matrix} \right.$

Видно, что наша конструкция $f^{*}(p)$ детектирует наклон исходной функции f(x). Когда аргумент p функции $f^{*}(p)$ совпадает с наклоном исходной функции f(x), индикатор срабатывает, и мы получаем конечное число на выходе. Более того, это число совпало с величиной, на которую нужно сдвинуть функцию вниз, чтобы «обернуть» исходную f снизу (и тогда $\infty$ при $p \neq k$ логична, ведь прямые и гарантированно пересекутся). Эта идея собрана в неравенстве Фенхеля — Юнга:

$f^{}(p) = \sup_{x}\left( px - f(x) \right)\ \ \Rightarrow \ \ f(x) \geq px - f^{*}(p)\ \ \forall x,p$

Прямые $px - f^{*}(p)$ целиком лежат под графиком – это как раз набор касательных, оборачивающих заданную функцию f снизу. Продолжим эксперименты: что если наклон меняется со временем?

$f(x) = x^{2}\ \ \Rightarrow \ \ f^{*}(p) = \sup_{x}\left( - x^{2} + px \right)$

Внутри супремума – парабола ветвями вниз. Теперь связалось, почему мы вообще рассматриваем выпуклые функции: супремум искать проще. В данном случае супремум достигается в вершине параболы $\frac{p}{2}$ . Подставим $x = \frac{p}{2}$ , тогда

$\ f^{*}(p) = px - f\left( \frac{p}{2} \right) = p*\frac{p}{2} - \left( \frac{p}{2} \right)^{2} = \frac{p^{2}}{4}$

Нашей исходной параболе $f(x) = x^{2}$ сопоставляется парабола $f^{*}(p) = \frac{p^{2}}{4}$ , лежащая ниже исходной. А что если провернуть наше преобразование в обратную сторону?

$f(x) = \frac{x^{2}}{4}\ \ \Rightarrow \ \ f^{*}(p) = \sup_{x}\left( - \frac{x^{2}}{4} + px \right)$

Аналогичные рассуждения про вершину параболы приведут нас к тому, что $f^{*}(p) = p^{2}$ . Получается, конкретно в этом примере мы провели преобразование дважды и вернулись к исходной функции $f^{**} = f$ .

Всегда ли это так? Оказывается, нет: только для выпуклых и замкнутых – это теорема Фенхеля-Моро:

$f^{**} \leq f$ всегда. Фенхель-Юнг: $f(x) \geq px - f^{*}(p)\ \ \forall x,p\ \ \Rightarrow \ \ f(x) \geq \sup_{p}{px - f^{*}(p)} = f^{**}(x)$
$f^{**}$ всегда выпукла и замкнута, по построению
довольно техническое доказательство, по которому $f^{**} = conv\ f$ (выпуклая оболочка f, т.е. множество всех хорд с концами в надграфике). По определению через выпуклый надграфик, выпуклая функция совпадает со своей выпуклой оболочкой, поэтому для выпуклых $f^{**} = f$

Заметим отдельно, что, опираясь на выпуклость, мы подбирали супремум по нулям производной: $f^{*}(p) = \left( px - f(x) \right)\left. \ \right|_{p = f^{'}(x)} = \left( px - f(x) \right)\left. \ \right|_{x = \left( f^{'} \right)^{- 1}(p)}$ . Такая замена требует обратимости $f^{'}$ , т.е. строгой выпуклости – это нам позже пригодится.

Подумаем, что именно делает наше преобразование. По графикам ранее было видно, что кривые как будто отражаются относительно некоторой границы между ними – что это за граница? Угадаем неподвижную точку:

$f(x) = \frac{x^{2}}{2},\ \ f^{*}(p) = \sup_{x}\left( - \frac{x^{2}}{2} + px \right) = \left( - \frac{x^{2}}{2} + px \right)\left. \ \right|_{p = f^{'}(x)} = \frac{p^{2}}{2}$

Предположим, что существует другая неподвижная точка $g = g^{*},\ \ g \neq f$ . Фенхель-Юнг: $px \leq g(x) + g^{*}(p) = g(x) + g(p)$ для $\forall x,p$ . В частности, для получим $\frac{x^{2}}{2} < g(x)$ . Но тогда:

$g^{*}(p) = \sup_{x}\left( px - g(x) \right) < \sup_{x}\left( px - \frac{x^{2}}{2} \right) = \frac{p^{2}}{2}$

т.е. $g \neq g^{*}$ , противоречие. Значит, $\frac{x^{2}}{2}$ – единственная неподвижная точка. По сути, проводя наше преобразование, мы «отражали» функции относительно «зеркала» $\frac{x^{2}}{2}$ .

Конструкция, которую мы все это время рассматривали, называется преобразованием Лежандра. Ключевая идея: выпуклая функция = верхняя огибающая своих нижних оценок.

Физика / механика

Существует три классических подхода для описания физических систем: по Ньютону, Лагранжу и Гамильтону. Ньютон работает с силами – нам это сейчас не интересно. Сфокусируемся на лагранжевом и гамильтоновом формализме, покажем связь между ними и место, где возникает преобразование Лежандра.

Рассмотрим задачу кинематики: по какой траектории полетит камень в поле гравитации (для простоты – бросим вертикально)?

Можно ввести координату , потенциальную энергию , кинетическую энергию $T = \frac{mv^{2}}{2} = \frac{m{\dot{q}}^{2}}{2}$ . В контексте выпуклых функций сразу бросается в глаза выпуклость кинетической энергии по скоростям (спойлер: это пригодится).

Полная энергия сохраняется, т.е. $0 = \frac{d}{dt}\left( mgq + \frac{m{\dot{q}}^{2}}{2} \right) = mg\dot{q} + m\dot{q}\ddot{q} = m\dot{q}(g + \ddot{q})$ , причем масса положительна, а скорость зануляется только в одной точке (вершине траектории). Остается $\ddot{q} = - g$ . Решив этот дифур, получим знакомое $q = q_{0} + v_{0}t - \frac{gt^{2}}{2}$ .

Вместо скорости $\dot{q}$ , можно описать систему через изменение импульса: $\frac{dp}{dt} = - mg,\ \ p = m\dot{q} = mv_{0} - mgt$ . Сократив, получим $\dot{q} = v_{0} - gt$ . Решив этот дифур, получим такой же ответ: $q = q_{0} + v_{0}t - \frac{gt^{2}}{2}$ .

Вывод: два описания системы (координаты и скорости vs координаты и импульсы) эквивалентны в рассмотренной нами выпуклой задаче.

Оказывается, оба подхода опираются на более глубокую теорию и тесно связаны

Лагранж описывает систему с помощью обобщенных координат и скоростей $\dot{q}$

Вводится кинетическая энергия $T\left( q,\dot{q} \right)$ и потенциальная энергия $U\left( q,\dot{q} \right)$ , их разность называется лагранжианом . Интеграл лагранжиана по времени (вдоль траектории) называется действием:

$S\lbrack q\rbrack = \int_{t_{1}}^{t_{2}}{L\left( q,\dot{q},t \right)dt}$

Это функционал, сопоставляющий траектории некоторое число. Постулируется Принцип наименьшего действия: природа стремится минимизировать этот функционал, поэтому траектории точек – экстремали. Почему это верно – отдельный глубокий вопрос.

Возмутим траекторию (требования к возмущению – гладкая функция, обнуляющаяся на концах):

$q_{\varepsilon}(t) = q(t) + \varepsilon \cdot \eta(t),\ \ {\dot{q}}_{\varepsilon}(t) = \dot{q}(t) + \varepsilon \cdot \dot{\eta}(t),\ \ \eta\left( t_{1} \right) = \eta\left( t_{2} \right) = 0$

Действие вдоль возмущенной траектории:

$S\left\lbrack q_{\varepsilon} \right\rbrack = S(\varepsilon) = \int_{t_{1}}^{t_{2}}{L\left( q + \varepsilon\eta,\dot{q} + \varepsilon\dot{\eta},t \right)dt}$

Проварьируем: продифференцируем по $\varepsilon$ , найдем экстремум:

$0 = \frac{\delta S}{\delta q} = \frac{dS}{d\varepsilon} = \int_{t_{1}}^{t_{2}}{\frac{dL\left( q + \varepsilon\eta,\dot{q} + \varepsilon\dot{\eta},t \right)}{d\varepsilon}dt} = \int_{t_{1}}^{t_{2}}{\left\lbrack \frac{\partial L}{\partial q}\frac{\partial q_{\varepsilon}}{\partial\varepsilon} + \frac{\partial L}{\partial\dot{q}}\frac{\partial{\dot{q}}_{\varepsilon}}{\partial\varepsilon} \right\rbrack dt} = \int_{t_{1}}^{t_{2}}{\left\lbrack \frac{\partial L}{\partial q}\eta \right\rbrack dt} + \int_{t_{1}}^{t_{2}}{\left\lbrack \frac{\partial L}{\partial\dot{q}}\dot{\eta} \right\rbrack dt}$

Второе слагаемое – проинтегрируем по частям, вспомнив, что $\eta\left( t_{1} \right) = \eta\left( t_{2} \right) = 0$ :

$0 = \frac{dS}{d\varepsilon} = \int_{t_{1}}^{t_{2}}{\left\lbrack \frac{\partial L}{\partial q}\eta \right\rbrack dt} + \left\lbrack \frac{\partial L}{\partial\dot{q}}\eta \right\rbrack_{t_{1}}^{t_{2}} - \int_{t_{1}}^{t_{2}}{\frac{d}{dt}\left( \frac{\partial L}{\partial\dot{q}} \right)\eta dt} = \int_{t_{1}}^{t_{2}}{\left\lbrack \frac{\partial L}{\partial q} - \frac{d}{dt}\left( \frac{\partial L}{\partial\dot{q}} \right) \right\rbrack\eta dt}$

Мы рассматривали выражения для произвольного $\eta$ , а по основной лемме вариационного исчисления это значит, что квадратная скобка тождественно равна нулю. Получим уравнение Эйлера-Лагранжа:

$\frac{\partial L}{\partial q} - \frac{d}{dt}\left( \frac{\partial L}{\partial\dot{q}} \right) = 0$

В лагранжевом формализме вводится также обобщенный импульс $p = \frac{\partial L}{\partial\dot{q}}$ . Можно подставить систему из школьной механики и проверить соответствие с привычным определением: $p = \frac{\partial}{\partial v}\left( \frac{1}{2}mv^{2} - U(x) \right) = mv$

Гамильтон описывает систему с помощью обобщенных координат и импульсов

Вводится гамильтониан $H = p\dot{q} - L$ . Сразу интересное наблюдение: это выражение уже похоже на содержимое супремума в преобразовании Лежандра, но подойдем к этой связи осторожнее. Пока отметим: в классической механике типично, что кинетическая энергия $T = \frac{mv^{2}}{2} = \frac{p\dot{q}}{2}$ , поэтому $H = p\dot{q} - L = 2T - (T - U) = T + U$ , т.е. полная энергия.

Распишем и проварьируем действие в терминах гамильтониана (дважды и независимо по q,p):

$S\lbrack q,p\rbrack = \int_{t_{1}}^{t_{2}}{\left\lbrack p\dot{q} - H(q,p) \right\rbrack dt},\ \ 0 = \frac{dS}{d\varepsilon} = \int_{t_{1}}^{t_{2}}{\left\lbrack \dot{q} - \frac{\partial H}{\partial p} \right\rbrack\xi dt} - \int_{t_{1}}^{t_{2}}{\left\lbrack \dot{p} + \frac{\partial H}{\partial q} \right\rbrack\eta dt}$

Получим уравнения Гамильтона:

$\dot{q} = \frac{\partial H}{\partial p},\ \ \dot{p} = - \frac{\partial H}{\partial q}$

Оказывается, Гамильтон – это лежандров образ Лагранжа. Рассмотрим $H = \sup_{\dot{q}}\left( p\dot{q} - L \right)$ . На самом деле, мы так и вводили выше, но не заметили, как неявно учли, что находимся в супремуме выражения. Действительно,

$0 = \frac{d}{d\dot{q}}\left( p\dot{q} - L \right) = p - \frac{\partial L}{\partial\dot{q}}$ , т.е. $p = \frac{\partial L}{\partial\dot{q}}$

Мы получили буквально определение импульса, только теперь оно не постулируется, а следует из поиска супремума.

Теперь продифференцируем и подставим определение импульса:

$dH = \dot{q}dp + pd\dot{q} - \frac{\partial L}{\partial q}dq - \frac{\partial L}{\partial\dot{q}}d\dot{q} - \frac{\partial L}{\partial t}dt = \dot{q}dp - \frac{\partial L}{\partial q}dq - \frac{\partial L}{\partial t}dt + \left( p - \frac{\partial L}{\partial\dot{q}} \right)d\dot{q} = \dot{q}dp - \frac{\partial L}{\partial q}dq - \frac{\partial L}{\partial t}dt$

При этом сам . Сопоставим частные производные:

$dH = \frac{\partial H}{\partial p}dp + \frac{\partial H}{\partial q}dq + \frac{\partial H}{\partial t}dt\ \ \ \Rightarrow \ \ \dot{q} = \frac{\partial H}{\partial p},\ \ \frac{\partial H}{\partial q} = - \frac{\partial L}{\partial q}$

Вспомним ур-е Эйлера-Лагранжа – на решении:

$\frac{\partial L}{\partial q} = \frac{d}{dt}\left( \frac{\partial L}{\partial\dot{q}} \right) = \dot{p},\ \ т.е.\dot{p} = - \frac{\partial H}{\partial q}$

Видна первая большая польза от преобразования Лежандра: мы свели ОДУ 2го порядка к паре ОДУ 1го порядка. Плюс существуют методы, которые проще использовать в гамильтоновом формализме (а в квантовой механике в принципе вводится именно оператор Гамильтона), поэтому преобразование Лежандра возникает часто.

Почему это вообще работает? И пара комментариев

Преобразование Лежандра корректно определено, если отображение $\dot{q} \rightarrow p$ локально обратимо (мы ведь выражаем $\dot{q}$ через из определения импульса $p = \frac{\partial L}{\partial\dot{q}}$ ). Достаточное условие – невырожденность гессиана $\frac{\partial^{2}L}{\partial{\dot{q}}^{2}}$ , но более физично условие положительной определенности – т.е. выпуклости лагранжиана по $\dot{q}$ . Почему более физично, на простом примере: $L = \frac{mv^{2}}{2} - U(x),\ \ \ \frac{\partial^{2}L}{\partial{\dot{q}}^{2}} = m > 0$ . Оказывается, гессиан лагранжиана описывает массы и их аналоги, что часто оказывается положительным. Сразу отметим: $p = \frac{\partial L}{\partial\dot{q}}$ и $\dot{q} = \frac{\partial H}{\partial p}$ взаимно обратные отображения, поэтому гессианы и – обратные матрицы (из теоремы об обратной функции). В частности, гессиан описывает величины с размерностью обратных масс (подвижностей)

Есть более глубокий взгляд: постулируется ряд симметрий, из которых выводятся суждения о виде

однородность пространства $\rightarrow T$ зависит только от $\dot{q}$
изотропность пространства $\rightarrow T$ квадратичная форма от $\dot{q}$
из принципа относительности Галилея следует, что законы механики одинаковы в любой ИСО, а значит может зависеть только от координат

Получается, что выпуклость часто гарантируется на фундаментальном уровне. Но не всегда: иногда физика допускает невыпуклые лагранжианы, и тогда эквивалентность Лагранжа vs Гамильтона теряется.

Дополнительно рассмотрим случай, когда какой-то $\dot{q}$ нельзя выразить через свой из $p = \frac{\partial L}{\partial\dot{q}}$ (как мы делали). Это означает $\frac{dp}{d\dot{q}} = \frac{\partial^{2}L}{\partial{\dot{q}}^{2}} = 0$ . Дважды проинтегрируем: $L = C_{1}(q) \cdot \dot{q} + C_{2}(q)$ . Подставим:

$H = \sup_{\dot{q}}\left( p\dot{q} - L \right) = \sup_{\dot{q}}\left( \left( p - C_{1}(q) \right) \cdot \dot{q} - C_{2}(q) \right) = \left\{ \begin{matrix} - C_{2}(q),\ \ p = C_{1}(q) \\ \infty,\ \ p \neq C_{1}(q) \end{matrix} \right.$

Это буквально наш старый пример

$f(x) = kx + b\ \ \Rightarrow \ \ f^{}(p) = \sup_{x}\left( (p - k) \cdot x - b \right) = \left\{ \begin{matrix} - b,\ \ p = k \\ \infty,\ \ p \neq k \end{matrix} \right.$

Получившийся результат означает, что энергия конечна лишь при специфичном условии на импульсы, т.е. импульсы не меняются свободно. Это называется первичной связью, и в своей основе – это лежандров образ линейной по скорости части лагранжиана. По названию ясно, что есть и вторичные — они устроены сложнее

Задачу можно свести к еще более удобному виду, если перейти в координаты, в которых гамильтониан $\equiv 0$ . Решение нам сейчас не интересно, в итоге получится уравнение Гамильтона-Якоби (пригодится, но не скоро, S — вдоль экстремали):

$\frac{\partial S}{\partial t} + H\left( q,\frac{\partial S}{\partial q},t \right) = 0,$

Итог: преобразование Лежандра позволяет переписать задачу в виде, в котором ее иногда проще решать.

Оптимизация

В этом разделе очень много применений преобразования Лежандра, сфокусируемся на нескольких ключевых. Ключевая задача всего раздела — найти минимум заданной функции.

Чаще преобразование Лежандра записывают в векторном виде: $f^{*}(p) = \sup_{x}\left( \left\langle x,p \right\rangle - f(x) \right)$ . Для гладких функций: $f^{*}(p) = \left\langle x,p \right\rangle - f(x)$ при $x = (\nabla f)^{- 1}(p)$ . Отсюда же следует: $\nabla f^{*} = (\nabla f)^{- 1}$ . Кроме того, $\nabla^{2}f^{*} = \left( \nabla^{2}f \right)^{- 1}$ (если гессиан обратим — по теореме об обратной функции, что мы пронаблюдали на примере лагранжиана и гамильтониана)

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

Условная оптимизация

Рассмотрим задачу $\min_{x}{f(x)}$ при . Составляем лагранжиан (сразу отметим: это не тот же лагранжиан, который был в физике, к сожалению, тут коллизия имен – лишь совпадение):

$L(x,\lambda) = f(x) + \lambda^{T}(b - Ax)$

Что произошло: для допустимых x (т.е. таких, что второе слагаемое зануляется) можно взять любой множитель $\lambda$ – все равно в скобке ноль. Но теперь давайте рассмотрим не только допустимые x, а вообще все x: свободы больше, значит минимум только опустится:

$\inf_{допустимые\ x}{f(x)} = \inf_{допустимые\ x}\left( f(x) + \lambda^{T}(b - Ax) \right) \geq \inf_{x}\left( f(x) + \lambda^{T}(b - Ax) \right)$

Обозначим $g(\lambda) = \inf_{x}\left( f(x) + \lambda^{T}(b - Ax) \right)$ – это семейство оценок снизу на нашу функцию f. Из этого семейства логично взять лучшую: $\sup_{\lambda}{g(\lambda)}$ – это все еще оценка снизу, но максимальная (наиболее «тугая»).

Перейдя к двойственной задаче, мы буквально провели преобразование Лежандра. Оказывается, преобразование Лежандра всегда было задачей условной оптимизации с линейными ограничениями:

$g(\lambda) = \inf_{x}{L(x,\lambda)} = \lambda^{T}b + \inf_{x}\left\lbrack f(x) - \left( A^{T}\lambda \right)^{T}x \right\rbrack = \lambda^{T}b - \sup_{x}\left\lbrack \left( A^{T}\lambda \right)^{T}x - f(x) \right\rbrack = \lambda^{T}b - f^{*}\left( A^{T}\lambda \right)$

Поскольку все время $g(\lambda)$ по построению были оценками снизу, $\sup_{\lambda}g \leq \inf f$ . Причем в невыпуклом случае зазор ожидается, потому что двойственная задача восстанавливает лишь $f^{**} = conv\ f$ .

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

Часто минимум функции нельзя/сложно/долго искать аналитически, либо нас интересует любой достаточно хороший локальный минимум. Тогда применяют итеративные алгоритмы, в частности, популярный в ML градиентный спуск. Ключевая идея – подгоняем аргумент функции в ту сторону, куда функция убывает:

$x_{k + 1} = x_{k} - \eta \cdot \nabla f\left( x_{k} \right)$

Эту инженерную интуицию можно обосновать аналитически. Аппроксимируем функцию f ее касательной в текущей точке. Попробуем перейти в ее минимум, но минимум касательной, ясно, на бесконечности, поэтому добавим штраф за отход от исходной строчки – норму шага с регулируемым коэффициентом $\eta$ :

$f(x) \approx f\left( x_{k} \right) + \left\langle \nabla f\left( x_{k} \right),\ x - x_{k} \right\rangle$ $x_{k + 1} = \arg{\min\left\lbrack f\left( x_{k} \right) + \left\langle \nabla f\left( x_{k} \right),\ x - x_{k} \right\rangle + \frac{1}{2\eta}\left\| x - x_{k} \right\|^{2} \right\rbrack} = \arg{\min_{x}\left\lbrack \eta\left\langle \nabla f\left( x_{k} \right),x \right\rangle + \frac{1}{2}\left\| x - x_{k} \right\|^{2} \right\rbrack}$

Видно, что градиентный спуск эквивалентен задаче оптимизации с евклидовым штрафом $\left\| x - x_{k} \right\|^{2}$

Продифференцируем выражение под аргминимумом по x:

$\nabla f\left( x_{k} \right) + \frac{1}{\eta}\left( x - x_{k} \right) = 0\ \ \Rightarrow \ \ x_{k + 1} = x_{k} - \eta \cdot \nabla f\left( x_{k} \right)$

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

Сопряженное пространство — это множество всех непрерывных линейных функционалов, заданных на исходном пространстве. Преобразование Лежандра переводило нас именно туда. В обоих пространствах существует метрика — симметричная положительно определенная матрица , на основе которой вводится скалярное произведение: $\left\langle u,v \right\rangle_{M} = u^{T}Mv$ . В сопряженном пространстве живет $\partial f$ – “сырой” ковектор-градиент, а в исходном — $\nabla f$ , причем $\left\langle \nabla f,v \right\rangle_{M} = \partial f(v)$ , или $\nabla f = M^{- 1}\partial f$ . В Евклидовой метрике , поэтому $\nabla f = \partial f$ , как мы и писали в формуле градиентного спуска.

Мы ранее уже сталкивались с ситуацией, когда сопряжение не изменяет функцию: $\psi = \frac{1}{2}x^{2}$ . Вернемся к ней, а точнее, рассмотрим $\psi(x) = \frac{1}{2}x^{T}Mx$ . Тогда $\nabla\psi(x) = Mx,\ \ \ \nabla\psi^{*}(p) = M^{- 1}p$ – это буквально наше правило $\nabla\psi^{*} = (\nabla\psi)^{- 1}$ из начала раздела. Как видно, сама матрица $M = \nabla^{2}\psi$ . В физике у нас уже возникал обратимый гессиан $\frac{\partial^{2}L}{\partial{\dot{q}}^{2}} = m > 0$ , масса/матрица инерции. Это та самая матрица , просто в другом контексте.

Возвращаемся к градиентному спуску. В более общем случае $M \neq I$ будет:

$x_{k + 1} = x_{k} - \eta \cdot \nabla f\left( x_{k} \right) = x_{k} - \eta \cdot \nabla\psi^{*}(\partial f)$

Просто мы не заметили, как неявно выбрали самосопряженную $\psi$ – ту самую из начала статьи.

Зеркальный спуск

Суммируем, в чем была проблема. Мы складываем объект из исходного пространства (x) с объектом из сопряженного (градиент). Мы незаметно отождествили пространства и $X^{*}$ , хотя в общем случае они различны. На практике это означает, что мы не учли геометрию задачи. Классический пример – минимизация на вероятностном симплексе: $x_{i} \geq 0,\ \ \sum_{}^{}x_{i} = 1$ . Евклидов штраф на нем всё время утыкается в границы, так что хотелось бы заложить в метод учет геометрии пространств

Эту проблему решает алгоритм зеркального спуска:

берем строго выпуклую гладкую функцию $\psi$ (distance-generating function, «зеркало»)
переводим точку в сопряженное пространство
делаем градиентный спуск в сопряженном пространстве – там, где на самом деле живут градиенты
возвращаемся в исходное пространство

Именно за счет переходов между исходным пространством и сопряженным алгоритм получил название «зеркальный» спуск:

$x_{k + 1} = \nabla\psi^{*}\left( \nabla\psi\left( x_{k} \right) - \eta \cdot \nabla f\left( x_{k} \right) \right)$

Градиентный спуск – частный случай зеркального при $\psi(x) = \frac{1}{2}x^{T}x$ .

У зеркального спуска есть эквивалентная запись (можно убедиться, взяв градиент выражения под argmax):

$x_{k + 1} = \arg{\min_{x}\left\lbrack \eta\left\langle \nabla f\left( x_{k} \right),x \right\rangle + D_{\psi}\left( x,x_{k} \right) \right\rbrack},\ \ где\ D_{\psi}(x,y) = \psi(x) - \psi(y) - \left\langle \nabla\psi(y),x - y \right\rangle$

Функция $D_{\psi}$ называется дивергенцией Брэгмана – по сути, это аналог расстояния. Для $\psi = \frac{1}{2}\left\| \ \cdot \ \right\|^{2}$ будет $D_{\psi}(x,y) = \frac{1}{2}\left\| x - y \right\|^{2}$ , и вся конструкция схлопывается в обычный градиентный спуск.

Вернемся к примеру с минимизацией на вероятностном симплексе $x_{i} \geq 0,\ \ \sum_{}^{}x_{i} = 1$ . Мы упоминали, что евклидов штраф не учитывает геометрию пространства, но какой подошел бы лучше? Оказывается, в качестве зеркала лучше всего подходит отрицательная энтропия $\psi = \sum_{}^{}{x_{i}\log x_{i}}$ .

$D_{\psi}(x,y) = \psi(x) - \psi(y) - \left\langle \nabla\psi(y),x - y \right\rangle = \sum_{}^{}{x_{i}\log x_{i}} - \sum_{}^{}{y_{i}\log y_{i}} - \sum_{}^{}{\left( \log y_{i} + 1 \right)\left( x_{i} - y_{i} \right)} = \sum_{}^{}{x_{i}\log x_{i}} - \sum_{}^{}{x_{i}\log y_{i}} - \sum_{}^{}x_{i} + \sum_{}^{}y_{i} = \sum_{}^{}{x_{i}\log\frac{x_{i}}{y_{i}}} = KL\left( x||y \right)$

Оказывается, дивергенция Брэгмана для отрицательной энтропии становится KL-дивергенцией – естественным аналогом расстояния для распределений.

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

Линейное приближение очень грубое, и возникает идея приближать квадратично. Ключевая идея: на каждом шаге строим касательную параболу, и переходим в ее минимум:

$f(x) \approx f\left( x_{k} \right) + \left\langle \nabla f\left( x_{k} \right),\ x - x_{k} \right\rangle + \frac{1}{2}\left( \ x - x_{k} \right)^{T}\nabla^{2}f\left( x_{k} \right)\left( \ x - x_{k} \right)$ $x_{k + 1} = \arg{\min{\left\lbrack f\left( x_{k} \right) + \left\langle \nabla f\left( x_{k} \right),\ x - x_{k} \right\rangle + \frac{1}{2}\left( \ x - x_{k} \right)^{T}\nabla^{2}f\left( x_{k} \right)\left( \ x - x_{k} \right) \right\rbrack = x_{k} - \left\lbrack \nabla^{2}f\left( x_{k} \right) \right\rbrack^{- 1}\partial f}}$

Ясно, что если сама оптимизируемая функция квадратичная, то она совпадет со своим квадратичным приближением, и мы перейдем в оптимум всего за 1 шаг.

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

Видно, что метод Ньютона использует в качестве метрики гессиан $M = \nabla^{2}f$ и эквивалентен зеркальному спуску с зеркалом $\psi(x) = \frac{1}{2}x^{T}\nabla^{2}f\left( x_{k} \right)x$ . При этом при переходе из сопряженного пространства обратно к иксам, градиент домножается как раз на обратный гессиан (как мы ранее получали связь между гессианом прямой и сопряженной функций, если гессиан невырожден)

Метрику/гессиан $M = \nabla^{2}f$ можно аппроксимировать. Такие методы называют квазиньютоновскими. Например, если М приближается матрицей ранга 2, метод называется BFGS.

Итог: преобразование Лежандра позволяет перейти к двойственной задаче, где удобно строить гарантированные оценки на функции или геометрия пространства более наглядна.

Статфиз / термодинамика

Вернемся немного к физике. Обсудим одно ложное утверждение, которое я слышал в контексте термодинамики: “преобразование Лежандра меняет интенсивные и экстенсивные величины”.

Что такое «интенсивные» и «экстенсивные»? Пусть система увеличилась в $\lambda$ раз. Тогда $g(\lambda X) = \lambda g(X)$ – экстенсивная (объем V, энтропия S), а $g(\lambda X) = g(X)$ – интенсивная (температура T, давление P). Вообще, функции со свойством $g(\lambda X) = \lambda^{k}g(X)$ называются однородными порядка . Такие функции удобны тем, что масштабируются предсказуемо.

Возьмем внутреннюю энергию . Ее переменные — экстенсивные, а производные $\frac{\partial U}{\partial S} = T,\ \ \frac{\partial U}{\partial V} = - P,\ \ \frac{\partial U}{\partial N} = \mu$ – интенсивные (их вычисление – чисто технический момент). Преобразование Лежандра дает свободную энергию Гельмгольца: $F(T,V,N) = \inf_{S}(U - TS)$ , и экстенсивная S действительно меняется на интенсивную T. Но возьмем теперь все возможные величины в расчете на 1 частицу: $U = Nu(s,v),\ \ F = Nf(T,v)$ . Важно, что s стала интенсивной. Тогда $f(T,v) = \inf_{s}\left( u(s,v) - Ts \right)$ , и Лежандр произвел обмен двух интенсивных величин. Тот факт, что часто лежандрова пара состоит из экстенсивной + интенсивной величины, это просто следствие однородности взятых функций.

Мы затронули интересный момент о сопряженности термодинамических потенциалов. Оказывается, эта ветка открывает ряд ценных связей

Энтропия

Статистический смысл энтропии появляется так. Пусть в ящике с вымышленной перегородкой находится несколько частиц. Суммарное кол-во частиц – макросостояние ящика, а распределение частиц по ячейкам – микросостояния. Всего микросостояний, реализующих данное макросостояние, — $\Omega$ штук, все они независимы и равновероятны (с вероятностью $p = \frac{1}{\Omega}$ ). Энтропия должна быть функцией от $\Omega$ – потому что просто больше не от чего оттолкнуться. При этом мы знаем, что энтропия аддитивна ( $S = S_{1} + S_{2}$ ), а количеств микросостояний мультипликативно, по аналогии с вероятностями независимых событий ( $\Omega = \Omega_{1}\Omega_{2}$ ). Единственная функция, которая переводит произведение в сумму ( $S\left( \Omega_{1}\Omega_{2} \right) = S\left( \Omega_{1} \right) + S\left( \Omega_{2} \right)$ ) – это логарифм, т.е. $S = k\ln\Omega$ . Подставив выражение в закон идеального газа, мы получим, что – константа Больцмана. Сам же ящик не обязан описывать геометрические положения частиц: он может, например, иллюстрировать распределение частиц по энергиям. В итоге $S = k\ln{\Omega(E)}$

Выразим общую энтропию через равные вероятности $p(E) = \frac{1}{\Omega(E)}$ :

$S = k\ln{\Omega(E)} = k \cdot \left( - \frac{\Omega}{\Omega} \right) \cdot \left( - \ln\Omega \right) = - k\sum_{i = 1}^{\Omega}{\frac{1}{\Omega}\ln\frac{1}{\Omega}} = - k\sum_{i = 1}^{\Omega}{p\ln p}$

Это похоже на энтропию Шеннона (об этом – позже)

Второе начало термодинамики гласит, что в изолированной системе энтропия растет, пока не достигнет максимума в термодинамическом равновесии. По виду последней функции ясно, что энтропия максимальна на равномерном распределении. Часто 2е начало искаженно пересказывают как «вселенная стремится к беспорядку», хотя точнее было бы «вселенную вероятнее застать в наиболее вероятном ее состоянии, и таким состоянием является то, что часто соответствует равномерному разбросу». На самом деле, рост энтропии – это просто следствие того факта, что равномерное распределение соответствует наибольшему числу микросостояний среди всех распределений. Эта идея лежит в основе принципа максимума энтропии

Свободная энергия

Посмотрим на вероятность $p_{i}$ микросостояния с энергией $E_{i}$ в системе из ящика и резервуара. Все состояния равновероятны, поэтому вероятность пропорциональна числу способов его реализовать: $p_{i}\ \sim\ \Omega\left( E - E_{i} \right)$ , причем энергия ящика много меньше энергии резервуара ( $E \gg E_{i}$ ). Разложив определение энтропии по Тейлору, получим $\frac{S}{k} = \ln{\Omega\left( E - E_{i} \right)} \approx \ln{\Omega(E)} - \frac{\partial\ln\Omega}{\partial E}E_{i}$ . По определению температуры, $\frac{\partial S}{\partial E} = \frac{1}{T}$ , поэтому $\frac{S}{k} \approx \ln{\Omega(E)} - \frac{E_{i}}{kT}$ . Тогда $\ln p_{i}\sim\ln{\Omega(E)} - \frac{E_{i}}{kT}$ . Взяв экспоненту и учтя константу, получим распределение Гиббса:

$p_{i} = \frac{1}{Z}e^{- E_{i}/kT} = \frac{e^{- E_{i}/kT}}{\sum_{i}^{}e^{- E_{i}/kT}}$

Это похоже на софтмакс с температурой (об этом – позже)

Подставим вероятности в формулу энтропии (вспомнив определение внутренней энергии):

$S = - k\sum_{}^{}{p_{i}\ln p_{i}} = - k\sum_{}^{}{p_{i}\left( - \ln Z - \frac{E_{i}}{kT} \right)} = k\ln Z + \frac{U}{T}\ \ \Rightarrow \ \ U - TS = - kT\ln Z$

А оказалось, что – это определение свободной энергии, поэтому:

$F = - \frac{\ln Z}{\beta},\ \ где\ Z = \sum_{i}^{}e^{- \beta E_{i}}\ и\ \ \beta = \frac{1}{kT}\ \ \ \Rightarrow \ \ \ F = - kT\ln\left\lbrack \sum_{i}^{}e^{- E_{i}/kT} \right\rbrack$

Это похоже на log-sum-exp (об этом – позже)

Связь энтропии и свободной энергии

Проинтегрируем по уровням энергии через плотность состояний:

$Z(\beta) = \int_{}^{}{\Omega(E)e^{- \beta E}dE} = \int_{}^{}e^{S(E)/k\ - \beta E}dE,\ \ где\ S,E\ \sim\ N$

При огромном числе частиц ( $N \rightarrow \infty$ ) работает ЦПТ. Поэтому подынтегральная функция концентрируется вблизи максимума, вырождается в дельта-функцию, и вклад в интеграл дает лишь одно значение подынтегральной функции – в ее максимуме (это на пальцах метод перевала):

$Z(\beta) = \exp\left\lbrack \sup_{E}\left( \frac{S(E)}{k}\ - \beta E \right) \right\rbrack$

Продифференцировав $\frac{\partial}{\partial E}\left( \frac{S(E)}{k}\ - \beta E \right) = 0$ , получим термодинамическое определение температуры: $\frac{\partial S}{\partial E} = \frac{1}{T}$ . Как уже было в разделе про Лагрàнжа и Гамильтона с определением импульса, тут определение температуры тоже возникло из поиска экстремума. Мы нашли максимум, потому что – вогнутая функция ( выпуклая), реализовав принцип максимума энтропии. Это логично: наличие установившейся температуры – понятный признак термодинамического равновесия, где энтропия максимальна.

Вернемся к определению свободной энергии:

$- \beta F(\beta) = \ln Z = \sup_{E}\left( \frac{S(E)}{k}\ - \beta E \right) = \sup_{E}\left( - \beta \cdot E\ - \left( - \frac{S(E)}{k} \right) \right) = \sup_{E}\left( p \cdot E - f(E) \right)$

Теперь, когда функция $f(E) = - \frac{S(E)}{k}$ стала выпуклая, проявилось преобразование Лежандра (через sup). Сама лежандрова пара – это функция $\ln Z$ по переменной $- \beta$ и функция $- \frac{S}{k}$ по переменной .

Такая пара может показаться громоздкой, но на самом деле она скрывает красивую связь: мы показывали, что $\ln Z$ – это log-sum-exp, а $- \frac{S(E)}{k}$ – это $- \sum_{}^{}{p\ln p}$ . Итог: log-sum-exp и негэнтропия лежандрово сопряжены

Итог: преобразование Лежандра в статфизе порождает много полезных для статистики и ML идей

Статистика

Если задача тервера понять по параметрам модели, какие данные она генерирует, то задача статистики – по данным получить оценку параметров. Заметная часть статистики основана на преобразовании Лежандра

Энтропия, softmax, logsumexp

Вспомним раздел про статфиз (энтропия, Гиббс, свободная энергия). Тем, кто сталкивался со статистикой/ML, формулы в разделе про статфиз уже заспойлерили 3 вещи:

энтропия имеет вид $- \sum_{}^{}{p\ln p}$ , как энтропия Шеннона
распределение гиббса похоже на softmax с температурой
свободная энергия имеет вид log-sum-exp

Возьмем негэнтропию Шеннона на вероятностном симплексе. Рассмотрим дискретные вероятности (в непрерывном случае действия аналогичны, но с интегралами вместо сумм). Плюс договоримся использовать натуральный логарифм вместо двоичного:

$- H(p) = \sum_{i}^{}{p_{i}\ln p_{i}},\ \ p \in \Delta = \left\{ p:\ \ \ p_{i} \geq 0\ \ ,\ \sum_{i}^{}p_{i} = 1 \right\}$

Посчитаем ее градиент:

$\nabla\left( - H(p) \right) = \begin{bmatrix} \ln p_{1} + 1 & \ldots & \ln p_{n} + 1 \end{bmatrix} = \theta$

Градиент негэнтропии не имеет какого-то устоявшегося названия, но с точностью до константы – это логиты / логарифмические шансы: $\theta_{i} = \ln\frac{p_{i}}{p_{base}} = \ln p_{i} - \ln p_{base}$ , причем точное равенство при $p_{base} = \frac{1}{e}$ .

Посчитаем гессиан и проверим, что функция выпуклая:

$\nabla^{2}\left( - H(p) \right) = \begin{pmatrix} \frac{1}{p_{1}} & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & \frac{1}{p_{n}} \end{pmatrix} \succ 0$

Гессиан негэнтропии работает с величинами, обратными вероятностям. По аналогии с лежандровой парой лагранжиан-гамильтониан в физике, и как мы выводили по теореме об обратной функции, можно предположить, что гессиан сопряженной функции как-то работает с прямыми вероятностями. Будет ли он формально обратным – не факт, надо проверять по теореме об обратной функции

Посчитаем преобразование Лежандра от негэнтропии:

$\psi^{*}(\theta) = \sup_{p \in \Delta}\left( \left\langle \theta,p \right\rangle - \sum_{i}^{}{p_{i}\ln p_{i}} \right)$

Это задача условной оптимизации на симплексе. Решать умеем с помощью функции Лагранжа из раздела про методы оптимизации:

$L = \sum_{i}^{}{p_{i}\theta_{i}} - \sum_{i}^{}{p_{i}\ln p_{i}} + \lambda\left( 1 - \sum_{i}^{}p_{i} \right),\ \ \frac{\partial L}{\partial p_{i}} = \theta_{i} - \ln p_{i} - 1 - \lambda = 0\ \ \Rightarrow \ \ p_{i} = e^{\theta_{i} - 1 - \lambda}$

Отдельно отметим, что это полностью соответствует определению преобразования Лежандра, т.к. множитель при $\lambda$ зануляется на всех из рассматриваемого пространства. Саму $\lambda$ подбираем из нормировки $\sum_{i}^{}p_{i} = 1$ :

$p_{i} = \frac{e^{\theta_{i}}}{\sum_{i}^{}e^{\theta_{i}}}\ \ \ \Rightarrow \ \ \ p = softmax\ \theta$

Подставим в $\psi^{*}$ :

$\psi^{*}(\theta) = \sum_{i}^{}{\theta_{i}p_{i}} - \sum_{i}^{}{p_{i}\ln p_{i}} = \sum_{i}^{}{p_{i}\theta_{i}} - \sum_{i}^{}{p_{i}\left( \theta_{i} - \ln{\sum_{i}^{}e^{\theta_{i}}} \right) = \sum_{i}^{}{p_{i}\ln{\sum_{i}^{}e^{\theta_{i}}}}} = \ln{\sum_{i}^{}e^{\theta_{i}}}$ $\left( - H(p) \right)^{*} = logsumexp(\theta)$

Продифференцируем, чтобы проверить наши догадки:

$\nabla logsumexp(\theta) = \begin{bmatrix} \frac{e^{\theta_{1}}}{\sum_{i}^{}e^{\theta_{i}}} & \ldots & \frac{e^{\theta_{n}}}{\sum_{i}^{}e^{\theta_{i}}} \end{bmatrix} = softmax\ \theta = p$

Градиент logsumexp уже имеет устоявшееся название, собственно, softmax. Посчитаем гессиан:

$\nabla^{2}logsumexp(\theta) = \backslash громоздкая\ выкладка\backslash\text{=}diag(p) - pp^{T} = \begin{pmatrix} p_{1} - p_{1}^{2} & \cdots & - p_{1}p_{n} \\ \vdots & \ddots & \vdots \\ - p_{n}p_{1} & \cdots & p_{n} - p_{n}^{2} \end{pmatrix} \succcurlyeq 0$

Гессиан вырожден, обратимости нет (симплекс и $\mathbb{R}^{n}$ имеют разную размерность: условие нормировки «съедает» одну степень свободы). Но мы почти угадали, что возникнут вероятности (правда, с комбинированной размерностью).

Наконец, все наши действия можно обобщить на непрерывный случай:

$\psi(p) = - H\lbrack p\rbrack = \int_{}^{}{p(x)\ln{p(x)}dx},\ \ \nabla\left( - H\lbrack p\rbrack \right) = \ln{p(x)} + 1,\ \ \nabla^{2}\left( - H\lbrack p\rbrack \right) = \frac{\delta(x - y)}{p(x)}$

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

$\psi^{*}(\theta) = \ln{\int_{}^{}{e^{\theta(x)}dx}},\ \ \nabla\psi^{*}(\theta) = p(x),\ \ \nabla^{2}\psi^{*}(\theta) = p(x)\delta(x - y) - p(x)p(y)$

Получили как градиент — саму плотность, а как гессиан – оператор, похожий по структуре на матрицу выше

Принцип максимума энтропии и ОМП

В блоке про энтропию мы провели такое преобразование Лежандра:

$\psi^{*}(\theta) = \sup_{p \in \Delta}\left( \left\langle \theta,p \right\rangle - \sum_{i}^{}{p_{i}\ln p_{i}} \right)$

И оказалось, что связанная преобразованием Лежандра пара переменных – вероятности и параметры $\theta$ .

Что означает сам супремум: мы максимизировали энтропию при фиксированном матоже параметра $\mathbb{E}_{p}\theta = \left\langle \theta,p \right\rangle = \sum_{i}^{}{p_{i}\theta_{i}}$ . В этом ограничении, softmax – распределение максимальной энтропии. Оказывается, к максимизации энтропии можно подойти с другой стороны, и начнем с пары прикладных задач.

Пример 1: монетка $X_{i}\ \sim\ Bernoulli(\mu)$ . Подбросим ее раз, посчитаем долю появлений решки $m = \frac{1}{n}\sum_{}^{}x_{i}$ . Как оценить по нашим данным, какое реальное распределение скорее всего могло эти данные породить?

Формально, мы ищем $\widehat{\mu} = \arg{\max_{\mu}{p\left( x \middle| \mu \right)}} = \arg{\max_{\mu}{\ln\left( p\left( x \middle| \mu \right) \right)}}$ , т.к. логарифм строго возрастает и не сдвигает аргмаксимум. Вероятность $p\left( x \middle| \mu \right)$ называют правдоподобием, функцию $l(\mu) = \ln\left( p\left( x \middle| \mu \right) \right)$ – логарифмом правдоподобия, ее градиент – скором, а саму $\widehat{\mu}$ – оценкой максимального правдоподобия (ОМП). Посчитаем:

$l(\mu) = \ln{p\left( x \middle| \mu \right)} = \ln{\prod_{i = 1}^{n}{p\left( x_{i} \middle| \mu \right)}} = \sum_{i = 1}^{n}{\ln{p\left( x_{i} \middle| \mu \right)}} = \sum_{i = 1}^{n}{\ln\left( \mu^{x_{i}}(1 - \mu)^{1 - x_{i}} \right)} = \ln\mu \cdot \sum_{}^{}x_{i} + \ln(1 - \mu) \cdot \left( n - \sum_{}^{}x_{i} \right) = mn\ln\mu + (n - mn)\ln(1 - \mu)$ $0 = \frac{d}{d\mu}l(\mu) = \frac{m}{\mu}n - \frac{1 - m}{1 - \mu}n\ \ \ \Rightarrow \ \ \frac{\widehat{\mu}}{1 - \widehat{\mu}} = \frac{m}{1 - m}\ \ \ \Rightarrow \ \ \ \ \widehat{\mu} = m$

Наиболее правдоподобна ситуация – когда монетка выпадает решкой с вероятностью, равной нашей экспериментально посчитанной доле решек, что логично, если о монетке не допускать ничего дополнительно.

Заметим, что в выводе у нас возникала конструкция $\frac{\widehat{\mu}}{1 - \widehat{\mu}}$ . Плюс мы в начале раздела заметили, что параметр $\theta$ , участвующий в сопряженной паре, был логитом. Хочется подогнать задачу под уже знакомый нам вид: возьмем параметр $\theta = \ln\frac{\mu}{1 - \mu}$ . Оптимальный $\widehat{\theta} = \ln\frac{m}{1 - m}$

Как бы мы решали эту задачу раньше, методом максимальной энтропии? Поищем такое распределение , которое дает максимальную энтропию при условии, что его матожидание равно наблюдаемому среднему: $\mathbb{E}_{p}X = \sum p(x)x = p(1) = m$ . Запишем задачу условной оптимизации, она же – преобразование Лежандра в пространстве вероятностей (два множителя Лагранжа — $\theta$ и $\lambda$ ):

$L = - \sum_{}^{}{p\ln p} + \theta\left( \sum_{}^{}{px} - m \right) + \lambda\left( 1 - \sum_{}^{}p \right),\ \ 0 = \frac{\partial L}{\partial p} = - \ln p - 1 - \theta x - \lambda\ \ \ \Rightarrow \ \ p(x) \propto e^{\theta x}$ $p(1) = m = \frac{e^{\widehat{\theta} \cdot 1}}{e^{\widehat{\theta} \cdot 1} + e^{\widehat{\theta} \cdot 0}} = \frac{1}{1 + e^{- \widehat{\theta}}} = \sigma(\widehat{\theta})$

Эта функция, обратная логиту, называется сигмоидой. Для категориального мы получали похожий результат — softmax. Из этого равенства можно выразить оптимальный $\widehat{\theta} = \ln\frac{m}{1 - m}$ – мы получили ровно такой же, как для ОМП. К тому же, укрепилась наша интуиция о записи параметра в виде логита.

Пример 2: нормальное распределение $X_{i}\ \sim\ N\left( \mu,\sigma^{2} \right),\ \ m_{1} = \frac{1}{n}\sum_{}^{}x_{i},\ \ m_{2}^{2} = \frac{1}{n}\sum_{}^{}x_{i}^{2}$ . Теперь у нас 2 параметра + 2 статистики, которые мы считаем в эксперименте. Обозначим вектор статистик как $T(x) = \begin{pmatrix} x & x^{2} \end{pmatrix}^{T}$

ОМП:

$l\left( \mu,\sigma^{2} \right) = - \frac{n}{2}\ln\left( 2\pi\sigma^{2} \right) - \frac{1}{2\sigma^{2}}\sum_{}^{}\left( x_{i} - \mu \right)^{2},\ \ \partial_{\mu} = \partial_{\sigma^{2}} = 0\ \ \Rightarrow \ \ \begin{pmatrix} \widehat{\mu} \\ {\widehat{\sigma}}^{2} \end{pmatrix} = \begin{pmatrix} m_{1} \\ m_{2}^{2} - m_{1}^{2} \end{pmatrix}$

Принцип максимума энтропии (при условиях $\mathbb{E}_{p}X = m_{1},\ \ \mathbb{E}_{p}X^{2} = m_{2}^{2}$ , со множителями $\theta_{1},\theta_{2},\lambda$ ):

$L = - \int_{}^{}{p\ln p} + \theta_{1}\left( m_{1} - \int_{}^{}{pxdx} \right) + \theta_{2}\left( m_{2}^{2} - \int_{}^{}{px^{2}dx} \right) + \lambda\left( 1 - \int_{}^{}{pdx} \right)$ $0 = \frac{\delta L}{\delta p} = - \ln p - 1 - \theta_{1}x - \theta_{2}x^{2} - \lambda\ \ \ \Rightarrow \ \ p(x) \propto e^{\theta_{1}x + \theta_{2}x^{2}}$

Значит, — гауссиана, а оптимальные $\begin{pmatrix} {\widehat{\theta}}_{1} & {\widehat{\theta}}_{2} \end{pmatrix} = \begin{pmatrix} \frac{\widehat{\mu}}{{\widehat{\sigma}}^{2}} & - \frac{1}{2{\widehat{\sigma}}^{2}} \end{pmatrix}$

На этих примерах оказалось, что метод максимальной энтропии и ОМП – прямая и Лежандрово двойственная задачи. Всегда ли это так?

Экспоненциальный класс, информация Фишера и ковариация

Что общего можно заметить в обоих примерах:

мы не работали с отдельными элементами выборки. Проводя эти эксперименты на практике, мы могли вообще не записывать отдельные исходы. Нам было достаточно лишь хранить 1-2 статистики и обновлять их итеративно – и как будто в этих статистиках было достаточно информации, чтобы описать всё распределение.
плотности выражались через экспоненту, если записывать их через правильные (с точки зрения преобразования Лежандра) параметры. Именно эти «естественные» параметры сопряжены с наблюдаемыми средними значениями тех «достаточных» статистик

Оказывается, существует удобный класс распределений, обобщающий наши наблюдения – он называется экспоненциальным. Все плотности в нем записываются как:

$p\left( x \middle| \theta \right) = h(x)\exp{\left( \theta T(x) - f(\theta) \right),\ \ где\ f(\theta) = \ln{\int_{}^{}{h(x)\exp{\left( \theta T(x) \right)dx}}}}$

Названия компонент: $\theta$ – натуральный параметр, – достаточная статистика, – базовая мера. Показатели экспонент линейны по $\theta$ , а нормировка $f(\theta)$ имеет знакомый нам вид logsumexp и выпукла (вне экспоненциального класса, кстати, выпуклость f не гарантируется):

$\nabla f = \mathbb{E}_{\theta}T(X) = \mu,\ \ \nabla^{2}f = Var\ T(X) \geq 0$

При расчете правдоподобия логарифм как раз полезен тем, что «съедает» экспоненту. Для 1 измерения:

$\frac{1}{n}l(\theta) = const + \theta\widehat{\mu} - f(\theta),\ \ где\ \widehat{\mu} = \frac{1}{n}\sum_{}^{}{T\left( x_{i} \right)}$

Сразу видно, что $\mathbb{E}\nabla l(\theta)\mathbb{= E}T - \mu = 0.$ А max по $\theta$ – буквально определение преобразования Лежандра:

$f^{*}\left( \widehat{\mu} \right) = \max_{\theta}\left( \theta\widehat{\mu} - f(\theta) \right),\ \ \theta = \nabla f^{*}(\mu) = (\nabla f)^{- 1}(\mu),\ \ \nabla^{2}f^{*} = \left( Var\ T(X)\ \right)^{- 1}$

В методе максимума энтропии мы минимизируем негэнтропию при условии $\mathbb{E}_{p}T = \widehat{\mu}$ . Оптимум дает $p(x) \propto he^{\theta T}$ – в экспоненциальном классе, а двойственная функция и преобразование Лежандра эквивалентны:

$\min_{p}{KL\left( p||h \right)}при\ \mathbb{E}_{p}T = \widehat{\mu}$

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

гессиан негэнтропии оказался равен его матрице Фишера (по определению, $I(\theta) = Var\ \nabla l(\theta)\mathbb{= - E}\nabla^{2}l(\theta)$ )
гессиан logsumexp оказался равен его матрице ковариации (по определению, $Cov\ X = \mathbb{E}XX^{T}\mathbb{- E}X\mathbb{E}X^{T}$ )

Вспомним также три известных нам смысла гессиана:

аналитический: вторая производная, локальная мера выпуклости/кривизны функции
физический: аналоги массы/инерции и обратным к ним
геометрический: гессиан кодирует геометрию (метрику) пространства, на котором применяется функция

Выше мы показали, матрица Фишера (в средних параметрах) и матрица ковариации достаточной статистики (в натуральных параметрах) связаны как гессианы в преобразовании Лежандра. Из него же следует теорема Крамера-Рао. Для несмещенной оценки $\widehat{\theta}$ ( $\mathbb{E}\widehat{\theta} = \theta$ ) по КБШ:

$1 = Cov\ \left( \widehat{\theta},\nabla l(\theta) \right)^{2} \leq Var\ \widehat{\theta} \cdot Var\ \nabla l(\theta) = Var\ \widehat{\theta} \cdot I(\theta)\ \ \Rightarrow \ \ Var\ \widehat{\theta} \geq {I(\theta)}^{- 1}$

Эта теорема означает, что никакая несмещенная оценка не может быть точнее границы ${I(\theta)}^{- 1}$ . Физическая аналогия ковариаций – массы, информации Фишера – подвижности (обратные массы). Крамер-Рао означает, что чем жёстче статистика связана с параметром (больше масса), тем меньше подвижности и тем меньше система «дрожит», а значит — тем точнее можно оценить ее реальные параметры по шумным наблюдениям.

Осталось поговорить про геометрию пространств. В разделе про оптимизацию мы могли по-разному выбирать метрику M, например, в методе Ньютона мы брали гессиан функции. В статистике естественная метрика – это как раз информация Фишера, т.е. тот же гессиан (логарифма нормировки). Честный градиент с учетом геометрии этого пространства называется натуральным:

$\widetilde{\nabla}L = {I(\theta)}^{- 1}\nabla L$

Используя в задачах оптимизации натуральный градиент, мы запускаем в точности метод зеркального спуска с зеркалом – логарифмом нормировки. Этот метод активно используется в современном ML, например, популярный метод Adam аппроксимирует $I(\theta)$ диагональной матрицей (плюс пара инженерных доработок). Подробнее про применения в ML поговорим позже.

Проверка гипотез

Вернемся к задаче с монеткой $X_{i}\ \sim\ Bernoulli(\mu)$ . Пусть нам удалось сузить представление о параметре модели до двух альтернатив: $H_{0}:X\ \sim\ P_{0}$ – монетка честная ( $\mu = \frac{1}{2}$ ) vs $H_{1}:X\ \sim\ P_{1}$ – монетка фальшивая ( $\mu = q$ ). Нужно по наблюдениям определить, фальшивая ли монетка, и чем мы рискуем, принимая решение. Введем функцию $\phi(x)$ – вероятность отвергнуть $H_{0}$ по наблюдению :

Ошибки	$H_{0}$ верна	$H_{1}$ верна
Отвергаем $H_{0}$ ( $\phi = 1$ )	I рода (ложная тревога), $\alpha = \mathbb{E}_{0}\phi$	Все ок
Не отвергаем $H_{0}$ ( $\phi = 0$ )	Все ок	II рода (пропуск), $\beta = \mathbb{E}_{1}(1 - \phi)$

На практике выбирают $\alpha_{0}$ — допустимый уровень ошибки I рода (уровень значимости), гарантируют, что ложных тревог не больше этого порога ( $\alpha \leq a_{0})$ , и в этих ограничениях минимизируют $\beta$ – ошибку II рода (максимизируют мощность $1 - \beta$ ). Это условная оптимизация:

$L(\phi,\lambda) = \mathbb{E}_{1}\phi - \lambda\left( \mathbb{E}_{0}\phi - \alpha_{0} \right) = \lambda\alpha_{0} + \sum_{x}^{}{\phi(x)\left( P_{1}(x) - \lambda P_{0}(x) \right)}$

Максимум $\max_{\phi}L$ возьмем поточечно: поставим $\phi(x) = 1$ там, где скобка положительна, т.е. $\frac{P_{1}(x)}{P_{0}(x)} > \lambda$ . Заметим, что $P_{0}$ и $P_{1}$ – это правдоподобия (вероятность выборки при условиях гипотез на параметры). Мы построили критерий отношения правдоподобия как тот, что дает максимальную можность при заданном уровне значимости – это лемма Неймана-Пирсона. Перейдем к логарифму правдоподобия:

$s(x) = \ln\frac{P_{1}(x)}{P_{0}(x)},\ \ S_{n} = \frac{1}{n}\sum_{}^{}{s(x_{i})},\ \ отвергаем\ H_{0}\ при\ S_{n} > \gamma$

$S_{n}$ – достаточная статистика, нам достаточно копить только ее вместо сырых данных. По ЗБЧ:

$Под\ H_{0}:\ \mathbb{E}_{0}s = \int_{}^{}{P_{0}(x)\ln\frac{P_{1}(x)}{P_{0}(x)}dx} = - KL(P_{0}|\left| P_{1} \right) \leq 0$ $Под\ H_{1}:\ \mathbb{E}_{1}s = \int_{}^{}{P_{1}(x)\ln\frac{P_{1}(x)}{P_{0}(x)}dx} = KL(P_{1}|\left| P_{0} \right) \geq 0$

Оценим ошибку I рода $\alpha_{n} = P_{0}\left( S_{n} \geq \gamma \right)$ . Для любого $\lambda \geq 0$ индикатор мажорируется экспонентой (неравенство Чернова):

$\alpha_{n} \leq e^{- \lambda n\gamma}\mathbb{E}_{0}\left( e^{\lambda\sum_{}^{}{s\left( x_{i} \right)}} \right) = e^{- \lambda n\gamma}\left( \mathbb{E}_{0}e^{\lambda s} \right)^{n} = e^{- n\left( \lambda\gamma - \psi(\lambda) \right)},\ \ \psi(\lambda) = \ln{\mathbb{E}_{0}\left( e^{\lambda s} \right)} = \ln{\sum_{x}^{}{P_{0}(x)^{1 - \lambda}P_{1}(x)^{\lambda}}}$

Мы получили целое семейство оценок сверху, из которых берем лучшую. Преобразование Лежандра $\psi^{*}$ от лог-статсуммы / производящей функции моментов $\psi$ называется функцией уклонения:

$\alpha_{n} \leq e^{- n\psi^{*}(\gamma)},\ \ \psi^{*}(\gamma) = \sup_{\lambda}\left( \lambda\gamma - \psi(\lambda) \right) = KL\left( P_{\lambda^{*}(\gamma)}||P_{0} \right)$

Заметим, что $\psi$ – это тот же logsumexp нормировщик f из экспоненциального класса:

$P_{\lambda}(x) = \frac{P_{0}(x)^{1 - \lambda}P_{1}(x)^{\lambda}}{e^{\psi(\lambda)}} = P_{0}(x)e^{\lambda s(x) - \psi(\lambda)}$

Мы уже работали с подобной конструкцией в разделе Статфиз: – энергия, $\lambda$ – обратная температура, $\psi$ – свободная энергия, $\psi^{*}$ — энтропия, $P_{\lambda}$ – распределение Гиббса (софтмакс с температурой $\frac{1}{\lambda}$ ). Наклон $\lambda^{*}$ ищем так, чтобы редкое значение $\gamma$ стало типичным средним, применяем метод перевала.

Вероятность получить данное или более экстремальное значение статистики при верной $H_{0}$ называют p-value. Его считают по данным, и если $pvalue \leq \alpha_{0}$ , то отклонение статистики считают неправдоподобно большим и $H_{0}$ отвергают. Оказалось, что:

$pvalue = P_{0}\left( S_{n} \geq \gamma \right)\ \propto \ e^{- n\psi^{*}(\gamma)},\ \ KL\left( {\widehat{P}}_{n}||P_{0} \right) \approx - \frac{1}{n}\ln(pvalue) \rightarrow \psi^{*}(\gamma)$

Минус логарифм p-value – это аналог энтропии и мера нашего удивления данным (насколько сильно отличаются эмпирическое и предполагаемое распределения). Само p-value – свойство хвоста распределения: когда у статистики есть производящая функция, её хвост управляется лежандровым образом этой функции

Итог: преобразование Лежандра пронизывает статистику, предоставляя много ценных идей и методов

ML

Байесовский вывод

Ключевая работа Байеса была опубликована посмертно – именно она заложила основы байесовской статистики и значительной части современного ML. Байес исследовал такой мысленный эксперимент:

экспериментатор с завязанными глазами кидает белый шар на бильярдный стол и хочет узнать, куда он попал ( $\theta$ – относительная позиция шара на столе, 0 – левый край, 1 — правый)
для этого экспериментатор кидает на стол другой шар, а ассистент подсказывает, куда тот приземлился относительно белого. Эта информация помогает скорректировать представление о положении белого шара
процедура повторяется n раз, пока экспериментатор не достигнет достаточной уверенности в ответе

$P\left( белый\ на\ \theta\ |\ справа\ k\ цветных \right) = \frac{P\left( справа\ k\ цветных\ |\ белый\ на\ \theta \right) \cdot P(белый\ на\ \theta)}{P(справа\ k\ цветных)}$

$P\left( белый\ на\ \theta\ |\ справа\ k\ цветных \right)$ – апостериор (представление о параметре, уточненное по данным)

$P(белый\ на\ \theta)$ – априор (текущая гипотеза), в начале предполагаем, что шар приземляется равномерно, т.е. имеет распределения $Beta(1,1) \equiv 1$

$P\left( справа\ k\ цветных\ |\ белый\ на\ \theta \right)$ – знакомое нам правдоподобие, у нас оно равно $C_{n}^{k}\theta^{k}(1 - \theta)^{n - k}$ – это биномиальное распределение (сумма n независимых испытаний Бернулли)

$P(справа\ k\ цветных)$ – называется evidence. По сути, нам интересен только числитель, а знаменатель — это нормировочная константа, которую на практике часто оценивают или подгоняют

Теорема Байеса объединяет определение условной вероятности и формулу полной вероятности:

$P\left( \theta\ |\ k \right) = \frac{P(\theta,k)}{P(k)} = \frac{P\left( k\ |\ \theta \right) \cdot P(\theta)}{\sum_{\theta}^{}{P\left( k\ |\ \theta \right) \cdot P(\theta)}}$

В нашем случае:

$P\left( \theta\ |\ k \right) = \frac{C_{n}^{k}\theta^{k}(1 - \theta)^{n - k} \cdot 1}{(n + 1)^{- 1}}\ \propto \ \theta^{k}(1 - \theta)^{n - k},\ \ это\ Beta(k + 1,n - k + 1)$

Что здесь можно заметить:

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

Если априор не меняет вид апостериора, то априор называется сопряженным к нему. Слово «сопряженный» намекает на преобразование Лежандра. Пусть дано распределение из экспоненциального класса:

$p\left( x \middle| \theta \right) = h(x)\exp{\left( \theta T(x) - f(\theta) \right),\ \ где\ f(\theta) = \ln{\int_{}^{}{h(x)\exp{\left( \theta T(x) \right)dx}}}}$

Правдоподобие для все выборки размера n:

$p\left( X \middle| \theta \right) = \prod_{}^{}{h\left( x_{i} \right)} \cdot \exp\left( \theta\sum_{}^{}{T\left( x_{i} \right)} - nf(\theta) \right)$

Чтобы получить сопряженный априор, мы требуем, чтобы апостериорное распределение сохранило форму. Единственный способ это сделать — задать априор в виде:

$p\left( \theta \middle| \alpha,\nu \right) = \exp\left( \theta\alpha - \nu f(\theta) - g(\alpha,\nu) \right)$

Где $\alpha$ – сумма достаточных статистик из полученных на данный момент (или гипотетических) наблюдений, $\nu$ – их количество (в гипотетическом случае – мера уверенности в априоре), $g(\alpha,\nu)$ – нормировочная константа. Мы хотим так подобрать эту константу, чтобы априор был корректным, т.е. интегрировался в единицу:

$g(\alpha,\nu) = \ln{\int_{}^{}{\exp\left( \theta\alpha - \nu f(\theta) \right)d\theta}}$

Решается уже знакомым нам методом перевала:

$g(\alpha,\nu) \approx \sup_{\theta}\left( \theta\alpha - \nu f(\theta) \right) = \nu \cdot \sup_{\theta}\left( \theta \cdot \frac{\alpha}{\nu} - f(\theta) \right)$

Оказалось, чтобы нормировать априор, требуется взять как преобразование Лежандра от :

$g(\alpha,\nu) = \nu \cdot f^{*}\left( \frac{\alpha}{\nu} \right)$

Регуляризация в линейной регрессии

Рассмотрим задачу: даны числовые метки и признаки . Предположим, что метки порождены линейной моделью с нормальным шумом как $y_{i}\ \sim\ N\left( \theta^{T}x_{i},\sigma^{2} \right)$ , и попробуем подобрать оптимальные веса $\widehat{\theta}$ – параметры модели. Сразу отметим, что часто модель записывают как $y_{i} = \theta^{T}x_{i} + b_{i}$ , но это эквивалентно добавлению признака из всех единиц и включению $b_{i}$ в вектор весов $\theta$ – так что будем писать в простом виде

Что означает «оптимальные»? Например, полученные как ОМП. Альтернативно — можно минимизировать среднюю ошибку модели, например, MSE — среднеквадратичную:

$MSE(\theta) = \frac{1}{2}\left\| Y - X\theta \right\|_{2}^{2} = \frac{1}{2}(Y - X\theta)^{T}(Y - X\theta) = \frac{1}{2}\sum_{i = 1}^{n}\left( y_{i} - \theta^{T}x_{i} \right)^{2}$

Итог: оптимум можно найти аналитически, градиентным спуском и как ОМП. Рассмотрим все варианты

1) Аналитический подход (метод наименьших квадратов, МНК):

$\frac{d}{d\theta}MSE = X^{T}X\theta - X^{T}Y = 0\ \ \Rightarrow \ \ \widehat{\theta} = \left( X^{T}X \right)^{- 1}X^{T}Y$

Возникает проблема при обращении матрицы $X^{T}X$ . Если она вырождена, обратить нельзя вовсе. Если она близка к вырожденной – обратить можно, но сами веса $\theta$ становятся огромными, и модель оказывается чересчур чувствительна ко входным данным – это называется переобучением. Инженерное решение: давайте добавим к MSE-лоссу штраф за размер весов – например, за сумму их квадратов, т.е. $L_{2}$ -норму:

$MSE(\theta) = \frac{1}{2}\left\| Y - X\theta \right\|_{2}^{2} + \frac{\lambda}{2}\left\| \theta \right\|_{2}^{2},\ \ \frac{d}{d\theta} = 0\ \ \Rightarrow \ \ \widehat{\theta} = \left( X^{T}X + \lambda I \right)^{- 1}X^{T}Y$

Теперь обращаемая матрица отделилась от вырожденной, сделав модель стабильнее.

2) Подход через максимум правдоподобия (ОМП):

Вспомним, что данные генерируются как $y_{i}\ \sim\ N\left( \theta^{T}x_{i},\sigma^{2} \right)$ . Для выборки из n независимых наблюдений:

$\ln{p\left( Y\ |\ X,\theta \right)} = \ln{\prod_{i = 1}^{n}{\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\left( - \frac{\left( y_{i} - \theta^{T}x_{i} \right)^{2}}{2\sigma^{2}} \right)}} = const - \frac{1}{2\sigma^{2}}\sum_{i = 1}^{n}\left( y_{i} - \theta^{T}x_{i} \right)^{2}$

Мы получили буквально выражение для $MSE(\theta)$ – только с минусом (ищем минимум ошибки vs максимум правдоподобия). Продолжим решение:

$\sum_{i = 1}^{n}\left( y_{i} - \theta^{T}x_{i} \right)^{2} = Y^{T}Y - 2\theta^{T}X^{T}Y + \theta^{T}X^{T}X\theta$ $\nabla_{\theta}\ln{p\left( Y\ |\ X,\theta \right)} = - \frac{1}{\sigma^{2}}\left( X^{T}X\theta - X^{T}Y \right) = 0\ \ \Rightarrow \ \ \widehat{\theta} = \left( X^{T}X \right)^{- 1}X^{T}Y$

Мы нашли в точности оценку по МНК – с той же проблемой возможной вырожденности $X^{T}X$ .

Пока что мы не учли всю силу теоремы Байеса. Что мы можем сказать о весах модели априорно, пока мы их еще не оценили? В целом, мы можем оттолкнуться от любого априора, и достаточное количество наблюдений его в конце концов исправят. Но раз мы можем выбирать, то давайте выберем распределение, которое сильно упростит расчеты: сопряженный априор. Оказывается, сопряженное к нормальному – тоже нормальное. Действительно: мы в самом начале получали, что $\psi(x) = \frac{1}{2}x^{2}$ – неподвижная точка преобразования Лежандра, а логарифм плотности нормального распределения имеет ровно такой же вид. Пусть:

$\theta\ \sim\ N\left( \theta_{0},\Sigma_{0} \right),\ \ p(\theta) = \frac{1}{(2\pi)^{d/2}\left| \Sigma_{0} \right|^{1/2}\ }\exp\left( - \frac{1}{2}\left( \theta - \theta_{0} \right)^{T}\Sigma_{0}^{- 1}\left( \theta - \theta_{0} \right) \right)$

Байесовский апдейт:

$p\left( \theta\ |\ X,Y \right) = \frac{p\left( Y\ |\ X,\theta \right) \cdot p(\theta)}{p\left( Y\ |\ X \right)}\ \propto \ \exp\left( - \frac{1}{2\sigma^{2}}(Y - X\theta)^{T}(Y - X\theta) - \frac{1}{2}\left( \theta - \theta_{0} \right)^{T}\Sigma_{0}^{- 1}\left( \theta - \theta_{0} \right) \right)$

Апостериорные ковариация и среднее:

$\Sigma_{n} = \left( \Sigma_{0}^{- 1} + \frac{1}{\sigma^{2}}X^{T}X \right)^{- 1},\ \ \theta_{n} = \Sigma_{n}\left( \Sigma_{0}^{- 1}\theta_{0} + \frac{1}{\sigma^{2}}X^{T}Y \right)$

Нормальное распределение принадлежит экспоненциальному классу, и его нормировочная константа равна $f(\eta) = \frac{1}{2}\eta^{T}\Sigma\eta$ , где $\eta = \Sigma^{- 1}\theta$ – натуральный параметр. Преобразование Лежандра: $f^{*}(\theta) = \frac{1}{2}\theta^{T}\Sigma^{- 1}\theta$ . Оказалось, при байесовском апдейте, апостериорная точность (обратная ковариация) суммируется: $\Sigma_{n}^{- 1} = \Sigma_{0}^{- 1} + \frac{1}{\sigma^{2}}X^{T}X$ – т.е. буквально происходит сложение информации Фишера в сопряженном пространстве.

Посчитав $\theta_{n}$ (т.е. проведя апдейт по всему датасету), получим уточненную оценку оптимального параметра $\widehat{\theta}$ :

$\widehat{\theta} = \left( \Sigma_{0}^{- 1} + \frac{1}{\sigma^{2}}X^{T}X \right)^{- 1}\left( \Sigma_{0}^{- 1}\theta_{0} + \frac{1}{\sigma^{2}}X^{T}Y \right)$

Мы выбирали априор сами, пусть он будет удобным, например $\Sigma_{0} = \sigma_{0}^{2}I$ и $\theta_{0} = 0$ . Тогда:

$\widehat{\theta} = \left( X^{T}X + \frac{\sigma^{2}}{\sigma_{0}^{2}}I \right)^{- 1}X^{T}Y$

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

$\widehat{\theta} = \left( X^{T}X + \lambda I \right)^{- 1}X^{T}Y,\ \ \lambda = \frac{\sigma^{2}}{\sigma_{0}^{2}}$

Наконец, отметим, что связь MSE-лосса и правдоподобия чуть глубже. Пусть вектор ошибки $\varepsilon = Y - X\theta$ , тогда $MSE(\theta) = \frac{1}{2}\left\| \varepsilon \right\|_{2}^{2}$ . Перейдем к двойственной задаче через преобразование Лежандра по $\varepsilon$ . Снова вспомним, что $\psi(x) = \frac{1}{2}x^{2}$ – неподвижная точка преобразования Лежандра, поэтому $MSE^{*}(\lambda) = \frac{1}{2}\left\| \lambda \right\|_{2}^{2},$ где $\lambda$ – двойственная переменная. Решение двойственной задачи имеет тот же вид. Более того, теперь решение без регуляризации $\widehat{\theta} = \left( X^{T}X \right)^{- 1}X^{T}Y$ явно выходит как частный случай, при ограничении $X^{T}\lambda = 0$

3) Итеративный подход:

$\theta_{t + 1} = \theta_{t} - \eta\nabla MSE\left( \theta_{t} \right),\ \ где\ \nabla MSE\left( \theta_{t} \right) = X^{T}X\theta_{t} - X^{T}Y$

Градиентный спуск останавливается, когда выходит на плато, т.е. $0 = \nabla MSE(\theta) = X^{T}X\theta - X^{T}Y$ , откуда снова получается решение $\widehat{\theta} = \left( X^{T}X \right)^{- 1}X^{T}Y$ .

Возникает новая проблема — теперь с учетом геометрии пространства весов: признаки могут иметь разный масштаб, поэтому градиентный спуск очень долго сходится. Когда-то для учета корректной метрики мы придумали зеркальный спуск – применим его здесь. Мы уже знаем, что MSE – это отрицательный логарифм правдоподобия, поэтому посчитаем натуральный градиент:

$I(\theta)\mathbb{= - E}\nabla^{2}l(\theta) = \nabla\frac{1}{\sigma^{2}}\left( X^{T}X\theta - X^{T}Y \right) = \frac{1}{\sigma^{2}}X^{T}X\ \ \ \ \Rightarrow \ \ \ \ \widetilde{\nabla}\ = \left( I(\theta) \right)^{- 1}\nabla\ = \sigma^{2}\left( X^{T}X \right)^{- 1}$

Возьмем удобный размер шага $\eta = \frac{1}{\sigma^{2}}$ :

$\theta_{t + 1} = \theta_{t} - \eta\sigma^{2}\left( X^{T}X \right)^{- 1}\left( X^{T}X\theta_{t} - X^{T}Y \right) = \theta_{t} - \theta_{t} + \left( X^{T}X \right)^{- 1}X^{T}Y = \left( X^{T}X \right)^{- 1}X^{T}Y$

То есть просто учтя геометрию пространства, мы всего за 1 шаг перешли в оптимальное решение. Почему так произошло? Вспомним метод Ньютона: берем квадратичное приближение функции, переходим в его минимум. Но MSE сама по себе квадратична, т.е. равна своему квадратичному приближению, и перейдя в минимум приближения, мы перешли в минимум всей функции. Градиентный спуск по MSE с натуральным градиентом эквивалентен методу Ньютона.

Остается два финальных вопроса: почему мы выбрали именно $L_{2}$ -норму весов и квадратичный MSE-лосс? На первый мы уже ответили – это эквивалентно предположению о том, что априорно веса распределены нормально. Мы могли использовать $L_{1}$ -норму весов (сумму модулей) – это эквивалентно предположению о том, что априорно веса распределены по Лапласу. Оказывается, ответ на второй вопрос – аналогичный: MSE-лосс кодирует априорное предположение, что данные генерируются с нормальным шумом, а MAE-лосс (сумма модулей ошибок) эквивалентен предположению о лаплассовском шуме. Функция модуля хоть и не гладкая в нуле, но выпуклая, поэтому везде тоже применимо преобразование Лежандра – с парой оговорок.

Логистическая регрессия, нейрон, SVM

Рассмотрим задачу классификации: даны дискретные (или бинарные) метки и числовые признаки . Предположим, что метки порождены из категориального распределения (или Бернулли). Попробуем подобрать оптимальные веса $\widehat{\theta}$ – параметры модели. В бинарном случая типична модель логистической регрессии $y = \sigma\left( \theta^{T}x \right)$ . По сути, мы говорим, что, с помощью линейной регрессии мы моделируем логарифмические шансы (логит – натуральный параметр Бернулли), а потом применяем сигмоиду и отображаем в корректные вероятности классов. Логарифм правдоподобия категориального (или Бернулли) – это кросс-энтропия: $- \sum_{cls}^{}{y_{cls}\ln p_{cls}}$ . Связки этой задачи с физикой и преобразованием Лежандра мы детально изучили раздела Статистика.

Классическая нейросеть состоит из набора слоев, каждый слой – из нейронов, а нейрон – это функция вида $y = \sigma\left( \theta^{T}x \right)$ , где $\sigma$ – функция активации (типично это сигмоида, и тогда нейрон совпадает с логистической регрессией). Обучение происходит по шагам:

forward-pass: нейрон получает на вход точку . Линейная операция переводит его в логит $z = \theta^{T}x$ – элемент сопряженного к вероятностям пространства, активация (как $\nabla\psi(z)$ ) переводит логит в первичное пространство вероятностей/средних. Считается лосс – лежандров зазор
backward-pass: расчет градиентов весов по цепочному правилу, от лосса назад по всем слоям. Градиент лосса по логиту $\frac{\partial L}{\partial z} = \widehat{y} - y$ лежит в первичном пространстве вероятностей/средних, рассчет градиентов переводит его в сопряженное (как $\nabla\psi(y)$ ). За одну итерацию веса модели обновляются на их градиент. Получается, forward и backward pass – это полный круг преобразований Лежандра.

Вернемся к задаче бинарной классификации и подойдем к ней с другой стороны. Пусть у нас есть линейно разделимая выборка $\left( x_{i},y_{i} \right)$ , где метки $y_{i} = \pm 1$ . Ищем гиперплоскость $w^{T}x + b = 0$ (здесь уже удобно разделить коэффициенты и смещение), которая разделит классы с максимальным зазором:

Support Vector Machines (SVMs): Important Derivations | Towards Data Science

Прямая задача:

$\min{\frac{1}{2}\left\| w \right\|_{2}^{2}},\ \ при\ y_{i}\left( w^{T}x_{i} + b \right) \geq 1\ \forall i$

Чтобы учесть ограничения, введем множители Лагранжа $\alpha_{i} \geq 0$ :

$L(w,b,\alpha) = \frac{1}{2}\left\| w \right\|_{2}^{2} - \sum_{}^{}{\alpha_{i}\left( y_{i}\left( w^{T}x_{i} + b \right) - 1 \right)}$ $\frac{\partial L}{\partial w} = w - \sum_{}^{}{\alpha_{i}y_{i}x_{i}} = 0,\ \ \frac{\partial L}{\partial b} = - \sum_{}^{}{\alpha_{i}y_{i}} = 0\ \ \ \ \Rightarrow \ \ \ w = \sum_{}^{}{\alpha_{i}y_{i}x_{i}},\ \ \sum_{}^{}{\alpha_{i}y_{i}} = 0$

Подставляем обратно в лагранжиан, после всех выкладок получаем двойственную задачу:

$\max_{\alpha}L = \max_{\alpha}\left( \sum_{i}^{}\alpha_{i} - \frac{1}{2}\sum_{i,j}^{}{\alpha_{i}\alpha_{j}y_{i}y_{j}x_{i}^{T}x_{j}} \right),\ \ при\ \sum_{i}^{}{\alpha_{i}y_{i}} = 0\ и\ \alpha_{i} \geq 0$

С точки зрения Лежандра, произошло преобразование знакомой нам стационарной точки:

$f(w) = \frac{1}{2}\left\| w \right\|_{2}^{2},\ \ f^{*}(\lambda) = \sup_{w}\left( \lambda^{T}w - \frac{1}{2}\left\| w \right\|_{2}^{2} \right) = \frac{1}{2}\left\| \lambda \right\|_{2}^{2}$

Но в нашей задаче ограничения $y_{i}\left( w^{T}x_{i} + b \right) \geq 1$ линейны, поэтому появляется $\lambda = \sum_{}^{}{\alpha_{i}y_{i}x_{i}}$ , а преобразование Лежандра порождает билинейную форму $\sum_{i,j}^{}{\alpha_{i}\alpha_{j}y_{i}y_{j}x_{i}^{T}x_{j}}$ . Мы перешли из пространства весов в пространство множителей Лагранжа $\alpha$ , где ограничения стали линейными ( $\sum_{}^{}{\alpha_{i}y_{i}} = 0$ ) и неотрицательными ( $\alpha_{i} \geq 0$ ). А когда мы нашли $\alpha$ , веса восстановим по формуле $w = \sum_{}^{}{\alpha_{i}y_{i}x_{i}}$ , а смещение – из $y_{i}\left( w^{T}x_{i} + b \right) = 1$ . Если данные не разделимы линейно, то вводятся штраф за нахлест с коэффициентом.

Более того, в двойственной задаче признаки $x_{j}$ входят только через скалярное произведение $x_{i}^{T}x_{j}$ . Можно заменить $x_{i}^{T}x_{j}$ на произвольную положительно определенную ядерную функцию $K\left( x_{i},x_{j} \right) = \psi\left( x_{i} \right)^{T}\psi\left( x_{j} \right)$ , которая задает метрику и дивергенцию Брегмана в пространстве признаков:

$\max_{\alpha}\left( \sum_{i}^{}\alpha_{i} - \frac{1}{2}\sum_{i,j}^{}{\alpha_{i}\alpha_{j}y_{i}y_{j}K\left( x_{i},x_{j} \right)} \right)$

Этот метод называется методом опорных векторов (SVM). Метод наглядно показывает пользу преобразования Лежандра: нам удалось перейти к задаче в таких координатах, что размерность перестала зависеть от числа признаков (теперь – только от числа объектов), ограничения стали простыми (линейными), а вся нелинейность ушла в ядра (скалярные произведения в пространстве признаков)

Генеративные модели, ELBO

Вернемся к линейной регрессии: что там на самом деле происходило? В начале мы сделали ряд априорных предположений:

о данных — наблюдения независимы и порождены некоторой моделью
об этой модели — таргет зависит от признаков линейно (с весами и некоторым шумом)
об этом шуме – нормальность и гомоскедастичность, $\varepsilon\ \sim\ N\left( 0,\sigma^{2}I \right)$
и об этих весах — тоже нормальность и гомоскедастичность, $\theta_{0}\ \sim\ N\left( 0,\sigma_{0}^{2}I \right)$

Потом мы задумались: какая модель $f_{\widehat{\theta}}(x)$ в этих предположениях наиболее правдоподобно описывает наблюдения? Оказалось:

$\widehat{\theta} = \arg{\max_{\theta}{P\left( Y\ |\ X,\theta \right)}} = \left( X^{T}X + \frac{\sigma^{2}}{\sigma_{0}^{2}}I \right)^{- 1}X^{T}Y$

Построенная модель является дискриминативной: это означает, что по данным нам признакам мы умеем восстанавливать наиболее правдоподобный таргет . О самом распределении признаков мы ничего не знаем. Покрутим определение условной вероятности:

$P\left( Y\ |\ X,\theta \right) = \frac{P\left( X,Y\ |\ \theta \right)}{P\left( X\ |\ \theta \right)}\ \propto \ P\left( X,Y\ |\ \theta \right)$

Предположим, что мы хотим сгенерировать реалистичную картинку, например, изображения котов — нам нужна модель $P\left( img,\ кот\ |\ \theta \right)$ . Дискриминативные модели нам не помогут: они могут лишь сопоставлять данной картинке некоторое число или метку класса: $P\left( кот\ |\ img,\ \theta \right)$ — «на данной картинке кот». Умея генерировать котов, мы можем легко построить их классификатор: $P\left( кот\ |\ img,\ \theta \right) \propto P\left( img,кот\ |\ \theta \right)$ . Но вот обратная задача принципиально сложнее.

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

Все картинки лежат в пространстве 3D-матриц размерности height x width x rgb. Большинство элементов этого пространства выглядят как шумный экран телевизора, но некоторое подмножество из них – это картинки котов. Картинки именно котов генерируются сложным распределением, зависящем от – это скрытые случайные величины (латенты), кодирующие “смысл” картинки. Научимся генерировать латенты и восстанавливать из них картинки котов с помощью модели с параметрами $\theta$ . Теорема Байеса:

$p_{\theta}\left( z\ |\ x \right) = \frac{p_{\theta}(x,z)}{p_{\theta}(x)} = \frac{p_{\theta}\left( x\ |\ z \right)q(z)}{\int_{}^{}{p_{\theta}\left( x\ |\ z \right)q(z)dz}} = \frac{p_{\theta}(x,z)}{\mathbb{E}_{z\ \sim\ q(z)}\left\lbrack p_{\theta}\left( x\ |\ z \right) \right\rbrack}$

Модель $p_{\theta}\left( x\ |\ z \right)$ называется декодер, $p_{\theta}\left( z\ |\ x \right)$ – энкодер, – априор на латенты (не зависит от параметров модели $\theta$ ). Вспомним, что знаменатель формулы Байеса $p_{\theta}(x)$ называется evidence: в нашем примере с картинками – это генератор всех возможных картинок, похожих на наш датасет . Мы хотим, чтобы картинки были максимально правдоподобными, т.е. ищем:

$\max_{\theta}{\ln{p_{\theta}(x)}} = \max_{\theta}{\ln{\int_{}^{}{p_{\theta}(x,z)dz}}} = \max_{\theta}{\ln{\int_{}^{}{e^{\ln{p_{\theta}(x,z)}}dz}}}$

Это выражение — встречавшийся нам много раз logsumexp. Мы знаем, что задача поиска ОМП – двойственная к максимизации энтропии, а logsumexp – это преобразование Лежандра от негэнтропии $\int_{}^{}{q(z)\ln{q(z)dz}}$ . Осталось только подставить выражение $\ln{p_{\theta}(x,z)}$ вместо двойственной координаты:

$\ln{p_{\theta}(x)} = \max_{q}\left( \int_{}^{}{q(z)\ln{p_{\theta}(x,z)}dz} - \int_{}^{}{q(z)\ln{q(z)}dz} \right) = \max_{q}\int_{}^{}{q(z)\ln\frac{p_{\theta}(x,z)}{q(z)}dz}$

Выражение под максимумом называется ELBO – Evidence Lower BOund. Название отсылает нас к разговору о преобразования Лежандра — как поиску наиболее «тугой» нижний оценки на выпуклую функцию. Более того, подогнав формулу к знакомому виду logsumexp, мы получили наглядное выражение для лежандрова зазора:

$ELBO\lbrack q\rbrack = \int_{}^{}{q(z)\ln\frac{p_{\theta}(x,z)}{q(z)}dz} = \int_{}^{}{q(z)\ln\frac{p_{\theta}\left( z\ |\ x \right)p_{\theta}(x)}{q(z)}dz} = \int_{}^{}{q(z)\ln\frac{p_{\theta}\left( z\ |\ x \right)}{q(z)}dz} + \int_{}^{}{q(z)\ln{p_{\theta}(x)}dz} = \left( \int_{}^{}{q(z)dz} \right)\ln{p_{\theta}(x)} - \int_{}^{}{q(z)\ln\frac{q(z)}{p_{\theta}\left( z\ |\ x \right)}dz} = \ln{p_{\theta}(x)} - KL(q(z)\ ||\ p_{\theta}\left( z\ |\ x \right))$

Поскольку $\ln{p_{\theta}(x)}$ не зависит от q, то максимизируя $ELBO\lbrack q\rbrack$ , мы минимизируем KL-дивергенцию между априорным и апостериорным распределениями латентов. Всегда $KL \geq 0$ , причем равенство достигается при $q(z) = p_{\theta}\left( z\ |\ x \right)$ . Более того, мы не можем напрямую считать правдоподобие $\ln{p_{\theta}(x)} = \ln{\int_{}^{}{p_{\theta}\left( x\ |\ z \right)q(z)dz}}$ (например, если оценивать по Монте-Карло, оценка получится смещенной, т.к. логарифм снаружи интеграла). А считать совместное ELBO — можем (логарифм занесли внутрь интеграла, подбираем сами). Снова польза от преобразования Лежандра: перешли в такие координаты, в которых задача решается проще.

Во всех учебных материалах, которые я видел, ELBO выводят методом «хитрых единиц». Ясно, что так делают, чтобы не углубляться в Лежандра и не отвлекаться от хода курса ради одного вывода. Но такое домножение на «хитрые единицы» скрывает более фундаментальную идею и крайне неинтуитивно:

$\ln{p_{\theta}(x)} = \left( \int_{}^{}{q(z)dz} \right) \cdot \ln\left( \frac{p_{\theta}(x,z)}{p_{\theta}\left( z\ |\ x \right)} \cdot \frac{q(z)}{q(z)} \right) = \int_{}^{}{q(z)}\ln\frac{p_{\theta}(x,z)}{q(z)}dz + \int_{}^{}{q(z)\ln\frac{q(z)}{p_{\theta}\left( z\ |\ x \right)}dz} = ELBO\lbrack q\rbrack + KL(q(z)\ ||\ p_{\theta}\left( z\ |\ x \right))$

Теперь остается вопрос: считать и максимизировать ELBO удобнее теоретически – но как?

Пример 1: кластеризация. Пусть у нас есть множество точек, и наша задача – разделить их на группы «похожих». Предположим, что точки порождены из смеси гауссиан, взвешенных категориальным распределением: $p_{\theta}(x) = \sum_{k = 1}^{K}{\pi_{k}N\left( x\ |\ \mu_{k},\Sigma_{k} \right)}$ . Каждый кластер порожден одной конкретной гауссианой, но мы не знаем, какой именно – это наш латент $z\ \sim\ Cat(\pi)$ . Вытянув некоторый латент, модель сопоставляет ему нужную гауссиану и вытягивает точку из нее: $x\ \sim\ N\left( \mu_{z},\Sigma_{z} \right)$ .

Идея: давайте зафиксируем текущие параметры модели $\theta_{t} = \left\{ \pi,\mu,\Sigma \right\}_{t}$ . При них обновим апостериор: с какой вероятностью точка $x_{i}$ порождена компонентой (мягкие принадлежности точек кластерам):

$r_{ik} = p_{\theta_{t}}\left( z_{i} = k\ |\ x_{i} \right) = \frac{\pi_{k}N\left( x_{i}\ |\ \mu_{k},\Sigma_{k} \right)}{\sum_{j = 1}^{K}{\pi_{j}N\left( x_{i}\ |\ \mu_{j},\Sigma_{j} \right)}} = \frac{\pi_{k}\exp\left( - \left( x_{i} - \mu_{k} \right)^{2}/2\sigma^{2} \right)}{\sum_{j}^{}{\pi_{j}\exp\left( - \left( x_{i} - \mu_{j} \right)^{2}/2\sigma^{2} \right)}}$

Обратим внимание: $r_{ik}$ – это softmax с температурой $2\sigma^{2}$ (и весами $\pi$ ). Устремив температуру к нулю, получим жесткую принадлежность кластерам – алгоритм k-means.

Теперь зафиксируем $r_{ik}$ и посчитаем ELBO:

$ELBO = \sum_{i,k}^{}{r_{ik}\ln{p_{\theta}(x,z)}} - \sum_{i,k}^{}{r_{ik}\ln r_{ik}} = \sum_{i,k}^{}{r_{ik}\ln\left( \pi_{k}N\left( x_{i}\ |\ \mu_{k},\Sigma_{k} \right) \right)} - \sum_{i,k}^{}{r_{ik}\ln r_{ik}}$

Максимизируем ELBO по параметрам модели $\theta = \left\{ \pi,\mu,\Sigma \right\}$ :

$\pi_{k} = \frac{\sum_{i}^{}r_{ik}}{n},\ \ \mu_{k} = \frac{\sum_{i}^{}{r_{ik}x_{i}}}{\sum_{i}^{}r_{ik}},\ \ \Sigma_{k} = \frac{\sum_{i}^{}{r_{ik}\left( x_{i} - \mu_{k} \right)\left( x_{i} - \mu_{k} \right)^{T}}}{\sum_{i}^{}r_{ik}}$

Повторяем до сходимости. Спустя много шагов, мы получим модель, умеющую относить точки к кластерам.

Заметим связи:

для каждой точки $x_{i}$ энергия компоненты $\eta_{i,k} = \ln\pi_{k} + \ln{N\left( x_{i}\ |\ \mu_{k},\Sigma_{k} \right)}$ – имеет вид $\psi\left( \eta_{i} \right) = logsumexp$ , а принадлежность $r_{ik} = \left\lbrack \nabla\psi\left( \eta_{i} \right) \right\rbrack_{k}$
первый шаг фиксирует натуральные параметры и переводит моменты/средние сразу в оптимум $\nabla\psi\left( \eta_{i} \right) = p_{\theta_{t}}\left( z_{i} = k\ |\ x_{i} \right)$ . Это шаг зеркального спуска подъема до упора
второй шаг приравнивает достаточные статистики модели ко взвешенным эмпирическим и двигает натуральные параметры. Это шаг с натуральным градиентом, $\nabla^{2}\psi^{*} = Фишер^{- 1}$
полный алгоритм – это покоординатный подъем ELBO: по одним его параметрам + отдельно по другим

Это называется EM-алгоритм:

Е-шаг: зафиксируем параметры модели, при них обновим апостериор $p_{\theta_{t}}\left( z\ |\ x \right)$ . Поскольку мы выставляем $q = p_{\theta_{t}}\left( z\ |\ x \right)$ , то зазор $KL\left( q(z)\ ||\ p_{\theta}\left( z\ |\ x \right) \right) = 0$ , граница ELBO касается $\ln{p_{\theta}(x)}$
М-шаг: фиксируем латент, считаем ELBO и максимизируем по параметрам. Отсюда $\ln{p_{\theta_{t}}(x)} \leq \ln{p_{\theta_{t + 1}}(x)}$ , т.е. правдоподобие растет монотонно

Пример 2: генерация картинок. Пусть у нас есть датасет из картинок, задача — научиться генерировать похожие

Идея: возьмем автоэнкодер — нейросеть из двух частей, энкодера и декодера. Сначала энкодер отобразит картинку в пространство малой размерности, а потом декодер попробует ее оттуда восстановить:

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

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

Аналог расстояния между распределениями – это KL-дивергенция. В лоссе два слагаемых: первое – качество реконструкции, второе – регуляризация за форму распределения латентов.

$- Loss = ELBO = \mathbb{E}_{z\ \sim\ q(z)}\left\lbrack \ln{p_{\theta}\left( x\ |\ z \right)} \right\rbrack - KL(q\left( z\ |\ x \right)\ ||\ N(0,I))$

Такая модель называется вариационный автоэнкодер (VAE) – и работает уже заметно лучше. Более того, архитектура VAE позволяет плавно перемещаться по латентному пространству:

Энкодер выдает моментную сторону латентного апостериора и учится приближать $\nabla\psi$ – отображение из энергии $\ln{p_{\theta}(x,z)}$ в апостериор. Декодер параметризует натуральную сторону выходного распределения и учится аппроксимировать энергию. Само преобразование Лежандра происходит в ELBO

Пример 3: диффузионные модели. Это принципиально другой подкласс генеративных моделей, опирающийся на теорию оптимального управления и стох.дифуры – темы намного более сложные для погружения. Поэтому рассмотрим очень по верхам и лишь покажем, где там Лежандр.

Снова хотим генерировать картинки. Идея:

берем картинку и постепенно добавляем к ней белый шум, пока картинка не исчезнет (прямой процесс)
научим модель восстанавливать исходную картинку из шума (обратный процесс). Обучаем снова через максимизацию ELBO. Обученной модели подаем на вход только шум, а на выходе — реалистичная картинка.

В непрерывном пределе, прямой процесс – СДУ: . Обращение времени – снова диффузия: $dx = \left\lbrack f - g^{2}\nabla_{x}\ln{p_{t}(x)} \right\rbrack dt + g(t)d\overline{w},\ \ \ \frac{dx}{dt} = f - \frac{1}{2}g^{2}\nabla_{x}\ln{p_{t}(x)}$ . Все неизвестное здесь – скор $\nabla_{x}\ln{p_{t}(x)}$ , его и учит нейросеть. Сгенерировать картинку = проинтегрировать обратное уравнение, двигаясь по скору. Посмотрим на генерацию как на задачу оптимального управления: к шумовому СДУ добавим рулящий снос , чтобы к концу процедуры прийти на данные, с квадратичной ценой руления $l(u) = \frac{1}{2}u^{2}$ . Функция ценности подчиняется уравнению Гамильтона-Якоби-Белмана:

$\partial_{t}V + \inf_{u}\left\lbrack \frac{1}{2}u^{2} + (f + gu)\nabla V \right\rbrack + \frac{1}{2}g^{2}\Delta V = 0$

Мы уже встречали уравнение ГЯ в физике, а ГЯБ – его стохастичный аналог. Внутренний инфимум – это буквально преобразование Лежандра от цены руления , взятое в точке $- g\nabla V$ . Решаем, подставляем, получаем «управляющий гамильтониан»: $H(x,p) = fp - \frac{1}{2}g^{2}p^{2},\ \ \ \partial_{t}V + H(x,\nabla V) + \frac{1}{2}g^{2}\Delta V = 0$ . Замена $V = - \ln\psi$ отождествляет $\nabla V$ со скором, решение: $u^{*} = - g\nabla V = g\nabla\ln\psi = g\nabla\ln{p_{t}(x)}$ . Подставляем в управляемое СДУ и получаем диффузионную модель.

Я опустил очень много деталей, но факт в том, что преобразование Лежандра лежит в основе модели

Как-то так. Опять же, приветствую любые комментарии и указания на ошибки/неточности

ссылка на оригинал статьи https://habr.com/ru/articles/1080734/