QR分解(Householder変換による)

[←]　[↑]　[→]

QR分解のアルゴリズム

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

【Householder変換】

/* AをRに入れる */

for j=1 to n

for i=1 to m

r_ij = a_ij

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

for j=1 to n

tmp = 0

for i=j+1 to m

tmp = tmp + |r_ij|²

v_j= r_jj + sgn(r_jj) sqrt(|r_jj|² + tmp)

tmp = sqrt(|v_j|² + tmp)

v_j= v_j /tmp

for i=j+1 to m

v_i= r_ij /tmp

for k=j to n

tmp = 0

for i=j to m

tmp = tmp + v_i^*r_ik

for i=j to m

r_ik = r_ik - 2 tmp v_i

tmp = 0

for i= j to m

tmp = tmp + v_i^* b_i

for i=j to m

b_i = b_i - 2 tmp v_i

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

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

◇

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

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

このHをA(m行n列、m≧n)にかけて

1行1列目の要素a₁₁の下のすべての要素a_1i(j=2 to m)をゼロにするには

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

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

このH₁をAにかければ、対角要素a₁₁の下のすべての要素a_1i(j=2 to m)をゼロにできる。

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

H_n … H₁A = R

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

A = QR , Q = (H_n … H₁)^-1 = (H₁^-1) … (H_n^-1)= H₁ … H_n

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

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

Rx = H_n … H₁ b

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

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

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

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

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

vを次のように決めて

RにHをかける

bにも同じHをかける

【Householder変換】

/* AをRに入れる */

for j=1 to n

for i=1 to m

r_ij = a_ij

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

for j=1 to n

tmp = 0

for i=j+1 to m

tmp = tmp + |r_ij|²

v_j= r_jj + sgn(r_jj) sqrt(|r_jj|² + tmp)

tmp = sqrt(|v_j|² + tmp)

v_j= v_j /tmp

for i=j+1 to m

v_i= r_ij /tmp

for k=j to n

tmp = 0

for i=j to m

tmp = tmp + v_i^*r_ik

for i=j to m

r_ik = r_ik - 2 tmp v_i

tmp = 0

for i= j to m

tmp = tmp + v_i^* b_i

for i=j to m

b_i = b_i - 2 tmp v_i

[←]　[↑]　[→]