QR分解(Gram-Schmidtの直交化法による)

[←]　[↑]　[→]

---

QR分解のアルゴリズム A=QR
【修正Gram-Schmidt , ループkj順】 for k=1 to n     /* AをQに入れる */     for i=1 to m         qik = aik     /* QとRの計算 */     for j=1 to k-1         rjk = 0         for i=1 to m             rjk = rjk + qij^* qik …(★)         for i=1 to m             qik = qik - rjkqij     rkk = 0     for i=1 to m         rkk = rkk + qik^* qik     rkk = sqrt(rkk)     for i=1 to m         qik = qik/rkk	【修正Gram-Schmidt , ループjk順】 /* AをQに入れる */ for j=1 to n     for i=1 to m         qij = aij /* QとRの計算 */ for j=1 to n     rjj = 0     for i=1 to m         rjj = rjj + qij^* qij     rjj = sqrt(rjj)     for i=1 to m         qij = qij/rjj     for k=j+1 to n         rjk = 0         for i=1 to m             rjk = rjk + qij^* qik …(★)         for i=1 to m             qik = qik - rjkqij
(注)オリジナルのGram-Schmidt法では(★)でqik の代わりにaik を使う。 rjk = rjk + qij^* aik
連立1次方程式Ax=bを解くには /* c=Q^†bを計算する*/ for j=1 to n    cj = 0    for i=1 to m       cj = cj + qij^* bi /* Rx=cをxについて解く(後退代入) */ xn = cn / rnn for i=n-1 to 1    tmp = 0    for j=i+1 to n       tmp = tmp + rij xj    xi = (ci - tmp)/ rii

一般のm行n列(m≧n)の行列Aを、列ベクトルが正規直交であるQと上三角行列Rに分解するアルゴリズムを導く。

◇

アルゴリズムを導く方針は、A=QRと分解できたとして要素ごとに式を書き下し、これをq,rについて解く。

簡単のためn=3の場合を書き下すと

a₁ = r₁₁ q₁

a₂ = r₁₂ q_{1 +} r₂₂ q₂

a₃ = r₁₃ q_{1 +} r₂₃ q_{2 +} r₃₃ q₃

これをq_jについて解くと

q₁ = a₁/r₁₁

q₂ =(a₂ - r₁₂ q₁)/ r₂₂

q₃ =(a₃ - r₁₃ q₁ - r₂₃ q₂)/ r₃₃

q_jが互いに直交してかつノルム1になるようにr_ijを決めると

r₁₁ = ||a₁||

r₁₂ = q₁^†a₂ , r₂₂ = ||a₂ - r₁₂ q₁||

r₁₃ = q₁^†a₃ , r₂₃ = q₂^†a₃ , r₃₃ = ||a₃ - r₁₃ q₁ - r₂₃ q₂||

これは、任意のベクトルa_jから正規直交ベクトルq_jを作るGram-Schmidtの直交化法と同じである。

v₁ = a₁ , q₁ = v₁/||v₁||

v₂ =(a₂ - (q₁^†a₂)q₁) , q₂ = v₂/||v₂||

v₃ =(a₃ - (q₁^†a₃)q₁ - (q₂^†a₃)q₂) , q₃ = v₃/||v₃||

この手順を一般化すると次のようになる。

・ステップk (k=1 to n)において q1からqk-1と直交するqkを求める for k=1 to n vk = ak     for j=1 to k-1         rjk = qj^†ak        vk = vk - rjk qj rkk = ||vk|| qk = vk /rkk 変数を節約するためvの代わりにqを使うとよい。

【古典Gram-Schmidt , ループkj順】 for k=1 to n     /* AをQに入れる */     for i=1 to m         qik = aik     /* QとRの計算 */     for j=1 to k-1         rjk = 0         for i=1 to m             rjk = rjk + qij^*aik         for i=1 to m             qik = qik - rjkqij     rkk = 0     for i=1 to m         rkk = rkk + qik^* qik     rkk = sqrt(rkk)     for i=1 to m         qik = qik/rkk

2重ループの順序を次のように入れ替えてもよい。

・ステップk (k=1 to n)において q1からqk-1と直交するqkを求める for j=1 to n vj = aj for j=1 to n rjj = ||vj|| qj = vj /rjj     for k=j+1 to n         rjk = qj^†ak         vk = vk - rjk qj 変数を節約するためvの代わりにqを使うとよい。

【古典Gram-Schmidt , ループjk順】 /* AをQに入れる */ for j=1 to n     for i=1 to m         qij = aij /* QとRの計算 */ for j=1 to n     rjj = 0     for i=1 to m         rjj = rjj + qij^* qij     rjj = sqrt(rjj)     for i=1 to m         qij = qij/rjj     for k=j+1 to n         rjk = 0         for i=1 to m             rjk = rjk + qij^* aik         for i=1 to m             qik = qik - rjkqij

◇◇

上記のGram-Schmidt法では、q₁からq_k-1と直交するv_kを

v_k = a_k

for j=1 to k-1  v_k = v_k - (q_j^†a_k) q_j  …(1)

∴　v_k = a_k - (q₁^†a_k) q₁ - (q₂^†a_k) q₂ - ... - (q_k-1^†a_k) q_k-1  …(1')

と計算しているので、例えばq₁とq₂の直交性の誤差がq₂とv_kの直交性に伝搬してしまう、ということが起きる。

すなわち、たまたまa₁とa₂が直交に近いと

v₂ = a₂ - (q₁^†a₂) q₁   where q₁ = a₁/||a₁||

の内積q₁^†a₂が0に近いため、引き算で桁落ちしてq₁とq₂の直交性に誤差が生じる。

q₂ ^† q₁ =ε

このとき(1')のv_kとq₂の内積をとるとわかるように

q₂ ^† v_k= - (q₁^†a_k)ε

とq₂とv_kの直交性に誤差が伝搬してしまう。

そこで(1)を次の(2)のように変えると

v_k = a_k

for j=1 to k-1  v_k = v_k - (q_j^†v_k) q_j  …(2)

∴　v_k = a_k - (q₁^†a_k) q₁ - (q₂^†(a_k - (q₁^†a_k) q₁)) q₂ - ... - (q_k-1^†a_k) q_k-1  …(2')

この(2')のv_kとq₂の内積をとるとわかるように

q₂ ^† v_k= 0

とq₁とq₂の直交性の誤差はキャンセルされる。

このように、オリジナルのGram-Schmidt法の(1)を(2)に修正して、q_jどうしの直交性の誤差がq_jとv_kの直交性に伝搬しないようにしたものを修正Gram-Schmidt法という。

---

[←]　[↑]　[→]