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

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

Геометрия

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

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

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

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

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

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

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

В-третьих, геометрически это означает, что наклон 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.$

Получилась некоторая «вывернутая дельта-функция»: во всех точках , в нуле 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 снизу (и тогда при логична, ведь прямые и гарантированно пересекутся). Эта идея собрана в неравенстве Фенхеля — Юнга:

$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 - fleft( 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)$ для . В частности, для получим $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)=mgdot{q} + mdot{q}ddot{q}=mdot{q}(g + ddot{q})$ , причем масса положительна, а скорость зануляется только в одной точке (вершине траектории). Остается . Решив этот дифур, получим знакомое $q=q_{0} + v_{0}t - frac{gt^{2}}{2}$ .

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

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

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

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

Вводится кинетическая энергия и потенциальная энергия , их разность называется лагранжианом . Интеграл лагранжиана по времени (вдоль траектории) называется действием:

$Slbrack qrbrack=int_{t_{1}}^{t_{2}}{Lleft( q,dot{q},t right)dt}$

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

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

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

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

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

$Sleftlbrack q_{varepsilon} rightrbrack=S(varepsilon)=int_{t_{1}}^{t_{2}}{Lleft( q + varepsiloneta,dot{q} + varepsilondot{eta},t right)dt}$

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

$0=frac{delta S}{delta q}=frac{dS}{dvarepsilon}=int_{t_{1}}^{t_{2}}{frac{dLleft( q + varepsiloneta,dot{q} + varepsilondot{eta},t right)}{dvarepsilon}dt}=int_{t_{1}}^{t_{2}}{leftlbrack frac{partial L}{partial q}frac{partial q_{varepsilon}}{partialvarepsilon} + frac{partial L}{partialdot{q}}frac{partial{dot{q}}_{varepsilon}}{partialvarepsilon} rightrbrack dt}=int_{t_{1}}^{t_{2}}{leftlbrack frac{partial L}{partial q}eta rightrbrack dt} + int_{t_{1}}^{t_{2}}{leftlbrack frac{partial L}{partialdot{q}}dot{eta} rightrbrack dt}$

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

$0=frac{dS}{dvarepsilon}=int_{t_{1}}^{t_{2}}{leftlbrack frac{partial L}{partial q}eta rightrbrack dt} + leftlbrack frac{partial L}{partialdot{q}}eta rightrbrack_{t_{1}}^{t_{2}} - int_{t_{1}}^{t_{2}}{frac{d}{dt}left( frac{partial L}{partialdot{q}} right)eta dt}=int_{t_{1}}^{t_{2}}{leftlbrack frac{partial L}{partial q} - frac{d}{dt}left( frac{partial L}{partialdot{q}} right) rightrbracketa dt}$

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

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

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

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

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

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

$Slbrack q,prbrack=int_{t_{1}}^{t_{2}}{leftlbrack pdot{q} - H(q,p) rightrbrack dt}, 0=frac{dS}{dvarepsilon}=int_{t_{1}}^{t_{2}}{leftlbrack dot{q} - frac{partial H}{partial p} rightrbrackxi dt} - int_{t_{1}}^{t_{2}}{leftlbrack dot{p} + frac{partial H}{partial q} rightrbracketa dt}$

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

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

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

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

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

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

$dH=dot{q}dp + pddot{q} - frac{partial L}{partial q}dq - frac{partial L}{partialdot{q}}ddot{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}{partialdot{q}} right)ddot{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}{partialdot{q}} right)=dot{p}, т.е.dot{p}=- frac{partial H}{partial q}$

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

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

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

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

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

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

