Сегодня мы поговорим о матричных разложениях как общем инструменте
Основные матричные разложения в вычислительной линейной алгебре:
Мы уже вводили QR-разложение некоторое время назад, но теперь мы обсудим его более подробно.
В вычислительной линейной алгебре нам нужно решать различные задачи, например:
Для этого мы представляем матрицу как сумму и/или произведение матриц с более простой структурой, чтобы решать упомянутые задачи быстрее / в более стабильной форме.
Что такое более простая структура?
Мы уже встречали несколько классов матриц со структурой.
Для плотных матриц наиболее важными классами являются
План сегодняшней лекции - обсудить разложения одно за другим и указать: - Как вычислить конкретное разложение - Когда разложение существует - Что делается в реальной жизни (LAPACK).
A = P L U,
где P - матрица перестановки, L - нижняя треугольная матрица, U - верхняя треугольная.
A^{-1} f = (L U)^{-1} f = U^{-1} L^{-1} f,
и это сводится к решению двух линейных систем
L y = f, \quad U x = y
с нижней и верхней треугольными матрицами соответственно.
Вопрос: какова сложность решения этих линейных систем?
Если матрица эрмитова положительно определенная, т.е.
A = A^*, \quad (Ax, x) > 0, \quad x \ne 0,
то она может быть представлена в виде
A = RR^*,
где R - нижняя треугольная матрица.
Это нам понадобится для QR разложения.
A = Q R,
где Q - ортогональная (унитарная) по столбцам матрица, а R - верхняя треугольная матрица.
Размеры матриц: Q имеет размер n \times m, R имеет размер m \times m если n\geq m. Смотрите наш постер для визуализации QR разложения
QR разложение определено для любой прямоугольной матрицы.
Это разложение играет ключевую роль во многих задачах: - Вычисление ортогональных базисов в линейном пространстве - Используется в предварительной обработке для SVD - QR-алгоритм для вычисления собственных векторов и собственных значений (один из 10 самых важных алгоритмов 20-го века) основан на QR разложении - Решение переопределенных систем линейных уравнений (задача линейных наименьших квадратов)
\Vert A x - b \Vert_2 \rightarrow \min_x,
где A имеет размер n \times m, n \geq m.
A = Q R,
и используем уравнение для псевдообратной матрицы в случае матрицы A полного ранга:
x = A^{\dagger}b = (A^*A)^{-1}A^*b = ((QR)^*(QR))^{-1}(QR)^*b = (R^*Q^*QR)^{-1}R^*Q^*b = R^{-1}Q^*b.
таким образом, x можно найти из
R x = Q^*b
Теорема. Каждая прямоугольная матрица размера n \times m имеет QR разложение.
Существует несколько способов доказать это и вычислить его:
Если у нас есть представление вида
A = QR,
тогда A^* A = ( Q R)^* (QR) = R^* (Q^* Q) R = R^* R, матрица A^* A называется матрицей Грама, и её элементы являются скалярными произведениями столбцов матрицы A.
(A^* A y, y) = (Ay, Ay) = \Vert Ay \Vert^2 > 0, \quad y\not = 0.
Следовательно, A^* A = R^* R всегда существует.
Тогда матрица Q = A R^{-1} является унитарной:
(A R^{-1})^* (AR^{-1})= R^{-*} A^* A R^{-1} = R^{-*} R^* R R^{-1} = I.
Когда матрица n \times m не имеет полного ранга по столбцам, она называется матрицей неполного ранга.
QR разложение, однако, также существует.
Для любой матрицы неполного ранга существует последовательность матриц полного ранга по столбцам A_k такая, что A_k \rightarrow A (почему?).
Каждая A_k может быть разложена как A_k = Q_k R_k.
Множество всех унитарных матриц компактно, поэтому существует сходящаяся подпоследовательность Q_{n_k} \rightarrow Q (почему?), и Q^* A_k \rightarrow Q^* A = R, которая является треугольной.
A^* A = R^* R,
и
Q = A R^{-1}.
import numpy as np
n = 20
r = 9
a = [[1.0 / (i + j + 0.5) for i in range(r)] for j in range(n)]
a = np.array(a)
q, Rmat = np.linalg.qr(a)
e = np.eye(r)
print('Built-in QR orth', np.linalg.norm(np.dot(q.T, q) - e))
gram_matrix = a.T.dot(a)
Rmat1 = np.linalg.cholesky(gram_matrix)
#q1 = np.dot(a, np.linalg.inv(Rmat1.T))
q1 = np.linalg.solve(Rmat1, a.T).T
print('Via Cholesky:', np.linalg.norm(np.dot(q1.T, q1) - e))Built-in QR orth 7.668461578388947e-16
Via Cholesky: 0.9966594169095181Грам-Шмидт: 1. q_1 := a_1/\Vert a_1 \Vert 2. q_2 := a_2 - (a_2, q_1) q_1, \quad q_2 := q_2/\Vert q_2 \Vert 3. q_3 := a_3 - (a_3, q_1) q_1 - (a_3, q_2) q_2, \quad q_3 := q_3/\Vert q_3 \Vert 4. И так далее
Заметим, что преобразование от Q к A имеет треугольную структуру, так как из k-го вектора мы вычитаем только предыдущие. Это следует из того факта, что произведение треугольных матриц является треугольной матрицей.
Метод Грама-Шмидта может быть очень неустойчивым (т.е. полученные векторы не будут ортогональными, особенно если q_k имеет малую норму).
Это называется потерей ортогональности.
Существует решение, называемое модифицированным методом Грама-Шмидта. Вместо того, чтобы делать
q_k := a_k - (a_k, q_1) q_1 - \ldots - (a_k, q_{k-1}) q_{k-1}
мы делаем это пошагово. Сначала устанавливаем q_k := a_k и ортогонализуем последовательно:
q_k := q_k - (q_k, q_1)q_1, \quad q_k := q_{k} - (q_k,q_2)q_2, \ldots
В точной арифметике это одно и то же. В арифметике с плавающей точкой это абсолютно разные вещи!
Обратите внимание, что сложность составляет \mathcal{O}(nm^2) операций
R = Q^* A,
и нам нужно найти определенную ортогональную матрицу Q, которая приводит матрицу к верхнетреугольной форме.
- Для простоты мы будем искать матрицу n \times n такую, что
Q^* A = \begin{bmatrix} * & * & * \\ 0 & * & * \\ 0 & 0 & * \\ & 0_{(n-m) \times m} \end{bmatrix}
H_1 A = \begin{bmatrix} * & * & * \\ 0 & * & * \\ 0 & * & * \\ 0 & * & * \end{bmatrix}
Затем,
H_2 H_1 A = \begin{bmatrix} * & * & * \\ 0 & * & * \\ 0 & 0 & * \\ 0 & 0 & * \end{bmatrix},
где
H_2 = \begin{bmatrix} 1 & 0 \\ 0 & H'_2, \end{bmatrix}
и H'_2 - это матрица Хаусхолдера размера 3 \times 3.
Наконец,
H_3 H_2 H_1 A = \begin{bmatrix} * & * & * \\ 0 & * & * \\ 0 & 0 & * \\ 0 & 0 & 0 \end{bmatrix},
Вы можете попробовать реализовать это самостоятельно, это просто.
В реальности, поскольку это разложение плотной матрицы, алгоритм следует реализовывать в терминах блоков (почему?).
Вместо использования преобразования Хаусхолдера, мы используем блочное преобразование Хаусхолдера вида
H = I - 2UU^*,
где U^* U = I.
QR-разложение также может быть использовано для вычисления (численного) ранга матрицы, см. Rank-Revealing QR Factorizations and the Singular Value Decomposition, Y. P. Hong; C.-T. Pan
Это делается с помощью так называемого разложения, выявляющего ранг.
Оно основано на представлении
P A = Q R,
где P - матрица перестановки (она переставляет столбцы), а R имеет блочную форму
R = \begin{bmatrix} R_{11} & R_{12} \\ 0 & R_{22}\end{bmatrix}.
Цель состоит в том, чтобы найти P такую, что норма R_{22} мала, тогда вы можете определить численный ранг, взглянув на неё.
Оценка: \sigma_{r+1} \leq \Vert R_{22} \Vert_2 (проверьте, почему).
LU и QR разложения могут быть вычислены с помощью прямых методов за конечное количество операций.
А как насчет формы Шура и SVD?
Они не могут быть вычислены с помощью прямых методов (почему?) и могут быть вычислены только с помощью итерационных методов.
Хотя итерационные методы все равно имеют ту же сложность \mathcal{O}(n^3) в арифметике с плавающей точкой благодаря быстрой сходимости.
A = Q T Q^*,
где T - верхнетреугольная матрица, а Q - унитарная матрица, и это разложение дает собственные значения матрицы (они находятся на диагонали T).
QR алгоритм был независимо предложен в 1961 году Кублановской и Фрэнсисом.
Не путайте QR алгоритм и QR разложение!
QR разложение - это представление матрицы, тогда как QR алгоритм использует QR разложение для вычисления собственных значений!
A = Q T Q^*,
и перепишем его в форме
Q T = A Q.
Слева мы можем видеть QR разложение матрицы AQ.
Мы можем использовать это для вывода итерации с фиксированной точкой для формы Шура, также известной как QR алгоритм.
Мы можем записать итерационный процесс
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} = Q^*_{k+1} A Q_{k+1} = ( Q_{k+1}^* A = R_{k+1} Q^*_k) = R_{k+1} \widehat{Q}_k.
Таким образом, мы приходим к стандартной форме QR алгоритма.
Итоговые формулы записываются в классической форме QRRQ:
Повторяем итерации, пока A_k не станет достаточно треугольной (например, норма поддиагональной части достаточно мала).
Утверждение
Матрицы A_k унитарно подобны A
A_k = Q^*_{k-1} A_{k-1} Q_{k-1} = (Q_{k-1} \ldots Q_1)^* A (Q_{k-1} \ldots Q_1)
и произведение унитарных матриц является унитарной матрицей.
Сложность каждого шага составляет \mathcal{O}(n^3), если выполняется общее QR разложение.
Мы надеемся, что A_k будет очень близка к треугольной матрице для достаточно большого k.
import numpy as np
n = 40
a = [[1.0/(i - j + 0.5) for i in range(n)] for j in range(n)]
a = np.array(a)
niters = 10000
for k in range(niters):
q, rmat = np.linalg.qr(a)
a = rmat.dot(q)
print('Leading 3x3 block of a:')
print(a[:4, :4])Leading 3x3 block of a:
[[ 3.08733278e-01 3.08447671e+00 1.39093779e-02 -9.15979870e-02]
[-3.05714157e+00 2.90031089e-01 1.49436447e-01 -4.41658920e-02]
[-2.25478613e-21 -1.27296314e-20 5.50815465e-01 3.04169557e+00]
[-1.45411659e-20 -7.68894785e-21 -3.00877077e+00 5.09624202e-01]]Сходимость QR алгоритма происходит от наибольших собственных значений к наименьшим.
Для одного собственного значения требуется как минимум 2-3 итерации.
Каждый шаг состоит из одной QR факторизации и одного произведения матриц, в результате сложность \mathcal{O}(n^3).
В: означает ли это, что общая сложность \mathcal{O}(n^4)?
О: к счастью, нет.
Мы можем ускорить QR алгоритм, используя сдвиги, так как A_k - \lambda I имеет те же векторы Шура.
Мы обсудим эти детали позже
Всегда ли сходится QR алгоритм?
Приведите пример.
Для матрицы A = \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix}
мы имеем A_k = A.
В предыдущей лекции мы рассмотрели степенной метод, который представляет собой A^k v – приближение собственного вектора.
QR алгоритм вычисляет (неявно) QR-разложение матрицы A^k:
A^k = A \cdot \ldots \cdot A = Q_1 R_1 Q_1 R_1 \ldots = Q_1 Q_2 R_2 Q_2 R_2 \ldots (R_2 R_1) = \ldots (Q_1 Q_2 \ldots Q_k) (R_k \ldots R_1).
A = U \Sigma V^*.
A^* A = V^* \Sigma^2 V,
и/или
AA^* = U^* \Sigma^2 U
с помощью QR алгоритма, но это плохая идея (см. матрицы Грама).