Простой итерационный алгоритм сингулярного разложения

Материал из MachineLearning.

Перейти к: навигация, поиск

Простой итерационный алгоритм сингулярного разложения матриц допускает простую высокопараллельную (в том числе, нейросетевую) реализацию. Сингулярное разложение матриц (англ. Singular value decomposition) необходимо для решения многих задач анализа данных. В том числе, анализ главных компонент сводится к сингулярному разложению матрицы центрированных данных.

Идея сингулярного разложения матрицы данных

Если \operatorname{X}=\left\{x_1,..., x_m \right\}^T — матрица, составленная из векторов-строк центрированных данных, то  C = \frac{1}{m-1}\operatorname{X}^T\operatorname{X} и задача о спектральном разложении ковариационной матрицы  C превращается в задачу о сингулярном разложении матрицы данных \operatorname{X}.

Число  \sigma \geq 0 называется сингулярным числом матрицы \operatorname{X} тогда и только тогда, когда существуют правый и левый сингулярные векторы: такие m-мерный вектор-строка b_{\sigma} и n-мерный вектор-столбец a_{\sigma} (оба единичной длины), что выполнено два равенства:

\operatorname{X} a_{\sigma} = \sigma b_{\sigma}^T \;;\, \, b_{\sigma} \operatorname{X}= \sigma a_{\sigma}^T.

Пусть p= \operatorname{rang} \operatorname{X} \leq \min\{n,m\}ранг матрицы данных. Сингулярное разложение матрицы данных \operatorname{X} — это её представление в виде

\operatorname{X}= \sum_{l=1}^p \sigma_l b_l^T a_l^T \;; \;\operatorname{X}^T= \sum_{l=1}^p \sigma_l a_l b_l \; \left(x_{ij}=\sum_{l=1}^p \sigma_l b_{li}a_{lj}\right),

где \sigma_l > 0 — сингулярное число, a_l=(a_{li}), \, i=1,... n — соответствующий правый сингулярный вектор-столбец, а b_l=(b_{li}), \, i=1,... m — соответствующий левый сингулярный вектор-строка (l=1,...p). Правые сингулярные векторы-столбцы a_l, участвующие в этом разложении, являются векторами главных компонент и собственными векторами эмпирической ковариационной матрицы C =\frac{1}{m-1}\operatorname{X} ^T \operatorname{X} , отвечающими положительным собственным числам  \lambda_l=\frac{1}{m-1}\sigma_l^2 > 0 .

Хотя формально задачи сингулярного разложения матрицы данных и спектрального разложения ковариационной матрицы совпадают, алгоритмы вычисления сингулярного разложения напрямую, без вычисления спектра ковариационной матрицы, более эффективны и устойчивы [1]. Это следует из того, что задача сингулярного разложения матрицы \operatorname{X} лучше обусловлена, чем задача разложения матрицы C =\frac{1}{m-1}\operatorname{X} ^T \operatorname{X} : для ненулевых собственных и сингулярных чисел

\frac{\lambda_{\max}}{\lambda_{\min}}= \left(\frac{\sigma_{\max}}{\sigma_{\min}}\right)^2.

Теория сингулярного разложения была создана Дж. Дж. Сильвестром (англ. J. J. Sylvester) в 1889 г. и изложена во всех подробных руководствах по теории матриц [1].

Простой итерационный алгоритм сингулярного разложения

Основная процедура — поиск наилучшего приближения произвольной m \times n матрицы X=(x_{ij}) матрицей вида b \otimes a = (b_i a_j) (где bm-мерный вектор, а an-мерный вектор) методом наименьших квадратов:

F(b, a) = \frac{1}{2}\sum_{i=1}^m \sum_{j=1}^n (x_{ij} - b_i a_j )^2 \to \min

Решение этой задачи дается последовательными итерациями по явным формулам. При фиксированном векторе a=(a_j) значения b=(b_i) , доставляющие минимум форме F(b, a) , однозначно и явно определяются из равенств \partial F/ \partial b_i = 0 :

\frac{\partial F}{\partial b_i} = - \sum_{j=1}^n (x_{ij} - b_i a_j )a_j = 0; \;\; b_i = \frac{\sum_{j=1}^n x_{ij}  a_j}{\sum_{j=1}^n a_j^2 }\, .

Аналогично, при фиксированном векторе b =(b_ i) определяются значения a=(a_j) :

a_j = \frac{\sum_{i=1}^m b_i x_{ij} }{\sum_{i =1}^m b_i ^2 }\, .

B качестве начального приближения вектора a возьмем случайный вектор единичной длины, вычисляем вектор b, далее для этого вектора b вычисляем вектор a и т. д. Каждый шаг уменьшает значение F(b, a) . В качестве критерия остановки используется малость относительного уменьшения значения минимизируемого функционала F(b, a) за шаг итерации (\Delta F / F ) или малость самого значения F.

В результате для матрицы X=(x_{ij}) получили наилучшее приближение матрицей P_1 вида b^1 \otimes a^1 = (b_i^1  a_j^1) (здесь верхним индексом обозначен номер итерации). Далее, из матрицы X вычитаем полученную матрицу P_1, и для полученной матрицы уклонений X_1=X-P_1 вновь ищем наилучшее приближение P_2 этого же вида и т. д., пока, например, норма X_k не станет достаточно малой. В результате получили итерационную процедуру разложения матрицы X в виде суммы матриц ранга 1, то есть X=P_1+P_2+... +P_q  \;  (P_l = b^l \otimes a^l) . Полагаем  \sigma_l = \|a^l\| \|b^l\| и нормируем векторы  a^l \, , \, b^l: a^l:= a^l/ \| a^l\|; \, \, b^l:= b^l/ \| b^l\|. В результате получена аппроксимация сингулярных чисел  \sigma_l и сингулярных векторов (правых —  a^l и левых — b^l).

К достоинствам этого алгоритма относится его исключительная простота и возможность почти без изменений перенести его на данные с пробелами[1], а также взвешенные данные.

Существуют различные модификации базового алгоритма, улучшающие точность и устойчивость. Например, векторы главных компонент a^l при разных l должны быть ортогональны «по построению», однако при большом числе итерации (большая размерность, много компонент) малые отклонения от ортогональности накапливаются и может потребоваться специальная коррекция a^l на каждом шаге, обеспечивающая его ортогональность ранее найденным главным компонентам.

Для квадратных симметричных положительно определённых матриц описанный алгоритм превращается в метод прямых итераций для поиска собственных векторов.

Примечания


Личные инструменты