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

Алгоритм вычисления собственных значений — алгоритм, позволяющий определить собственные значения и собственные векторы заданной матрицы. Создание эффективных и устойчивых алгоритмов для этой задачи является одной из ключевых задач вычислительной математики.

Если задана 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 имеет простой^{[note 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^{[note 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}n
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^{[note 3]}. Нормальные, эрмитовы и вещественные симметричные матрицы имеют ряд полезных свойств:

• Каждый корневой собственный вектор нормальной матрицы является простым собственным вектором.
• Любая нормальная матрица подобна диагональной, поскольку её нормальная жорданова форма является диагональной матрицей.
• Собственные вектора, соответствующие различным собственным значениям нормальной матрицы, ортогональны.
• Для любой нормальной матрицы A C ⁿ имеет ортонормальный базис, состоящий из собственных векторов матрицы A. Соответствующая матрица собственных векторов является унитарной.
• Собственные значения эрмитовой матрицы являются вещественными числами, поскольку (λ - λ)v = (A^* - A)v = (A - A)v = 0 для ненулевого собственного вектора v.
• Если матрица A вещественна, существует ортонормальный базис для R ⁿ, состоящий из собственных векторов матрицы A, в том и только в том случае, когда A симметрична.

Возможно как для вещественных, так и для комплексных матриц иметь все собственные значения вещественными, не будучи при этом эрмитовой матрицей. Например, вещественная треугольная матрица имеет все свои собственные значения на диагонали, но, в общем случае, не симметрична.

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

Для задачи решения системы линейных уравнений Av = b, где A является обратимой, число обусловленности κ(A^-1, b) определяется выражением ||A||_op||A^-1||_op, где || ||_op — операторная норма, подчинённая обычной евклидовой норме на C ⁿ. Поскольку это число не зависит от b и является тем же самым как для A, так и для A^-1, оно обычно называется числом обусловленности κ(A) матрицы A. Это значение κ(A) является также абсолютным значением отношения наибольшего собственного значения матрицы A к её наименьшему собственному значению. Если A является унитарной, то ||A||_op = ||A^-1||_op = 1, так что κ(A) = 1. В общем случае для матриц часто сложно вычислить операторную норму. По этой причине обычно используют другие нормы матрицы для оценки числа обусловленности.

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

Было показано, что число обусловленности задачи вычисления собственного подпространства нормальной матрицы A, соответствующего собственному значению λ, обратно пропорционально минимальному расстоянию между λ и другими, отличными от λ, собственными значениями матрицы A^[3]. В частности, задача определения собственного подпространства для нормальных матриц хорошо обусловлена для изолированных собственных значений. Если собственные значения не изолированы, лучшее, на что мы можем рассчитывать, это определение подпространства всех собственных векторов близлежащих собственных значений.

Любой нормированный многочлен является характеристическим многочленом сопровождающей матрицы, поэтому алгоритм для вычисления собственных значений можно использовать для нахождения корней многочленов. Теорема Абеля — Руффини показывает, что любой такой алгоритм для размерности большей 4 должен либо быть бесконечным, либо вовлекать функции более сложные, чем элементарные арифметические операции или дробные степени. По этой причине алгоритмы, вычисляющие точно собственные значения за конечное число шагов, существуют только для специальных классов матриц. В общем случае алгоритмы являются итеративными, дающими на каждой итерации очередное приближение к решению.

Некоторые алгоритмы дают все собственные значения, другие дают несколько значений или даже всего одно, однако и эти алгоритмы можно использовать для вычисления всех собственных значений. Как только собственное значение λ матрицы A определено, его можно использовать либо для приведения алгоритма к получению другого собственного значения, либо для сведения задачи к такой, которая не имеет λ в качестве решения.

Приведение обычно осуществляется сдвигом: A заменяется на A - μE для некоторой константы μ. Собственное значение, найденное для A - μE, должно быть добавлено к μ, чтобы получить собственное значение матрицы A. Например, в степенном методе μ = λ. Итерация степенного метода находит самое большое по абсолютной величине значение, так что даже если λ является приближением к собственному значению, итерация степенного метода вряд ли найдёт его во второй раз. И наоборот, методы, основанные на обратных итерациях находят наименьшее собственное значение, так что μ выбирается подальше от λ в надежде оказаться ближе к какому-нибудь другому собственному значению.

Приведение можно совершить путём сужения матрицы 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 растущая сложность формул корней^[англ.] делает этот подход менее привлекательным.

Для матриц 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.

• Список алгоритмов вычисления собственных значений^[англ.]

• Дж. Голуб, Ч. Ван Лоун. Матричные вычисления. — Москва: «Мир», 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. — Jan 1991. — Т. 12, вып. 1. — С. 41—48. — doi:10.1137/0612005.