QR分解(Householder変換による)

[] [] []


QR分解のアルゴリズム

A=QR ただしQ自身でなくQbを求める

Householder変換】

/* ARに入れる */

for j=1 to n

    for i=1 to m

        rij = aij

/* RQbの計算 */

for j=1 to n

    tmp = 0

    for i=j+1 to m

        tmp = tmp + |rij|2

    vj = rjj + sgn(rjj) sqrt(|rjj|2 + tmp)

    tmp = sqrt(|vj|2 + tmp)

    for i=j+1 to m

        vj = rjj /tmp

    for k=j to n

        tmp = 0

       for i=j to m

            tmp = tmp + vi* rik

       for i=j to m

            rik = rik - 2 tmp vi

    tmp = 0

   for i= j to m

        tmp = tmp + vi* bi

   for i=j to m

        bi = bi - 2 tmp vi

連立1次方程式Ax=bを解くには

/* Rx= Qbxについて解く(後退代入) */

xn = bn / rnn

for i=n-1 to 1

   tmp = 0

   for j=i+1 to n

      tmp = tmp + rij xj

   xi = (bi - tmp)/ rii

前項ではGivens回転を使って三角化したが、鏡映変換(ハウスホルダー変換)を使っても三角化できる

法線v の平面に対する鏡映変換Hは次のように表せる。

∵ベクトルa から平面におろした垂線を2倍に伸ばしてa に足せば、a の鏡映になる

このHA(mn列、m≧n)にかけて

11列目の要素a11の下のすべての要素a1i(j=2 to m)をゼロにするには

と、a1e方向に鏡映変換するように平面の法線vを選べばよいから

となる。v 2とおりあるが、なるべく桁落ちしないようにv が大きい方を選んで(v 2乗するので全体の符号は変えてもよい)

このH1Aにかければ、対角要素a11の下のすべての要素a1i(j=2 to m)をゼロにできる。

さらに、この手順をすべての対角要素ajj (j=1,...,n)について繰り返せばAを上三角行列Rにできる。

Hn … H1A = R

Hjはユニタリ行列かつエルミート行列なので、Hの逆行列はHそのものだから

A = QR , Q = (Hn … H1) -1 = (H1 -1) … (Hn -1)= H1 … Hn

QR分解できたことになる。

なおAx=bを解くだけならばQを陽に求める必要はなく、bHをかけていけばよい。

Rx = Hn … H1 b

以上の手順を整理すると次のようになる。

手順

コード

Aを作業領域Rに入れる。

j(=1 to n)のステップにおいて

対角要素rjjの下すべてrij (i=j+1 to m)をゼロにする行列Hjを求める

j列目の列ベクトルに対して

vを次のように決めて

RHをかける

bにも同じHをかける

Householder変換】

/* ARに入れる */

for j=1 to n

    for i=1 to m

        rij = aij

/* RQbの計算 */

for j=1 to n

    tmp = 0

    for i=j+1 to m

        tmp = tmp + |rij|2

    vj = rjj + sgn(rjj) sqrt(|rjj|2 + tmp)

    tmp = sqrt(|vj|2 + tmp)

    for i=j+1 to m

        vj = rjj /tmp

    for k=j to n

        tmp = 0

       for i=j to m

            tmp = tmp + vi* rik

       for i=j to m

            rik = rik - 2 tmp vi

    tmp = 0

   for i= j to m

        tmp = tmp + vi* bi

   for i=j to m

        bi = bi - 2 tmp vi


[] [] []