Разбираемся в ML без воды: от базы до Attention. Часть 12: Понижение размерности и PCA. data science.. data science. pca.. data science. pca. svd.. data science. pca. svd. Машинное обучение.. data science. pca. svd. Машинное обучение. понижение размерности.. data science. pca. svd. Машинное обучение. понижение размерности. проклятие размерности.

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

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

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

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

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

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

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

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

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

Во-вторых, данные становятся “редкими”: представьте квадрат со стороной 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 фактически вращает старые оси координат (признаки). Он находит такое направление, вдоль которого данные вытянуты сильнее всего, и делает его первой, самой важной осью — первой главной компонентой.

PCA

PCA

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

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

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

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

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

Итак, прямую можно задать единичным вектором ||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_comega, получим

text{Var}(z)=frac{1}{n}(X_comega)^TX_comega=frac{1}{n}omega^TX_c^TX_comega

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

max_{||omega||=1}omega^TSigmaomega.

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

mathcal{L}(omega, lambda)=omega^TSigmaomega - lambda(omega^Tomega - 1)

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

2Sigmaomega - 2lambdaomega=0

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

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

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

  • Вычислить ковариационную матрицу для центрированных данных;

  • Найти собственные векторы ковариационной матрицы;

  • Отсортировать их по собственным значениям;

  • Взять первые k (число компонент);

  • Умножать данные на эти векторы и получить сжатый датасет размера mathbb{R}^{n times k} .

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

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

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

X_c=USV^T,

Где

  • U — ортогональная матрица размера n times n. Её столбцы — это собственные векторы матрицы XX^{T}.

  • V — ортогональная матрица размера d times d. Её столбцы — это собственные векторы матрицы X^{T}X

  • S — прямоугольная диагональная матрица размера 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

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

Sigma=VS^TISV^T=VS^TSV^T

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

Sigma=VLambda V^T

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

Sigma V=VLambda V^{T}V

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

Sigma V=V Lambda

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

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

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

Метод

Сложность по времени

Сложность по памяти

Проблема

1. Через ковариацию

O(n cdot d^2 + d^3)

O(d^2)

Кубическая сложность

2. Полный SVD

O(n cdot d^2)

O(n^2 + d^2)

Часто n большое

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

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

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

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

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

v=Sigma q=X^T_c(X_cq_1)

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

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

alpha_1=q_1^TSigma 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^TSigma q_1 -  alpha_1(q_1^Tq_1)

Так как q_{1} единичный (q_1^T q_1=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:

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

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

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}.

Запишем систему уравнений в общем виде для k векторов. Соберем векторы 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} слева. Так как все векторы q строго перпендикулярны друг другу, то Q_k^T Q_k=I_k (единичная матрица размера 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)

O(d^2)

2. Полный SVD

O(n cdot d^2)

O(n^2+d^2)

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

O(n cdot d cdot k)

O(n cdot k + d cdot k)

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

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

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

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

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

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

Автор: ysrgsyn

Источник