О чём эта запись?

В Москве мы проводили для студентов МГУ, ВШЭ, МФТИ, ИППИ совместный семинар по статистическому кластерному анализу. Недавно я делал доклад на тему алгоритма суммы-произведения и графа сомножителей (sum-product algorithm and factor graphs). Я никогда не занимался этой темой, но удачным образом доклад породил интересную дискуссию про то, как используются байесовские методы и статистическая физика в обработке изображений. Поверхностный поиск выдал мне несколько книг по теме (обе на русском): (1) случайное поле Гиббса и (2) случайные поля и марковские цепи в обработке изображений. Первая книга написана полностью на языке теоретической физики, вторая больше про обработку изображений и про байесовскую статистику.

Текущая заметка основана на статье «Factor Graphs and the Sum-Product Algorithm (2001)» за авторством Kschischang, Frey, Loeliger. Изначально мой интерес в этой статье был основан на том, что мы хотели понять, как работает известный алгоритм кластеризации «Affinity Propagation». Авторы Frey и Dueck придумали этот алгоритм и опубликовали про него несколько работ. Алгоритм Affinity Propagation весьма схож с алгоритмом В.Г.Спокойного Adaptive Weight Clustering, над которым мы тогда работали, и мы хотели узнать «врага в лицо».

Оказалось, что идея графа сомножителей является куда более общей, чем просто один алгоритм кластеризации, и содержит красивые математические идеи. Чем отличается эта заметка от оригинальной статьи? Я постараюсь заострить внимание на моментах, которые вызвали наибольший интерес в течение семинара, и пропущу доказательства и технические подробности. Стоит заметить, что статья написана очень хорошим и понятным языком, так что мне почти ничего не нужно делать. Часть материала просто скопирована дословно.

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

Алгоритм суммы-произведения является обобщением нескольких других известных алгоритмов: Forward-Backward Algorithm, алгоритма Витерби для марковских цепей, фильтра Калмана, быстрого преобразования Фурье и нескольких алгоритмов из теории линейного кодирования. Это довольно длинный список, но я постараюсь объяснить, как довольно несложная модель передачи сообщений в двудольном графе может объединять все эти понятия. Данная заметка покрывает только первую часть доклада, поэтому мы обсудим только одно приложение.

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

Давайте начнём с такой задачи. Пусть \( x_1, \ldots, x_n \) это набор переменных, где каждая переменная принимает значения только в конечном алфавите \( \Sigma \). Пусть \( g(x_1, \ldots, x_n) \) некоторая функция, принимающая вещественные значения с некоторой дополнительной внутренней структурой, которую мы обсудим позже. Мы хотим вычислить так называемые маргинальные функции: \[ g_i(x_i) \overset{def}= \sum_{x_1} \cdots \sum_{x_{i-1}} \sum_{x_{i+1}} \cdots \sum_{x_{n}} g(x_1, \ldots, x_n), \] где значение \( g_i(a) \) получается суммированием функции \( g(x_1, \ldots, x_n) \) по всем значениях переменных с условием \( x_i = a \). Авторы оригинальной статьи придумали специальное обозначение для сумм такого рода: \[ g_i (x_i) = \sum_{\sim \{ x_i \}} g(x_1, \ldots, x_n) \] Например, если \( h \) это функция от трёх переменных, то \[ h_2(x_2) = \sum_{\sim \{ x_2 \}} h(x_1, x_2, x_3) = \sum_{x_1} \sum_{x_3} h(x_1, x_2, x_3). \] Вообще говоря, необходимо экспоненциальное число шагов для того, чтобы вычислить такую функцию. Нам нужны некоторые дополнительные предположения про устройство функции \( g \). Давайте предположим, что функция \( g(x_1, \ldots, x_n) \) раскладывается в произведение нескольких локальных функций, то есть \[ g(x_1, \ldots, x_n) = \prod_{j \in J} f_j (X_j), \] где каждое \( X_j \) является подмножеством \( \{ x_1, \ldots, x_n \} \), и \( J \) является дискретным множеством индексов. Ниже мы рассмотрим пример.

Почему эта задача является интересной? Почему мы рассматриваем именно такую структуру? Часто ли на практике встречаются функции такого вида?

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

Определение.

