Разбираемся в ML без воды: от базы до Attention. Часть 12: Понижение размерности и PCA / Хабр

В предыдущей части мы разобрали градиентный бустинг — финального босса в классическом обучении с учителем. Мы научились строить мощные ансамбли, которые выжимают максимум из табличных данных. Кажется, что на этом можно ставить точку и прыгать в современный мир нейросетей и Deep Learning.

Но до этого момента мы жили в идеальной теплице: у нас всегда была разметка (тот самый target, который нужно предсказать), а количество признаков в таблицах было разумным. В реальности все иначе. Данных часто слишком много, в них куча шума, а правильных ответов никто не разметил.

В этой части мы закроем очередную проблему в классическом ML — столкнемся лицом к лицу с проклятием размерности (curse of dimensionality). Поймем, как сжимать многомерные пространства, не теряя важный смысл, и как заставить машину самостоятельно группировать объекты в кластеры, вообще не имея готовых классов.

Проклятие размерности

Хорошо ли иметь много столбцов в датасете? На первый взгляд, кажется, да. Чем больше признаков, тем больше потенциально полезной информации может содержаться в данных. Закономерностей найдется больше, модель (хоть дерево, хоть линейная регрессия) получит больше информации для поиска закономерностей. Одним словом — профит. И это, в целом, истина.

Но у монетки есть и обратная сторона, по имени проклятием размерности (Curse of Dimensionality). Этот термин придумал математик Ричард Беллман для динамического программирования, но он идеально вписался и в мире машинного обучения. Разберем в чем его суть:

При увеличении количества признаков возникают несколько неприятных моментов.

Первое и, пожалуй, самое очевидное — вычисления становятся дороже: чем больше матрица признаков, тем больше памяти нужно для её хранения. Кроме того, нужно выполнить больше вычислений при работе с этой матрицей.

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

Во-вторых, данные становятся "редкими": представьте квадрат со стороной 1.

В 1D это отрезок длины 1.
В 2D — квадрат площадью 1.
В 10D — гиперкуб объёма 1.

Теперь возьмём маленькую область со стороной 0.1.

В 1D её "объём" равен 0.1, или 10% от всего пространства;
В 2D: 0.01, или же 1%;
В 10D: 10 $^{-10}$ , то есть 0....01%.

Получается, общее пространство вокруг разрастается до гигантских масштабов, а область с нашими данными превращается в микроскопическую точку. Вокруг остается один пустой "воздух". Чтобы заполнить этот 10-мерный куб данными с той же плотностью, что и обычный отрезок, нам понадобится экспоненциально больше строк в датасете. Откуда найти столько данных? — непонятно.

Из этого вытекает третья проблема: расстояния теряют свой смысл.

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

$\frac{d_{\max}-d_{\min}}{d_{\min}}\to 0$

Самая дальняя точка оказывается почти на таком же расстоянии, как и самая близкая. То есть понятие "сосед" становится менее информативным.

Также, есть и четвертая проблема: возможность переобучения:

Если признаков много, модель получает огромную свободу.

Например у нас 100 объектов и 100500 признаков. Модель может найти гиперплоскость, которая идеально разделит обучающую выборку, но на новых данных всё развалится.

Поэтому чем выше размерность, тем больше данных требуется. А данные, увы, не всегда валяются под соседним кустом.

Понижение размерности

Как мы увидели, проблем у высокой размерности достаточно много. И один из классических способов борьбы с ними звучит почти гениально: раз высокая размерность создаёт проблемы, давайте её понизим.

Однако делать это нужно со вкусом. Самый грубый подход иногда действительно работает: мы можем без сожалений удалить столбец с группой крови предыдущих владельцев из датасета, который используем для предсказания цен на дома. Иными словами, можно и нужно избавляться от признаков, которые вообще не связаны с поставленной задачей.

И да, формально это тоже является понижением размерности, активно используется на практике и называется отбором признаков (feature selection). Минус такого подхода заключается в том, что в конечном счёте мы своими руками избавляемся от части информации.

Но есть и другой подход — можно составить новые признаки, используя старые. Простым примером является создание столбца "площадь" из столбцов "длина" и "ширина". Этот метод тоже имеет свое название — извлечение признаков (feauture extraction).

Посмотрим как это работает, на примере ветерана машинного обучения — метода PCA.

Principal component analysis (PCA)

Итак, Метод главных компонент (Principal component analysis, PCA) — это способ преобразовать исходные признаки так, чтобы получить новое пространство, в котором информация "упакована" более компактно.

Возьмем датасет с автомобилями. Если у всех машин пробег равен ровно 100500 км, этот признак, очевидно, бесполезен. Но если пробег варьируется от 5000 до 100000 км, разброс огромен, и этот признак уже может быть информативным. PCA ищет такие направления в данных, вдоль которых этот разброс максимален, и объявляет их новыми осями — главными компонентами.

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

