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);