Графом сомножителей называется двудольный граф, который представляет структуру факторизации. Первая доля содержит вершины, соответствующие переменны, вторая доля — функциям. Между вершинами, соответствующими переменной \( x_i \) и функции \( f_j \) (последняя называется вершиной-сомножителем) есть ребро тогда и только тогда, когда \( x_i \) входит в функцию \( f_j \) как аргумент.

Пример.

Пусть \[ g(x_1, x_2, x_3, x_4, x_5) = f_A(x_1) f_B(x_2) f_C(x_1, x_2, x_3) f_D(x_3, x_4) f_E(x_3, x_5). \] Тогда соответствующий граф сомножителей изображён на рисунке:

Дерево выражения

Если граф сомножителей не содержит циклов, то его можно представить в виде дерева.

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

Описанная ситуация может быть использована для того, чтобы упростить вычисления используя закон дистрибутивности. Пусть мы хотим найти \( g_1(x_1) \) из предыдущего примера. Тогда сперва мы вычисляем маргинальные суммы для вершин \( x_4, x_5 \), затем «передаём сообщения» от \( f_D, f_E \) к вершине \( x_3 \) (позже мы поймём, что именно означает передать сообщение), затем вычисляем маргинальные суммы для \( x_3 \), и так далее.

Более формально, используя закон дистрибутивности, маргинальная сумма для \( g_1(x_1) \) может быть представлена в виде \[ g_1 (x_1) = f_A (x_1) \Big( \sum_{x_2} f_B(x_2) \Big( \sum_{x_3} f_C(x_1, x_2, x_3) \Big( \sum_{x_4} f_D(x_3, x_4) \Big) \Big( \sum_{x_5} f_E(x_3, x_5) \Big) \Big) \Big) \] Более того, используя обозначение исключающего суммирования \( \sum_{\sim \{ x \}} \) вместо \(\sum_{x}\), та же самая сумма может быть представлена в виде \[ g_1(x_1) = f_A (x_1) \sum_{\sim \{x_1 \}} \Big( f_B(x_2)f_C(x_1, x_2, x_3)
\Big( \sum_{\sim \{x_3\}} f_D(x_3, x_4) \Big) \Big( \sum_{\sim \{x_3\}} f_E(x_3, x_5) \Big) \Big) \]

Проверка корректности выражения выше может занять некоторое время.

Упражнение.

Напишите аналогичную декомпозицию для \( g_3 (x_3) \).

Предложение.

Если граф сомножителей не содержит циклов, то граф сомножителей не только кодирует представление глобальной функции, но также даёт алгоритм для быстрого вычисления маргинальных функций.

Упражнение.

Верно ли, что вычисления ускоряются по сравнению с методом «грубой силы»?

Вычисление маргинальных функций

До сих пор, дерево было «подвешено» за вершину \( x_1 \), мы мы могли вычислить маргинальную сумму \( g_1(x_1) \). Если мы хотим вычислить, скажем, маргинальную сумму \( g_3(x_3) \), текущая техника позволяет только подвесть заново это дерево за другую вершину и проделать всю процедуру разложения заново. Сперва мы поймём, как сделать процесс автоматическим для одной маргинальной суммы, а затем модифицируем процедуру таким образом, чтобы считать все маргинальные суммы.

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

Что является сообщением? Рассмотрим ребро \( \{x, f\} \), соединяющее переменную \( x\) с вершиной-сомножителем \( f\). Независимо от того, куда смотрит ребро, в качестве сообщения мы будем рассматривать некоторую функцию от одной переменной \( x\), скажем, \( \varphi_{ \{x, f\} }(x) \). Каждая вершина получает сообщения, и затем отсылает сообщения.

Как уже говорилось, мы инициализируем процесс в листьях дерева. В течение инициализации, каждая листовая вершина, соответствующая переменной \( x\), посылает тождественную функцию (\( \varphi(x) \equiv 1 \)) в качестве сообщения своему родителю, а каждая листовая вершина-сомножитель \( f\) посылает описание \( f \) своему родителю. Заметим, что если \( f\) является листовой вершиной, то функция \( f\) зависит ровно от одной переменной.

Введём следующие правила обновления состояний:

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

Вершина-сомножитель \( f \) с родителем \( x \) формирует произведение \( f \) с сообщениями, полученными от потомков, а затем действует на результат оператором \( \sum_{\sim \{ x \}} \) исключающего суммирования.

