Pull to refresh

Препарируем t-SNE

Reading time 10 min
Views 79K
Работая над статьей «Глубокое обучение на R...», я несколько раз встречал упоминание t-SNE — загадочной техники нелинейного снижения размерности и визуализации многомерных переменных (например, здесь), был заинтригован и решил разобраться во всем в деталях. t-SNE это t-distributed stochastic neighbor embedding. Русский вариант с «внедрением соседей» в некоторой мере звучит нелепо, поэтому дальше буду использовать английский акроним.


Алгоритм t-SNE, который также относят к методам множественного обучения признаков, был опубликован в 2008 году (ссылка 1 в конце статьи) голландским исследователем Лоуренсом ван дер Маатеном (сейчас работает в Facebook AI Research) и чародеем нейронных сетей Джеффри Хинтоном. Классический SNE был предложен Хинтоном и Ровейсом в 2002 (ссылка 2). В статье 2008 года описывается несколько «трюков», которые позволили упростить процесс поиска глобальных минимумов, и повысить качество визуализации. Одним из них стала замена нормального распределения на распределение Стьюдента для данных низкой размерности. Кроме того, была сделана удачная реализация алгоритма (в статье есть ссылка на MatLab), которая потом портировалась в другие популярные среды.

Немного математики


Начнем с «классического» SNE и сформулируем задачу. У нас есть набор данных с точками, описываемыми многомерной переменной с размерностью пространства существенно больше трех. Необходимо получить новую переменную, существующую в двумерном или трехмерном пространстве, которая бы в максимальной степени сохраняла структуру и закономерности в исходных данных. SNE начинается с преобразования многомерной евклидовой дистанции между точками в условные вероятности, отражающие сходство точек. Математически это выглядит следующим образом (формула 1):



Эта формула показывает, насколько точка Xj близка к точке Xi при гауссовом распределении вокруг Xi с заданным отклонением σ. Сигма будет различной для каждой точки. Она выбирается так, чтобы точки в областях с большей плотностью имели меньшую дисперсию. Для этого используется оценка перплексии:



Где H(Pi) – энтропия Шеннона в битах (формула 2):



В данном случае перплексия может быть интерпретирована как сглаженная оценка эффективного количества «соседей» для точки Xi. Она задается в качестве параметра метода. Авторы рекомендуют использовать значение в интервале от 5 до 50. Сигма определяется для каждой пары Xi и Xj при помощи алгоритма бинарного поиска.

Для двумерных или трехмерных «коллег» пары Xi и Xj, назовем их для ясности Yi и Yj, не представляет труда оценить условную вероятность, используя ту же формулу 1. Стандартное отклонение предлагается установить в 1/√2:



Если точки отображения Yi и Yj корректно моделируют сходство между исходными точками высокой размерности Xi и Xj, то соответствующие условные вероятности pj|i и qj|i будут эквивалентны. В качестве очевидной оценки качества, с которым qj|i отражает pj|i, используется дивергенция или расстояние Кульбака-Лейблера. SNE минимизирует сумму таких расстояний для всех точек отображения при помощи градиентного спуска. Функция потерь для данного метода будет определяться формулой 3:



При этом градиент выглядит на удивление просто:



Авторы предлагают следующую физическую аналогию для процесса оптимизации: Все точки отображения соединены пружинами. Жесткость пружины, соединяющей точки i и j зависит от разности между сходством двух точек в многомерном пространстве и двух точек в пространстве отображения. В этой аналогии, градиент — это результирующая сила, действующая на точку в пространстве отображения. Если систему «отпустить», через какое-то время она придет в равновесие, это и будет искомое распределение. Алгоритмически, поиск равновесия предлагается делать с учетом моментов:



где η – параметр, определяющий скорость обучения (длину шага), а α – коэффициент инерции. Использование классического SNE позволяет получить неплохие результаты, но может быть связано с трудностями в оптимизации функции потерь и проблемой скученности (в оригинале – crowding problem). t-SNE если и не решает эти проблемы совсем, то существенно облегчает. Функция потерь t-SNE имеет два принципиальных отличая. Во-первых, у t-SNE симметричная форма сходства в многомерном пространстве и более простой вариант градиента. Во-вторых, вместо гауссова распределения для точек из пространства отображения используется t-распределение (Стьюдента), «тяжелые» хвосты которого облегчают оптимизацию и решают проблему скученности.

