Householder変換

寄稿:東條遼平

ベクトルの変換

例えば,

img116.png

のようなベクトルがあったときに,対称行列Pを考え,

img117.png

のような計算をし,

img118.png

のように第一要素以外がゼロであるようなベクトルに変換することができます.

これはx'とxが対象となるようなある超平面の法線ベクトルをvとしたときに,

img119.png

として計算された行列Pを使えばいいのですが,次の図を見るとわかりやすいかと思います.

householder1.png

便宜上2次元で表記していますが,超平面といっても高次元で平面のようなものを 考えているだけですので問題はありません.

そこでベクトルvの単位ベクトルとして

img120.png

を考えると,

img121.png

より,ベクトルの足し算を考えると

img122.png

ここで行列を導入します.

img123.png

ここで単位行列やらなにやらが湧いてでてきましたが, 式の意味が変わらないからこう変形しているだけです.こうすることで,

img124.png

とでき,

img119.png

とおくと

img117.png

となります.

ところでこの変換後のベクトルx'ですが,要素は何でも良いわけではなく, 変換前のxと大きさが同じでなくては困ります.そのためx'の第一要素は |x|となり,残りがゼロなのですが,符号は融通がききます. というのも,2次元以上ですので第一成分の符号が正だろうが負だろうが ベクトルvを調節すれば条件を満たす超平面が存在できるからです. そのため毎回正でも負でも数学的には問題ないのですが, PCで実装するとなると丸め誤差というものが発生するため ゼロに近い値というのは演算を進めて行くと正確さに欠けてしまいます. ベクトルvを求めるときにxとx'で引き算を行いますので, 第一成分の絶対値がゼロから離れるようにxとx'の第一成分は符号が異なる方が望ましいという訳です.

ここでベクトルx,x'を

img125.png
img126.png

とすると,ベクトルvは定数aを使って

img127.png

とできます.次に

img128.png
img129.png

とできますが,ベクトルvは向きが大事なのであって,大きさはいくらでも構いません. ですが,後々行列の三重対角化などを行うときに|v|=1の方が都合が良いので

img130.png

としておきます.

ソースコード

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に置き換えられ,返り値は 変換後のベクトルの第一成分となります. プログラム自体はかなりシンプルになっています.

Valid XHTML 1.1! home > コンピュータ > プログラミング >
リロード   新規 編集 凍結 差分 添付 複製 改名   トップ 一覧 検索 最終更新 バックアップ   ヘルプ   最終更新のRSS
Modified by 物理のかぎプロジェクト PukiWiki 1.4.5_1 Copyright © 2001-2005 PukiWiki Developers Team. License is GPL.
Based on "PukiWiki" 1.3 by yu-jiPowered by PHP 5.2.17HTML convert time to 0.256 sec.