Пусть \( n(v) \) обозначает множество соседей заданной вершины \( v\). Вычисление сообщений можно изобразить следующим образом:

от переменной к локальной функции

\[ \mu_{x \to f} (x) = \prod_{h \in n(x) \backslash \{f\}} \mu_{h \to x} (x) \]

от локальной функции к переменной

\[ \mu_{f \to x} (x) = \sum_{\sim \{x\}} \Big( f(X) \prod_{y \in n(f) \backslash \{x\}} \mu_{y \to f} (y) \Big) \]

Согласно этой схеме, вычисление заканчивается в корневой вершине \( x_i \), потому что дерево подвешено за вершину \( x_i \). Конечная маргинальная функция \( g_i (x_i) \) получается как произведение всех сообщений, полученных в \( x_i \).

Упражнение.

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

Описанная процедура хороша, если мы хотим только одну маргинальную сумму. Давайте подумаем, как можно модифицировать её для нескольких маргиналов.

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

Теперь мы снова инициируем процесс передачи сообщений в листовых вершинах. И снова, каждая вершина остаётся в состоянии ожидания до тех пор, пока не получит сообщения от всех смежных вершин кроме одной. Как только сообщения получены, эта вершина может вычислить сообщение и послать его по оставшемуся ребру своему соседу. После отправки сообщения, данная вершина возвращается в состояние ожидания, и ждёт обратного сообщения по этому ребру. Как только сообщение получено, вершина может вычислить и послать сообщения остальным соседям. Алгоритм завершается как только по каждому ребру было передано по два сообщения, по одному в каждом направлении.

Пример процедуры передачи сообщений

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

Если подробнее, сообщения генерируются таким образом:

Шаг 1.

\begin{align} \mu_{A \to 1} (x_1) &= \sum_{\sim \{x_1 \}} f_A(x_1) = f_A (x_1) \\ \mu_{B \to 2} (x_2) &= \sum_{\sim \{x_2 \}} f_B(x_2) = f_B (x_2) \\ \mu_{4 \to D} (x_4) &= 1 \\ \mu_{5 \to E} (x_5) &= 1. \end{align}

Шаг 2.

\begin{align} \mu_{1 \to С} (x_1) &= \mu_{A \to 1} (x_1) \\ \mu_{2 \to С} (x_2) &= \mu_{B \to 2} (x_2) \\ \mu_{D \to 3} (x_3) &= \sum_{\sim \{x_3\}} \mu_{4 \to D} (x_4) f_D (x_3, x_4) \\ \mu_{E \to 3} (x_3) &= \sum_{\sim \{x_3\}} \mu_{5 \to E} (x_5) f_E (x_3, x_5) . \end{align}

Шаг 3.

\begin{align} \mu_{C \to 3} (x_3) &= \sum_{\sim \{x_3\}} \mu_{1 \to C} (x_1) \mu_{2 \to C} (x_2) f_C (x_1, x_2, x_3) \\ \mu_{3 \to C} (x_3) &= \mu_{D \to 3} (x_3) \mu_{E \to 3} (x_3) \end{align}

Шаг 4.

\begin{align} \mu_{C \to 1} (x_1) &= \sum_{\sim \{x_1\}} \mu_{3 \to C} (x_3) \mu_{2 \to C} (x_2) f_C (x_1, x_2, x_3) \\ \mu_{C \to 2} (x_2) &= \sum_{\sim \{x_2\}} \mu_{3 \to C} (x_3) \mu_{1 \to C} (x_1) f_C (x_1, x_2, x_3) \\ \mu_{3 \to D} (x_3) &= \mu_{C \to 3} (x_3) \mu_{E \to 3} (x_3) \\ \mu_{3 \to E} (x_3) &= \mu_{C \to 3} (x_3) \mu_{D \to 3} (x_3). \end{align}

Шаг 5.

\begin{align} \mu_{1 \to A} (x_1) &= \mu_{C \to 1} (x_1) \\ \mu_{2 \to B} (x_2) &= \mu_{C \to 2} (x_2) \\ \mu_{D \to 4} (x_4) &= \sum_{\sim \{x_4\}} \mu_{3 \to D} (x_3) f_D (x_3, x_4) \\ \mu_{E \to 4} (x_5) &= \sum_{\sim \{x_5\}} \mu_{3 \to E} (x_3) f_E (x_3, x_5) \\ \end{align}

