В поиске собственных значений (матриц)

Как найти собственные числа и собственные значения матрицы? Методы, излагаемые в курсе линейной алгебры, основанные на определении — применимы ли они к реальным данным? Существует ли простой алгоритм поиска этих величин, который можно понять, а не просто поверить?

Теория

Кратко обрисуем предмет нашего интереса: пусть A — квадратная матрица размера n на n, тогда ненулевой вектор x, удовлетворяющий равенству $Ax=\lambda x$ , называется собственным вектором, а число $\lambda$ — собственным значением. Интуиция, стоящая за определением, говорит нам следующее: линейный оператор, выраженный матрицей A, действует в направлении вектора X как сжатие/растяжение.

И ещё немного теоретических отступлений

Несложно заметить, что я действительно обрисовываю теорию очень краткими мазками: более того, здесь и всюду дальше я буду считать, что элементы матрицы — вещественные числа из $\mathbb{R}$ , и только при редких оказиях я буду упоминать комплексные числа $\mathbb{C}$ . Тот факт, что линейное пространство и соответствующие линейные операторы и матрицы могут быть построены над произвольным полем — здесь нам не нужен, бесконечномерные линейные пространства — тем более. Если у вас есть какие-то вопросы к этому параграфу, но не к предыдущему — просто продолжайте читать дальше.

Непосредственно из матричных операций и определения детерминанта получается характеристическое — оно же вековое — уравнение: $\det(A - \lambda I)$ , раскрываемое в монический многочлен степени n относительно лямбда. Согласно основной теореме алгебры и теореме Безу, такое уравнение всегда имеет n возможно комплексных корней (вопрос кратности корней и связанных с этим неприятностей — рассмотрим позднее).

Наивный метод поиска предлагает нам решить характеристическое уравнение, подставить полученные корни в систему $(A - \lambda I) x = 0$ , решить полученную однородную систему линейных уравнений, получив таким образом собственные вектора x.

Казалось бы — а что тут плохого? Тут всё отлично, если у вас n=2, потому что все мы со школы умеем решать уравнения второй степени, таким образом мы легко найдём собственные числа, и, наверное, систему из двух линейных уравнений — тоже решим. Более того, уже за нас всё давно решили, даже в аналитическом (общем) виде, смотреть здесь. В случае n=3 — если это пример с «подогнанными» числами, скорее всего, получится угадать собственные значения и решить системы. Если числа вполне произвольные, то, в принципе, есть формула корней кубического уравнения, в ещё более ужасном виде она существует для уравнений 4-й степени, для степеней 5 и выше — её нет, иначе говоря, единственный рабочий вариант — применение численных методов. Вот к этому мы и перейдём.

Почему «прямой» метод не прямой

Может возникнуть идея прямой работы через определение, с уже численными методами решения характеристического уравнения. На самом деле при попытке такого подхода упускается ещё один шаг: получение характеристического уравнения. Действительно, чтобы получить его, нужно раскрыть детерминант, причём в символьном виде, если опять же следовать по его определению, то сложность его вычисления будет n!.. Таким образом, сначала нужно применить менее наивный метод нахождения характеристического полинома. И только после этого приступать к его численным решениям.

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

▍ Идея метода

Предположим, для матрицы A существует собственный базис $\{e_i\}_i, i=1..n$ — базис, состоящий из собственных векторов, имеющий размерность n. Возьмём произвольный (случайный) вектор x, рассмотрим $inline$ , казалось бы, какое-то странное выражение, и непонятно, что дальше с ним делать. Но не зря был упомянут собственный базис: с одной стороны, x раскладывается в нём в следующую сумму $\sum_{n=1}^{n} x_i e_i$ , а с другой — собственные вектора $inline$ , «проходя» через умножение на A, дают выражения $\{\lambda_i e_i\}_i, i=1..n$ . Из вышесказанного получаем $A^k x = \sum_{i=1}^{i=n} \lambda_i^k x_i e_i$ . Теперь сделаем ещё одно предположение: одно из собственных чисел строго больше всех остальных (по модулю). Тогда становится достаточно ясно, что одно из слагаемых будет превосходить все другие, при росте k в выражении $inline$ . Сначала покажем, как «вытащить» собственный вектор, позже вернёмся к собственному числу.

