Участник:USKoroleva/t-sne

Свойства и структура алгоритмов

Общее описание алгоритма

T-распределенное стохастическое вложение соседей (t-Distributed Stochastic Neighbor Embedding) — это алгоритм, предназначенный для визуализации высокоразмерных данных. Он относится к классу алгоритмов нелинейного снижения размерности (non-linear dimensionality reduction). Основная задача t-SNE — преобразование многомерных данных в пространство более низкой размерности (обычно 2D или 3D) таким образом, чтобы сохранить структуру локальных соседств. Другими словами, алгоритм предназначен для того, чтобы точки данных, которые были близки друг к другу в исходном высокоразмерном пространстве, оставались близкими в новом пространстве, а далекие точки оставались далекими.

Алгоритм работает с множеством точек в высокоразмерном пространстве. Каждая точка может представлять собой любой вектор признаков (например, вектор word2vec, гистограмму пикселей изображения, набор генетических маркеров).

В отличие от линейных методов, таких как PCA (Principal Component Analysis), t-SNE способен улавливать сложные, нелинейные зависимости между данными.

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

В основе t-SNE лежит вероятностный подход. Он не вычисляет прямое расстояние между точками, а вместо этого моделирует распределения вероятностей, которые отражают схожести между точками.

Математическое описание алгоритма

Алгоритм работает в два этапа:

Построение распределения вероятностей в высокоразмерном пространстве.

Если дан набор из n объектов высокой размерности $$\mathbf{x}_1, \mathbf{x}_2, ..., \mathbf{x}_n,$$ то для каждой пары точек данных исходного пространства вычисляются вероятности $$p_{j|i}$$ вероятность того, что объект i выберет объект j в качестве соседа. Эта вероятность пропорциональна гауссовскому распределению, центрированному вокруг точки i. $$ p_{j|i} = \frac{\exp\left(-||\mathbf{x}_i - \mathbf{x}_j||^2 / 2\sigma_i^2\right)}{\sum_{k \neq i} \exp\left(-||\mathbf{x}_i - \mathbf{x}_k||^2 / 2\sigma_i^2\right)} $$ $$\sigma_i \text{— дисперсия гауссова распределения для объекта i}$$ Дисперсия подбирается для каждого объекта так, чтобы энтропия распределения достигла заданного значения (параметр перплексии), оно определяет количество ближайших соседей, которые точка будет учитывать при построении распределения.

Затем условные вероятности симметризуются, чтобы получить совместное распределение p в высокоразмерном пространстве: $$ p_{ij} = \frac{p_{j|i} + p_{i|j}}{2n} $$ Чем больше полученное значение, тем более похожи объекты в исходном пространстве.

Построение распределения вероятностей и оптимизация в пространстве низкой размерности.

Алгоритм t-SNE стремится получить отображение $$\mathbf{y}_{1},\dots,\mathbf{y}_{N}$$ в D-мерное пространство $$\mathbf{y}_{i}\in \mathbb{R}^{d}$$ которое отражает похожести, насколько это возможно.

В низкоразмерном пространстве для объектов определяется аналогичное распределение схожести. Однако вместо гауссова распределения используется t-распределение Стьюдента с одной степенью свободы(распределение Коши). Это ключевое отличие, которое решает проблему "скукоживания". $$ q_{ij} = \frac{(1 + \|\mathbf{y}_i - \mathbf{y}_j\|^2)^{-1}}{\sum_{k \neq l} (1 + \|\mathbf{y}_k - \mathbf{y}_l\|^2)^{-1}} $$ t-распределение имеет более тяжелые "хвосты", что позволяет точкам в низкоразмерном пространстве "разбегаться" и более эффективно моделировать умеренные расстояния.

Алгоритм итеративно обновляет положения точек в низкоразмерном пространстве, минимизируя расхождение Кульбака-Лейблера между двумя распределениями P и Q $$ KL(P \mid \mid Q) = \sum_{i \neq j} p_{ij} \log\frac{p_{ij}}{q_{ij}} $$ Минимизация осуществляется с помощью метода градиентного спуска.

Результатом оптимизации является отображение, которое отражает похожесть между объектами пространства высокой размерности.

Вычислительное ядро алгоритма

Вычислительное ядро №1: Построение матрицы схожестей в высокоразмерном пространстве

Вычисление попарных расстояний между всеми точками: O(N²)

