SVD-разложение прямоугольной матрицы

Сингулярным разложением матрицы A размером MxN называется её представление в виде A = U W V^T, где U - ортогональная матрица размером MxM, V - ортогональная матрица размером NxN, W - матрица размером MxN, на главной диагонали которой находятся неотрицательные числа, расположенные в порядке убывания, а все внедиагональные элементы равны нолю.

С учетом свойств матрицы W, большей частью состоящей из нулей, для получения матрицы A требуется не M столбцов матрицы U, а лишь первые min(M,N) столбцов (в примере выше - три столбца), аналогично, лишь первые min(M,N) строк матрицы V^T влияют на результат произведения. Эти столбцы и строки называются левыми и правыми сингулярными векторами.

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

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

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

Алгоритм

Выше был описан принцип работы современных алгоритмов сингулярного разложения: приведение матрицы к двухдиагональной форме с последующей её диагонализацией QR-алгоритмом. Эта простая схема вполне работоспособна, однако в неё можно внести дополнение, заметно повышающее быстродействие программы. Схема улучшенного алгоритма, описанная ниже, практически полностью позаимствована из пакета LAPACK (подпрограмма xGESVD).

Для квадратных и близких к ним по форме матриц описанный выше алгоритм работает достаточно хорошо. В случае, если входная матрица прямоугольная и имеет сильно вытянутую форму (такие матрицы часто встречаются в задачах статистики), то вместо приведения её к бидиагональной форме можно использовать LQ или QR разложение (в зависимости от того, какая из сторон матрицы больше), чтобы представить её в виде произведения прямоугольной ортогональной матрицы Q и квадратной (верхне- или нижнетреугольной) матрицы, после чего бидиагонализировать квадратную матрицу заметно меньшего размера, чем исходная матрица. В зависимости от соотношения числа строк и столбцов возможно увеличение скорости в 2-6 раз по сравнению с вариантом, в котором бидиагонализации подвергается исходная прямоугольная матрица.

Вторым улучшением, носящим более технический характер, является оптимизация вычисления матрицы U. При вычислении матрицы U осуществляются операции со столбцами матрицы, которые, при принятом во многих языках программирования способе хранения матриц по строкам, оказываются неоптимальными в плане использования кэша процессора. В случае, если доступна дополнительная память, можно вычислять матрицу U^T, осуществляя операции со строками матрицы, что намного эффективнее, особенно в случае большой матрицы U, после чего транспонировать полученный результат.

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

Алгоритм, описанный выше, очень хорош, намного лучше своих предшественников. Однако, за все надо платить - вместе со всеми используемыми подпрограммами алгоритм занимает около трех тысяч строк. С одной стороны, я не в восторге от таких размеров, с другой - этот алгоритм объективно лучше. Старый алгоритм в стиле LINPACK, который был здесь раньше, теперь находится в разделе устаревших алгоритмов. Если компактность реализации критична, можете использовать его, но я не рекомендую это делать - новый алгоритм надежнее и быстрее (как алгоритмически, так и за счет использования векторных операций).

Описание подпрограммы

Подпрограмма RMatrixSVD осуществляет SVD-разложение прямоугольной матрицы размером MxN. На выходе подпрограмма возвращает массив сингулярных значений, упорядоченных по убыванию, и, по требованию, матрицы U и V^T, причем возможно как возвращение только левых и правых сингулярных векторов, так и полных матриц размером MxM и NxN (в зависимости от параметров UNeeded и VTNeeded). Обратите внимание, что в переменной VT алгоритм возвращает не матрицу V, а матрицу V^T.