Householder変換 をテンプレートにして作成

開始行:
RIGHT:寄稿：東條遼平 * ベクトルの変換 [#g2445ea2] 例えば、 #ref(img116.png,nolink) のようなベクトルがあったときに、対称行列Pを考え、 #ref(img117.png,nolink) のような計算をし、 #ref(img118.png,nolink) のように第一要素以外がゼロであるようなベクトルに変換する... これはx'とxが対象となるようなある超平面の法線ベクトルをv... #ref(img119.png,nolink) として計算された行列Pを使えばいいのですが、次の図を見ると... #ref(householder1.png,nolink) 便宜上2次元で表記していますが、超平面といっても高次元で平... 考えているだけですので問題はありません。 そこでベクトルvの単位ベクトルとして #ref(img120.png,nolink) を考えると、 #ref(img121.png,nolink) より、ベクトルの足し算を考えると #ref(img122.png,nolink) ここで行列を導入します。 #ref(img123.png,nolink) ここで単位行列やらなにやらが湧いてでてきましたが、 式の意味が変わらないからこう変形しているだけです。こうす... #ref(img124.png,nolink) とでき、 #ref(img119.png,nolink) とおくと #ref(img117.png,nolink) となります。 ところでこの変換後のベクトルx'ですが、要素は何でも良いわ... 変換前のxと大きさが同じでなくては困ります。そのためx'の第... |x|となり、残りがゼロなのですが、符号は融通がききます。 というのも、2次元以上ですので第一成分の符号が正だろうが負... ベクトルvを調節すれば条件を満たす超平面が存在できるからで... そのため毎回正でも負でも数学的には問題ないのですが、 PCで実装するとなると丸め誤差というものが発生するため ゼロに近い値というのは演算を進めて行くと正確さに欠けてし... ベクトルvを求めるときにxとx'で引き算を行いますので、 第一成分の絶対値がゼロから離れるようにxとx'の第一成分は符... ここでベクトルx、x'を #ref(img125.png,nolink) #ref(img126.png,nolink) とすると、ベクトルvは定数aを使って #ref(img127.png,nolink) とできます。次に #ref(img128.png,nolink) #ref(img129.png,nolink) とできますが、ベクトルvは向きが大事なのであって、大きさは... ですが、後々行列の三重対角化などを行うときに|v|=1の方が都... #ref(img130.png,nolink) としておきます。 * ソースコード [#q9f447d2] double Householder(Vector *x) { double norm; double weight; int i, j; norm = sqrt(GetDotProduct(x, x)); if(norm){ if(x->v[0] < 0) norm = - norm; x->v[0] += norm; weight = 1/(sqrt(2*norm*x->v[0])); for(i=0; i<x->size; i++){ x->v[i] *= weight; } } return -norm; } GetDotProductは内積を取る関数です。 引数として変換前のベクトルを渡し、ベクトルvに置き換えられ... 変換後のベクトルの第一成分となります。 プログラム自体はかなりシンプルになっています。 - &ref(main.c); - &ref(calculation.c); - &ref(calculation.h);

終了行:
RIGHT:寄稿：東條遼平 * ベクトルの変換 [#g2445ea2] 例えば、 #ref(img116.png,nolink) のようなベクトルがあったときに、対称行列Pを考え、 #ref(img117.png,nolink) のような計算をし、 #ref(img118.png,nolink) のように第一要素以外がゼロであるようなベクトルに変換する... これはx'とxが対象となるようなある超平面の法線ベクトルをv... #ref(img119.png,nolink) として計算された行列Pを使えばいいのですが、次の図を見ると... #ref(householder1.png,nolink) 便宜上2次元で表記していますが、超平面といっても高次元で平... 考えているだけですので問題はありません。 そこでベクトルvの単位ベクトルとして #ref(img120.png,nolink) を考えると、 #ref(img121.png,nolink) より、ベクトルの足し算を考えると #ref(img122.png,nolink) ここで行列を導入します。 #ref(img123.png,nolink) ここで単位行列やらなにやらが湧いてでてきましたが、 式の意味が変わらないからこう変形しているだけです。こうす... #ref(img124.png,nolink) とでき、 #ref(img119.png,nolink) とおくと #ref(img117.png,nolink) となります。 ところでこの変換後のベクトルx'ですが、要素は何でも良いわ... 変換前のxと大きさが同じでなくては困ります。そのためx'の第... |x|となり、残りがゼロなのですが、符号は融通がききます。 というのも、2次元以上ですので第一成分の符号が正だろうが負... ベクトルvを調節すれば条件を満たす超平面が存在できるからで... そのため毎回正でも負でも数学的には問題ないのですが、 PCで実装するとなると丸め誤差というものが発生するため ゼロに近い値というのは演算を進めて行くと正確さに欠けてし... ベクトルvを求めるときにxとx'で引き算を行いますので、 第一成分の絶対値がゼロから離れるようにxとx'の第一成分は符... ここでベクトルx、x'を #ref(img125.png,nolink) #ref(img126.png,nolink) とすると、ベクトルvは定数aを使って #ref(img127.png,nolink) とできます。次に #ref(img128.png,nolink) #ref(img129.png,nolink) とできますが、ベクトルvは向きが大事なのであって、大きさは... ですが、後々行列の三重対角化などを行うときに|v|=1の方が都... #ref(img130.png,nolink) としておきます。 * ソースコード [#q9f447d2] double Householder(Vector *x) { double norm; double weight; int i, j; norm = sqrt(GetDotProduct(x, x)); if(norm){ if(x->v[0] < 0) norm = - norm; x->v[0] += norm; weight = 1/(sqrt(2*norm*x->v[0])); for(i=0; i<x->size; i++){ x->v[i] *= weight; } } return -norm; } GetDotProductは内積を取る関数です。 引数として変換前のベクトルを渡し、ベクトルvに置き換えられ... 変換後のベクトルの第一成分となります。 プログラム自体はかなりシンプルになっています。 - &ref(main.c); - &ref(calculation.c); - &ref(calculation.h);

ページ名: