intel Math Karnel Library (以降MKL)を用いて、一般固有値問題の対角化を行います。
特に、登場する行列が実対称行列である場合について考えます。

PDF版はこちらをどうぞ。
https://slpr.sakura.ne.jp/qp/supplement_data/diagonalize/GeneralizedEVP_symmetrixMatrix.pdf

概要

本稿で取り扱う一般固有値問題は、下記の固有値問題です。
\(\displaystyle
A\mathbf{x}=\lambda B \mathbf{x}
\)
ここで、\(A\)は\(N\times N\)の実対称行列、\(B\)は\(N\times N\)の実対称正定値行列（行列式がゼロではない）、\(\lambda, \mathbf{x}\)はこの問題の解であり、固有値、固有ベクトルを意味します。

MKLでは、下記の順で一般固有値問題の固有値・固有ベクトルを得ます。

1. 一般固有値問題を標準形式の固有値問題に縮退（変形）
2. 標準形式の固有値問題を解く

ここでいう標準形式の固有値問題とは、下記の形を持つ固有値問題です。
\(\displaystyle
A\mathbf{x}=\lambda \mathbf{x}
\)
ここで、\(A\)は\(N\times N\)の実対称行列です。

一般固有値問題から標準形式の固有値問題への式変形は、単純には両辺に左から\(B\)の逆行列\(B^{-1}\)を書ければ良いです。
ただし対称行列の逆行列は対称行列となりますが、この方法では\(B^{-1}A\)の計算を実施すると対称行列では無くなってしまうので、あまり良い方法ではありません。
また、数値計算においては逆行列をあらわに求めることはあまり意味がなく、内部的には例えばコレスキー因子分解を利用して標準形式の固有値問題に焼き直しています。
詳細は本稿のpdf、または
桂田祐史著, 『一般化固有値問題』http://nalab.mind.meiji.ac.jp/~mk/labo/text/generalized-eigenvalue-problem.pdf
をご参照ください

プログラム

MKLを使用して実対称行列を対象にする場合、?potrfと?sygstを用いて、一般固有値問題から標準形式に変形します。

変形後、通常の固有値問題となりますので、?syevdなりで得れば良いです。

下記のプログラムで一般固有値問題を解くことができます。

https://slpr.sakura.ne.jp/qp/supplement_data/diagonalize/main.f90

MKLにリンクして実行すると下記の結果を得ます。

Tweet