Краткий обзор предыдущей лекции

  • QR разложение и алгоритм Грама-Шмидта
  • Разложение Шура и QR-алгоритм (основы)

План на сегодня

Сегодня мы поговорим о:

  • Алгоритмах для симметричных задач на собственные значения
    • QR алгоритм (более подробно)
    • Разделяй и властвуй
    • Метод бисекции
  • Алгоритмах вычисления SVD

Вычисление формы Шура

  • Напомним, что мы пытаемся избежать сложности \mathcal{O}(n^3) для каждой итерации.

  • Идея состоит в том, чтобы придать матрице более простую структуру, чтобы каждый шаг QR алгоритма стал дешевле.

  • В случае произвольной матрицы мы можем использовать форму Хессенберга.

Форма Хессенберга

Матрица A называется матрицей в форме Хессенберга, если

a_{ij} = 0, \quad \mbox{если } i \geq j+2.

H = \begin{bmatrix} * & * & * & * & * \\ * & * & * & * & * \\ 0 & * & * & * & *\\ 0 & 0 & * & * & *\\ 0 & 0 & 0 & * & * \\ \end{bmatrix}.

Приведение произвольной матрицы к форме Хессенберга

  • Применяя отражения Хаусхолдера, мы можем привести любую матрицу к форме Хессенберга

U^* A U = H

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

  • Вычислительная сложность такого приведения составляет \mathcal{O}(n^3) операций.

  • В форме Хессенберга вычисление одной итерации QR алгоритма требует \mathcal{O}(n^2) операций (например, используя вращения Гивенса, как?), и форма Хессенберга сохраняется при QR итерации (проверьте почему).

Симметричный (эрмитов) случай

  • В симметричном случае у нас A = A^*, тогда H = H^* и верхняя форма Хессенберга становится трехдиагональной матрицей.

  • Далее мы будем говорить о случае симметричной трехдиагональной формы.

  • Любая симметричная (эрмитова) матрица может быть приведена к трехдиагональной форме с помощью отражений Хаусхолдера.

  • Ключевой момент заключается в том, что трехдиагональная форма сохраняется при QR алгоритме, и стоимость одного шага может быть снижена до \mathcal{O}(n)!

QR алгоритм: итерации

  • Итерации QR алгоритма имеют следующий вид:

A_k = Q_k R_k, \quad A_{k+1} = R_k Q_k.

  • Если A_0 = A является трехдиагональной симметричной матрицей , эта форма сохраняется при QR алгоритме.

Давайте посмотрим..

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

#Generate a random tridiagonal matrix

n = 20
rng = np.random.RandomState(0)
d = rng.normal(size=(n,))

rng = np.random.RandomState(1)
sub_diag = rng.normal(size=(n-1,))

mat = np.diag(d) + np.diag(sub_diag, -1) + np.diag(sub_diag, 1)
mat1 = np.abs(mat)
mat1 = mat1/np.max(mat1.flatten())
plt.spy(mat)
q, r = np.linalg.qr(mat)
plt.figure()
b = r.dot(q)
b[abs(b) <= 1e-12] = 0
plt.spy(b)
#plt.figure()
#plt.imshow(np.abs(r.dot(q)))
b[0, :]
Array([-1.47168069, -0.75157939,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ],      dtype=float64)

Трехдиагональная форма

  • В трехдиагональной форме вам не нужно вычислять матрицу Q: вам нужно вычислить только трехдиагональную часть, которая появляется после итераций

A_k = Q_k R_k, \quad A_{k+1} = R_k Q_k,

в случае, когда A_k = A^*_k и также является трехдиагональной.

  • Такая матрица определяется \mathcal{O}(n) параметрами; вычисление QR более сложное, но можно вычислить A_{k+1} напрямую без вычисления Q_k.

  • Это называется неявным QR-шагом.

Теорема о неявной QR-итерации

Все неявные QR-алгоритмы основаны на следующей теореме:

Пусть

Q^* A Q = H

является неприводимой верхней матрицей Хессенберга. Тогда первый столбец матрицы Q определяет все остальные её столбцы. Его можно найти из уравнения

A Q = Q H.

