前書き

記事の概要

今回はQR法(固有値を求めるアルゴリズム)に欠かせない QR分解 について解説していきたいと思います。
この記事には

数式によるアルゴリズムの解説
QR分解に関する数学的な話
プログラム実装例

が含まれています。それぞれ自分に必要な部分を参照してもらえればいいなと思っています。

前提知識

大学1年生で勉強する「線形代数学」の基本的な事項は既知として書いていきます。
キーワードは 直交行列 / 三角行列 / グラム・シュミットの直交化法 です。
もし記事中でわからない言葉が出てきた場合には、線形代数学の教科書を読んだり他のサイトを見たりして調べてください。

導入

QR分解って何？

$n$ 次正方行列 $A$ が与えられたとき、それを $n$ 次直交行列 $Q$ と $n$ 次上三角行列 $R$ の積に分解することを QR分解 と言います。より一般的には $m \times n$ 行列 $A$ を $m$ 次直交行列 $Q$ と $m \times n$ 行列 $R$ に分解するのがQR分解です。
行列 $A$ がQR分解可能であるための必要十分条件は $\text{rank}(A) = n$ です。 $m = n$ 、すなわち $A$ が正方行列である場合には $A$ が正則であるときに限ってQR分解を実行することができます。

QR分解ができると何が嬉しいの？

冒頭で触れたように、QR分解は QR法 と呼ばれるアルゴリズムの基礎をなしています。QR法は固有値を求めるためのアルゴリズムで、他のアルゴリズム(ヤコビ法、べき乗法、逆反復法など)と比較して数値的に安定しているのが特徴です。

QR分解することができれば行列式や逆行列を簡単に求めることができます。直交行列の行列式は1になるので

$$
\det(A) = \det(Q R) = \det(Q) \det(R) = \det(R) = \prod_{i=1}^n R_{ii}
$$

というように行列式を求めることができます。また

$$
A^{-1} = (Q R)^{-1} = R^{-1} Q^{-1} = R^{-1} Q^T
$$

というように逆行列を求めることができます。上三角行列の逆行列は $O(n^2)$ で比較的簡単に求めることができるため、 $A^{-1}$ を同じ計算量オーダーで計算することができます。

理論

アルゴリズム

QR分解を行うための方法はいくつかありますが、ここでは グラム・シュミットの直交化法 / ハウスホルダー変換 による方法の2つを紹介したいと思います。

グラム・シュミットの直交化法

$n$本の線型独立な列ベクトル $\{ {\bf a_1}, \dots, {\bf a_n} \}$ が与えられたとき、それと同じ空間を張るような(正規)直交基底 $\{ {\bf e_1}, \dots, {\bf e_n} \}$ を抽出するアルゴリズムを グラム・シュミットの直交化法 と言います。以下の式に従って順に ${\bf e_k}$ を構成していきます。(この記事ではこの式について詳しく説明することはしません。)

\begin{eqnarray*}
{\bf \tilde{e_k}} &=& {\bf a_k} - \sum_{i=1}^{k-1} \langle {\bf a_k}, {\bf e_i} \rangle \cdot {\bf e_i} \\
{\bf e_k} &=& \frac{{\bf \tilde{e_k}}}{\| {\bf \tilde{e_k}} \|}
\end{eqnarray*}

ここで $\langle \cdot, \cdot \rangle, \| \cdot \|$ はそれぞれ標準内積、(ベクトルの)ノルムを表しています。

各$k$についてこの式を書き下してみましょう。

\begin{eqnarray*}
{\bf e_1} &=& \| {\bf \tilde{e_1}} \|^{-1} {\bf a_1} \\
{\bf e_2} &=& \| {\bf \tilde{e_2}} \|^{-1} ({\bf a_2} - \langle {\bf a_2}, {\bf e_1} \rangle {\bf e_1}) \\
{\bf e_3} &=& \| {\bf \tilde{e_3}} \|^{-1} ({\bf a_3} - \langle {\bf a_3}, {\bf e_1} \rangle {\bf e_1} - \langle {\bf a_3}, {\bf e_2} \rangle {\bf e_2}) \\
&\dots& \\
{\bf e_n} &=& \| {\bf \tilde{e_k}} \|^{-1} ({\bf a_n} - \langle {\bf a_n}, {\bf e_1} \rangle {\bf e_1} - \langle {\bf a_n}, {\bf e_2} \rangle {\bf e_2} \dots - \langle {\bf a_n}, {\bf e_{n-1}} \rangle {\bf e_{n-1}})
\end{eqnarray*}

