Степенной метод и определения

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

• Также отметим, что для Эрмитовых матриц \(\left(A=A^*\right)\) собственные значения всегда действительны (докажите!).

• Степенной метод - простейший метод вычисления максимального по модулю собственного значения.

Наивный способ вычисления собственных чисел.

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

\[ p(\lambda)=\operatorname{det}(A-\lambda I) \]

1. Вычислить коэффициенты многочлена

2. Найти его корни Это хорошая идея? Посмотрим на небольшой пример

import numpy as np
import matplotlib.pyplot as plt
n = 40
a = [[1.0 / (i - j + 0.5) for i in range(n)] for j in range(n)]
a = np.array(a)
ev = np.linalg.eigvals(a)
#There is a special numpy function for chacteristic polynomial
cf = np.poly(a)
ev_roots = np.roots(cf)
#print('Coefficients of the polynomial:', cf)
#print('Polynomial roots:', ev_roots)
plt.scatter(ev_roots.real, ev_roots.imag, marker='x', label='np.roots')
b = a + 0.01 * np.random.randn(n, n)
ev_b = np.linalg.eigvals(b)
plt.scatter(ev.real, ev.imag, marker='o', label='np.linalg.eigvals')
#plt.scatter(ev_b.real, ev_b.imag, marker='o', label='Малое возмущение')
plt.legend(loc='best')
plt.xlabel('Real part')
plt.ylabel('Imaginary part') 

Text(0, 0.5, 'Imaginary part')

Мораль:

• Не делайте так, если только у вас нет серьёзной причины

• Поиск корней полинома - очень плохо обусловленная задача (задача может быть обусловлена не так плохо, но с использованием другого базиса в пространстве многочленов). Заметим, что матрицей Грама для мономов

\[ h_{i j}=\int_{0}^{1} x^{i} x^{j} d x=\frac{1}{i+j+1} \]

является матрица Гильберта, у которой сингулярные числа экспоненциально убывают. Таким образом, мономы почти линейно зависимы.

Круги Гершгорина

• Есть интересная теорема, которая часто помогает локализовать собственные значения матрицы \(A = (a_{ij})\).

• Она называется теоремой Гершгорина.

• Она утверждает, что все собственные значения \(\lambda_i, i=1, \ldots, n\) находятся внутри объединения кругов Гершгорина \(C_i\), где \(C_i\) - окружность на комплексной плоскости с центром в \(a_{i i}\) и радиусом

\[ r_i=\sum_{j \neq i}\left|a_{i j}\right| \]

Более того, если круги не пересекаются, то они содержат по одному собственному значению внутри каждого круга.

Доказательство.

Сначала покажем, что если матрица \(A\) обладает строгим диагональным преобладанием, то есть

\[ \left|a_{i i}\right|>\sum_{j \neq i}\left|a_{i j}\right|, \]

тогда такая матрица невырождена. Разделим диагональную и недиагональную часть и получим

\[ A=D+S=D\left(I+D^{-1} S\right) \]

где \(\left\|D^{-1} S\right\|_1<1\). Поэтому, в силу теоремы о ряде Неймана, матрица \(I+D^{-1} S\) обратима и, следовательно, \(A\) также обратима. \(\Large \blacksquare\)

Теперь докажем утверждение теоремы от противного:

• если любое из собственных чисел лежит вне всех кругов, то матрица \((A-\lambda I)\) обладает свойством строгого диагонального преобладания

• поэтому она обратима

• это означает, что если \((A-\lambda I) x=0\), то \(x=0\).

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
n = 3
fig, ax = plt.subplots(1, 1)
a = [[5, 1, 1], [1, 0, 0.5], [2, 0, 10]]
#a = [[1.0 / (i - j + 0.5) for i in xrange(n)] for j in xrange(n)]
a = np.array(a)
#a = np.diag(np.arange(n))
a = a + 2 * np.random.randn(n, n)
#u = np.random.randn(n, n)
#a = np.linalg.inv(u).dot(a).dot(u)
xg = np.diag(a).real
yg = np.diag(a).imag
rg = np.zeros(n)
ev = np.linalg.eigvals(a)
for i in range(n):
    rg[i] = np.sum(np.abs(a[i, :])) - np.abs(a[i, i])
    crc = plt.Circle((xg[i], yg[i]), radius=rg[i], fill=False)
    ax.add_patch(crc)
plt.scatter(ev.real, ev.imag, color='r', label="Eigenvalues")
plt.axis('equal')
plt.legend()
ax.set_title('Eigenvalues and Gershgorin circles')
fig.tight_layout()

Замечание: Существуют более сложные фигуры, под названием овалы Cassini, которые содержат спектр

\[ M_{i j}=\left\{z \in \mathbb{C}:\left|a_{i i}-z\right| \cdot\left|a_{j j}-z\right| \leq r_{i} r_{j}\right\}, \quad r_{i}=\sum_{l \neq i}\left|a_{i l}\right| \]

https://en.wikipedia.org/wiki/Cassini_oval

Степенной метод

Задача на собственые значения

