Алгоритм обчислення власних значень

Однією з найважливіших задач обчислювальної математики є створення ефективних і стійких алгоритмів знаходження власних значень матриці. Ці алгоритми обчислення власних значень можуть також знаходити власні вектори.

Якщо задано n × n квадратну матрицю A над дійсними або комплексними числами, власне значення λ і відповідний йому кореневий вектор v - це пара, що задовольняє рівності^[1]

${\displaystyle \left(A-\lambda E\right)^{k}{\mathbf {v} }=0,}$

де v ненульовий n × 1 вектор-стовпець, E — n × n одинична матриця, k — додатне ціле, а λ і v можуть бути комплексними, навіть якщо A дійсне. Якщо k = 1, вектор просто називається власним вектором. У цьому випадку Av = λv. Будь-яке власне значення λ матриці A має простий^{[прим. 1]} власний вектор, відповідний йому так, що якщо k — найменше ціле, за якого (A - λE)^k v = 0 для кореневого вектора v, то (A - λE)^k-1 v буде простим власним вектором. Значення k завжди можна взяти меншим або рівним n. Зокрема, (A - λE)ⁿ v = 0 для всіх кореневих векторів v відповідних λ.

Для будь-якого власного значення λ матриці A ядро ker(A - λE) складається з усіх власних векторів, відповідних λ, (разом з 0) і називається власним підпростором числа λ, а векторний підпростір ker((A - λE)ⁿ) складається з усіх кореневих векторів (доповнений нульовим вектором) і називається кореневим підпростором. Геометрична кратність значення λ є розмірністю його власного підпростору. Алгебрична кратність значення λ є розмірністю його кореневого підпростору. Подальші терміни пов'язані з рівністю

${\displaystyle p_{A}\left(z\right)={\rm {det}}\left(zE-A\right)=\prod _{i=1}^{k}(z-\lambda _{i})^{\alpha _{i}},}$

Тут det — визначник, λ_i — всі різні власні значення матриці A, а α_i — відповідні алгебричні кратності. Функція p_A(z) — це характеристичний многочлен матриці A. Отже, алгебрична кратність є кратністю власних значень як коренів характеристичного многочлена. Оскільки будь-який власний вектор є кореневим вектором, геометрична кратність менша або дорівнює алгебричній кратності. Сума алгебричних кратностей дорівнює n степеню характеристичного многочлена. Рівняння p_A(z) = 0 називають характеристичним рівнянням, оскільки його корені є точно власними значеннями матриці A. В теоремі Гамільтона — Келі сама матриця A задовольняє тому самому рівнянню: p_A(A) = 0^{[прим. 2]}. Як наслідок, стовпці матриці ${\displaystyle \textstyle \prod _{i\neq j}(A-\lambda _{i}E)^{\alpha _{i}}}$ мають бути або 0, або кореневими векторами, відповідними власному значенню λ_j, оскільки вони знищуються матрицею ${\displaystyle \textstyle (A-\lambda _{j}E)^{\alpha _{j}}.}$

Будь-який набір кореневих векторів різних власних значень лінійно незалежний, так що базис для всього C ⁿ можна вибрати з набору кореневих векторів. Точніше цей базис {v_i}ⁿ_i=1 можна вибрати і вибудувати так, що

• якщо v_i і v_j мають одне і те саме власне значення, то це виконуватиметься для будь-якого v_k при k між i і j;
• якщо v_i не є простим власним вектором і якщо λ_i — його власне значення, то (A - λ_iE )v_i = v_i-1 (зокрема v₁ має бути простим власним вектором).

Якщо ці базисні вектори розташувати як стовпці матриці V = [ v₁ v₂ … v_n ], то V можна використовувати для перетворення матриці A в її нормальну жорданову форму:

${\displaystyle V^{-1}AV={\begin{bmatrix}\lambda _{1}&\beta _{1}&0&\ldots &0\\0&\lambda _{2}&\beta _{2}&\ldots &0\\0&0&\lambda _{3}&\ldots &0\\\vdots &\vdots &\vdots &\ddots &\vdots \\0&0&0&\ldots &\lambda _{n}\end{bmatrix}},}$

