QR分解 (キューアールぶんかい、英 : QR decomposition, QR factorization )とは、m × n 実行列 A を、 m 次直交行列 Q と m × n 上三角行列 R との積への分解により表すこと、またはそう表した表現をいう。このような分解は常に存在する。
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 が存在する。
もしA がn 個の線形独立 な列を持つなら、Q の最初のn 列はA の列空間 の正規直交基底 をなす。より一般的に、1 ≤ k ≤ n の任意のk について、Q の最初のk 列はA の最初のk 列の線型包 をなす[3] 。A の任意の列k がQ の最初のk 列にのみ依存するということは、R が三角行列であることから明らかである[3] 。
矩形行列 より一般的に、m ≥ n である複素m ×n 行列A を、m ×m ユニタリ行列 Q とm ×n 上三角行列R に分解することができる。m ×n 上三角行列の下から(m −n )行はすべてゼロであるため、R や、R とQ 両方の分割を簡単に行うことができる。
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}} ここで、R 1 はn ×n 上三角行列、0 は(m − n )×n 零行列、Q 1 はm ×n 行列、Q 2 はm ×(m − n )行列で、Q 1 とQ 2 は両方直交する列を持つ。
Q 1 R 1 をGolub & Van Loan (1996 , §5.2)はAの薄い(thin )QR分解と呼び、 Trefethen & Bauは軽減(reduced )QR分解と呼んでいる[3] 。もしA が最大階数 n であり、R 1 の対角成分を正にするならば、R 1 とQ 1 は一意に定まる。しかし一般的にQ 2 はそうではない。R 1 はA * A (A が実行列の場合A T A に等しい)のコレスキー分解 の上三角部分に等しい。
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 の最初の列を選んで得られるハウスホルダー行列Q 1 にA を乗算する。この結果、行列Q 1 A は左の列が(最初の行を除き)ゼロになる。
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 ′ (Q 1 A から最初の列と最初の行を除いたもの)に繰り返すと、ハウスホルダー行列Q ′2 が得られる。Q ′2 はQ 1 より小さいということに注意すること。A ′の代わりにQ 1 A で計算したいため、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分解アルゴリズムである。しかしながら、新しい零要素を生成する毎回の鏡映変化において行列Q とR 両方の行列全体を書き換えるため、ハウスホルダー変換は必要とするメモリ帯域幅が多く、並列化できない。(ただし鏡映変換のブロック化に基づくブロック化ハウスホルダ変換を使えばメモリ帯域幅への要求を低減することができる)
ギブンス回転の使用 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分解を半単純リー群に一般化している。
脚注 [脚注の使い方 ]
^ a b c L. N. Trefethen and D. Bau, Numerical Linear Algebra (SIAM, 1997). ^ Strang, Gilbert (2019). Linear Algebra and Learning from Data (1 ed.). Wellesley: Wellesley Cambridge Press. p. 143. ISBN 978-0-692-19638-0 ^ 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. MR 3024913. 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++での実装