QR-алгоритм

QR-алгоритм — это численный метод в линейной алгебре, предназначенный для решения полной проблемы собственных значений, то есть отыскания всех собственных чисел и собственных векторов матрицы. Был разработан в конце 1950-х годов независимо В. Н. Кублановской и Дж. Фрэнсисом^[англ.].

Алгоритм

Пусть A — вещественная матрица, для которой мы хотим найти собственные числа и векторы. Положим A₀=A. На k-м шаге (начиная с k = 0) вычислим QR-разложение A_k=Q_kR_k, где Q_k — ортогональная матрица (то есть Q_k^T = Q_k⁻¹), а R_k — верхняя треугольная матрица. Затем мы определяем A_k+1 = R_kQ_k.

Заметим, что

[math]\displaystyle{ 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]

то есть все матрицы A_k являются подобными, то есть их собственные значения равны.

Пусть все диагональные миноры матрицы A не вырождены. Тогда последовательность матриц A_k при [math]\displaystyle{ k \to \infty, }[/math] сходится по форме к клеточному правому треугольному виду, соответствующему клеткам с одинаковыми по модулю собственными значениями.^[1]

Для того, чтобы получить собственные векторы матрицы, нужно перемножить все матрицы Q_k.

Алгоритм считается вычислительно устойчивым, т. к. производится ортогональными преобразованиями подобия.

Доказательство для симметричной положительно определённой матрицы

Будем считать, что собственные числа положительно-определённой матрицы A упорядочены по убыванию:

[math]\displaystyle{ \lambda_1 \gt \lambda_2 \gt ... \gt \lambda_n \gt 0. }[/math]

Пусть

[math]\displaystyle{ \Lambda=\mathrm{diag}\left(\lambda_{1},...,\lambda_{n}\right), }[/math]

а S — матрица, составленная из собственных векторов матрицы A. Тогда, матрица A может быть записана в виде спектрального разложения

[math]\displaystyle{ A=S\Lambda S^{T}. }[/math]

Найдём выражение для степеней исходной матрицы через матрицы Q_k и R_k. С одной стороны, по определению QR-алгоритма:

[math]\displaystyle{ 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]\displaystyle{ A^{k}=Q_{1}\cdot...\cdot Q_{k}\cdot R_{k}\cdot...\cdot R_{1} }[/math]

Введя следующие обозначения:

[math]\displaystyle{ S_{k}=Q_{1}\cdot...\cdot Q_{k}, }[/math]

[math]\displaystyle{ T_{k}=R_{k}\cdot...\cdot R_{1}, }[/math]

получим

[math]\displaystyle{ A^{k}=S_{k}T_{k}. }[/math]

С другой стороны:

[math]\displaystyle{ A^{k}=S\Lambda^{k}S^{T}. }[/math]

Приравнивая правые части последних двух формул, получим:

[math]\displaystyle{ S\Lambda^{k}S^{T}=S_{k}T_{k}. }[/math]

Предположим, что существует LU-разложение матрицы S^T:

[math]\displaystyle{ S^{T}=LU, }[/math]

тогда

[math]\displaystyle{ S\Lambda^{k}LU=S_{k}T_{k}. }[/math]

Умножим справа на обратную к U матрицу, а затем — на обратную к Λ^k:

[math]\displaystyle{ S\Lambda^{k}L=S_{k}T_{k}U^{-1}, }[/math]

[math]\displaystyle{ S\Lambda^{k}L\Lambda^{-k}=S_{k}T_{k}U^{-1}\Lambda^{-k}. }[/math]

Можно показать, что

[math]\displaystyle{ \Lambda^{k}L\Lambda^{-k}\to \mathrm{diag}\left(l_{11},...,l_{nn}\right)=L^{\prime}. }[/math]

При [math]\displaystyle{ k \to \infty }[/math] без ограничения общности можно считать, что на диагонали матрицы L стоят единицы, поэтому

[math]\displaystyle{ S_{k}T_{k}U^{-1}\Lambda^{-k}\to S. }[/math]

Обозначим

[math]\displaystyle{ P_{k}=T_{k}U^{-1}\Lambda^{-k}, }[/math]

причём матрица P_k является верхнетреугольной, как произведение верхнетреугольных и диагональных матриц.

Таким образом, мы доказали, что

[math]\displaystyle{ S_{k}P_{k}\to S }[/math].

Из единственности QR-разложения следует, что если произведение ортогональной и треугольной матрицы сходится к ортогональной матрице, то треугольная матрица сходится к единичной матрице. Из сказанного следует, что

