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

[] [] []


QR分解のアルゴリズム

A=QR

【修正Gram-Schmidt , ループkj順】

for k=1 to n

    /* AQに入れる */

    for i=1 to m

        qik = aik

    /* QRの計算 */

    for j=1 to k-1

        rjk = 0

        for i=1 to m

            rjk = rjk + qij * qik ()

        for i=1 to m

            qik = qik - rjk qij

    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順】

/* AQに入れる */

for j=1 to n

    for i=1 to m

        qij = aij

/* QRの計算 */

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 - rjk qij

()オリジナルのGram-Schmidt法では()qik の代わりにaik を使う。rjk = rjk + qij * aik

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

/* c=Qbを計算する*/

for j=1 to n

   cj = 0

   for i=1 to m

      cj = cj + qij* bi

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

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

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

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

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

a1 = r11 q1

a2 = r12 q1 + r22 q2

a3 = r13 q1 + r23 q2 + r33 q3

これをqjについて解くと

q1 = a1/r11

q2 =(a2 - r12 q1)/ r22

q3 =(a3 - r13 q1 - r23 q2)/ r33

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

r11 = ||a1||

r12 = q1a2 , r22 = ||a2 - r12 q1||

r13 = q1a3 , r23 = q2a3 , r33 = ||a3 - r13 q1 - r23 q2||

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

v1 = a1 , q1 = v1/||v1||

v2 =(a2 - (q1a2)q1) , q2 = v2/||v2||

v3 =(a3 - (q1a3)q1 - (q2a3)q2) , q3 = v3/||v3||

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

手順

コード

・ステップk (k=1 to n)において

q1からqk-1と直交するqkを求める

for k=1 to n

    vk = ak

    for j=1 to k-1

        rjk = qjak

       vk = vk - rjk qj

   rkk = ||vk||

    qk = vk /rkk

変数を節約するためvの代わりにqを使うとよい。

【古典Gram-Schmidt , ループkj順】

for k=1 to n

    /* AQに入れる */

    for i=1 to m

        qik = aik

    /* QRの計算 */

    for j=1 to k-1

        rjk = 0

        for i=1 to m

            rjk = rjk + qij * aik

        for i=1 to m

            qik = qik - rjk qij

    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 = qjak

        vk = vk - rjk qj

変数を節約するためvの代わりにqを使うとよい。

【古典Gram-Schmidt , ループjk順】

/* AQに入れる */

for j=1 to n

    for i=1 to m

        qij = aij

/* QRの計算 */

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 - rjk qij

◇◇

上記のGram-Schmidt法では、q1からqk-1と直交するvk

vk = ak

for j=1 to k-1  vk = vk - (qjak) qj  (1)

∴ vk = ak - (q1ak) q1 - (q2ak) q2 - ... - (qk-1ak) qk-1  (1')

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

すなわち、たまたまa1a2が直交に近いと

v2 = a2 - (q1a2) q1   where q1 = a1/||a1||

の内積q1a20に近いため、引き算で桁落ちしてq1q2の直交性に誤差が生じる。

q2 q1 =ε

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

q2 vk= - (q1ak)ε

q2vkの直交性に誤差が伝搬してしまう。

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

vk = ak

for j=1 to k-1  vk = vk - (qjvk) qj  (2)

∴ vk = ak - (q1ak) q1 - (q2(ak - (q1ak) q1)) q2 - ... - (qk-1ak) qk-1  (2')

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

q2 vk= 0

q1q2の直交性の誤差はキャンセルされる。

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


[] [] []