1. ホーム
  2. プログラマーのための線形代数
  3. 行列の分解
  4. QR分解

QR分解

QR分解は、行列の分解テクニックの一つで、線型最小二乗問題や固有値問題を解くときに使われます。当ページでは、このQR分解について知っておきたいことをまとめて確認することができます。具体的には以下のことがわかります

当ページでわかること

  • QR分解とは
  • QR分解のやり方
  • PythonでQR分解

それでは見ていきましょう。

目次

QR分解とは

QR 分解は \(m \times n\) の行列 \(A\)を \(m \times m\) の直交行列 \(Q\) と \(m \times n\) の上三角行列 \(R\) の積に分解するものです。なお直交行列とは、それぞれの列、またはそれぞれの行がお互いに直交した単位ベクトル(長さが 1 のベクトル)である行列のことです。

たとえば以下の行列 A は、ご覧の通り、直交行列 Q と上三角行列 \(R\) の積で表すことが可能です。

\[\begin{eqnarray}
\overset{A}{
\begin{pmatrix}
1 & 2 \\
3 & 4 \\
5 & 6
\end{pmatrix}
}
=
\overset{Q}{
\begin{pmatrix}
-0.17 & 0.90 & 0.40 \\
-0.51 & 0.28 & -0.82 \\
-0.85 & -0.35 & 0.40
\end{pmatrix}}
\overset{R}{
\begin{pmatrix}
-5.92 & -7.44 \\
0 & 0.83 \\
0 & 0
\end{pmatrix}}
\end{eqnarray}\]

このようにQR分解は、LU分解と違って、もともとの行列が正方行列でなくても可能であり、もともとの行列を直交行列と上三角行列に分解するものです。このQR分解を行うことの利点は次の通りです。

QR分解の利点

  • 線型最小二乗問題を解けるようになる。
  • 固有値問題を解けるようになる。

QR分解のやり方

QR分解のやり方には次の 3 つの方法があります。

  • グラム・シュミット分解
  • ハウスホルダー変換
  • ギブンス回転

グラム・シュミット分解は、直交行列の数値の精度が不安定ですが、もっとも計算がしやすいという利点があります。そのため、線形代数ライブラリが利用できない場合に、自分自身でQR分解のプログラムを実装する際に便利です。

ハウスホルダー変換は、数値的にもっとも安定した方法ですが、計算に手間がかかります。そのためこの方法をプログラムに実装するのは動作が重くなるため向いていません。

ギブンス回転は、数値的に安定していますが、アルゴリズムをプログラムとして実装するのに最も手間がかかります。ただし、この方法で実装したプログラムは動作が高速です。そのため厳密な数値計算ライブラリでは、この計算方法でQR分解の関数が実装されている場合が多いです。

ここではグラム・シュミット分解についてのみ解説することにします。

グラム・シュミット分解

あらためて、QR 分解は \(m \times n\) の行列 \(A\)を、\(m \times m\) の直交行列 \(Q\) と \(m \times n\) の上三角行列 \(R\) の積に分解するものです。

それでは以下の行列を例にやってみましょう。

\[\begin{eqnarray}
A=
\begin{pmatrix}
1 & 1 & -1\\
1 & -1 & 1\\
1 & 2 & 3
\end{pmatrix}
\end{eqnarray}\]

まずはグラム・シュミットの正規直交化法を使って、\(A\) を直交行列にします。 \(A\) の列ベクトル \(a_1,a_2,a_3\) が 3 次元の場合、直交行列 \(Q\) のそれぞれの列ベクトル \(q_1, q_2, q_3\) は次のように求められます(ベクトル \(a\) の\(L_2\)ノルムを \(\|a_1 \|\) とします。そしてベクトル \(a\) と \(q\) の内積が \(a \cdot q\) です)

  • \(q_1 = \dfrac{a_1}{\|a_1\|}\)
  • \(v_2 = a_2-(a_2 \cdot q_1)q_1\), \(q_2=\dfrac{v_2}{\|v_2\|}\)
  • \(v_3 =a_3-(a_3 \cdot q_1)q_1-(a_3 \cdot q_2)q_2\), \(q_3= \dfrac{q_3}{|q_3|}\)

実際に行列 \(A\) の直交行列を求めてみましょう。

まず、\(q_1\) は次のように求められます。

