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

У Google миллиарды веб-страниц, и каждой нужен рейтинг важности. Рейтинг — это компоненты собственного вектора матрицы ссылок размером в миллиарды строк: PageRank — гигантский степенной метод на матрице переходов цепи Маркова, крутящийся на обычном железе. Умножай вектор на матрицу, нормируй, повторяй — и спектр сам проявится. Здесь разбираем эту идею в миниатюре: от итераций в столбик до отношения Рэлея, сдвигов и моста к QR-алгоритму.

Почему характеристический многочлен — тупик#

В курсе линейной алгебры собственные значения ищут через характеристическое уравнение: для матрицы 2×2 это квадратное уравнение, для 3×3 — кубическое. Логичный ход — продолжить и для больших матриц. Упираемся в стену: начиная с пятой степени формул корней не существует в принципе (теорема Абеля—Руффини), а характеристическое уравнение степени 1000 — это 1001 коэффициент, миллион слагаемых в старших и никакого шанса на честное решение.

Даже если многочлен как-то выписать, корни катастрофически чувствительны к коэффициентам. Классический пример Уилкинсона: у многочлена $(x-1)(x-2)\cdots (x-20)$ все корни — целые числа от 1 до 20. Сдвиньте коэффициент при $x^{19}$ на $2^{-23}\approx 0,00000012$ — и пара вещественных корней 10 и 11 превращается в комплексные $10,10\pm 0,64i$. Коэффициенты характеристического многочлена считаются через определители с округлениями — спектр оказывается испорчен ещё до решения уравнения. Отсюда железное правило численной линейной алгебры: не превращай матрицу в многочлен, работай с ней напрямую.

Степенной метод: одно умножение за итерацию#

Степенной метод игнорирует многочлен целиком. Нужна одна операция — умножение матрицы на вектор. Стартуем с любого вектора $x_0$ и повторяем:

Почему это вообще работает? Разложим стартовый вектор по собственным векторам: $x_0=c_1{v}_1+c_2{v}_2+\dots$ Каждое умножение тянет из слагаемого свой множитель, и после $k$ шагов $A^k{x}_0=c_1{\lambda }_1^k{v}_1+c_2{\lambda }_2^k{v}_2+\dots$ Если $|{\lambda }_1|$ строго больше остальных модулей, слагаемое с ${\lambda }_1^k$ обгоняет всех: нормированный вектор поворачивается к $v_1$, а норма $\Vert y_k\Vert$ тянется к $|{\lambda }_1|$. Из каждой итерации выжимают три вещи:

• оценка по норме: $\hat{\lambda }=\Vert y_k\Vert$ — длина вектора до нормирования
• оценка по компонентам: $\hat{\lambda }=(y_k{)}_i/(x_k{)}_i$ — частное соответствующих компонент
• критерий остановки: $|{\hat{\lambda }}_{k+1}-{\hat{\lambda }}_k|<\epsilon$ — приращение оценки застопорилось; тот же принцип, что в методе Ньютона

Практическая деталь: в стартовом векторе должна быть хоть капля $v_1$, то есть $c_1\ne 0$. Если $x_0$ точно ортогонален главному направлению, метод уползёт к следующему собственному вектору. В машинной арифметике округления обычно сами подмешивают нужную каплю, но рассчитывать на это как на гарантию не стоит.

Прогоняем руками: матрица с корнями 5 и 2#

Достаточно теории — посчитаем. Берём матрицу, которую уже щупали в интерактиве:

Собственные значения ${\lambda }_1=5$ и ${\lambda }_2=2$, векторы $v_1=(1;1)$ и $v_2=(1;-2)$. Проверьте умножением: $A(1;1)=(5;5)$ и $A(1;-2)=(2;-4)=2(1;-2)$. Отношение модулей $5/2=2,5$ обещает сходимость к пятёрке со скоростью «ошибка в 2,5 раза меньше за шаг». Стартуем с $x_0=(1;0)$.

