Вычисление SVD на компьютере: алгоритмы и литература

[B@R5uk]

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

В книге Lloyd N. Trefethen, David Bau III - Numerical linear algebra (1997) в главе 5 Eigenvalues, в лекции 31 Computing the SVD подробно расписывается более-менее современный алгоритм, состоящий из двух этапов. Первый заключается в том, что исходная матрица путём последовательности преобразований Хаусхольдера приводится к бидиагональному виду. На втором этапе неким волшебным методом эта матрица раскладывается непосредственно на U, Σ и V. Авторы подробно расписывают сколько флопов нужно для работы первой фазы и как это всё можно чуть-чуть улучшить финтом ушами через QR-разложение, но ничего не говорят про алгоритм второй фазы, кроме того, что он является методом итерационного приближения со сверхлинейной сходимостью и что это обеспечивает его время работы на одну степень n быстрее, чем время первого этапа. В итоге, я понимаю, как реализовать первую часть, но не весь алгоритм целиком.

Буду очень благодарен качественной литературе по теме.

[zgemm]

Если первую часть осилили, то со второй - разберетесь (я не могу посоветовать где почитать, так как везде криво). Там есть два алгоритма:

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

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

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

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

Попробуйте по этим идеям снова перечитать Трефентена, думаю, тогда будет существенно понятнее.

-- 27.11.2025, 23:57 --

B@R5uk в сообщении #1710820 писал(а):

Хочется самому реализовать точный и достаточно быстрый вариант хотя бы для матриц общего вида.

отдельно прокомментирую то, что вы хотите быстрый вариант. Это не тривиально, если вы это хотите сделать на компьютерах, у которых скорость доступа к памяти отличается от скорости доступа к регистрам, а такие компьютеры - это практически все, кроме разве что какой-нибудь atmega. А вот тут нужно очень хорошо чувствовать BLASы, и постоянно держать в голове теорему старшего Воеводина про умножение больших матриц. Не зря сам код SVD в лапаке содержит сотни тысяч строк. То есть я могу советовать, и постараюсь помочь, но реально это очень не тривиально и быстрые алгоритмы не всегда получаются именно те, у которых минимальное число флопов.

[Евгений Машеров]

Такой алгоритм я видел в книге Уилкинсона и Райнша "Справочник алгоритмов на языке АЛГОЛ. Линейная алгебра", и даже перекладывал его на С. "Вторая часть" это QR-алгоритм для нахождения собственных значений, только доработанный для двух-, а не трёхдиагональной матрицы.
Есть параграф про вычисление SVD в Деммель, "Вычислительная линейная алгебра".

[zgemm]

У Вилкинсона нет метода деления двухдиагональной на две подзадачи, а у Деммеля примерно также не понятно как у Трефентена написано. Они в общем-то и на конференциях в то время парой ходили, и мыслили схоже. Мне кажется, что даже Тыртышников чуть понятнее в своем кратком курсе написал про зануление поддиагонального элемента, но я очень хорошо помню, что те, кто по его учебнику учились, на экзамене не могли внятно ответить как это делается, а я любил издеваться, задавая студенту такой вопрос :)

[B@R5uk]

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

Евгений Машеров, спасибо за наводку! Обе книжки качнул, судя по беглому просмотру, это то, что мне нужно. Буду осваивать.

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

[zgemm]

B@R5uk в сообщении #1711057 писал(а):

zgemm, мне скорее важнее точность, а не скорость.

Смотря что вы вкладываете в понятие "точность" - если вам нужна точность порядка где-то половины машинной, то да, то, что на Алголе у Вилкинсона - именно то, что Вам надо. Если же то, что вкладывают современные разработчики Лапака, то надо начинать с Деммеля, так как Трефентен вроде это толком не обсуждал, а при Вилкинсоне это еще не открыли.

(Оффтоп)

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

В настоящий момент Лапак используется в таком огромном числе кодов начиная от CFD и до ИИ, и с 1981 реально вылизывался, что даже понять все нюансы алгоритмов, что есть в современном лапаковском SVD не работая в группе Донгарры думаю, даже за год 100% занятости не реально. Я сильно был в теме до 2008, и по себе знаю сколько там на этот алгоритм сил угрохано, но с тех пор там еще много что было добавлено.

B@R5uk в сообщении #1711057 писал(а):

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

а никак с разреженными, от слова совсем. Любая бидиагонализация сразу приводит к плотной квадратной матрице слева и справа от бидиагональной, и дальнейшие итерации требуют кубической зависимости от минимальной размерности матрицы. То есть даже теоретические изыскания в этом направлении могут только на проценты ускорить суммарное решение этой задачи, поэтому за эти проценты никто гнаться не стал.

А вот если вам нужны, например, только сингулярные значения из какой-то области, или например, самые малые или самые большие значения - то тут есть простор где разгуляться, и надо читать в сторону Крыловских итерационных методов. У того же Деммеля кажется тоже было что-то. Довольно много есть в смежной документации на cuSparse от Nvidia.

[Евгений Машеров]

Прежде всего - это не наилучшие книжки, а просто те, которые я читал и в которых есть ответ на вопрос о вычислении SVD. Думается, есть и получше, просто я с ходу не назову. Что до разреженности - то тут нужны специфические алгоритмы, что-то вроде Ланцоша, использующие в качестве "базовой операции" умножение матрицы на вектор, что при сильной разреженности может быть намного быстрее. Только доработанные для сингулярного разложения прямоугольной матрицы, а не разложения по собственным значениям квадратной, связанной с прямоугольной, как $A=B^TB$. Или можно свести задачу о сингулярных значениях к задаче о собственных способом, описанным у Трефетена на странице 235. Памяти надо больше, чем при получении квадратной матрицы умножением, зато собственные числа такой матрицы не будут квадратами сингулярных исходной, так что влияние ошибок меньше.