\[ A x=\lambda x, \quad\|x\|_2=1 \text { для устойчивости. } \]

может быть записана как итерации с неподвижной точкой, которые называются степенным методом и дают максимальное по модулю собственное значение матрицы \(A\).

Степенной метод имеет вид

\[ x_{k+1}=A x_k, \quad x_{k+1}:=\frac{x_{k+1}}{\left\|x_{k+1}\right\|_2} . \]

и \(x_{k+1} \rightarrow v_1\), где \(A v_1=\lambda_1 v_1\) и \(\lambda_1\) максимальное по модулю собственное значение, и \(v_1-\) соответствующий собственный вектор. На \((k+1)\)-ой итерации приближение для \(\lambda_1\) может быть найдено следующим образом

\[ \lambda^{(k+1)}=\left(A x_{k+1}, x_{k+1}\right), \]

Заметим, что \(\lambda^{(k+1)}\) не требуется для \((k+2)\)-ой итерации, но может быть полезно для оценки ошибки на каждой итерации: \(\left\|A x_{k+1}-\lambda^{(k+1)} x_{k+1}\right\|\).

Метод сходится со скоростью геометричекой прогрессии, с константой \(q=\left|\frac{\lambda_2}{\lambda_1}\right|<1\), где \(\lambda_1>\lambda_2 \geq \cdots \geq \lambda_n\). Это означает, что сходимость может быть сколь угодно медленной при близких значениях у \(\lambda_1\) и \(\lambda_2\).

Анализ сходимости для эрмитовой матрицы

• Рассмотрим степенной метод более подробно для случая эрмитовой матрицы

• Через несколько слайдов вы увидите, что любая эрмитова матрица диагонализуема, поэтому существует ортонормированный базис из собственных векторов \(v_1, \ldots, v_n\) такой что \(A v_i=\lambda_i v_i\).

• Разложим \(x_0\) в этом базисе с коэффициентами \(c_i\):

\[ x_0=c_1 v_1+\cdots+c_n v_n . \]

• Поскольку \(v_i\) - собственные векторы, выполнены следующие равенства

\[\begin{split} \begin{aligned} x_1 &=\frac{A x_0}{\left\|A x_0\right\|}=\frac{c_1 \lambda_1 v_1+\cdots+c_n \lambda_n v_n}{\left\|c_1 \lambda_1 v_1+\cdots+c_n \lambda_n v_n\right\|} \\ & \vdots \\ x_k &=\frac{A x_{k-1}}{\left\|A x_{k-1}\right\|}=\frac{c_1 \lambda_1^k v_1+\cdots+c_n \lambda_n^k v_n}{\left\|c_1 \lambda_1^k v_1+\cdots+c_n \lambda_n^k v_n\right\|} \end{aligned} \end{split}\]

• Получаем следующее выражение

\[ x_k=\frac{c_1}{\left|c_1\right|}\left(\frac{\lambda_1}{\left|\lambda_1\right|}\right)^k \frac{v_1+\frac{c_2}{c_1} \frac{\lambda_2^k}{\lambda_1^k} v_2+\cdots+\frac{c_n}{c_1} \frac{\lambda_n^k}{\lambda_1^k} v_n}{\left\|v_1+\frac{c_2}{c_1} \frac{\lambda_2^k}{\lambda_1^k} v_2+\cdots+\frac{c_n}{c_1} \frac{\lambda_n^k}{\lambda_1^k} v_n\right\|} \]

которое сходится к \(v_1\) при \(\left|\frac{c_1}{\left|c_1\right|}\left(\frac{\lambda_1}{\left|\lambda_1\right|}\right)^k\right|=1\) и \(\left(\frac{\lambda_2}{\lambda_1}\right)^k \rightarrow 0\) если \(\left|\lambda_2\right|<\left|\lambda_1\right|\).

Что необходимо помнить про степенной метод.

• Степенной метод даёт оценку для максимального по модулю собственного числа или спектрального радиуса матрицы

• Одна итерация требует одного умножения матрицы на вектор. Если можно умножить вектор на матрицу за \(\mathcal{O}(n)\) (напрмиер, она разреженная), тогда степенной метод можно использовать для больших \(n\)

• Сходимость может быть медленной

• Для грубой оценки максимального по модулю собственного значения и соответствующего вектора достаточно небольшого числа итераций

• Вектор решения лежит в Крыловском подпространстве \(\left\{x_0, A x_0, \ldots, A^k x_0\right\}\) и имеет вид \(\mu A^k x_0\), где \(\mu\) нормировочная постоянная.

Для какой матрицы легко найти весь спектр?

Существует класс матриц, для которого собственные числа можно найти очень легко, - это треугольные матрицы

\[\begin{split} A=\left(\begin{array}{ccc} \lambda_1 & * & * \\ 0 & \lambda_2 & * \\ 0 & 0 & \lambda_3 \end{array}\right) \end{split}\]

Собственные числа матрицы \(A-\lambda_1, \lambda_2, \lambda_3\). Почему? Потому что детерминант имеет вид

