Ctrl+K

QR алгоритм

QR алгоритм#

  • QR алгоритм был предложен в 1961 г. независимо В. Н. Кублановской и Ј. Francis’ом. Статью про историю этого алгоритма и его авторов можно прочитать тут. https://www.atm.org.uk/write/MediaUploads/Resources/Mid_Plenary_FrancisGolub.pdf

  • Не путайте QR алгоритм и QR разложение!

  • QR разложение - это представление матрицы в виде произведения двух матриц, а QR алгоритм использует QR разложение для вычисления разложения Шypa.

  • Является «золотым стандартом» для вычисления собственных значений

Теория QR алгоритма#

Рассмотрим выражение

\[ A=Q T Q^* \]

и перепишем его в виде

\[ Q T=A Q \]

Слева замечаем QR разложение матрицы \(A Q\).

Запишем следующий итерационный процесс

\[ Q_{k+1} R_{k+1}=A Q_k, \quad Q_{k+1}^* A=R_{k+1} Q_k^* \]

Введём новую матрицу

\[ A_k=Q_k^* A Q_k=Q_k^* Q_{k+1} R_{k+1}=\widehat{Q}_k R_{k+1} \]

тогда аппроксимация для \(A_{k+1}\) имеет вид

\[ A_{k+1}=Q_{k+1}^* A Q_{k+1}=\left(Q_{k+1}^* A=R_{k+1} Q_k^*\right)=R_{k+1} \widehat{Q}_k . \]

Итак, мы получили стандартную форму записи QR алгоритма. Финальные формулы обычно записывают в QRRQ-форме:

  1. Инициализируем \(A_0=A\).

  2. Вычислим QR разложение матрицы \(A_k: A_k=Q_k R_k\).

  3. Обновим аппроксимацию \(A_{k+1}=R_k Q_k\).

Продолжаем итерации пока \(A_k\) не станет достаточно треугольной (например, норма подматрицы под главной диагональю не станет достаточно мала).

Утверждение Матрицы \(A_k\) унитарно подобны матрице \(A\)

\[ A_k=Q_{k-1}^* A_{k-1} Q_{k-1}=\left(Q_{k-1} \ldots Q_1\right)^* A\left(Q_{k-1} \ldots Q_1\right) \]

а произведение унитарных матриц - унитарная матрица. Сложность одной итерации \(\mathcal{O}\left(n^3\right)\), если используется \(\mathrm{QR}\) разложение для общего случая. Мы ожидаем, что \(A_k\) будет очень близка к треугольной матрице для достаточно большого \(k\).

import numpy as np
n = 4
a = [[1.0/(i + j + 0.5) for i in range(n)] for j in range(n)]
niters = 200
for k in range(niters):
    q, rmat = np.linalg.qr(a)
    a = rmat.dot(q)
print('Leading 3x3 block of a:')
display(a[:3, :3])

Leading 3x3 block of a:

array([[ 2.41052440e+000,  2.69817894e-018, -4.95361754e-017],
       [ 2.42500623e-168,  3.49984625e-001,  1.75163388e-017],
       [ 0.00000000e+000,  6.56745067e-273,  1.53236733e-002]])

Реализация QR-алгоритма#

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

def qr_algorithm(A, num_iter, eps):
    T = A.copy()
    U = np.eye(A.shape[0])
    conv = [(T, U)]
    for i in range(num_iter):
        Q, R = np.linalg.qr(T)
        T = R @ Q
        U = U @ Q
        conv.append((T, U))
        if np.sum(np.abs(np.tril(T, k=-1))) < eps:
            break
    return T, U, conv[1:]

n = 7
A = np.random.randn(n, n)
# A = A.T @ A
A = A + A.T
true_eigvals, true_eigvec = np.linalg.eig(A)
print(true_eigvals)

[-4.90241018  6.15117923 -3.18614889 -2.09635507 -0.34667853  3.92398239
  2.59969224]

T, U, conv = qr_algorithm(A, 2000, 1e-6)
print(np.linalg.norm(A - U @ T @ U.T))

7.410718418419184e-14