Дополнительно рассмотрим случай, когда какой-то нельзя выразить через свой из $p=frac{partial L}{partialdot{q}}$ (как мы делали). Это означает $frac{dp}{ddot{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( pdot{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.$

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

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

$frac{partial S}{partial t} + Hleft( q,frac{partial S}{partial q},t right)=0,$

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

Оптимизация

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

Чаще преобразование Лежандра записывают в векторном виде: $f^{*}(p)=sup_{x}left( leftlangle x,p rightrangle - f(x) right)$ . Для гладких функций: $f^{*}(p)=leftlangle x,p rightrangle - 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 (т.е. таких, что второе слагаемое зануляется) можно взять любой множитель – все равно в скобке ноль. Но теперь давайте рассмотрим не только допустимые 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}leftlbrack f(x) - left( A^{T}lambda right)^{T}x rightrbrack=lambda^{T}b - sup_{x}leftlbrack left( A^{T}lambda right)^{T}x - f(x) rightrbrack=lambda^{T}b - f^{*}left( A^{T}lambda right)$

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

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

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

$x_{k + 1}=x_{k} - eta cdot nabla fleft( x_{k} right)$

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

$f(x) approx fleft( x_{k} right) + leftlangle nabla fleft( x_{k} right), x - x_{k} rightrangle$ $x_{k + 1}=arg{minleftlbrack fleft( x_{k} right) + leftlangle nabla fleft( x_{k} right), x - x_{k} rightrangle + frac{1}{2eta}left| x - x_{k} right|^{2} rightrbrack}=arg{min_{x}leftlbrack etaleftlangle nabla fleft( x_{k} right),x rightrangle + frac{1}{2}left| x - x_{k} right|^{2} rightrbrack}$

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

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

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

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

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

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

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

$x_{k + 1}=x_{k} - eta cdot nabla fleft( x_{k} right)=x_{k} - eta cdot nablapsi^{*}(partial f)$

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

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

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

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

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

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

$x_{k + 1}=nablapsi^{*}left( nablapsileft( x_{k} right) - eta cdot nabla fleft( x_{k} right) right)$

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

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

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

$x_{k + 1}=arg{min_{x}leftlbrack etaleftlangle nabla fleft( x_{k} right),x rightrangle + D_{psi}left( x,x_{k} right) rightrbrack}, где D_{psi}(x,y)=psi(x) - psi(y) - leftlangle nablapsi(y),x - y rightrangle$

Функция $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) - leftlangle nablapsi(y),x - y rightrangle=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}logfrac{x_{i}}{y_{i}}}=KLleft( x||y right)$

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

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

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

$f(x) approx fleft( x_{k} right) + leftlangle nabla fleft( x_{k} right), x - x_{k} rightrangle + frac{1}{2}left( x - x_{k} right)^{T}nabla^{2}fleft( x_{k} right)left( x - x_{k} right)$ $x_{k + 1}=arg{min{leftlbrack fleft( x_{k} right) + leftlangle nabla fleft( x_{k} right), x - x_{k} rightrangle + frac{1}{2}left( x - x_{k} right)^{T}nabla^{2}fleft( x_{k} right)left( x - x_{k} right) rightrbrack=x_{k} - leftlbrack nabla^{2}fleft( x_{k} right) rightrbrack^{- 1}partial f}}$

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

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

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

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

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

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

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

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

