前編では、固有値解析を含むアプリケーションプログラムにおいては、既存のMathlibの全面利用ではなく、自ら構築したプログラムスキームの中に固有値算法を取り込むことを勧めたところで終わった。もっとも、マトリックスの次元Nが大きい場合の話であるが。

しかし、お勧めするとは言っても、前編で挙げたp«NかつNが大次元という課題をカバーする一般固有値解析の算法といえば、悲しいことに選択肢に限りがある。実際問題、次の三つの算法があるだけではないかと思う-浅学な筆者ゆえ、他に存在することを知らずにいるかもしれないが、その場合はお許しくだされ。いずれも初期ベクトルを設定してからスタートする反復計算型の算法である。しかし、固有値、固有ベクトルの求め方が三者三様である。

インバース・パワー法：

固有ベクトルを求めてから、固有値を求める

サブスペース法：

通常は固有値、固有ベクトルを同時に求める

ランチョス法：

固有値を求めてから固有ベクトルを求める

以下、この三つの算法の概略説明および筆者がプログラム開発で利用してきた感想を述べてみたい。

[インバース・パワー法]

最小固有値一つを求める算法で、日本語で逆反復法とも呼ばれているものである。理論といい、算法といい理解し易い方法なので、初中級レベルの数値計算プログラムの練習問題には格好の材料と言えるかもしれない。

孤立した固有値ならば、ものの見事に早く固有値を探し出す優れた算法だが、問題は重複固有値の存在、近接固有値が存在する場合である(図1)。この二つが存在する場合は、初期ベクトルの設定にもよるが、うまく収束しないことを覚悟する必要がある。

重複固有値、近接固有値というのは、例えば、3次のマトリックスを想定した場合、固有方程式を固有値変数λで展開した三次の固有多項式をイメージすれば分かりが早い。三次関数の根がすなわち固有値であるから、重複固有値とは、図1左の場合で、数学書では縮退固有値とも呼んでいる。弾性構造物の振動解析でいえば、対称性を持つ構造では、対称度数に相応した同値の、すなわち重複固有値が存在するため、このようなときには結果の評価に注意が必要である。

近接固有値の方は、同じく図1右にあるように、根の値に違いがあるにはあるのだが、お互い非常に近接した隣接固有値のことをいう。

数値計算の結果には、もちろん数値誤差を持つのが当然なので、物理問題によっては、微小な差異を持つ複数の固有値どうしが一体、重複固有値なのか近接固有値なのか、判断に苦しむ場面に出くわすこともある。

 

また、2番目以降の固有値を求めたい場合の使用では、シフト処理(原点移動)というテクニックを利用することになる。シフト処理の詳細については、固有値算法を説く適当な教科書を見ていただくとして、ここでは省く。

本法を一言で評価するとすれば、「きれいな花には刺がある」というところだろうか。物理的評価ができる経験者が注意深く使用するならまだしも、間違っても初心者がお手軽に使用するものではないと思う-つい、お手軽に使いたい簡便性の魅力を持っているのだが。

 

[サブスペース法]

元の固有ベクトルをリッツベクトルという直交ベクトルで近似展開することから始まる。イメージ的には、リッツベクトルを一種の直交関数とみなせば、ちょうど任意関数をフーリエ級数展開で近似できるという数学処理に似てなくもない。しかもN次元であったベクトルが、リクエストする固有値の数p個のベクトルで展開されることで、N次元マトリックスA、Bよりも格段に小さいp次元マトリックスA_S、B_Sの固有値問題に帰着されることになる-実際の計算では、A_S、B_Sの次元はpよりもやや大きい数が設定されるが。

このとき、いかにA、Bが粗マトリックスおよびバンドマトリックスの状態であってもA_S、B_Sは小さな密マトリックスになる。変換後は、小サイズになった標準固有値解析なので、それ対象の固有値算法の選択幅は広くなるのだが、ここの算法では、安定な算法として定番のヤコビ法が採用されることが多いようだ。もちろん、この段階の固有値計算でMathLibを利用するのも一法である。冒頭で、本法が、固有値、固有ベクトルを同時に求めると言ったが、この由来は実にヤコビ法にある、という舞台裏もある。

一方、最初に仮定したリッツベクトルがそのままというわけではなく、逆反復法を利用して改善処理を繰り返すことも必要である。この意味で、本法は、同時逆反復法とも呼ばれるのだが、サブスペース法の名前の由来は、元のベクトル空間から、リッツベクトル空間という部分空間(サブスペース)に帰着させるところから来ている。