де λ_i — власне значення, β_i = 1 якщо (A - λ_i+1)v_i+1 = v_i і β_i = 0 в інших випадках.

Якщо W є оборотною матрицею і λ — власне значення матриці A з відповідним кореневим вектором v, то (W ^-1AW - λE )^k W ^-kv = 0. Отже, λ є власним значенням матриці W ^-1AW з відповідним кореневим вектором W ^-kv. Таким чином, подібні матриці мають однакові власні значення.

Ермітово-спряжена матриця M^* до комплексної матриці M — це траспонована матриця із заміною всіх елементів спряженими значеннями: M ^* = M ^T. Квадратну матрицю A називають нормальною, якщо вона комутує з ермітово-спряженою: A^*A = AA^*. Матрицю називають ермітовою, якщо вона дорівнює своїй спряженій: A^* = A. Всі ермітові матриці нормальні. Якщо A має тільки дійсні елементи, то спряжена до неї — це просто транспонована матриця і вона буде ермітовою в тому і тільки в тому випадку, коли вона симетрична. Якщо застосувати це до стовпців, спряженість можна використовувати для визначення канонічного скалярного добутку в C ⁿ: w • v = w^* v^{[прим. 3]}. Нормальні, ермітові і дійсні симетричні матриці мають низку корисних властивостей:

• Кожен кореневий власний вектор нормальної матриці є простим власним вектором.
• Будь-яка нормальна матриця подібна до діагональної, оскільки її нормальна жорданова форма є діагональною матрицею.
• Власні вектори, відповідні різним власним значенням нормальної матриці, ортогональні.
• Для будь-якої нормальної матриці A C ⁿ має ортонормальний базис, що складається з власних векторів матриці A. Відповідна матриця власних векторів є унітарною.
• Власні значення ермітової матриці є дійвсними числами, оскільки (λ — λ)v = (A^* — A)v = (A — A)v = 0 для ненульового власного вектора v.
• Якщо матриця A дійсна, існує ортонормальний базис для R ⁿ, що складається з власних векторів матриці A, в тому і тільки в тому випадку, коли A симетрична.

Як дійсні, так і комплексні матриці, які не є при цьому ермітовими матрицями, можуть мати всі власні значення дійсними. Наприклад, дійсна трикутна матриця має всі свої власні значення на діагоналі, але, в загальному випадку, не симетрична.

Будь-яку задачу обчислювальної математики можна розглядати як обчислення деякої функції ƒ від деякого аргументу x. Число обумовленості κ(ƒ, x) задачі — це відношення відносної похибки результату обчислення до відносної похибки параметра функції і залежить як від функції, так і від параметра. Число обумовленості описує, наскільки зростає похибка під час обчислень. Десятковий логарифм цього числа говорить про кількість знаків, які ми втрачаємо відносно вихідних даних. Число обумовленості стосується найкращого сценарію, відбиваючи нестабільність самої задачі, незалежно від способу розв'язування. Ніякий алгоритм не може дати результату кращого, ніж визначений числом обумовленості, хіба що випадково. Однак погано розроблений алгоритм може дати істотно гірші результати. Наприклад, як буде згадано нижче, задача знаходження власних значень нормальної матриці завжди добре обумовлена, проте задача знаходження коренів многочлена може бути дуже погано обумовленою^[en]. Такі алгоритми обчислення власних значень, які працюють шляхом знаходження коренів характеристичного многочлена, можуть виявитися погано обумовленими, навіть якщо сама задача добре обумовлена.

Для задачі розв'язування системи лінійних рівнянь Av = b, де A є оборотною, число обумовленості κ(A^-1, b) визначається виразом ||A||_op||A⁻¹||_op, де || ||_op — операторна норма, підпорядкована звичайній евклідовій нормі на C ⁿ. Оскільки це число не залежить від b і є тим самим як для A, так і для A^-1, його зазвичай називають числом обумовленості κ(A) матриці A. Це значення κ(A) є також абсолютним значенням відношення найбільшого власного значення матриці A до її найменшого власного значення. Якщо A є унітарною, то ||A||_op = ||A⁻¹||_op = 1, так що κ(A) = 1. У загальному випадку для матриць часто складно обчислити операторну норму. З цієї причини для оцінення числа обумовленості зазвичай використовують інші норми матриці.