\[
q_1
=
\dfrac{(1, 1, 1)}{\sqrt{1^2+1^2+1^2}}
=
\dfrac{1}{\sqrt{3}}(1, 1, 1)
\]

次に \(q_2\) は次のように求められます(内積の計算は省きます)

\[\begin{eqnarray}
v_2 = a_2 – (a_2 \cdot q_1)q_1
&=&
(1, -1, 2)

\left(\dfrac{2}{\sqrt{3}}\right)
\dfrac{1}{\sqrt{3}}(1,1,1)\\
&=&
\dfrac{1}{3}(1,-5,4)
\end{eqnarray}\]

\[\begin{eqnarray}
q_2
=
\dfrac{v_2}{\|v_2\|}
&=&
\dfrac
{
\dfrac{1}{3}(1,-5,4)
}
{
\dfrac{1}{3}\sqrt{(1^2+\{-5\}^2+4^2)}
}\\
&=&
\dfrac{1}{\sqrt{42}}(1, -5, 4)
\end{eqnarray}\]

最後に \(q_3\) は次のように求められます。

\[\begin{eqnarray}
v_3 = a_3-(a_3 \cdot q_1)q_1-(a_3 \cdot q_2)q_2
&=&
a_3-\sqrt{3}q_1 – \dfrac{6}{\sqrt{42}}q_2\\
&=&
\dfrac{1}{7}(-3,1,2)
\end{eqnarray}\]

\[\begin{eqnarray}
q_3
=
\dfrac{v_3}{\|v_3\|}
&=&
\dfrac
{
\dfrac{1}{7}(-3,1,2)
}
{
\dfrac{1}{7}\sqrt{(-3^2+1^2+2^2)}
}\\
&=&
\dfrac{1}{\sqrt{14}}(-3, 1, 2)
\end{eqnarray}\]

以上のことから行列 \(A\) の直交行列 \(Q\) は以下の通りになります。

\[\begin{eqnarray}
Q
=
\begin{pmatrix}
\frac{1}{\sqrt{3}} & \frac{1}{\sqrt{42}} & \frac{-3}{\sqrt{14}}\\
\frac{1}{\sqrt{3}} & \frac{-5}{\sqrt{42}} & \frac{1}{\sqrt{14}}\\
\frac{1}{\sqrt{3}} & \frac{4}{\sqrt{42}} & \frac{2}{\sqrt{14}}\
\end{pmatrix}
\end{eqnarray}\]

これで直交行列を求められました。次に上三角行列 \(R\) を求めるのですが、その際は、直交行列の次の性質を使います(\(I\) は単位行列です)

\[
Q^T \cdot Q = Q \cdot Q^T = I
\]

この性質より、\(R\) は次のように求められます。

\[Q^TA = Q^TQR = R\]

\[\begin{eqnarray}
R
=
Q^TA
=
\begin{pmatrix}
\frac{3}{\sqrt{3}} & \frac{2}{\sqrt{3}} & \frac{6}{\sqrt{3}}\\
0 & \frac{14}{\sqrt{42}} & \frac{3}{\sqrt{42}}\\
0 & 0 & \frac{5}{\sqrt{14}}\
\end{pmatrix}
\end{eqnarray}\]

以上がQR分解です。

PythonでQR分解

Python では QR 分解は、NumPy の linalg.qr() 関数として実装されています。デフォルトでは、この関数は、次元削減された Q と R 行列を返します。mode 引数に “complete” と渡すと、\(n \times m\) の \(Q\) と \(m \times n\) の \(R\) を得られます。ただし、多くのアプリケーションでは、そこまでは必要とされていません。

以下がコード例です。

In [1]:
# NumPy と SciPy のインポート
import numpy as np
from scipy.linalg import qr

# 長方行列を定義
A = np.array([
    [1,2,],
    [3,4,],
    [5,6]])
print(A)
[[1 2]
 [3 4]
 [5 6]]
In [2]:
# QR分解
Q,R = qr(A)
print(Q, "complete")
[[-0.16903085  0.89708523  0.40824829]
 [-0.50709255  0.27602622 -0.81649658]
 [-0.84515425 -0.34503278  0.40824829]] complete
In [3]:
print(R)
[[-5.91607978 -7.43735744]
 [ 0.          0.82807867]
 [ 0.          0.        ]]
In [4]:
# 元々の行列の再作成
B= Q @ R 
print(B)
[[1. 2.]
 [3. 4.]
 [5. 6.]]