На данный момент в наших рассуждениях мы всё ещё оперируем значениями, которых мы не знаем, всё, что нам по-настоящему известно — это матрица A и вектор $inline$ , но на самом деле этого достаточно. Последнее недостающее звено — деление $inline$ на его норму $\lVert \vec{A^kx} \rVert$ , где под последней мы будем понимать квадратный корень из суммы квадратов компонентов, иначе говоря — Евклидову норму. Действительно, пусть $\lambda_1$ — наибольшее собственное число, вытащив $x_1 \lambda_1^k$ одновременно из числителя и знаменателя дроби $\frac{A^kx}{\lVert A^kx \rVert}$ $= \frac{x_1 \lambda_1^k (e_1 + \frac{\sum_{i=2}^{i=n} \lambda_i^k x_i}{x_1 \lambda_1^k})}{x_1 \lambda_1^k \sqrt{1 + \frac{\sum_{i=2}^{i=n} (\lambda_i^k x_i)^2}{(x_1 \lambda_1^k)^2})}}$ , мы получим числитель, стремящийся к $inline$ , и знаменатель, стремящийся к 1.

Имея собственный вектор, можно получить соответствующее собственное число, используя отношение Рэлея:

$\frac{((Ax),x)}{(x,x)} = \lambda \frac{(x,x)}{(x,x)} = \lambda$

▍ Об условиях и ограничениях

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

Итак, я почти уверен, что в любом курсе линейной алгебры затрагивается тема собственного базиса, но вот то, что он не обязан существовать — вполне возможно, часто пропускается или «зажёвывается». Не погружаясь глубоко в детали: для кратных корней характеристического уравнения может не существовать достаточно (то есть столько, сколько есть кратность корня) линейно-независимых собственных векторов, в связи с чем вводится понятие алгебраических и геометрических кратностей. Пожалуй, я ограничусь примером, рассмотрев матрицу: $\begin{pmatrix} 1 & 1\\ 0 & 1 \end{pmatrix}$ , чьё собственное значение, имеющее кратность два, равно одному и собственный вектор $\begin{pmatrix} 1\\ 0 \end{pmatrix}$ — предлагаю проверить читателям самим, для развития некоторой интуиции. Но на самом деле, это не просто пример «плохой» матрицы, без собственного базиса, это пример Жордановой клетки, где на диагонали стоит собственное число (единица), и единицы (для любого значения собственного числа) расположены в верхнем треугольнике над ней.

На Хабре есть очень обстоятельная статья на эту тему, здесь же я замечу только, что существует некий аналог собственного базиса (причём часть там будет собственными векторами, а часть так называемыми присоединёнными) в котором любая матрица представляется в виде блоков — клеток Жордана, причём наше обычное разложение по собственному базису можно рассматривать как частный случай, с очень простой клеткой Жордана. И самое главное: Жорданова матрица нильпотентна, то есть при умножении саму на себя, единицы будут сдвигаться вправо, а значит не позже, чем через n итераций, для нас не будет никакой разницы между собственным базисом и этим базисом, а значит, метод может работать вне зависимости от диагнолизируемости матрицы A.

Пункт насчёт существования максимального по модулю собственного значения, на самом деле означает, что максимальное значение не может быть кратным корнем, и оно же не может парой комплексно-сопряжённых комплексных корней, потому что их модуль одинаков, в обеих случаях последовательность $\frac{A^kx}{\lVert A^kx \rVert}$ не сойдётся.
Парочка примеров на эту тему, единичная матрица $\begin{pmatrix} 1 & 0\\ 0 & 1 \end{pmatrix}$ , которая уже представлена в диагональном виде, в тоже время очевидно, что $\frac{A^kx}{\lVert A^kx \rVert}$ для любого сколько угодно большого k мы опять получим вектор x, пусть и нормированный.

