LQ分解とは

LQ分解とは，任意の$m\times n$行列$\boldsymbol{A}$を

\boldsymbol{A}=\boldsymbol{L}\boldsymbol{Q}

または

\boldsymbol{A}=\boldsymbol{L}_{r}\boldsymbol{Q}_{r}

と分解する操作のことです．ただし，$\boldsymbol{L}$は$m\times n$下三角行列，$\boldsymbol{Q}$は$n\times n$正規直交行列，$\boldsymbol{L}_{r}$は$\boldsymbol{L}$の左$m\times r$部分行列，$\boldsymbol{Q}_{r}$は$\boldsymbol{Q}$の上$r\times n$部分行列，$r$は$\boldsymbol{A}$のランク（階数）です．
後者はいわゆる階数分解（任意の行列を列フルランクな行列と行フルランクな行列の積の形に分解する操作のこと）の一種となっています．

正規直交行列とは，それを構成する基底ベクトルが互いに正規直交性，すなわち$\boldsymbol{Q}^{\mathrm{T}}\boldsymbol{Q}=\boldsymbol{Q}\boldsymbol{Q}^{\mathrm{T}}=\boldsymbol{1}$を満たす行列のことです（ただし$\boldsymbol{1}$は単位行列で，文脈に合わせて任意のサイズをとれるものとします）．

階数分解における列フルランクな行列は像空間写像行列，行フルランクな行列は内部空間への写像行列としての意味を持ちます．後者が正規直交性を持つならば，その直交補空間は核空間（零空間） となります．
LQ分解で現れる$\boldsymbol{Q}$から$\boldsymbol{Q}_{r}$を除いた下$(n-r)\times n$部分行列の転置は，そのまま核空間写像行列（零空間写像行列） となります．
このように，行列により写像される空間の性質を理解する上でLQ分解は重宝します．

LQ分解の具体的な計算方法は幾つか知られていますが，本記事では，高速なGram=Schmidtの直交化による求め方と，行列の退化に対し強いHouseholder法による求め方を紹介します．

なお，LQ分解と同種の計算にQR分解があります．これは，任意の$m\times n$行列$\boldsymbol{A}$を

\boldsymbol{A}=\boldsymbol{Q}\boldsymbol{R}

と分解する操作のことです．ただし，このときの$\boldsymbol{Q}$は$m\times m$正規直交行列，$\boldsymbol{R}$は$m\times n$上三角行列です．
上三角行列を転置すると下三角行列になるので，$\boldsymbol{A}^{\mathrm{T}}$をLQ分解した結果を転置すれば，同じものが得られます．すなわちQR分解とLQ分解は本質的に同じ操作ですので，本記事ではこれ以上の説明は省きます．

Gram=Schmidtの直交化によるLQ分解

Gram=Schmidtの直交化

単位ベクトル$\boldsymbol{e}_{1}$があるとして，それと線形独立なあるベクトル$\boldsymbol{a}$を元に，$\boldsymbol{e}_{1}$に直交する単位ベクトル$\boldsymbol{e}_{2}$を作ることを考えます．任意次元のベクトルを考えるので，外積は使えないものとします．

$\boldsymbol{a}$を$\boldsymbol{e}_{1}$が張る直線に射影したベクトルは$\left(\boldsymbol{a}^{\mathrm{T}}\boldsymbol{e}_{1}\right)\boldsymbol{e}_{1}$ですから，$\boldsymbol{a}$からこれを引けば，$\boldsymbol{e}_{1}$に直交するベクトル

\boldsymbol{n}_{2}=\boldsymbol{a}-\left(\boldsymbol{a}^{\mathrm{T}}\boldsymbol{e}_{1}\right)\boldsymbol{e}_{1}

となります．念のため，$\boldsymbol{n}_{2}^{\mathrm{T}}\boldsymbol{e}_{1}=0$であることは明らかです．
$\boldsymbol{e}_{2}$は，$\boldsymbol{n}_{2}$を正規化して

\boldsymbol{e}_{2}=\frac{\boldsymbol{n}_{2}}{\|\boldsymbol{n}_{2}\|}

とすれば良いでしょう．

同様に，二つの互いに直交する単位ベクトル$\boldsymbol{e}_{1}$，$\boldsymbol{e}_{2}$があるとして，それらと線形独立なあるベクトル$\boldsymbol{a}$を元に，$\boldsymbol{e}_{1}$，$\boldsymbol{e}_{2}$のどちらにも直交する単位ベクトル$\boldsymbol{e}_{3}$を作ることを考えます．

$\boldsymbol{a}$を$\boldsymbol{e}_{1}$と$\boldsymbol{e}_{2}$が張る平面に射影したベクトルは$\left(\boldsymbol{a}^{\mathrm{T}}\boldsymbol{e}_{1}\right)\boldsymbol{e}_{1}+\left(\boldsymbol{a}^{\mathrm{T}}\boldsymbol{e}_{2}\right)\boldsymbol{e}_{2}$ですから，$\boldsymbol{a}$からこれを引けば，$\boldsymbol{e}_{1}$と$\boldsymbol{e}_{2}$の両方に直交するベクトル

