Участник: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 Информационный граф

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

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