Русская Википедия:QR-алгоритм
QR-алгоритм — это численный метод в линейной алгебре, предназначенный для решения полной проблемы собственных значений, то есть отыскания всех собственных чисел и собственных векторов матрицы. Был разработан в конце 1950-х годов независимо В. Н. Кублановской и Шаблон:Нп3.
Алгоритм
Пусть A — вещественная матрица, для которой мы хотим найти собственные числа и векторы. Положим A0=A. На k-м шаге (начиная с k = 0) вычислим QR-разложение Ak=QkRk, где Qk — ортогональная матрица (то есть QkT = Qk−1), а Rk — верхняя треугольная матрица. Затем мы определяем Ak+1 = RkQk.
Заметим, что
- <math> A_{k+1} = R_k Q_k = Q_k^{-1} Q_k R_k Q_k = Q_k^{-1} A_k Q_k = Q_k^{T} A_k Q_k, </math>
то есть все матрицы Ak являются подобными, то есть их собственные значения равны.
Пусть все диагональные миноры матрицы A не вырождены. Тогда последовательность матриц Ak при <math> k \to \infty,</math> сходится по форме к клеточному правому треугольному виду, соответствующему клеткам с одинаковыми по модулю собственными значениями.[1]
Для того, чтобы получить собственные векторы матрицы, нужно перемножить все матрицы Qk.
Алгоритм считается вычислительно устойчивым, т. к. производится ортогональными преобразованиями подобия.
Доказательство для симметричной положительно определённой матрицы
Будем считать, что собственные числа положительно-определённой матрицы A упорядочены по убыванию:
- <math> \lambda_1 > \lambda_2 > ... > \lambda_n > 0. </math>
Пусть
- <math> \Lambda=\mathrm{diag}\left(\lambda_{1},...,\lambda_{n}\right),</math>
а S — матрица, составленная из собственных векторов матрицы A. Тогда, матрица A может быть записана в виде спектрального разложения
- <math>A=S\Lambda S^{T}. </math>
Найдём выражение для степеней исходной матрицы через матрицы Qk и Rk. С одной стороны, по определению QR-алгоритма:
- <math> A^{k}=A_{1}^{k}=\left(Q_{1}R_{1}\right)^{k}=Q_{1}\left(R_{1}Q_{1}\right)^{k-1}R_{1}=Q_{1}A_{2}^{k-1}R_{1}. </math>
Применяя это соотношение рекурсивно, получим:
- <math> A^{k}=Q_{1}\cdot...\cdot Q_{k}\cdot R_{k}\cdot...\cdot R_{1} </math>
Введя следующие обозначения:
- <math> S_{k}=Q_{1}\cdot...\cdot Q_{k}, </math>
- <math> T_{k}=R_{k}\cdot...\cdot R_{1}, </math>
получим
- <math> A^{k}=S_{k}T_{k}. </math>
С другой стороны:
- <math> A^{k}=S\Lambda^{k}S^{T}. </math>
Приравнивая правые части последних двух формул, получим:
- <math>S\Lambda^{k}S^{T}=S_{k}T_{k}.</math>
Предположим, что существует LU-разложение матрицы ST:
- <math>S^{T}=LU,</math>
тогда
- <math>S\Lambda^{k}LU=S_{k}T_{k}.</math>
Умножим справа на обратную к U матрицу, а затем — на обратную к Λk:
- <math>S\Lambda^{k}L=S_{k}T_{k}U^{-1},</math>
- <math>S\Lambda^{k}L\Lambda^{-k}=S_{k}T_{k}U^{-1}\Lambda^{-k}.</math>
Можно показать, что
- <math>\Lambda^{k}L\Lambda^{-k}\to \mathrm{diag}\left(l_{11},...,l_{nn}\right)=L^{\prime}.</math>
При <math> k \to \infty</math> без ограничения общности можно считать, что на диагонали матрицы L стоят единицы, поэтому
- <math> S_{k}T_{k}U^{-1}\Lambda^{-k}\to S. </math>
Обозначим
- <math>P_{k}=T_{k}U^{-1}\Lambda^{-k},</math>
причём матрица Pk является верхнетреугольной, как произведение верхнетреугольных и диагональных матриц.
Таким образом, мы доказали, что
- <math>S_{k}P_{k}\to S</math>.
Из единственности QR-разложения следует, что если произведение ортогональной и треугольной матрицы сходится к ортогональной матрице, то треугольная матрица сходится к единичной матрице. Из сказанного следует, что
- <math>S_{k}=Q_{1}\cdot...\cdot Q_{k}\to S.</math>
То есть матрицы Sk сходятся к матрице собственных векторов матрицы A.
Так как
- <math> A_{k+1}=Q_{k}^{T}A_{k}Q_{k}=...=\left(Q_{k}^{T}\cdot...\cdot Q_{1}^{T}\right)A_{1}\left(Q_{1}\cdot...\cdot Q_{k}\right)=\left(Q_{1}\cdot...\cdot Q_{k}\right)^{T}A\left(Q_{1}\cdot...\cdot Q_{k}\right),</math>
то
- <math> A_{k+1}=S_{k}^T A S_{k}.</math>
Переходя к пределу, получим:
- <math> \lim_{k\to\infty} A_{k}=\lim_{k\to\infty} A_{k+1} = S^T A S = S^T S \Lambda S^T S = \Lambda. </math>
Итак, мы доказали, что QR-алгоритм позволяет решить полную проблему собственных значений для симметричной положительно-определённой матрицы.
Реализация QR-алгоритма
При определенных условиях последовательность матриц <math>A_k</math> сходится к треугольной матрице, разложению Шура матрицы <math>A</math>. В этом случае собственные числа треугольной матрицы располагаются на ее диагонали, и задача нахождения собственных чисел считается решенной. В тестах на сходимость не практично требовать точных нулей в нулевой части матрицы, но можно воспользоваться теоремой о кругах Гершгорина, задающей пределы ошибок.
В исходном состоянии матрицы (без дополнительных преобразований) стоимость итераций относительно высока. Стоимость алгоритма можно уменьшить, сначала приведя матрицу <math>A</math> к форме верхней матрицы Хессенберга (стоимость получения которой методом, основанном на преобразовании Хаусхолдера, оценивается как <math>{10 \over 3}n^3 + O(n^2)</math> арифметических операций), и затем используя конечную последовательность ортогональных преобразований подобия. Данный алгоритм чем-то похож на двухсторонюю QR-декомпозицию. (В обычной QR-декомпозиции матрица отражений Хаусхолдера умножается на исходную только слева, тогда как в случае использования формы Хессенберга матрица отражений умножается на исходную и слева, и справа.) Нахождение QR-декомпозиции верхней матрицы Хессенберга оценивается как <math>6n^2+O(n)</math> арифметических операций. Из-за того, что форма матрицы Хессенберга почти верхнетреугольная (у неё только один поддиагональный элемент не равен нулю), удается сразу снизить число итераций требуемых для схождения QR-алгоритма.
В случае, если исходная матрица симметричная, верхняя матрица Хессенберга также симметричная и поэтому является трехдиагональной. Этим же свойством обладает вся последовательность матриц <math>A_k</math>. В этом случае стоимость процедуры оценивается как <math>{4 \over 3}n^3+O(n^2)</math> арифметических операций с использованием метода, основанного на преобразовании Хаусхолдера. Нахождение QR-разложения симметричной трехдиагональной матрицы оценивается как <math>O(n)</math> операций.
Скорость сходимости зависит от степени разделения собственных чисел, и в практической реализации в явном или неявном виде используюся «сдвиги» для усиления разделения собственных чисел и для ускорения сходимости. В типичном виде для симметричных матриц QR-алгоритм точно находит одно собственное число (уменьшая размерность матрицы) за одну или две итерации, что делает этот подход как эффективным, так и надежным.
Реализация QR-алгоритма в неявном виде
В современной вычислительной практике QR-алгоритм реализуется с использованием его неявной версии, что упрощает добавление множественных «сдвигов». Исходно матрица приводится к форме верхней матрицы Хессенберга <math>A_0=QAQ^T</math>, так же как и в явной версии. Затем, на каждом шаге первая колонка <math>A_k</math> преобразуется через малоразмерное преобразование подобия Хаусхолдера к первой колонке <math>p(A_k)</math> (или <math>p(A_k)e_1</math>), где <math>p(A_k)</math> — это полином степени <math>r</math>, который определяет стратегию «сдвигов» (обычно <math>p(x)=(x-\lambda)(x-\bar{\lambda})</math>, где <math>\lambda</math> и <math>\bar{\lambda}</math> — это два собственных числа остаточной подматрицы <math>A_k</math> размера 2×2, это так называемый неявный двойной сдвиг). Затем последовательные преобразования Хаусхолдера размерности <math>r+1</math> производятся с целью вернуть рабочую матрицу <math>A_k</math> к форме верхней матрицы Хессенберга.
Примечания
Ссылки
