QR分解

QR分解(キューアールぶんかい、: QR decomposition, QR factorization)とは、m × n 実行列 Aを、 m直交行列 Qm × n 上三角行列 R との積への分解により表すこと、またはそう表した表現をいう[1]。このような分解は常に存在する[2]

QR分解は線型最小二乗問題を解くために使用される。また、固有値問題の数値解法の1つであるQR法の基礎となっている。

定義

正方行列

すべての実正方行列 A直交行列Qと上三角行列(別名右三角行列)Rを用いて

A = Q R {\displaystyle A=QR}

と分解できる。もしA正則ならば、Rの対角成分が正になるような因数分解は一意に定まる。

もしAが複素正方行列ならば、Qユニタリ行列 (つまり Q Q = Q Q = I {\displaystyle Q^{*}Q=QQ^{*}=I} )となるような分解A = QRが存在する。

もしAn個の線形独立な列を持つなら、Qの最初のn列はA列空間正規直交基底をなす。より一般的に、1 ≤ k ≤ nの任意のkについて、Qの最初のk列はAの最初のk列の線型包をなす[3]。Aの任意の列kQの最初のk列にのみ依存するということは、Rが三角行列であることから明らかである[3]

矩形行列

より一般的に、m ≥ nである複素m×n行列Aを、m×mユニタリ行列Qm×n上三角行列Rに分解することができる。m×n上三角行列の下から(mn)行はすべてゼロであるため、Rや、RQ両方の分割を簡単に行うことができる。

A = Q R = Q [ R 1 0 ] = [ Q 1 , Q 2 ] [ R 1 0 ] = Q 1 R 1 {\displaystyle 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}}

ここで、R1n×n上三角行列、0は(m − nn零行列、Q1m×n行列、Q2m×(m − n)行列で、Q1Q2は両方直交する列を持つ。

Q1R1Golub & Van Loan (1996, §5.2)はAの薄い(thin)QR分解と呼び、 Trefethen & Bauは軽減(reduced)QR分解と呼んでいる[3]。もしAが最大階数nであり、R1の対角成分を正にするならば、R1Q1は一意に定まる。しかし一般的にQ2はそうではない。R1A* A (Aが実行列の場合ATAに等しい)のコレスキー分解の上三角部分に等しい。

QL・RQ・LQ分解

同様に、Lを下(lower)三角行列として、QL、RQ、LQ分解を定義することができる。

QR分解の計算

QR分解を計算する手法として、グラム・シュミット分解ハウスホルダー変換ギブンス回転などがある。それぞれ利点と欠点がある。

グラム・シュミットの正規直交化法の使用

グラム・シュミットの正規直交化法を最大階数行列の列 A = [ a 1 , , a n ] {\displaystyle A=\left[{\boldsymbol {a}}_{1},\ldots ,{\boldsymbol {a}}_{n}\right]} に適用することを考える。内積 v , w = v T w {\displaystyle \langle {\boldsymbol {v}},{\boldsymbol {w}}\rangle ={\boldsymbol {v}}^{\textsf {T}}{\boldsymbol {w}}} (複素ベクトルの場合 v , w = v w {\displaystyle \langle {\boldsymbol {v}},{\boldsymbol {w}}\rangle ={\boldsymbol {v}}^{*}{\boldsymbol {w}}} )とする。

射影の定義より、

proj u a = u , a u , u u {\displaystyle \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}}}

したがって、