Для задачі обчислення власних значень Бауер і Файк довели^[en], що якщо λ є власним значенням діагоналізовної n × n матриці A з матрицею власних векторів V, то абсолютна похибка обчислення λ обмежена добутом κ(V) і абсолютною похибкою в A: ${\displaystyle \;\;|\delta \lambda |\leqslant \kappa _{op}(V)\|\delta A\|_{op}}$^[2]. Як наслідок, число обумовленості для обчислення λ дорівнює κ(λ, A) = κ(V) = ||V ||_op ||V ⁻¹||_op. Якщо матриця A нормальна, то V є унітарною і κ(λ, A) = 1. Таким чином, задача обчислення власних значень нормальних матриць добре обумовлена.

Показано, що число обумовленості задачі обчислення власного підпростору нормальної матриці A, відповідного власного значення λ, обернено пропорційне мінімальній відстані між λ та іншими, відмінними від λ, власними значеннями матриці A^[3]. Зокрема, задача визначення власного підпростору для нормальних матриць добре обумовлена для ізольованих власних значень. Якщо власні значення не ізольовані, найкраще, на що можна розраховувати, це визначення підпростору всіх власних векторів довколишніх власних значень.

Будь-який нормований многочлен є характеристичним многочленом супутньої матриці, тому алгоритм для обчислення власних значень можна використовувати для знаходження коренів многочленів. Теорема Абеля — Руффіні показує, що будь-який такий алгоритм для розмірності більшої від 4 має або бути нескінченним, або залучати функції складніші, ніж елементарні арифметичні операції або дробові степені. З цієї причини алгоритми, які обчислюють точно власні значення за скінченне число кроків, існують тільки для спеціальних класів матриць. У загальному випадку алгоритми є ітеративними і дають на кожній ітерації чергове наближення до розв'язку.

Деякі алгоритми дають усі власні значення, інші дають кілька значень або навіть всього одне, однак і ці алгоритми можна використовувати для обчислення всіх власних значень. Як тільки власне значення λ матриці A визначено, його можна використовувати або для зведення алгоритму до отримання іншого власного значення, або для зведення задачі до такої, для якої λ не є розв'язком.

Зведення зазвичай виконують зсувом: A замінюється на A - μE для деякої константи μ. Щоб отримати власне значення матриці A, потрібно власне значення, знайдене для A - μE, додати до μ. Наприклад, у степеневому методі μ = λ. Ітерація степеневого методу знаходить найбільше за абсолютною величиною значення, так що навіть якщо λ є наближенням до власного значення, ітерація степеневого методу навряд чи знайде його вдруге. І навпаки, методи, засновані на зворотних ітераціях знаходять найменше власне значення, так що μ вибирається подалі від λ в надії опинитися ближче до якогось іншого власного значення.

Зведення можна виконати звуженням матриці A до простору стовпців матриці A - λE. Оскільки A - λE вироджена, простір стовпців має меншу розмірність. Алгоритм обчислення власних значень можна тоді застосувати до цієї звуженої матриці. Процес можна продовжувати, поки не буде знайдено всі власні значення.

Якщо алгоритм знаходження власних значень не дає власного вектора, поширеним є застосування алгоритму, заснованого на зворотній ітерації, з прирівнюванням μ до найближчої апроксимації власного значення. Це швидко приводить до власного вектора найближчого до μ власного значення. Для невеликих матриць альтернативою є використання стовпцевого підпростору добутку A - λ́E для кожного з решти власних значень λ́.