Бинарный поиск дисперсий σᵢ для каждой точки для достижения заданной перплексии

Расчет и симметризация условных вероятностей

Вычислительное ядро №2: Минимизация расхождения Кульбака-Лейблера

Расчет попарных расстояний в низкоразмерном пространстве: O(N²)

Вычисление градиента

Обновление координат точек

Оба этапа имеют квадратичную сложность O(N²) относительно числа точек. Основные операции — вычисление попарных расстояний и вероятностных распределений.

Макроструктура алгоритма

Построение графа сходств в высокоразмерном пространстве

1. Расстояния: Рассчитать матрицу попарных расстояний для исходных данных.

2. Вероятности: На основе расстояний и заданной перплексии построить матрицу совместных вероятностей P, которая представляет собой взвешенный и симметризованный граф сходств.

Оптимизация низкоразмерного представления

1. Инициализация: Случайно разместить точки в целевом пространстве.

2. Цикл оптимизации (повторять T раз):

• Рассчитать сходства Q: По текущим позициям точек в низкоразмерном пространстве рассчитать новую матрицу сходств Q с помощью тяжелых хвостов распределения Стьюдента.

• Рассчитать градиент: Оценить направление улучшения, вычислив градиент между P и Q.

• Обновить позиции: Сдвинуть точки в направлении, минимизирующем расхождение между P и Q.

Схема реализации последовательного алгоритма

Входные параметры и данные: X - матрица высокоразмерных данных размерности n × D

perplexity - целевая перплексия (обычно между 5 и 50)

T - количество итераций

η - скорость обучения

d - размерность выходного пространства (обычно 2 или 3)

# 1. Предобработка данных (нормализация)
X = (X - np.mean(X, axis=0)) / np.std(X, axis=0)
# 2. Инициализация низкоразмерного представления
Y = np.random.normal(0, 1e-4, (n, d))
Y_prev = np.copy(Y)  # для хранения предыдущего состояния (импульс)
# 3. Вычисление попарных квадратов расстояний в высокоразмерном пространстве
D_high = np.zeros((n, n))
for i in range(n):
    for j in range(i+1, n):
        dist_sq = np.sum((X[i] - X[j])**2)
        D_high[i, j] = dist_sq
        D_high[j, i] = dist_sq

Вычисление совместных вероятностей P

def binary_search_perplexity(distances_i, target_perplexity, tol=1e-5, max_iter=50):
   Бинарный поиск для нахождения оптимальной дисперсии
   sigma_min, sigma_max = -np.inf, np.inf
   sigma = 1.0  # начальное предположение
   
   for iter in range(max_iter):
       # Вычисление условных вероятностей с текущей сигмой
       p_cond = np.exp(-distances_i * sigma)
       p_cond = p_cond / np.sum(p_cond)
       
       # Вычисление текущей перплексии
       entropy = -np.sum(p_cond * np.log2(p_cond + 1e-8))
       current_perplexity = 2 ** entropy
       
       # Корректировка границ поиска
       if abs(current_perplexity - target_perplexity) < tol:
           break
           
       if current_perplexity > target_perplexity:
           sigma_min = sigma
           sigma = sigma_max if np.isfinite(sigma_max) else sigma * 2
       else:
           sigma_max = sigma
           sigma = sigma_min if np.isfinite(sigma_min) else sigma / 2
   
   return p_cond
# Основной цикл вычисления P
P = np.zeros((n, n))
for i in range(n):
   p_conditional = binary_search_perplexity(D_high[i], perplexity)
   P[i] = p_conditional
# Симметризация
P = (P + P.T) / (2 * n)

Градиентный спуск

# Инициализация импульса
gain = np.ones((n, d))  # адаптивная скорость обучения
for t in range(T):
   # 1. Вычисление попарных расстояний в низкоразмерном пространстве
   dist_low = np.zeros((n, n))
   for i in range(n):
       for j in range(i+1, n):
           dist_sq = np.sum((Y[i] - Y[j])**2)
           dist_low[i, j] = dist_sq
           dist_low[j, i] = dist_sq  
   # 2. Вычисление матрицы Q (распределение Стьюдента)
   Q = 1.0 / (1.0 + dist_low)
   np.fill_diagonal(Q, 0)  # обнуляем диагональ
   Q = Q / np.sum(Q)
   Q = np.maximum(Q, 1e-12)   
   # 3. Вычисление градиента
   grad = np.zeros((n, d))
   for i in range(n):
       for j in range(n):
           if i != j:
               # Градиент
               factor = 4.0 * (P[i, j] - Q[i, j]) * Q[i, j] * (1.0 + dist_low[i, j])
               grad[i] += factor * (Y[i] - Y[j])    
   # 4. Обновление с импульсом и адаптивной скоростью обучения
   # Простой импульс
   Y_new = Y - eta * grad + momentum * (Y - Y_prev)   
   # 5. Обновление состояний для следующей итерации
   Y_prev = np.copy(Y)
   Y = Y_new

