QR分解

テンプレート:参照方法 テンプレート:ページ番号 テンプレート:翻訳直後 QR分解（キューアールぶんかい、テンプレート:Lang-en-short）とは、m × n 実行列 Aを、 m 次直交行列 Q と m × n 上三角行列 R との積への分解により表すこと、またはそう表した表現をいうテンプレート:Sfn。このような分解は常に存在するテンプレート:Sfn。

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

定義

正方行列

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

$A=QR$

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

もしAが複素正方行列ならば、Qがユニタリ行列 (つまり$Q^{*}Q=Q{Q}^{*}=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\left[\begin{array}{c}{R}_1\\ 0\end{array}\right]=\left[\begin{array}{c}{Q}_1,Q_2\end{array}\right]\left[\begin{array}{c}{R}_1\\ 0\end{array}\right]=Q_1{R}_1$

ここで、R₁はn×n上三角行列、0は(m − n)×n零行列、Q₁はm×n行列、Q₂はm×(m − n)行列で、Q₁とQ₂は両方直交する列を持つ。

Q₁R₁をテンプレート:HarvtxtはAの薄い(thin)QR分解と呼び、 Trefethen & Bauは軽減(reduced)QR分解と呼んでいる。もしAが最大階数nであり、R₁の対角成分を正にするならば、R₁とQ₁は一意に定まる。しかし一般的にQ₂はそうではない。R₁はA^* A (Aが実行列の場合A^TAに等しい)のコレスキー分解の上三角部分に等しい。

QL・RQ・LQ分解

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

QR分解の計算

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

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

テンプレート:Details グラム・シュミットの正規直交化法を最大階数行列の列$A=[{𝒂}_1,\dots ,{𝒂}_n]$に適用することを考える。内積$⟨𝒗,𝒘⟩={𝒗}^{\mathsf{\text{T}}}𝒘$ (複素ベクトルの場合 $⟨𝒗,𝒘⟩={𝒗}^{*}𝒘$)とする。

射影の定義より、

${\mathrm{proj}}_{𝒖}𝒂=\frac{⟨𝒖,𝒂⟩}{⟨𝒖,𝒖⟩}𝒖$

したがって、

$\begin{array}{cccc}{𝒖}_1& ={𝒂}_1,& {𝒆}_1& =\frac{{𝒖}_1}{\Vert {𝒖}_1\Vert }\\ {𝒖}_2& ={𝒂}_2-{\mathrm{proj}}_{{𝒖}_1}{𝒂}_2,& {𝒆}_2& =\frac{{𝒖}_2}{\Vert {𝒖}_2\Vert }\\ {𝒖}_3& ={𝒂}_3-{\mathrm{proj}}_{{𝒖}_1}{𝒂}_3-{\mathrm{proj}}_{{𝒖}_2}{𝒂}_3,& {𝒆}_3& =\frac{{𝒖}_3}{\Vert {𝒖}_3\Vert }\\ & ⋮& & ⋮\\ {𝒖}_k& ={𝒂}_k-\sum _{j=1}^{k-1}{\mathrm{proj}}_{{𝒖}_j}{𝒂}_k,& {𝒆}_k& =\frac{{𝒖}_k}{\Vert {𝒖}_k\Vert }\end{array}$

ここで${𝒂}_i$を新しく計算された正規直交基底上に表すことができ、$⟨{𝒆}_i,{𝒂}_i⟩=\Vert {𝒖}_i\Vert$であるから、

$\begin{array}{cc}{𝒂}_1& =⟨{𝒆}_1,{𝒂}_1⟩{𝒆}_1\\ {𝒂}_2& =⟨{𝒆}_1,{𝒂}_2⟩{𝒆}_1+⟨{𝒆}_2,{𝒂}_2⟩{𝒆}_2\\ {𝒂}_3& =⟨{𝒆}_1,{𝒂}_3⟩{𝒆}_1+⟨{𝒆}_2,{𝒂}_3⟩{𝒆}_2+⟨{𝒆}_3,{𝒂}_3⟩{𝒆}_3\\ & ⋮\\ {𝒂}_k& =\sum _{j=1}^k⟨{𝒆}_j,{𝒂}_k⟩{𝒆}_j\end{array}$

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

$A=QR$

ただし、

$Q=[{𝒆}_1,\dots ,{𝒆}_n],$

$R=\left[\begin{array}{cccc}⟨{𝒆}_1,{𝒂}_1⟩& ⟨{𝒆}_1,{𝒂}_2⟩& ⟨{𝒆}_1,{𝒂}_3⟩& \dots \\ 0& ⟨{𝒆}_2,{𝒂}_2⟩& ⟨{𝒆}_2,{𝒂}_3⟩& \dots \\ 0& 0& ⟨{𝒆}_3,{𝒂}_3⟩& \dots \\ ⋮& ⋮& ⋮& \ddots \end{array}\right]$

例

$A=\left[\begin{array}{ccc}12& -51& 4\\ 6& 167& -68\\ -4& 24& -41\end{array}\right]$

の分解を考える。

正規直交行列$Q$に対して、

$Q^{\mathsf{\text{T}}}Q=I$

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

$\begin{array}{cc}U=\left[\begin{array}{ccc}{𝒖}_1& {𝒖}_2& {𝒖}_3\end{array}\right]& =\left[\begin{array}{ccc}12& -69& -58/5\\ 6& 158& 6/5\\ -4& 30& -33\end{array}\right];\\ Q=\left[\begin{array}{ccc}\frac{{𝒖}_1}{\Vert {𝒖}_1\Vert }& \frac{{𝒖}_2}{\Vert {𝒖}_2\Vert }& \frac{{𝒖}_3}{\Vert {𝒖}_3\Vert }\end{array}\right]& =\left[\begin{array}{ccc}6/7& -69/175& -58/175\\ 3/7& 158/175& 6/175\\ -2/7& 6/35& -33/35\end{array}\right]\end{array}$

したがって、

$\begin{array}{cc}{Q}^{\mathsf{\text{T}}}A& =Q^{\mathsf{\text{T}}}QR=R;\\ R& =Q^{\mathsf{\text{T}}}A=\left[\begin{array}{ccc}14& 21& -14\\ 0& 175& -70\\ 0& 0& 35\end{array}\right]\end{array}$

RQ分解との関係

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

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

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

利点と欠点

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

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

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

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

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

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

$\alpha =-\exp (i\arg x_k)\Vert 𝒙\Vert$

テンプレート:Harvとして、さらに以下のQの導出において転置を共役転置に読み替えること。

ここで、${𝒆}_1$をベクトル(1 0 … 0)^T、||·||をユークリッド距離、$I$をm×m単位行列とし、

$\begin{array}{cc}𝒖& =𝒙-\alpha {𝒆}_1,\\ 𝒗& =\frac{𝒖}{\Vert 𝒖\Vert },\\ Q& =I-2𝒗{𝒗}^{\mathsf{\text{T}}}\end{array}$

あるいは、$A$が複素行列ならば、

$Q=I-2𝒗{𝒗}^{*}$

とおく。

$Q$はm×mハウスホルダー行列であり、

$Q𝒙={\left[\begin{array}{cccc}\alpha & 0& \cdots & 0\end{array}\right]}^{\mathsf{\text{T}}}$

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

$Q_{1}A=\left[\begin{array}{cccc}{\alpha }_1& \star & \dots & \star \\ 0& & & \\ ⋮& & A^{\prime }& \\ 0& & & \end{array}\right]$

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

$Q_k=\left[\begin{array}{cc}{I}_{k-1}& 0\\ 0& {Q_k}^{\prime }\end{array}\right]$

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

$R=Q_t\cdots Q_2{Q}_{1}A$

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

$Q=Q_1^{\mathsf{\text{T}}}{Q}_2^{\mathsf{\text{T}}}\cdots Q_t^{\mathsf{\text{T}}}$

とすると、 $A=QR$は$A$のQR分解である。

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

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

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

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

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

と表せる。

例

$A=\left[\begin{array}{ccc}12& -51& 4\\ 6& 167& -68\\ -4& 24& -41\end{array}\right]$

の分解を考える。

まず、行列Aの最初の列、ベクトル${𝒂}_1={\left[\begin{array}{ccc}12& 6& -4\end{array}\right]}^{\mathsf{\text{T}}}$を、$\Vert {𝒂}_1\Vert {𝒆}_1={\left[\begin{array}{ccc}1& 0& 0\end{array}\right]}^{\mathsf{\text{T}}}$に変換する鏡映を見つける必要がある。

今、

$\begin{array}{cc}𝒖& =𝒙-\alpha {𝒆}_1,\\ 𝒗& =\frac{𝒖}{\Vert 𝒖\Vert }\end{array}$

とする。

ここで、

$\alpha =14$であり、$𝒙={𝒂}_1={\left[\begin{array}{ccc}12& 6& -4\end{array}\right]}^{\mathsf{\text{T}}}$

であるから、

$𝒖={\left[\begin{array}{ccc}-2& 6& -4\end{array}\right]}^{\mathsf{\text{T}}}=\left[\begin{array}{c}2\end{array}\right]{\left[\begin{array}{ccc}-1& 3& -2\end{array}\right]}^{\mathsf{\text{T}}}$ and $𝒗=\frac{1}{\sqrt{14}}{\left[\begin{array}{ccc}-1& 3& -2\end{array}\right]}^{\mathsf{\text{T}}}$であり、

$\begin{array}{cc}{Q}_1=& I-\frac{2}{\sqrt{14}\sqrt{14}}\left[\begin{array}{c}-1\\ 3\\ -2\end{array}\right]\left[\begin{array}{ccc}-1& 3& -2\end{array}\right]\\ =& I-\frac{1}{7}\left[\begin{array}{ccc}1& -3& 2\\ -3& 9& -6\\ 2& -6& 4\end{array}\right]\\ =& \left[\begin{array}{ccc}6/7& 3/7& -2/7\\ 3/7& -2/7& 6/7\\ -2/7& 6/7& 3/7\end{array}\right]\end{array}$

である。

今、

$Q_{1}A=\left[\begin{array}{ccc}14& 21& -14\\ 0& -49& -14\\ 0& 168& -77\end{array}\right]$

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

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

$A^{\prime }=M_{11}=\left[\begin{array}{cc}-49& -14\\ 168& -77\end{array}\right]$

に再び適用する。

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

$Q_2=\left[\begin{array}{ccc}1& 0& 0\\ 0& -7/25& 24/25\\ 0& 24/25& 7/25\end{array}\right]$

を得る。

今、

$Q=Q_1^{\mathsf{\text{T}}}{Q}_2^{\mathsf{\text{T}}}=\left[\begin{array}{ccc}6/7& -69/175& 58/175\\ 3/7& 158/175& -6/175\\ -2/7& 6/35& 33/35\end{array}\right]$

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

$\begin{array}{cc}Q& =Q_1^{\mathsf{\text{T}}}{Q}_2^{\mathsf{\text{T}}}=\left[\begin{array}{ccc}0.8571& -0.3943& 0.3314\\ 0.4286& 0.9029& -0.0343\\ -0.2857& 0.1714& 0.9429\end{array}\right]\\ R& =Q_2{Q}_{1}A=Q^{\mathsf{\text{T}}}A=\left[\begin{array}{ccc}14& 21& -14\\ 0& 175& -70\\ 0& 0& -35\end{array}\right]\end{array}$

を得た。

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

利点と欠点

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

ギブンス回転の使用

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

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

例

$A=\left[\begin{array}{ccc}12& -51& 4\\ 6& 167& -68\\ -4& 24& -41\end{array}\right]$

の分解を考える。

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

$\begin{array}{cc}{G}_1& =\left[\begin{array}{ccc}\cos (\theta )& 0& -\sin (\theta )\\ 0& 1& 0\\ \sin (\theta )& 0& \cos (\theta )\end{array}\right]\\ & \approx \left[\begin{array}{ccc}0.94868& 0& -0.31622\\ 0& 1& 0\\ 0.31622& 0& 0.94868\end{array}\right]\end{array}$

ここで$G_{1}A$の結果は${𝒂}_{31}$要素がゼロである。

$G_{1}A\approx \left[\begin{array}{ccc}12.64911& -55.97231& 16.76007\\ 6& 167& -68\\ 0& 6.64078& -37.6311\end{array}\right]$

同様にしてそれぞれ非対角要素$a_{21}$・$a_{32}$要素がゼロであるようなギブンス行列$G_2$・$G_3$を構成し、三角行列$R$を構成する。直交行列$Q^{\mathsf{\text{T}}}$はすべてのギブンス行列の積$Q^{\mathsf{\text{T}}}=G_3{G}_2{G}_1$で表される。したがって、$G_3{G}_2{G}_{1}A=Q^{\mathsf{\text{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の対角要素とすると、

$|\det (A)|=|\det (R)|=\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\left[\begin{array}{c}R\\ O\end{array}\right],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⟺A=QR{P}^{\mathsf{\text{T}}}$

列のピボットはAが(ほぼ)階数落ちである、またはその疑いがある場合に便利である。また、数値的精度を向上させることもできる。通常、Rの対角成分が非増加、つまり$\left|r_{11}\right|\ge \left|r_{22}\right|\ge \dots \ge \left|r_{nn}\right|$となるようにPを選ぶ。この手段により特異値分解よりも低い計算コストでAの(数値的な)階数を求めることができ、いわゆるテンプレート:仮リンクの基礎となっている。

線形逆問題への利用

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

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

ここで、$R_1^{-1}$はガウスの消去法で計算でき、${\left(R_1^{\mathsf{\text{T}}}\right)}^{-1}b$はテンプレート:仮リンクを用いることで直接計算できる。後者の手法の方が数値的精度が高く、計算量も少ないという利点がある。

ノルム$\Vert A\hat{x}-b\Vert$を最小にするような過決定($m\ge 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^{\mathsf{\text{T}}}b\right)$と表せる。劣決定の場合と同様に、$R_1$の逆行列を直接計算しなくても、テンプレート:仮リンクを用いることで早く正確に$\hat{x}$を求めることができる。($Q_1$と$R_1$は数値ライブラリによっては高速な(economic)QR分解として実装されている。)

一般化

テンプレート:仮リンクはQR分解を半単純リー群に一般化している。

参考文献

和文

• テンプレート:Cite book
• テンプレート:Cite book
• テンプレート:Cite book

英文

• テンプレート:Cite book
• テンプレート:Cite book
• テンプレート:Cite book
• テンプレート:Cite book
• テンプレート:Cite book
• テンプレート:Citation.

関連項目

• QR法
• テンプレート:仮リンク
• 固有値分解
• LU分解
• 特異値分解

外部リンク

• Online Matrix Calculator 行列のQR分解の実行
• LAPACK users manual QR分解の計算のサブルーチンの詳細
• Mathematica users manual QR分解の計算のルーチンの詳細と例
• ALGLIB LAPACKのC++、C#、Delphiなどへの移植
• Eigen::QR QR分解のC++での実装