Оскільки власними значеннями трикутної матриці є діагональні елементи, в загальному випадку не існує скінченного методу, подібного до виключень Гауса, для зведення матриці до трикутної форми, зберігаючи при цьому власні значення, однак можна отримати щось близьке до трикутної матриці. Верхня матриця Гессенберга — це квадратна матриця, в якій усі елементи нижче від першої піддіагоналі дорівнюють нулю. Нижня матриця Гессенберга — це квадратна матриця, в якій усі члени вище від першої наддіагоналі дорівнюють нулю. Матриці, які є як нижніми, так і верхніми матрицями Гессенберга — це тридіагональні матриці. Матриці Гессенберга і тридіагональні матриці є початковими для багатьох алгоритмів обчислення власних значень, оскільки нульові значення зменшують складність задачі. Існує кілька методів зведення матриць до матриць Гессенберга з тими ж власними значеннями. Якщо початкова матриця симетрична або ермітова, то кінцева матриця буде тридіагональною.

Якщо потрібні тільки власні значення, не потрібно обчислювати матрицю подібності, оскільки перетворена матриця має ті самі власні значення. Якщо також потрібні і власні вектори, матриця подібності необхідна для перетворення власних векторів матриці Гессенберга на власні вектори початкової матриці.

Ітеративні алгоритми розв'язують задачу обчислення власних значень побудовою послідовностей, що збігаються до власних значень. Деякі алгоритми дають також послідовності векторів, що збігаються до власних векторів. Найчастіше послідовності власних значень виражаються через послідовності подібних матриць, які збігаються до трикутної або діагональної форми, що дозволяє потім просто отримати власні значення. Послідовності власних векторів виражаються через відповідні матриці подібності.

Не існує простих алгоритмів прямого обчислення власних значень для матриць у загальному випадку, але для багатьох спеціальних класів матриць власні значення можна обчислити прямо. Це:

Оскільки визначник трикутної матриці є добутком її діагональних елементів, то для трикутної матриці T ${\displaystyle \mathrm {det} \left(\lambda E-T\right)={\prod }_{i}\left(\lambda -T_{ii}\right)}$. Отже, власні значення T — це її діагональні елементи.

Якщо p — будь-який многочлен і p(A) = 0, то власні значення матриці A задовольняють тому ж рівнянню. Якщо вдається розкласти многочлен p на множники, то власні значення A містяться серед його коренів.

Наприклад, проєктор — це квадратна матриця P, що задовольняє рівнянню P² = P. Коренями відповідного скалярного поліноміального рівняння λ² = λ будуть 0 і 1. Таким чином, проєктор має 0 і 1 власними значеннями. Кратність власного значення 0 — це дефект P, тоді як кратність 1 — це ранг P.

Інший приклад — матриця A, що задовольняє рівнянню A² = α²E для деякого скаляра α. Власні значення мають бути рівними ±α. Оператори проєктування

${\displaystyle P_{+}={\frac {1}{2}}\left(E+{\frac {A}{\alpha }}\right)}$

${\displaystyle P_{-}={\frac {1}{2}}\left(E-{\frac {A}{\alpha }}\right)}$

задовольняють рівностям

${\displaystyle AP_{+}=\alpha P_{+}\quad AP_{-}=-\alpha P_{-}}$

і

${\displaystyle P_{+}P_{+}=P_{+}\quad P_{-}P_{-}=P_{-}\quad P_{+}P_{-}=P_{-}P_{+}=0.}$

Простори стовпців матриць P₊ і P_- є підпросторами власних векторів матриці A, які відповідають +α і -α, відповідно.

Для розмірностей від 2 до 4 існують формули з використанням радикалів, які можна використовувати для обчислення власних значень. Для матриць 2х2 і 3х3 це звичайна практика, але для матриць 4х4 через зростання складності формул коренів^[en] цей підхід менш привабливий.

Для матриць 2×2

${\displaystyle A={\begin{bmatrix}a&b\\c&d\end{bmatrix}},}$

характеристичний многочлен дорівнює

${\displaystyle {\rm {det}}{\begin{bmatrix}\lambda -a&-b\\-c&\lambda -d\end{bmatrix}}=\lambda ^{2}\,-\,\left(a+d\right)\lambda \,+\,\left(ad-bc\right)=\lambda ^{2}\,-\,\lambda \,{\rm {tr}}(A)\,+\,{\rm {det}}(A).}$

Власні значення можна знайти як корені квадратного рівняння:

${\displaystyle \lambda ={\frac {{\rm {tr}}(A)\pm {\sqrt {{\rm {tr}}^{2}(A)-4{\rm {det}}(A)}}}{2}}.}$