Последовательная сложность алгоритма

Основные операции:

1. Арифметические операции с вещественными числами
2. Вычисления exp и log
3. Операции с памятью
4. Сложность по фазам:

• Матрица расстояний в высокоразмерном пространстве:

Для каждой из O(n²) пар точек выполняется O(D) операций

• Совместные вероятности P: O(n²) с доминирующими exp/log

На каждой итерации:

Вычисление exp для n элементов: O(n)

Суммирование: O(n)

Вычисление энтропии: O(n)

Итого на точку: O(n)

Общая сложность фазы: O(n²) операций + O(n²) вычислений exp

• Градиентный спуск (T итераций): O(T·n²d)

Итеративный градиентный спуск

Детализация на одну итерацию:

Расстояния в низкоразмерном пространстве: O(n²d)

Вычисление Q: O(n²) (включая нормализацию)

Вычисление градиента: O(n²d)

Обновление позиций: O(nd)

Итого на итерацию: O(n²d)

Общая последовательная сложность: O(n²(D + T·d))

Последовательная сложность алгоритма

Основные операции:

1. Арифметические операции с вещественными числами
2. Вычисления exp и log
3. Операции с памятью
4. Сложность по фазам:

• Матрица расстояний в высокоразмерном пространстве:

Для каждой из O(n²) пар точек выполняется O(D) операций

• Совместные вероятности P: O(n²) с доминирующими exp/log

На каждой итерации:

Вычисление exp для n элементов: O(n)

Суммирование: O(n)

Вычисление энтропии: O(n)

Итого на точку: O(n)

Общая сложность фазы: O(n²) операций + O(n²) вычислений exp

• Градиентный спуск (T итераций): O(T·n²d)

Итеративный градиентный спуск

Детализация на одну итерацию:

Расстояния в низкоразмерном пространстве: O(n²d)

Вычисление Q: O(n²) (включая нормализацию)

Вычисление градиента: O(n²d)

Обновление позиций: O(nd)

Итого на итерацию: O(n²d)

Общая последовательная сложность: O(n²(D + T·d))

Информационный граф

Граф представляет собой направленную ациклическую структуру с 10 ярусами (уровнями) вычислений, организованными в виде конвейера. Каждый ярус соответствует логическому этапу алгоритма, причем операции в пределах одного яруса могут выполняться параллельно, а между ярусами существуют строгие последовательные зависимости данных.

Ярус 0 (Входные данные): Содержит два узла: матрицу исходных данных X размером 3×d (где d - размерность исходного пространства) и параметр perplexity. Эти данные являются исходными для всех последующих вычислений.

Ярус 1 (Вычисление попарных расстояний): Состоит из 6 параллельных узлов, каждый из которых вычисляет квадрат евклидова расстояния между парой точек: dist_0_1, dist_0_2, dist_1_0, dist_1_2, dist_2_0, dist_2_1. Каждая операция имеет вид ‖xᵢ - xⱼ‖² и выполняется независимо от других, что позволяет достичь максимального параллелизма на этом этапе.

Ярус 2 (Определение параметров σ): Содержит 3 параллельных узла sigma_0, sigma_1, sigma_2. Каждый узел выполняет бинарный поиск значения σᵢ для соответствующей точки i, обеспечивающего заданную перплексию. Эти вычисления также независимы друг от друга, но каждый узел sigma_i зависит от всех расстояний, исходящих из точки i.

Ярус 3 (Условные вероятности): Включает 6 узлов вычисления условных вероятностей p_j|i = exp(-‖xᵢ - xⱼ‖²/2σᵢ²)/Σ. Узлы организованы в шахматном порядке для минимизации перекрытий: p_1|0, p_2|0, p_0|1, p_2|1, p_0|2, p_1|2. Каждое вычисление зависит от соответствующего расстояния dist_i_j и параметра sigma_i.