Посмотрим как это работает с математической точки зрения.

Возьмем датасет (матрицу) размерности $n \times d$ , где n — количество объектов (строк), а — количество признаков (столбцов).

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

$\mu_{j} = \frac{1}{n}\sum_{i=1}^{n}x_{ij}$

И отнимаем её от всех значений $x_j := x_j - \mu_j$ . Тем самым центр превращается в начало координат. Теперь надо думать на какую прямую лучше всего спроецировать наши точки. Чтобы не путаться в обозначениях, напишем полученное как , где — центрированное.

Итак, прямую можно задать единичным вектором $||\omega|| = 1, \ \omega \in \mathbb{R}^d$ . А проекции точки на прямую имеют вид $z_i = x_{i}\omega$ , где $i = 1, \dots, n$ . А мы хотим найти такое направление, при котором разброс этих проекций был максимальным, то есть задача выглядит так:

$\max_{||\omega||=1}\text{Var}(z)$

И вот здесь нам как раз пригодится центрирование, ибо дисперсия $\text{Var}(z)$ раскроется просто как

$\text{Var}(z) = \frac{1}{n}\sum_{i=1}^{n}z_i^2 = \frac{1}{n} z^Tz$

Подставляя $z = X_c\omega$ , получим

$\text{Var}(z) = \frac{1}{n}(X_c\omega)^TX_c\omega = \frac{1}{n}\omega^TX_c^TX_c\omega$

Вспоминаем ковариационную матрицу: $\Sigma = X^TX$ . Получаем итоговую задачу PCA:

$\max_{||\omega||=1}\omega^T\Sigma\omega.$

Дальше все достаточно естественно — люди мы простые, видим условную задачу, вводим Лагранжиан (подробнее об этом было в части 8).

$\mathcal{L}(\omega, \lambda) = \omega^T\Sigma\omega - \lambda(\omega^T\omega - 1)$

По стандартному методу, берем производную по $\omega$ , приравниваем к нулю:

$2\Sigma\omega - 2\lambda\omega = 0$

Отсюда получаем $\Sigma\omega = \lambda\omega$ . И это не просто уравнение, это классическое определение собственного вектора. Получается, что $\omega$ собственный вектор матрицы $\Sigma$ , а $\lambda$ — собственное значение.

Теперь мы получили главную идею: $\omega^T\Sigma\omega = \lambda$ , если $\omega$ — собственный вектор. А значит для задачи максимизации мы просто берем собственный вектор с наибольшим собственным значением.

То есть, весь PCA это:

Вычислить ковариационную матрицу для центрированных данных;
Найти собственные векторы ковариационной матрицы;
Отсортировать их по собственным значениям;
Взять первые (число компонент);
Умножать данные на эти векторы и получить сжатый датасет размера $\mathbb{R}^{n \times k}$ .

Казалось бы, всё отлично, но раз мы занимаемся понижением размерности, то фич у нас многовато. А из этого следует, что ковариационная матрица будет огромной. Именно поэтому, на практике почти всегда обходятся без вычисления ковариационной матрицы, а используют метод под названием SVD.

Singular Value Decomposition / Сингулярное разложение

Суть SVD заключается в том, что мы можем разложить матрицу на вращение, сжатие и на ещё одно вращение.

Где

— ортогональная матрица размера $n \times n$ . Её столбцы — это собственные векторы матрицы $XX^{T}$ .
— ортогональная матрица размера $d \times d$ . Её столбцы — это собственные векторы матрицы $X^{T}X$
— прямоугольная диагональная матрица размера $n \times d$ . На её главной диагонали стоят неотрицательные числа, упорядоченные по убыванию $(s_1 \ge s_2 \ge \dots \ge 0)$ . Это сингулярные числа, которые равны корням из собственных значений $(s_i = \sqrt{\lambda_i})$ .

Если подставить такое разложение в ковариационную матрицу $\Sigma = X^TX$ , получим

$\Sigma = (USV^T)^T(USV^T) = VS^TU^TUSV^T$

Вспомним, что матрица называется ортогональной, если , где — единичная матрица. Поэтому,

$\Sigma = VS^TISV^T = VS^TSV^T$

Так как матрица диагональная, то , а само произведение дает нам диагональную матрицу, с элементами $s_{ij} = \lambda_{i}$ , при , и $s_{ij} = 0$ , при $i\ne j$ . Обозначим эту матрицу собственных значений как $\Lambda$ . Стало быть,

$\Sigma = V\Lambda V^T$