Сходимость QR-алгоритма

  • Сходимость QR-алгоритма является очень деликатным вопросом (см. Е. Е. Тыртышников, “Краткое введение в численный анализ” для подробностей).

Резюме. Если у нас есть разложение вида

A = X \Lambda X^{-1}, \quad A = \begin{bmatrix}A_{11} & A_{12} \\ A_{21} & A_{22}\end{bmatrix}

и

\Lambda = \begin{bmatrix} \Lambda_1 & 0 \\ 0 & \Lambda_2 \end{bmatrix}, \quad \lambda(\Lambda_1)=\{\lambda_1,\dots,\lambda_m\}, \ \lambda(\Lambda_2)=\{\lambda_{m+1},\dots,\lambda_r\},

и существует разрыв между собственными значениями \Lambda_1 и \Lambda_2 (|\lambda_1|\geq \dots \geq |\lambda_m| > |\lambda_{m+1}| \geq\dots \geq |\lambda_r| >0), тогда блок A^{(k)}_{21} матрицы A_k в QR-итерации стремится к нулю со скоростью

\Vert A^{(k)}_{21} \Vert \leq C q^k, \quad q = \left| \frac{\lambda_{m+1}}{\lambda_{m}} \right |,

где m - размер \Lambda_1.

Таким образом, нам нужно увеличить разрыв! Это можно сделать с помощью QR-алгоритма со сдвигами.

QR-алгоритм со сдвигами

A_{k} - s_k I = Q_k R_k, \quad A_{k+1} = R_k Q_k + s_k I

Скорость сходимости для версии со сдвигом тогда равна

\left| \frac{\lambda_{m+1} - s_k}{\lambda_{m} - s_k} \right |,

где \lambda_m - m-е по величине собственное значение матрицы по модулю. - Если сдвиг близок к собственному значению, то скорость сходимости лучше. - Существуют различные стратегии выбора сдвигов. - Введение сдвигов - это общая стратегия для улучшения сходимости итерационных методов нахождения собственных значений. - В следующих слайдах мы проиллюстрируем, как выбрать сдвиг на более простом алгоритме.

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

  • Вспомним степенной метод для вычисления собственных значений.

x_{k+1} := A x_k, \quad x_{k+1} := \frac{x_{k+1}}{\Vert x_{k+1} \Vert}.

  • Он сходится к собственному вектору, соответствующему наибольшему по модулю собственному значению.

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

  • Попробуем использовать стратегию сдвигов. Если мы сдвигаем матрицу как

A := A - \lambda_k I

то соответствующее собственное значение становится маленьким (а нам нужно большое). - Это не то, что мы хотели!

Метод обратной итерации и метод итерации Рэлея

  • Чтобы сделать малое собственное значение большим, нам нужно обратить матрицу, и это дает нам метод обратной итерации

x_{k+1} = (A - \lambda I)^{-1} x_k,

где \lambda - сдвиг, который является приближением к собственному значению, которое мы хотим найти.

  • Как и для степенного метода, сходимость является линейной.

  • Для ускорения сходимости можно использовать метод итерации Рэлея (обратная итерация с адаптивными сдвигами), который задается выбором адаптивного сдвига:

x_{k+1} = (A - \lambda_k I)^{-1} x_k,

\lambda_k = \frac{(Ax_k, x_k)}{(x_k, x_k)}

  • В симметричном случае A = A^* сходимость локально кубическая, а в остальных случаях - локально квадратичная.

Сингулярные значения и собственные значения (1)

  • Теперь поговорим о сингулярных значениях и собственных значениях.

  • SVD

A = U \Sigma V^*

существует для любой матрицы.

  • Это также можно рассматривать как приведение данной матрицы к диагональной форме с помощью двусторонних унитарных преобразований:

\Sigma = U^* A V.

  • С помощью двусторонних преобразований Хаусхолдера мы можем привести любую матрицу к бидиагональной форме B (как?).

Сингулярные значения и собственные значения (2)

  • Неявный QR алгоритм (со сдвигами) дает способ вычисления собственных значений (и формы Шура).

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

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

  1. Работать с трехдиагональной матрицей

T = B^* B

  1. Работать с расширенной матрицей