[math]\displaystyle{ S_{k}=Q_{1}\cdot...\cdot Q_{k}\to S. }[/math]

То есть матрицы S_k сходятся к матрице собственных векторов матрицы A.

Так как

[math]\displaystyle{ 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]\displaystyle{ A_{k+1}=S_{k}^T A S_{k}. }[/math]

Переходя к пределу, получим:

[math]\displaystyle{ \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]\displaystyle{ A_k }[/math] сходится к треугольной матрице, разложению Шура матрицы [math]\displaystyle{ A }[/math]. В этом случае собственные числа треугольной матрицы располагаются на ее диагонали, и задача нахождения собственных чисел считается решенной. В тестах на сходимость не практично требовать точных нулей в нулевой части матрицы, но можно воспользоваться теоремой о кругах Гершгорина, задающей пределы ошибок.

В исходном состоянии матрицы (без дополнительных преобразований) стоимость итераций относительно высока. Стоимость алгоритма можно уменьшить, сначала приведя матрицу [math]\displaystyle{ A }[/math] к форме верхней матрицы Хессенберга (стоимость получения которой методом, основанном на преобразовании Хаусхолдера, оценивается как [math]\displaystyle{ {10 \over 3}n^3 + O(n^2) }[/math] арифметических операций), и затем используя конечную последовательность ортогональных преобразований подобия. Данный алгоритм чем-то похож на двухсторонюю QR-декомпозицию. (В обычной QR-декомпозиции матрица отражений Хаусхолдера умножается на исходную только слева, тогда как в случае использования формы Хессенберга матрица отражений умножается на исходную и слева, и справа.) Нахождение QR-декомпозиции верхней матрицы Хессенберга оценивается как [math]\displaystyle{ 6n^2+O(n) }[/math] арифметических операций. Из-за того, что форма матрицы Хессенберга почти верхнетреугольная (у неё только один поддиагональный элемент не равен нулю), удается сразу снизить число итераций требуемых для схождения QR-алгоритма.

В случае, если исходная матрица симметричная, верхняя матрица Хессенберга также симметричная и поэтому является трехдиагональной. Этим же свойством обладает вся последовательность матриц [math]\displaystyle{ A_k }[/math]. В этом случае стоимость процедуры оценивается как [math]\displaystyle{ {4 \over 3}n^3+O(n^2) }[/math] арифметических операций с использованием метода, основанного на преобразовании Хаусхолдера. Нахождение QR-разложения симметричной трехдиагональной матрицы оценивается как [math]\displaystyle{ O(n) }[/math] операций.

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

Реализация QR-алгоритма в неявном виде

В современной вычислительной практике QR-алгоритм реализуется с использованием его неявной версии, что упрощает добавление множественных «сдвигов». Исходно матрица приводится к форме верхней матрицы Хессенберга [math]\displaystyle{ A_0=QAQ^T }[/math], так же как и в явной версии. Затем, на каждом шаге первая колонка [math]\displaystyle{ A_k }[/math] преобразуется через малоразмерное преобразование подобия Хаусхолдера к первой колонке [math]\displaystyle{ p(A_k) }[/math] (или [math]\displaystyle{ p(A_k)e_1 }[/math]), где [math]\displaystyle{ p(A_k) }[/math] — это полином степени [math]\displaystyle{ r }[/math], который определяет стратегию «сдвигов» (обычно [math]\displaystyle{ p(x)=(x-\lambda)(x-\bar{\lambda}) }[/math], где [math]\displaystyle{ \lambda }[/math] и [math]\displaystyle{ \bar{\lambda} }[/math] — это два собственных числа остаточной подматрицы [math]\displaystyle{ A_k }[/math] размера 2×2, это так называемый неявный двойной сдвиг). Затем последовательные преобразования Хаусхолдера размерности [math]\displaystyle{ r+1 }[/math] производятся с целью вернуть рабочую матрицу [math]\displaystyle{ A_k }[/math] к форме верхней матрицы Хессенберга.

Примечания

↑ Численные методы / Н. С. Бахвалов, Н. П. Жидков, Г. М. Кобельков. — 3-е изд. — М.: БИНОМ, Лаборатория знаний, 2004. — С. 321. — 636 с. — ISBN 5-94774-175-X.

Ссылки

Заметки Питера Олвера о ортогональных базисах и QR-алгоритме (англ.)

[1] Численные методы / Н. С. Бахвалов, Н. П. Жидков, Г. М. Кобельков. — 3-е изд. — М.: БИНОМ, Лаборатория знаний, 2004. — С. 321. — 636 с. — ISBN 5-94774-175-X.

[1]