Остается финальный шаг: умножим обе части этого равенства на матрицу справа:

$\Sigma V=V\Lambda V^{T}V$

Так как матрица V ортогональна, то и , а значит

$\Sigma V= V \Lambda$

Фактически, это то же самое уравнение $\Sigma \omega = \lambda \omega$ , к которому мы пришли через Лагранжиан.

Матрица , которую SVD берет напрямую из наших центрированных данных $X_{c}$ , состоит из тех самых собственных векторов ковариационной матрицы, которые мы искали для максимизации дисперсии. То есть таким образом мы пропустили вычисление $X^{T}X$ .

Посмотрим на разницу:

Метод	Сложность по времени	Сложность по памяти	Проблема
1. Через ковариацию	$O(n \cdot d^2 + d^3)$		Кубическая сложность
2. Полный SVD	$O(n \cdot d^2)$		Часто большое

Как говорится, хрен редьки не слаще. Однако, всё это время мы брали все собственные вектора, а затем оставили из них первые . Для решения этой проблемы существует метод Ланцоша.

Метод Ланцоша

В чем суть? Раз мы не обязаны найти все собственные векторы огромной матрицы, мы можем спроецировать её на подпространство гораздо меньшей размерности k, где и найдем нужное нам количество.

Пусть у нас есть ковариационная матрица $\Sigma = X_c^T X_c$ размера $d \times d$ . Мы физически не можем её вычислить и хранить в памяти. Наша цель — найти её первый собственный вектор $\omega _{1}$ (главную компоненту).

Возьмем случайный единичный вектор размера и умножим его на ковариационную матрицу справа. Назовем результат как.

$v = \Sigma q = X^T_c(X_cq_1)$

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

Теперь измерим, какая часть от вектора $\Sigma q_1$ сонаправлена с нашим исходным вектором ․ Для этого найдем их скалярное произведение и обозначим полученное число как $\alpha_{1}$ :

$\alpha_1 = q_1^T\Sigma q_1$

Также вычтем из $\Sigma q_1$ его проекцию на $q_{1}$ , чтобы получить строго перпендикулярный остаток. Назовем этот остаток вектором $r_{1}$ :

$r_1 = \Sigma q_1 - \alpha_1q_1$

Посмотрим на следующее произведение:

$q_1^Tr_1 = q_1^T(\Sigma q_1 - \alpha_1q_1) = q_1^T\Sigma q_1 - \alpha_1(q_1^Tq_1)$

Так как $q_{1}$ единичный (), а $q_1^T \Sigma q_1 = \alpha_1$ , то получаем:

$q_1^Tr_1 = \alpha_1 - \alpha_1 = 0$

Как и ожидалось, вектор $r_{1}$ перпендикулярен $q_{1}$ .

Обозначим длину вектора $r_{1}$ как $\beta_1$ ( $\beta_1 = \Vert r_1 \Vert$ ). Отнормируем его, чтобы получить второй единичный вектор:

$q_{2}=\frac{r_{1}}{\beta_{1}}\implies r_{1}=\beta_{1}q_{2}$

Подставим это обратно в наше уравнение:

$\beta_{1}q_{2}=\Sigma q_{1}-\alpha_{1}q_{1}\implies \Sigma q_{1}=\alpha_{1}q_{1}+\beta_{1}q_{2}$

Теперь у нас есть два взаимно перпендикулярных вектора $q_{1}$ и $q_{2}$ . Умножим матрицу $\Sigma$ на наш новый вектор $q_{2}$ : $\Sigma q_{2}$

Аналогично предыдущему шагу, этот вектор будет как-то проецироваться на уже имеющиеся у нас направления $q_{1}$ и $q_{2}$ . Найдем эти проекции:

Проекция на $q_{2}$ : $\alpha_2 = q_2^T \Sigma q_2$
Проекция на $q_{1}$ : $q_1^T \Sigma q_2$

По свойствам транспонирования $q_1^T \Sigma q_2 = (\Sigma q_1)^T q_2$ . Подставим сюда найденное нами выражение для $\Sigma q_1$ :

$(\Sigma q_{1})^{T}q_{2}=(\alpha_{1}q_{1}+\beta_{1}q_{2})^{T}q_{2}=\alpha_{1}(q_{1}^{T}q_{2})+\beta_{1}(q_{2}^{T}q_{2})$

Итак, $q_{1}$ и $q_{2}$ перпендикулярны, длина $q_{2}$ равна 1. Нетрудно заметить, что наше выражение превращается просто в число $\beta_{1}$ .

Красота симметрии: проекция нового вектора $\Sigma q_2$ на самый первый вектор $q_{1}$ оказалась в точности равна числу $\beta _{1}$ , которое мы посчитали на прошлом шаге!

