Обзор существующих алгоритмов

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

В основе наиболее популярных современных алгоритмов сингулярного разложения лежит приведение матрицы ортогональным преобразованием к двухдиагональной форме (достаточно простая задача, решающаяся за конечное число операций) и последующая её диагонализация итеративным QR-алгоритмом. Разные алгоритмы обычно отличаются тем, как они осуществляют итерации QR-алгоритма. Первоначально широкое распространение получило семейство алгоритмов, основанных на алгоритме Голуба-Кахана-Рейнча ("Golub Kahan Reinsch algorithm"), именно этот метод реализован в LINPACK/EISPACK. Достоинствами этого метода являются его простота и компактность реализации - исходный код занимает лишь 300 строк. Но это единственные его достоинства. Хотя обычно алгоритм успешно справляется с задачей, в некоторых сложных случаях его сходимость и точность оставляют желать лучшего.

В библиотеке LAPACK реализовано два алгоритма сингулярного разложения, исправляющих недостатки алгоритма из LINPACK. Первый алгоритм (подпрограммы xBDSQR, xGESVD), исходный код которого лег в основу алгоритма, представленного на этой странице, характеризуется более высокой точностью, чем его аналог из LINPACK, была значительно улучшена сходимость, благодаря чему новый алгоритм вытеснил своего предшественника.

Следует отметить такую важную особенность нового алгоритма, как повышенную точность нахождения малых сингулярных значений двухдиагональной матрицы. Реализованный в LINPACK вариант QR-итерации находил все сингулярные значения с одинаковой абсолютной погрешностью, при этом самое большое сингулярное значение находилось с точностью, близкой к машинной. Для решения систем линейных уравнений, обращения матриц и других задач этого вполне достаточно - в этих задачах важна именно абсолютная точность сингулярного разложения. Однако, в некоторых задачах наиболее интересны малые сингулярные значения, относительная погрешность которых слишком велика. Например, если первое сингулярное значение равно 1 и найдено с абсолютной погрешностью 10 ^-15 (около 15 верных цифр), то сингулярное значение, равное 10 ^-10 и найденное с той же абсолютной погрешностью, будет иметь лишь 5 верных цифр. Новый алгоритм находит все сингулярные значения с одинаковой относительной погрешностью, так что в данном примере оба значения будут иметь 15 верных значащих цифр.

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

Как выше упоминалось, в пакете LAPACK реализовано два алгоритма сингулярного разложения. Второй алгоритм, в англоязычной литературе называющийся "divide-and-conquer", разбивает задачу SVD-разложения большой двухдиагональной матрицы на несколько задач меньшего размера, к которым применяется обычный QR-алгоритм. Этот алгоритм характеризуется более высокой скоростью работы для больших матриц, чем обычный QR-алгоритм. В частности, SVD-разложение квадратной матрицы алгоритмом "divide-and-conquer" при N=100 в 2-4 раза быстрее, чем разложение обычным алгоритмом (включая время, требуемое на приведение матрицы к двухдиагональной форме), а при N=1000 в 6-7 раз быстрее. Однако, этот алгоритм технически настолько сложен, что в ближайшее время вряд ли будет включен в ALGLIB. Если для вашей задачи критична скорость SVD-разложения больших матриц, то имеет смысл обратиться к LAPACK.