Что такое «интенсивные» и «экстенсивные»? Пусть система увеличилась в раз. Тогда – экстенсивная (объем V, энтропия S), а – интенсивная (температура 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 частицу: . Важно, что s стала интенсивной. Тогда $f(T,v)=inf_{s}left( u(s,v) - Ts right)$ , и Лежандр произвел обмен двух интенсивных величин. Тот факт, что часто лежандрова пара состоит из экстенсивной + интенсивной величины, это просто следствие однородности взятых функций.

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

Энтропия

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

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

$S=kln{Omega(E)}=k cdot left( - frac{Omega}{Omega} right) cdot left( - lnOmega right)=- ksum_{i=1}^{Omega}{frac{1}{Omega}lnfrac{1}{Omega}}=- ksum_{i=1}^{Omega}{pln p}$

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

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

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

Посмотрим на вероятность $p_{i}$ микросостояния с энергией $E_{i}$ в системе из ящика и резервуара. Все состояния равновероятны, поэтому вероятность пропорциональна числу способов его реализовать: $p_{i} sim Omegaleft( E - E_{i} right)$ , причем энергия ящика много меньше энергии резервуара ( $E gg E_{i}$ ). Разложив определение энтропии по Тейлору, получим $frac{S}{k}=ln{Omegaleft( E - E_{i} right)} approx ln{Omega(E)} - frac{partiallnOmega}{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}simln{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=- ksum_{}^{}{p_{i}ln p_{i}}=- ksum_{}^{}{p_{i}left( - ln Z - frac{E_{i}}{kT} right)}=kln Z + frac{U}{T} Rightarrow U - TS=- kTln Z$

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

$F=- frac{ln Z}{beta}, где Z=sum_{i}^{}e^{- beta E_{i}} и beta=frac{1}{kT} Rightarrow F=- kTlnleftlbrack sum_{i}^{}e^{- E_{i}/kT} rightrbrack$

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

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

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

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

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

$Z(beta)=expleftlbrack sup_{E}left( frac{S(E)}{k} - beta E right) rightrbrack$

Продифференцировав $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). Сама лежандрова пара – это функция по переменной и функция $- frac{S}{k}$ по переменной .

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

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

Статистика

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

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

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

энтропия имеет вид $- sum_{}^{}{pln 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}$

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

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

Градиент негэнтропии не имеет какого-то устоявшегося названия, но с точностью до константы – это логиты / логарифмические шансы: $theta_{i}=lnfrac{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( leftlangle theta,p rightrangle - sum_{i}^{}{p_{i}ln p_{i}} right)$

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

$L=sum_{i}^{}{p_{i}theta_{i}} - sum_{i}^{}{p_{i}ln p_{i}} + lambdaleft( 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}$

Отдельно отметим, что это полностью соответствует определению преобразования Лежандра, т.к. множитель при зануляется на всех из рассматриваемого пространства. Саму подбираем из нормировки $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 громоздкая выкладкаbackslashtext{=}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)=- Hlbrack prbrack=int_{}^{}{p(x)ln{p(x)}dx}, nablaleft( - Hlbrack prbrack right)=ln{p(x)} + 1, nabla^{2}left( - Hlbrack prbrack right)=frac{delta(x - y)}{p(x)}$

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

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

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

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

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

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

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

Что означает сам супремум: мы максимизировали энтропию при фиксированном матоже параметра $mathbb{E}_{p}theta=leftlangle theta,p rightrangle=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}{pleft( x middle| mu right)}}=arg{max_{mu}{lnleft( pleft( x middle| mu right) right)}}$ , т.к. логарифм строго возрастает и не сдвигает аргмаксимум. Вероятность называют правдоподобием, функцию – логарифмом правдоподобия, ее градиент – скором, а саму – оценкой максимального правдоподобия (ОМП). Посчитаем:

$l(mu)=ln{pleft( x middle| mu right)}=ln{prod_{i=1}^{n}{pleft( x_{i} middle| mu right)}}=sum_{i=1}^{n}{ln{pleft( x_{i} middle| mu right)}}=sum_{i=1}^{n}{lnleft( mu^{x_{i}}(1 - mu)^{1 - x_{i}} right)}=lnmu cdot sum_{}^{}x_{i} + ln(1 - mu) cdot left( n - sum_{}^{}x_{i} right)=mnlnmu + (n - mn)ln(1 - mu)$ $0=frac{d}{dmu}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=lnfrac{mu}{1 - mu}$ . Оптимальный $widehat{theta}=lnfrac{m}{1 - m}$

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

