QR分解

**QR分解**（キューアールぶんかい、）とは、*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&nbsp;≤&nbsp;*k*&nbsp;≤&nbsp;*n*の任意の*k*について、*Q*の最初の*k*列は*A*の最初の*k*列の線型包をなす*。A*の任意の列*k*が*Q*の最初の*k*列にのみ依存するということは、*R*が三角行列であることから明らかである。

### 矩形行列
より一般的に、*m*&nbsp;≥&nbsp;*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 


ここで、*R*1は*n*×*n*上三角行列、*0*は(*m*&nbsp;−&nbsp;*n*)×*n*零行列、*Q*1は*m*×*n*行列、*Q*2は*m*×(*m*&nbsp;−&nbsp;*n*)行列で、*Q*1と*Q*2は両方直交する列を持つ。

*Q*1*R*1をはAの薄い(*thin*)QR分解と呼び、 Trefethen & Bauは軽減(*reduced*)QR分解と呼んでいる。もし*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 = \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分解との関係
RQ分解は行列*A*を上三角行列*R*(別名右三角)と直交行列*Q*に変換する。QR分解との違いはこれらの行列の順番だけである。

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

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

#### 利点と欠点

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

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

ハウスホルダー変換 (または*ハウスホルダー鏡映*)はベクトルを取り、平面または超平面に関する鏡映をする変換である。この演算は*m*×*n*行列A (*m*&nbsp;≥&nbsp;*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**の最初の列を選んで得られるハウスホルダー行列*Q*1に*A*を乗算する。この結果、行列*Q*1*A*は左の列が(最初の行を除き)ゼロになる。
: Q_1A = \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 = \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*&nbsp;−&nbsp;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}
      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{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*の(数値的な)階数を求めることができ、いわゆるの基礎となっている。

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

次元がm \times nで階数がmであるような行列Aに対して、劣決定(m )線形問題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 \\
                           0
  \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分解を半単純リー群に一般化している。

## 脚注

## 参考文献
### 和文

* 
* 
* 

### 英文

* 
* 
* 
* 
* 
* .

## 関連項目
* QR法
* 極分解
* 固有値分解
* LU分解
* 特異値分解

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