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

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

[←]　[↑]　[→]

QR分解のアルゴリズム

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

Aの要素は実数とする。

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

/* AをRに入れる */

for j=1 to n

for i=1 to m

r_ij = a_ij

/* RとQ^†bの計算 */

for j=1 to n

for i=j+1 to m

d = sqrt(r_jj² + r_ij²)

c = r_jj/d , s = r_ij/d

r_jj = d , r_ij= 0

for k=j+1 to n

tmp = c r_jk + s r_ik

r_jk = -s r_jk + c r_ik

r_ik = tmp

tmp = c b_j + s b_i

b_j = -s b_j + c b_i

b_i = tmp

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

/* AをRに入れる */

for j=1 to n

for i=1 to m

r_ij = a_ij

/* RとQ^†bの計算 */

for j=1 to n

for i=m-1 to j

d = sqrt(r_ij² + r_i+1j²)

c = r_ij/d , s = r_i+1j/d

r_ij = d , r_i+1j= 0

for k=j+1 to n

tmp = c r_ik + s r_i+1k

r_i+1k = -s r_ik + c r_i+1k

r_ik = tmp

tmp = c b_i + s b_i+1

b_i+1 = -s b_i + c b_i+1

b_i = tmp

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

/* AをRに入れる */

for j=1 to n

for i=1 to m

r_ij = a_ij

/* RとQ^†bの計算 */

for j=1 to n

for i=m-1 to j

r_ij = d , r_i+1j= 0

for k=j+1 to n

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

/* Rx= Q^†bをxについて解く(後退代入) */

x_n = b_n / r_nn

for i=n-1 to 1

tmp = 0

for j=i+1 to n

tmp = tmp + r_ijx_j

x_i = (b_i - tmp)/ r_ii

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

◇

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

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

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

i行j列目の要素をゼロにするには

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

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

G_N … G₁A = R

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

A = QR , Q = G₁^T …G_N^T

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

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

Rx = G_N … G₁ b

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

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

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

・j列(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回転、トップダウン順】

/* AをRに入れる */

for j=1 to n

for i=1 to m

r_ij = a_ij

/* RとQ^†bの計算 */

for j=1 to n

for i=j+1 to m

d = sqrt(r_jj² + r_ij²) , c = r_jj/d , s = r_ij/d

r_jj = d , r_ij= 0

for k=j+1 to n

tmp = c r_jk + s r_ik

r_jk = -s r_jk + c r_ik

r_ik = tmp

tmp = c b_j + s b_i

b_j = -s b_j + c b_i

b_i = tmp

上記はトップダウン順(i=j+1 to m)なので離れた行どうし

の演算になっているが、ボトムアップ順(i=m-1 to j)に

すれば隣どうしの行の演算になる

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

/* RとQ^†bの計算 */

for j=1 to n

for i=m-1 to j

d = sqrt(r_ij² + r_i+1j²) , c = r_ij/d , s = r_i+1j/d

r_ij = d , r_i+1j= 0

for k=j+1 to n

tmp = c r_ik + s r_i+1k

r_i+1k = -s r_ik + c r_i+1k

r_ik = tmp

tmp = c b_i + s b_i+1

b_i+1 = -s b_i + c b_i+1

b_i = tmp

◇◇

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

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

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

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

j列の回転結果は

j列の右側を回転する

bを回転する

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

/* RとQ^†bの計算 */

for j=1 to n

for i=m-1 to j

　

r_ij = d , r_i+1j= 0

for k=j+1 to n

　　　　

◇◇◇

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

[←]　[↑]　[→]