\boldsymbol{n}_{3}=\boldsymbol{a}-\left(\boldsymbol{a}^{\mathrm{T}}\boldsymbol{e}_{1}\right)\boldsymbol{e}_{1}-\left(\boldsymbol{a}^{\mathrm{T}}\boldsymbol{e}_{2}\right)\boldsymbol{e}_{2}

となります．これも念のため，$\boldsymbol{n}_{3}^{\mathrm{T}}\boldsymbol{e}_{1}=\boldsymbol{n}_{3}^{\mathrm{T}}\boldsymbol{e}_{2}=0$であることは明らかです．$\boldsymbol{e}_{3}$は，$\boldsymbol{n}_{3}$を正規化して

\boldsymbol{e}_{3}=\frac{\boldsymbol{n}_{3}}{\|\boldsymbol{n}_{3}\|}

とすれば良いです．

これを一般化すれば，互いに直交する単位ベクトル$\boldsymbol{e}_{1},\cdots,\boldsymbol{e}_{k-1}$，およびそれらと線形独立なベクトル$\boldsymbol{a}$から，$\boldsymbol{e}_{1},\cdots,\boldsymbol{e}_{k-1}$全てに直交する単位ベクトル$\boldsymbol{e}_{k}$が，

\displaylines{
\boldsymbol{n}_{k}=\boldsymbol{a}
-\left(\boldsymbol{a}^{\mathrm{T}}\boldsymbol{e}_{1}\right)\boldsymbol{e}_{1}
-\cdots
-\left(\boldsymbol{a}^{\mathrm{T}}\boldsymbol{e}_{k-1}\right)\boldsymbol{e}_{k-1}
\\
\boldsymbol{e}_{k}=\frac{\boldsymbol{n}_{k}}{\|\boldsymbol{n}_{k}\|}
}

という計算によって作れます．

今，互いに線形独立な$m$個の$n$次元ベクトル$\boldsymbol{a}_{1},\cdots,\boldsymbol{a}_{m}$があるとします．上記を踏まえれば，まず$\boldsymbol{a}_{1}$に平行な単位ベクトル$\boldsymbol{e}_{1}$を

\displaylines{
\boldsymbol{n}_{1}=\boldsymbol{a}_{1}
\\
\boldsymbol{e}_{1}=\frac{\boldsymbol{n}_{1}}{\|\boldsymbol{n}_{1}\|}
}

と定め，そこから$k=2,\cdots,m$について順次

\displaylines{
\boldsymbol{n}_{k}=\boldsymbol{a}_{k}
-\left(\boldsymbol{a}_{k}^{\mathrm{T}}\boldsymbol{e}_{1}\right)\boldsymbol{e}_{1}
-\cdots
-\left(\boldsymbol{a}_{k}^{\mathrm{T}}\boldsymbol{e}_{k-1}\right)\boldsymbol{e}_{k-1}
\\
\boldsymbol{e}_{k}=\frac{\boldsymbol{n}_{k}}{\|\boldsymbol{n}_{k}\|}
}

とすることで，互いに直交する$m$個の単位ベクトル$\boldsymbol{e}_{1},\cdots,\boldsymbol{e}_{m}$を作ることができます．これがGram=Schmidtの直交化です．

Gram=Schmidtの直交化とLQ分解との関係

簡単のため$m=n$であるとして，上記の一連の式を変形すれば，

\begin{cases}
\boldsymbol{a}_{1}=\|\boldsymbol{n}_{1}\|\boldsymbol{e}_{1}
\\
\boldsymbol{a}_{2}=\left(\boldsymbol{a}_{2}^{\mathrm{T}}\boldsymbol{e}_{1}\right)\boldsymbol{e}_{1}+\|\boldsymbol{n}_{2}\|\boldsymbol{e}_{2}
\\
\vdots
\\
\boldsymbol{a}_{m}=
\left(\boldsymbol{a}_{m}^{\mathrm{T}}\boldsymbol{e}_{1}\right)\boldsymbol{e}_{1}
+\cdots
+\left(\boldsymbol{a}_{m}^{\mathrm{T}}\boldsymbol{e}_{n-1}\right)\boldsymbol{e}_{m-1}
+\|\boldsymbol{n}_{m}\|\boldsymbol{e}_{m}
\end{cases}

すなわち

\begin{bmatrix}
\boldsymbol{a}_{1} & \boldsymbol{a}_{2} & \cdots & \boldsymbol{a}_{m}
\end{bmatrix}
=
\begin{bmatrix}
\boldsymbol{e}_{1} & \boldsymbol{e}_{2} & \cdots & \boldsymbol{e}_{m}
\end{bmatrix}
\begin{bmatrix}
\|\boldsymbol{n}_{1}\| & \boldsymbol{a}_{2}^{\mathrm{T}}\boldsymbol{e}_{1} & \cdots & \boldsymbol{a}_{m}^{\mathrm{T}}\boldsymbol{e}_{1} \\
0 & \|\boldsymbol{n}_{2}\| & & \boldsymbol{a}_{m}^{\mathrm{T}}\boldsymbol{e}_{2} \\
\vdots & \ddots & \ddots & \vdots \\
0 & \cdots & 0 & \|\boldsymbol{n}_{m}\|
\end{bmatrix}