Якщо визначити ${\displaystyle \textstyle {\rm {gap}}\left(A\right)={\sqrt {{\rm {tr}}^{2}(A)-4{\rm {det}}(A)}}}$ як відстань між двома власними значеннями, легко обчислити

${\displaystyle {\frac {\partial \lambda }{\partial a}}={\frac {1}{2}}\left(1\pm {\frac {a-d}{{\rm {gap}}(A)}}\right),\qquad {\frac {\partial \lambda }{\partial b}}={\frac {\pm c}{{\rm {gap}}(A)}}}$

з подібними формулами для c і d. Звідси випливає, що обчислення добре обумовлене, якщо власні значення ізольовані.

Власні вектори можна отримати, використовуючи теорему Гамільтона — Келі. Якщо λ₁, λ₂ — власні значення, то (A - λ₁E )(A - λ₂E ) = (A - λ₂E )(A - λ₁E ) = 0, так що стовпці (A - λ₂E ) обнуляються матрицею (A - λ₁E ) і навпаки. Припускаючи, що жодна з матриць не дорівнює нулю, стовпці кожної матриці мають містити власні вектори для іншого власного значення (якщо ж матриця нульова, то A є добутком одиничної матриці на константу і будь-який ненульовий вектор є власним).

Наприклад, нехай

${\displaystyle A={\begin{bmatrix}4&3\\-2&-3\end{bmatrix}},}$

тоді tr(A) = 4 - 3 = 1 і det(A) = 4(-3) - 3(-2) = -6, так що характеристичне рівняння

${\displaystyle 0=\lambda ^{2}-\lambda -6=(\lambda -3)(\lambda +2),}$

а власні значення рівні 3 і -2. Тепер

${\displaystyle A-3E={\begin{bmatrix}1&3\\-2&-6\end{bmatrix}}}$, ${\displaystyle \qquad A+2E={\begin{bmatrix}6&3\\-2&-1\end{bmatrix}}}$

В обох матрицях стовпці відрізняються скалярними коефіцієнтами, так що можна вибирати будь-який стовпець. Так, (1, -2) можна використовувати як власний вектор, відповідного власному значенню -2, а (3, -1) — як власний вектор для власного значення 3, що легко можна перевірити множенням на матрицю A.

Якщо A — матриця 3 х 3, то характеристичним рівнянням буде:

${\displaystyle {\rm {det}}\left(\alpha E-A\right)=\alpha ^{3}-\alpha ^{2}{\rm {tr}}(A)-\alpha {\frac {1}{2}}\left({\rm {tr}}(A^{2})-{\rm {tr}}^{2}(A)\right)-{\rm {det}}(A)=0.}$

Це рівняння можна розв'язати за допомогою методів Кардано або Лагранжа, але афінне перетворення матриці A істотно спрощує вираз і веде прямо до тригонометричного розв'язку. Якщо A = pB + qE, то A і B мають однакові власні вектори і β є власним значенням матриці B тоді й лише тоді, коли α = pβ + q є власним значенням матриці A. Якщо покласти ${\displaystyle \textstyle q={\rm {tr}}(A)/3}$ і ${\displaystyle \textstyle p=\left({\rm {tr}}\left((A-qE)^{2}\right)/6\right)^{1/2}}$, одержимо

${\displaystyle {\rm {det}}\left(\beta E-B\right)=\beta ^{3}-3\beta -{\rm {det}}(B)=0.}$

Підстановка β = 2cos θ і деяке спрощення з використанням тотожності cos 3θ = 4cos³ θ - 3cos θ зводить рівняння до cos 3θ = det(B) / 2. Отже,

${\displaystyle \beta =2{\rm {cos}}\left({\frac {1}{3}}{\rm {arccos}}\left({\rm {det}}(B)/2\right)+{\frac {2k\pi }{3}}\right),\quad k=0,1,2.}$

Якщо det(B) є комплексним числом або більше 2 за абсолютною величиною, арккосинус для всіх трьох величин k слід брати з однієї гілки. Проблема не виникає, якщо A дійсна і симетрична, приводячи до простого алгоритму:^[12]

