« 統計学復習メモ8: 回帰直線と主成分分析 | メイン | EclipseにTomcatプラグイン導入 »

 Maximaでの固有ベクトル計算について

筆者のMaxima 5.17+sbclの環境だと、標準のeigenパッケージを使って固有ベクトルを求めようとすると、エラーが出て求まらないことが多い。例えば、

load(eigen);
m: matrix([0.123,0.456,0.789],[0.456,0.789,0],[0.789,0,0]);
eivects(m);

とやると、異様に時間がかかる上、
algsys failure: the eigenvector(s) for the1th eigenvalue will be missing.
algsys failure: the eigenvector(s) for the2th eigenvalue will be missing.
algsys failure: the eigenvector(s) for the3th eigenvalue will be missing.

というメッセージに続いて、虚数混じりのおぞましい値が固有値として出てきて、固有ベクトルは出力されない。そもそも固有値が間違っている。

固有ベクトルを英語でeigenvectorと言うことをたまたま覚えてたので、Maximaのヘルプでeigenvectorを検索して、出てきたのがeigenパッケージのeigenvectors()/eivects()だったので、それを使ったのだが、それが間違いだったようだ。Maximaのヘルプのeigenvectorsの項にも、固有値が正しく求められないことがあると、はっきり書いてある。対処方法も書いてあるが、とても面倒である。

インターネットで調べてみると、Maximaのメーリングリストのアーカイブに辿り着いた。どうやら2007年4月くらいからこの問題は原因も含めて知られているが、まだ解決していないらしい。

数式を含まない、純粋な数値行列なら、他の方法があることがわかった。
(1) eigens_by_jacobi()
標準関数として存在する。対称行列に限るという制約があるが、例えば共分散行列の固有ベクトルを求めるなら何の問題も無い。先に言ってくれという感じである。

(%i1) m: matrix([0.123,0.456,0.789],[0.456,0.789,0],[0.789,0,0]);
(%i2) eigens_by_jacobi(m);
(%o2) [[- .7932213800442011, 1.180843735860609, .5243776441835921],
[ .6946378770918785 .5975383492309946 .4005323219261915 ]
[ ]
[ - .2001963037214474 .6953728293011738 - .6902014693159894 ]]
[ ]
[ - .6909411405362675 .3992549930407618 .6026572747810012 ]

eivects()と同様、固有値の行列と固有ベクトルの行列が一緒に返される。

(2) LAPACKパッケージのdgeev()
昔からある最強の線形演算ライプラリであるLAPACKのLisp版がMaximaに同梱されている。最初に使用する時にコンパイルがなされるが、これが環境依存であり、Lispの実行環境によってはコンパイルに失敗してしまうらしい。無事にコンパイルに成功すると、dgeev()で最強の計算結果が得られる。

(%i1) load(lapack);
(%i2) m: matrix([0.123,0.456,0.789],[0.456,0.789,0],[0.789,0,0]);
(%i3) dgeev(m,true,true);
(%o3) [[- .7932213800442016, 1.180843735860611, .5243776441835922],
[ .6946378770918787 - .5975383492309948 .4005323219261915 ]
[ ]
[ - .2001963037214475 - .6953728293011737 - .6902014693159894 ],
[ ]
[ - .6909411405362672 - .3992549930407617 .6026572747810011 ]

この例は対称行列だが、もちろん非対称の行列でも問題ない。

筆者はgclとsbclで試したが、sbclではコンパイルが終了しなかったが、暴走したMaximaを強制終了すると、次からはコンパイル無しでロードできるようになった。(テスト段階で暴走していた?)

なお、root権限でインストールしたMaximaのLAPACKパッケージをコンパイルするには、root権限が必要である。(UNIX系のOSなら1回rootでMaximaを起動してload(lapack);とやれば良い)

トラックバック

このエントリーのトラックバックURL:
http://ynomura.dip.jp/cgi-bin/mt/mt-tb.cgi/116

コメント (2)

yuji:

eigen_by_jacobi関数、役に立ちました。
ありがとうございます。

一つ質問があるんですが、数式を含んだ状態で固有値を求めた際、虚数を含まないようにする命令はありますか?
A:matrix([-z,-1,0],[-1,-z/2,-0.25],[0,-0.25,-z/2]);
B:eigenvalues(A);
として固有値を求めると、対称行列なのに虚数が表われて、よく分からない数式になってしまいます。

ynomura:

eigenvalues()はeigenパッケージにあり、まだeigenvectors()と共にバグがありますので、不正な答えが出力されることがあります。
また、3x3の大きさの行列の固有値は3次方程式の解であり、3次方程式の解の公式は(実数解でも一時的に)虚数が含まれますので、行列に変数が残っていると、実数解しかあり得ないとわかってても、虚数抜きで表現するのは困難だと思われます。(ちなみにMathematicaだと、こういう場合は「カルダノの公式」が適用されて、やはり実数解が虚数を使って表現されるようです)

固有値は固有方程式の解であることを使って、Maximaで
solve(determinant(A-diagmatrix(3,x))=0,x);
とすると、虚数は入りますが、正しい答えが得られます。

もし、zに何か実数が代入されてから計算されても良いのであれば、
A(z):=matrix([-z,-1,0],[-1,-z/2,-0.25],[0,-0.25,-z/2]); B(z):=eigens_by_jacobi(A(z));
としてから
B(3)
などとすると、実数表示の解が得られます。(ここでも残念ながら、eigenvalues()を使うと虚数を含む解の公式表示になってしまうようです)

コメント投稿フォーム

※投稿されたコメントはオーナーが承認するまで表示されません。


Powered by Movable Type 3.35