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()
直しました。
ありがとうございます。
ChatGPTに「QR分解」を説明してもらいました
[Q] QR分解とは何ですか? どのような目的で使用されますか? 初心者にもわかるように説明してください。
[A] QR分解は、線型代数学において行列を2つの行列の積の形に分解する手法です。QR分解では、与えられた行列を正規直交行列(Q行列)と上三角行列(R行列)の積に分解します。具体的に説明すると、与えられた行列をAとします。その行列AをQR分解すると、A = QRとなります。ここで、Qは正規直交行列(列ベクトル同士の内積が0でない場合は1、0である場合は0となる行列)であり、Rは上三角行列(対角線よりも上の要素が0となる行列)です。QR分解は、主に以下のような目的で使用されます:
1. 線型方程式の解の求め方の改善: QR分解を用いることで、線型方程式の解を効率的に求めることができます。QR分解を行った後、上三角行列Rを用いて連立方程式を解くことが容易になります。
2. 最小二乗法の最適化問題の解法: 最小二乗法は、データに対してモデルを当てはめる際に使用されます。QR分解は、最小二乗法の最適なパラメータを求める手法の1つとして使用されます。
3. 固有値と固有ベクトルの計算: QR分解は、行列の固有値と固有ベクトルを計算するための手法としても利用されます。QR分解によって上三角行列Rを得ることができ、固有値はRの対角要素に対応し、固有ベクトルはQの列ベクトルとなります。初心者にも分かりやすく説明すると、QR分解は行列を2つの特殊な形の行列に分解する方法です。この分解によって、線型方程式の解や最小二乗法の問題の解、固有値と固有ベクトルを求めることができます。QR分解は、数値計算や統計解析、機械学習などさまざまな分野で広く使用される重要な手法です。