NaLab02:行列・ベクトルの演習

投稿者: | 2018-05-06

今回のB4計算機ゼミでは,行列とベクトルを使った数値実験の演習のために「べき乗法」(Power Method)の課題が出た.

べき乗法とは

べき乗法は,固有値問題を取り扱う際に用いられる数値解法のひとつであり,ある行列Aと初期ベクトルを与えることにより,絶対値最大の固有値とそれに対応する固有ベクトルを1つだけ求める手法である.

べき乗法のアルゴリズム

べき乗法のアルゴリズムは以下の通りである.

・初期ベクトルx(0)を選び,x(k)=x(k)/r(k)と正規化する.ここで,反復回数をk=0と初期化しておく.また,r(k)をk番目の時のxのノルムとする.【*】
x(k+1)=Ax(k)を計算する.
・r(k+1)=r(k+1) / r(k)を計算する.
・|r(k) ー r(k+1)| ≦ |r(k+1)|を収束条件とする.
・これを満たさない場合は,k = k+1 として,【*】からの反復を繰り返す.

プログラムの部品

べき乗法は,上記のアルゴリズムで実行可能だが,プログラムは部品を関数(ルーチン)化しておき,今後も利用できるようにしておくことが望ましい.今回は行列とベクトル積,定数倍(除算)することがあるので,それぞれをルーチン化しておいた.

ただし,mainで宣言してmallocした動的メモリ領域を利用しているため,関数にもポインタ渡しで動作させていることに気を付けたい.
代入先のベクトル要素もあらかじめ初期化されていることが前提.

行列とベクトルの積

行列・ベクトル積は,行列の列数とベクトルの長さが一致しなければならない.

ベクトルのノルム・ベクトルの内積

よくよく見てみれば,中身は同じ動作をしていることに気づく.

ベクトルの定数倍(除算ver.)

シンプルに.

こんな部品さえ作っておけば,今後いちからプログラムを書かなくても,部品をくっつけてプログラムを構築していけばよい.

べき乗法のプログラム

プログラムすべてを表すと,以下のようになる.
ここでは,行列をランダムで生成し,対称行列になるよう構築している.
また,最終表示は,反復回数・絶対値最大の固有値・固有ベクトルを表示させ,固有ベクトルに関しては,1つめの要素で正規化して表示させるようにもしてみた.(固有ベクトルは任意の定数倍でもよいことを考慮している.手計算で求めるときに任意定数を置くことを思い出してほしい.)

留意点

アルゴリズム中で正規化をしないと,オーバーフローもしくはアンダーフローしてしまうため,計算が途中で終了してしまう.正規化はそれを避けるためも改善策である.
また,初期ベクトルの置き方によっては,収束せずに発散してしまうこともある.

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です