plt.spy(T, markersize=10, precision=1e-6)

<matplotlib.lines.Line2D at 0x7f925d409760>
../../_images/2816ac8d41c42bdc4b16cd02f497ec2b44f8b95e0618493d28495d83cb5cc4d0.png 
plt.figure(figsize=(10, 7))
conv_qr = np.array([np.sum(np.abs(np.tril(T, k=-1))) for T, U in conv])
plt.plot(conv_qr)
plt.yscale("log")
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.xlabel("\# iterations", fontsize=20)
plt.ylabel("Sum of lower triangular part in modulos", fontsize=20)
plt.grid(True)
../../_images/f572cdec5b8fdcb72967889321238dd4966f5728b0be82adf337ad6ee00f6a34.png 
U.round(4), true_eigvec.round(4)

(array([[-0.4441, -0.3348, -0.1244,  0.4476, -0.0476, -0.6016, -0.3326],
        [ 0.6661,  0.1351,  0.3278, -0.1021, -0.2134, -0.5688, -0.2261],
        [-0.3834,  0.0422,  0.6656, -0.2172, -0.1887,  0.2568, -0.5093],
        [ 0.3495, -0.2065,  0.0116,  0.6566, -0.4027,  0.4643, -0.1616],
        [ 0.282 , -0.4769, -0.281 , -0.2959,  0.424 ,  0.1729, -0.5629],
        [-0.0987,  0.3416, -0.5925, -0.2812, -0.5804,  0.0219, -0.3257],
        [-0.0258, -0.6935,  0.0624, -0.3798, -0.4876, -0.0511,  0.3603]]),
 array([[-0.3348, -0.4441,  0.4476,  0.6016, -0.3326,  0.1244,  0.0476],
        [ 0.1351,  0.6661, -0.1021,  0.5688, -0.2261, -0.3278,  0.2134],
        [ 0.0422, -0.3834, -0.2172, -0.2568, -0.5093, -0.6656,  0.1887],
        [-0.2065,  0.3495,  0.6566, -0.4643, -0.1616, -0.0116,  0.4027],
        [-0.4769,  0.282 , -0.2959, -0.1729, -0.5629,  0.281 , -0.424 ],
        [ 0.3416, -0.0987, -0.2812, -0.0219, -0.3257,  0.5925,  0.5804],
        [-0.6935, -0.0258, -0.3798,  0.0511,  0.3603, -0.0624,  0.4876]]))

Ускорение QR алгоритма#

Матрица \(A\) имеет верхне-гессенбергову форму, если \(a_{i j}=0\), при \(i \geq j+2\).

\[\begin{split} H=\left[\begin{array}{lllll} * & * & * & * & * \\ * & * & * & * & * \\ 0 & * & * & * & * \\ 0 & 0 & * & * & * \\ 0 & 0 & 0 & * & * \end{array}\right] \end{split}\]

С помощью отражений Хаусхолдера можно привести любую матрицу к верхне-гессенберговой форме:

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

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

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

  • Также верхне-гессенбергова форма матрицы сохраняется после выполнения одной итерации QR алгоритма (проверим ниже).

Матрица Хаусхолдера - это матрица вида

\[ H \equiv H(v)=I-2 v v^*, \]

где \(v-n \times 1\) вектор и \(v^* v=1\). Покажите, что \(H\) унитарна и Эрмитова. Это также матрица отражения относительно плоскости с нормалью \(v\):

\[ H x=x-2\left(v^* x\right) v \]

image.png

Свойства матрицы Хаусхолдера#

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

\[\begin{split} H\left[\begin{array}{c} \times \\ \times \\ \times \\ \times \end{array}\right]=\left[\begin{array}{l} \times \\ 0 \\ 0 \\ 0 \end{array}\right] \end{split}\]

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

Пусть \(e_1=(1,0, \ldots, 0)^T\), тогда нам надо найти такой вектор \(v\) что

\[ H x=x-2\left(v^* x\right) v=\alpha e_1, \]

где \(\alpha\) неизвестная константа. В силу унитарной инвариантности \(\|\cdot\|_2\) мы получим