T = \begin{bmatrix} 0 & B \\ B^* & 0 \end{bmatrix}

  • Случай 1 хорош, если вы не формируете T напрямую!

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

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

Рассмотрены: - QR алгоритм: “золотой стандарт” вычисления собственных значений

Следующие слайды: - Метод “разделяй и властвуй” - Метод бисекции - Метод Якоби

Разделяй и властвуй

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

T = \begin{bmatrix} T'_1 & B \\ B^{\top} & T'_2 \end{bmatrix}

  • Мы можем записать матрицу T как

T = \begin{bmatrix} T_1 & 0 \\ 0 & T_2 \end{bmatrix} + b_m v v^*

где vv^* - матрица ранга 1, v = (0,\dots,0,1,1,0,\dots,0)^T.

  • Предположим, мы уже разложили T_1 и T_2:

T_1 = Q_1 \Lambda_1 Q^*_1, \quad T_2 = Q_2 \Lambda_2 Q^*_2

  • Тогда (проверьте),

\begin{bmatrix} Q^*_1 & 0 \\ 0 & Q^*_2 \end{bmatrix} T\begin{bmatrix} Q_1 & 0 \\ 0 & Q_2 \end{bmatrix} = D + \rho u u^{*}, \quad D = \begin{bmatrix} \Lambda_1 & 0 \\ 0 & \Lambda_2\end{bmatrix}

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

Диагональная плюс низкоранговая матрица

Вычисление собственных значений матрицы

D + \rho u u^*

является непростой задачей.

Характеристический многочлен имеет вид

\det(D + \rho uu^* - \lambda I) = \det(D - \lambda I)\det(I + \rho (D - \lambda I)^{-1} uu^*) = 0.

Тогда (докажите!!)

\det(I + \rho (D - \lambda I)^{-1} uu^*) = 1 + \rho \sum_{i=1}^n \frac{|u_i|^2}{d_i - \lambda} = 0

Подсказка: найдите \det(I + w u^*), используя тот факт, что \text{det}(C) = \prod_{i=1}^n\lambda_i(C) и \text{trace}(C) = \sum_{i=1}^n \lambda_i.

Характеристическое уравнение

1 + \rho \sum_{i=1}^n \frac{|u_i|^2}{d_i - \lambda} = 0

Как найти корни?

import numpy as np
import matplotlib.pyplot as plt

lm = np.array([1, 2, 3, 4])
M = len(lm)
D = np.array(lm)
a = np.min(lm)
b = np.max(lm)
t = np.linspace(-1, 6, 1000)
u = 0.5 * np.ones(M)
rho = 1
def fun(lam):
    return 1 + rho * np.sum(u**2/(D - lam))
res = [fun(lam) for lam in t]
plt.figure(figsize=(10,8))
plt.plot(t, res, 'k')
#plt.plot(np.zeros_like(t))
plt.ylim([-6, 6])
plt.tight_layout()
plt.yticks(fontsize=24)
plt.xticks(fontsize=24)
plt.grid(True)

  • Функция имеет только один корень на отрезке [d_i, d_{i+1}]

  • Мы, кстати, доказали теорему о перемежаемости Коши (что происходит с собственными значениями при возмущении ранга 1)

Как найти корень

  • Метод Ньютона не сработает (нарисуйте картинку с касательной линией).

  • Заметим, что метод Ньютона - это просто аппроксимация функции f(\lambda) линейной функцией.

  • Гораздо лучшей аппроксимацией является гипербола:

f(\lambda) \approx c_0 + \frac{c_1}{d_i - \lambda} + \frac{c_2}{d_{i+1} - \lambda}.

  • Чтобы подобрать коэффициенты, нам нужно вычислить f(\lambda) и f'(\lambda) в конкретной точке.

  • После этого аппроксимацию можно получить, решив квадратное уравнение

Важные проблемы

  • Во-первых, стабильность: этот метод был заброшен на долгое время из-за нестабильности вычисления собственных векторов.

  • В рекурсии нам нужно вычислить собственные векторы матрицы D + \rho uu^*.

  • Точное выражение для собственных векторов просто (давайте проверим!)