Пример похитрее, оператор вращения вектора на плоскости — тоже линейный, матрицу его приводить не нужно — мы и так знаем, что делает умножение вектора на него — вращает вектор, соответственно, наша последовательность никогда не сойдётся к собственному (и любому другому) вектору. Более того, это в принципе невозможно: из этого же определения ясно, что собственные вектора такого оператора — комплексные, кроме случая углов поворота, кратным 180-ти градусам — наша последовательность не может выдать комплекснозначный результат.

Обещанная ремарка: есть исчезающе малый шанс выбрать начальный вектор x не удачно, а именно так, что в его разложении в собственном/жордановом базисе его проекция на искомый $inline$ равняется нулю, иначе говоря, x ортогонален $inline$ .

▍ Имплементация

Посмотрим результаты работы метода, вооружившись для начала матрицей размера n=2 — в этом случае легко наглядно показать результаты работы метода и его сходимость к собственному вектору. Будем использовать язык python с библиотеками numpy и matplotlib, как удобные средства оперирования матрицами и создания визуализаций.

import numpy as np from matplotlib.animation import FuncAnimation import matplotlib.pyplot as plt from numpy import linalg as LA  # given matrix A = np.array([[4, 2], [3, 4]])  eigenvalues, eigenvectors = LA.eig(A) # an eigenvector corresponding to the max eigenvalue eigenvector = -eigenvectors[:,0]  print("eigenvalues", eigenvalues) print("eigenvector",eigenvector)  k = 10 X_origin = [0, 0] Y_origin = [0, 0]  # Starting vector #X = np.random.rand(A.shape[0]) # initial guess vector X = np.array([-0.7, 0.7])  fig, ax = plt.subplots() plt.title("Convergence of the power iteration method") plt.grid() plt.xlim(-1, 1) plt.ylim(-1, 1) U = [eigenvector[0], X[0]] V = [eigenvector[1], X[1]]  q = ax.quiver(X_origin, Y_origin, U, V, color=['b', 'r'], angles='xy', scale_units='xy', scale=1)  def animate(i):     # Power it     global X     if i:          X = np.matmul(A, X)         # Normalize it         X = X / np.linalg.norm(X, ord=2)      U = [eigenvector[0], X[0]]     V = [eigenvector[1], X[1]]     q.set_UVC(U, V)  ani = FuncAnimation(fig, animate, frames=k, repeat=False, interval=1000)  plt.show()  l = np.dot(np.matmul(A, X), X) / np.dot(X, X)  print("lamda = ", l) print("X = ", X)

Синим цветом обозначен искомый собственный вектор, красным — начальное приближение и результаты итераций. Видно, что, несмотря на начальный угол, близкий к 90 градусам, вектор приближений очень быстро поворачивается до искомого вектора, начиная визуально совпадать с ним с 5-й итерации.

Посмотрим вывод:

eigenvalues [6.44948974 1.55051026]
eigenvector [-0.63245553 -0.77459667]
lamda = 6.449515757348442
X = [-0.63247569 -0.77458021]

Точность в собственном векторе — на уровне 5-го знака после запятой, в собственном значении — на уровне 4-го знака.

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

import numpy as np import matplotlib.pyplot as plt from numpy import linalg as LA  np.random.seed(42) n = 50 A = np.random.randint(-100, 100, size=(n, n)) # make it symmetric  A = (A + A.T) / 2 eigenvalues, eigenvectors = LA.eig(A) eigenvalues_abs = np.abs(eigenvalues) index = np.argmax(eigenvalues_abs)  lamda1 = eigenvalues[index]  eigenvector = eigenvectors[:,index] X = np.random.rand(A.shape[0]) k = 40   diff = np.zeros(k) diff_x = np.linspace(1, k, num=k)  for i in range(k):     X = np.matmul(A, X)     # Normalize it     X = X / np.linalg.norm(X, ord=2)          diff[i] = np.linalg.norm(X - eigenvector)  eigenvalues_abs[index] = 0 index = np.argmax(eigenvalues_abs)  lamda2 = eigenvalues[index]  print("rate = ", lamda2 / lamda1) print("diff", diff)  fig, ax = plt.subplots() ax.plot(diff_x, diff) plt.title("Convergence of the power iteration method") plt.grid() plt.show()

