t-SNE с нуля (ft. NumPy)

t-SNE с нуля (ft. NumPy)

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

Эта статья предназначена для читателей, которые заинтересованы в понимании t-SNE посредством перевода математики из оригинальной статьи — Лоренса ван дер Маатена и Джеффри Хинтона — в реализацию кода на python.[1] Я нахожу, что такого рода упражнения достаточно проливают свет на внутреннюю работу статистических алгоритмов / моделей и действительно проверяют ваше базовое понимание относительно этих алгоритмов / моделей. Как минимум, успешная реализация всегда приносит большое удовлетворение!

Эта статья будет доступна лицам с любым уровнем взаимодействия с t-SNE. Однако, стоить заметить, что эта статья не будет каким-либо полноценным руководством:

• Строго концептуальное введение и исследование t-SNE, поскольку существует множество других замечательных ресурсов, которые делают это; тем не менее, я сделаю всё возможное, чтобы связать математические уравнения с их интуитивными / концептуальными аналогами на каждом этапе реализации.
• Всестороннее обсуждение применений и плюсов /минусов t-SNE, а также прямые сравнения t-SNE с другими методами уменьшения размерности. Я, однако, кратко коснусь этих тем на протяжении всей этой статьи, но ни в коем случае не буду подробно их освещать.

Без лишних слов давайте начнём с краткого введения в t-SNE.

Краткое введение в t-SNE

Стохастическое вложение соседей с t-распределением (t-SNE) – это инструмент уменьшения размерности, который в основном используется в наборах данных с пространством пространственных объектов большой размерности и позволяет визуализировать данные или спроецировать их в пространство меньшей размерности (обычно 2-D). Это особенно полезно для визуализации нелинейно разделяемых данных, в которых линейные методы, такие как анализ главных компонент (PCA), потерпели бы неудачу. Обобщение линейных структур уменьшения размерности (таких как PCA) на нелинейные подходы (такие как t-SNE) также известно как многообразное обучение. Эти методы могут быть чрезвычайно полезны для визуализации и понимания базовой структуры многомерного нелинейного набора данных, а также могут быть полезны для распутывания и группировки воедино наблюдений, которые похожи в многомерном пространстве. Для получения дополнительной информации о t-SNE и других разнообразных методах обучения документация scikit-learn является отличным ресурсом. Кроме того, чтобы прочитать о некоторых интересных областях применения t-SNE, на странице Википедии освещаются некоторые из этих областей со ссылками на работу.

Давайте начнём с разбиения названия стохастического вложение соседей с t-распределением на его компоненты. t-SNE – это расширение стохастического вложения соседей (SNE), представленное 6 годами ранее в этой статье Джеффри Хинтоном и Сэмом Роуэйсом. Итак, давайте начнём с этого. Стохастическая часть названия происходит от того факта, что целевая функция не является выпуклой и, следовательно, при разных инициализациях могут получаться разные результаты. Вложение соседей подчёркивает природу алгоритма — оптимальное отображение точек в исходном многомерном пространстве в соответствующее низкоразмерное пространство при наилучшем сохранении структуры точек ”окрестности”. SNE состоит из следующих (упрощённых) шагов:

1. Получение матрицы подобия между точками в исходном пространстве: Вычислите условные вероятности для каждой точки данных j относительно каждой точки данных i. Эти условные вероятности вычисляются в исходном многомерном пространстве с использованием гауссова уравнения с центром в точке i и принимают следующую интерпретацию: вероятность того, что i выберет j в качестве своего соседа в исходном пространстве. Это создаёт матрицу, которая представляет сходства между точками.
2. Инициализация: Выберите случайные начальные точки в пространстве нижнего измерения (скажем, 2-D) для каждой точки данных в исходном пространстве и вычислите новые условные вероятности аналогично описанным выше в этом новом пространстве.
3. Отображение: Итеративно улучшайте точки в пространстве меньшей размерности до тех пор, пока расхождения Кульбака-Лейблера между всеми условными вероятностями не будут сведены к минимуму. По сути, мы сводим к минимуму различия в вероятностях между матрицами подобия двух пространств, чтобы обеспечить наилучшее сохранение сходства при отображении исходного набора данных в низкоразмерный датасетах.