\[ \|x\|_2=\|H x\|_2=\left\|\alpha e_1\right\|_2=|\alpha| . \]

и

\[ \alpha=\pm\|x\|_2 \]

Также, можем выращить \(v\) из равенства \(x-2\left(v^* x\right) v=\alpha e_1\):

\[ v=\frac{x-\alpha e_1}{2 v^* x} \]

Умножив последнее выражение на \(x^*\) получим

\[\begin{split} \begin{gathered} x^* x-2\left(v^* x\right) x^* v=\alpha x_1 ; \\ \|x\|_2^2-2\left(v^* x\right)^2=\alpha x_1 \\ \left(v^* x\right)^2=\frac{\|x\|_2^2-\alpha x_1}{2} . \end{gathered} \end{split}\]

Итак, \(v\) существует и равна

\[ v=\frac{x \pm\|x\|_2 e_1}{2 v^* x}=\frac{x \pm\|x\|_2 e_1}{\pm \sqrt{2\left(\|x\|_2^2 \mp\|x\|_2 x_1\right)}} \]

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

  • Если матрица \(A\) симметричная (эрмитова), то \(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\) симметричная трёхдиагональная матрица , то эта форма сохраняется.

  • Давайте это проверим!

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style('white')

#Generate a random tridiagonal matrix
n = 20
d = np.random.randn(n)
sub_diag = np.random.randn(n-1)

mat = np.diag(d) + np.diag(sub_diag, -1) + np.diag(sub_diag, 1)
plt.spy(mat)
plt.title("Original matrix", fontsize=16)
q, r = np.linalg.qr(mat)
plt.figure()
b = r.dot(q)
b[abs(b) <= 1e-12] = 0
plt.spy(b)
plt.title("After one iteration of QR algorithm", fontsize=16)
#plt.figure()
#plt.imshow(np.abs(r.dot(q)))
b[0, :]

array([-0.96810585,  1.51121923,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ])
../../_images/da99cde8d103e37ac199681341acb72e33028afb454ad70a7dfe4a245fcbc16c.png ../../_images/9add18aff186a0665ca0c3e4aab0cb23aad0aee9d4420256516454449850ff12.png

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

  • Работая с трёхдиагональной формой, вам не нужно вычислять матрицу \(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-шаr.

Теорема о неявном QR алгоритме

  • Bce реализации неявного QR алгоритма основаны на следующей теореме

Теорема.

Пусть

\[ Q^* A Q=H \]

верхне-гессенбергова форма матрицы. Тогда первый столбец матрицы \(Q\) определяет все остальные её столбцы. Он может быть найден из следующего уравнения

\[ A Q=Q H \]

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

  • Сходимость QR алгоритма - непростой вопрос (см. Е.Е. Тыртышников «Краткий курс численного анализа»)! Итог. Если у нас есть разложение вида

\[\begin{split} A=X \Lambda X^{-1}, \quad A=\left[\begin{array}{ll} A_{11} & A_{12} \\ A_{21} & A_{22} \end{array}\right] \end{split}\]

и

\[\begin{split} \Lambda=\left[\begin{array}{cc} \Lambda_1 & 0 \\ 0 & \Lambda_2 \end{array}\right], \quad \lambda\left(\Lambda_1\right)=\left\{\lambda_1, \ldots, \lambda_m\right\}, \lambda\left(\Lambda_2\right)=\left\{\lambda_{m+1}, \ldots, \lambda_r\right\} \end{split}\]

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

\[ \left\|A_{21}^{(k)}\right\| \leq C q^k, \quad q=\left|\frac{\lambda_{m+1}}{\lambda_m}\right|, \]

где \(m\) размер матрицы \(\Lambda_1\)

Таким образом, нам нужно увеличить зазор между \(\Lambda_1\) и \(\Lambda_2\). Это можно сделать с помощью QR алгоритма со сдвигами.

Модификации 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\)-ое большее по модулю собственное значение. Если сдвиг близок к собственному вектору, сходимость более быстрая.

  • Существуют различные стратегии выбора сдвигов.

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

Содержание