今回は、少し数値解析の話をしよう。工学問題を解析する数値解析のプログラミングをしていると、時折、下のようなマトリックス３重積の場面に出くわすことがある。

この A、B、C がすべて大きさの違う長方形マトリックスというケースは数学での行列演算の練習問題ならば考えられるが、工学解析ではあまり考えられない。工学問題では Bが正方マトリックスであることが多いと思う。

マトリックス３重積を実際に計算する場合、前から乗算していくのと、後ろから乗算していくのでは、乗算回数に差が出ることに注意が必要である。今、A が（L, N）の、C が（N,M）のマトリックスであるとすると――もちろん、このとき B は（N, N）のマトリックスである―― L＞M のときは、後ろから、L＜M のときは、前から乗算を進めるのが乗算回数で有利となる。

L＝M の場合は、どちら側から始めても乗算回数は同じである。どちらかといえば、工学問題ではこのケースが多いので、A、B、C の各マトリックスに何か特徴がなければ、これ以上話題がなくて話が止まってしまう。

だが、その心配はご無用。特徴があるのである。工学問題の多くでは B が対称マトリックスになっている。さらに、扱っている物理量がテンソル的なものであるとき、A が C の転置マトリックスになっているということが多い。こうなると、俄然、話が続くのである。

有限要素法で頻発するこのマトリックス３重積、見方によっては B の変換式であり、これを合同変換と呼ぶことは、既に第38話で紹介済みである。

一方、対称マトリックスを三角分解する方法の１つに、コレスキー分解というアルゴリズムがある。マトリックスを下三角マトリックスと上三角マトリックスの積の形に変換する上に、２つのマトリックスが互いに転置マトリックスの関係になっている分解である。すなわち、B に適用すれば下の通りである。

Bの分解式をマトリックス３重積に代入すれば、次のように整理される。

ということは、UC というマトリックスを求めておけば、自分自身の転置マトリックスとの間の積でマトリックス３重積が求まることになる。

ここで、B を（N×N）の、C を（N×M）のマトリックスとして、元のマトリックス３重積をそのまま乗算した場合と、コレスキー分解を利用した場合で乗算回数を比較してみる。

２番目の式の最後の項は、コレスキー分解時の乗算回数である。もし、N が大きな数であるときは、三乗の項で効いてくるから、この方法は悲劇を迎える。しかし、幸いにも、多くの問題で N は小さい数値で、しかも、N＜M であることが多い。有限要素法の要素剛性マトリックスを求める場面でいえば、ソリッド要素で、N＝６である。M の方は、要素の構成節点数で変化し、低次六面体要素で、M＝24、高次六面体要素で、M＝60となる。

上の比較式を眺めると、コレスキー分解を利用した場合、次の回数だけ乗算が節約できることが分かる。

N、M に上の数字を代入してみれば分かる通り、両者に大差が生じるわけではない。しかし、マトリックス３重積がループの中にある場合は、その差が増幅される。

要素剛性マトリックスの計算では、マトリックス積が図51-1に示す位置にある。この図では、コレスキー分解は１回済ませておけばいいことも分かる。通常、コレスキー分解を利用する場合は、このように１回で済むことも多い。それで、節約回数の内の N³ / 6 の項は無視してよい。結局、要素剛性マトリックスの計算においては、乗算の節約回数が次の通りとなる。

ここで、具体的に数字を計算してみよう。仮に要素数を１万要素として、六面体要素に対して数字を入れてみると、低次要素で、INT＝８ゆえ、約2,900万回、高次要素で、INT＝27ゆえ（フル積分の場合）、約２億４千万回の節約回数となる。節約回数だけを見ていると、効果大に思えてしまうが、必須の NEL×INT×NM² という乗算回数が頑として残るので、比率で評価すると、低次要素で9.4%、高次要素で４%の乗算減数となり、ちょっと、がっかりの数値に思う人もいるかもしれない。

この例で見るとおり、工学の数値解析でマトリックス３重積に出くわした際の、コレスキー分解法の評価を下せば、既に、コレスキー分解のツールを持っている人は、是非、利用すべきである。初級プログラマーとの違いを見せつけるのが、こういうテクニックである。一方、新たに単発的に必要となった人までが、わざわざ、コレスキー分解のプログラムを作成してまで使用するほどのことでもないというところだろうか。

次に、有限要素法に現れるもう１つの代表的なマトリックス３重積として、要素剛性マトリックスの座標変換をちょっと考えてみよう。しかし、この局面で、コレスキー分解法を適用しようとしたら、途端に計算破綻を発見することになるだろう。残念ながら、この場合には適用不可なのである。その理由は、コレスキー分解そのものにある。

コレスキー分解のアルゴリズムは、計算過程において対角項の平方根計算が必要であり、それが分母位置にも来る。したがって、対角項は常に正値でないといけない。換言すれば、分解対象のマトリックスが正定値でないといけない。ところが、要素剛性マトリックスには拘束条件が導入されていないので、マトリックスが特異マトリックスとなっている（第36話参照）。これが理由で適用不可なのである。

最後に、コレスキー分解についてもう少し紹介しておこう。連立方程式の解法の１つとしてコレスキー分解が利用されることも多い。と言うよりも、むしろ連立方程式の解法のアルゴリズムとして生まれたと言っていいものだ。

係数マトリックスが三角分解さえされていれば、次のように未知変数を一旦、作業変数に置き換えさえすれば、後は簡単な連続代入の計算で済んでしまう。

このアルゴリズムは、上で言ったように、分解過程で対角項に負値になることを許さないので、非線形解析での利用が出来ないなど、やや適用限界が狭くなるが、それでも通常の数値解析分野では広く使われているはずである。

コレスキー分解を開発したコレスキー（Cholesky；仏1875－1918）はフランスの陸軍少佐だったそうです。有名なエコール・ポリテクニックの卒業生で数学に強かったので、測量部隊に配属になり、測量データの整理からコレスキー分解を考案したそうです。他書では皆無と思われるコレスキー少佐の興味深い経歴が次の本で紹介されている。

“反復法の数理” 藤野清次／張紹良著、1996年、朝倉書店。

2008年12月 記