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なら１回rootでMaximaを起動してload(lapack);とやれば良い）