В качестве альтернативы минимизации суммы дивергенций Кульбака-Лейблера между условными вероятностями pi|j и qi|j предлагается минимизировать одиночную дивергенцию между совместной вероятностью P в многомерном пространстве и совместной вероятностью Q в пространстве отображения:



где pii и qii = 0, pij = pji, qij = qji для любых i и j, а pij определяется по формуле:



где n — количество точек в наборе данных. Градиент для симметричного SNE получается существенно проще, чем для классического:



Проблема скученности заключается в том, что расстояние между двумя точками в пространстве отображения, соответствующими двум среднеудаленным точкам в многомерном пространстве, должно быть существенно больше, нежели расстояние, которое позволяет получить гауссово распределение. Проблему решают хвосты Стьюдента. В t-SNE используется t-распределение с одной степенью свободы. Совместная вероятность для пространства отображения в этом случае будет определяться формулой 4:



А соответствующий градиент – выражением 5:



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

Алгоритм


Алгоритм t-SNE в упрощенном виде можно представить следующим псевдокодом:

Data: набор данных X = {x1, x2, …, xn},
параметр функции потерь: перплексия Perp,
Параметры оптимизации: количество итераций T, скорость обучения η, момент α(t).
Result: представление данных Y(T) = {y1, y2, …, yn} (в 2D или 3D).
begin
	вычислить попарное сходство pj|i c перплексией Perp (используя формулу 1)
	установить pij = (pj|i + pi|j)/2n
	инициализировать Y(0) = {y1, y2, …, yn} точками нормального распределения (mean=0, sd=1e-4)
	for t = 1 to T do
 		вычислить сходство точек в пространстве отображения qij (по формуле 4)
		вычислить градиент δCost/δy (по формуле 5)
		установить Y(t) = Y(t-1) + ηδCost/δy + α(t)(Y(t-1) - Y(t-2))
	end
end


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

Второй трюк менее очевидный — «раннее гиперусиление» (в оригинале – «early exaggeration»). Заключается он в умножении в начале оптимизации всех pij на некоторое целое число, например на 4. Смысл в том, чтобы для больших pij получить бόльшие qij. Это позволит для кластеров в исходных данных получить плотные и широко разнесенные кластеры в пространстве отображения.

Реализация на R


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

Для экспериментов предлагается использовать следующие параметры:
Количество итераций: 1000;
Перплексия: 40;
Раннее гиперусиление: х4 для первых 50 итераций;
Момент α: 0.5 для первых 250 итераций, 0.8 для оставшихся 750;
Скорость обучения η: первоначально 100 (в оригинальном коде — 500) и изменяется после каждой итерации в соответствие со схемой, описанной Робертом Джекобсом в 1988 (ссылка 3).

Алгоритм можно разбить на три логические части:
  1. Вычисление перплексии и гауссового ядра для вектора текущих значений сигма (вернее gamma, я использовал переменную, обратную квадрату сигма)
  2. Вычисление попарного сходства pij в многомерном пространстве при помощи бинарного поиска для заданной перплексии
  3. Вычисление попарного сходства для пространства отображения, функции потерь и градиента


Проиллюстрирую только самые важные конструкции. Прежде всего, необходимо рассчитать матрицу квадратов евклидовой дистанции для исходного набора данных. Для этого лучше использовать функцию rdist() из пакета fields. Она работает в разы быстрее, чем популярная dist() из базового пакета stats. Возможно, это по тому, что она оптимизирована именно для евклидовой дистанции. Далее нужно инициализировать пространство отображения. Это можно сделать при помощи rnorm() с нулевым мат. ожиданием и стандартным отклонением 1е-4.

Расчет логарифма перплексии и гауссового ядра выглядит следующим образом:

P <- exp(-Di * gamma) # kernel calc    
sumP <- sum(P)
PP <- log(sumP) + gamma * sum(Di * P) / sumP