結論から言えば、アプリケーションプログラムで採用する固有値解析

の算法としては一番無難ではないか、と筆者は感じている。重複固有値、近接固有値が存在していても、安定的に解を求めてくれるからである。

欠点と言えば、初期ベクトルとして、ライバルのランチョス法が1本用意すればいいのに対して、p本(実際は＋α)用意する必要があることである。しかも、それらの初期ベクトルはお互い一次従属であってはならないので、勝手な設定は出来ない。実際、初期ベクトルの設定によっては、反復計算がうまく収束しない経験もする。アプリケーションプログラムの開発者は、うまくいかなかったときの、初期ベクトル変更機能も用意しておくべきだと思う。

余談となるが、サブスペース法の創始者は、昔汎用有限要素法解析ソフトで有名だったSAPやADINAの開発に携わったK.J.Batheである。彼の1971年の論文で世に提案されたのが嚆矢のようだ。おそらく日本のこの種のアプリケーションプログラムでは、一番多く採用されている実績のある固有値算法だと推測する。

 

[ランチョス法]

サブスペース法の最大のライバル的算法である。大きな問題点を持つゆえからか、一度は停滞気味だったのだが、近年、市販有限要素法解析ソフトで導入されだして脚光を浴びて復活した固有値算法である。

問題点というのは、数値の扱いにデリケートなところがあり、次々生み出す近似ベクトルの計算において、発生した数値誤差がベクトル同士の直交性を崩しやすいという欠点が指摘されていたのだ。Batheなどは自分の著書で、やはり、自分が開発したサブスペース法への思い入れが強いのか、この点を痛烈に批判している。復活ならしめた理由に、反復計算過程でベクトル同士の再直交処理が組み込まれたアルゴリズムにある。

またランチョス法の復活を押した理由は、先にも言った1本の初期ベクトルの準備でいいという有利さもあるのだが、最大の理由は、計算速度がサブスペース法に比べて有利だという点にあると思う。有利と言っても、中規模クラスの計算では、さほどの差があるわけではなく、大規模クラスの計算において、計算スピード面に関してランチョス法がサブスペース法に優ってくるようだ。

ところで、固有値算法では、いくつもの変換作業がつきものだが、前編で述べた一般固有値解析から、標準固有値解析への変換作業もそのひとつである。また、固有値方程式での元の係数マトリックスを3重対角マトリックスへ変換して固有値を求めやすくする手法も主流の固有値算法の一つである。ランチョス法はこの二つの変換を一挙にやり遂げてしまう見事な算法なのだ。しかも、この3重対角マトリックスの大きさが元の次元Nではなくて、Nに比べて極端に小さくなる(3重対角マトリックのサイズはユーザーが指定した数nになる)ことが本法のミソである。

3重対角マトリックスへの変換では、ランチョス(本エッセイ41話参照)が生みの親なのかどうかまで分からないが、彼が関係したことは間違いないだろうランチョス変換という方法が使用されている。このためこの固有値算法はランチョス法と名付けられている。ランチョス変換の理論背景は魅惑的で、そのアルゴリズムは比較的簡単で、そして結果は驚くべきものである。

反復過程で次々生み出される直交ベクトルから計算されるαi、βi（i=1～n）からなる3重対角マトリックス(図3)の固有値が元の係数マトリックスでのそれと同じ、というなんという数学的幸運さか、思わず拍手したくなる。

さて、3重対角マトリックスで採られる固有値算法は、二分法が定番のようだ。この固有値算法はMathLibをまるまる利用してもいいし、数学書を参考に自己開発でもいいのではないだろうか。

ランチョス法の欠点と言えば、上で言った誤差問題を除いても、やや厄介なことが一つある。nの決定に少し苦労するのである。ユーザーにとってnの値は自動で決まるのではなく、外から与える必要がある。要求する固有値数の精度に対して、どの程度のnを設定すればいいのか、判断に困ってしまう。無難を念じて、大きめの数値を設定すれば、計算スピードが早いというこの算法のメリットを減殺してしまうことになり、ランチョス法を採用するという意味がなくなってしまう。通常は、要求する固有値数と同程度以上の数値を設定することが多いだろうが、その辺の匙加減がエンジニアの感に委ねられることになる。この問題点があるため、数値解析の大家戸川隼人先生も初心者には勧められないと断言されている。

筆者の体験談を追記すれば、対称構造物の振動解析のような重複固有値が生じる問題では、リクエストする固有値数よりも結構上乗せした数値をnに設定しないと、うまく固有値が求まらないというケースがあった。

2011年11月記