Участник:Asenin/Многомерное шкалирование

Автор: Сенин Александр Николаевич, студент ММП ВМК МГУ (417)

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

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

Для неподготовленного читателя вкратце опишем, как обычно ставятся задачи машинном обучении. В классическом машинном обучении имеется множество объектов, для каждого объекта предполагается некоторый ответ (целевая переменная). Считаем, что существует некоторая зависимость между объектами и ответами, в общем случае она неизвестна. Общая задача машинного обучения и состоит в том, чтобы восстановить эту зависимость: для каждого объекта предсказать соответствующую ему целевую переменную. Обычно, у нас уже есть некоторые знания об этой зависимости, чаще всего они выражаются в совокупности прецедентов: пар (объект, целевая переменная). Такая совокупность называется обучающей выборкой. Предполагается, что в конечном счете мы для любого объекта будем уметь возвращать целевую переменную. Существенная часть алгоритмов машинного обучения сужает понятие объекта до конечного вектора - признакового описания, получает матрицу объекты-признаки, каждая строка такой матрицы соответствует объекту, и каждой строке соответствует целевая переменная.

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

Постановка получилась сильно общей, вполне естественно существование более узких постановок и разных подходов к решению. Резюмируем общую задачу: по вещественной матрице попарных различий найти для каждого объекта такое положение в пространстве размерности [math]p[/math] (чаще всего [math]p=2,3[/math]), что попарные различия будут лучше всего сохранены. Тогда общее описание алгоритма: получаем на вход квадратную вещественную матрицу [math]D[/math] размера [math]N x N[/math], возвращаем [math]N[/math] векторов размерности [math]p[/math] - координаты точек, которые лучше всего описывают наши объекты.

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

Конкретизируем подход к решению задачи многомерного шкалирования. Будем рассматривать задачу классического многомерного шкалирования (Classical Multidimensional Scaling, cMDS).

Пусть входная матрица различий - матрица попарных расстояний для евклидовой метрики [math]D = (d_{ij})[/math]. Здесь предполагаем, что данные нам различия во входной матрице - расстояния, с точки зрения математики, причем посчитанные евклидовой метрикой. Это допущение позволит найти конфигурацию, которая идеально точно воспроизводит расстояния на парах, но про размерность полученных координат мы не сможем ничего сказать. Однако, даже если метрика была не евклидовой, или даже не была формально математическим расстоянием, а также в случае, если требуются координаты фиксированной небольшой размерности, мы все равно получим результат.

Задача классического MDS (cMDS) - найти [math]X = (x_1, ..., x_N)^T[/math], т.ч. [math]d_{ij} = ||x_i - x_j||_2[/math].

Решение не единственно: [math]X^* = X + c^T[/math] тоже решение, т.к. [math]d_{ij} = ||x_i - x_j|| = ||(x_i + c) - (x_j + c)||[/math]. Ищем центрированную конфигурацию [math]\overline{x} = 0[/math]. Матрица [math]D[/math] евклидова, т.е. [math]\exists \{x_i\}_{i=1}^{N} \in \ R^p[/math], т.ч. [math]d_{ij}^2 = (x_i - x_j)^T(x_i - x_j)[/math]

Идея: восстановить [math]\{x_i\}_{i=1}^{N} \in R^p[/math], при условии [math]\overline{x} = 0[/math]. Решение: восстановим матрицу Грама [math]B = (b_{ij}),[/math] где [math]b_{ij} = x_i^T x_j[/math]

Обозначим [math]X = (x_1, ..., x_N)^T \Rightarrow B = XX^T[/math]

Таким образом, если получим матрицу Грама, то сможем по ее спектральному разложению получить [math]X[/math].

Восстановление матрицы Грама в cMDS:

[math]d_{ij}^2 = (x_i - x_j)^T(x_i - x_j) = x_i^Tx_i + x_j^Tx_j - 2x_i^Tx_j = b_{ii} + b_{jj} - 2b_{ij} \\ \overline{x} = 0 \Rightarrow \sum_{i=1}^N b_{ij} = 0 \\ \dfrac{1}{N} \sum\limits_{i=1}^{N} d_{ij}^2 = \dfrac{1}{N} \sum\limits_{i=1}^{N} b_{ii} + b_{jj} \\ \dfrac{1}{N} \sum\limits_{j=1}^{N} d_{ij}^2 = b_{ii} + \dfrac{1}{N} \sum\limits_{j=1}^{N} b_{jj} \\ \dfrac{1}{N^2} \sum\limits_{i, j=1}^{N} d_{ij}^2 = \dfrac{2}{N} \sum\limits_{i=1}^{N} b_{ii} \\ b_{ij} = -\dfrac{1}{2}(d_{ij}^2 - d_{i \bullet}^2 - d_{\bullet j}^2 + d_{\bullet \bullet}^2)[/math]