t-SNE улучшает SNE двумя основными способами:

1. Это сводит к минимуму расхождения Кульбака-Лейблера между совместными вероятностями, а не условными. Авторы называют это “симметричным SNE”, поскольку их подход гарантирует, что совместные вероятности p_ij = p_ji. Это приводит к гораздо лучшему функционированию функции затрат, которую легче оптимизировать.
2. Он вычисляет сходства между точками, используя распределение Стьюдента-t с одной степенью свободы (также распределение Коши), а не гауссово в низкоразмерном пространстве (шаг 2 выше). Здесь мы можем видеть, откуда берётся буква “t” в t-SNE. Это улучшение помогает устранить “проблему скученности”, выделенную авторами, и ещё больше улучшить задачу оптимизации. Эту “проблему скученности” можно представить следующим образом: представьте, что у нас есть 10-мерное пространство, объема пространства, доступного в 2-мерном пространстве, будет недостаточно для точного захвата этих умеренно непохожих точек по сравнению с объемом пространства для соседних точек относительно объема пространства, доступного в 10-мерном пространстве. Проще говоря, просто представьте, что вы берёте трёхмерное пространство и проецируете его в 2-D. В трёхмерном пространстве будет гораздо больше общего пространства для моделирования сходств относительно 2-D проекции. Распределение Стьюдента-t помогает облегчить эту проблему, поскольку имеет более тяжелые хвосты, чем нормальное распределение.

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

Внедрение с Scratch

Давайте теперь перейдем к пониманию t-SNE посредством реализации оригинальной версии алгоритма, представленной в статье Лоуренса ван дер Маатена и Джеффри Хинтона. Сначала мы начнём с пошаговой реализации приведённого ниже алгоритма, который будет охватывать 95% основного алгоритма. Авторы отмечают два дополнительных улучшения: 1) Раннее преувеличение и 2) Адаптивная скорость обучения. Мы обсудим только добавление ранних преувеличений, поскольку это наиболее способствует интерпретации внутренней работы реальных алгоритмов, поскольку адаптивная скорость обучения направлена на повышение скорости сходимости.

1. Входные и выходные данные

Следуя оригинальной статье, мы будем использовать общедоступный набор данных MNIST из OpenML с изображениями рукописных цифр от 0 до 9.[2] Мы также произвольно выберем 1000 изображений из набора данных и уменьшим размерность набора данных с помощью анализа основных компонентов (PCA) и сохраним 30 компонентов. Оба они предназначены для увеличения вычислительного времени алгоритма, поскольку код здесь оптимизирован не для скорости, а скорее для интерпретируемости и обучения.

from sklearn.datasets import fetch_openml
from sklearn.decomposition import PCA
import pandas as pd

# Fetch MNIST data
mnist = fetch_openml('mnist_784', version=1, as_frame=False)
mnist.target = mnist.target.astype(np.uint8)

X_total = pd.DataFrame(mnist["data"])
y_total = pd.DataFrame(mnist["target"])

X_reduced = X_total.sample(n=1000)
y_reduced = y_total.loc[X_total.index]

# PCA to keep 30 components
X = PCA(n_components=30).fit_transform(X_reduced) 

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

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

Что касается выходных данных, напомним, что мы ищем низкоразмерное отображение исходного набора данных X. В этом примере мы будем отображать исходное пространство в двумерное пространство. Таким образом, нашим новым результатом будет 1000 изображений, которые теперь представлены в 2-мерном пространстве, а не в исходном 30-мерном пространстве:

2. Вычисление сходства X в исходном пространстве

Теперь, когда у нас есть исходные данные, первым шагом является вычисление попарных сходств в исходном многомерном пространстве. То есть, для каждого изображения i мы вычисляем вероятность того, что я выбрал бы изображение j в качестве его соседа в исходном пространстве для каждого j. Эти вероятности вычисляются с помощью нормального распределения, центрированного вокруг каждой точки, а затем нормализуются для суммирования до 1. Математически мы имеем:

Обратите внимание, что в нашем случае с n = 1000 эти вычисления приведут к получению матрицы оценок сходства 1000 x 1000. Мы устанавливаем p = 0 всякий раз, когда i = j b / c, мы моделируем попарное сходство. Однако вы можете заметить, что мы не упомянули, как определяется σ. Это значение определяется для каждого наблюдения i с помощью поиска по сетке на основе заданной пользователем желаемой запутанности распределений. Мы поговорим об этом непосредственно ниже, но давайте сначала посмотрим, как мы будем программировать eq. (1):

def get_original_pairwise_affinities(X:np.array([]), 
                                     perplexity=10):

    '''
    Function to obtain affinities matrix.
    '''

    n = len(X)

    print("Computing Pairwise Affinities....")

    p_ij = np.zeros(shape=(n,n))
    for i in range(0,n):
        
        # Equation 1 numerator
        diff = X[i]-X
        σ_i = grid_search(diff, i, perplexity) # Grid Search for σ_i
        norm = np.linalg.norm(diff, axis=1)
        p_ij[i,:] = np.exp(-norm**2/(2*σ_i**2))

        # Set p = 0 when j = i
        np.fill_diagonal(p_ij, 0)
        
        # Equation 1 
        p_ij[i,:] = p_ij[i,:]/np.sum(p_ij[i,:])

    # Set 0 values to minimum numpy value (ε approx. = 0) 
    ε = np.nextafter(0,1)
    p_ij = np.maximum(p_ij,ε)

    print("Completed Pairwise Affinities Matrix. \n")

    return p_ij

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

где H(P) – энтропия Шеннона для P.

В нашем случае мы установим perplexity = 10 и установим, что пространство поиска определяется как [0,01 * стандартное отклонение норм для разницы между изображениями i и j, 5 * стандартное отклонение норм для разницы между изображениями i и j], разделённое на 200 равных шагов. Зная это, мы можем определить нашу функцию grid_search() следующим образом:

def grid_search(diff_i, i, perplexity):

    '''
    Helper function to obtain σ's based on user-specified perplexity.
    '''

    result = np.inf # Set first result to be infinity

    norm = np.linalg.norm(diff_i, axis=1)
    std_norm = np.std(norm) # Use standard deviation of norms to define search space

    for σ_search in np.linspace(0.01*std_norm,5*std_norm,200):

        # Equation 1 Numerator
        p = np.exp(-norm**2/(2*σ_search**2)) 

        # Set p = 0 when i = j
        p[i] = 0 

        # Equation 1 (ε -> 0) 
        ε = np.nextafter(0,1)
        p_new = np.maximum(p/np.sum(p),ε)
        
        # Shannon Entropy
        H = -np.sum(p_new*np.log2(p_new))
        
        # Get log(perplexity equation) as close to equality
        if np.abs(np.log(perplexity) - H * np.log(2)) < np.abs(result):
            result = np.log(perplexity) - H * np.log(2)
            σ = σ_search
    
    return σ

Учитывая эти функции, мы можем вычислить аффинное преобразование через p_ij = get_original_pairwise_affinities(X), где мы получаем следующую матрицу:

Обратите внимание, что диагональным элементам задаётся значение ε ≈ 0 по конструкции (всякий раз, когда i = j). Напомним, что ключевым расширением алгоритма t-SNE является вычисление совместных вероятностей, а не условных вероятностей. Это вычисляется просто следующим образом:

Таким образом, мы можем определить новую функцию:

def get_symmetric_p_ij(p_ij:np.array([])):

    '''
    Function to obtain symmetric affinities matrix utilized in t-SNE.
    '''
        
    print("Computing Symmetric p_ij matrix....")

    n = len(p_ij)
    p_ij_symmetric = np.zeros(shape=(n,n))
    for i in range(0,n):
        for j in range(0,n):
            p_ij_symmetric[i,j] = (p_ij[i,j] + p_ij[j,i]) / (2*n)
    
    # Set 0 values to minimum numpy value (ε approx. = 0)
    ε = np.nextafter(0,1)
    p_ij_symmetric = np.maximum(p_ij_symmetric,ε)

    print("Completed Symmetric p_ij Matrix. \n")

    return p_ij_symmetric