$L=- sum_{}^{}{pln p} + thetaleft( sum_{}^{}{px} - m right) + lambdaleft( 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}=lnfrac{m}{1 - m}$ – мы получили ровно такой же, как для ОМП. К тому же, укрепилась наша интуиция о записи параметра в виде логита.

Пример 2: нормальное распределение $X_{i} sim Nleft( 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}$

ОМП:

$lleft( mu,sigma^{2} right)=- frac{n}{2}lnleft( 2pisigma^{2} right) - frac{1}{2sigma^{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_{}^{}{pln p} + theta_{1}left( m_{1} - int_{}^{}{pxdx} right) + theta_{2}left( m_{2}^{2} - int_{}^{}{px^{2}dx} right) + lambdaleft( 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 статистики и обновлять их итеративно – и как будто в этих статистиках было достаточно информации, чтобы описать всё распределение.
плотности выражались через экспоненту, если записывать их через правильные (с точки зрения преобразования Лежандра) параметры. Именно эти «естественные» параметры сопряжены с наблюдаемыми средними значениями тех «достаточных» статистик

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

$pleft( 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}}}}$

Названия компонент: – натуральный параметр, – достаточная статистика, – базовая мера. Показатели экспонент линейны по , а нормировка имеет знакомый нам вид 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 + thetawidehat{mu} - f(theta), где widehat{mu}=frac{1}{n}sum_{}^{}{Tleft( x_{i} right)}$

Сразу видно, что А max по – буквально определение преобразования Лежандра:

$f^{*}left( widehat{mu} right)=max_{theta}left( thetawidehat{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}{KLleft( 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}Xmathbb{E}X^{T}$ )

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

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

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

$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 аппроксимирует диагональной матрицей (плюс пара инженерных доработок). Подробнее про применения в ML поговорим позже.

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

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

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

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

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

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

$s(x)=lnfrac{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)lnfrac{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)lnfrac{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)$ . Для любого индикатор мажорируется экспонентой (неравенство Чернова):

$alpha_{n} leq e^{- lambda ngamma}mathbb{E}_{0}left( e^{lambdasum_{}^{}{sleft( x_{i} right)}} right)=e^{- lambda ngamma}left( mathbb{E}_{0}e^{lambda s} right)^{n}=e^{- nleft( lambdagamma - 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^{*}$ от лог-статсуммы / производящей функции моментов называется функцией уклонения:

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

Заметим, что – это тот же 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)}$

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

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

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

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

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

ML

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

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

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

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

– апостериор (представление о параметре, уточненное по данным)

– априор (текущая гипотеза), в начале предполагаем, что шар приземляется равномерно, т.е. имеет распределения

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

– называется evidence. По сути, нам интересен только числитель, а знаменатель - это нормировочная константа, которую на практике часто оценивают или подгоняют

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

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

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

$Pleft( 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)$

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

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

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

$pleft( 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:

$pleft( X middle| theta right)=prod_{}^{}{hleft( x_{i} right)} cdot expleft( thetasum_{}^{}{Tleft( x_{i} right)} - nf(theta) right)$

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

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

$g(alpha,nu)=ln{int_{}^{}{expleft( thetaalpha - nu f(theta) right)dtheta}}$

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

$g(alpha,nu) approx sup_{theta}left( thetaalpha - 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 Nleft( theta^{T}x_{i},sigma^{2} right)$ , и попробуем подобрать оптимальные веса – параметры модели. Сразу отметим, что часто модель записывают как $y_{i}=theta^{T}x_{i} + b_{i}$ , но это эквивалентно добавлению признака из всех единиц и включению $b_{i}$ в вектор весов – так что будем писать в простом виде

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

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

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

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

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

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

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

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

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

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

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

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

$sum_{i=1}^{n}left( y_{i} - theta^{T}x_{i} right)^{2}=Y^{T}Y - 2theta^{T}X^{T}Y + theta^{T}X^{T}Xtheta$ $nabla_{theta}ln{pleft( Y | X,theta right)}=- frac{1}{sigma^{2}}left( X^{T}Xtheta - 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 Nleft( theta_{0},Sigma_{0} right), p(theta)=frac{1}{(2pi)^{d/2}left| Sigma_{0} right|^{1/2} }expleft( - frac{1}{2}left( theta - theta_{0} right)^{T}Sigma_{0}^{- 1}left( theta - theta_{0} right) right)$

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

$pleft( theta | X,Y right)=frac{pleft( Y | X,theta right) cdot p(theta)}{pleft( Y | X right)} propto expleft( - frac{1}{2sigma^{2}}(Y - Xtheta)^{T}(Y - Xtheta) - 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}Sigmaeta$ , где $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}=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-лосса и правдоподобия чуть глубже. Пусть вектор ошибки , тогда $MSE(theta)=frac{1}{2}left| varepsilon right|_{2}^{2}$ . Перейдем к двойственной задаче через преобразование Лежандра по . Снова вспомним, что $psi(x)=frac{1}{2}x^{2}$ – неподвижная точка преобразования Лежандра, поэтому $MSE^{*}(lambda)=frac{1}{2}left| lambda right|_{2}^{2},$ где – двойственная переменная. Решение двойственной задачи имеет тот же вид. Более того, теперь решение без регуляризации $widehat{theta}=left( X^{T}X right)^{- 1}X^{T}Y$ явно выходит как частный случай, при ограничении $X^{T}lambda=0$

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

$theta_{t + 1}=theta_{t} - etanabla MSEleft( theta_{t} right), где nabla MSEleft( theta_{t} right)=X^{T}Xtheta_{t} - X^{T}Y$

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

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

$I(theta)mathbb{=- E}nabla^{2}l(theta)=nablafrac{1}{sigma^{2}}left( X^{T}Xtheta - 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} - etasigma^{2}left( X^{T}X right)^{- 1}left( X^{T}Xtheta_{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

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

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

forward-pass: нейрон получает на вход точку . Линейная операция переводит его в логит $z=theta^{T}x$ – элемент сопряженного к вероятностям пространства, активация (как ) переводит логит в первичное пространство вероятностей/средних. Считается лосс – лежандров зазор
backward-pass: расчет градиентов весов по цепочному правилу, от лосса назад по всем слоям. Градиент лосса по логиту $frac{partial L}{partial z}=widehat{y} - 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}}$ . Мы перешли из пространства весов в пространство множителей Лагранжа , где ограничения стали линейными ( $sum_{}^{}{alpha_{i}y_{i}}=0$ ) и неотрицательными ( $alpha_{i} geq 0$ ). А когда мы нашли , веса восстановим по формуле $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}$ на произвольную положительно определенную ядерную функцию $Kleft( x_{i},x_{j} right)=psileft( x_{i} right)^{T}psileft( x_{j} right)$ , которая задает метрику и дивергенцию Брегмана в пространстве признаков:

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

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

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

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

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

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

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

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

$Pleft( Y | X,theta right)=frac{Pleft( X,Y | theta right)}{Pleft( X | theta right)} propto Pleft( X,Y | theta right)$

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

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

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

$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)}leftlbrack p_{theta}left( x | z right) rightrbrack}$

Модель $p_{theta}left( x | z right)$ называется декодер, $p_{theta}left( z | x right)$ – энкодер, – априор на латенты (не зависит от параметров модели ). Вспомним, что знаменатель формулы Байеса $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)lnfrac{p_{theta}(x,z)}{q(z)}dz}$

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