1. 1
   Умножаем: $A{x}_0=(4\cdot 1+1\cdot 0;2\cdot 1+3\cdot 0)=(4;2)$
2. 2
   Норма до нормирования: $\Vert y_0\Vert =\sqrt{4^2+2^2}=\sqrt{20}\approx 4,472$ — первая оценка ${\lambda }_1$
3. 3
   Нормируем: $x_1=(0,894;0,447)$ — направление уже развернулось от оси к биссектрисе
4. 4
   Вторая итерация: $A{x}_1=(4,025;3,130)$, норма $\sqrt{26}\approx 5,099$, нормируем: $x_2=(0,789;0,614)$
5. 5
   Угол между $x_2$ и $v_1=(1;1)$ — около $7^{\circ }$, хотя на старте было $45^{\circ }$

Предел виден невооружённым глазом: $x_k\to (1;1)/\sqrt{2}\approx (0,707;0,707)$ — нормированный собственный вектор, а оценки $\lambda$ тают к $5$. Заметьте две честные подробности. Норма немонотонна: $5,099\to 5,091\to 5,044$ — матрица несимметрична, путь к пределу виляет. И Рэлей на первом шаге случайно попал точно в $5,000$ — совпадение, а не закономерность: после второго шага ошибка вернулась к типичному затуханию.

Ошибка тает в 2,5 раза за итерацию#

Скорость степенного метода зашита в спектре: ошибка примерно пропорциональна $({\lambda }_2/{\lambda }_1{)}^k$ — отношению модулей первого соперника и лидера. Здесь $(2/5{)}^k=0,4^k$: 0,4; 0,16; 0,064; 0,026… Углы в таблице послушно делятся на 2,5: $18,4^{\circ }\to 7,1^{\circ }\to 2,8^{\circ }\to 1,1^{\circ }$. Перевод в знаки: $0,4^k=0,1$ при $k\approx 2,5$, то есть каждые два-три шага метод выигрывает один верный десятичный знак.

Медленно? По меркам Гаусса — да, но сравнение нечестное. Итерация стоит $n^2$ умножений для плотной матрицы, а для разреженной — числа ненулевых элементов, которых может быть $5n$ вместо $n^2$. Память — один вектор, $O(n)$. Именно поэтому метод жив на гигантских разреженных задачах, где плотные алгоритмы не помещаются ни во время, ни в диски: PageRank с демпфированием 0,85 жмёт ошибку в $0,85^k$ раз — за сотню итераций в $10^7$ раз, и этого хватает на весь веб.

Отношение Рэлея и сдвиг: точность и наведение#

Норма — грубая оценка собственного числа. Существует точный скалярный функционал:

Если $x$ отклонён от собственного вектора на угол $\theta$, то ошибка $r(x)-{\lambda }_1\sim {\theta }^2$ — второй порядок против первого у нормы. В таблице это видно: на шаге 2 угол $7,1^{\circ }$, и Рэлей даёт $5,077$ — ошибка $0,077$ против $0,099$ по норме; на шестом шаге Рэлей уже $5,003$. Отношение Рэлея — это ещё и способ остановиться красиво: когда нужен $\lambda$ с запасом точности, итерации завершают по нему, а не по норме.

Вторая работа скаляра $\sigma$ — наведение. Матрица $A-\sigma E$ имеет собственные значения ${\lambda }_i-\sigma$: векторы те же, спектр сдвинут целиком. Подбирая $\sigma$, можно сделать наибольшим по модулю нужное число. Пример: сдвиг $\sigma =1,5$ даёт для матрицы-героя числа $3,5$ и $0,5$ — теперь доминирует $3,5$, то есть то самое $\lambda =2$, к которому исходный метод не собирался. Бонус: отношение модулей улучшается с $0,4$ до $0,5/3,5=1/7\approx 0,14$ — сходимость почти втрое резвее. Нашли $3,5$ — прибавьте сдвиг: $\lambda =3,5+1,5=2$.

Обратный степенной метод: где взять вектор#

Практика ставит задачу наоборот: приближение $\lambda$ найдено — сдвигом, QR, из постановки задачи, — а нужен собственный вектор. Ответ изящен: итерировать с обратной матрицей. Собственные значения $A^{-1}$ суть $1/{\lambda }_i$, собственные векторы — те же. Для матрицы-героя:

У $A^{-1}$ два собственных числа: $0,5$ (вектор $(1;-2)$) и $0,2$ (вектор $(1;1)$). Обратный степенной метод сходится к вектору наименьшего по модулю $\lambda$ матрицы $A$. Связка «сдвиг плюс обращение» $(A-\sigma E{)}^{-1}$ — главный инструмент всего семейства: чем ближе $\sigma$ к искомому $\lambda$, тем сильнее доминирует нужное число и тем быстрее сходимость; при почти точном сдвиге вектор проявляется за одну-две итерации. Обратную матрицу никто не строит явно — на каждом шаге решают систему $(A-\sigma E)y=x$ итерационно, теми же методами, что и обычные СЛАУ.

От степенного метода — к QR-алгоритму#

Одну мысль стоит унести с собой: возведение в степень само ищет собственные направления. Считаем степени матрицы-героя: $A^2=\left(\begin{array}{cc}18& 7\\ 14& 11\end{array}\right)$, $A^4=\left(\begin{array}{cc}422& 203\\ 406& 219\end{array}\right)$, $A^8=\left(\begin{array}{cc}260502& 130123\\ 260246& 130379\end{array}\right)$. Отношение второй компоненты к первой внутри столбцов: $0,78$ и $1,57$, затем $0,96$ и $1,08$, затем $0,999$ и $1,002$. Столбцы $A^k$ сползают друг к другу — оба поворачиваются к направлению $(1;1)$. Матрица, возведённая в большую степень, забывает всё, кроме доминирующего направления.

QR-алгоритм — промышленное продолжение этой идеи. На каждом шаге матрицу раскладывают в произведение ортогональной $Q$ и треугольной $R$ и переставляют множители: $A_{m+1}=RQ$. Ортогональность заменяет жёсткое нормирование степенного метода, и процесс сходится сразу ко всем собственным значениям — фактически диагонализация без единого характеристического многочлена. Современные реализации (приведение к Хессенбергу плюс сдвиги) находят весь спектр плотной матрицы за $O(n^3)$, а векторы дочитывают уже знакомыми обратными итерациями со сдвигами.

Когда степенной метод не сходится#

Общий знаменатель всех болезней один: метод смотрит только на модули. Сравнялись два модуля — выделенного направления больше нет, и тянуть не к чему. Сдвиг лечит почти всё, потому что меняет расклад модулей, не трогая векторы; именно поэтому в промышленных алгоритмах сдвиги стоят на каждом шаге.

• Степенной метод: умножай на $A$, нормируй, повторяй — сходится к наибольшему по модулю $\lambda$, если спектр честно разделён
• Оценки числа: норма $\Vert y_k\Vert$, частное компонент, отношение Рэлея — точнее всех, ошибка $\sim {\theta }^2$
• Скорость — $({\lambda }_2/{\lambda }_1{)}^k$; сдвиг $\sigma$ перестраивает спектр и наводит метод на нужный корень
• Обратный степенной метод достаёт векторы; пара «QR + обратные итерации» — сердце numpy.linalg.eig и LAPACK
• PageRank — тот же степенной метод на матрице миллиардов строк; ошибка сжимается в $0,85^k$ раз за итерацию

Прогоните степенной метод на своей матрице из задачника и сверьте с ответом: посчитать три итерации в столбик — пять минут, а понимание приходит именно из ручного счёта. Потренироваться можно в тренажёре численных методов. Теория собственных пар разобрана в уроке линейной алгебры, а словарные статьи про собственный вектор и диагонализацию матрицы пригодятся, если разложение $x_0=c_1{v}_1+c_2{v}_2$ вызывает вопросы.