u 1 = a 1 , e 1 = u 1 u 1 u 2 = a 2 proj u 1 a 2 , e 2 = u 2 u 2 u 3 = a 3 proj u 1 a 3 proj u 2 a 3 , e 3 = u 3 u 3 u k = a k j = 1 k 1 proj u j a k , e k = u k u k {\displaystyle {\begin{aligned}{\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{aligned}}}

ここで a i {\displaystyle {\boldsymbol {a}}_{i}} を新しく計算された正規直交基底上に表すことができ、 e i , a i = u i {\displaystyle \left\langle {\boldsymbol {e}}_{i},{\boldsymbol {a}}_{i}\right\rangle =\left\|{\boldsymbol {u}}_{i}\right\|} であるから、

a 1 = e 1 , a 1 e 1 a 2 = e 1 , a 2 e 1 + e 2 , a 2 e 2 a 3 = e 1 , a 3 e 1 + e 2 , a 3 e 2 + e 3 , a 3 e 3 a k = j = 1 k e j , a k e j {\displaystyle {\begin{aligned}{\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{aligned}}}

これは行列の形に書くことができ、

A = Q R {\displaystyle A=QR}

ただし、

Q = [ e 1 , , e n ] , {\displaystyle Q=\left[{\boldsymbol {e}}_{1},\ldots ,{\boldsymbol {e}}_{n}\right],}
R = [ e 1 , a 1 e 1 , a 2 e 1 , a 3 0 e 2 , a 2 e 2 , a 3 0 0 e 3 , a 3 ] {\displaystyle 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 = [ 12 51 4 6 167 68 4 24 41 ] {\displaystyle A={\begin{bmatrix}12&-51&4\\6&167&-68\\-4&24&-41\end{bmatrix}}}

の分解を考える。

正規直交行列 Q {\displaystyle Q} に対して、

Q T Q = I {\displaystyle Q^{\textsf {T}}\,Q=I}

が常に成り立つから、グラム・シュミット法により以下の手順で Q {\displaystyle Q} を計算できる。

U = [ u 1 u 2 u 3 ] = [ 12 69 58 / 5 6 158 6 / 5 4 30 33 ] ; Q = [ u 1 u 1 u 2 u 2 u 3 u 3 ] = [ 6 / 7 69 / 175 58 / 175 3 / 7 158 / 175 6 / 175 2 / 7 6 / 35 33 / 35 ] {\displaystyle {\begin{aligned}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{aligned}}}

したがって、

Q T A = Q T Q R = R ; R = Q T A = [ 14 21 14 0 175 70 0 0 35 ] {\displaystyle {\begin{aligned}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{aligned}}}

RQ分解との関係

RQ分解は行列Aを上三角行列R(別名右三角)と直交行列Qに変換する。QR分解との違いはこれらの行列の順番だけである。

QR分解はグラム・シュミットの正規直交化法をA最初の列から最後の列の順に適用する。

RQ分解はグラム・シュミットの正規直交化法をA最後の行から最初の行の順に適用する。

利点と欠点

グラム・シュミットの正規直交化法は本来、数値的に不安定である。射影の応用として直交化との幾何学的な類似性があるが、直交化自体は数値的な差異が生じやすい。しかしながら、実装が簡単という大きな利点があり、外部線形代数ライブラリが利用できない場合のプロトタイピングに便利なアルゴリズムである。

ハウスホルダー変換の使用

QR分解のためのハウスホルダー変換: 目標はベクトル x {\displaystyle x} を同じ長さかつ e 1 {\displaystyle e_{1}} の共線であるベクトルに変換する線形変換を見つけることである。直交射影(グラム・シュミット法)を使うこともできるが、ベクトル x {\displaystyle x} e 1 {\displaystyle e_{1}} が直交に近い場合、数値的に不安定である。代わりに、ハウスホルダー変換により点線を通して鏡映する(点線は x {\displaystyle x} e 1 {\displaystyle e_{1}} のなす角の二等分線)。この変換による最大角は45度である。

ハウスホルダー変換 (またはハウスホルダー鏡映)はベクトルを取り、平面または超平面に関する鏡映をする変換である。この演算はm×n行列 A {\displaystyle A} (m ≥ n)のQR変換の計算に使うことができる。

Qは一つの座標を除いたすべての座標が未知でもベクトルを鏡映するために使用できる。

スカラαに対して x = | α | {\displaystyle \|{\boldsymbol {x}}\|=|\alpha |} を満たすような A {\displaystyle A} の任意の実m次元列ベクトル x {\displaystyle {\boldsymbol {x}}} を考える。もしこのアルゴリズムが浮動小数点演算を用いて実装されている場合、桁落ち(による誤差の拡大)を防ぐため、行列Aの最終的な上三角部分においてすべての要素が0である後のピボット座標 x k {\displaystyle x_{k}} に対し、α x {\displaystyle {\boldsymbol {x}}} k番目の座標の逆符号とする。複素行列の場合には、

α = exp ( i arg x k ) x {\displaystyle \alpha =-\exp(i\arg x_{k})\|{\boldsymbol {x}}\|}

(Stoer & Bulirsch 2002, p. 225)として、さらに以下のQの導出において転置を共役転置に読み替えること。

ここで、 e 1 {\displaystyle {\boldsymbol {e}}_{1}} をベクトル(1 0 … 0)T、||·||をユークリッド距離 I {\displaystyle I} m×m単位行列とし、

u = x α e 1 , v = u u , Q = I 2 v v T {\displaystyle {\begin{aligned}{\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{aligned}}}

あるいは、 A {\displaystyle A} が複素行列ならば、

Q = I 2 v v {\displaystyle Q=I-2{\boldsymbol {v}}{\boldsymbol {v}}^{*}}

とおく。

Q {\displaystyle Q} m×mハウスホルダー行列であり、

Q x = [ α 0 0 ] T   {\displaystyle Q{\boldsymbol {x}}={\begin{bmatrix}\alpha &0&\cdots &0\end{bmatrix}}^{\textsf {T}}\ }

これによりm×n行列Aを上三角の形に漸次変換できる。まず、xの最初の列を選んで得られるハウスホルダー行列Q1Aを乗算する。この結果、行列Q1Aは左の列が(最初の行を除き)ゼロになる。

Q 1 A = [ α 1 0 A 0 ] {\displaystyle Q_{1}A={\begin{bmatrix}\alpha _{1}&\star &\dots &\star \\0&&&\\\vdots &&A'&\\0&&&\end{bmatrix}}}

この操作をA′ (Q1Aから最初の列と最初の行を除いたもの)に繰り返すと、ハウスホルダー行列Q2が得られる。Q2Q1より小さいということに注意すること。A′の代わりにQ1Aで計算したいため、A′を左上に拡張し、ひとつの1を埋める必要がある。つまり、一般的には

Q k = [ I k 1 0 0 Q k ] {\displaystyle Q_{k}={\begin{bmatrix}I_{k-1}&0\\0&Q_{k}'\end{bmatrix}}}

とする。 t {\displaystyle t} 回このプロセスを繰り返すと、 t = min ( m 1 , n ) {\displaystyle t=\min(m-1,n)} のとき、

R = Q t Q 2 Q 1 A {\displaystyle R=Q_{t}\cdots Q_{2}Q_{1}A}

は上三角行列である。そこで、

Q = Q 1 T Q 2 T Q t T {\displaystyle Q=Q_{1}^{\textsf {T}}Q_{2}^{\textsf {T}}\cdots Q_{t}^{\textsf {T}}}

とすると、 A = Q R {\displaystyle A=QR} A {\displaystyle A} のQR分解である。

この鏡映変換を用いた計算方法は先述のグラム・シュミット法よりも数値的安定性がある。

下表にサイズnの正方行列を仮定したときの、ハウスホルダー変換によるQR分解のk番目のステップにおける計算量を示す。

演算 k番目のステップにおける計算量
乗算 2 ( n k + 1 ) 2 {\displaystyle 2(n-k+1)^{2}}
加算 ( n k + 1 ) 2 + ( n k + 1 ) ( n k ) + 2 {\displaystyle (n-k+1)^{2}+(n-k+1)(n-k)+2}
除算 1 {\displaystyle 1}
平方根 1 {\displaystyle 1}

これらの数をn − 1ステップ(サイズnの正方行列であるため)まで合計して、このアルゴリズムの(浮動小数点演算の観点からの)複雑さは

2 3 n 3 + n 2 + 1 3 n 2 = O ( n 3 ) {\displaystyle {\frac {2}{3}}n^{3}+n^{2}+{\frac {1}{3}}n-2=O\left(n^{3}\right)}

と表せる。

A = [ 12 51 4 6 167 68 4 24 41 ] {\displaystyle A={\begin{bmatrix}12&-51&4\\6&167&-68\\-4&24&-41\end{bmatrix}}}

の分解を考える。

まず、行列Aの最初の列、ベクトル a 1 = [ 12 6 4 ] T {\displaystyle {\boldsymbol {a}}_{1}={\begin{bmatrix}12&6&-4\end{bmatrix}}^{\textsf {T}}} を、 a 1 e 1 = [ 1 0 0 ] T {\displaystyle \left\|{\boldsymbol {a}}_{1}\right\|\;{\boldsymbol {e}}_{1}={\begin{bmatrix}1&0&0\end{bmatrix}}^{\textsf {T}}} に変換する鏡映を見つける必要がある。

今、

u = x α e 1 , v = u u {\displaystyle {\begin{aligned}{\boldsymbol {u}}&={\boldsymbol {x}}-\alpha {\boldsymbol {e}}_{1},\\{\boldsymbol {v}}&={\frac {\boldsymbol {u}}{\|{\boldsymbol {u}}\|}}\end{aligned}}}

とする。

ここで、

α = 14 {\displaystyle \alpha =14} であり、 x = a 1 = [ 12 6 4 ] T {\displaystyle {\boldsymbol {x}}={\boldsymbol {a}}_{1}={\begin{bmatrix}12&6&-4\end{bmatrix}}^{\textsf {T}}}

であるから、

u = [ 2 6 4 ] T = [ 2 ] [ 1 3 2 ] T {\displaystyle {\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 v = 1 14 [ 1 3 2 ] T {\displaystyle {\boldsymbol {v}}={1 \over {\sqrt {14}}}{\begin{bmatrix}-1&3&-2\end{bmatrix}}^{\textsf {T}}} であり、
Q 1 = I 2 14 14 [ 1 3 2 ] [ 1 3 2 ] = I 1 7 [ 1 3 2 3 9 6 2 6 4 ] = [ 6 / 7 3 / 7 2 / 7 3 / 7 2 / 7 6 / 7 2 / 7 6 / 7 3 / 7 ] {\displaystyle {\begin{aligned}Q_{1}={}&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{aligned}}}

である。

今、

Q 1 A = [ 14 21 14 0 49 14 0 168 77 ] {\displaystyle Q_{1}A={\begin{bmatrix}14&21&-14\\0&-49&-14\\0&168&-77\end{bmatrix}}}

を見てみると、すでにほぼ三角行列である。あとは(3, 2)要素を零にするだけでよい。

(1, 1)における小行列を取り、同じプロセスを

A = M 11 = [ 49 14 168 77 ] {\displaystyle A'=M_{11}={\begin{bmatrix}-49&-14\\168&-77\end{bmatrix}}}

に再び適用する。

先述のメソッドと同様にして、このプロセスの次のステップが正しく動作するために、1で直和を取ることにより、ハウスホルダー変換

Q 2 = [ 1 0 0 0 7 / 25 24 / 25 0 24 / 25 7 / 25 ] {\displaystyle Q_{2}={\begin{bmatrix}1&0&0\\0&-7/25&24/25\\0&24/25&7/25\end{bmatrix}}}

を得る。

今、

Q = Q 1 T Q 2 T = [ 6 / 7 69 / 175 58 / 175 3 / 7 158 / 175 6 / 175 2 / 7 6 / 35 33 / 35 ] {\displaystyle 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}}}

または、有効数字四桁で、

Q = Q 1 T Q 2 T = [ 0.8571 0.3943 0.3314 0.4286 0.9029 0.0343 0.2857 0.1714 0.9429 ] R = Q 2 Q 1 A = Q T A = [ 14 21 14 0 175 70 0 0 35 ] {\displaystyle {\begin{aligned}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{aligned}}}

を得た。

行列Qは直交行列であり、Rは上三角行列であるため、A = QRは求めるべきQR分解である。

利点と欠点

ハウスホルダー変換の使用は、R行列のゼロを生成するメカニズムに鏡映を利用しているため、最もシンプルでかつ数値的に安定したQR分解アルゴリズムである。しかしながら、新しい零要素を生成する毎回の鏡映変化において行列QR両方の行列全体を書き換えるため、ハウスホルダー変換は必要とするメモリ帯域幅が多く、並列化できない。(ただし鏡映変換のブロック化に基づくブロック化ハウスホルダ変換を使えばメモリ帯域幅への要求を低減することができる)

ギブンス回転の使用

QR分解はギブンス回転を使用しても計算できる。各回転により行列の亜対角要素がゼロになり、R行列を構成できる。すべてのギブンス回転を結合することで直交行列Qを構成できる。

実際には、行列全体を構成して乗算をするようなギブンス回転は行われない。代わりに、疎な要素を計算するような無駄な計算をしない、疎なギブンス行列乗算と同等なあるギブンス回転の手順が採られる。そのギブンス回転の手順は少しの非対角成分をゼロにするだけで済み、ハウスホルダー変換よりも容易に並列化できる。

A = [ 12 51 4 6 167 68 4 24 41 ] {\displaystyle A={\begin{bmatrix}12&-51&4\\6&167&-68\\-4&24&-41\end{bmatrix}}}

の分解を考える。

まず、左下隅の要素、 a 31 = 4 {\displaystyle {\boldsymbol {a}}_{31}=-4} をゼロにする回転行列を構成する必要がある。この行列 G 1 {\displaystyle G_{1}} はギブンス回転で求めることができる。まずベクトル [ 12 4 ] {\displaystyle {\begin{bmatrix}12&-4\end{bmatrix}}} を、X軸に沿って回転させる。このベクトルは角度 θ = arctan [ ( 4 ) 12 ] {\displaystyle \theta =\arctan \left[{-(-4) \over 12}\right]} を持つ。直交ギブンス回転行列 G 1 {\displaystyle G_{1}} を次のように作る。

G 1 = [ cos ( θ ) 0 sin ( θ ) 0 1 0 sin ( θ ) 0 cos ( θ ) ] [ 0.94868 0 0.31622 0 1 0 0.31622 0 0.94868 ] {\displaystyle {\begin{aligned}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{aligned}}}

ここで G 1 A {\displaystyle G_{1}A} の結果は a 31 {\displaystyle {\boldsymbol {a}}_{31}} 要素がゼロである。

G 1 A [ 12.64911 55.97231 16.76007 6 167 68 0 6.64078 37.6311 ] {\displaystyle G_{1}A\approx {\begin{bmatrix}12.64911&-55.97231&16.76007\\6&167&-68\\0&6.64078&-37.6311\end{bmatrix}}}

同様にしてそれぞれ非対角要素 a 21 {\displaystyle a_{21}} a 32 {\displaystyle a_{32}} 要素がゼロであるようなギブンス行列 G 2 {\displaystyle G_{2}} G 3 {\displaystyle G_{3}} を構成し、三角行列 R {\displaystyle R} を構成する。直交行列 Q T {\displaystyle Q^{\textsf {T}}} はすべてのギブンス行列の積 Q T = G 3 G 2 G 1 {\displaystyle Q^{\textsf {T}}=G_{3}G_{2}G_{1}} で表される。したがって、 G 3 G 2 G 1 A = Q T A = R {\displaystyle G_{3}G_{2}G_{1}A=Q^{\textsf {T}}A=R} であり、QR分解は A = Q R {\displaystyle A=QR} である。

利点と欠点

ギブンス回転によるQR分解は、アルゴリズムを完全に動かすのに必要な行の順序を決定するのが簡単ではないため、実装に最も手間がかかる。しかしながら、新しいゼロ要素 a i j {\displaystyle a_{ij}} がゼロになる予定の要素の行(i)とその上の行(j)にしか影響しないという特筆すべき利点がある。これにより、ギブンス回転アルゴリズムはハウスホルダー変換手法よりも帯域幅効率が良く、容易に並列化できる。

行列式や固有値の積との関係

QR分解を正方行列の行列式の絶対値を求めるのに利用できる。ある行列が A = Q R {\displaystyle A=QR} と分解できるとする。このとき

det ( A ) = det ( Q ) det ( R ) {\displaystyle \det(A)=\det(Q)\cdot \det(R)}

である。

Qはユニタリであるため、 | det ( Q ) | = 1 {\displaystyle |\det(Q)|=1} である。したがって、 r i i {\displaystyle r_{ii}} Rの対角要素とすると、

| det ( A ) | = | det ( R ) | = | i r i i | {\displaystyle \left|\det(A)\right|=\left|\det(R)\right|=\left|\prod _{i}r_{ii}\right|}

となる。

さらに、行列式は固有値の積に等しいため、 λ i {\displaystyle \lambda _{i}} A {\displaystyle A} の固有値とすると、

| i r i i | = | i λ i | {\displaystyle \left|\prod _{i}r_{ii}\right|=\left|\prod _{i}\lambda _{i}\right|}

となる。

QR分解の定義を非正方行列に導入し、固有値を特異値に置き換えることで、上記性質を非正方行列 A {\displaystyle A} に拡張することができる。

非正方行列AのQR分解を

A = Q [ R O ] , Q Q = I {\displaystyle A=Q{\begin{bmatrix}R\\O\end{bmatrix}},\qquad Q^{*}Q=I}

とする。ただし、 O {\displaystyle O} は零行列、 Q {\displaystyle Q} はユニタリ行列。

特異値分解と行列式の性質から、 σ i {\displaystyle \sigma _{i}} A {\displaystyle A} の特異値として、

| i r i i | = i σ i {\displaystyle \left|\prod _{i}r_{ii}\right|=\prod _{i}\sigma _{i}}

A {\displaystyle A} R {\displaystyle R} の特異値は同じであるが、複素固有値が異なる場合があることに注意すること。しかしながら、Aが正方ならば、下記は真である。

i σ i = | i λ i | {\displaystyle {\prod _{i}\sigma _{i}}=\left|{\prod _{i}\lambda _{i}}\right|}

結論として、QR分解を使うことによって行列の固有値や特異値の積を効率よく計算することができる。

列のピボット

ピボットQRは列のピボット[4]の新しいステップにおいて、それぞれ初めに残りの列で最も大きいものを取るという点で、通常のグラム・シュミット法とは異なっている。したがって、置換行列Pを次のように導入する。

A P = Q R A = Q R P T {\displaystyle AP=QR\quad \iff \quad A=QRP^{\textsf {T}}}

列のピボットはAが(ほぼ)階数落ちである、またはその疑いがある場合に便利である。また、数値的精度を向上させることもできる。通常、Rの対角成分が非増加、つまり | r 11 | | r 22 | | r n n | {\displaystyle \left|r_{11}\right|\geq \left|r_{22}\right|\geq \ldots \geq \left|r_{nn}\right|} となるようにPを選ぶ。この手段により特異値分解よりも低い計算コストでAの(数値的な)階数を求めることができ、いわゆるRank Revealing QR分解(英語版)の基礎となっている。

線形逆問題への利用

行列の逆行列を直接求めるのに比べ、QR分解を利用した逆問題の解法は、条件数が減少していることからも分かるように、数値的に安定している[5]

次元が m × n {\displaystyle m\times n} で階数が m {\displaystyle m} であるような行列 A {\displaystyle A} に対して、劣決定( m < n {\displaystyle m<n} )線形問題 A x = b {\displaystyle Ax=b} を解くためには、まず A {\displaystyle A} の転置行列のQR分解 A T = Q R {\displaystyle A^{\textsf {T}}=QR} を求める。ただし、Qは直交行列(つまり、 Q T = Q 1 {\displaystyle Q^{\textsf {T}}=Q^{-1}} )であり、Rは R = [ R 1 0 ] {\displaystyle R={\begin{bmatrix}R_{1}\\0\end{bmatrix}}} という特殊な形である。ここで、 R 1 {\displaystyle R_{1}} m × m {\displaystyle m\times m} 正方右三角行列、零行列は ( n m ) × m {\displaystyle (n-m)\times m} 次元である。計算すると、この逆問題の解を次のように表すことができる。 x = Q [ ( R 1 T ) 1 b 0 ] {\displaystyle x=Q{\begin{bmatrix}\left(R_{1}^{\textsf {T}}\right)^{-1}b\\0\end{bmatrix}}}

ここで、 R 1 1 {\displaystyle R_{1}^{-1}} ガウスの消去法で計算でき、 ( R 1 T ) 1 b {\displaystyle \left(R_{1}^{\textsf {T}}\right)^{-1}b} 前方置換法(英語版)を用いることで直接計算できる。後者の手法の方が数値的精度が高く、計算量も少ないという利点がある。

ノルム A x ^ b {\displaystyle \|A{\hat {x}}-b\|} を最小にするような過決定( m n {\displaystyle m\geq n} )問題 A x = b {\displaystyle Ax=b} の解 x ^ {\displaystyle {\hat {x}}} を求めるためには、まず A {\displaystyle A} のQR分解 A = Q R {\displaystyle A=QR} を求める。 Q 1 {\displaystyle Q_{1}} を直交行列 Q {\displaystyle Q} 全体のうち最初の n {\displaystyle n} 列を含む m × n {\displaystyle m\times n} 行列、 R 1 {\displaystyle R_{1}} を先述の通りに置くと、この問題の解は x ^ = R 1 1 ( Q 1 T b ) {\displaystyle {\hat {x}}=R_{1}^{-1}\left(Q_{1}^{\textsf {T}}b\right)} と表せる。劣決定の場合と同様に、 R 1 {\displaystyle R_{1}} の逆行列を直接計算しなくても、後方置換法(英語版)を用いることで早く正確に x ^ {\displaystyle {\hat {x}}} を求めることができる。( Q 1 {\displaystyle Q_{1}} R 1 {\displaystyle R_{1}} は数値ライブラリによっては高速な(economic)QR分解として実装されている。)

一般化

岩澤分解はQR分解を半単純リー群に一般化している。

脚注

[脚注の使い方]
  1. ^ Golub & Van Loan 2013, 5.2 The QR Factorization.
  2. ^ Golub & Van Loan 2013, Theorem 5.2.1 (QR Factorization).
  3. ^ a b c L. N. Trefethen and D. Bau, Numerical Linear Algebra (SIAM, 1997).
  4. ^ Strang, Gilbert (2019). Linear Algebra and Learning from Data (1 ed.). Wellesley: Wellesley Cambridge Press. p. 143. ISBN 978-0-692-19638-0 
  5. ^ Parker, Geophysical Inverse Theory, Ch1.13.

参考文献

和文

英文

  • Golub, Gene H.; Van Loan, Charles F. (2013) (英語). Matrix computations (Fourth ed.). Johns Hopkins University Press. ISBN 978-1421407944. MR3024913. https://books.google.co.jp/books?id=X5YfsuCWpxMC 
  • Trefethen, Lloyd N.; Bau III, David (1997) (英語). Numerical Linear Algebra. Society for Industrial and Applied Mathematics. ISBN 9780898713619 
  • Horn, Roger A.; Johnson, Charles R. (1985). “Section 2.8.” (英語). Matrix Analysis. Cambridge University Press. ISBN 0-521-38632-2 
  • Press, William H.; Teukolsky, Saul A.; Vetterling, William T.; Flannery, Brian P. (2007). “Section 2.10. QR Decomposition” (英語). Numerical Recipes: The Art of Scientific Computing (3rd ed.). New York: Cambridge University Press. ISBN 978-0-521-88068-8 
  • Stoer, Josef; Bulirsch, Roland (2002) (英語). Introduction to Numerical Analysis. Texts in applied mathematics, 12 (3rd ed.). Springer. ISBN 0-387-95452-X 
  • Golub, Gene H.; Van Loan, Charles F. (1996), Matrix Computations (3rd ed.), Johns Hopkins, ISBN 978-0-8018-5414-9 .

関連項目

外部リンク

  • Online Matrix Calculator 行列のQR分解の実行
  • LAPACK users manual QR分解の計算のサブルーチンの詳細
  • Mathematica users manual QR分解の計算のルーチンの詳細と例
  • ALGLIB LAPACKのC++、C#、Delphiなどへの移植
  • Eigen::QR QR分解のC++での実装
連立一次方程式
ベクトル
ベクトル空間
計量ベクトル空間
行列線型写像
演算・操作
不変量
クラス
行列式
多重線型代数
数値線形代数
基本的な概念
ソフトウェア
ライブラリ
反復法・技法
人物
行列値関数
その他
カテゴリ カテゴリ