となっていることが分かります．これを転置すれば，

\begin{bmatrix}
\boldsymbol{a}_{1}^{\mathrm{T}} \\
\boldsymbol{a}_{2}^{\mathrm{T}} \\
\vdots \\
\boldsymbol{a}_{m}^{\mathrm{T}}
\end{bmatrix}
=
\begin{bmatrix}
 \|\boldsymbol{n}_{1}\| & 0 & \cdots & 0 \\
 \boldsymbol{a}_{2}^{\mathrm{T}}\boldsymbol{e}_{1} & \|\boldsymbol{n}_{2}\| & \ddots & \vdots \\
 \vdots & & \ddots & 0 \\
 \boldsymbol{a}_{m}^{\mathrm{T}}\boldsymbol{e}_{1} & \boldsymbol{a}_{m}^{\mathrm{T}}\boldsymbol{e}_{2} & \cdots & \|\boldsymbol{n}_{m}\| \\
\end{bmatrix}
\begin{bmatrix}
\boldsymbol{e}_{1}^{\mathrm{T}} \\
\boldsymbol{e}_{2}^{\mathrm{T}} \\
\vdots \\
\boldsymbol{e}_{m}^{\mathrm{T}}
\end{bmatrix}

となります．右辺の左側行列は下三角行列になっていますから，これはLQ分解された形です．
実際には，前節の手順によって$\boldsymbol{e}_{k}$（$k=1,\cdots,m$）を順次求めていきながら，求まった$\boldsymbol{e}_{k}$に対し行列$\boldsymbol{L}$の$k$列目成分を埋めていく，というやり方になります．

上記は

\begin{bmatrix}
\boldsymbol{a}_{1}^{\mathrm{T}} \\
\boldsymbol{a}_{2}^{\mathrm{T}} \\
\vdots \\
\boldsymbol{a}_{m}^{\mathrm{T}}
\end{bmatrix}
=
\begin{bmatrix}
 \boldsymbol{a}_{1}^{\mathrm{T}}\boldsymbol{e}_{1} & 0 & \cdots & 0 \\
 \boldsymbol{a}_{2}^{\mathrm{T}}\boldsymbol{e}_{1} & \boldsymbol{a}_{2}^{\mathrm{T}}\boldsymbol{e}_{2} & \ddots & \vdots \\
 \vdots & & \ddots & 0 \\
 \boldsymbol{a}_{m}^{\mathrm{T}}\boldsymbol{e}_{1} & \boldsymbol{a}_{m}^{\mathrm{T}}\boldsymbol{e}_{2} & \cdots & \boldsymbol{a}_{m}^{\mathrm{T}}\boldsymbol{e}_{m} \\
\end{bmatrix}
\begin{bmatrix}
\boldsymbol{e}_{1}^{\mathrm{T}} \\
\boldsymbol{e}_{2}^{\mathrm{T}} \\
\vdots \\
\boldsymbol{e}_{m}^{\mathrm{T}}
\end{bmatrix}

と綺麗に書くこともできるのですが，計算手順としては，各$k$列の計算においてまず$\boldsymbol{n}_{k}$を求めるので，これのノルムを用いて対角成分を与える方が効率が良いです．

実装

一般的には$m=n$とは限りません．$\boldsymbol{e}_{k}$は，$\boldsymbol{a}_{k}$の数$m$を超えて作ることは出来ないので，Gram=Schmidtの直交化により行えるLQ分解は，必然的に

\boldsymbol{A}=\boldsymbol{L}_{r}\boldsymbol{Q}_{r}

の形をとることになります．改めて，$r$は行列$\boldsymbol{A}$のランクです．

$r$は，$\boldsymbol{A}$の行ベクトル$\boldsymbol{a}_{1}^{\mathrm{T}},\cdots,\boldsymbol{a}_{m}^{\mathrm{T}}$のうち互いに線形独立なものの個数と言い換えられます．少なくとも，$m$と$n$のうち小さい方の値よりも大きくなることはありませんが，値を事前に知ることは一般的には出来ません．したがって，上記の計算の反復過程でランクを数えていくことになります．
具体的には，ある$k$において$\|\boldsymbol{n}_{k}\|=0$となったら，それは$\boldsymbol{a}_{k}$が$\boldsymbol{e}_{1},\cdots,\boldsymbol{e}_{k-1}$に従属するということですので，$r$をインクリメントせず，$\boldsymbol{a}_{k}^{\mathrm{T}}$に対応する単位直交ベクトルのインデックスもインクリメントしない，という処理が必要になります．

これらを考慮したZMにおける実装を示します．