Строим по [math]D[/math] матрицу Грама [math]B[/math]:

[math]b_{ij} = -\dfrac{1}{2}(d_{ij}^2 - d_{i \bullet}^2 - d_{\bullet j}^2 + d_{\bullet \bullet}^2)[/math]

[math]B = C_N A C_N[/math], где [math]A = -\dfrac{1}{2} D^2[/math] поэлементно, [math]C_n = E - \dfrac{1}{n} \textbf{1} \textbf{1}^T[/math]

Выше получено, что [math]B = XX^T,[/math] [math]X \in R^{N \times p} \Rightarrow B[/math] симметричная, неотрицательно определенная, ранга [math]rg{B} = rg{XX^T} = rg{X} = p[/math]

[math]B[/math] имеет [math]p[/math] положительных собственных значений и [math]n-p[/math] нулевых собственных значений.

[math]B = \Gamma \Lambda \Gamma^T[/math], где [math]\Lambda = diag(\lambda_1, ..., \lambda_p)[/math] и [math]\Gamma = (\gamma_1, ..., \gamma_p)[/math] матрица собственных векторов.

Итоговый результат:

[math]X = \Gamma \Lambda^{\frac{1}{2}}[/math]

Важное замечание: таким образом, мы найдем конфигурацию точек некоторой фиксированной, вообще говоря, нами не выбираемой, размерности [math]p[/math], которая может быть сильно больше целевых для задач визуализации [math]p=2,3[/math] (и многих других задач). Существуют теоремы, которые показывают, что можно смотреть на величину собственного значения, как на показатель полезности соответствующего собственного вектора, с точки зрения сохранения набольшей информации. Другими словами, наибольший собственный вектор описывает направление наибольшего разброса данных, то есть проецирование данных на это направление сохранит больше всего информации. Тогда, если полученное с помощью алгоритма представление данных [math]X[/math] обладает неудовлетворительно большой размерностью, то смело можем оставить в матрице [math]\Gamma[/math] только то число собственных векторов, какую размерность точек хотим получить, причем будем оставлять вектора, отвечающие наибольшим собственным значениям.

Рассмотрим теперь случай, когда расстояния были посчитаны произвольной, необязательно евклидовой метрикой. Существуют соответстующая теорема, которая обещает нахождение конфигурации с точными равенствами попарных расстояний для евклидовой метрики, но для неевклидовой метрики существование такой конфигурации обычно не гарантируется. Этому случаю соответствуют, например, существование отрицательных собственных значений в спектральном разложении. Здесь решение аналогичное: находим матрицу Грама, находим ее спектральное разложение, а далее оставляем все те же [math]p[/math] собственных векторов, отвечающих наибольшим собственным значениям.

Далее будем рассматривать конкретное приложение классического многомерного шкалирования во многих задачах машинного обучения, а именно понижение размерности. Эта процедура применятся к данным для снижения вычислительной сложности последующих алгоритмов, которые будут с этими данными оперировать. В общей постановке многомерного шкалирования считается, что у нас от данных остаются только попарные расстояния, однако, для нашей цели этот лишний шаг избыточен: если у нас уже известны полные данные, нам нет смысла дополнительно считать матрицу попарных расстояний (обычно огромного размера), мы можем сразу посчитать матрицу Грама, а по ней найти спектральное разложение, оставить в разложении [math]p[/math] собственных векторов, отвечающих наибольшим собственным значениям, и спроецировать на них (с предварительным умножением на корни собственных значений) наши данные. Другими словами, используем всю доступную нам информацию: вместо вычисления матрицы попарных расстояний, восстановления по ней матрицы Грама, будем сразу же считать матрицу Грама.

Наше соображение о том, что в спектральном разложении нужно брать не все собственные векторы, чтобы получить представление нужной размерности, позволяет упростить задачу: нам не требуется строить полное спектральное разложение, а именно находить все собственные векторы и собственные значений, достаточно найти лишь несколько наибольших собственных значений и соответствующих им собственных векторов.

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

