三重対角化
をテンプレートにして作成
home
>
サイトマップ
開始行:
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)
というように変換します。このとき変換行列は超平面の法線ベ...
#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行3...
要素の値はHouseholder法によって既に値はわかっているので、...
右下の2×2の小行列を計算すれば良いことになります。また、...
変換前の右下2×2の小行列の両側から、変換行列の右下2×2...
これによって大幅に演算コストを下げる事ができます。
以上だけでも随分違うのですが、まだ行列の両側から行列を掛...
なのでできるだけまとめて計算できるようにしてみます。
#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->w...
for(j=i+2; j<m->width; j++){
m->a[i*m->width + j] = m->a[j*m->width...
}
}
FreeVector(v);
}
Householder変換するベクトルの次元は毎回1つずつ減っていき...
ループする度にデクリメントしています。 Householderの関数...
Householder変換後のベクトルの第一成分で、また、これは渡し...
大きさなのでこれがゼロであれば既にその行と列は3重対角化...
そのため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-...
}
}
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-...
g->v[i-offset]...
}
}
FreeVector(g);
}
殆どひたすら
#ref(img148.png,nolink)
を求めています。
#ref(img144.png,nolink)
のAが行列全体ではなくて値が変更される小行列を指しているた...
添字を底上げしていますが、殆ど上の式そのままです。
- &ref(main.c);
- &ref(calculation.c);
- &ref(calculation.h);
終了行:
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)
というように変換します。このとき変換行列は超平面の法線ベ...
#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行3...
要素の値はHouseholder法によって既に値はわかっているので、...
右下の2×2の小行列を計算すれば良いことになります。また、...
変換前の右下2×2の小行列の両側から、変換行列の右下2×2...
これによって大幅に演算コストを下げる事ができます。
以上だけでも随分違うのですが、まだ行列の両側から行列を掛...
なのでできるだけまとめて計算できるようにしてみます。
#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->w...
for(j=i+2; j<m->width; j++){
m->a[i*m->width + j] = m->a[j*m->width...
}
}
FreeVector(v);
}
Householder変換するベクトルの次元は毎回1つずつ減っていき...
ループする度にデクリメントしています。 Householderの関数...
Householder変換後のベクトルの第一成分で、また、これは渡し...
大きさなのでこれがゼロであれば既にその行と列は3重対角化...
そのため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-...
}
}
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-...
g->v[i-offset]...
}
}
FreeVector(g);
}
殆どひたすら
#ref(img148.png,nolink)
を求めています。
#ref(img144.png,nolink)
のAが行列全体ではなくて値が変更される小行列を指しているた...
添字を底上げしていますが、殆ど上の式そのままです。
- &ref(main.c);
- &ref(calculation.c);
- &ref(calculation.h);
ページ名:
home
>
Modified by
物理のかぎプロジェクト
PukiWiki 1.4.5_1
Copyright © 2001-2005
PukiWiki Developers Team
. License is
GPL
.
Based on "PukiWiki" 1.3 by
yu-ji
Powered by PHP 5.3.29HTML convert time to 0.002 sec.