Теперь снова вычтем из $\Sigma q_2$ обе эти проекции, чтобы получить следующий перпендикулярный остаток $r_{2}$ :

$r_{2}=\Sigma q_{2}-\alpha_{2}q_{2}-\beta_{1}q_{1}$

Измеряем длину остатка $\beta_2 = \Vert r_2 \Vert$ и находим следующий единичный вектор :

$q_3 = \frac{r_2}{\beta_2} \implies r_2 = \beta_2q_3$

Отсюда получим

$\Sigma q_2 = \beta_1q_1 + \alpha_2q_2 + \beta_2q_3$

Аналогично для :

$\Sigma q_n = \beta_{n-1}q_{n-1} + \alpha_nq_n + \beta_{n}q_{n+1}$

То есть вектор $\Sigma q_n$ зависит только от $q_{n-1}$ и $q_{n+1}$ .

Запишем систему уравнений в общем виде для векторов. Соберем векторы $q_1, \dots, q_k$ в прямоугольную матрицу $Q_{k}$ размера $d \times k$ .

Тогда всю эту цепочку уравнений можно записать как одно матричное равенство:

$\Sigma Q_{k}=Q_{k}T_{k}$

Матрица $T_{k}$ заполняется нашими числами $\alpha$ и $\beta$ согласно уравнению. Поскольку каждый шаг зависит только от двух соседей, числа ложатся строго по трем диагоналям, а все остальные элементы матрицы нулевые:

$T_{k}=\left(\begin{matrix}\alpha {1}&\beta {1}&0&0\\ \beta {1}&\alpha {2}&\beta {2}&0\\ 0&\beta {2}&\alpha {3}&\beta {3}\\ 0&0&\beta {3}&\alpha {4}\end{matrix}\right)$

Умножим уравнение $\Sigma Q_k = Q_k T_k$ на $Q_{k}^{T}$ слева. Так как все векторы строго перпендикулярны друг другу, то (единичная матрица размера $k \times k$ ). Получаем:

$Q_{k}^{T}\Sigma Q_{k}=T_{k}$

Мы с нуля вывели, как превратить гигантскую ковариационную матрицу в крошечную трехдиагональную матрицу размера $k \times k$ .

Теперь компьютеру не составит труда собственные векторы $s_{i}$ этой матрицы $T_{k}$ :

$T_{k}s_{i}=\theta_{i}s_{i}$

Но давайте поймем куда делись наши "омеги и лямбды".

Подставим вместо $T_{k}$ выражение $Q_k^T \Sigma Q_k$ :

$(Q_{k}^{T}\Sigma Q_{k})s_{i}=\theta_{i}s_{i}$

Умножим на $Q_{k}$ слева:

$Q_k Q_k^T \Sigma Q_k s_i = \theta_i Q_k s_i$

Поскольку вектор $\Sigma Q_k s_i$ по определению метода Ланцоша целиком лежит внутри нашего подпространства $Q_{k}$ , операция проекции $Q_{k}Q_{k}^{T}$ оставляет его неизменным и полностью исчезает. Получим:

$\Sigma (Q_{k}s_{i})=\theta_{i}(Q_{k}s_{i})$

Сравниваем это с нашим Лагранжианом $\Sigma \omega = \lambda \omega$ :

Наши искомые омеги (собственные векторы) — это просто $\omega_i = Q_k s_i$ .
Наши лямбды (собственные значения) — это $\lambda_i = \theta _{i}$ .

Заключение

Напоследок, добавим этот метод в нашу таблицу для сравнения.

Метод	Сложность по времени	Сложность по памяти
1. Через ковариацию	$O(n \cdot d^2 + d^3)$
2. Полный SVD	$O(n \cdot d^2)$
3. Метод Ланцоша	$O(n \cdot d \cdot k)$	$O(n \cdot k + d \cdot k)$

Думаю, комментарии излишни.

На если запустить алгоритм в таком виде, то через 3–4 шага из-за микроскопических погрешностей в вычислениях свойство ортогональности векторов начинает теряться (в численных методах это называют loss of orthogonality). Новые векторы перестают быть перпендикулярными самым первым, и алгоритм ломается.

Именно поэтому на практике (включая тот же scikit-learn и с его ARPACK) никто не запускает метод Ланцоша "в лоб". К нему обязательно прикручивают так называемую процедуру переортогонализации.

Её смысл в том, что алгоритм заставляют принудительно проверять и выравнивать перпендикулярность нового вектора не только по отношению к двум соседям, а вообще ко всем ранее найденным векторам.

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

На практике встречаются и другие методы, но мы пока остановимся на этих. А также существуют и другие методы понижения размерности в целом. О них и поговорим в дальнейшем!