これを各 ${\bf a_k}$ について解くと

\begin{eqnarray*}
{\bf a_1} &=& \| {\bf \tilde{e_1}} \| {\bf e_1} \\
{\bf a_2} &=& \langle {\bf a_2}, {\bf e_1} \rangle {\bf e_1} + \| {\bf \tilde{e_2}} \| {\bf e_2} \\
{\bf a_3} &=& \langle {\bf a_3}, {\bf e_1} \rangle {\bf e_1} + \langle {\bf a_3}, {\bf e_2} \rangle {\bf e_2} + \| {\bf \tilde{e_3}} \| {\bf e_3} \\
&\dots& \\
{\bf a_n} &=& \langle {\bf a_n}, {\bf e_1} \rangle {\bf e_1} + \langle {\bf a_n}, {\bf e_2} \rangle {\bf e_2} \dots + \langle {\bf a_n}, {\bf e_{n-1}} \rangle {\bf e_{n-1}} + \| {\bf \tilde{e_n}} \| {\bf e_n}
\end{eqnarray*}

となります。行列の形式で表現すれば

\begin{eqnarray*}
A &=& Q R \\
Q &\equiv& \left( {\bf e_1}, \dots, {\bf e_n} \right) \\
R &\equiv& \left(
\begin{array}{ccccc}
\| {\bf \tilde{e_1}} \| & \langle {\bf a_2}, {\bf e_1} \rangle & \dots & \dots & \langle {\bf a_n}, {\bf e_1} \rangle \\
0 & \| {\bf \tilde{e_2}} \| & \langle {\bf a_3}, {\bf e_2} \rangle & \dots & \langle {\bf a_n}, {\bf e_2} \rangle \\
\vdots & \ddots & \ddots & \ddots & \vdots \\
\vdots & & \ddots & \| {\bf \tilde{e_{n-1}}} \| & \langle {\bf a_n}, {\bf e_{n-1}} \rangle \\
0 & \dots & \dots & 0 & \| {\bf \tilde{e_n}} \|
\end{array}
\right)
\end{eqnarray*}

とすることができます。 $Q, R$ の定義から明らかなように、 $Q$ は正規直交基底を並べたものなので直交行列になっており、 $R$ は上三角行列になっています。

このように、グラム・シュミットの直交化法を利用すると、自然にQR分解を実行することができます。
なお、 $m > n$ の場合にはさらに $m - n$ 個の線形独立なベクトルをとることができますが、それを含めた $m$ 本の線形独立なベクトルを並べて $Q$ とします。 $R$ については第 $n$ 行までを上述のアルゴリズムによって計算し、第 $n+1$ 行以降はすべて0で埋めれば良いです。

ハウスホルダー変換による方法

与えられた $n$ 次元ベクトル ${\bf x}$ を $n$ 次元上のある超平面に関して鏡映変換する(対称移動させる)操作を ハウスホルダー変換 と言います。超平面はそれに直交する法線ベクトル ${\bf n}$ で特徴づけることができるので、この変換 $h$ は

$$
h({\bf x}) = {\bf x} - 2 \langle {\bf x}, {\bf n} \rangle {\bf n}
$$

と書くことができます。ハウスホルダー変換は、2次元空間においては「直線に関するベクトルの対称移動」、3次元では「平面に関するベクトルの対称移動」となっているので、そこから類推してもらえると良いかと思います。

この変換の表現行列 $H$ は

$$
H {\bf x} = (I_n - 2 {\bf n} {\bf n}^T) {\bf x}
$$

となります。( $I_n$ は $n$ 次単位行列を表しています。)この $H$ を ハウスホルダー行列 と呼びます。

ハウスホルダー行列は様々な性質を持っています。
最も重要なのは、ハウスホルダー行列は直交行列であり、かつ対称行列であるということです。これは $Q$ を構成する際になくてはならない性質となっています。この性質は容易に示すことができます。
他にも ${\bf y} = H {\bf x}$ なる ${\bf x}, {\bf y}$ に関して ${\bf x} = H {\bf y}$ が成立するという性質も重要です。これはハウスホルダー変換がベクトルの鏡映変換であるということを考えれば直感的にも明らかだと思います。

