Какие быстрые алгоритмы существуют для вычисления усеченного SVD?

Возможно, здесь не по теме, но уже существует несколько ( один , два ) связанных вопросов.

Поиски в литературе (или поиск в Google по усеченным алгоритмам SVD) обнаруживают множество статей, которые используют усеченные SVD по-разному, и утверждают (разочаровывающе, часто без цитирования), что существуют быстрые алгоритмы для его вычисления, но никто кажется, указывает на то, что эти алгоритмы.

Единственное, что я могу найти, - это один рандомизированный алгоритм , используемый в библиотеке redSVD .

То, что я хотел бы видеть, - это набор точных и неточных алгоритмов, подходящих для понимания того, как работают системы (но, конечно, не обязательно для их фактической реализации!).

У кого-нибудь есть хорошая рекомендация для такого рода вещей?

Говоря в широком смысле, существует два подхода для вычисления разложения по собственным значениям или по сингулярным значениям. Один из подходов заключается в диагонализации матрицы, и это, по сути, приводит к одновременному получению всего разложения по собственным значениям / сингулярным значениям (всего спектра собственных значений), см. Краткий обзор здесь: Какие эффективные алгоритмы для вычисления разложения по сингулярным значениям (SVD)? Альтернативой является использование итерационного алгоритма, который выдает один (или несколько) собственных векторов за раз. Итерации могут быть остановлены после того, как желаемое количество собственных векторов было вычислено.

Я не думаю, что есть итерационные алгоритмы специально для SVD. Это связано с тем, что можно вычислить SVD матрицы B размером , выполнив собственное разложение квадратной симметричной ( n + m ) × ( n + m ) матрицы A = ( 0 B B ⊤ 0 ) . Поэтому вместо того , чтобы спросить , что алгоритмы вычисления усеченного СВД, вы должны спросить , что итерационный алгоритм вычислений eigendecomposition: алгоритм усеченного СВД ≈ итерационный алгоритм eigendecomposition $n\times m$$\mathbf B$$(n+m)\times(n+m)$

Самый простой итерационный алгоритм называется степенной итерацией и действительно очень прост:

1. Инициализировать случайный .$\mathbf x$
2. Обновить .$\mathbf x \gets \mathbf A\mathbf x$
3. Нормализовать .$\mathbf x \gets \mathbf x / \|\mathbf x\|$
4. Перейти к шагу № 2, если не сходятся.

Все более сложные алгоритмы в конечном счете основаны на идее мощной итерации, но становятся довольно сложными. Необходимая математика задается подпространствами Крылова . Алгоритмы - итерация Арнольди (для квадратных несимметричных матриц), итерация Ланцоша (для квадратных симметричных матриц) и их варианты, такие как, например, «неявно перезапущенный метод Ланцоша» и еще много чего.

Вы можете найти это, например, в следующих учебниках:

1. Голуб и Ван Лоан, Матричные вычисления
2. Trefethen & Bau, Численная линейная алгебра
3. Деммель, Прикладная числовая линейная алгебра
4. Саад, Численные методы для задач больших собственных значений

Все разумные языки программирования и пакеты статистики (Matlab, R, Python numpy, как вы его называете) используют одни и те же библиотеки Фортрана для выполнения разложения по собственным / сингулярным значениям. Это LAPACK и ARPACK . ARPACK означает ARnoldi PACKage, и это все об итерациях Арнольди / Ланцоша. Например, в Matlab есть две функции для SVD: svdвыполняет полную декомпозицию через LAPACK, и svdsвычисляет заданное число единичных векторов через ARPACK, и это на самом деле просто оболочка для eigsвызова в матрице с квадратными размерами.

Обновить

$\mathbf B$$\mathbf A$$\mathbf A$$\mathbf B$$\mathbf A$

Для этих методов также есть библиотека Фортрана, она называется PROPACK :

Программный пакет PROPACK содержит набор функций для вычисления сингулярного разложения больших и разреженных или структурированных матриц. Процедуры SVD основаны на алгоритме бидиагонализации Ланцоша с частичной реортогонализацией (BPRO).

Тем не менее, PROPACK кажется гораздо менее стандартным, чем ARPACK, и изначально не поддерживается в стандартных языках программирования. Он написан Расмусом Ларсеном, у которого есть большая 90-страничная статья 1998 года « Бидиагонализация Ланцоша» с частичной реортогонализацией с хорошим обзором. Благодаря @MichaelGrant через этот поток вычислительной науки SE .

Среди более свежих работ наиболее популярными, по-видимому, являются Baglama & Reichel, 2005, Augmented неявно перезапускает методы бидиагонализации Lanczos , что, вероятно, соответствует уровню техники. Спасибо @Dougal за предоставленную ссылку в комментариях.

Обновление 2

Действительно, в обзорном документе, который вы сами цитировали, есть совершенно другой подход, описанный вами: Halko et al. 2009, Нахождение структуры со случайностью: вероятностные алгоритмы построения приближенных матричных разложений . Я не знаю достаточно об этом, чтобы комментировать.

Я просто наткнулся на поток с помощью быстрого поиска SVD, поэтому я пытаюсь сам разобраться, но, возможно, вам стоит взглянуть на адаптивное перекрестное приближение (ACA).

$M$$M=\sum_{i=0}^k U_i\cdot V^T_i$$N\times N$$O(N)$

Опять же, это зависит от вашей проблемы, если это работает. Во многих случаях, с которыми я лично сталкиваюсь, ACA является очень полезным числовым инструментом.

Примечание. Я хотел написать это как комментарий, но поскольку я только создал эту учетную запись, у меня недостаточно репутации для комментариев ... Но публикация работает.

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

1. Прежде чем что-либо делать, подгоните простую модель (например, глобальное среднее + значения столбца и константы строки), и только после того, как вы это сделаете, вам следует перейти к использованию усеченного SVD для подгонки остатков.
2. Инициализируйте случайный вектор длины k (где это урезанный ранг) для каждой строки и столбца (для каждого фильма и пользователя в случае Netflix).
3. Удерживайте векторы строк фиксированными и обновите векторы столбцов, чтобы минимизировать ошибки относительно известных записей в матрице. Процедура приведена в коде Matlab в статье.
4. Держите векторы столбцов фиксированными и обновляйте векторы строк аналогичным образом.
5. Повторите 3 и 4, пока вы не сходитесь или не получаете достаточно хорошие результаты.