おっぱいそん!

pythonを数値計算に使うときのテクニックとかをまとめていきたい。

2つの行列が相似であることの確認

2つの正方行列A,Bが与えられた時に、この2つの行列が相似、つまり、ある正則行列P(Pはユニタリ行列でなくてもよい)が存在して、
A = P B P^(-1)
と書けるかどうかを、数値的に確かめたい。
(結論から言うと、僕はどうすればいいか分かりません。いい解決法を知っている方は教えて下さい)

解析的に

解析的には、2つの正方行列A, Bのジョルダン標準系を計算して、比べれば良い。
これは、「A, Bのジョルダン標準形が一致すること(ただし、ジョルダンブロックの並べ方の不定性があるので、適当に並べ替える)」と「A, Bが相似であること」は同値であるからである。
(証明は、ジョルダン標準形の一意性から簡単に示せるので略。)
(ちなみに、ユニタリ同値であるかどうかを確かめたい時にも、ジョルダン標準形が使え、ジョルダン標準形が一致することに加えて、A, Bをジョルダン標準形にする正則行列P_A*P_B^(-1)がユニタリーがどうかを確かめればいい。)

数値的に

しかしながら、ジョルダン標準形と数値計算は相性が悪い。
というのも、ある考えている行列に縮退した固有値があり、固有ベクトルが1次独立でなければ、対応するジョルダンブロックの非対角成分に1が現れるが、数値計算においてとても小さい誤差が入り、固有値の縮退がとけると、その瞬間に固有ベクトルは1次独立となり、非対角成分がゼロになってしまうからである。つまり、数値計算における小さい誤差で、ジョルダン標準形の値が急激に変わるので、ジョルダン標準形は数値計算的に安定でないのだ。
参考:http://en.wikipedia.org/wiki/Jordan_normal_form#Numerical_analysis
実際、ジョルダン標準形は数値的に安定でないため、numpyなどに関数は用意されていません。代数計算をしてくれるSympyにはジョルダン標準形を求めてくれる関数があるみたいです。
参考:http://stackoverflow.com/questions/20312216/compute-jordan-normal-form-of-matrix-in-python-numpy


そこで、数値計算においては他の方法を使いたい。

まず、特性多項式(charactaristic polynominal)であるが、これは相似であることの必要条件であるが、十分条件でない。つまり、相似であれば、特性多項式が一致する(証明はdetの性質より簡単)ことは言えるが、特性多項式が一致したからといって、相似であるとは限らないので、使えない。

http://math.stackexchange.com/questions/14075/how-do-i-tell-if-matrices-are-similar
上の知恵袋の記事や、Wikipediaの記事によると、Schur decompositionなるものが数値的には安定で使えると書かれているが、Schur decompositionは分解した際の、上三角行列が一意的でない(http://ja.wikipedia.org/wiki/%E3%82%B7%E3%83%A5%E3%83%BC%E3%83%AB%E5%88%86%E8%A7%A3)ので十分条件にしかなっていないと思う。つまり、2つの行列A,BをSchur decompositionをした時に、上三角行列が一致していれば相似であるが、逆は言えない。この上三角行列をさらに何らかの一意的な標準系みたいなものに持っていければ、使えるのだが...
参考:http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.linalg.schur.html

行列の標準形にはジョルダン標準形以外のものもあって、Frobenius normal formやWeyr canonical formというものもあるらしいのですが(Wikipediaによると、標準形は一意的なので使えそう)、数値的に安定かどうかは知りません。また、numpyなどに関数が用意されていないようです。
参考:http://en.wikipedia.org/wiki/Canonical_form#Linear_algebra

ということで、行列の相似を数値的に確かめる方法を調べた結果、諦めました…
いい方法を知っている方は教えて下さい。