では、ハウスホルダー変換によるQR分解を見ていきましょう。
大雑把に説明すると、ハウスホルダー変換を $A$ の各列に対して順に適用していき、最終的に上三角行列にすることで $Q, R$ を求めます。(あとで詳しく説明しますが、変換過程で作用させた行列の積が $Q$ となり、最終的に得られた上三角行列が $R$ となっています。)

まず、 $A = ({\bf a_1}, \dots, {\bf a_n})$ の第1列について考えます。ハウスホルダー変換によって $A$ を上三角行列にするために、まず

\begin{eqnarray*}
H {\bf a_1} =
\left(
\begin{array}{c}
\| {\bf a_1} \| \\
0 \\
\vdots \\
0
\end{array}
\right)
\end{eqnarray*}

を満たすような変換行列 $H$ を考えます。¹ ${\bf x} \equiv {\bf a_1}, {\bf y} \equiv (\| a_1 \|, 0, \dots, 0)^T$ と書くことにすると、そのような $H$ は

$$
H = I - \frac{({\bf x} - {\bf y}) ({\bf x} - {\bf y})^T}{\| {\bf x} - {\bf y} \|^2}
$$

です。(導出は折りたたんであるので、必要ならばクリックして展開してください。)

*ハウスホルダー行列 $H$ の導出*

ハルスホルダー行列の定義式 $H = I - 2 {\bf v} {\bf v}^T$ を上の条件に代入すると

\begin{eqnarray*}
H {\bf x} = {\bf y}
&\Leftrightarrow& (I - 2 {\bf v} {\bf v}^T) {\bf x} = {\bf y} \\
&\Leftrightarrow& ({\bf v}^T {\bf x}) {\bf v} = {\bf x} - {\bf y}
\end{eqnarray*}

となり、 ${\bf v} = c ({\bf x} - {\bf y})$ と表すことができます。 ${\bf v}$ は単位ベクトルなので、定数 $c$ は一意に定まって

$$
{\bf v} = \frac{{\bf x} - {\bf y}}{\| {\bf x} - {\bf y} \|}
$$

となります。

$A$ に対してこのような $H$ を左から作用させると

\begin{eqnarray*}
H A = H ({\bf a_1}, \dots, {\bf a_n}) = 
\left(
\begin{array}{cc}
\| a_1 \| & \ast \\
{\bf 0} & B^{(1)}
\end{array}
\right)
\end{eqnarray*}

となり、第1列の第2行以降の成分を0にすることができます。ここで、 $B^{(1)}$ は $n-1$ 次正方行列です。
第2列以降についても帰納的に実行することができます。第 $k$ 列まで変換が完了していると仮定して、第 $k+1$ 列の変換を考えましょう。すなわち、 $k$ 列目まで変換された結果得られる行列を $A^{(k)}$ と書くことにすると、

\begin{eqnarray*}
A^{(k)} =
\left(
\begin{array}{cccc}
a_{11}^{(k)} & \dots & a_{1k}^{(k)} & \ast \\
& \ddots & \vdots & \ast \\
& & a_{kk}^{(k)} & \ast \\
\huge{O} & & & B^{(k)} \\
\end{array}
\right)
\end{eqnarray*}

となっていると仮定します。このとき、先ほどと同様に $n-k$ 次正方行列 $B^{(k)}$ に対してその1列目の第2行以降を0にするような変換 $G^{(k+1)}$ を施せば良いです。 $A^{(k)}$ が与えられたとき、その小行列 $B^{(k)}$ に変換 $G^{(k+1)}$ が適用されるような $n$ 次正方行列の変換 $H^{(k+1)}$ は

\begin{eqnarray*}
H^{(k+1)} =
\left(
\begin{array}{cc}
I_k & {\bf 0}^T \\
{\bf 0} & G^{(k+1)}
\end{array}
\right)
\end{eqnarray*}

となります。( $I_k$ は $k$ 次単位行列です。)この変換によって小行列 $B^{(k+1)}$ が $G^{(k+1)}$ によって変換され、その結果 $A^{(k+1)}$ の第 $k+1$ 列の第 $k+2$ 列以降の成分が0になることが確認できると思います。この変換の仕方を踏まえると