/* LQ decomposition based on Gram=Schmidt's method. (destructive) */
int zMatDecompLQ_GramSchmidt_DST(zMat m, zMat l, zMat q)
{
  int i, j, rank;
  double *mp, r;

  zMatZero( l );
  for( rank=0, i=0; i<zMatRowSizeNC(m); i++ ){
    mp = zMatRowBuf(m,i);
    for( j=0; j<rank; j++ ){
      r = zRawVecInnerProd( mp, zMatRowBuf(q,j), zMatColSizeNC(m) );
      zRawVecCatDRC( mp, -r, zMatRowBuf(q,j), zMatColSizeNC(q) );
      zMatSetElemNC( l, i, j, r );
    }
    if( zIsTol( ( r = zRawVecNorm(mp,zMatColSizeNC(m)) ), ZM_LQ_DECOMP_GRAMSCHMIDT_TOL ) ) continue;
    zRawVecDiv( mp, r, zMatRowBuf(q,rank), zMatColSizeNC(m) );
    zMatSetElemNC( l, i, rank, r );
    if( rank < zMatColSizeNC(m) ) rank++;
  }
  return rank;
}

/* LQ decomposition based on Gram=Schmidt's method. */
int zMatDecompLQ_GramSchmidt(const zMat m, zMat l, zMat q)
{
  zMat mcp;
  int rank;

  if( !( mcp = zMatClone( m ) ) ){
    ZALLOCERROR();
    return -1;
  }
  rank = zMatDecompLQDST( mcp, l, q );
  zMatFree( mcp );
  return rank;
}

Householder法によるLQ分解

Householder変換

Householder変換とは，あるベクトルを，ある超平面（鏡面）を挟んで鏡像対称なベクトル（鏡像ベクトル）に写像する操作です．

元のベクトルを$\boldsymbol{a}$，鏡面の単位法線ベクトルを$\boldsymbol{\nu}$とおきましょう．鏡面は原点を通るものとします．鏡面上に$\boldsymbol{a}$を射影したベクトル$\boldsymbol{m}$は

\boldsymbol{m}=\boldsymbol{a}-\left(\boldsymbol{a}^{\mathrm{T}}\boldsymbol{\nu}\right)\boldsymbol{\nu}

と求まります．
$\boldsymbol{a}$から$\boldsymbol{m}$に向かうベクトルを2倍すれば$\boldsymbol{a}$と鏡像対称な位置$\boldsymbol{a}^{\prime}$に到達するので，

\boldsymbol{a}^{\prime}=
\boldsymbol{a}-2\left(\boldsymbol{a}^{\mathrm{T}}\boldsymbol{\nu}\right)\boldsymbol{\nu}
=\left(\boldsymbol{1}-2\boldsymbol{\nu}\boldsymbol{\nu}^{\mathrm{T}}\right)\boldsymbol{a}

を得ます．よって，原点を通り$\boldsymbol{\nu}$を単位法線ベクトルとする鏡面によって任意のベクトルを鏡像変換する行列は

\boldsymbol{P}=\boldsymbol{1}-2\boldsymbol{\nu}\boldsymbol{\nu}^{\mathrm{T}}

ということになります．これが正規直交行列であることは，

\boldsymbol{P}^{\mathrm{T}}\boldsymbol{P}
=
\left(\boldsymbol{1}-2\boldsymbol{\nu}\boldsymbol{\nu}^{\mathrm{T}}\right)^{\mathrm{T}}
\left(\boldsymbol{1}-2\boldsymbol{\nu}\boldsymbol{\nu}^{\mathrm{T}}\right)
=
\boldsymbol{1}-2\boldsymbol{\nu}\boldsymbol{\nu}^{\mathrm{T}}-2\boldsymbol{\nu}\boldsymbol{\nu}^{\mathrm{T}}+4\boldsymbol{\nu}\boldsymbol{\nu}^{\mathrm{T}}\boldsymbol{\nu}\boldsymbol{\nu}^{\mathrm{T}}
=\boldsymbol{1}

となることから確認できます（$\boldsymbol{\nu}^{\mathrm{T}}\boldsymbol{\nu}=1$であることを用いました）．

さらに都合が良いことに，

\boldsymbol{P}^{\mathrm{T}}
=\left(\boldsymbol{1}-2\boldsymbol{\nu}\boldsymbol{\nu}^{\mathrm{T}}\right)^{\mathrm{T}}
=\boldsymbol{1}-2\boldsymbol{\nu}\boldsymbol{\nu}^{\mathrm{T}}
=\boldsymbol{P}

ですので，この$\boldsymbol{P}$は対称行列でもあります．これは，「鏡像の鏡像は自分自身である」ということを考えても明らかです．

逆に，ベクトル$\boldsymbol{a}$を別のベクトル$\boldsymbol{a}^{\prime}$に鏡像変換するような鏡面の法線ベクトル$\boldsymbol{\nu}$を見つけることを考えましょう．ただし，$\|\boldsymbol{a}^{\prime}\|=\|\boldsymbol{a}\|$である必要があります．

\boldsymbol{a}^{\prime}=\boldsymbol{a}-2\left(\boldsymbol{a}^{\mathrm{T}}\boldsymbol{\nu}\right)\boldsymbol{\nu}
\qquad\Leftrightarrow\qquad
2\left(\boldsymbol{a}^{\mathrm{T}}\boldsymbol{\nu}\right)\boldsymbol{\nu}=\boldsymbol{a}-\boldsymbol{a}^{\prime}

ですから，$\boldsymbol{\nu}$は$\boldsymbol{a}-\boldsymbol{a}^{\prime}$に平行であることが分かります．$\boldsymbol{\nu}$は単位ベクトルなので，

\boldsymbol{\nu}=\frac{\boldsymbol{a}-\boldsymbol{a}^{\prime}}{\|\boldsymbol{a}-\boldsymbol{a}^{\prime}\|}

となります（$\boldsymbol{a}^{\prime}=\boldsymbol{a}$である場合は，変換する必要が無いので除外して考えて構いません）．

これを応用して，Householder変換によってベクトル$\boldsymbol{a}$をある単位ベクトル$\boldsymbol{e}$で張られる直線上に写像する鏡面の法線ベクトル$\boldsymbol{\nu}$を，簡単に見つけることができます．

鏡像ベクトル$\boldsymbol{a}^{\prime}$は$\boldsymbol{a}$とノルムが等しくなりますから，

\boldsymbol{a}^{\prime}=\pm\|\boldsymbol{a}\|\boldsymbol{e}

として上記の計算を適用すれば良いのです．
複号となっているのは，同一直線上に鏡像の候補が二つあるためですが，桁落ちを避けるために，$\boldsymbol{e}^{\mathrm{T}}\boldsymbol{a}<0$であるときは$\boldsymbol{e}$と同じ向きの，$\boldsymbol{e}^{\mathrm{T}}\boldsymbol{a}\geq 0$であるときは$\boldsymbol{e}$と逆向きのベクトルをそれぞれ選んで，$\boldsymbol{a}^{\prime}$が$\boldsymbol{a}$からなるべく離れるようにするのが好ましいです．よって，

\boldsymbol{a}^{\prime}=-\mathrm{sgn}\left(\boldsymbol{e}^{\mathrm{T}}\boldsymbol{a}\right)\|\boldsymbol{a}\|\boldsymbol{e}

とすれば良いでしょう．

Householder変換により行列を下三角行列化する

行列$\boldsymbol{A}$を，右から変換行列をかける（のと等価な操作をする）ことで下三角行列化する，というのがHouseholder法によるLQ分解の基本的な考え方です．変換は$\boldsymbol{A}$の列ベクトルではなく行ベクトルに対してなされることになります．$\boldsymbol{A}$のある行ベクトル$\boldsymbol{a}^{\mathrm{T}}$を別の行ベクトル$\boldsymbol{a}^{\prime\mathrm{T}}$に鏡像変換する操作は

\boldsymbol{a}^{\prime\mathrm{T}}
=\boldsymbol{a}^{\mathrm{T}}-2\left(\boldsymbol{a}^{\mathrm{T}}\boldsymbol{\nu}\right)\boldsymbol{\nu}^{\mathrm{T}}
=\boldsymbol{a}^{\mathrm{T}}\left(\boldsymbol{1}-2\boldsymbol{\nu}\boldsymbol{\nu}^{\mathrm{T}}\right)
(=\boldsymbol{a}^{\mathrm{T}}\boldsymbol{P})

となります．写像行列$\boldsymbol{P}$が対称であるため，左右が入れ替わっただけになっています．

これに基づいて，まず$\boldsymbol{A}$の1行目ベクトル$\boldsymbol{a}_{1}^{\mathrm{T}}=[a_{1,1}~\cdots~a_{1,n}]$を$\boldsymbol{e}_{1}^{\mathrm{T}}=[1~0~\cdots~0]$で張られる直線上に写像することを考えます．写像先のベクトルは

\boldsymbol{a}_{1}^{\prime\mathrm{T}}=[s_{1}~0~\cdots~0]
\qquad(s_{1}=-\mathrm{sgn}(a_{1,1})\|\boldsymbol{a}_{1}\|)

となり，鏡面の法線ベクトル$\boldsymbol{\nu}_{1}$は

\boldsymbol{u}_{1}=\begin{bmatrix}
a_{1,1}-s_{1} & a_{1,2} & \cdots & a_{1,n}
\end{bmatrix}
,\qquad
\boldsymbol{\nu}_{1}=\boldsymbol{u}_{1}/\|\boldsymbol{u}_{1}\|

と求まります．また，残りの行は

\boldsymbol{a}_{i}^{\prime\mathrm{T}}
=\boldsymbol{a}_{i}^{\mathrm{T}}-2\left(\boldsymbol{a}_{i}^{\mathrm{T}}\boldsymbol{\nu}_{1}\right)\boldsymbol{\nu}_{1}^{\mathrm{T}}
\qquad\mbox{($i=2,\cdots,m$)}

と変換されます．ただし，$\boldsymbol{a}_{i}^{\mathrm{T}}$（$i=2,\cdots,m$）は$\boldsymbol{A}$の$i$行目ベクトルです．