- Получаем на вход матрицу объекты-признаки [math]X[/math].

- Строим матрицу Грама попарных скалярных произведений [math]XX^T[/math].

- Находим [math]p[/math] ее наибольших собственных значений и отвечающие им собственные векторы.

- Организуем собственные векторы в столбцы матрицы [math]\Gamma[/math], корни собственных значений располагаем на диагонали матрицы [math]\Lambda^{\frac{1}{2}}[/math], причем в порядке невозрастания собственных значений.

- Возвращаем восстановленную матрицу [math]X = \Gamma \Lambda^{\frac{1}{2}}[/math].

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

- Получаем на вход матрицу попарных расстояний [math]D[/math].

- Возводим матрицу поэлементно в квадрат [math]D^2[/math].

- Находим [math]A = -\dfrac{1}{2} D^2[/math].

- Находим матрицу Грама [math]B = C_N A C_N[/math], где [math]C_n = E - \dfrac{1}{n} \textbf{1} \textbf{1}^T[/math].

В подавляющем большинстве задач задана именно выборка [math]X[/math], поэтому далее будем реализовывать параллельно именно первый алгоритм.

Конкретизируем теперь некоторые шаги последовательного алгоритма. Матрица Грама с точки зрения вычислений определяется вполне однозначно, то же самое и с вычислением [math]X = \Gamma \Lambda^{\frac{1}{2}}[/math] на последнем шаге. Основной вопрос в нахождении собственных векторов и собственных значений. Для нахождения всего спектра матрицы можно использовать и метод вращения Якоби, и алгоритм [math]QR[/math], но в нашей задаче эти методы избыточны, нам не требуется весь спектр, нам требуется лишь несколько наибольших его элементов. Здесь нам очень удачно подойдет степенной метод нахождения собственных значений. Это итерационный алгоритм поиска собственного значения с максимальной абсолютной величиной и одного из соответствующих собственных векторов для произвольной матрицы [math]A[/math].

- Берем некоторый начальный вектор [math]r_0[/math].

- Итеративно пересчитываем [math]r_{k+1} = \dfrac{Ar_k}{||Ar_k||}[/math] до сходимости, так находим собственный вектор.

- Вычисляем [math]\mu_k = \dfrac{r_k^T Ar_k}{r_k^T r_k}[/math] - соответствующее собственное значение.

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

- Вычисляем [math]A_1 = A - \lambda r_k r_k^T[/math] - матрицу, сохраняющую все собственные значения матрицы [math]A[/math], кроме [math]\lambda[/math]. В качестве [math]\lambda[/math] кладем найденное до этого наибольшее по модулю собственное значение.

- Приеняем предыдущие шаги для нахождения следующего по модулю собственного значения, и т.д.

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

Самый частый фрагмент алгоритма - итеративное нахождение собственных значений степенным методом.

- Берем некоторый начальный вектор [math]r_0[/math].

- Итеративно пересчитываем [math]r_{k+1} = \dfrac{Ar_k}{||Ar_k||}[/math] до сходимости, так находим собственный вектор.

- Вычисляем [math]\mu_k = \dfrac{r_k^T Ar_k}{r_k^T r_k}[/math] - соответствующее собственное значение.

По сути, ключевая операция здесь - умножение матрицы на вектор.

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

Алгоритм на довольно высоком уровне абстракции будет выглядеть следующим образом:

- Получаем на вход матрицу объекты-признаки [math]X[/math].

- Строим матрицу Грама попарных скалярных произведений [math]XX^T[/math].

- Находим [math]p[/math] ее наибольших собственных значений и отвечающие им собственные векторы.

- Организуем собственные векторы в столбцы матрицы [math]\Gamma[/math], корни собственных значений располагаем на диагонали матрицы [math]\Lambda^{\frac{1}{2}}[/math], причем в порядке невозрастания собственных значений.

- Возвращаем восстановленную матрицу [math]X = \Gamma \Lambda^{\frac{1}{2}}[/math].

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

1. Получаем на вход матрицу объекты-признаки [math]X[/math].

2. Строим матрицу Грама попарных скалярных произведений [math]XX^T[/math].

3. Находим [math]p[/math] ее наибольших собственных значений и отвечающие им собственные векторы, а именно:

3(а). Берем некоторый начальный вектор [math]r_0[/math].