(D - \alpha_i I)^{-1}u, где \alpha_i - вычисленный корень.

  • Причина нестабильности:
    • если \alpha_i и \alpha_{i+1} близки, то соответствующие собственные векторы коллинеарны, но они должны быть ортогональными
    • если \alpha_i и \alpha_{i+1} близки, то они близки к d_i, поэтому матрицы D - \alpha_i I и D - \alpha_{i+1} I близки к сингулярным

Теорема Лёвнера

  • Решение пришло в виде использования странной теоремы Лёвнера:

Если \alpha_i и d_i удовлетворяют теореме о перемежаемости

d_n < \alpha_n < \ldots < d_{i+1} < \alpha_{i+1} \ldots

То существует вектор \widehat{u} такой, что \alpha_i являются точными собственными значениями матрицы

\widehat{D} = D + \widehat{u} \widehat{u}^*.

  • Таким образом, сначала вычисляются собственные значения, затем вычисляется \widehat{u} и только потом собственные векторы.

Разделяй и властвуй и Быстрый Мультипольный Метод

В вычислениях метода “разделяй и властвуй” нам необходимо вычислять суммы вида

f(\lambda) = 1 + \rho \sum_{i=1}^n \frac{|u_i|^2}{(d_i - \lambda)},

и делать это как минимум для n точек.

  • Сложность тогда составляет \mathcal{O}(n^2), как и для QR алгоритма.

  • Можем ли мы сделать её \mathcal{O}(n \log n)?

  • Ответ - да, но нам придётся заменить точные вычисления приближёнными с помощью Быстрого Мультипольного Метода.

Давайте объясним немного подробнее…

Еще несколько алгоритмов

  • Совершенно другой подход основан на бисекции.

  • Для матрицы A ее инерция определяется как тройка

(\nu, \zeta, \pi),

где \nu - количество отрицательных, \zeta - нулевых и \pi - положительных собственных значений.

  • Если X невырожденная, то

Inertia(A) = Inertia(X^* A X)

Бисекция с помощью метода Гаусса

  • Для заданного z мы можем выполнить исключение Гаусса:

A - zI = L D L^*,

и инерцию диагональной матрицы тривиально вычислить.

  • Таким образом, если мы хотим найти все собственные значения в интервале a, b

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

  • Иллюстрация: если Inertia(A)=(5,0,2) и после сдвига Inertia(A-zI)=(4,0,3), z\in[a,b], то мы знаем, что \lambda(A)\in[a,z].

Метод Якоби

  • Вспомним, что такое вращения Якоби (Гивенса)

  • В плоскости они соответствуют ортогональной матрице 2 \times 2 вида

\begin{pmatrix} \cos \phi & \sin \phi \\ -\sin \phi & \cos \phi \end{pmatrix},

и в n-мерном случае мы выбираем две переменные i и j и выполняем вращение.

Метод Якоби (продолжение)

  • Идея метода Якоби заключается в минимизации суммы квадратов внедиагональных элементов:

\Gamma(A) = \mathrm{off}( U^* A U), \quad \mathrm{off}^2(X) = \sum_{i \ne j} \left|X_{ij}\right|^2 = \|X \|^2_F - \sum\limits_{i=1}^n x^2_{ii}.

путем применения последовательных вращений Якоби U для обнуления внедиагональных элементов.

  • Когда “опорный элемент” выбран, его легко исключить.

  • Основной вопрос заключается в том, каков порядок проходов, которые мы должны сделать (т.е. в каком порядке исключать элементы).

  • Если мы всегда исключаем наибольшие внедиагональные элементы, метод имеет квадратичную сходимость.

  • На практике используется циклический порядок (т.е., (1, 2), (1, 3), \ldots, (2, 3), \ldots).

Метод Якоби: сходимость

  • Чтобы показать сходимость метода Якоби, нам нужно доказать, что после каждого вращения сумма квадратов внедиагональных элементов уменьшается:

\text{off}(B) < \text{off}(A),

