QR-алгоритм

Нехай 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.

Зауважимо, що

${\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},}$

тобто всі матриці A_k є подібними, тобто їхні власні значення рівні.

Нехай усі діагональні мінори матриці A не вироджені. Тоді послідовність матриць A_k при ${\displaystyle k\to \infty ,}$ збігається за формою до клітинного правого трикутного вигляду, відповідного клітинам з однаковими за модулем власними значеннями.[1]

Для того, щоб отримати власні вектори матриці, потрібно перемножити всі матриці Q_k.

Алгоритм вважається обчислювально стійким, оскільки проводиться ортогональними перетвореннями подібності.

Будемо вважати, що власні числа додатноозначеної матриці A впорядковано за спаданням:

${\displaystyle \lambda _{1}>\lambda _{2}>...>\lambda _{n}>0.}$

Нехай

${\displaystyle \Lambda =\mathrm {diag} \left(\lambda _{1},...,\lambda _{n}\right),}$

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

${\displaystyle A=S\Lambda S^{T}.}$

Знайдемо вираз для степенів початкової матриці через матриці Q_k і R_k. З одного боку, за визначенням QR-алгоритму:

${\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}.}$

Застосовуючи це співвідношення рекурсивно, отримаємо:

${\displaystyle A^{k}=Q_{1}\cdot ...\cdot Q_{k}\cdot R_{k}\cdot ...\cdot R_{1}}$

Увівши такі позначення:

${\displaystyle S_{k}=Q_{1}\cdot ...\cdot Q_{k},}$

${\displaystyle T_{k}=R_{k}\cdot ...\cdot R_{1},}$

маємо

${\displaystyle A^{k}=S_{k}T_{k}.}$

З іншого боку:

${\displaystyle A^{k}=S\Lambda ^{k}S^{T}.}$

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

${\displaystyle S\Lambda ^{k}S^{T}=S_{k}T_{k}.}$

Припустимо, що існує LU-розклад матриці S^T:

${\displaystyle S^{T}=LU,}$

тоді

${\displaystyle S\Lambda ^{k}LU=S_{k}T_{k}.}$

Помножимо праворуч на обернену до U матрицю, а потім — на обернену до Λ^k:

${\displaystyle S\Lambda ^{k}L=S_{k}T_{k}U^{-1},}$

${\displaystyle S\Lambda ^{k}L\Lambda ^{-k}=S_{k}T_{k}U^{-1}\Lambda ^{-k}.}$

Можна показати, що

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

При ${\displaystyle k\to \infty }$ без обмеження загальності можна вважати, що на діагоналі матриці L стоять одиниці, тому

${\displaystyle S_{k}T_{k}U^{-1}\Lambda ^{-k}\to S.}$

Позначимо

${\displaystyle P_{k}=T_{k}U^{-1}\Lambda ^{-k},}$

причому матриця P_k є верхньотрикутною, як добуток верхньотрикутних і діагональних матриць.

Таким чином, ми довели, що

${\displaystyle S_{k}P_{k}\to S}$.

З єдиності QR-розкладу випливає, що якщо добуток ортогональної і трикутної матриці збігається до ортогональної матриці, то трикутна матриця збігається до одиничної матриці. Зі сказаного випливає, що

${\displaystyle S_{k}=Q_{1}\cdot ...\cdot Q_{k}\to S.}$

Тобто матриці S_k збігаються до матриці власних векторів матриці A.

Оскільки

${\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),}$

то

${\displaystyle A_{k+1}=S_{k}^{T}AS_{k}.}$

Переходячи до границі, отримаємо:

${\displaystyle \lim _{k\to \infty }A_{k}=\lim _{k\to \infty }A_{k+1}=S^{T}AS=S^{T}S\Lambda S^{T}S=\Lambda .}$

Отже, ми довели, що QR-алгоритм дозволяє розв'язати повну проблему власних значень для симетричної додатноозначеної матриці.

За певних умов послідовність матриць ${\displaystyle A_{k}}$ збігається до трикутної матриці, розкладу Шура матриці ${\displaystyle A}$. В цьому випадку власні числа трикутної матриці розташовуються на її діагоналі, і задача знаходження власних чисел вважається розв'язаною. У тестах на збіжність не практично вимагати точних нулів у нульовій частині матриці, але можна скористатися теоремою Гершгоріна про кола, що задає межі помилок.

У початковому стані матриці (без додаткових перетворень) вартість ітерацій відносно висока. Вартість алгоритму можна зменшити, спочатку звівши матрицю ${\displaystyle A}$ до форми верхньої матриці Гессенберга (вартість отримання якої методом, заснованим на перетворенні Хаусхолдера, оцінюється як ${\displaystyle {10 \over 3}n^{3}+O(n^{2})}$ арифметичних операцій), і потім скориставшись скінченною послідовністю ортогональних перетворень подібності. Цей алгоритм нагадує двобічну QR-декомпозицію. (У звичайній QR-декомпозиції матриця відображень Хаусхолдера множиться на початкову тільки зліва, тоді як у разі використання форми Хессенберга матриця відображень множиться на початкову і зліва, і справа.) Знаходження QR-декомпозиції верхньої матриці Хессенберга оцінюється як ${\displaystyle 6n^{2}+O(n)}$ арифметичних операцій. Оскільки форма матриці Хессенберга майже верхньотрикутна (у ній тільки один піддіагональний елемент не дорівнює нулю), вдається зразу знизити число ітерацій необхідних для збігання QR-алгоритму.

У випадку, якщо початкова матриця симетрична, верхня матриця Хессенберга також симетрична і тому є тридіагональною. Цю ж властивість має вся послідовність матриць ${\displaystyle A_{k}}$. У цьому випадку вартість процедури оцінюється як ${\displaystyle {4 \over 3}n^{3}+O(n^{2})}$ арифметичних операцій з використанням методу, заснованого на перетворенні Хаусхолдера. Знаходження QR-розкладу симетричної тридіагональної матриці оцінюється як ${\displaystyle O(n)}$ операцій.

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

У сучасній обчислювальній практиці QR-алгоритм реалізують із використанням його неявної версії, що спрощує додавання множинних «зсувів». Початково матриця зводиться до форми верхньої матриці Хессенберга ${\displaystyle A_{0}=QAQ^{T}}$, так само як і в явній версії. Потім, на кожному кроці перша колонка ${\displaystyle A_{k}}$ перетворюється через малорозмірне перетворення подоби Хаусхолдера до першої колонки ${\displaystyle p(A_{k})}$ (або ${\displaystyle p(A_{k})e_{1}}$), де ${\displaystyle p(A_{k})}$ — це поліном степеня ${\displaystyle r}$, який визначає стратегію «зсувів» (зазвичай ${\displaystyle p(x)=(x-\lambda )(x-{\bar {\lambda }})}$, де ${\displaystyle \lambda }$ і ${\displaystyle {\bar {\lambda }}}$ — це два власних числа залишкової підматриці ${\displaystyle A_{k}}$ розміру 2х2, це так званий неявний подвійний зсув). Потім виконують послідовні перетворення Хаусхолдера розмірності ${\displaystyle r+1}$ з метою повернути робочу матрицю ${\displaystyle A_{k}}$ до форми верхньої матриці Хессенберга.

• Нотатки Пітера Олвера про ортогональні базиси і QR-алгоритм(англ.)