QR分解(Givens回転による、実数の場合)

[] [] []


QR分解のアルゴリズム

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

Aの要素は実数とする。

Givens回転、トップダウン順】

/* ARに入れる */

for j=1 to n

    for i=1 to m

        rij = aij

/* RQbの計算 */

for j=1 to n

    for i=j+1 to m

        d = sqrt(rjj2 + rij2)

        c = rjj/d , s = rij/d

        rjj = d , rij = 0

        for k=j+1 to n

            tmp = c rjk + s rik

            rjk = -s rjk + c rik

            rik = tmp

        tmp = c bj + s bi

        bj = -s bj + c bi

        bi = tmp

Givens回転、ボトムアップ順】

/* ARに入れる */

for j=1 to n

    for i=1 to m

        rij = aij

/* RQbの計算 */

for j=1 to n

    for i=m-1 to j

        d = sqrt(rij2 + ri+1j2)

        c = rij/d , s = ri+1j/d

        rij = d , ri+1j = 0

        for k=j+1 to n

            tmp = c rik + s ri+1k

            ri+1k = -s rik + c ri+1k

            rik = tmp

        tmp = c bi + s bi+1

        bi+1 = -s bi + c bi+1

        bi = tmp

Givens回転、ボトムアップ順、CORDIC

/* ARに入れる */

for j=1 to n

    for i=1 to m

        rij = aij

/* RQbの計算 */

for j=1 to n

    for i=m-1 to j

    

        rij = d , ri+1j = 0

        for k=j+1 to n

          

    

連立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

前項のGram-Schmidt法では、三角行列Rをかけて列ベクトルを直交化していったので直交性に誤差が生じてしまう。そこで逆に、列ベクトルがちゃんと正規直交しているQをかけて三角化すればよく、こちらの方が数値計算の安定性がよいらしい。直交行列と言えば回転か鏡映だから、ここでは2次元の回転(ギブンス回転)を使ったアルゴリズムを導く

本項では簡単のため回転行列が実数で済むように行列Aの要素は実数とし、次項で複素数に拡張する

2次元の回転行列(mm列のうちj,ij,i列の4要素以外は単位行列になっている)を

とおき、これをA(mn列、m≧n)にかけて

ij列目の要素をゼロにするには

s,cを決めればよい。これをi(i=j+1,...,m)について繰り返せばajjの下をすべてゼロにできる。

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

GN … G1A = R

Gの逆行列はGを転置したものだから

A = QR , Q = GNT … G1T

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

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

Rx = GN … G1 b

ここでGをかける回数Nは、2次元の回転でAの要素を1個づつゼロにしているので次のとおり。

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

手順

コード

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

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

i行目(i=j+1 to m)をゼロにする行列G(ij)を求める

G(ij)は次の4要素以外は単位行列であるので

G(ij)Rにかけるとi行目とj行目以外は変わらない。

またj列の左側は既にゼロになっており、j列の結果は

とわかっているので、j列の右側だけ計算すればよい

bにも同じG(ij)をかける

Givens回転、トップダウン順】

/* ARに入れる */

for j=1 to n

    for i=1 to m

        rij = aij

/* RQbの計算 */

for j=1 to n

    for i=j+1 to m

        d = sqrt(rjj2 + rij2) , c = rjj/d , s = rij/d

        rjj = d , rij = 0

        for k=j+1 to n

            tmp = c rjk + s rik

            rjk = -s rjk + c rik

            rik = tmp

        tmp = c bj + s bi

        bj = -s bj + c bi

        bi = tmp

上記はトップダウン順(i=j+1 to m)なので離れた行どうしの演算になっているが、ボトムアップ順(i=m-1 to j)にすれば隣どうしの行の演算になる

Givens回転、ボトムアップ順】

/* RQbの計算 */

for j=1 to n

    for i=m-1 to j

        d = sqrt(rij2 + ri+1j2) , c = rij/d , s = ri+1j/d

         rij = d , ri+1j = 0

         for k=j+1 to n

             tmp = c rik + s ri+1k

             ri+1k = -s rik + c ri+1k

             rik = tmp

        tmp = c bi + s bi+1

        bi+1 = -s bi + c bi+1

        bi = tmp

◇◇

回転行列の要素を求めるとき、平方根とわり算が必要になるが、もし固定小数点でよければ加減算だけで計算するCORDICという方法があるCORDICには、半径と偏角を求めるVectorモードと、回転結果を求めるRotationモードがある。

これらを使えば平方根とわり算を使わずに次のように計算できる。

手順

コード

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

i行目(i=m-1 to j)をゼロにする行列G(ij)のパラメータを求める

j列の回転結果は

j列の右側を回転する

bを回転する

Givens回転、ボトムアップ順】

/* RQbの計算 */

for j=1 to n

    for i=m-1 to j

    

        rij = d , ri+1j = 0

        for k=j+1 to n

        

    

◇◇◇

Givens回転によるアルゴリズムは、ソフトウェアでは実装は簡単だが計算が重くなるのであまり使われず、CORDICを実装しやすいFPGAで使われるようである。


[] [] []