$ELBOlbrack qrbrack=int_{}^{}{q(z)lnfrac{p_{theta}(x,z)}{q(z)}dz}=int_{}^{}{q(z)lnfrac{p_{theta}left( z | x right)p_{theta}(x)}{q(z)}dz}=int_{}^{}{q(z)lnfrac{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)lnfrac{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, то максимизируя , мы минимизируем KL-дивергенцию между априорным и апостериорным распределениями латентов. Всегда , причем равенство достигается при $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 lnleft( frac{p_{theta}(x,z)}{p_{theta}left( z | x right)} cdot frac{q(z)}{q(z)} right)=int_{}^{}{q(z)}lnfrac{p_{theta}(x,z)}{q(z)}dz + int_{}^{}{q(z)lnfrac{q(z)}{p_{theta}left( z | x right)}dz}=ELBOlbrack qrbrack + KL(q(z) || p_{theta}left( z | x right))$

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

Пример 1: кластеризация. Пусть у нас есть множество точек, и наша задача – разделить их на группы «похожих». Предположим, что точки порождены из смеси гауссиан, взвешенных категориальным распределением: $p_{theta}(x)=sum_{k=1}^{K}{pi_{k}Nleft( x | mu_{k},Sigma_{k} right)}$ . Каждый кластер порожден одной конкретной гауссианой, но мы не знаем, какой именно – это наш латент . Вытянув некоторый латент, модель сопоставляет ему нужную гауссиану и вытягивает точку из нее: $x sim Nleft( 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}Nleft( x_{i} | mu_{k},Sigma_{k} right)}{sum_{j=1}^{K}{pi_{j}Nleft( x_{i} | mu_{j},Sigma_{j} right)}}=frac{pi_{k}expleft( - left( x_{i} - mu_{k} right)^{2}/2sigma^{2} right)}{sum_{j}^{}{pi_{j}expleft( - left( x_{i} - mu_{j} right)^{2}/2sigma^{2} right)}}$

Обратим внимание: $r_{ik}$ – это softmax с температурой $2sigma^{2}$ (и весами ). Устремив температуру к нулю, получим жесткую принадлежность кластерам – алгоритм 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}lnleft( pi_{k}Nleft( x_{i} | mu_{k},Sigma_{k} right) right)} - sum_{i,k}^{}{r_{ik}ln r_{ik}}$

Максимизируем ELBO по параметрам модели :

$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}=lnpi_{k} + ln{Nleft( x_{i} | mu_{k},Sigma_{k} right)}$ – имеет вид $psileft( eta_{i} right)=logsumexp$ , а принадлежность $r_{ik}=leftlbrack nablapsileft( eta_{i} right) rightrbrack_{k}$
первый шаг фиксирует натуральные параметры и переводит моменты/средние сразу в оптимум $nablapsileft( 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)$ , то зазор $KLleft( 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: генерация картинок. Пусть у нас есть датасет из картинок, задача - научиться генерировать похожие

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

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

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

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

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

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

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

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

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

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

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

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

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

В непрерывном пределе, прямой процесс – СДУ: . Обращение времени – снова диффузия: $dx=leftlbrack f - g^{2}nabla_{x}ln{p_{t}(x)} rightrbrack dt + g(t)doverline{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}leftlbrack frac{1}{2}u^{2} + (f + gu)nabla V rightrbrack + frac{1}{2}g^{2}Delta V=0$

Мы уже встречали уравнение ГЯ в физике, а ГЯБ – его стохастичный аналог. Внутренний инфимум – это буквально преобразование Лежандра от цены руления , взятое в точке . Решаем, подставляем, получаем «управляющий гамильтониан»: $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$ . Замена отождествляет со скором, решение: $u^{*}=- gnabla V=gnablalnpsi=gnablaln{p_{t}(x)}$ . Подставляем в управляемое СДУ и получаем диффузионную модель.

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

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

Автор: Anatolywww

Источник ^[1]

Сайт-источник PVSM.RU: https://www.pvsm.ru

Путь до страницы источника: https://www.pvsm.ru/fizika/457673

Ссылки в тексте:

[1] Источник: https://habr.com/ru/articles/1080734/?utm_source=habrahabr&utm_medium=rss&utm_campaign=1080734

Нажмите здесь для печати.