SVD-разложение двухдиагональной матрицы

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

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

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

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

Подпрограмма RMatrixBDSVD осуществляет SVD-разложение двухдиагональной матрицы, заданной массивами D и E (главная диагональ и побочная диагональ). В зависимости от параметра IsUpper матрица трактуется, как верхняя или нижняя двухдиагональная.

Алгоритм сохраняет полученные сингулярные значения в массиве D, замещая главную диагональ. Матрицы преобразования U и V^T выводятся следующим способом. В алгоритм передаются переменные U и VT. Переменная U содержит матрицу с NRU строками и N столбцами, и умножается справа на матрицу преобразования U. Переменная VT содержит матрицу с N строками и NCVT столбцами, и умножается слева на матрицу преобразования V^T.

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

Параметр IsFractionalAccuracyRequired контролирует погрешность, с которой находятся сингулярные значения. Если он равен True, то сингулярные значения находятся с одинаковой относительной погрешностью, иначе они находятся с одинаковой абсолютной погрешностью. Обычно имеет смысл устанавливать это параметр в False (что несколько увеличивает скорость работы), поскольку в большинстве практических применений повышенная точность нахождения малых сингулярных значений не нужна.