Алгоритм

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

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

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

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

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