Завершение.

\begin{align} g_1(x_1) &= \mu_{A \to 1} (x_1) \mu_{C \to 1} (x_1)\\ g_2(x_2) &= \mu_{B \to 2} (x_2) \mu_{C \to 2} (x_2)\\ g_3(x_3) &= \mu_{C \to 3} (x_3) \mu_{D \to 3} (x_3) \mu_{E \to 3} (x_3)\\ g_4(x_4) &= \mu_{D \to 4} (x_4) \\ g_5(x_5) &= \mu_{E \to 5} (x_5) \end{align}

Некоторые замечания

Трюк с полукольцом

Можно заметить, что маргинальные суммы используют две операции: \( \{+, \times\} \) и эти две операции удовлетворяют закону дистрибутивности. Фактически, можно использовать любые две операции, которые удовлетворяют закону дистрибутивности, как это выполнено в случае коммутативного полукольца. Например, если мы хотим найти максимум \( \max f(x_1, \ldots, x_n) \), можно использовать набор операций \( \{\max, \times\} \), где \( \max \) теперь играет роль суммы. Это позволяет сформулировать новый оптимизационный алгоритм в терминах графов сомножителей!

Графы сомножителей с циклами

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

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

Важно подчеркнуть, что теперь мы не инициализируем вершины, а вместо этого инициализируем рёбра. Теперь, в течение инициализации мы предполагаем, что каждое ребро заряжено сообщением с единицей.

Упражнение.

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

Итак, граф сомножителей полученный из функции \( g \), содержит циклы. Какие гарантии можно получить для такого алгоритма? Оказывается, что каждая итерация алгоритма суммы-произведения минимизирует так называемую дивергенцию Кульбака-Лейблера. То есть, существует полуинвариант, который минимизируется нашей процедурой. Давайте поймём, что происходит, шаг за шагом.

Метод Бете

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

Предложение (Аппроксимация среднего поля)

Обозначим через \( N(a) \) множество аргументов функции с индексом \( a\), и через \( d_i \) мощность множества соседей \( x_i \).

Предположим, что \( p(\boldsymbol x) \) задаёт вероятностное распределение в виде произведения, где \[ p(\boldsymbol x) = \prod_{a = 1}^{M} q_a(\boldsymbol x_{N(a)}) \] и соответствующие граф не имеет циклов. Тогда это распределение можно представить в виде \[ p(\boldsymbol x) = \dfrac{\prod_{a=1}^M p_a(\boldsymbol x_{N(a)})}{\prod_{i=1}^N (p_i(x_i))^{d_i - 1}},
\] где функции \( p_n (x_n) \) и \( p_{a}(\boldsymbol x_{N(a)}) \) подчиняются ограничениям \begin{align} \sum_{x_n} p_n (x_n) &= 1 \quad \forall n \in [1, N], \\ \sum_{\boldsymbol x_{N(a)}} p_a(\boldsymbol x_{N(a)}) &= 1,\\ \sum_{\boldsymbol x_{N(a) \backslash n}} p_a(\boldsymbol x_{N(a)}) &= p_n(x_n) \quad \forall n \in [1, N]. \end{align}

Упражнение.

Найдите разложение Бете для функции из нашего примера, предполагая, что наша функция задаёт вероятностное распределение.

Напомню, что все переменные \( x_1, \ldots, x_n \) принимают значения лишь в конечном алфавите \( \Sigma \). Это означает, что любая функция от одной переменной \( x_i \) может быть представлена в виде вектора длины \( |\Sigma| \). Это значит, что разложение Бете может быть параметризовано с помощью фиксированного набора параметров.

Теорема.

Если граф сомножителей, соответствующий \(p(\boldsymbol x)\) имеет циклы, то шаги алгоритма суммы-произведения эквивалентны координатному спуску для дивергенции Кульбака-Лейблера (KL-дивергенции) между \( p(\boldsymbol x)\) и распределением типа Бете.

Доказательство и точную формулировку можно найти в приложении A диссертации Delbert Dueck.

Теперь мы можем перейти к приложениям.

Моделирование систем с помощью графов сомножителей

Affinity propagation.

Affinity propagation (распространение сродства) это алгоритм кластеризации. Его цель состоит в том, чтобы максимизировать функционал net similarity.