Ярус 4 (Симметризация вероятностей): Содержит 3 узла p_01, p_02, p_12, которые выполняют симметризацию: pᵢⱼ = (pⱼ|ᵢ + pᵢ|ⱼ)/2N. Каждый узел зависит от двух соответствующих условных вероятностей (например, p_01 зависит от p_1|0 и p_0|1).

Ярус 5 (Инициализация): Единственный узел Y_init, создающий случайную инициализацию матрицы низкоразмерных представлений Y размером 3×2.

Ярус 6 (Низкоразмерные сходства): 3 параллельных узла q_01, q_02, q_12, вычисляющих t-распределенные сходства: qᵢⱼ = 1/(1+‖yᵢ - yⱼ‖²). Каждый узел зависит от Y_init и соответствующих координат yᵢ и yⱼ.

Ярус 7 (Вычисление градиентов): 3 параллельных узла grad_0, grad_1, grad_2, вычисляющих градиенты функции стоимости КЛ-дивергенции: ∂C/∂yᵢ = 4Σⱼ(pᵢⱼ - qᵢⱼ)(yᵢ - yⱼ)qᵢⱼ. Каждый градиент grad_i зависит от всех симметризованных вероятностей pᵢⱼ и низкоразмерных сходств qᵢⱼ, связанных с точкой i.

Ярус 8 (Обновление координат): 3 параллельных узла update_0, update_1, update_2, выполняющих обновление по правилу: yᵢ = yᵢ - η·∂C/∂yᵢ + α·Δyᵢ(t-1). Каждый узел зависит от соответствующего градиента grad_i.

Ярус 9 (Выходные данные): Финальный узел Y_final, содержащий полученное низкоразмерное представление данных размером 3×2.

Типы дуг и зависимостей:

Голубые сплошные линии - потоки данных от входных параметров к вычислительным узлам

Красные линии - зависимости вероятностей pᵢⱼ при вычислении градиентов

Зеленые линии - переход градиентов к обновлению координат

Серые пунктирные линии - остальные вычислительные зависимости

Критические пути: Основной критический путь проходит через цепочку: X → dist_0_1 → sigma_0 → p_1|0 → p_01 → grad_0 → update_0 → Y_final. Аналогичные пути существуют для других точек.

Ресурс параллелизма алгоритма

Ярусно-параллельная форма и параллельная сложность

1. Фаза инициализации и вычисления матрицы расстояний

Ярус 1: Вычисление попарных квадратов расстояний для всех (i,j), i < j

       - Степень параллелизма: O(n²)
       - Высота яруса: O(D) операций
       - Тип операций: арифметические (+, -, *)

2. Фаза вычисления матрицы P (совместных вероятностей)

Для каждой точки i независимо (параллелизм по i):

   Ярус 2.1: Бинарный поиск сигмы (последовательно, K итераций)
       - Степень параллелизма: O(n) (по точкам i)
       - Высота на точку: O(K·n) операций
       - K ≈ 50 итераций на практике
   
   Ярус 2.2: Симметризация P = (P + Pᵀ)/(2n)
       - Степень параллелизма: O(n²)
       - Высота яруса: O(1) операций

3. Фаза градиентного спуска (T итераций)

Ярус 3.1: Вычисление dist_low[i,j] = ||y_i - y_j||² для всех (i,j)

       - Степень параллелизма: O(n²)
       - Высота яруса: O(d) операций

Ярус 3.2: Вычисление Q[i,j] = (1 + dist_low[i,j])⁻¹ для всех (i,j)

       - Степень параллелизма: O(n²)
       - Высота яруса: O(1) операций (деление)

Ярус 3.3: Нормализация Q (суммирование и деление)

       - Редукционная операция: O(log n) ярусов
       - Высота: O(log n) операций

Ярус 3.4: Вычисление градиента для каждой точки i

       - Степень параллелизма: O(n)
       - Высота яруса: O(nd) операций

Ярус 3.5: Обновление позиций Y

       - Степень параллелизма: O(n)
       - Высота яруса: O(d) операций

Общая параллельная сложность: O(D + K·n + T·(d + log n))

Иерархическая структура параллелизма

Уровень 1: Массовый параллелизм (степень O(n²))

Вычисление попарных расстояний

Симметризация матриц

