QR分解
QR分解(キューアールぶんかい、QR decomposition, QR factorization)とは、m × n 実行列 Aを、 m 次直交行列 Q と m × n 上三角行列 R との積への分解により表すこと、またはそう表した表現をいう。このような分解は常に存在する。
QR分解は線型最小二乗問題を解くために使用される。また、固有値問題の数値解法の1つであるQR法の基礎となっている。
すべての実正方行列 Aは直交行列Qと上三角行列(別名右三角行列)Rを用いて
A = QR
と分解できる。もしAが正則ならば、Rの対角成分が正になるような因数分解は一意に定まる。
もしAが複素正方行列ならば、Qがユニタリ行列 (つまりQ^* Q = QQ^* = I)となるような分解A = QRが存在する。
もしAがn個の線形独立な列を持つなら、Qの最初のn列はAの列空間の正規直交基底をなす。より一般的に、1 ≤ k ≤ nの任意のkについて、Qの最初のk列はAの最初のk列の線型包をなす。Aの任意の列kがQの最初のk列にのみ依存するということは、Rが三角行列であることから明らかである。
より一般的に、m ≥ nである複素m×n行列Aを、m×mユニタリ行列Qとm×n上三角行列Rに分解することができる。m×n上三角行列の下から(m−n)行はすべてゼロであるため、Rや、RとQ両方の分割を簡単に行うことができる。
A = QR = Q \begin{bmatrix} R_1 \\ 0 \end{bmatrix}
= \begin{bmatrix} Q_1, Q_2 \end{bmatrix} \begin{bmatrix} R_1 \\ 0 \end{bmatrix}
= Q_1 R_1
ここで、R1はn×n上三角行列、0は(m − n)×n零行列、Q1はm×n行列、Q2はm×(m − n)行列で、Q1とQ2は両方直交する列を持つ。
Q1R1をはAの薄い(thin)QR分解と呼び、 Trefethen & Bauは軽減(reduced)QR分解と呼んでいる。もしAが最大階数nであり、R1の対角成分を正にするならば、R1とQ1は一意に定まる。しかし一般的にQ2はそうではない。R1はA* A (Aが実行列の場合ATAに等しい)のコレスキー分解の上三角部分に等しい。
同様に、Lを下(lower)三角行列として、QL、RQ、LQ分解を定義することができる。
QR分解を計算する手法として、グラム・シュミット分解、ハウスホルダー変換、ギブンス回転などがある。それぞれ利点と欠点がある。
グラム・シュミットの正規直交化法を最大階数行列の列A = \left[\boldsymbol{a}_1, \ldots, \boldsymbol{a}_n\right]に適用することを考える。内積\langle\boldsymbol{v}, \boldsymbol{w}\rangle = \boldsymbol{v}^\textsf{T} \boldsymbol{w} (複素ベクトルの場合 \langle\boldsymbol{v}, \boldsymbol{w}\rangle = \boldsymbol{v}^* \boldsymbol{w})とする。
射影の定義より、
\operatorname{proj}_{\boldsymbol{u}\boldsymbol{a} =
\frac{\left\langle\boldsymbol{u}, \boldsymbol{a}\right\rangle}{\left\langle\boldsymbol{u}, \boldsymbol{u}\right\rangle}{\boldsymbol{u}
したがって、
\begin{align}
\boldsymbol{u}_1 &= \boldsymbol{a}_1, &
\boldsymbol{e}_1 &= {\boldsymbol{u}_1 \over \|\boldsymbol{u}_1\〕 \\
\boldsymbol{u}_2 &= \boldsymbol{a}_2 - \operatorname{proj}_{\boldsymbol{u}_1}\,\boldsymbol{a}_2, &
\boldsymbol{e}_2 &= {\boldsymbol{u}_2 \over \|\boldsymbol{u}_2\〕 \\
\boldsymbol{u}_3 &= \boldsymbol{a}_3 - \operatorname{proj}_{\boldsymbol{u}_1}\,\boldsymbol{a}_3 - \operatorname{proj}_{\boldsymbol{u}_2}\,\boldsymbol{a}_3, &
\boldsymbol{e}_3 &= {\boldsymbol{u}_3 \over \|\boldsymbol{u}_3\〕 \\
& \vdots &
& \vdots \\
\boldsymbol{u}_k &= \boldsymbol{a}_k - \sum_{j=1}^{k-1}\operatorname{proj}_{\boldsymbol{u}_j}\,\boldsymbol{a}_k,&
\boldsymbol{e}_k &= {\boldsymbol{u}_k \over \|\boldsymbol{u}_k\〕
\end{align}
ここで\boldsymbol{a}_iを新しく計算された正規直交基底上に表すことができ、\left\langle\boldsymbol{e}_i, \boldsymbol{a}_i\right\rangle = \left\|\boldsymbol{u}_i\right\|であるから、
\begin{align}
\boldsymbol{a}_1 &= \langle\boldsymbol{e}_1, \boldsymbol{a}_1\rangle \boldsymbol{e}_1 \\
\boldsymbol{a}_2 &= \langle\boldsymbol{e}_1, \boldsymbol{a}_2\rangle \boldsymbol{e}_1
+ \langle\boldsymbol{e}_2, \boldsymbol{a}_2\rangle \boldsymbol{e}_2 \\
\boldsymbol{a}_3 &= \langle\boldsymbol{e}_1, \boldsymbol{a}_3\rangle \boldsymbol{e}_1
+ \langle\boldsymbol{e}_2, \boldsymbol{a}_3\rangle \boldsymbol{e}_2
+ \langle\boldsymbol{e}_3, \boldsymbol{a}_3\rangle \boldsymbol{e}_3 \\
&\vdots \\
\boldsymbol{a}_k &= \sum_{j=1}^k \langle \boldsymbol{e}_j, \boldsymbol{a}_k \rangle \boldsymbol{e}_j
\end{align}
これは行列の形に書くことができ、
A = QR
ただし、
Q = \left[\boldsymbol{e}_1, \ldots, \boldsymbol{e}_n\right],
R = \begin{bmatrix}
\langle\boldsymbol{e}_1, \boldsymbol{a}_1\rangle &
\langle\boldsymbol{e}_1, \boldsymbol{a}_2\rangle &
\langle\boldsymbol{e}_1, \boldsymbol{a}_3\rangle & \ldots \\
0 &
\langle\boldsymbol{e}_2, \boldsymbol{a}_2\rangle &
\langle\boldsymbol{e}_2, \boldsymbol{a}_3\rangle & \ldots \\
0 &
0 &
\langle\boldsymbol{e}_3, \boldsymbol{a}_3\rangle & \ldots \\
\vdots &
\vdots &
\vdots &
\ddots
\end{bmatrix}
A = \begin{bmatrix}
12 & -51 &   4 \\
6 & 167 & -68 \\
-4 &  24 & -41
\end{bmatrix}
の分解を考える。
正規直交行列Qに対して、
Q^\textsf{T} \,Q = I
が常に成り立つから、グラム・シュミット法により以下の手順でQを計算できる。
\begin{align}
U = \begin{bmatrix}
\boldsymbol u_1 & \boldsymbol u_2 & \boldsymbol u_3
\end{bmatrix} &=
\begin{bmatrix}
12 & -69 & -58/5 \\
6 & 158 &   6/5 \\
-4 &  30 & -33
\end{bmatrix}; \\
Q = \begin{bmatrix}
\dfrac{\boldsymbol u_1}{\|\boldsymbol u_1\〕 &
\dfrac{\boldsymbol u_2}{\|\boldsymbol u_2\〕 &
\dfrac{\boldsymbol u_3}{\|\boldsymbol u_3\〕
\end{bmatrix} &=
\begin{bmatrix}
6/7 & -69/175 & -58/175 \\
3/7 & 158/175 &   6/175 \\
-2/7 &   6/35  & -33/35
\end{bmatrix}
\end{align}
したがって、
\begin{align}
Q^\textsf{T} A &= Q^\textsf{T}Q\,R = R; \\
R &= Q^\textsf{T}A =
\begin{bmatrix}
14 &  21 & -14 \\
0 & 175 & -70 \\
0 &   0 &  35
\end{bmatrix}
\end{align}
RQ分解は行列Aを上三角行列R(別名右三角)と直交行列Qに変換する。QR分解との違いはこれらの行列の順番だけである。
QR分解はグラム・シュミットの正規直交化法をAの最初の列から最後の列の順に適用する。
RQ分解はグラム・シュミットの正規直交化法をAの最後の行から最初の行の順に適用する。
グラム・シュミットの正規直交化法は本来、数値的に不安定である。射影の応用として直交化との幾何学的な類似性があるが、直交化自体は数値的な差異が生じやすい。しかしながら、実装が簡単という大きな利点があり、外部線形代数ライブラリが利用できない場合のプロトタイピングに便利なアルゴリズムである。
ハウスホルダー変換 (またはハウスホルダー鏡映)はベクトルを取り、平面または超平面に関する鏡映をする変換である。この演算はm×n行列A (m ≥ n)のQR変換の計算に使うことができる。
Qは一つの座標を除いたすべての座標が未知でもベクトルを鏡映するために使用できる。
スカラαに対して\|\boldsymbol{x}\| = |\alpha|を満たすようなAの任意の実m次元列ベクトル\boldsymbol{x}を考える。もしこのアルゴリズムが浮動小数点演算を用いて実装されている場合、桁落ちを防ぐため、行列Aの最終的な上三角部分においてすべての要素が0である後のピボット座標x_kに対し、αを\boldsymbol{x}のk番目の座標の逆符号とする。複素行列の場合、
\alpha = -\exp(i \arg x_k) \|\boldsymbol{x}\|
とし、以下のQの導出において転置を共役転置に読み替えること。
ここで、\boldsymbol{e}_1をベクトル(1 0 … 0)T、||·||をユークリッド距離、Iをm×m単位行列とし、
\begin{align}
\boldsymbol{u} &= \boldsymbol{x} - \alpha\boldsymbol{e}_1, \\
\boldsymbol{v} &= {\boldsymbol{u} \over \|\boldsymbol{u}\〕, \\
Q &= I - 2 \boldsymbol{v}\boldsymbol{v}^\textsf{T}
\end{align}
あるいは、Aが複素行列ならば、
Q = I - 2\boldsymbol{v}\boldsymbol{v}^*
とおく。
Qはm×mハウスホルダー行列であり、
Q\boldsymbol{x} = \begin{bmatrix} \alpha & 0 & \cdots & 0 \end{bmatrix}^\textsf{T} \
これによりm×n行列Aを上三角の形に漸次変換できる。まず、xの最初の列を選んで得られるハウスホルダー行列Q1にAを乗算する。この結果、行列Q1Aは左の列が(最初の行を除き)ゼロになる。
Q_1A = \begin{bmatrix}
\alpha_1 & \star & \dots & \star \\
0 &       &       &       \\
\vdots &       &    A' &       \\
0 &       &       &
\end{bmatrix}
この操作をA′ (Q1Aから最初の列と最初の行を除いたもの)に繰り返すと、ハウスホルダー行列Q′2が得られる。Q′2はQ1より小さいということに注意すること。A′の代わりにQ1Aで計算したいため、A′を左上に拡張し、ひとつの1を埋める必要がある。つまり、一般的には
Q_k = \begin{bmatrix}
I_{k-1} & 0    \\
0  & Q_k'
\end{bmatrix}
とする。
t回このプロセスを繰り返すと、t = \min(m - 1, n)のとき、
R = Q_t \cdots Q_2 Q_1 A
は上三角行列である。であるからして、
Q = Q_1^\textsf{T} Q_2^\textsf{T} \cdots Q_t^\textsf{T}
とすると、
A = QRはAのQR分解である。
このメソッドは先述のグラム・シュミット法よりも数値的安定性がある。
下表にサイズnの正方行列を仮定したときの、ハウスホルダー変換によるQR分解のk番目のステップにおける計算量を示す。
これらの数をn − 1ステップ(サイズnの正方行列であるため)まで合計して、このアルゴリズムの(浮動小数点演算の観点からの)複雑さは
\frac{2}{3}n^3 + n^2 + \frac{1}{3}n - 2 = O\left(n^3\right)
と表せる。
A = \begin{bmatrix}
12 & -51 &   4 \\
6 & 167 & -68 \\
-4 &  24 & -41
\end{bmatrix}
の分解を考える。
まず、行列Aの最初の列、ベクトル\boldsymbol{a}_1 = \begin{bmatrix} 12 & 6 & -4 \end{bmatrix}^\textsf{T}を、\left\|\boldsymbol{a}_1\right\| \;\boldsymbol{e}_1 = \begin{bmatrix} 1 & 0 & 0\end{bmatrix}^\textsf{T} に変換する鏡映を見つける必要がある。
今、
\begin{align}
\boldsymbol{u} &= \boldsymbol{x} - \alpha\boldsymbol{e}_1, \\
\boldsymbol{v} &= \frac{\boldsymbol{u}{\|\boldsymbol{u}\〕
\end{align}
とする。
ここで、
\alpha = 14であり、\boldsymbol{x} = \boldsymbol{a}_1 = \begin{bmatrix} 12 & 6 & -4 \end{bmatrix}^\textsf{T}
であるから、
\boldsymbol{u} = \begin{bmatrix} -2 & 6 & -4 \end{bmatrix}^\textsf{T} = \begin{bmatrix} 2 \end{bmatrix}\begin{bmatrix} -1 & 3 & -2 \end{bmatrix}^\textsf{T} and \boldsymbol{v} = {1 \over \sqrt{14}\begin{bmatrix} -1 & 3 & -2 \end{bmatrix}^\textsf{T} であり、
\begin{align}
={} &I - {2 \over \sqrt{14}\sqrt{14}
\begin{bmatrix} -1 \\ 3 \\ -2 \end{bmatrix}
\begin{bmatrix} -1 &  3 &  -2 \end{bmatrix} \\
={} &I - {1 \over 7}\begin{bmatrix}
1 & -3 &  2 \\
-3 &  9 & -6 \\
2 & -6 &  4
\end{bmatrix} \\
={} &\begin{bmatrix}
6/7 &  3/7 & -2/7 \\
3/7 & -2/7 &  6/7 \\
-2/7 &  6/7 &  3/7 \\
\end{bmatrix}
\end{align}
である。
今、
Q_1A = \begin{bmatrix}
14 &  21 & -14 \\
0 & -49 & -14 \\
0 & 168 & -77
\end{bmatrix}
を見てみると、すでにほぼ三角行列である。あとは(3, 2)要素をゼロにするだけでよい。
(1, 1)における小行列を取り、同じプロセスを
A' = M_{11} = \begin{bmatrix}
-49 & -14 \\
168 & -77
\end{bmatrix}
に再び適用する。
先述のメソッドと同様にして、このプロセスの次のステップが正しく動作するために、1で直和を取ることにより、ハウスホルダー変換
Q_2 = \begin{bmatrix}
1 &     0 &  0 \\
0 & -7/25 & 24/25 \\
0 & 24/25 &  7/25
\end{bmatrix}
を得る。
今、
Q = Q_1^\textsf{T} Q_2^\textsf{T} = \begin{bmatrix}
6/7 & -69/175 & 58/175 \\
3/7 & 158/175 & -6/175 \\
-2/7 &   6/35  & 33/35
\end{bmatrix}
または、有効数字四桁で、
\begin{align}
Q &= Q_1^\textsf{T} Q_2^\textsf{T} = \begin{bmatrix}
0.8571 & -0.3943 &  0.3314 \\
0.4286 &  0.9029 & -0.0343 \\
-0.2857 &  0.1714 &  0.9429
\end{bmatrix} \\
R &= Q_2 Q_1 A = Q^\textsf{T} A = \begin{bmatrix}
14 &  21 & -14 \\
0 & 175 & -70 \\
0 &   0 & -35
\end{bmatrix}
\end{align}
を得た。
行列Qは直交行列であり、Rは上三角行列であるため、A = QRは求めるべきQR分解である。
ハウスホルダー変換の使用は、R行列のゼロを生成するメカニズムに鏡映を利用しているため、最もシンプルな数値的に安定したQR分解アルゴリズムである。しかしながら、新しいゼロ要素を生成するすべての鏡映においてQ、R行列両方の全体が変化するため、ハウスホルダー変換は帯域幅が重く、並列化できない。
QR分解はギブンス回転を使用しても計算できる。各回転により行列の亜対角要素がゼロになり、R行列を構成できる。すべてのギブンス回転を結合することで直交行列Qを構成できる。
実際には、行列全体を構成して乗算をするようなギブンス回転は行われない。代わりに、疎な要素を計算するような無駄な計算をしない、疎なギブンス行列乗算と同等なあるギブンス回転の手順が採られる。そのギブンス回転の手順は少しの非対角成分をゼロにするだけで済み、ハウスホルダー変換よりも容易に並列化できる。
A = \begin{bmatrix}
12 & -51 &   4 \\
6 & 167 & -68 \\
-4 &  24 & -41
\end{bmatrix}
の分解を考える。
まず、左下隅の要素、\boldsymbol{a}_{31} = -4をゼロにする回転行列を構成する必要がある。この行列G_1はギブンス回転で求めることができる。まずベクトル\begin{bmatrix} 12 & -4 \end{bmatrix}を、X軸に沿って回転させる。このベクトルは角度\theta = \arctan\left[ {-(-4) \over 12}\right] を持つ。直交ギブンス回転行列G_1を次のように作る。
\begin{align}
G_1 &= \begin{bmatrix}
\cos(\theta) & 0 & -\sin(\theta) \\
0 & 1 &             0 \\
\sin(\theta) & 0 &  \cos(\theta)
\end{bmatrix} \\
&\approx \begin{bmatrix}
0.94868 & 0 & -0.31622 \\
0       & 1 &  0       \\
0.31622 & 0 &  0.94868
\end{bmatrix}
\end{align}
ここでG_1Aの結果は\boldsymbol{a}_{31}要素がゼロである。
G_1A \approx \begin{bmatrix}
12.64911 & -55.97231 &  16.76007 \\
6       & 167       & -68       \\
0       &   6.64078 & -37.6311
\end{bmatrix}
同様にしてそれぞれ非対角要素a_{21}・a_{32}要素がゼロであるようなギブンス行列G_2・G_3を構成し、三角行列Rを構成する。直交行列Q^\textsf{T}はすべてのギブンス行列の積Q^\textsf{T} = G_3 G_2 G_1で表される。したがって、G_3 G_2 G_1 A = Q^\textsf{T} A = Rであり、QR分解はA = QRである。
ギブンス回転によるQR分解は、アルゴリズムを完全に動かすのに必要な行の順序を決定するのが簡単ではないため、実装に最も手間がかかる。しかしながら、新しいゼロ要素a_{ij}がゼロになる予定の要素の行(i)とその上の行(j)にしか影響しないという特筆すべき利点がある。これにより、ギブンス回転アルゴリズムはハウスホルダー変換手法よりも帯域幅効率が良く、容易に並列化できる。
QR分解を正方行列の行列式の絶対値を求めるのに利用できる。ある行列がA = QRと分解できるとする。このとき
\det(A) = \det(Q) \cdot \det(R)
である。
Qはユニタリであるため、|\det(Q)| = 1である。したがって、r_{ii}をRの対角要素とすると、
\left|\det(A)\right| = \left|\det(R)\right| = \left|\prod_{i} r_{ii}\right|
となる。
さらに、行列式は固有値の積に等しいため、\lambda_iをAの固有値とすると、
\left|\prod_{i} r_{ii}\right| = \left|\prod_{i} \lambda_{i}\right|
となる。
QR分解の定義を非正方行列に導入し、固有値を特異値に置き換えることで、上記性質を非正方行列Aに拡張することができる。
非正方行列AのQR分解を
A = Q \begin{bmatrix} R \\ O \end{bmatrix}, \qquad Q^* Q = I
とする。ただし、Oは零行列、Qはユニタリ行列。
特異値分解と行列式の性質から、\sigma_iをAの特異値として、
\left|\prod_i r_{ii}\right| = \prod_i\sigma_{i}
AとRの特異値は同じであるが、複素固有値が異なる場合があることに注意すること。しかしながら、Aが正方ならば、下記は真である。
{\prod_i \sigma_i} = \left|{\prod_i \lambda_i}\right|
結論として、QR分解を使うことによって行列の固有値や特異値の積を効率よく計算することができる。
ピボットQRは列のピボットの新しいステップにおいて、それぞれ初めに残りの列で最も大きいものを取るという点で、通常のグラム・シュミット法とは異なっている。したがって、置換行列Pを次のように導入する。
AP = QR\quad \iff\quad A = QRP^\textsf{T}
列のピボットはAが(ほぼ)階数落ちである、またはその疑いがある場合に便利である。また、数値的精度を向上させることもできる。通常、Rの対角成分が非増加、つまり\left|r_{11}\right| \ge \left|r_{22}\right| \ge \ldots \ge \left|r_{nn}\right|となるようにPを選ぶ。この手段により特異値分解よりも低い計算コストでAの(数値的な)階数を求めることができ、いわゆるRank Revealing QR分解の基礎となっている。
行列の逆行列を直接求めるのに比べ、QR分解を利用した逆問題の解法は、条件数が減少していることからも分かるように、数値的に安定している。
次元がm \times nで階数がmであるような行列Aに対して、劣決定(m < n)線形問題Ax = bを解くためには、まずAの転置行列のQR分解A^\textsf{T} = QRを求める。ただし、Qは直交行列(つまり、Q^\textsf{T} = Q^{-1})であり、RはR = \begin{bmatrix} R_1 \\ 0 \end{bmatrix}という特殊な形である。ここで、R_1はm \times m正方右三角行列、零行列は(n-m) \times m次元である。計算すると、この逆問題の解を次のように表すことができる。
x = Q \begin{bmatrix}
\left(R_1^\textsf{T}\right)^{-1}b \\
\end{bmatrix}
ここで、R_1^{-1}はガウスの消去法で計算でき、\left(R_1^\textsf{T}\right)^{-1} bは前方置換法を用いることで直接計算できる。後者の手法の方が数値的精度が高く、計算量も少ないという利点がある。
ノルム\|A \hat{x} - b\|を最小にするような過決定(m \geq n)問題Ax = bの解\hat{x}を求めるためには、まずAのQR分解A = QRを求める。Q_1を直交行列Q全体のうち最初のn列を含むm \times n行列、R_1を先述の通りに置くと、この問題の解は\hat{x} = R_1^{-1} \left(Q_1^\textsf{T} b\right) と表せる。劣決定の場合と同様に、R_1の逆行列を直接計算しなくても、後方置換法を用いることで早く正確に\hat xを求めることができる。(Q_1とR_1は数値ライブラリによっては高速な(economic)QR分解として実装されている。)
岩澤分解はQR分解を半単純リー群に一般化している。