QR法
をテンプレートにして作成
home
>
サイトマップ
開始行:
RIGHT:寄稿:東條遼平
* 高速な対称行列の対角化 [#gd8fa78d]
実対称行列をJacobi法よりも高速に対角化し、固有値を求める...
Jacobi法は非対角成分を一つ選び、その成分をゼロにするよう...
もぐら叩きのような方法を取るため、どこかをゼロにしてもま...
あまり効率がよくありません。その点、QR法は事前に三重対角...
三重対角化された状態を保ったまま副対角成分をゼロにする方...
また、右下2×2の小行列の固有値を対角成分から先に引いてか...
対角化への収束が早くなります。
* 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列目...
行列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ができます。この行列は先程ゼ...
ゼロで無くなっているのですが、かなりゼロに近い値に変わっ...
これを1セットとして何度か繰り返し、4行3列の要素が許容誤...
4行目、4列目には手を着けず、左上3×3の小行列だと考えて繰り...
すると次は3行2列が許容誤差以下になり、最終的に2行1列が許...
対角化が終了します。このように次元がだんだんと小さくなっ...
より高速に演算が行えます。
* 問題点 [#gbbc2773]
以上でQR法により対角化することができますが、このまま実装...
それはだんだんと対角化が進み、対角成分が固有値に近付いて...
副対角成分の差が大きくなり、副対角成分の値がゼロに近付き...
そのため、次のように、全ての対角成分からある定数分だけ引...
後から足し直す事にします。
#ref(img171.png,nolink)
#ref(img172.png,nolink)
実際にこの上の式を下の式に代入するとわかりますが、
先程までの式と全く同じになります。また、この定数の値をど...
扱っている行列Aの右下2×2の小行列の固有値が良いようです。
2×2なら解の公式よりすぐに求める事ができますし、対角化が...
副対角成分がゼロに近付いても、同じく対角成分も2×2の小行...
固有値を引いているので一緒にゼロに近付きます。そのため
高速なまま対角化がすすめられます。
* ソースコード [#w612cb07]
Tridiag(m);
nowsize = m->width;
q = CreateMatrix(nowsize, nowsize);
v = CreateVector(nowsize);
まず対角化したい行列mを三重対角化しています。
処理が進ごとに次元が下がって行きますので現在の次元をnowsi...
行列qは先程の行列Q、ベクトルvはテンポラリ−として使ってい...
while(nowsize > 1){
double det, b, u;
if(fabs(m->a[(nowsize-1)*m->width + nowsize-2]...
nowsize--;
continue;
}
u = GetEigenvalue22(m->a[(nowsize-2)*m->width ...
m->a[(nowsize-2)*m->width ...
m->a[(nowsize-1)*m->width ...
m->a[(nowsize-1)*m->width ...
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-...
}
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は現在扱っている行列の次元なのでこれが1になると対...
当り前なのでwhile文の中からいきます。
if(fabs(m->a[(nowsize-1)*m->width + nowsize-2]) &l...
nowsize--;
continue;
}
この部分は例えば現在4×4の行列であれば
4行3列、つまり下側の副対角成分の一番下の部分が許容誤差...
左上の3×3の小行列のみを考えるため次元を1つ小さくします。
u = GetEigenvalue22(m->a[(nowsize-2)*m->width + no...
m->a[(nowsize-2)*m->width + no...
m->a[(nowsize-1)*m->width + no...
m->a[(nowsize-1)*m->width + no...
for(i=0; i<nowsize; i++){
m->a[i*m->width + i] -= u;
}
GetEigenvalue22の関数に右下2×2の小行列の4つの値を渡し...
値に近い方を返しています。そしてこの値をそれぞれの対角成...
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);
終了行:
RIGHT:寄稿:東條遼平
* 高速な対称行列の対角化 [#gd8fa78d]
実対称行列をJacobi法よりも高速に対角化し、固有値を求める...
Jacobi法は非対角成分を一つ選び、その成分をゼロにするよう...
もぐら叩きのような方法を取るため、どこかをゼロにしてもま...
あまり効率がよくありません。その点、QR法は事前に三重対角...
三重対角化された状態を保ったまま副対角成分をゼロにする方...
また、右下2×2の小行列の固有値を対角成分から先に引いてか...
対角化への収束が早くなります。
* 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列目...
行列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ができます。この行列は先程ゼ...
ゼロで無くなっているのですが、かなりゼロに近い値に変わっ...
これを1セットとして何度か繰り返し、4行3列の要素が許容誤...
4行目、4列目には手を着けず、左上3×3の小行列だと考えて繰り...
すると次は3行2列が許容誤差以下になり、最終的に2行1列が許...
対角化が終了します。このように次元がだんだんと小さくなっ...
より高速に演算が行えます。
* 問題点 [#gbbc2773]
以上でQR法により対角化することができますが、このまま実装...
それはだんだんと対角化が進み、対角成分が固有値に近付いて...
副対角成分の差が大きくなり、副対角成分の値がゼロに近付き...
そのため、次のように、全ての対角成分からある定数分だけ引...
後から足し直す事にします。
#ref(img171.png,nolink)
#ref(img172.png,nolink)
実際にこの上の式を下の式に代入するとわかりますが、
先程までの式と全く同じになります。また、この定数の値をど...
扱っている行列Aの右下2×2の小行列の固有値が良いようです。
2×2なら解の公式よりすぐに求める事ができますし、対角化が...
副対角成分がゼロに近付いても、同じく対角成分も2×2の小行...
固有値を引いているので一緒にゼロに近付きます。そのため
高速なまま対角化がすすめられます。
* ソースコード [#w612cb07]
Tridiag(m);
nowsize = m->width;
q = CreateMatrix(nowsize, nowsize);
v = CreateVector(nowsize);
まず対角化したい行列mを三重対角化しています。
処理が進ごとに次元が下がって行きますので現在の次元をnowsi...
行列qは先程の行列Q、ベクトルvはテンポラリ−として使ってい...
while(nowsize > 1){
double det, b, u;
if(fabs(m->a[(nowsize-1)*m->width + nowsize-2]...
nowsize--;
continue;
}
u = GetEigenvalue22(m->a[(nowsize-2)*m->width ...
m->a[(nowsize-2)*m->width ...
m->a[(nowsize-1)*m->width ...
m->a[(nowsize-1)*m->width ...
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-...
}
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は現在扱っている行列の次元なのでこれが1になると対...
当り前なのでwhile文の中からいきます。
if(fabs(m->a[(nowsize-1)*m->width + nowsize-2]) &l...
nowsize--;
continue;
}
この部分は例えば現在4×4の行列であれば
4行3列、つまり下側の副対角成分の一番下の部分が許容誤差...
左上の3×3の小行列のみを考えるため次元を1つ小さくします。
u = GetEigenvalue22(m->a[(nowsize-2)*m->width + no...
m->a[(nowsize-2)*m->width + no...
m->a[(nowsize-1)*m->width + no...
m->a[(nowsize-1)*m->width + no...
for(i=0; i<nowsize; i++){
m->a[i*m->width + i] -= u;
}
GetEigenvalue22の関数に右下2×2の小行列の4つの値を渡し...
値に近い方を返しています。そしてこの値をそれぞれの対角成...
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);
ページ名:
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.