RIGHT:寄稿：東條遼平 * 高速な対称行列の対角化 [#gd8fa78d] 実対称行列をJacobi法よりも高速に対角化し、固有値を求めることができる方法です。 Jacobi法は非対角成分を一つ選び、その成分をゼロにするような直交行列を掛ける事で、 もぐら叩きのような方法を取るため、どこかをゼロにしてもまたどこかがゼロでなくなり、 あまり効率がよくありません。その点、QR法は事前に三重対角化することで、 三重対角化された状態を保ったまま副対角成分をゼロにする方法を取ります。 また、右下２×２の小行列の固有値を対角成分から先に引いてから変換するため、 対角化への収束が早くなります。 * QR法の考え [#a099ffc4] 行列Aに対し、直交行列Pを次々と左から掛けて上三角行列R(対角成分より下が全てゼロ)を 作ったとします。 #ref(img151.png,nolink) このとき、単位行列に右から同じものを掛けていったものを #ref(img152.png,nolink) とすると、Qは直交行列であり、 #ref(img153.png,nolink) とのように行列Aを直交行列と上三角行列に分解することが出来ます。 このとき、行列QとRを入れ換えた、 #ref(img157.png,nolink) という新たな行列を考えたとき、 #ref(img154.png,nolink) と変換出来ることから、この新たな行列の固有値は元の行列Aの固有値と等しくなります。 つまりこれを何度も繰り返すことで、対角化すればその対角成分を調べる事で、 行列Aの固有値を知ることができる事になります。一般的に書くと次のようになります。 #ref(img155.png,nolink) と分解し、新たな行列として #ref(img158.png,nolink) #ref(img156.png,nolink) とします。 * 流れ [#p4d45b0a] #ref(img159.png,nolink) という三重対角化行列を考えると、これを上三角行列Rに変換しなければなりません。 しかし、上三角行列といっても対角成分よりも下側が全てゼロであればいいので、 実際には副対角成分の片側だけをゼロにすればいいだけです。 行列Aでいうと2行1列、3行2列、4行3列の3つです。 それではまず2行1列の値をゼロにしてみます。 このとき掛ける直交行列の見つけ方はJacobi法と同じで回転行列を考えます。 #ref(img160.png,nolink) という行列を左から掛けると、 #ref(img161.png,nolink) ここで2行1列の要素がゼロになるようにsinとcosを選びます。 #ref(img162.png,nolink) #ref(img163.png,nolink) #ref(img164.png,nolink) #ref(img165.png,nolink) #ref(img166.png,nolink) また、このとき1行目と2行目しか値が変化しない事を利用すると行列の掛け算の演算コストを 大幅に下げる事ができます。Qも同じく右から掛けるので1列目と2列目しか値が変化しません。 行列A、行列Qのそれぞれの新しい値は次のようになります。 #ref(img167.png,nolink) #ref(img168.png,nolink) すると行列Aは次のようになります。 #ref(img169.png,nolink) 今度は3行2列の要素をゼロにするのですが、直交行列Pは次のようになります。 #ref(img170.png,nolink) このように4行3列までゼロにしたら直交行列Qと上三角行列Rへの分解が終了となりますので、 #ref(img157.png,nolink) を計算して新たな三重対角行列Aができます。この行列は先程ゼロにした要素が ゼロで無くなっているのですが、かなりゼロに近い値に変わっています。 これを１セットとして何度か繰り返し、4行3列の要素が許容誤差以下になれば 4行目、4列目には手を着けず、左上3×3の小行列だと考えて繰り返して行きます。 すると次は3行2列が許容誤差以下になり、最終的に2行1列が許容誤差以下になって 対角化が終了します。このように次元がだんだんと小さくなって行くことで より高速に演算が行えます。 * 問題点 [#gbbc2773] 以上でQR法により対角化することができますが、このまま実装すると不都合な点があります。 それはだんだんと対角化が進み、対角成分が固有値に近付いて行くと対角成分と 副対角成分の差が大きくなり、副対角成分の値がゼロに近付きにくくなるのです。 そのため、次のように、全ての対角成分からある定数分だけ引いておき、 後から足し直す事にします。 #ref(img171.png,nolink) #ref(img172.png,nolink) 実際にこの上の式を下の式に代入するとわかりますが、 先程までの式と全く同じになります。また、この定数の値をどう選ぶかですが、 扱っている行列Aの右下２×２の小行列の固有値が良いようです。 ２×２なら解の公式よりすぐに求める事ができますし、対角化が進んで 副対角成分がゼロに近付いても、同じく対角成分も２×２の小行列の 固有値を引いているので一緒にゼロに近付きます。そのため 高速なまま対角化がすすめられます。 * ソースコード [#w612cb07] Tridiag(m); nowsize = m->width; q = CreateMatrix(nowsize, nowsize); v = CreateVector(nowsize); まず対角化したい行列mを三重対角化しています。 処理が進ごとに次元が下がって行きますので現在の次元をnowsizeに保持します。 行列qは先程の行列Q、ベクトルvはテンポラリ−として使っています。 while(nowsize > 1){ double det, b, u; if(fabs(m->a[(nowsize-1)*m->width + nowsize-2]) < eps){ nowsize--; continue; } u = GetEigenvalue22(m->a[(nowsize-2)*m->width + nowsize-2], m->a[(nowsize-2)*m->width + nowsize-1], m->a[(nowsize-1)*m->width + nowsize-2], m->a[(nowsize-1)*m->width + nowsize-1]); for(i=0; i<nowsize; i++){ m->a[i*m->width + i] -= u; } q->height = nowsize; q->width = nowsize; InitMatrix(q); NextAQ(m, q, nowsize); for(i=0; i<nowsize; i++){ double sum=0; int k; for(j=0; j<nowsize; j++){ for(k=i; k<nowsize; k++){ sum += m->a[i*m->width + k]*q->a[k*q->width + j]; } v->v[j] = sum; sum = 0; } for(j=0; j<nowsize; j++){ m->a[i*m->width + j] = v->v[j]; } } for(i=0; i<nowsize; i++){ m->a[i*m->width + i] += u; } } これがQR法の心臓部の全てですが、長いのでちょっとずつ見て行きます。 nowsizeは現在扱っている行列の次元なのでこれが１になると対角化が終了するのは 当り前なのでwhile文の中からいきます。 if(fabs(m->a[(nowsize-1)*m->width + nowsize-2]) < eps){ nowsize--; continue; } この部分は例えば現在４×４の行列であれば ４行３列、つまり下側の副対角成分の一番下の部分が許容誤差以下になっていれば 左上の３×３の小行列のみを考えるため次元を１つ小さくします。 u = GetEigenvalue22(m->a[(nowsize-2)*m->width + nowsize-2], m->a[(nowsize-2)*m->width + nowsize-1], m->a[(nowsize-1)*m->width + nowsize-2], m->a[(nowsize-1)*m->width + nowsize-1]); for(i=0; i<nowsize; i++){ m->a[i*m->width + i] -= u; } GetEigenvalue22の関数に右下２×２の小行列の４つの値を渡し、その固有値のうち、一番右下の 値に近い方を返しています。そしてこの値をそれぞれの対角成分から一時引いておきます。 q->height = nowsize; q->width = nowsize; InitMatrix(q); NextAQ(m, q, nowsize); qのサイズを更新してInitMatrixで行列qを単位行列に初期化します。 NextAQが終わるとmは上三角行列R、qは直交行列Qになります。 これより後は新たなAを求めるためにRQの順でかけ算し、 先程引いておいたμを足し直しているだけです。 NextAQの内部は少しごちゃごちゃしてしまいましたが、 sinとcosを求めて、 #ref(img167.png,nolink) #ref(img168.png,nolink) を計算し、渡された行列mが上三角行列になるまで計算しているに過ぎません。 - &ref(main.c); - &ref(calculation.c); - &ref(calculation.h);