[zgemm]

Евгений Машеров в сообщении #1711077 писал(а):

Или можно свести задачу о сингулярных значениях к задаче о собственных способом, описанным у Трефетена на странице 235.

только там умалчивается, что точнее, чем корень из машинной точности таким способом получить нельзя.

Причем еще в 1986 Саад предложил для этого GMRES, который может как линейный системы решать, так и сингулярные числа и векторы получать, и без того, что точность ограничена корнем из машинной, но только пока ТС не определится нужны ли ему все, или только какие-то - это нет смысла обсуждать - я писал причину выше.

[Евгений Машеров]

Именно так, как на странице 235 - это точность порядка машинной. Но за счёт дублирования матриц, а не умножения (где именно и случается потеря точности).

[B@R5uk]

Ну, кроме общего повышения квалификации и освоения SVD у меня сейчас стоит такая подзадача.

Есть почти квадратная матрица и она разрежена: число ненулевых элементов пропорционально её стороне. "Неквадратность" у матрицы возможна двоякая: как переопределённость, так и недоопределённость. Мне надо эту матрицу обращать (точнее решать систему с ней), но с обоих сторон (прямую и транспонированную матрицу). Это, очевидно, может приводить к проблемам, поэтому нужна регуляризация. SVD в этом плане решает все задачи одним выстрелом, даже регуляризация делается принятием осведомлённого волевого решения (выбором того, какие именно сингулярные значения считать нулевыми).

Это всего лишь стоящая предо мной подзадача, мне ещё предстоит выяснить даст ли направление, раскрывающееся её решением то, что мне нужно. Если дело выгорит, то можно будет задуматься в серьёз и об эффективности алгоритмов. Пока лишь хочу освоить общий случай и ознакомиться с тем, что на этом поприще уже проделано. Все предложенные выше книжки довольно старые. Неужели после 10-х годов ничего нового в этом направлении не придумывали и не публиковали? Или же это на столько хорошо проработанная тема, что уже и нечего?

[zgemm]

B@R5uk
Ээээ, с SVD вы совсем в другую сторону смотрите, чем принято.

Вам надо итерационно решать СЛАУ с регуляризацией, а именно

$$\min_x ||Ax - b||_2 + \lambda ||x||_2$$

и

$$\min_{x_t} ||A^T x_t - b_t||_2 + \lambda ||x_t||_2$$

что в общем-то приводит к классическим крыловским методам, и в общем я бы начал с сопряженного градиента, который в Вашем случае сойдется (при ненулевой лямбде) за $n^2$ арифметических операций.

Куча всего есть в спарс-паках, cuSparse, но в общем-то сопряженный градиент спокойно можно (в вызовах BLAS) в 20 строк запрограммировать.

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

-- 01.12.2025, 20:45 --

Если у вас действительно в каждой строчке и столбце в среднем 1-2 ненулевых элемента, еще можно смотреть в сторону разреженного и точного LU с минимальным заполнением, причем так, чтобы пока есть куда разлагать, то у вас есть эти LU факторы, а как кончится - у вас там дополнение по Шуру. Его уже или итерационно (если большое) или через SVD, если маленькое.

[vpb]

А вот есть такая книжка А.Н.Малышев, Введение в вычислительную линейную алгебру, так там третья глава называется "Сингулярное разложение и решение систем линейных уравнений". Как я себе представяляю, Малышев --- это студенческий аналог книжки Годунов, Антонов, Кирилюк, Костин, Гарантированная точность решения линейных систем, но та, хотя математически строгая, не очень читабельная.

З.Ы. Пардон, не точно привел название. "Гарантированная точность решения систем линейных уравнений в евклидовых пространствах".

[B@R5uk]

vpb, спасибо. Качну, чтобы было, читну. Но это тоже довольно старенькие книжки.

[B@R5uk]

А где можно поподробней почитать про алгоритмы, выполняющие rank-1 update для SVD (вместе с матрицами левых и правых векторов)? Что-то не помню, если это у уже предложенных книжках.

[Schrodinger's cat]

Всем добрый день.
Если еще актуально. Пару лет назад занимался данной проблемой.
1 Деммель мне вообще никак не помог. Написанное там не воспринималось моим мозгом. (это мой личный опыт, у вас может быть иначе)
2 Достаточно просто понять "простой итерационный" алгоритм.
https://rutube.ru/video/868c6edaf1de6356ed50a681ccd0962d/
У автора небольшая помарка, я в комментах на ютубе об этом ему сказал, исправил он видео или нет не знаю. Догадаться не сложно.
Работает очень долго, но вычисляет ))
3 Вращениями Якоби. Зануляются под диагональные элементы. Реализация легко гуглится. Об этом уже сказали.
Работает тоже долго, но вычисляет ))
4 Хаусхолдером с QR... уже не помню почему но реализацию в коде не довел до конца. Но сомневаюсь что по скорости это будет сильно лучше Якоби, чтобы побить скорость Intel MKL.

Главное что я для себя подчерпнул из этого опыта это то, что соревноваться с Intel MKL просто бесполезно. Мои реализации в 20 раз медленнее! Сомневаюсь что у вас получится хотя бы приблизиться к этой скорости если не использовать оптимизации по кэшам процессора и кучу прочих ухищрений.