\[ \operatorname{det}(A-\lambda I)=\left(\lambda-\lambda_1\right)\left(\lambda-\lambda_2\right)\left(\lambda-\lambda_3\right) . \]

Таким образом, вычисление собственных значений для треугольной матрицы - простая задача. Теперь на помощь приходят унитарные матрицы. Пусть \(U\) унитарная матрица, то есть \(U^* U=I\). Тогда выполнены следующие равенства

\[ \operatorname{det}(A-\lambda I)=\operatorname{det}\left(U\left(U^* A U-\lambda I\right) U^*\right)=\operatorname{det}\left(U U^*\right) \operatorname{det}\left(U^* A U-\lambda I\right)=\operatorname{det}\left(U^* A U-\lambda I\right), \]

где мы используем свойства детерминанта от произведения матриц, \(\operatorname{det}(A B)=\operatorname{det}(A) \operatorname{det}(B)\). Это означает, что матрицы \(U^* A U\) и \(A\) имеют одинаковые характеристические многочлены, и, следовательно, одинаковые собственные значения. Если мы приведём матрицу \(A\) к верхнетреугольному виду \(T\) с помощью унитарной матрицы \(U: U^* A U=T\), мы решили задачу! Умножим слева и справа на \(U\) и \(U^*\) соответственно, получим нужное нам разложение:

\[ A=U T U^* . \]

• Это знаменитое разложение Шура.

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

• Разложение Шура показывает, почему нам нужны матричные разложения: они представляют матрицу в виде произведения трёх матриц подходящей структуры.

Теорема Шура.

Теорема: Каждая матрица \(A \in \mathbb{C}^{n \times n}\) может быть представлена в виде формы Шура \(A=U T U^*\), где \(U\) унитарная, а \(T\) верхнетреугольная. Набросок доказательства.

1. Каждая матрица имеет как минимум один ненулевой собственный вектор (для корня характеристического многочлена матрица \((A-\lambda I)\) вырождена и имеет нетривиальное ядро). Пусть

\[ A v_1=\lambda_1 v_1, \quad\left\|v_1\right\|_2=1 \]

2. Пусть \(U_1=\left[v_1, v_2, \ldots, v_n\right]\), где \(v_2, \ldots, v_n\) любые векторы ортогональные \(v_1\). Тогда

\[\begin{split} U_1^* A U_1=\left(\begin{array}{cc} \lambda_1 & * \\ 0 & A_2 \end{array}\right) \end{split}\]

где \(A_2\) матрица размера \((n-1) \times(n-1)\). Это называется блочнотреугольной формой. Теперь мы можем проделать аналогичную процедуру для матрицы \(A_2\) и так далее.

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

Нормальные матрицы.

Важное приложение теоремы Шура связано с так называемыми нормальными матрицами. Определение. Матрица \(A\) называется нормальной матрицей, если

\[ A A^*=A^* A . \]

Q: Какие примеры нормальных матриц вы можете привести?

Примеры: эрмитовы матрицы, унитарные матрицы.

Нормальные матрицы

Tеорема: \(A\) - нормальная матрица, тогда и только тогда, когда \(A=U \Lambda U^*\), где \(U\) унитарна и \(\Lambda\) диагональна.

Набросок доказательства:

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

В другую сторону доказательство более сложное. Рассмотрим форму Шура для матрицы \(A\). Тогда \(A A^*=A^* A\) означает, что \(T T^*=T^* T\).

Сравнив поэлементно матрицы в левой и правой части, сразу становится видно, что единственная верхнетреугольная матрица, для которой такое равенство будет выполнено - это диагональная матрица!

Важное следствие

Любая нормальная матрица - унитарно диагонализуема. Это означает, что она может быть приведена к диагональному виду с помощью унитарной матрицы U. Другими словами, каждая нормальная матрица имеет ортогональный базис из собственных векторов.

Псевдоспектр.

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

• Тогда, если \(А\) эрмитова матрица, отношение Релея определяется как

\[ R_{A}(x)=\frac{(A x, x)}{(x, x)} \]

и максимальное собственное значение равно максимальному значению \(R_{A}(x)\), аналогично для минимального собственного значения.

• Таким образом, мы можем использовать методы оптимизации для поиска этих экстремальных собственных значений.

А теперь приведём понятие, которое является обобщением собственных чисел

• Для динамических систем с матрицей \(А\), спектр может много сообщить о поведении системы (например, о её устойчивости)

• Однако для не нормальных матриц, спектр может быть неустойчивым относительно малых возмущений матрицы

• Для измерения подобных возмущений было разработана концепция псевдоспектра.

Псевдоспектр.

Рассмотрим объединение всех возможных собственных значений для всевозможных возмущений матрицы \( A .\)

\[ \Lambda_{\epsilon}(A)=\left\{\lambda \in \mathbb{C}: \exists E, x \neq 0:(A+E) x=\lambda x, \quad\|E\|_{2} \leq \epsilon .\right\} \]

Для малых \(Е\) и нормальных \(А\) это круги вокруг собственных значений. Для не нормальных матриц структура может сильно отличаться. Подробности можно найти тут: http://www.cs.ox.ac.uk/pseudospectra/