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)

 

$\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]$

 

コメント

  1. HNaito より:

    下記は誤植と思われますので、ご確認ください。
    (6) 式で、a2= → a3=
    (7) 式の 2 行 3 列目の成分で、q1^T a3 → q2^T a3
    SymPyプログラムで、sympy.init_printing() → init_printing()