Последнее выражение получается, если формулу (1) для pj|i подставить формулу (2) для перплексии. Пусть n будет количеством точек в исходных данных. Di – вектор длины n-1, с квадратами попарных дистанций для строки исходного данных. Получается он при помощи конструкции Di < — D[i, -i], где D – матрица квадратов дистанций n x n, и если i = 1, то Di – это первая строка матрицы D без первого столбца. Ядро, P – вектор, длины n-1, а логарифм перплексии – вещественное число (за счет стоящей перед поэлементным произведением суммы). В оригинальном коде используется натуральный логарифм, я решил поступить также.

Для поиска gamma существуют следующие ограничения – количество попыток (делений) — не более 50, порог разницы между логарифмами заданной и вычисленной перплексии 1e-5. Одна итерация бинарного поиска выглядит следующим образом:

            if (perpDiff > 0){
                   gammamin = gamma[i]
                   if (is.infinite(gammamax)) gamma[i] = gamma[i] * 2
                   else gamma[i] = (gamma[i] + gammamax)/2
            } else{
                   gammamax = gamma[i]
                   if (is.infinite(gammamin))  gamma[i] = gamma[i]/ 2
                   else gamma[i] = ( gamma[i] + gammamin) / 2
            }

gammamin и gammamax изначально устанавливаются в – и + Inf.

Функция потерь считается как сумма по строкам поэлементного произведения двух матриц n x n:

cost =  sum(apply(P * log(P/Q),1,sum)) 

P – симметричная версия совместной вероятности для многомерного пространства, рассчитывается как P = .5 * (P + t(P)) из условной попарной вероятности, t – оператор транспонирования матриц. Q – совместная вероятность для пространства отображения:

            num = 1/(1 + (rdist(Y))^2)
            diag(num)=0 # Set diagonal to zero
            Q = num / sum(num) # Normalize to get probabilities

Y – инициализированная матрица точек n x 2 для двумерного пространства отображения. Q получается матрицей n x n за счет использования rdist().

Градиент рассчитывается в векторизованной форме сразу для всех Yi:

L = (P - Q) * num;
grads = 4 * (diag(apply(L, 1, sum)) - L) %*% Y

Здесь grads это матрица, размером n x 2, которая образуется за счет матричного произведения. Конструкция в круглых скобках — матрица размера n x n, на главной диагонали которой построчные суммы матрицы L, а все остальные элементы взяты со знаком минус. Матричное произведение позволяет умножить на разность Yi и Yj и сразу посчитать суммы по j (формула 4). Правда градиент получается обратного знака.

Для поиска оптимальных значений Y используется что-то похожее на эвристику Роберта Джекобса «дельта-бар-дельта». Смысл её в том чтобы адаптивно изменять шаг обучения. Если градиент не меняет свой знак на очередной итерации, шаг увеличивается линейно (к коэффициенту усиления добавляется 0.2), если знак меняется — то уменьшается экспоненциально (коэффициент умножается на 0.8). С учетом того, что градиент у нас обратного знака:

gains = (gains + .2) * abs(sign(grads) != sign(incs)) 
     + gains * .8 * abs(sign(grads) == sign(incs))		
gains[gains < min_gain] = min_gain
          
incs = momentum * incs - gains  * epsilon * grads
 Y = Y + incs


Эксперименты


Я не мог не попробовать повторить эксперимент, описанный в статье (ссылка 1). Идея этого эксперимента – в визуализации структуры кластеров для набора изображений MNIST. Одна из частей этого набора содержит 60000 изображений рукописных цифр от 0 до 9, 28х28 пикселей. Набор данных можно загрузить с сайта гуру нейронных сетей, Яна Лекуна.

С загрузкой в R придется повозиться, но помогает функция readBin(), которая позволяет читать (и записывать) бинарные данные. 20 случайных цифр:



Как и в оригинальном эксперименте я использовал не весь набор, а только 6000 изображений, выбранных случайным образом при помощи функции sample(). 1000 итераций заняла чуть менее полутора часов на x64 c Intel Core i7 2.4 GHz. Результат на картинке:



