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

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

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

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

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

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

t-распределенное стохастическое соседнее встраивание ( t-SNE) — это инструмент уменьшения размерности, который в основном используется в наборах данных с большим пространственным пространством признаков и позволяет визуализировать данные вниз или спроецировать их в пространство меньшей размерности (обычно 2-мерное). Д). Это особенно полезно для визуализации нелинейно разделимых данных, когда линейные методы, такие как анализ главных компонентов (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. Отображение: итеративно улучшайте точки в пространстве меньшего размера, пока расхождения Кульбака-Лейблера между всеми условными вероятностями не будут минимизированы. По сути, мы минимизируем различия в вероятностях между матрицами сходства двух пространств, чтобы обеспечить наилучшее сохранение сходства при сопоставлении исходного набора данных высокой размерности с набором данных низкой размерности.

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

Реализация с нуля

Давайте теперь перейдем к пониманию t-SNE путем реализации оригинальной версии алгоритма, представленной в статье Лоренса ван дер Маатен и Джеффри Хинтона. Сначала мы начнем с пошаговой реализации алгоритма 1 ниже, который покроет 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. В этом примере мы будем отображать исходное пространство в двумерное пространство. Таким образом, наш новый вывод будет состоять из 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

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

В нашем случае мы установим недоумение = 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 σ

Обратите внимание, диагональные элементы установлены равными ε ≈ 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

Теперь мы завершили первый основной шаг в 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))

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

Опять же, мы устанавливаем q = 0 всякий раз, когда i = j . Обратите внимание, что это уравнение отличается от уравнения. (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

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

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

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

В питоне имеем:

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

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

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

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

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

Заключение

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

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

Ресурсы

[1] ван дер Маатен, LJP; Хинтон, Г. Э. Визуализация многомерных данных с использованием t-SNE. Журнал исследований машинного обучения 9: 2579–2605, 2008 г.

[2] ЛеКун и др. (1999): Лицензия на набор рукописных цифр (изображений) MNIST: CC BY-SA 3.0

Получите доступ ко всему коду через этот репозиторий GitHub: https://github.com/jakepenzak/Blog-Posts

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