QR分解

QR分解

QR分解

 $m\times n$ 行列 $A$ を $m\times m$ 直交行列 $Q$ と $m\times n$ 上三角行列 $R$ の積に分解することを QR分解 (QR decomposition) とよび、数値線型代数において重要な役割を担います。NumPy には QR分解を実行する numpy.linalg.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}_2=\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}_1^T\boldsymbol{a}_3\\
0&0&\parallel\boldsymbol{u}_3\parallel
\end{bmatrix}\tag{7}\]
 現在ではハウスホルダー変換やギブンス回転などの高度な手法が確立されているので、QR分解にグラム・シュミット分解が用いられることはほとんどありませんが、中級プログラミングの練習として良い題材なので、意欲のある人は実装に挑戦してみてください。
 

numpy.linalg.qr()

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

# SLA_063

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()メソッドが備わっています。

# SLA_064

from sympy import *

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

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

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

display(Q, R)
$\displaystyle \left[\begin{matrix}\frac{2}{3} & \frac{2}{3} & - \frac{1}{3}\\\frac{2}{3} & - \frac{1}{3} & \frac{2}{3}\\\frac{1}{3} & - \frac{2}{3} & - \frac{2}{3}\end{matrix}\right]$
$\displaystyle \left[\begin{matrix}3 & 2 & \frac{1}{3}\\0 & 1 & \frac{10}{3}\\0 & 0 & \frac{4}{3}\end{matrix}\right]$