Цвета задаются при помощи меток, MNIST – это размеченный набор данных. Результат можно существенно улучшить. Как советуют авторы, для исходного набора точек следует сделать анализ главных компонент (PCA) и выбрать первые 30:



Существует несколько реализаций PCA для R. Простейший вариант в две строки:

pca <- prcomp(X)
X <- pca$x[,1:30]

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

Техника t-SNE не была бы, наверное, так популярна, если бы не возможность визуализации распределенного представления слов. В этом эксперименте я использовал сравнительно небольшую модель (около 900 слов и 300-мерное пространство признаков), обученную при помощи word2vec.



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

Оказывается, с распределенным представлением слов, подготовленным на реальных данных, не всё так просто. Проблем две: как правило, очень большое количество кластеров и существенные различия по их размеру. Для визуализации можно сделать случайную выборку, однако результат будет, скорее всего, не репрезентативным. Авторы объясняют это на следующем примере: Если A, B и С — равноудаленные точки в многомерном пространстве, но между A и B есть еще точки, а между B и C нет, то алгоритм с большей вероятностью объединит A и B в один кластер нежели B и C.

Но бывает так, что важна именно наглядность, а не полнота картины. Тогда вместо случайной выборки точек можно попробовать отобрать интересные кластеры. Кластеризация встроена в оригинальный word2vec, кроме того можно использовать k-means. Возможно, этот вариант покажется не совсем «честным», но, на мой взгляд, имеет право на существование. Выборка данных для распределенного представления слов из пятнадцати кластеров выглядит уже не такой запутанной:



Для больших наборов данных есть специальная реализация t-SNE, в которой для вычисления сходства в многомерном пространстве используются деревья ближайших соседей и метод случайного блуждания (описана в статье 1). На практике можно использовать пакет Rtsne, обертку для C++ реализации Barnes-Hut-SNE. Rtsne позволяет обработать все 60000 изображений MNIST менее чем за час. Подробности можно посмотреть в публикации 6.

С Rtsne я еще не экспериментировал, а вот tsne (порт MatLab-овского кода для R) испытал. Визуализация получается гораздо хуже, чем мой вариант. Заглянув в архив пакета я обнаружил две вещи: Во-первых, при расчете гауссового ядра дистанция используется без квадарата, во-вторых, в выражении для перплексии стоит матричное произведение. В моей версии и в коде для MatLab оно поэлементное. Странно, что вообще получились видимые кластеры. Объясняю это «живучестью» алгоритма, но допускаю, что в чем-то я не до конца разобрался. Буду рад замечаниям и желаю всем впечатляющих многомерных визуализаций.

Формулы для этой статьи подготовлены при помощи онлайн редактора формул LATEX.

Ссылки:


  1. L.J.P. van der Maaten and G.E. Hinton. Visualizing High-Dimensional Data Using t-SNE. Journal of Machine Learning Research 9(Nov):2579-2605, 2008. PDF
  2. G.E. Hinton and S.T. Roweis. Stochastic Neighbor Embedding. In Advances in Neural Information. Processing Systems, volume 15, pages 833–840, Cambridge, MA, USA, 2002. The MIT Press. PDF
  3. R.A. Jacobs. Increased rates of convergence through learning rate adaptation. Neural Networks, 1: 295–307, 1988. PDF
  4. Эксперименты с t-SNE на Python: Алгоритм t-SNE. Иллюстрированный вводный курс (перевод оригинальной статьи на O’Reilly). datareview.info/article/algoritm-t-sne-illyustrirovannyiy-vvodnyiy-kurs
  5. Rtsne против tsne: www.codeproject.com/Tips/788739/Visualization-of-High-Dimensional-Data-using-t-SNE Здесь же доступен набор данных для экспериментов с распределенным представлением слов.
  6. L.J.P. van der Maaten. Accelerating t-SNE using Tree-Based Algorithms. Journal of Machine Learning Research 15(Oct):3221-3245, 2014. PDF
Tags:
Hubs:
+11
Comments 4
Comments Comments 4

Articles