皆さん、こんにちは。
今回は、エルミート行列に限れば全ての固有値と固有ベクトルを同時に求められる数値解析の手法「ヤコビ法」を紹介します。
※連立1次方程式を解くアルゴリズムにも「ヤコビ法」という手法がありましたが連立1次方程式の数値解析 ～ガウスの消去法、ヤコビ法、ガウス・ザイデル法～ - ちょぴん先生の数学部屋 (hatenablog.com)、それとは全く別物です。
これと区別するために、今回紹介する固有値解析の手法は「ヤコビの対角化法」と呼ばれることもあります。
1. 相似変換
本題に入る前に、これから複数回にわたって紹介する固有値問題の数値解析を考えるにあたって不可欠な概念である「相似変換」について説明しておきます。
相似変換とは、次のように行列Aを正則行列Pを使って行列Bに変換することです。
①の右辺は、Pの逆行列とP自身でAをサンドイッチすることでAと相似な行列Bを作っています。
これは行列Aを対角化するときの計算式と全く同じであり、行列の対角化は相似変換の特別な場合と言えます。
行列の対角化（詳細版） - ちょぴん先生の数学部屋 (hatenablog.com)
相似な２つの行列AとBは、次の値が一致します。
他にも一致するものはたくさんあるのですが、特に重要なのは「固有値が一致する」という事実です。
早速証明していきましょう。
示すに当たっては、次の行列式の性質を使います。
行列の積の行列式は、それぞれの行列の行列式の積と等しい。そして、逆行列の行列式は、元の逆行列の逆数になる。というものです。
2×2行列の場合について証明すると次のようになります。
これらの事実を使って、まずは相似な行列A,Bの行列式が一致することを証明すると、
という感じになります。
次に固有値が一致することを示すのですが、それには元となる「固有方程式」が一致することが示せれば十分です。
ということで、Bの固有方程式を変形するとAの固有方程式になることを示していきます。
このように、行列の固有値は相似変換をしても不変です。この性質を使って、固有値を調べたい行列Aに対して相似変換を繰り返して都合のいい形の行列にする、というのが、この後紹介していくヤコビ法、QR法、ランチョス法、ハウスホルダー変換の共通した根本の考え方になります。
但し、固有値は不変ですが残念ながら固有ベクトルの方は変換のたびに変化してしまいます。なので、固有ベクトルの方は一般には固有値を求めた後元の行列Aに立ち返って改めて求める必要があります。
そんな中、今回紹介するヤコビ法では行列Aを相似変換で最終的に「対角行列」そのものに変形するので、幸いにして固有ベクトルも同時に求めることができます。
2. ヤコビ法のコンセプト
それでは、ヤコビ法の概略を説明します。
最初に断っておきますが、ヤコビ法が使えるのはAがエルミート行列の場合だけです。
エルミート行列とユニタリ行列 ～対角化がより美しくなる特別な行列～ - ちょぴん先生の数学部屋 (hatenablog.com)
物理の問題を扱う際は多くの行列はエルミート行列になってくれるのであまり心配する必要はありませんが、一応注意です。
話を元に戻すと、ヤコビ法は、Aの絶対値が最大の非対角成分を0にするようなユニタリ行列を使ってAを何度も相似変換して最終的に対称行列Λにしていくアルゴリズムになります。
ユニタリ行列を使った相似変換は「ユニタリ変換」と呼ばれるので、以後この用語を使用します。
Aをユニタリ変換してA1が得られ、
A1をユニタリ変換してA2が得られ、
この工程をn回繰り返した結果対称行列Λが出来上がったとすると、
これまでの変換式を全て代入してまとめると、
となります。波線部分の行列を
とまとめると、⑬はAの対角化の計算式そのものになっています。
一般にユニタリ行列の積もユニタリ行列になるので、⑭のPもユニタリ行列です。
（証明は↓）
そして、エルミート行列Aは、正規直交化した固有ベクトルを並べてできるユニタリ行列を使って対角化できました。
ということは、⑭のPがAを対角化するユニタリ行列そのものなので、Pの各行がAの固有ベクトルとなります。固有値は言うまでもなく対角行列Λの対角成分です。
このようにして、ユニタリ変換を繰り返すことで対角化し固有値と固有ベクトルを同時に計算するのがヤコビ法となります。
※各工程でユニタリ行列を使うのは、相似変換の際に用いる逆行列がユニタリ行列の場合は転置するだけで求まるので圧倒的に計算コストが小さくできるからです。
3. ギブンス回転
以上がヤコビ法の全体の流れなのですが、では、各工程のユニタリ行列をどうやって求めるのか？次にこの話をしていきます。
概要の部分で、ヤコビ法で使うユニタリ変換は「Aの絶対値が最大の非対角成分を0にする」と説明しました。
これを繰り返すことでAの対角成分以外を全部0にして対角行列にしようという魂胆です。
では、どうやって非対角成分を0にするのか？ヤコビ法では、その方法として「回転行列」を使います。
具体的には、次のように単位行列の一部を回転行列に変えたものを変換に使うユニタリ行列にします。
このような回転行列を使ったユニタリ変換は「ギブンス回転」と呼ばれます。
さて、このユニタリ変換でAがどのように変化するかを見ていきましょう。
Aからs行s列目、s行t列目、t行s列目、t行目t列目だけ抜き出してできる2×2行列をCとします。
Cの非対角成分をユニタリ変換で0にすることが目標です。
今、ギブンス回転を作るのに使った回転行列R(θ)を
として、R(θ)による相似変換でCが
Bになったとします。回転角θをうまく選んでBの非対角成分が0になって欲しいわけです。
実際にBを計算すると、
となります。ここでAはエルミート行列なのでCもエルミート行列となり、非対角成分は互いに等しくなるためpと置き直しています。
（※Aの成分は全て実数だとして考えています）
ここで、Bの非対角成分を
とすると、㉑の結果から
と三角関数の和で表現できます。さらに、
と各係数をまとめると、q=0となるθの条件は、
となります。これをθについて解けば、tanの逆関数を用いて
となります。※対角成分の値が等しく中身の分母が0になる場合はθ=π/4とします。
㉖で計算されるθを使ったギブンス回転を行うことで、Aのs行t列目（とt行s列目)の成分を0に変えることができます。
Aのユニタリ変換のたびに、その都度ギブンス回転の行列を作って、それを使ってユニタリ変換します。
ギブンス回転を使うと確かに狙った成分を0に変えることはできるのですが、元々0だった別の成分が0でない値に変わってしまうことがあります。なので、ある成分を0にしたら別の成分が0でなくなって・・・というモグラ叩きになってしまいがちで、対角行列に至るまでにユニタリ変換の回数が嵩んでしまうのがヤコビ法の欠点です。
また、非対角成分のすべて（n×nのエルミート行列の場合、実質n(n-1)/2個の成分）を0にする必要があるため、行列のサイズnが大きくなるほどn^2のオーダーで0にすべき非対角成分が増えていくためそれだけ計算時間が爆増します。
なので、ヤコビ法は巨大なサイズの行列には不向きで、n=30程度以下の小サイズの行列の固有値解析に使用される傾向があります。
4. Pythonでの実装例
以上説明してきたヤコビ法のアルゴリズムをPythonで記述すると次のようになります。
非対角成分を完全に0にしようとすると計算量が莫大になってしまうので、ある程度の小さい値になったところで妥協しています。
実際に、前回も登場した次の行列
でプログラムを実行すると、
9回のユニタリ変換で、全ての固有値と固有ベクトルが正しく計算できています。
先ほどヤコビ法は「巨大なサイズの行列には不向き」と説明しました。しかし、有限要素法など何万ものサイズになる行列に対して固有値を計算したい場面が多々あります。
次回は、そんな巨大な行列に対しても固有値を計算できる「QR法」を紹介します。
