Вычисление 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]

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