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(m≧n)の行列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法という。


[] [] []