где B = U^*AU - матрица после применения вращения.

  • Рассмотрим, что происходит при вращении. Пусть p и q - индексы строк и столбцов, которые мы выбрали для вращения.

  • Вспомним, что \text{off}^2(A) = \sum_{i \neq j} |a_{ij}|^2 - сумма квадратов всех внедиагональных элементов.

  • После вращения:

    1. Норма Фробениуса матрицы не меняется: \|B\|_F = \|A\|_F (свойство унитарных преобразований)
    2. Диагональные элементы, кроме b_{pp} и b_{qq}, остаются теми же: b_{ii} = a_{ii} для i \neq p, q
    3. Элемент a_{pq} (и симметричный ему a_{qp}) становится нулем

  • Запишем сумму квадратов внедиагональных элементов:

\text{off}^2(B) = \|B\|^2_F - \sum_{i=1}^n b_{ii}^2

= \|A\|^2_F - \sum_{i \neq p,q} a_{ii}^2 - (b_{pp}^2 + b_{qq}^2)

  • Можно показать, что b_{pp}^2 + b_{qq}^2 = a_{pp}^2 + a_{qq}^2 + 2a_{pq}^2
  • Подставляя это выражение:

\text{off}^2(B) = \|A\|^2_F - \sum_{i \neq p,q} a_{ii}^2 - (a_{pp}^2 + a_{qq}^2 + 2a_{pq}^2)

= \|A\|^2_F - \sum_{i=1}^n a_{ii}^2 - 2a_{pq}^2

= \text{off}^2(A) - 2a_{pq}^2

  • Поскольку a_{pq}^2 > 0 (мы выбираем ненулевой элемент), получаем \text{off}^2(B) < \text{off}^2(A)

  • Если мы всегда выбираем наибольший внедиагональный элемент a_{pq} = \gamma, то:

    • Все внедиагональные элементы не превосходят \gamma: |a_{ij}| \leq \gamma
    • Всего внедиагональных элементов 2N = n(n-1)
    • Поэтому \text{off}^2(A) \leq 2N\gamma^2

  • Отсюда получаем: 2\gamma^2 \geq \frac{\text{off}^2(A)}{N}

  • После вращения: \text{off}^2(B) = \text{off}^2(A) - 2\gamma^2 \leq \text{off}^2(A) - \frac{\text{off}^2(A)}{N}

  • Это дает нам: \text{off}^2(B) \leq \text{off}^2(A) \cdot (1 - \frac{1}{N})

  • Или в терминах \text{off}(A): \text{off}(B) \leq \sqrt{1 - \frac{1}{N}} \cdot \text{off}(A)

  • После N шагов (полного цикла) мы уменьшаем \text{off}(A) примерно в e^{1/2} \approx 1.65 раз:

\left(1 - \frac{1}{N}\right)^{\frac{N}{2}} \approx e^{-\frac{1}{2}}

  • Это показывает линейную сходимость метода. При этом локально (вблизи решения) сходимость становится квадратичной.

Jacobi: итоги

Метод Якоби был первым численным методом для вычисления собственных значений, предложенным в 1846 году.

  • Большая константа
  • Очень точный (высокая относительная погрешность для малых собственных значений)
  • Хорошие возможности параллелизации

Итоги по этой части

  • Множество алгоритмов для вычисления решения задачи на собственные значения:
    • QR
    • Разделяй и властвуй
    • Метод бисекции
    • Метод Якоби

Следующая лекция

  • Мы начинаем разреженную и/или структурированную численную линейную алгебру.

Questions?

from IPython.core.display import HTML
def css_styling():
    styles = open("./styles/custom.css", "r").read()
    return HTML(styles)
css_styling()
---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
Cell In [2], line 5
      3     styles = open("./styles/custom.css", "r").read()
      4     return HTML(styles)
----> 5 css_styling()

Cell In [2], line 3, in css_styling()
      2 def css_styling():
----> 3     styles = open("./styles/custom.css", "r").read()
      4     return HTML(styles)

File ~/miniconda3/envs/tensorflow/lib/python3.9/site-packages/IPython/core/interactiveshell.py:282, in _modified_open(file, *args, **kwargs)
    275 if file in {0, 1, 2}:
    276     raise ValueError(
    277         f"IPython won't let you open fd={file} by default "
    278         "as it is likely to crash IPython. If you know what you are doing, "
    279         "you can use builtins' open."
    280     )
--> 282 return io_open(file, *args, **kwargs)

FileNotFoundError: [Errno 2] No such file or directory: './styles/custom.css'