QR-алгоритм

QR-алгоритм — це чисельний метод у лінійній алгебрі, призначений для розв'язування повної задачі власних значень, тобто відшукання всіх власних значень і власних векторів матриці. Розробили в кінці 1950-х років незалежно В. М. Кублановська і Джон Френсіс^[en].

Алгоритм

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

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

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 при $k\to \infty ,$ збігається за формою до клітинного правого трикутного вигляду, відповідного клітинам з однаковими за модулем власними значеннями.^[1]

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

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

Доказ для симетричної додатноозначеної матриці

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

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

Нехай

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

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

A=S\Lambda S^{T}.

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

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}.

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

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

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

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

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

маємо

A^{k}=S_{k}T_{k}.

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

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

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

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

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

S^{T}=LU,

тоді

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

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

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

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

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

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

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

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

Позначимо

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

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

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

S_{k}P_{k}\to S

.

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

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

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

Оскільки

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),

то

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

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

\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-алгоритм дозволяє розв'язати повну проблему власних значень для симетричної додатноозначеної матриці.

Реалізація QR-алгоритму

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

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

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

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

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

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