% Given a real symmetric 3x3 matrix A, compute the eigenvalues

p1 = A(1,2)^2 + A(1,3)^2 + A(2,3)^2
if (p1 == 0) 
  % A is diagonal.
  eig1 = A(1,1)
  eig2 = A(2,2)
  eig3 = A(3,3)
else
  q = trace(A)/3
  p2 = (A(1,1) - q)^2 + (A(2,2) - q)^2 + (A(3,3) - q)^2 + 2 * p1
  p = sqrt(p2 / 6)
  B = (1 / p) * (A - q * E)    % E is the identity matrix
  r = det(B) / 2

  % In exact arithmetic for a symmetric matrix -1 <= r <= 1
  % but computation error can leave it slightly outside this range.
  if (r <= -1) 
   phi = pi / 3
  elseif (r >= 1)
   phi = 0
  else
   phi = acos(r) / 3
  end

  % the eigenvalues satisfy eig3 <= eig2 <= eig1
  eig1 = q + 2 * p * cos(phi)
  eig3 = q + 2 * p * cos(phi + (2*pi/3))
  eig2 = 3 * q - eig1 - eig3   % since trace(A) = eig1 + eig2 + eig3
end

Власні вектори також A можна отримати, скориставшись теоремою Гамільтона — Келі. Якщо α₁, α₂, α₃ — різні власні значення матриці A, то (A - α₁E)(A - α₂E)(A - α₃E) = 0. Тоді стовпці добутку будь-яких двох з цих матриць містять власні вектори третього власного значення. Однак якщо a₃ = a₁, то (A - α₁E)²(A - α₂E) = 0 і (A - α₂E)(A - α₁E)² = 0. Таким чином, кореневий власний підпростір α₁ натягнуий на стовпці A - α₂E, тоді як звичайний власний підпростір натягнутий на стовпці (A - α₁E)(A - α₂E). Звичайний власний підпростір α₂ натягнутий на стовпці (A - α₁E)².

Наприклад, нехай

${\displaystyle A={\begin{bmatrix}3&2&6\\2&2&5\\-2&-1&-4\end{bmatrix}}.}$

Характеристичне рівняння буде

${\displaystyle 0=\lambda ^{3}-\lambda ^{2}-\lambda +1=(\lambda -1)^{2}(\lambda +1),}$

з власними значеннями 1 (кратності 2) і −1. Обчислюємо

${\displaystyle A-E={\begin{bmatrix}2&2&6\\2&1&5\\-2&-1&-5\end{bmatrix}},\qquad A+E={\begin{bmatrix}4&2&6\\2&3&5\\-2&-1&-3\end{bmatrix}}}$,

а потім

${\displaystyle (A-E)^{2}={\begin{bmatrix}-4&0&-8\\-4&0&-8\\4&0&8\end{bmatrix}},\qquad (A-E)(A+E)={\begin{bmatrix}0&4&4\\0&2&2\\0&-2&-2\end{bmatrix}}}$.

Тоді (-4, -4, 4) є власним вектором для −1, а (4, 2, -2) є власним вектором для 1. Вектори (2, 3, -1) і (6, 5, -3) є кореневими векторами, відповідними значенню 1, будь-який з яких можна скомбінувати з (-4, -4, 4) і (4, 2, -2), утворюючи базис кореневих векторів матриці A.

• Список тем числового аналізу^[en]

• Дж. Голуб, Ч. Ван Лоун. Матричные вычисления. — М. : «Мир», 1999. — ISBN 5-03-002406-9.
• Б. Парлетт. Симметричная проблема собственных значений. — М. : «Мир», 1983.
• Дж. Х. Уилкинсон. Алгебраическая проблема собственных значений. — М. : «Наука» Главная редакция физико-математической литературы, 1970.

• Adam W. Bojanczyk, Adam Lutoborski. Computation of the Euler angles of a symmetric 3X3 matrix // SIAM Journal on Matrix Analysis and Applications. — 1991. — Т. 12, вип. 1 (1 січня). — С. 41—48. — DOI:10.1137/0612005.