Десятое занятие: реализация QR-разложения

Мы используем QR-разложение при вычислении собственных значений и решении задачи наименьших квадратов. Это важный компонент численной линейной алгебры.

«В численной линейной алгебре есть один алгоритм, который важнее других — это QR-разложение», — пишет Трефетен на странице 48.

Вспомним, что в QR-алгоритме, который мы видели на предыдущем занятии, используется QR-разложение, но не следует путать эти два понятия.

В NumPy

import numpy as np

np.set_printoptions(suppress=True, precision=4)

n = 5
A = np.random.rand(n, n)
npQ, npR = np.linalg.qr(A)

Проверяем, является ли Q ортогональной матрицей:

np.allclose(np.eye(n), npQ @ npQ.T), np.allclose(np.eye(n), npQ.T @ npQ)
# (True, True)

Проверяем, является ли R треугольной матрицей.

npR
array([[-0.8524, -0.7872, -1.1163, -1.2248, -0.7587],
       [ 0.    , -0.9363, -0.2958, -0.7666, -0.632 ],
       [ 0.    ,  0.    ,  0.4645, -0.1744, -0.3542],
       [ 0.    ,  0.    ,  0.    ,  0.4328, -0.2567],
       [ 0.    ,  0.    ,  0.    ,  0.    ,  0.1111]])

Когда вектор b проецируется на прямую a, его проекция p является частью вектора b вдоль прямой a. Рассмотрим интерактивную диаграмму из раздела «Проекция» на сайте «Погружение в математику».

Источник: погружение в математику

Вот как выглядит проекция вектора на плоскость:

Источник: линейная алгебра с точки зрения наименьших квадратов регрессии

При проецировании вектора b на прямую a его проекция p представляет собой некоторое кратное a. Пусть p = $\hat{x}a$, где $\hat{x}$ — скаляр.

Ортогональность

Ключевым моментом в проекции является ортогональность: линия от b до p перпендикулярна a.

Это означает: $a \cdot (b - \hat{x}a) = 0$

Поэтому: $\hat{x} = \frac{a \cdot b}{a \cdot a}$

Грам-Шмидт

Классический (неустойчивый) Грам-Шмидт

Для каждой колонки j вычисляется одна проекция: $v_j = P_ja_j$

где $P_j$ ортогонален пространству $q_1, ..., q_{j-1}$.

def cgs(A):
    m, n = A.shape
    Q = np.zeros([m,n], dtype=np.float64)
    R = np.zeros([n,n], dtype=np.float64)
    for j in range(n):
        v = A[:,j]
        for i in range(j):
            R[i,j] = np.dot(Q[:,i], A[:,j])
            v = v - (R[i,j] * Q[:,i])
        R[j,j] = np.linalg.norm(v)
        Q[:, j] = v / R[j,j]
    return Q, R

Q, R = cgs(A)

np.allclose(A, Q @ R)
# True

Проверяем, что Q является унитарной матрицей.

np.allclose(np.eye(len(Q)), Q.dot(Q.T))
# True

np.allclose(npQ, -Q)
# True

R
array([[ 0.02771,  0.02006, -0.0164 , ...,  0.00351,  0.00198,  0.00639],
       [ 0.     ,  0.10006, -0.00501, ...,  0.07689, -0.0379 , -0.03095],
       [ 0.     ,  0.     ,  0.01229, ...,  0.01635,  0.02988,  0.01442],
       ...,
       [ 0.     ,  0.     ,  0.     , ...,  0.     , -0.     , -0.     ],
       [ 0.     ,  0.     ,  0.     , ...,  0.     ,  0.     , -0.     ],
       [ 0.     ,  0.     ,  0.     , ...,  0.     ,  0.     ,  0.     ]])

Грам-Шмидт напоминает вам об итерации Арнольди, которая используется для преобразования матрицы в форму Хессенберга, поскольку она также является структурированным процессом ортогонализации.

Улучшенный Грам-Шмидт

Классический (неустойчивый) Грам-Шмидт: для каждой колонки j вычисляется одна проекция: $v_j = P_ja_j$

где $P_j$ ортогонален пространству $q_1, ..., q_{j-1}$.

Улучшенный Грам-Шмидт: для каждой колонки j вычисляются n-1 проекций: $P_j = P_{\perp q_{j-1}\cdots\perp q_{2}\perp q_{1}}$

import numpy as np
n = 3
A = np.random.rand(n,n).astype(np.float64)

def cgs(A):
    m, n = A.shape
    Q = np.zeros([m,n], dtype=np.float64)
    R = np.zeros([n,n], dtype=np.float64)
    for j in range(n):
        v = A[:,j]
        for i in range(j):
            R[i,j] = np.dot(Q[:,i], A[:,j])
            v = v - (R[i,j] * Q[:,i])
        R[j,j] = np.linalg.norm(v)
        Q[:, j] = v / R[j,j]
    return Q, R

def mgs(A):
    V = A.copy()
    m, n = A.shape
    Q = np.zeros([m,n], dtype=np.float64)
    R = np.zeros([n,n], dtype=np.float64)
    for i in range(n):
        R[i,i] = np.linalg.norm(V[:,i])
``` ```
Q[:, i] = V[:, i] / R[i, i]
for j in range(i, n):
    R[i, j] = np.dot(Q[:, i], V[:, j])
    V[:, j] = V[:, j] - R[i, j]*Q[:, i]
return Q, R

В данном фрагменте кода на языке Python происходит вычисление матриц $Q$ и $R$ с использованием матриц $V$ и $A$.

Для вычисления элементов матрицы $Q$ используется формула:

$Q_{ij} = \frac{V_{ij}}{R_{ii}}$, где $i$ — номер строки, а $j$ — номер столбца.

Затем для вычисления элементов матрицы $R$ используется следующая формула:

$R_{ij} = (Q * V)_{ij}$, где $i$ и $j$ имеют те же значения, что и в формуле для $Q$.

После этого элементы матрицы $V$ обновляются по формуле:

$V_{ij} = V_{ij} - (R * Q)_{ij}$.

Наконец, матрицы $Q$ и $R$ возвращаются из функции.

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