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)
$\displaystyle \left[\begin{matrix}3 & 2 & \frac{1}{3}\\0 & 1 & \frac{10}{3}\\0 & 0 & \frac{4}{3}\end{matrix}\right]$
下記は誤植と思われますので、ご確認ください。
(6) 式で、a2= → a3=
(7) 式の 2 行 3 列目の成分で、q1^T a3 → q2^T a3
SymPyプログラムで、sympy.init_printing() → init_printing()
直しました。
ありがとうございます。