Вычисления на уровне пар объектов (i,j)

Уровень 2: Крупноблочный параллелизм (степень O(n))

Бинарный поиск для разных точек

Вычисление градиентов по точкам

Обновление позиций точек

Уровень 3: Векторный параллелизм (степень O(d))

Операции внутри вычисления одного расстояния

Обновление d-мерных векторов

Скошенный параллелизм Алгоритм обладает вертикальным скошенным параллелизмом внутри итераций градиентного спуска:

Вычисления на ярусах 3.1-3.5 могут быть частично перекрыты

Потоковый параллелизм при конвейерной обработке итераций

Однако присутствует строгая последовательность итераций:

Y_{t} должен быть полностью вычислен перед началом вычислений Y_{t+1}

Это создает критический путь длиной O(T·(d + log n))

Входные и выходные данные алгоритма

Основные входные данные представляют собой плотную вещественную матрицу объектов-признаков X размерности n × D, где n — количество объектов (точек данных), обычно от 100 до 10 000 для стандартного варианта алгоритма (при бо́льших n используется аппроксимативный метод Barnes-Hut), а D — размерность исходного пространства, чаще всего лежащая в диапазоне от 50 до 1000, хотя может достигать и десятков тысяч. Данные, как правило, предварительно нормализуются (среднее значение по каждому признаку приводится к нулю, стандартное отклонение — к единице), что является критически важным для стабильной работы алгоритма. Объем входных данных составляет n × D × sizeof(float64) байт. Помимо матрицы X, алгоритм требует задания ключевых гиперпараметров. Наиболее важным из них является перплексия (perplexity) — вещественное число, обычно выбираемое в диапазоне от 5 до 50 (стандартное значение — 30), которое определяет баланс между сохранением локальной и глобальной структуры данных. Другим обязательным параметром является целевая размерность d (обычно 2 или 3 для визуализации).

Основным результатом работы алгоритма является плотная вещественная матрица сниженной размерности Y размерности n × d (практически всегда n × 2). Её объем составляет n × d × sizeof(float64) байт. Эта матрица представляет собой искомое низкоразмерное отображение, где каждая строка соответствует координатам одного объекта в двумерном (или трехмерном) пространстве. Важно понимать, что Y сохраняет лишь относительные, а не абсолютные расстояния: близость точек в Y указывает на схожесть соответствующих объектов в исходном высокомерном пространстве, однако глобальная геометрия (например, точные попарные расстояния) не сохраняется. Выходные данные обычно центрированы (среднее значение близко к нулю) и имеют контролируемый алгоритмом масштаб. Помимо основного результата, алгоритм часто возвращает значение финальной функции потерь, а также историю значений функции потерь по итерациям.

Свойства алгоритма

Вычислительная мощность, определяемая как отношение числа операций к объему перемещаемых данных, составляет примерно O(D + T·d) операций на единицу данных. Это умеренное значение указывает на то, что алгоритм может быть чувствителен к пропускной способности памяти, особенно на этапе вычисления попарных расстояний, где объем промежуточных данных O(n²) значительно превышает объем входных O(nD). С точки зрения устойчивости алгоритм условно устойчив. Для обеспечения численной стабильности требуются дополнительные меры, такие как добавление малых констант (например, 1e-12) для предотвращения деления на ноль и вычисления логарифма от нуля, а также ограничение вероятностей снизу. Алгоритм чувствителен к выбросам в данных, поскольку сфокусирован на сохранении локальных расстояний, и его результат сильно зависит от корректного выбора параметра перплексии (perplexity).

Сбалансированность алгоритма неоднородна. По типам операций преобладают арифметические действия над дорогостоящими вычислениями экспонент и логарифмов. С точки зрения параллельной нагрузки, этап вычисления попарных расстояний хорошо сбалансирован, представляя O(n²) однородных задач, в то время как этап бинарного поиска оптимальной сигмы для каждой точки может требовать разного числа итераций, создавая потенциальный дисбаланс. Важным свойством является неполная детерминированность алгоритма. Разные результаты при повторных запусках могут быть вызваны тремя факторами: случайной инициализацией низкоразмерного представления, неоднозначностью критерия остановки в бинарном поиске и влиянием порядка суммирований на накопление ошибок округления при параллельном выполнении. Детерминированность можно восстановить, фиксируя seed генератора случайных чисел и используя детерминированные алгоритмы редукции.

