第9話 固有値計算雑感 - 後編
前編では、固有値解析を含むアプリケーションプログラムにおいては、既存のMathlibの全面利用ではなく、自ら構築したプログラムスキームの中に固有値算法を取り込むことを勧めたところで終わりました。もっとも、マトリックスの次元Nが大きい場合の話ですが。
しかし、お勧めするとは言っても、前編で挙げたp≦NかつNが大次元という課題をカバーする一般固有値解析の算法といえば、選択肢に限りがあります。実際問題、次の三つの算法があるだけではないかと思います － 浅学な筆者ゆえ、他に存在することを知らずにいるかもしれませんが、その場合はお許しください。
いずれも初期ベクトルを設定してスタートする反復計算型の算法です。しかし、固有値、固有ベクトルの求め方が三者三様です。
以下、この三つの算法の概略説明および筆者がプログラム開発して利用した感想を述べてみます。
[インバース・パワー法]
最小固有値一つを求める算法で、日本語で逆反復法とも呼ばれています。理論といい、算法といい理解し易い方法なので、初中級レベルの数値計算プログラムの練習問題には格好の材料と言えるかもしれません。
孤立した固有値ならば、ものの見事に早く固有値を探し出す優れた算法ですが、問題は重複固有値の存在、近接固有値が存在する場合です(図1)。この二つが存在する場合は、初期ベクトルの設定にもよりますが、うまく収束しないことを覚悟する必要があります。
重複固有値、近接固有値というのは、例えば、3次のマトリックスを想定した場合、固有方程式を固有値変数λで展開した三次の固有多項式をイメージすれば分かりが早いです。三次関数の根がすなわち固有値ですから、重複固有値とは、図1左の場合で、数学書では縮退固有値とも呼んでいます。弾性構造物の振動解析でいえば、対称性を持つ構造では、対称度数に相応した同値の、すなわち重複固有値が存在するため、このようなときには結果の評価に注意が必要です。
近接固有値の方は、同じく図1右にあるように、根の値に違いがあるにはあるのですが、お互い非常に近接した隣接固有値のことをいいます。
数値計算の結果には、もちろん数値誤差を持つのが当然ですので、物理問題によっては、微小な差異を持つ複数の固有値どうしが一体、重複固有値なのか近接固有値なのか、判断に苦しむ場面に出くわすこともあります。
また、2番目以降の固有値を求めたい場合の使用では、シフト処理(原点移動)というテクニックを利用することになります。シフト処理の詳細については、固有値算法を説く適当な教科書を見ていただくとして、ここでは省きます。
本法を一言で評価するとすれば、「きれいな花には刺がある」というところでしょうか。物理的評価ができる経験者が注意深く使用するならまだしも、間違っても初心者がお手軽に使用するものではないと思います － つい、お手軽に使いたい簡便性の魅力を持っていますが。
[サブスペース法]
元の固有ベクトルをリッツベクトルという直交ベクトルで近似展開することから始まります。イメージ的には、リッツベクトルを一種の直交関数とみなせば、ちょうど任意関数をフーリエ級数展開で近似できるという数学処理に似てなくもないです。しかもN次元であったベクトルが、リクエストする固有値の数p個のベクトルで展開されることで、N次元マトリックスA、Bよりも格段に小さいp次元マトリックスAS、BSの固有値問題に帰着されることになります － 実際の計算では、AS、BSの次元はpよりもやや大きい数が設定されますが。
このとき、いかにA、Bが粗マトリックスおよびバンドマトリックスの状態であってもAS、BSは小さな密マトリックスになります。変換後は、小サイズになった標準固有値解析なので、それ対象の固有値算法の選択幅は広くなるのですが、ここの算法では、安定な算法として定番のヤコビ法が採用されることが多いようです。もちろん、この段階の固有値計算でMathLibを利用するのも一法です。冒頭で、本法が、固有値、固有ベクトルを同時に求めると言いましたが、この由来は実にヤコビ法にある、という舞台裏もあります。
一方、最初に仮定したリッツベクトルがそのままというわけではなく、逆反復法を利用して改善処理を繰り返すことも必要です。この意味で、本法は、同時逆反復法とも呼ばれるのですが、サブスペース法の名前の由来は、元のベクトル空間から、リッツベクトル空間という部分空間(サブスペース)に帰着させるところから来ています。
結論から言いますと、アプリケーションプログラムで採用する固有値解析の算法としては一番無難ではないか、と筆者は感じています。重複固有値、近接固有値が存在していても、安定的に解を求めてくれるからです。
欠点と言えば、初期ベクトルとして、ライバルのランチョス法が1本用意すればいいのに対して、p本(実際は＋α)用意する必要があることです。しかも、それらの初期ベクトルはお互い一次従属であってはならないので、勝手な設定は出来ません。実際、初期ベクトルの設定によっては、反復計算がうまく収束しない経験もします。アプリケーションプログラムの開発者は、うまくいかなかったときの、初期ベクトル変更機能も用意しておくべきだと思います。
余談となりますが、サブスペース法の創始者は、昔汎用有限要素法解析ソフトで有名だったSAPやADINAの開発に携わったK.J.Batheです。彼の1971年の論文で世に提案されたのが嚆矢のようです。おそらく日本のこの種のアプリケーションプログラムでは、一番多く採用されている実績のある固有値算法だと思います。
[ランチョス法]
サブスペース法の最大のライバル的算法です。大きな問題点を持つゆえからか、一度は停滞気味だったのですが、近年、市販有限要素法解析ソフトで導入されだして脚光を浴びて復活した固有値算法です。
問題点というのは、数値の扱いにデリケートなところがあり、次々生み出す近似ベクトルの計算において、発生した数値誤差がベクトル同士の直交性を崩しやすいという欠点が指摘されていたのです。Batheなどは自分の著書で、やはり、自分が開発したサブスペース法への思い入れが強いのか、この点を痛烈に批判しています。復活ならしめた理由に、反復計算過程でベクトル同士の再直交処理が組み込まれたアルゴリズムがあります。
またランチョス法の復活を押した理由は、先にも言いました1本の初期ベクトルの準備でいいという有利さもあるのですが、最大の理由は、計算速度がサブスペース法に比べて有利だという点にあると思います。有利と言っても、中規模クラスの計算では、さほどの差があるわけではなく、大規模クラスの計算においては、計算スピード面に関してランチョス法がサブスペース法に優ってきます。
ところで、固有値算法では、いくつもの変換作業がつきものですが、前編で述べた一般固有値解析から、標準固有値解析への変換作業もそのひとつです。また、固有値方程式での元の係数マトリックスを3重対角マトリックスへ変換して固有値を求めやすくする手法も主流の固有値算法の一つです。ランチョス法はこの二つの変換を一挙にやり遂げてしまう見事な算法なのです。しかも、この3重対角マトリックスの大きさが元の次元Nではなくて、Nに比べて極端に小さくなる(3重対角マトリックのサイズはユーザーが指定した数nになる)ことが本法のミソであります。
3重対角マトリックスへの変換では、ランチョス(続・有限要素法よもやま話27話参照)が生みの親なのかどうかまで分かりませんが、彼が関係したことは間違いないでしょうランチョス変換という方法が使用されます。このためこの固有値算法はランチョス法と名付けられています。ランチョス変換の理論背景は魅惑的で、そのアルゴリズムは比較的簡単で、そして結果は驚くべきものです。
反復過程で次々生み出される直交ベクトルから計算されるαi、βi（i=1～n）からなる3重対角マトリックス(図3)の固有値が元の係数マトリックスでのそれと同じ、というなんという数学的幸運さか、思わず拍手したくなりますよ。
さて、3重対角マトリックスで採られる固有値算法は、二分法が定番のようです。この固有値算法はMathLibをまるまる利用してもいいでしょうし、数学書を参考に自己開発でもいいのではないでしょうか。
ランチョス法の欠点と言えば、上で言った誤差問題を除いても、やや厄介なことが一つあります。nの決定に少し苦労するのです。ユーザーにとってnの値は自動で決まるのではなく、外から与える必要があります。要求する固有値数の精度に対して、どの程度のnを設定すればいいのか、判断に困ります。無難を念じて、大きめの数値を設定すれば、計算スピードが早いというこの算法のメリットを減殺してしまうことになり、ランチョス法を採用するという意味がなくなってしまいます。通常は、要求する固有値数と同程度以上の数値を設定することが多いでしょうが、その辺の匙加減がエンジニアの感に委ねられます。この問題点があるため、数値解析の大家戸川隼人先生も初心者には勧められないと断言されています。
筆者の体験談を追記すれば、対称構造物の振動解析のような重複固有値が生じる問題では、リクエストする固有値数よりも結構上乗せした数値をnに設定しないと、うまく固有値が求まらない、というケースがありました。
2011.11.01記