変換後の行列を$\boldsymbol{A}_{1}(=[\boldsymbol{a}_{1}^{\prime}~\cdots~\boldsymbol{a}_{n}^{\prime}]^{\mathrm{T}})$とおきましょう．形式的に$\boldsymbol{P}_{1}=\boldsymbol{1}-2\boldsymbol{\nu}_{1}\boldsymbol{\nu}_{1}^{\mathrm{T}}$とおくと，$\boldsymbol{A}\boldsymbol{P}_{1}=\boldsymbol{A}_{1}$かつ$\boldsymbol{P}_{1}=\boldsymbol{P}_{1}^{\mathrm{T}}$ですから，これは

\boldsymbol{A}=\boldsymbol{A}_{1}\boldsymbol{P}_{1}

と分解したのと同じことです．

以降の操作は再帰的に考えます．$\boldsymbol{A}$が

\displaylines{
\boldsymbol{A}=\boldsymbol{A}_{k}\boldsymbol{Q}_{k}
\\
\boldsymbol{A}_{k}=\left[\begin{array}{c|c}
\boldsymbol{L}_{k} & \boldsymbol{O} \\
\hline
* & \tilde{\boldsymbol{A}}_{k}
\end{array}\right]
}

という形になっているとしましょう．ただし，$\boldsymbol{Q}_{k}$は$n\times n$正規直交行列，$\boldsymbol{L}_{k}$は$k\times k$下三角行列，$\tilde{\boldsymbol{A}}_{k}$は$(m-k)\times (n-k)$行列です．$\boldsymbol{A}_{k}$の左下は計算に影響しないので，「気にしない」の意味で$*$としています．

先ほどと全く同じ考え方を適用することで，$\tilde{\boldsymbol{A}}_{k}$の1行目ベクトル$\tilde{\boldsymbol{a}}_{1}^{\mathrm{T}}=[\tilde{a}_{1,1}~\cdots~\tilde{a}_{1,n-k}]$を
$\tilde{\boldsymbol{e}}_{1}^{\mathrm{T}}=[1~0~\cdots~0]$で張られる直線上に写像することができます．すなわち，写像先のベクトルは

\tilde{\boldsymbol{a}}_{1}^{\prime\mathrm{T}}=[\tilde{s}_{1}~0~\cdots~0]
\qquad(\tilde{s}_{1}=-\mathrm{sgn}(\tilde{a}_{1,1})\|\tilde{\boldsymbol{a}}_{1}\|)

となり，鏡面の法線ベクトル$\tilde{\boldsymbol{\nu}}_{k}$は

\tilde{\boldsymbol{u}}_{k+1}=\begin{bmatrix}
\tilde{a}_{1,1}-\tilde{s}_{1} & \tilde{a}_{1,2} & \cdots & \tilde{a}_{1,n-k}
\end{bmatrix}
,\qquad
\tilde{\boldsymbol{\nu}}_{k+1}=\tilde{\boldsymbol{u}}_{k+1}/\|\tilde{\boldsymbol{u}}_{k+1}\|

と求まります．また，残りの行は

\tilde{\boldsymbol{a}}_{i}^{\prime\mathrm{T}}
=\tilde{\boldsymbol{a}}_{i}^{\mathrm{T}}-2\left(\tilde{\boldsymbol{a}}_{i}^{\mathrm{T}}\tilde{\boldsymbol{\nu}}_{k+1}\right)\tilde{\boldsymbol{\nu}}_{k+1}^{\mathrm{T}}
\qquad\mbox{($i=2,\cdots,m-k$)}

と変換されます．ただし，$\tilde{\boldsymbol{a}}_{i}^{\mathrm{T}}$（$i=2,\cdots,m-k$）は$\tilde{\boldsymbol{A}}_{k}$の$i$行目ベクトルです．

これもまた先ほどと同様に，変換後の行列を$\tilde{\boldsymbol{A}}_{k+1}$とおきましょう．形式的に$\tilde{\boldsymbol{P}}_{k+1}=\boldsymbol{1}-2\tilde{\boldsymbol{\nu}}_{k+1}\tilde{\boldsymbol{\nu}}_{k+1}^{\mathrm{T}}$とおくと，上記の操作は

\tilde{\boldsymbol{A}}_{k}=\tilde{\boldsymbol{A}}_{k+1}\tilde{\boldsymbol{P}}_{k+1}

と分解したのと同じことです．そして，この操作は元の$\boldsymbol{A}_{k}$の$1\sim k$行目と$1\sim k$列目に全く影響を与えないので，$\boldsymbol{\nu}_{k+1}=[0~\cdots~0~\tilde{\boldsymbol{\nu}}_{k+1}^{\mathrm{T}}]^{\mathrm{T}}$，$\boldsymbol{P}_{k+1}=\boldsymbol{1}-2\boldsymbol{\nu}_{k+1}\boldsymbol{\nu}_{k+1}^{\mathrm{T}}$とおくと，

\boldsymbol{A}_{k}=\boldsymbol{A}_{k+1}\boldsymbol{P}_{k+1}

と分解したのと同じこととも言えます．

