Многомерное шкалирование
Введение
Многомерное шкалирование - это методы визуализации многомерных данных. Пусть имеется $n$-мерное пространство признаков, в котором каждый объект описывается координатами $\mathbf{x}=\{x^1,...,x^n\}$ (номера признаков помечаем верхними индексами). В пространстве находится $N$ объектов с координатами $\mathbf{x}_1,...,\mathbf{x}_N$ (нижний индекс - номер объекта). Необходимо отобразить это пространство на $m$-мерное пространство в координатах $\mathbf{y}=\{y^1,...,y^m\}$, где $m$ существенно меньше $n$ (обычно $m=2$ или $m=3$). При этом исходная структура данных должна быть сохранена. В частности, объекты близкие в пространстве $\mathbf{x}$, должны быть близки и в пространстве $\mathbf{y}$.
Минимизация расстояний
В качестве меры близости объектов с координатами $\mathbf{a}$ и $\mathbf{b}$ будем использовать евклидово расстояние: $$ d(\mathbf{a},\mathbf{b}) = \|\mathbf{a}-\mathbf{b} \|= \sqrt{(\mathbf{a}-\mathbf{b})^2} = \sqrt{(a^1-b^1)^2+...+(a^n-b^n)^2} $$ и аналогично в $m$-мерном пространстве. Напомним, что евклидово расстояние подразумевает, что все признаки объектов имеют равную дисперсию, независимы и одинаково важны при классификации (дают равный вклад в расстояние между объектами).
Обозначим расстояния между объектами в обоих пространствах следующим образом (они симметричны по индексам $i,j$):
$$
D_{ij}=\bigr\|\mathbf{x}_i-\mathbf{x}_j \bigr\|,~~~~~~~~~~~d_{ij}=\bigr\|\mathbf{y}_i-\mathbf{y}_j \bigr\|.
$$
Критерий совпадения величин $D_{ij}$ и $d_{ij}$ задаётся той или иной функцией ошибок $Cost$.
Она зависит от положения объектов $\mathbf{y}_1$,...$\mathbf{y}_N$ в $m$-мерном пространстве (в $n$-мерном пространстве они заданы).
Минимизация функции ошибок проводится градиентным методом. Вычислим производную $Cost$ по $\mathbf{y}_p$.
От координат $p$-того объекта зависят расстояния $d_{p1},...,d_{pN}$ (в списке нет $d_{pp}$), поэтому:
$$
\frac{\partial Cost}{\partial \mathbf{y}_p} = \sum^{i \neq p}_{i}
\frac{\partial Cost}{\partial d_{pi}}\,\frac{\partial d_{pi}}{\partial \mathbf{y}_p}
=\sum^{i\neq p}_{i}
\frac{\partial Cost}{\partial d_{pi}}\,\frac{\mathbf{y}_p-\mathbf{y}_i}{d_{pi}}.
$$
В градиентном методе минимизации функции $Cost$ необходимо задать в $m$-мерном пространстве координаты объектов $\mathbf{y}_1,...,\mathbf{y}_N$ случайным образом.
Затем на каждой итерации их подправлять:
$$
\mathbf{y}_p(t+1) = \mathbf{y}_p(t) - g\,
\sum^{i\neq p}_{i} \,\frac{1}{d_{pi}}\,\frac{\partial C}{\partial d_{pi}}\, \big(\mathbf{y}_p(t)-\mathbf{y}_i(t)\big),
$$
где $g$ - скорость обучения.
Множитель при $\mathbf{y}_p(t)-\mathbf{y}_i(t)$ обычно интерпретируется как некая "сила", действующая на $p$-тую частицу
со стороны $i$-той частицы. Если производная $\partial C/\partial d_{pi}$ положительна - $p$-я частица притягивается к $i$-й,
а если отрицательна - отталкивается (см. рисунок справа для $\partial C/\partial d_{pi} < 0$).
Sammon mapping
Простейший выбор функции ошибок был предложен J.W. Sammon в 1969: $$ Cost = \sum^{k < q}_{k,q} \frac{\bigr(d_{kq} - D_{kq}\bigr)^2}{D_{kq}}. $$ Функция ошибок $Cost$ является суммой квадратичных отклонений расстояний между объектами в обоих пространствах. Такая сумма должна быть минимальной. При этом вес больших расстояний (см. знаменатель) меньше, чем вес маленьких расстояний. Знак $k < q$ над суммой подразумевает, что первый индекс суммирования всегда меньше второго: $$ \sum^{k < q}_{k,q} S_{kq}\equiv \frac{1}{2}\sum^{k \neq q}_{k,q} S_{kq} \equiv \sum^{N}_{k=1}\sum^{N}_{q=k+1} S_{qk} = S_{12}+S_{13}+...+S_{1N} + S_{23}+S_{24}+...+S_{2N} + ...+ S_{N-1,N}. $$ Несложно видеть (скажем для $p=1$), что $$ \frac{\partial C}{\partial d_{pi}} = 2\,\frac{d_{pi}-D_{pi}}{D_{pi}}. $$ Поэтому итерации градиентного метода выглядят следующим образом (множитель 2 включаем в $g$): $$ \mathbf{y}_p(t+1) = \mathbf{y}_p(t) - g\, \sum^{i\neq p}_{i} \frac{d_{pi}-D_{pi}}{D_{pi}d_{pi}}\, \big(\mathbf{y}_p(t)-\mathbf{y}_i(t)\big). $$ Если $d_{pi} > D_{pi}$ объекты притягиваются (уменьшая тем самым расстояние $d_{pi}$). В противном случае, они отталкиваются.
Стохастическое вложение соседей (SNE)
Минимизация разности расстояний в двух различных пространствах не учитывает скученности объектов.
Так, на рисунке справа можно предположить, что важнее сохранить
расстояние между объектами $A$ и $B$, чем между $B$ и $C$ (объект $C$, удалённый от сгущения может быть
ошибкой или случайным выбросом). Поэтому стоит использовать меру близости, учитывающую
густоту или разреженность объектов в пространстве.
В 2002 Hinton и Roweis предложили метод "Стохастического вложения соседей" (Stochastic Neighbor Embedding, SNE), в котором учитывается плотность распределения объектов. Позднее, в 2008 Maaten и Hinton его несколько модифицировали, назвав "t-SNE методом".
Идея SNE-метода состоит в записи распределения вероятностей объектов в $n$-мерном и $m$-мерном пространствах. Координаты $\mathbf{y}=\{y^1,...,y^m\}$ объектов в $m$-мерном пространстве подбираются таким образом, чтобы эти распределения были максимально похожи.
Пусть вероятность того, что в точке $\mathbf{x}$ окажется объект, соседний к объекту $\mathbf{x}_k$ (и отличный от него) равна: $$ P_k(\mathbf{x}) = \frac{1}{Z_k}\,e^{-\alpha_k\,\bigr(\mathbf{x}-\mathbf{x}_k\bigr)^2},~~~~~~~~Z_k=\sum^{i\neq k}_{i} e^{-\alpha_k\,\bigr(\mathbf{x}_i-\mathbf{x}_k\bigr)^2}, ~~~~~~~~~~~~~~~~P_k(\mathbf{x}_k)=0, $$ где $\alpha_k$ - некоторые константы, а нормировочный множитель $Z_k$ означает, что вероятности нормируются на единицу по всем объектам: $$ \sum^{i\neq k}_{i} P_k(\mathbf{x}_i) = 1. $$ Величина $P_k(\mathbf{x}_i)$ интерпретируется как условная вероятность обнаружить объект в $\mathbf{x}_i$, если есть объект в $\mathbf{x}_k$. Чем дальше точка $\mathbf{x}$ от $\mathbf{x}_k$, тем менее вероятнее там встретить соседа $k$-того объекта. Отметим, что в общем случае $P_i(\mathbf{x}_j)\neq P_j(\mathbf{x}_i)$.
Выбор параметров вероятностей SNE
Параметры $\alpha_k$ уникальны для каждого объекта. Они выбираются тем большими, чем выше плотность объектов вокруг $\mathbf{x}_k$. Для этого вычисляется энтропия распределения: $$ H_k = -\sum^N_{i=1} P_k(\mathbf{x}_i)\ln P_k(\mathbf{x}_i) $$ и параметры $\alpha_k$ подбираются так, чтобы $H_k$ была постоянна для всех $k$.
Чтобы пояснить такой выбор $\alpha_k$, поставим численный эксперимент. Сгенерим в n-мерном пространстве облако объектов ($N=1001$) с нулевым средним положением и нормально распределённым расстоянием от центра координат с волатильностями $\sigma=1,2,3,4,5$. Чем $\sigma$ больше, тем менее плотное распределение частиц возле начала координат. Зададим $\alpha = 1/(2\sigma^2)$ и вычислим энтропию относительно точки $\mathbf{x}=\mathbf{0}$ при разных размерностях пространства:
σ α H(n=10) H(n=100) 1 0.500 6.808 6.808 2 0.125 6.822 6.808 3 0.056 6.812 6.810 4 0.031 6.809 6.818 5 0.020 6.804 6.802
Как видно, она постоянна, поэтому выбор параметра $\alpha$ при фиксированной энтропии $H_k=const$, действительно, отражает плотность распределения частиц.
Отметим, что при энтропию быстрее вычислять при помощи эквивалентной формулы: $$ H_k = \ln Z_k+ \frac{\alpha_k}{Z_k}\,\sum^N_{i=1} e^{-\alpha_k\,\bigr(\mathbf{x}-\mathbf{x}_k\bigr)^2}~ (\mathbf{x}_i-\mathbf{x}_k)^2 $$
Функция ошибки SNE
В качестве распределения вероятностей в $m$-мерном пространстве выбирается $$ Q_k(\mathbf{x}) = \frac{1}{Z_k}\,e^{-\bigr(\mathbf{y}-\mathbf{y}_k\bigr)^2},~~~~~~~~Z_k=\sum^{i\neq k}_{i} e^{-\bigr(\mathbf{y}_i-\mathbf{y}_k\bigr)^2}, ~~~~~~~~~~~~~~~~Q_k(\mathbf{y}_k)=0. $$ В отличии от $n$-мерного пространства вероятности имеют одинаковый множитель при расстояниях в показателях экспонент. Это подразумевает, что в пространстве $m=2$ или $m=3$ объекты окажутся распределёнными с примерно равной плотностью (да и подбирать дополнительные параметры проблематично).
Выберем функцию ошибки (=функцию стоимости), следуя теории информации, в виде расстояния Кульбака-Лейблера между двумя распределениями вероятностей: $$ Cost = \sum^{i\neq j}_{i,j} P_i(\mathbf{x}_j)\,\ln \frac{P_i(\mathbf{x}_j)}{Q_i(\mathbf{y}_j)},~~~~~~~~~~~~\sum^{j\neq i}_iQ_i(\mathbf{y}_j) = 1. $$ В силу нормировочного условия на $Q_i(\mathbf{y}_j)$, эта величина достигает минимального значения (нуля) когда все $Q_i(\mathbf{y}_j)=P_i(\mathbf{x}_j)$.
Градиент SNE
Вычислим производную функции ошибок $Cost$ по расстояниям $d_{pq}$. При взятии производной от симметричного тензора необходимо учитывать формулу: $$ \frac{\partial d_{ij}}{\partial d_{pq}} = \delta_{ip}\delta_{jq}+\delta_{iq}\delta_{jp}, $$ где $\delta_{ip}$ - символ Кронекера, равный 1 при $i=p$ и 0 при $i\neq p$. Действительно, в соотвествии с этой формулой например, $\partial d_{12}/\partial d_{12}=\partial d_{12}/\partial d_{21}=1$, а, скажем, $\partial d_{12}/\partial d_{13}=0$. В связи с этим, имеем: $$ \frac{\partial Cost}{\partial d_{pq}} = -\sum^{i\neq j}_{i,j} \frac{P_i(\mathbf{x}_j)}{ Q_i(\mathbf{y}_j)}\,\frac{\partial Q_i(\mathbf{y}_j)}{\partial d_{pq}} = \sum^{i\neq j}_{i,j} P_i(\mathbf{x}_j)\, \Bigr[ 2\, d_{ij}\,(\delta_{ip}\delta_{jq}+\delta_{iq}\delta_{jp}) +\frac{1}{Z_i}\frac{\partial Z_i}{\partial d_{pq}} \Bigr] ~=~ 2d_{pq}\,\bigr[P_p(\mathbf{x}_q)+P_q(\mathbf{x}_p)\bigr] + \sum_i\frac{1}{Z_i}\,\frac{\partial Z_i}{\partial d_{pq}}, $$ где в последнем слагаемом сумма по $j$ положена единице в силу нормировочного условия на $P_i(\mathbf{x}_j)$. Производная нормировочного множителя равна (учитываем, что $d_{iq}\,\delta_{ip}\equiv d_{pq}\,\delta_{ip}$): $$ \frac{\partial Z_i}{\partial d_{pq}} = \sum^{j\neq i}_{j} \frac{\partial\, e^{-d^2_{ij}}}{\partial~ d_{pq}} = -2 \sum^{j\neq i}_{j} e^{-d^2_{ij}} ~ d_{ij}\,(\delta_{ip}\delta_{jq}+\delta_{iq}\delta_{jp}) = -2\,Z_i \,d_{pq} ~\Bigr(Q_i(\mathbf{y}_q) ~ \delta_{ip} + Q_i(\mathbf{y}_p) ~ \delta_{iq}\Bigr). $$ В результате "сила", действующая на $p$-й объект со стороны $q$-тых объектов, равна: $$ \frac{1}{d_{pq}}\,\frac{\partial Cost}{\partial d_{pq}} = 2\,\Bigr(P_p(\mathbf{x}_q)-Q_p(\mathbf{y}_q)+P_q(\mathbf{x}_p)-P_q(\mathbf{y}_p)\Bigr). $$ Соответственно, итерации по уточнению положения объектов в $m$-мерном пространстве имеют вид: $$ \mathbf{y}_p(t+1) = \mathbf{y}_p(t) - g\, \sum^{q\neq p}_{q} \Bigr(P_p(\mathbf{x}_q)-Q_p(\mathbf{y}_q)+P_q(\mathbf{x}_p)-P_q(\mathbf{y}_p)\Bigr)\, \big(\mathbf{y}_p(t)-\mathbf{y}_q(t)\big). $$
Метод t-SNE
В методе t-SNE, по сравнению с SNE, сделаны два изменения. Во-первых, от несимметричных "условных" вероятностей в $n$-мерном пространстве переходят к симметричным "совместным" вероятностям: $$ P_{ij}=P(\mathbf{x}_i, \mathbf{x}_j) = \frac{P_i(\mathbf{x}_j)+P_j(\mathbf{x}_i)}{2N} $$ Число объектов $N$ в знаменателе подразумевает следующую нормировку вероятностей: $$ \sum^{i\neq j}_{i,j} P_{ij} = 1. $$
Для вероятностей в $m$-мерном пространстве берётся распределение Коши или эквивалентно t-распределение Стьюдента с одной степенью свободы (отсюда приставка "t-" в названии метода): $$ Q_{ij}=Q(\mathbf{y}_i,\mathbf{y}_j) = \frac{1/Z}{1+(\mathbf{y}_i-\mathbf{y}_j)^2},~~~~~~~~~~~~Z = \sum^{i\neq j}_{i,j}\frac{1}{1+(\mathbf{y}_i-\mathbf{y}_j)^2}, ~~~~~~~~~~~~~~\sum^{i\neq j}_{i,j} Q_{ij} = 1. $$ Вероятности $Q_{ij}$, как и $P_{ij}$, симметричны: $Q_{ij}=Q_{ji}$.
Для функции ошибок по-прежнему выбирается расстояние Кульбака-Лейблера: $$ Cost = \sum^{i\neq j}_{i,j} P_{ij}\,\ln \frac{P_{ij}}{Q_{ij}}. $$
Градиент t-SNE
Вычисление градиента в t-SNE методе производится аналогично: $$ \frac{\partial Cost}{\partial d_{pq}} = -\sum^{i\neq j}_{i,j} \frac{P_{ij}}{ Q_{ij}}\,\frac{\partial Q_{ij} }{\partial d_{pq}} = \sum^{i\neq j}_{i,j} P_{ij} \, \Bigr[ \frac{2\, d_{ij}\,(\delta_{ip}\delta_{jq}+\delta_{iq}\delta_{jp})}{1+d^2_{ij}} +\frac{1}{Z}\frac{\partial Z}{\partial d_{pq}} \Bigr] = \frac{4\,d_{pq}\,P_{pq}}{1+d^2_{pq}} + \frac{1}{Z}\,\frac{\partial Z}{\partial d_{pq}}, $$ где в последнем слагаемом учтено нормировочное условие на $P_{ij}$. Ещё одно дифференцирование: $$ \frac{\partial Z}{\partial d_{pq}} = \sum^{j\neq i}_{i,j} \frac{\partial\, (1+d^2_{ij})^{-1}}{\partial~ d_{pq}} = -2 \sum^{j\neq i}_{i,j} \frac{ d_{ij}\,(\delta_{ip}\delta_{jq}+\delta_{iq}\delta_{jp})}{(1+d^2_{ij})^2} = -\frac{4\, d_{pq}}{(1+d^2_{pq})^2} $$ приводит к финальному, достаточно простому выражению для "силы": $$ \frac{1}{d_{pq}}\,\frac{\partial Cost}{\partial d_{pq}} = 4\,\frac{P_{pq}-Q_{pq}}{1+(\mathbf{y}_p-\mathbf{y}_q)^2}. $$ Итерации по уточнению положения объектов в $m$-мерном пространстве t-SNE метода имеют вид: $$ \mathbf{y}_p(t+1) = \mathbf{y}_p(t) - g\, \sum^{q\neq p}_{q} \frac{P_{pq}-Q_{pq}}{1+(\mathbf{y}_p-\mathbf{y}_q)^2}~\, \big(\mathbf{y}_p(t)-\mathbf{y}_q(t)\big) + m \cdot\bigr(\mathbf{y}_p(t)-\mathbf{y}_p(t-1)\bigr), $$ где последнее слагаемое является инерционным (метод тяжёлого шарика, $0 \le m\le 1$).
Параметры обучения
Авторы (Maaten и Hinton) рекомендуют в 2-мерном пространстве задавать случайное положение объектов с нормальным распределением с нулевым средним и маленькой $10^{-4}$ дисперсией (волатильность $10^{-2}$).
Мои эксперименты показывают, что в итерационную формулу стоит добавить ещё один параметр $\alpha$, а авторы включают также параметр $\beta$: $$ \mathbf{y}_p(t+1) = \mathbf{y}_p(t) - g\, \sum^{q\neq p}_{q} \frac{\beta\,P_{pq}-Q_{pq}}{1+\alpha\,(\mathbf{y}_p-\mathbf{y}_q)^2}~\, \big(\mathbf{y}_p(t)-\mathbf{y}_q(t)\big) + m \cdot\bigr(\mathbf{y}_p(t)-\mathbf{y}_p(t-1)\bigr), $$ На начальных стадиях обучения $\alpha$ выбирается малым ($\alpha=0.1-0.2$), чтобы объектам легче было "проталкиваться" к своим соседям через облако остальных объектов. Постепенно, этот параметр повышается до единицы. Параметр скорости обучения в методе t-SNE у меня достаточно высок ($g=1000-5000$), хотя авторы рекомендуют значение 100 с адаптивным увеличением по алгоритму Jacobs (1988). Параметр инерции, обычно, равен $m=0.5$. Авторы рекомендуют его повысить до $0.8$ после 250 итераций.
Авторы также дают ряд эзотерических советов. Например, умножать $P_{ij}$ на $\beta=4$ в начальных стадиях (50 итераций) оптимизации (early exaggeration), затем возвращать его к теоретическому значению $\beta=1$. Ещё один трюк, названный ранней компрессией (early compression) добавляет L2-штрафы в функцию ошибок, которые пропорциональны сумме квадратов расстояний объектов от начала координат. В формуле итерций это даёт дополнительный член $\gamma\,y_p(t)$. При этом предполагается, что когда расстояния между объектами невелики, кластерам легче проходить друг через друга (выше у меня для этого служит параметр $\alpha$)
Энтропия (в терминах натурального логарифма) обычно выбирается равной порядка 3.7 Чем меньше энтропия, тем более "ажурные" получаются кластеры.
Самостоятельно поиграться с параметрами обучения можно здесь.