QR分解

QR分解

　$m\times n$ 行列 $A$ を $m\times m$ 直交行列 $Q$ と $m\times n$ 上三角行列 $R$ の積に分解することを QR分解 (QR decomposition) とよび、数値線型代数において重要な役割を担います。NumPy には QR分解を実行する numpy.linalg.qr() が実装されています。

QR分解の手順

　最初に QR 分解の理論について簡単に学んでおきましょう。
　記述を簡単にするために、$3\times 3$ 正方行列 $A$ を考えます。
　前回記事で学んだグラム・シュミット直交化の手順を再掲します。
　
\[\begin{align*}&\boldsymbol{u}_1=\boldsymbol{a}_1,\quad \boldsymbol{q}_1=\frac{\boldsymbol{u}_1}{\parallel\boldsymbol{u}_1\parallel}\tag{1}\\[6pt]&\boldsymbol{u}_2=\boldsymbol{a}_2-(\boldsymbol{q}_1^T\boldsymbol{a}_2)\boldsymbol{q}_1,\quad \boldsymbol{q}_2=\frac{\boldsymbol{u}_2}{\parallel\boldsymbol{u}_2\parallel}\tag{2}\\[6pt]&\boldsymbol{u}_3=\boldsymbol{a}_3-(\boldsymbol{q}_1^T\boldsymbol{a}_3)\boldsymbol{q}_1-(\boldsymbol{q}_2^T\boldsymbol{a}_3)\boldsymbol{q}_2,\quad \boldsymbol{q}_3=\frac{\boldsymbol{u}_3}{\parallel\boldsymbol{u}_3\parallel}\tag{3}\end{align*}\]
　したがって、$\{\boldsymbol{a}_i\}$ は
　
\[\begin{align*}&\boldsymbol{a}_1=\parallel\boldsymbol{u}_1\parallel\boldsymbol{q}_1\tag{4}\\[6pt]&\boldsymbol{a}_2=\parallel\boldsymbol{u}_2\parallel\boldsymbol{q}_2(\boldsymbol{q}_1^T\boldsymbol{a}_2)\boldsymbol{q}_1\tag{5}\\[6pt]&\boldsymbol{a}_3=\parallel\boldsymbol{u}_3\parallel\boldsymbol{q}_3+(\boldsymbol{q}_1^T\boldsymbol{a}_3)\boldsymbol{q}_1+(\boldsymbol{q}_2^T\boldsymbol{a}_3)\boldsymbol{q}_2\tag{6}\end{align*}\]
となりますが、行列を使って以下のように表すことができます。
　
\[\begin{bmatrix}\boldsymbol{a}_1&\boldsymbol{a}_2&\boldsymbol{a}_3\end{bmatrix}=\begin{bmatrix}\boldsymbol{q}_1&\boldsymbol{q}_2&\boldsymbol{q}_3\end{bmatrix}\begin{bmatrix}\parallel\boldsymbol{u}_1\parallel&\boldsymbol{q}_1^T\boldsymbol{a}_2&\boldsymbol{q}_1^T\boldsymbol{a}_3\\0&\parallel\boldsymbol{u}_2\parallel&\boldsymbol{q}_2^T\boldsymbol{a}_3\\0&0&\parallel\boldsymbol{u}_3\parallel\end{bmatrix}\tag{7}\]
　現在ではハウスホルダー変換やギブンス回転などの高度な手法が確立されているので、QR分解にグラム・シュミット分解が用いられることはほとんどありませんが、中級プログラミングの練習として良い題材なので、意欲のある人は実装に挑戦してみてください。

numpy.linalg.qr()

　numpy.linalg.qr() は行列を受け取って QR分解 を実行します。

# numpy_qr_decomposition

# In[1]

import numpy as np
from scipy import linalg

# 小数点以下3桁で表示
# 指数表記は用いずに常に小数表示
np.set_printoptions(precision=3, suppress=True)

# 乱数シードを0に設定
np.random.seed(0)

# -9～9のランダム要素をもつ4×6行列を生成
A = np.random.randint(-9, 10, (4, 6))

# QR分解を実行
Q, R = np.linalg.qr(A)

print("行列A:\n{}\n".format(A))
print("Q:\n{}\n".format(Q))
print("R:\n{}".format(R))

'''
行列A:
[[ 3  6 -9 -6 -6 -2]
 [ 0  9 -5 -3  3 -8]
 [-3 -2  5  8 -4  4]
 [-1  0  7 -4  6  6]]

Q:
[[-0.688 -0.232 -0.021  0.687]
 [-0.000 -0.945 -0.060 -0.321]
 [ 0.688 -0.188 -0.335  0.616]
 [ 0.229 -0.133  0.940  0.214]]

R:
[[-4.359 -5.506 11.241  8.718  2.753  5.506]
 [ 0.000 -9.523  4.947  3.255 -1.487  6.478]
 [ 0.000  0.000  5.400 -6.130  6.927  4.826]
 [ 0.000  0.000  0.000  0.910 -6.263  4.943]]
'''

MutableDenseMatrix.QRdecomposition()

　SymPy の Matrix オブジェクトには QR分解を実行する QRdecomposition()メソッドが備わっています。

# sympy_qr_decomposition

# In[1]

from sympy import *

# LaTeX形式で結果を表示
# バージョンによっては必要ありません
init_printing()

A = Matrix([[2, 2,  2],
            [2, 1,  0],
            [1, 0, -3]])

# QR分解
Q, R = A.QRdecomposition()

display(Q, R)