Вводя p_ij сверху, мы имеем p_ij_symmetric = get_symmetric_p_ij(p_ij), где мы получаем следующую диаграмму сродств:

Теперь мы завершили первый главный шаг в t-SNE! Мы вычислили симметричную матрицу аффинности в исходном многомерном пространстве. Прежде чем мы перейдём непосредственно к этапу оптимизации, мы обсудим основные компоненты задачи оптимизации на следующих двух этапах, а затем объединим их в нашем цикле for.

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

Теперь мы хотим выбрать случайное начальное решение в пространстве меньшей размерности следующим образом:

def initialization(X: np.array([]),
                   n_dimensions = 2):

    return np.random.normal(loc=0,scale=1e-4,size=(len(X),n_dimensions))

где вызывая y0 = initialization(X), мы получаем случайное начальное решение:

Теперь мы хотим вычислить диаграмму сродства в этом пространстве меньшей размерности. Однако напомним, что мы делаем это, используя распределение Стьюдента-t с 1 степенью свободы:

Опять же, мы устанавливаем q = 0 всякий раз, когда i = j. Обратите внимание, что это уравнение отличается от уравнения eq. (1) в том смысле, что знаменатель находится над всеми i и, следовательно, симметричен по конструкции. Вводя это в код, мы получаем:

def get_low_dimensional_affinities(Y:np.array([])):
    '''
    Obtain low-dimensional affinities.
    '''

    n = len(Y)
    q_ij = np.zeros(shape=(n,n))

    for i in range(0,n):

        # Equation 4 Numerator
        diff = Y[i]-Y
        norm = np.linalg.norm(diff, axis=1)
        q_ij[i,:] = (1+norm**2)**(-1)

    # Set p = 0 when j = i
    np.fill_diagonal(q_ij, 0)

    # Equation 4 
    q_ij = q_ij/q_ij.sum()

    # Set 0 values to minimum numpy value (ε approx. = 0)
    ε = np.nextafter(0,1)
    q_ij = np.maximum(q_ij,ε)

    return q_ij

Здесь мы ищем диаграмму сродства 1000 x 1000, но теперь в пространстве более низких измерений. Вызывая q_ij = get_low_dimensional_affinities(y0), мы получаем:

4. Вычисление градиента функции затрат

Напомним, что наша функция затрат представляет собой расхождение Кулбека-Лейблера между совместными распределениями вероятностей в пространстве высокой размерности и пространстве низкой размерности:

Интуитивно мы хотим свести к минимуму разницу в матрицах сродства p_ij и q_ij, тем самым наилучшим образом сохранив структуру “окрестностей” исходного пространства. Задача оптимизации решается с помощью градиентного спуска, но сначала давайте рассмотрим вычисление градиента для приведённой выше функции затрат. Авторы выводят (см. приложение А к статье) градиент функции затрат следующим образом:

В python у нас есть:

def get_gradient(p_ij: np.array([]),
                q_ij: np.array([]),
                Y: np.array([])):
    '''
    Obtain gradient of cost function at current point Y.
    '''

    n = len(p_ij)

    # Compute gradient
    gradient = np.zeros(shape=(n, Y.shape[1]))
    for i in range(0,n):

        # Equation 5
        diff = Y[i]-Y
        A = np.array([(p_ij[i,:] - q_ij[i,:])])
        B = np.array([(1+np.linalg.norm(diff,axis=1))**(-1)])
        C = diff
        gradient[i] = 4 * np.sum((A * B).T * C, axis=0)

    return gradient  

Вводя соответствующие аргументы, мы получаем градиент при y0 через gradient = get_gradient(p_ij_symmetric,q_ij,y0) с соответствующим выводом:

Теперь у нас есть всё необходимое для решения задачи оптимизации!

5. Повторение и оптимизация низкоразмерного отображения

Чтобы обновить наше низкоразмерное отображение, мы используем градиентный спуск с импульсом, как описано авторами:

где η – наша скорость обучения, а α (t) – наш импульсный член в зависимости от времени. Скорость обучения определяет размер шага на каждой итерации, а член импульса позволяет алгоритму оптимизации набирать инерцию в плавном направлении пространства поиска, не увязая при этом в зашумлённых частях градиента. Мы установим η=200 для нашего примера и зафиксируем α(t)=0,5, если t < 250, и α(t)=0,8 в противном случае. У нас есть все компоненты, необходимые выше для вычисления в соответствии с правилом обновления, таким образом, теперь мы можем выполнить нашу оптимизацию за заданное количество итераций T (мы установим T = 1000).

Прежде чем мы перейдём к итерационной схеме, давайте сначала представим усовершенствование, которое авторы называют “ранним преувеличением”. Этот член является константой, которая масштабирует исходную диаграмму сродств p_ij. Оно уделяет больше внимания моделированию очень похожих точек (высокие значения в p_ij из исходного пространства) в новом пространстве на ранней стадии и, таким образом, формированию “кластеров” из очень похожих точек. Раннее преувеличение включается в начале схемы итерации (T<250), а затем отключается в противном случае. В нашем случае раннее преувеличение будет равно 4. Мы увидим это в действии на нашем рисунке ниже после реализации.

Теперь, собрав все части алгоритма воедино, мы получаем следующее:

def tSNE(X: np.array([]), 
        perplexity = 10,
        T = 1000, 
        η = 200,
        early_exaggeration = 4,
        n_dimensions = 2):
    
    n = len(X)

    # Get original affinities matrix 
    p_ij = get_original_pairwise_affinities(X, perplexity)
    p_ij_symmetric = get_symmetric_p_ij(p_ij)
    
    # Initialization
    Y = np.zeros(shape=(T, n, n_dimensions))
    Y_minus1 = np.zeros(shape=(n, n_dimensions))
    Y[0] = Y_minus1
    Y1 = initialization(X, n_dimensions)
    Y[1] = np.array(Y1)

    print("Optimizing Low Dimensional Embedding....")
    # Optimization
    for t in range(1, T-1):
        
        # Momentum & Early Exaggeration
        if t < 250:
            α = 0.5
            early_exaggeration = early_exaggeration
        else:
            α = 0.8
            early_exaggeration = 1

        # Get Low Dimensional Affinities
        q_ij = get_low_dimensional_affinities(Y[t])

        # Get Gradient of Cost Function
        gradient = get_gradient(early_exaggeration*p_ij_symmetric, q_ij, Y[t])

        # Update Rule
        Y[t+1] = Y[t] - η * gradient + α * (Y[t] - Y[t-1]) # Use negative gradient 

        # Compute current value of cost function
        if t % 50 == 0 or t == 1:
            cost = np.sum(p_ij_symmetric * np.log(p_ij_symmetric / q_ij))
            print(f"Iteration {t}: Value of Cost Function is {cost}")

    print(f"Completed Embedding: Final Value of Cost Function is {np.sum(p_ij_symmetric * np.log(p_ij_symmetric / q_ij))}")
    solution = Y[-1]

    return solution, Y

Вызывая решение, Y = tSNE(X), мы получаем следующий результат :

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

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

Заключение

И вот оно, мы написали его с нуля! Я надеюсь, что вы сочли это упражнение проливающим свет на внутреннюю работу системы и, как минимум, приносящим удовлетворение. Обратите внимание, что эта реализация предназначена не для оптимизации скорости, а скорее для понимания. Для повышения скорости вычислений и производительности были реализованы дополнения к алгоритму t-SNE, такие как варианты алгоритма Барнса-Хата (древовидные подходы), использующие PCA в качестве инициализации встраивания или использующие дополнительные расширения градиентного спуска, такие как адаптивные скорости обучения. Реализация в scikit-learn использует многие из этих улучшений.

Как всегда, я надеюсь, что вам понравилось читать это так же, как мне понравилось это писать.

Ресурсы

[1] van der Maaten, L.J.P.; Hinton, G.E. Visualizing High-Dimensional Data Using t-SNE. Journal of Machine Learning Research 9:2579–2605, 2008.

[2] LeCun et al. (1999): The MNIST Dataset Of Handwritten Digits (Images) License: CC BY-SA 3.0