これを$\boldsymbol{A}=\boldsymbol{A}_{k}\boldsymbol{Q}_{k}$に代入すれば，

\boldsymbol{A}=\boldsymbol{A}_{k+1}\boldsymbol{P}_{k+1}\boldsymbol{Q}_{k}

となります．形式的に

\boldsymbol{A}=\boldsymbol{A}_{k+1}\boldsymbol{Q}_{k+1}

と見なせば，変換行列$\boldsymbol{Q}_{k}$も

\boldsymbol{Q}_{k+1}=\boldsymbol{P}_{k+1}\boldsymbol{Q}_{k}

という再帰的構造を持つことが分かります．$\boldsymbol{P}_{k+1}$と$\boldsymbol{Q}_{k}$はともに正規直交行列なので，それらの積である$\boldsymbol{Q}_{k+1}$もまた正規直交行列となります．

ここで，$\boldsymbol{Q}_{k}$の下$(n-k)\times n$部分行列を$\tilde{\boldsymbol{Q}}_{k}$とおきましょう．

\boldsymbol{Q}_{k}=\left[\begin{array}{c}
* \\
\hline
\tilde{\boldsymbol{Q}}_{k}
\end{array}\right]

と表せるということです．$\boldsymbol{\nu}_{k+1}$の$1\sim k$番目成分が0であることより，

\boldsymbol{Q}_{k+1}=\left[\begin{array}{c}
* \\
\hline
(\boldsymbol{1}-2\tilde{\boldsymbol{\nu}}_{k+1}\tilde{\boldsymbol{\nu}}_{k+1}^{\mathrm{T}})\tilde{\boldsymbol{Q}}_{k}
\end{array}\right]

が成り立ちます．したがって，$\tilde{\boldsymbol{Q}}_{k}$の$i$列目ベクトルを$\tilde{\boldsymbol{q}}_{i}$（$i=1,\cdots,n$），$\tilde{\boldsymbol{Q}}_{k+1}$の$i$列目ベクトルを$\tilde{\boldsymbol{q}}_{i}^{\prime}$（$i=1,\cdots,n$）とそれぞれおけば，

\tilde{\boldsymbol{q}}_{i}^{\prime}=\tilde{\boldsymbol{q}}_{i}-2(\tilde{\boldsymbol{q}}_{i}^{\mathrm{T}}\tilde{\boldsymbol{\nu}}_{k+1})\tilde{\boldsymbol{\nu}}_{k+1}
\qquad\mbox{($i=1,\cdots,n$)}

となります．

なお，もう少し計算を工夫することもできます．

\|\tilde{\boldsymbol{u}}_{k+1}\|^{2}=\|\tilde{\boldsymbol{a}}_{1}-\tilde{s}_{1}\tilde{\boldsymbol{e}}_{1}\|^{2}
=\|\tilde{\boldsymbol{a}}_{1}\|^{2}-2\tilde{s}_{1}\tilde{\boldsymbol{a}}_{1}^{\mathrm{T}}\tilde{\boldsymbol{e}}_{1}+\tilde{s}_{1}^{2}\|\tilde{\boldsymbol{e}}_{1}\|^{2}
=2\tilde{s}_{1}(\tilde{s}_{1}-\tilde{a}_{1,1})

ですから，

\beta=1/\tilde{s}_{1}(\tilde{s}_{1}-\tilde{a}_{1,1})

とおけば，

\displaylines{
\tilde{\boldsymbol{a}}_{i}^{\prime\mathrm{T}}
=\tilde{\boldsymbol{a}}_{i}^{\mathrm{T}}-\beta\left(\tilde{\boldsymbol{a}}_{i}^{\mathrm{T}}\tilde{\boldsymbol{u}}_{k+1}\right)\tilde{\boldsymbol{u}}_{k+1}^{\mathrm{T}}
\qquad\mbox{($i=k+1,\cdots,m$)}
\\
\tilde{\boldsymbol{q}}_{i}^{\prime}
=\tilde{\boldsymbol{q}}_{i}-\beta\left(\tilde{\boldsymbol{q}}_{i}^{\mathrm{T}}\tilde{\boldsymbol{u}}_{k+1}\right)\tilde{\boldsymbol{u}}_{k+1}
\qquad\mbox{($i=1,\cdots,n$)}
}

となって，$\tilde{\boldsymbol{u}}_{k+1}$を正規化して$\tilde{\boldsymbol{\nu}}_{k+1}$を求める計算が無くなり，計算効率とメモリの使用効率を若干向上できます．$\tilde{s}_{1}$が$\tilde{a}_{1,1}$と異符号であるので，$\tilde{s}_{1}\neq 0$ならば$\beta$の分母は非零であることに注意して下さい．

実装

前節の内容に基づけば，

\boldsymbol{A}_{0}=\boldsymbol{A}
,\qquad
\boldsymbol{Q}_{0}=\boldsymbol{1}

とおき，$k=1,\cdots,l$（$l$は$m$，$n$のうち小さい方の値）について上記の計算を反復して$\boldsymbol{A}_{k}$と$\boldsymbol{Q}_{k}$を順次更新することで，最終的に下三角行列$\boldsymbol{L}=\boldsymbol{A}_{l}$と，変換行列$\boldsymbol{Q}=\boldsymbol{Q}_{l}$が得られます．