Информационная структура алгоритма характеризуется наличием вершин с высокой степенью исхода, что соответствует операциям рассылки данных. Например, координаты каждой точки Y_i на каждой итерации используются в O(n) операциях вычисления расстояний до других точек. Это предъявляет высокие требования к эффективной организации доступа к данным, возможно, через стратегии кэширования или широковещательную рассылку в распределенных системах. В графе алгоритма присутствуют "длинные" дуги, наиболее значимые из которых образуют критический путь длиной в T итераций, соединяя последовательные состояния Y_t и Y_{t+1}. Эта зависимость ограничивает потенциальное ускорение от параллелизации, так как итерации градиентного спуска должны выполняться строго последовательно.

Для специализированной аппаратной реализации информационный граф допускает компактные укладки. Одним из подходов является временнóй конвейер, где различные стадии обработки соседних итераций выполняются одновременно. Для фиксированной итерации вычислительные вершины, соответствующие независимым парам точек (i, j), могут быть расположены в виде двумерной решетки, что позволяет эффективно использовать параллельные ресурсы. При работе с очень большими наборами данных матрицы сходств могут быть разбиты на блоки, что открывает возможность для блочных вычислений и применения коллективных операций типа reduce-scatter. Ключевыми практическими проблемами реализации остаются квадратичный рост потребления памяти O(n²), чувствительность к подбору гиперпараметров и отсутствие гарантии сходимости к глобальному оптимуму функции стоимости.

Программная реализация алгоритма

Особенности реализации последовательного алгоритма

При проектировании параллельной программы для алгоритма t-SNE необходимо системно учесть несколько ключевых аспектов, формирующих её архитектуру. В основе лежит иерархически организованный ресурс параллелизма. На верхнем уровне выделяется макроструктура алгоритма: независимая обработка данных для вычисления матрицы расстояний и совместных вероятностей (степень параллелизма O(n²)) и итеративный градиентный спуск (последовательная цепочка длиной T). Внутри этих этапов параллелизм раскрывается через структуру вложенных циклов: внешний цикл по точкам данных (i) и внутренний по парам (i, j). Это создает комбинацию массового параллелизма на уровне операций с парами объектов и параллелизма конечного на уровне обработки отдельных точек или блоков данных.

Распределение операций между вычислительными потоками или процессами может осуществляться по нескольким схемам. Наиболее естественной для этапа вычисления расстояний является блочно-циклическое распределение пар (i, j) или строк/столбцов формируемых матриц. Для этапа градиентного спуска, где основная нагрузка приходится на вычисление градиента для каждой точки i, эффективно распределять точки между процессами. Распределение данных должно следовать за распределением вычислений, минимизируя необходимость в коммуникации. Матрицы P и Q, имеющие размер O(n²), могут храниться распределенно, блоками, либо реплицироваться, если это позволяют ресурсы памяти. Критически важным является локальное хранение векторов Y_i для каждой точки, за обработку которой отвечает данный процесс.

Оценка коммуникационных затрат показывает, что наиболее тяжелыми операциями являются редукции при нормализации матрицы Q (требуется глобальная сумма элементов) и, в распределенной реализации, обмены данными для вычисления попарных расстояний между точками, принадлежащим разным процессам. Общий объем пересылок для наивной реализации может достигать O(T·n²), что делает коммуникацию потенциальным узким местом. Для реализации с использованием MPI рекомендуется использовать гибридный подход MPI+OpenMP: MPI-процессы работают с крупными блоками данных, а внутри каждого процесса OpenMP-нити распараллеливают обработку своего блока. Для эффективной работы на GPU (CUDA) следует максимально задействовать массовый параллелизм, поручая каждому потоку обработку одной пары (i, j) на этапах вычисления расстояний и градиента. Ключевая оптимизация — использование разделяемой памяти для промежуточных результатов и минимизация передач между host и device. Для многоядерных CPU с общей памятью (OpenMP) основное внимание уделяется векторизации внутренних циклов (авто-векторизация или директивы SIMD) и обеспечению локальности данных для эффективного использования кеш-памяти, например, через переупорядочивание циклов. Таким образом, успешная параллельная реализация t-SNE требует тщательного согласования иерархии параллелизма алгоритма с иерархией памяти и вычислительных ресурсов целевой платформы.

Масштабируемость алгоритма и его реализации

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

Итеративный характер оптимизации встраивания с использованием градиентного спуска.

