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

Материал из 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 на каждом шаге, обеспечивающая его ортогональность ранее найденным главным компонентам.

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

Сингулярное разложение тензоров и тензорный метод главных компонент

Часто вектор данных имеет дополнительную структуру прямоугольной таблицы (например, плоское изображение) или даже многомерной таблицы - то есть тензора: x_{i_{1}i_{2}...i_{q}}, 1 \leq i_{j} \leq n_j. В этом случае также эффективно применять сингулярное разложение. Определение, основные формулы и алгоритмы переносятся практически без изменений: вместо матрицы данных имеем q+1-индексную величину \operatorname{X}=(x_{i_{0}i_{1}i_{2}...i_{q}}), где первый индекс i_{0}-номер точки (тензора) данных.

Основная процедура — поиск наилучшего приближения тензора x_{i_{0}i_{1}i_{2}...i_{q}} тензором вида a^0_{i_{0}} a^1_{i_{1}}a^2_{i_{2}}...a^q_{i_{q}} (где a^0=(a^0_{i_{0}})m-мерный вектор (m - число точек данных), a^l=(a^l_{i_{l}}) — вектор размерности n_l при l>0) методом наименьших квадратов:

F= \frac{1}{2}\sum_{i_{0}=1}^m \sum_{i_{1}=1}^{n_1}...\sum_{i_{q}=1}^{n_q} (x_{i_{0}i_{1}...i_{q}} - a^0_{i_{0}} a^1_{i_{1}}...a^q_{i_{q}})^2 \to \min

Решение этой задачи дается последовательными итерациями по явным формулам. Если заданы все векторы-сомножители кроме одного a^k_{i_{k}}, то этот оставшийся определяется явно из достаточных условий минимума.

a^k_{i_{k}}= \frac{\sum_{i_{0}=1}^m \sum_{i_{1}=1}^{n_1}... \sum_{i_{k-1}=1}^{n_{k-1}}\sum_{i_{k+1}=1}^{n_{k+1}}...\sum_{i_{q}=1}^{n_{q}}x_{i_{0}i_{1}...i_{k-1}i_{k}i_{k+1}...i_{q}} a^0_{i_{0}}a^{k-1}_{i_{k-1}}a^{k+1}_{i_{k+1}}...a^q_{i_{q}}}{\prod_{j\neq k}\|a^j\|^2 }\, .

B качестве начального приближения векторов a^l=(a^l_{i_{l}}) (l>0) возьмем случайные векторы единичной длины, вычислим вектор a^0, далее для этого вектора a^0 и данных векторов a^2 , a^3, ...вычисляем вектор a^1 и т. д. (циклически перебирая индексы) Каждый шаг уменьшает значение F(b, a). Алгоритм, очевидно, сходится. В качестве критерия остановки используется малость относительного уменьшения значения минимизируемого функционала F за цикл или малость самого значения F. Далее, из тензора \operatorname {X} вычитаем полученное приближение a^0_{i_{0}} a^1_{i_{1}}a^2_{i_{2}}...a^q_{i_{q}} и для остатка вновь ищем наилучшее приближение этого же вида и т. д., пока, например, норма очередного остатка не станет достаточно малой.

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


Примечания


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