2025-02-27

線形代数勉強：行列５分解 QR分解(2) Gram-Schmidtで軽く実装

（以下は自習資料です）

QR分解のコードを書いてみる。とりあえずグラム・シュミットの直交化を適用して正規直交行列 $ Q $ を求め、それを利用して $R$ を求める。

全コードはこちら

def qr_decomposition_gram_schmidt(
    A: ND,
) -> Tuple[ND, ND]:
    """与えられた正方行列aをGram-Schmidt法でQR分解する.

    Args:
        A (numpy.ndarray): 正方行列。

    Returns:
        Tuple[numpy.ndarray, numpy.ndarray]: 正規直交行列Qと上三角行列Rのタプル。
    """
    assert A.shape[0] == A.shape[1]
    ncols = A.shape[1]

    Q = numpy.zeros_like(A)
    R = numpy.zeros_like(A)

    for i in range(0, ncols):
        v_i = A[:, i].copy()  # A[:, i] を直接書き換えないようにコピー
        for j in range(0, i):
            v_i -= (Q[:, j] @ v_i) * Q[:, j]  # 直交化
        R[i, i] = norm(v_i)
        Q[:, i] = v_i / norm(v_i)  # 正規化
        for k in range(i + 1, ncols):
            R[i, k] = Q[:, i] @ A[:, k]

    return Q, R

なお、アルゴリズムの途中で、列ベクトルを求めるごとに毎回正規化していくコードになっているが、これを行わないと、正しい結果にならないようだ。理論上は正規化しなくても列ベクトル同士の直交性には影響しないはずだと思うが、ChatGPTによれば、列ベクトルごとのノルムの差が大きくなると計算に誤差がでるとのことのようだ。この問題を軽減するために修正Gram-Schmidt法というものがあるらしい。$R$の計算を同時にやってしまえばこれは不要なので実装しないことにした。