Сформулируем задачу следующим образом: задана последовательность наблюдений \( X_1, \ldots, X_n \), а также матрица схожести (это как расстояния только наоборот) \( s(i, j) \). Часто бывает так, что мы не можем наблюдать сами переменные \( X_1, \ldots, X_n \), но мы всегда можем наблюдать матрицу схожести.

В задаче метрической кластеризации, то есть когда \( X_j \in \mathbb R^d\), можно выбрать схожесть по правилу \( s(i, j) = - d(X_i, X_j) \) для \( i \neq j\), однако для \( i = j\) необходимо, чтобы \( s(X_i, X_j) \neq 0 \). Фактически, существует несколько возможных стратегий для того, чтобы определить матрицу схожести вершины с самой собой: \[ s(i, i) = -\lambda; \quad \text{or} \quad s(i, i) = \mathrm{Median}_{j}(s(i, j)) \]
Среди наблюдений \( c_i \in \{X_1, \ldots, X_n \} \) мы выбираем так называемые экземпляры \( c_1, c_2, \ldots, c_n \) так, чтобы сумма схожестей \[ \sum_{i = 1}^{n} s(i, c_i) \] была максимальной. При этом есть одно ограничение, которое препятствует «жадному» назначению экземпляров. \[ (c_i = k) \, \Rightarrow \, (c_k = k), \] Необходимо, чтобы назначенные экземпляры были корректными в том смысле, что если точка \( X_k \) является экземпляром для точки \( X_i \), то она должна являться экземпляром для самой себя. Это требование можно переформулировать следующим образом: множество точек разбито на непересекающиеся подмножества точек, где каждое подмножество имеет свой собственный единственный экземпляр.

Reference: Brendan J. Frey and Delbert Dueck, Clustering by Passing Messages Between Data Points, Science Feb. 2007
Reference: Brendan J. Frey and Delbert Dueck, “Clustering by Passing Messages Between Data Points”, Science Feb. 2007

Такая оптимизационная задача может быть проинтерпретирована как проблема размещения объектов (facility Location Problem), и она является NP-трудной.

Оказывается, что целевая функция вместе с ограничением корректности может быть преобразована в некоторую другую функцию без ограничений. А именно, мы рассматриваем \[ F(\boldsymbol c, \boldsymbol s) = \prod_{i=1}^{N} e^{s(i, c_i)} \prod_{k=1}^{N} f_k (\underbrace{c_1, \ldots, c_N}_{\boldsymbol c}) \] Второй член содержит ограничение корректности, определённое следующим образом: \[ f_k(\boldsymbol c) = \begin{cases} 0, & c_k \neq k, \, \exists i \colon c_i = k\\ 1, & \mathrm{otherwise} \end{cases} \]

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

Фактически, анализ механизма передачи сообщений для функции \( F(\boldsymbol c, \boldsymbol s) \) является довольно запутанным, необходимо рассматривать несколько случаев для корректных и некорректных конфигураций. Передача сообщений между \( c_i\) и \( f_j \) является наиболее важной, в то время как сообщения между \( s(i, c_i) \) и \( c_i \) можно легко исключить. Авторы алгоритма Affinity Propagation доказывают, что процедуру передачи сообщений, после некоторого преобразования обозначений, можно представить в следующем виде:

Affinity Propagation

  1. Set \( a(i, k) \) to zero (availabilities).
  2. Repeat until convergence:
    • \( \forall i, k \colon r(i, k) = s(i, k) - \max_{k’ \neq k} [s(i, k’) + a(i, k’)] \)
    • \( \forall i, k \colon a(i, k) = \begin{cases} \sum_{i’ \neq i} [r(i’, k)]_{+}, & k = i\\ [r(k,k) + [r(i’, k)]_{+}]_{-}, & k \neq i \end{cases} \)
  3. Output: cluster assignments.
    • \( \hat{\boldsymbol c} = (\hat c_{1}, \ldots, \hat c_{N}) \), where
    • \( \hat c_i = \arg\max_{k} [a(i,k) + r(i,k)] \)

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

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

Тем не менее, утверждается, что этот алгоритм хорошо себя ведёт на практике.

Заинтересованный читатель может обратиться к официальному FAQ, которое включает некоторые вероятностные интерпретации affinity propagation: FAQ Page.