В результате чего мы получим следующий график:

Тут можно понять, что норма разницы уменьшается, но как именно? Лучше перерисовать ошибку в логарифмическом масштабе:

Тут можно заметить, что после определённого шага ошибка становится «линейна», то есть ведёт себя как геометрическая прогрессия со знаменателем, равным отношению второго по модулю собственного числа к максимальному собственному числу (rate в скрипте).

Вывод работы скрипта, на котором можно проверить, что это действительно так:

rate = 0.9139630521574826
diff [1.33026844 1.28732128 1.24311358 1.19893046 1.15496257 1.1105661
1.06502103 1.01786565 0.96895095 0.91839837 0.86653833 0.81384951
0.7609016 0.70830146 0.65664456 0.60647469 0.55825471 0.51234941
0.46902015 0.42842897 0.39064931 0.3556805 0.32346342 0.29389557
0.26684447 0.24215888 0.21967773 0.19923697 0.18067448 0.16383357
0.14856524 0.13472953 0.12219622 0.110845 0.10056539 0.09125633
0.08282571 0.07518982 0.0682727 0.06200557]

Что могло бы быть дальше

Если рассматривать метод формально, то у него существует очевидное ограничение: нахождение только одной пары вектор/число за раз. Это легко преодолимо: если вместо исходной матрицы A, рассмотреть матрицу $A-\lambda E$ , где лямбда — найденное собственное значение, то лямбда «зануляется», позволяя искать следующее по величине значение.

Давайте вернёмся к матрице n=2 (для простоты) и проверим это:

A = np.array([[4, 2], [3, 4]])  eigenvalues, eigenvectors = LA.eig(A)  print("eigenvectors", eigenvectors) print("eigenvalues", eigenvalues)  k = 10  A = A - np.identity(2) * eigenvalues[0] X = np.random.rand(A.shape[0])  for i in range(k):     X = np.matmul(A, X)     # Normalize it     X = X / np.linalg.norm(X, ord=2)      print("calculated eigenvector = ", X) print("calculated eigenvalue=", np.dot(np.matmul(A, X), X) / np.dot(X, X) + eigenvalues[0])

Получаемый вывод:

eigenvectors [[ 0.63245553 -0.63245553]
[ 0.77459667 0.77459667]]
eigenvalues [6.44948974 1.55051026]
calculated eigenvector = [ 0.63245553 -0.77459667]
calculated eigenvalue= 1.5505102572168221

Всё работает, очевидно, что можно обобщить метод на случай n > 2, но… Это предполагает довольно много повторяющихся действий. В примере с n = 50 я уже употреблял приём «симметризации» матрицы, можно пойти дальше и сказать, что на самом деле, часто нас интересуют собственные значения и вектора именно симметричных матриц. Что, если нам изначально вооружиться n случайными, причём ортогональными (потому что искомые вектора — ортогональны) векторами, и попробовать реализовать тот же самый алгоритм? — На самом деле это не сработает, но это послужит своеобразным мостиком к QR алгоритму поиска собственных значений.

Также похожие идеи используются в методах, основанных на подпространстве Крылова.

Выводы

Степенной метод обладает рядом несомненных достоинств:

Простота теории, лежащей в его основе.
Простота реализации.

К недостаткам можно отнести:

Находит одну пару вектор/значение за раз.
Скорость сходимости может быть плохой при близких величинах собственных чисел.
Для произвольной матрицы без дополнительных проверок может не сойтись.

Источники

Telegram-канал со скидками, розыгрышами призов и новостями IT ?

ссылка на оригинал статьи https://habr.com/ru/articles/845652/