RIGHT:寄稿：東條遼平 * 三重対角化とは [#a4995ab5] 三重対角化とは、 #ref(img131.png,nolink) というように、対称行列を変換し、対角成分とその両隣の副対角成分以外を全てゼロにする 変換です。固有値を求めるときに、三重対角化を経由して対角化する方が一般的に高速になります。 イメージ的にも行列の成分でゼロが増えれば計算しなくてはならない量が減りそうです。 具体的には変換したい対称行列に何らかの行列を掛けて三重対角化する事になりますが、 どんな行列でも良いという訳ではなく、直交型の行列で、しかも右と左の両側から 掛けていかなければなりません。そうしないと変換後の行列の固有値が変わってしまいます。 つまり、式で書くと次のようにしなければなりません。Pは直交行列です。 #ref(img132.png,nolink) 例えば上の行列であれば、1行目と1列目を見てみると、行ベクトル、列ベクトル共に a、b、c、dというベクトルがk、l、0、0に変換されていると言えます。 もっというと、aという要素をそのまま置いておき(k=a)、 b、c、dというベクトルがl、0、0というベクトルに変わったと言えます。 ベクトルの第一成分以外をゼロにする変換行列は Householder変換 で求める事ができます。行列Aは対称行列なので行と列共に同じ変換行列で変換でき、 左から掛ければ列ベクトルが変換され、右から掛ければ行ベクトルが変換されます。 * 流れ [#b69ebcd8] #ref(img133.png,nolink) という行列Aを三重対角化する流れを説明します。これは対称行列です。まず、 #ref(img134.png,nolink) をHouseholder変換によって #ref(img135.png,nolink) というように変換します。このとき変換行列は超平面の法線ベクトルvを使って #ref(img119.png,nolink) と計算できますので、 #ref(img136.png,nolink) のような行列になります。これを行列Aの両側から掛けると、 #ref(img137.png,nolink) となります。続いて、 #ref(img138.png,nolink) をHouseholder変換によって #ref(img139.png,nolink) と変換し、変換行列は #ref(img140.png,nolink) となります。これを両側から掛けると #ref(img141.png,nolink) となり三重対角化が完了します。 * 問題点 [#j4162b6c] 以上で三重対角化が出来、Householder変換が出来ていれば 大して難しいものではありません。しかし、そもそも三重対角化をしようとしているのは 高速に行列を対角化するために経由させているので、三重対角化自体に 大きく時間が掛かかるような事は避けたいものです。そのため少し無駄を省く事にします。 まず行列の掛け算は演算コストがO(N3)になります。 これはN次元の行列の掛け算をするとN3回の演算が必要だという事です。 そのため次元の大きい行列の掛け算はできるだけ減らしたいということになります。 ここで、上の流れにおいて、 #ref(img137.png,nolink) から #ref(img141.png,nolink) に変換するとき、両側から #ref(img140.png,nolink) を掛ける事になりますが、このとき変換前と後で、1行目と1列目、2行2列目の成分は 変化していません。残りの部分は変化しますが、そのうち2行3列、2行4列、3行2列、4行2列の 要素の値はHouseholder法によって既に値はわかっているので、結局 右下の２×２の小行列を計算すれば良いことになります。また、変換後の右下２×２の小行列は 変換前の右下２×２の小行列の両側から、変換行列の右下２×２の小行列を掛けたものと同じです。 これによって大幅に演算コストを下げる事ができます。 以上だけでも随分違うのですが、まだ行列の両側から行列を掛けるという演算が残ります。 なのでできるだけまとめて計算できるようにしてみます。 #ref(img144.png,nolink) #ref(img148.png,nolink) とおくと、 #ref(img145.png,nolink) #ref(img150.png,nolink) より、 #ref(img142.png,nolink) #ref(img143.png,nolink) #ref(img146.png,nolink) #ref(img147.png,nolink) #ref(img149.png,nolink) * ソースコード [#sc73fd81] 三重対角化の部分は以下のようになります。 void Tridiag(Matrix *m) { Vector *v; int i, j; double tmp; v = CreateVector(m->width); for(i=0; i<m->width-2; i++){ v->size--; for(j=i+1; j<m->width; j++){ v->v[j-i-1] = m->a[i*m->width + j]; } if(!(tmp = Householder(v))) continue; Mult(m, v); m->a[i*m->width + i+1] = m->a[(i+1)*m->width + i] = tmp; for(j=i+2; j<m->width; j++){ m->a[i*m->width + j] = m->a[j*m->width + i] = 0; } } FreeVector(v); } Householder変換するベクトルの次元は毎回１つずつ減っていきますので、 ループする度にデクリメントしています。 Householderの関数の戻り値は、 Householder変換後のベクトルの第一成分で、また、これは渡したベクトルの 大きさなのでこれがゼロであれば既にその行と列は３重対角化されている事になります。 そのためcontinueしてループの先頭に戻ります。 Mult関数は、Householder変換の際、超平面の法線ベクトルvを受け取り、 行列mの両側からvを使って作られた変換行列を掛けます。このとき 右下の小行列のみ更新しますので、その後Householder変換した行列の成分を更新します。 static void Mult(Matrix* x, Vector* v) { int i, j, offset; double tmp; Vector *g; offset = x->width - v->size; g = CreateVector(v->size); for(i=offset; i<x->width; i++){ g->v[i-offset] = 0; for(j=offset; j<x->width; j++){ g->v[i-offset] += x->a[i*x->width + j]*v->v[j-offset]; } } tmp = GetDotProduct(g, v); for(j=0; j<g->size; j++){ g->v[j] = 2*(g->v[j] - v->v[j]*tmp); } for(i=offset; i<x->width; i++){ for(j=offset; j<x->width; j++){ x->a[i*x->width +j] -= (v->v[i-offset]*g->v[j-offset] + g->v[i-offset]*v->v[j-offset]); } } FreeVector(g); } 殆どひたすら #ref(img148.png,nolink) を求めています。 #ref(img144.png,nolink) のAが行列全体ではなくて値が変更される小行列を指しているためoffsetを使って 添字を底上げしていますが、殆ど上の式そのままです。 - &ref(main.c); - &ref(calculation.c); - &ref(calculation.h);