3(b). Итеративно пересчитываем [math]r_{k+1} = \dfrac{Ar_k}{||Ar_k||}[/math] до сходимости, так находим собственный вектор.

3(c). Вычисляем [math]\mu_k = \dfrac{r_k^T Ar_k}{r_k^T r_k}[/math] - соответствующее собственное значение.

3(d). Вычисляем [math]A_{new} = A - \lambda r_k r_k^T[/math] - матрицу, сохраняющую все собственные значения матрицы [math]A[/math], кроме [math]\lambda[/math]. В качестве [math]\lambda[/math] кладем найденное до этого наибольшее по модулю собственное значение. Повторяем степенной метод.

4. Организуем собственные векторы в столбцы матрицы [math]\Gamma[/math], корни собственных значений располагаем на диагонали матрицы [math]\Lambda^{\frac{1}{2}}[/math], причем в порядке невозрастания собственных значений.

5. Возвращаем восстановленную матрицу [math]X = \Gamma \Lambda^{\frac{1}{2}}[/math].

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

Пусть входная матрица имеет размер [math]n[/math] x [math]k[/math], столбцы суть признаки, строки суть объекты. Могут быть разные вариации взаимного отношения [math]n[/math] и [math]k[/math], но чаще всего [math]n \gt k[/math], причем в подавляющем числе случаев [math]n \gt \gt k[/math], бывает даже на порядки. То есть о матрице [math]X[/math] можно мыслить, как о прямоугольной.

Обсудим сложность (2) шага. На этом шаге мы считаем скалярные произведения строк (объектов). Для непосредственного вычисления матрицы Грама нам не требуется явно транспонировать матрицу. Достаточно вычислить [math]nk[/math] скалярных произвведений, каждое скалярное произведение - это [math]n[/math] умножений и [math]n-1[/math] сложение. Тогда последовательная сложность на этом шаге [math]O(n^2k)[/math].

Обсудим сложность (3b) шага. Число размерности искомой конфигурации точек [math]p[/math] обычно [math]2,3[/math] для задач визуализации, и много меньше [math]k[/math] для задач понижения размерности. Поэтому множитель числа итераций здесь уйдет в константу О большого. Далее нам нужно итеративно до сходимости пересчитывать [math]r_{k+1}[/math]. На одной итерации нужно одно умножение матрицы на вектор, и вычисление нормы полученного результата, затем деление каждой компоненты результата на норму. Важно: теперь матрица [math]A[/math] - суть матрица Грама, а значит имеет размеры [math]n[/math] x [math]k[/math]. Вычисление [math]Ar_k[/math] требует [math]n^2[/math] умножений и сложений. Вычисление нормы по сути скалярный квадрат, то есть [math]n[/math] умножений и [math]n-1[/math] сложение, деление на норму еще [math]n[/math] умножений. Итоговая сложность этого шага [math]O(n^2)[/math], основные вычислительные затраты идут на умножение матрицы на вектор. Если получится распараллелить алгоритм так, что можно считать сложность умножения матрицы на вектор линейной, то от всего шага получится добится линейной сложности. Более того, мы ищем собственные значения итерационно, а значит этот шаг будет повторятся сильно чаще остальных, это тоже нужно держать в уме при оценки сложности всего алгоритма, а также при выборе стратегии распараллеливания.

Шаг (3c) представляет собой два скалярных произведения (умножение матрицы на вектор получено на предыдущем шаге), то есть имеет линейную сложность [math]O(n)[/math].

На шаге [math](3d)[/math] требуется по вектору посчитать матрицу, это [math]n^2[/math] умножений, а затем вычесть результат из матрицы предыдущего шага, это еще [math]n^2[/math] сложений. Последовательная сложность этого шага [math]O(n^2).[/math]

Шаг (4) может быть выполнен грамотной организацией структур данных еще на предыдущих шагах, из вычислительных операций здесь только взятие [math]p[/math] корней, в общей сложности этот шаг учитывать не будем.

Шаг (5) по сути представляет собой умножение [math]p[/math] векторов длины [math]n[/math] на число, можем считать его сложность [math]O(n)[/math].

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

См. изображения с информационной структурой вычисления [math]XX^T[/math] и [math]Ar_k[/math]

Информационная структура алгоритма вычисления XX^T

Информационная структура алгоритма вычисления Ar_k

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