Наличие глобальных зависимостей: каждая точка вложенного пространства зависит от всех остальных через силы притяжения/отталкивания.

Эти особенности напрямую сказываются на масштабируемости как самого алгоритма, так и его реализаций.

1. Точки барьерной синхронизации и глобальных операций

При типичной параллельной реализации t-SNE (например, в библиотеках FIt-SNE, OpenTSNE, cuML, Multicore-TSNE) можно выделить следующие критические участки, требующие синхронизации:

Вычисление попарных расстояний. Хотя расстояния можно вычислять локально для подмножеств точек, их последующая нормализация в softmax требует глобальной операции редукции (сумма по всем парам). Это приводит к барьерной синхронизации и ограничивает слабую масштабируемость. Вычисление градиентов вложения: В классическом t-SNE градиент для каждой точки зависит от всех остальных (особенно при вычислении сил отталкивания). Это требует алл-редьюс (all-reduce) или алл-гезер (all-gather) операций на каждом шаге итерации.

Оптимизация с импульсом (momentum):

Обновление координат точек и скоростей требует синхронизации после каждого шага градиентного спуска. Эти точки синхронизации становятся узкими местами при увеличении числа процессоров, особенно в распределённых системах с высокой задержкой сети.

2. Оценка сильной и слабой масштабируемости Сильная масштабируемость (fixed problem size) Рассмотрим фиксированный набор данных (например, MNIST: N=70000 точек). При увеличении числа процессоров наблюдается:

Ускорение до 4–8 ядер (в зависимости от реализации и архитектуры CPU/GPU), благодаря параллелизации вычисления расстояний и градиентов. Насыщение производительности на ~16–32 ядрах: дальнейшее увеличение числа ядер не приводит к значимому ускорению из-за затрат на синхронизацию и ограниченной параллельности в фазе softmax. Отрицательное ускорение при переходе к мультиузловым системам без специализированной оптимизации (например, деревьев ближайших соседей или аппроксимаций, как в FIt-SNE). Слабая масштабируемость (scaling problem size with processors)

При увеличении N можно достичь более устойчивой масштабируемости, особенно при использовании аппроксимированных версий t-SNE:

FIt-SNE с использованием FFT-accelerated interpolation и блочного разбиения данных демонстрирует линейную слабую масштабируемость до сотен ядер. GPU-реализации (cuML) показывают сильную масштабируемость на одном узле (до 4–8 GPU), но страдают при мультиузловом распределении из-за частых all-reduce операций.

3. Влияние иерархии памяти и блочной структуры

Как и в случае умножения матриц, производительность t-SNE чувствительна к иерархии памяти:

При размерах задачи, умещающихся в L3 cache или оперативной памяти одного узла, наблюдаются пики эффективности. При выходе за пределы RAM одного узла и переходе к распределённой памяти (например, в кластере) эффективность падает из-за: обмена данными по сети, фрагментации данных, несбалансированной загрузки. Особенно это проявляется при вычислении матрицы сходства P: даже при блочном распределении строк/столбцов матрицы, операции типа softmax требуют глобального знания суммы всех элементов.

5. Выводы Алгоритм t-SNE по своей природе плохо масштабируем в классическом виде из-за квадратичной сложности и глобальных зависимостей. Масштабируемость реализации может быть значительно улучшена за счёт: аппроксимаций (деревья, интерполяция, случайная выборка), эффективного распределения данных (блочные схемы), минимизации глобальных синхронизаций (асинхронные градиенты, локальные softmax). Слабая масштабируемость оказывается значительно лучше, чем сильная, особенно при использовании современных оптимизированных библиотек (FIt-SNE, cuML). Иерархия памяти играет ключевую роль: оптимальная производительность достигается при таком размере задачи, который помещается в кэш или локальную память узла. Таким образом, максимальная производительность t-SNE достигается не просто увеличением числа процессоров, а балансом между объёмом данных, степенью аппроксимации и архитектурой вычислительной системы.

Литература

van der Maaten L., Hinton G. Visualizing Data using t-SNE // Journal of Machine Learning Research. 2008. Vol. 9. P. 2579–2605.

Erich Schubert, Michael Gertz. Intrinsic t-Stochastic Neighbor Embedding for Visualization and Outlier Detection // SISAP 2017 – 10th International Conference on Similarity Search and Applications. — 2017.