$\boldsymbol{Q}$は正方行列になりますので，Householder法により行えるLQ分解はGram=Schmidtの直交化によるそれと異なり，

\boldsymbol{A}=\boldsymbol{L}\boldsymbol{Q}

の形をとることになります．もちろん，$\boldsymbol{L}$と$\boldsymbol{Q}$から$\boldsymbol{L}_{r}$と$\boldsymbol{Q}_{r}$を取り出すことは可能です．

なお，反復計算の途中で$\tilde{s}_{1}=0$となった場合，そのときの行ベクトル$\boldsymbol{a}_{k}^{\mathrm{T}}$は，それ以前に現れた$\boldsymbol{a}_{1}^{\mathrm{T}},\cdots,\boldsymbol{a}_{k-1}^{\mathrm{T}}$に従属するということですので，Gram=Schmidtの直交化による方法と同様に，$\boldsymbol{A}$のランク$r$をインクリメントしないという処理が必要になります．

これらを考慮したZMにおける実装を示します．

/* LQ decomposition based on Householder method. (destructive) */
static int _zMatDecompLQ_Householder_DST(zMat m, zMat q)
{
  double s, ds, norm_inv, reflection;
  double *u;
  int i, j, size, colsize, rank;

  zMatIdentNC( q );
  size = zMatMinSize( m );
  for( rank=0, i=0; i<size; i++ ){
    u = &zMatElemNC(m,i,i);
    colsize = zMatColSizeNC(m) - i;
    if( zIsTiny( ( s = -zSgn( zMatElemNC(m,i,i) ) * zRawVecNorm( u, colsize ) ) ) ||
        zIsTiny( ( ds = s - zMatElemNC(m,i,i) ) ) ) continue;
    norm_inv = 1.0 / ( s * ds );
    *u = -ds; /* to use u temporarily as a reflection vector */
    for( j=0; j<zMatColSizeNC(q); j++ ){
      reflection = -norm_inv * zRawMatColInnerProd( zMatRowBufNC(q,i), u, zMatRowSizeNC(q)-i, zMatColSizeNC(q), j );
      zRawMatColCatDRC( zMatRowBufNC(q,i), reflection, u, zMatRowSizeNC(q)-i, zMatColSizeNC(q), j );
    }
    for( j=zMatRowSizeNC(m)-1; j>i; j-- ){
      reflection = -norm_inv * zRawVecInnerProd( &zMatElemNC(m,j,i), u, colsize );
      zRawVecCatDRC( &zMatElemNC(m,j,i), reflection, u, colsize );
    }
    *u = s;
    zRawVecZero( u + 1, colsize - 1 );
    rank++;
  }
  return rank;
}

/* LQ decomposition based on Householder method. */
int zMatDecompLQ_Householder(const zMat m, zMat l, zMat q)
{
  zMat ltmp, qtmp;
  int rank = -1;

  ltmp = zMatClone( m );
  qtmp = zMatAllocSqr( zMatColSizeNC(m) );
  if( !ltmp || !qtmp ){
    ZALLOCERROR();
    goto TERMINATE;
  }
  if( ( rank = _zMatDecompLQ_Householder_DST( ltmp, qtmp ) ) <= 0 ){
    ZRUNERROR( ZM_ERR_MAT_CANNOTDECOMPOSEZEROMAT );
    goto TERMINATE;
  }
  if( l ){
    zMatSetColSizeNC( l, rank );
    zMatGetNC( ltmp, 0, 0, l );
  }
  if( q ){
    zMatSetRowSizeNC( q, rank );
    zMatGetNC( qtmp, 0, 0, q );
  }

 TERMINATE:
  zMatFreeAtOnce( 2, ltmp, qtmp );
  return rank;
}

どちらを使うべきか

二つの方法の特徴をまとめます．

どちらの方法も，ベクトルの内積，線形結合，正規化という基本的な演算の組み合わせですし，反復方法もシンプルなので，簡単に実装することができます．

その上で，Gram=Schmidtの直交化による行列のLQ分解は，1行ずつ順番に単位直交ベクトルを求めていくことができるので，高速です．一方，与えられた行列の行ベクトルの中で，「線形独立ではあるけれど平行にとても近い組み合わせ」があったときに，数値的に不安定になるという難点があります．また，結果が階数分解表現でしか得られないのは制約になり得ます．

Householder法による行列のLQ分解は，全体を少しずつ変形していって最終的に全ての単位直交ベクトルの組を得る方法です．このため，高速さという点ではGram=Schmidtの直交化による方法よりも若干劣ります．一方，写像先のベクトルの方向を適切に選ぶことで数値的に不安定になることを避けられるのは安心感があります．正規直交行列が正方行列として得られるため，像空間写像行列と核空間写像行列が同時に求まるということも，有利に働く場面は多いでしょう．

実用に当たっては，これらのことを理解した上で適切に使い分けることになります．