Из обсуждения последовательной сложности алгоритма ясно следующее: просматривается последовательная структура на макроуровне, но на уровне каждой подзадачи оперируем с понятиями параллельной природы: матрицы, векторы, их совместные операции. По сложности выделяются первый этап вычисления матрицы Грама со сложностью [math]O(n^2 k)[/math] и умножение матрицы на вектор [math]O(n^2)[/math]. Причем, вторая операция имеет сильно больший приоритет - нам нужно умножать матрицу на вектор на каждой итерации степенного метода, а эти итерации будут повторяться до сходимости, то есть потенциально значительное число раз. Все остальные этапы по сути имеют линейную сложность [math]O(n)[/math], за исключением построения новой матрицы [math]A[/math], однако это операцию нужно провести всего [math]p[/math] раз ([math]p \lt \lt n[/math]) по сравнению с несравнимо большим число итераций для сходимости степенного метода.

Обсудим ресурса параллелизма вычисления [math]XX^T[/math]. По сути это [math]n^2[/math] скалярных произведений векторов длины [math]k[/math].

Для вычисления [math]XX^T[/math] нам потребуется выполнить следующие ярусы:

[math]n[/math] ярусов умножений и сложений, в каждом ярусе [math]nk[/math] операций

Таким образом, по высоте ЯПФ имеем линейную сложность, по ширине билинейную.

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

Обсудим ресурс параллелизма умножения матрицы на вектор [math]Ar_k[/math].

Для вычисления [math]Ar_k[/math] нам потребуется выполнить следующие ярусы:

[math]n[/math] ярусов умножений, в каждом ярусе [math]n[/math] операций

Таким образом, и по высоте ЯПФ, и по ширине ЯПФ имеем линейную сложность. Это важный результат, т.к. умножение матрицы на вектор - самая частая операция в нашей постановке алгоритма cMDS, потому что встречается на каждой итерации степенного метода.

Все остальные операции в последовательном алгоритме (см. шаги в предыдущем параграфе), как было выяснено, имеют линейную сложность, причем <<природного>> параллелизма в них нет. Тогда, если мы распараллелим вычисление матрицы Грама и умножение матрицы на вектор, то сложность алгоритма сведется к единразовому вычислению матрицы Грама с билинейной сложностью [math]O(nk)[/math], а затем вычислительное ядро алгоритма будет обладать линейной сложностью [math]O(n)[/math] (тут важно понимать, что число итераций степенного метода может несколько ухудшить эту оценку), [math]p[/math] раз потребуется вычислить [math]r_k r_k^T[/math] за [math]O(n^2)[/math] операций.

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

На вход поступает прямоугольная вещественная матрица размера [math]n[/math] на [math]k[/math], алгоритм возвращает матрицу размера [math]n[/math] на [math]p[/math], где [math]p[/math] - гиперпараметр, задаваемый пользователем.

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

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

Результаты запуска реализации cMDS на суперкомпьютере Ломоносов-2

Прокомментируем полученный результат. Действительно получилось значительно ускорить время работы алгоритма за счет распараллеливания самых частых операций вычислительного ядра алгоритма. Однако, сложность итогового алгоритма далека от линейной, причинами этому являются последовательная структура алгоритма на макроуровне (параллелизма удалось достичь лишь на фрагментах микроуровня алгоритма), а также число итераций степенного метода (помним, что на одной итерации за счет параллелизма мы можем достичь линейной сложности) все таки усложняет алгоритм, кроме того, дает вклад билинейная сложность вычисления [math]XX^T[/math] перед началом выполнения алгоритма. Увеличение сложности задачи (размера входной матрицы) закономерно приводит к значительному увеличению времени работы, однако параллелизм позволяет компенсировать это увеличение. Тем не менее, получен хороший результат: в случае матрицы [math]40000[/math] строк на [math]3000[/math] столбцов удалось сократить время работы с порядка 15 минут до около 1 минуты.

2.7 Существующие реализации алгоритма

Существующие реализации, использующие параллелизма на C++, автору неизвестны. В библиотеке sklearn для языка Python есть метод многомерного шкалирования, однако не гарантируется, что он решает задачу cMDS именно поиском собственных векторов матрицы Грама.

3 Литература

[1] Wikipedia contributors. (2020, April 13). Multidimensional scaling. In Wikipedia, The Free Encyclopedia. Retrieved 20:33, June 1, 2020, from https://en.wikipedia.org/w/index.php?title=Multidimensional_scaling&oldid=950612463