$$
A^{(n)} = H^{(n)} A^{(n-1)} = \dots = H^{(n)} \dots H^{(1)} A^{(0)} = H^{(n)} \dots H^{(1)} A
$$

が上三角行列になっています。
この上三角行列を $R \equiv A^{(n)}$ と定義します。すべての $k$ について $H^{(k)}$ は直交行列かつ対称行列であるため、

$$
A = (H^{(n)} \dots H^{(1)})^{-1} R = (H^{(1)})^{-1} \dots (H^{(n)})^{-1} R = H^{(1)} \dots H^{(n)} R
$$

と変形することができます。直交行列の積も直交行列であることから、 $Q \equiv H^{(1)} \dots H^{(n)}$ とすることで

$$
A = Q R
$$

というQR分解を得ることができます。

数学的な話

QR分解可能であるための条件

$n$ 次正方行列 $A$ がQR分解可能であることの必要十分条件は「 $A$ が正則である」ことです。これはグラム・シュミットの直交化法によって正規直交化するベクトルの組が1次独立であることからわかります。

QR分解の一意性(命題)

$n$ 次正方行列 $A$ がQR分解可能であるとき、その分解は一意的です。すなわち、 $A = Q R = Q^{\prime} R^{\prime}$ のように2通りにQR分解できたとすると、 $Q = Q^{\prime}, R = R^{\prime}$ が成立します。

QR分解の一意性(証明)

$n$ 次正方行列 $A$ がQR分解可能であるとき、その分解は一意的である。

この命題は $Q, R$ がともに正則であることに注意すれば容易に示すことができます。²
$A = Q R = Q^{\prime} R^{\prime}$ のように2通りにQR分解できたと仮定します。このとき、 $Q, R$ はともに正則なので

$$
(Q^{\prime})^{-1} Q = R^{\prime} R^{-1}
$$

となります。上三角行列どうしの積、および(正則な)上三角行列の逆行列が上三角行列になることから、右辺は上三角行列になります。一方、左辺は直交行列になっています。

等号から左辺および右辺は直交行列かつ上三角行列になります。このような行列は単位行列に限られることが知られています。(証明は折りたたんであります。)

*直交行列かつ上三角行列が単位行列になることの証明*

> $X$ を直交行列かつ上三角行列であるとします。

$X$ の逆行列を考えると、 $X$ は直交行列であることから $X^{-1} = X^T$ となっています。一方、 $X$ は上三角行列であることから $X^{-1}$ は上三角行列になっています。 $X$ および $X^T$ が上三角行列であるため、 $X$ (および $X^{-1}$) は対角行列です。
直交行列の各列のノルムが1であることを考慮すると、その対角成分は1となる必要があります。

以上より、直交行列かつ上三角行列は単位行列となります。

したがって

$$
(Q^{\prime})^{-1} Q = R^{\prime} R^{-1} = I_n \Leftrightarrow Q = Q^{\prime}, R = R^{\prime}
$$

が成立します。

実装

QR分解(グラム・シュミットの直交化法 / ハウスホルダー変換による方法)のPythonによる実装はGitHubにあります。

グラム・シュミットの直交化法において、行数 $m$ が列数 $n$ よりも大きい場合には、 $m$ 次元空間の標準基底を順番に見ていき、それまでに追加された基底と線型独立なものを探すようにしています。 while ループで実装していますが、このループは必ず停止することに注意してください。³

後書き

ラーニングピラミッドによれば「他者に伝える」ことで知識の定着率が上がるそうなので、これからも積極的に記事を書いていきたいと思います。

参考文献・リンク

直交変換(直交行列を表現行列とするような変換)はベクトルのノルムを保つため、変換後のベクトルの第2成分以降を0とすると、第1成分は一意に定まります。 ↩
グラム・シュミットの直交化法によるQR分解における $R$ の対角成分を見てみると、すべて $A$ の行ベクトルのノルムになっています。 $A$ は正則なので、それらが0になることはありません。 ↩
$m$ 個未満のベクトルの組で $m$ 次元空間の標準基底のすべてと線形従属になるようなものは存在しません。言い換えれば、 $m$ 個未満のベクトルで $m$ 次元空間を張ることはできません。そのため、標準基底のいずれかは線形独立なベクトルとして採用することができます。 ↩