Многомерное шкалирование

---

Введение

Многомерное шкалирование - это методы визуализации многомерных данных. Пусть имеется $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 Чем меньше энтропия, тем более "ажурные" получаются кластеры.

Самостоятельно поиграться с параметрами обучения можно здесь.

---