どうも、木村（@kimu3_slime）です。
今回は、Julia（SymPy）で固有値、固有ベクトル、対角化、ジョルダン標準形を求める方法を紹介します。
準備
SymPyを使うので、持っていなければインストールしておきましょう。
1 2 | using Pkg Pkg.add("SymPy") |
準備として、以下のコードを実行しておきます。
1 | using SymPy |
固有値、固有ベクトル、対角化、ジョルダン標準形を求める方法
簡単な例
1 | B = Sym[1 2; 2 4] |
1 | B.eigenvals() |
この行列の固有値と、その代数的重複度を求めましょう。「行列.eigenvals()」です。
1 2 3 | Dict{Any, Any} with 2 entries: 5 => 1 0 => 1 |
辞書型で結果が返ってきます。固有値5の代数的重複度は1、固有値0の代数的重複度は0という意味です。
中身を参照したいのですが、
1 | B.eigenvals()[5] |
1 | KeyError: key 5 not found |
なぜかキーが存在しないと言われてしまいます。「keys(B.eigenvals())」で確認すると含まれていますし、「5」の型（Int64）はAny型として扱われるはずなのですが……。SymPy側とJuliaのかみ合わせで何か問題があるのかもしれません。
キーとしてSym型を指定することで、
1 | B.eigenvals()[Sym(5)] |
\[ \begin{aligned}5\end{aligned} \]
対応する値が返ってきます。辞書型は「Dict{Any, Any}」と書いてあるのですが、実態としては「Dict{Sym,Any}」ということなのでしょう。
次の方法はすっきりしています。「行列.eigenvects()」で、固有値、代数的重複度、固有ベクトルをセットで求めてくれます。
1 | B.eigenvects() |
これらが固有値、固有ベクトルであることを計算でチェックしてみましょう。eigenvects()の入れ子構造に注意。
1 2 | B.eigenvects()[1][3][1] B*B.eigenvects()[1][3][1] |
\[ \begin{aligned}\left[ \begin{array}{r}-2\\1\end{array}\right]\end{aligned} \]
\[ \begin{aligned}\left[ \begin{array}{r}0\\0\end{array}\right]\end{aligned} \]
1 2 3 | B.eigenvects()[2][3][1] B*B.eigenvects()[2][3][1] B*B.eigenvects()[2][3][1] == 5*B.eigenvects()[2][3][1] |
\[ \begin{aligned}\left[ \begin{array}{r}\frac{1}{2}\\1\end{array}\right]\end{aligned} \]
\[ \begin{aligned}\left[ \begin{array}{r}\frac{5}{2}\\5\end{array}\right]\end{aligned} \]
それぞれ、固有値0,5に対応する固有ベクトルであることがわかりました。
「行列.diagonalize()」で、行列が対角化可能ならば、変換行列、対角行列を返します。
1 | B.diagonalize() |
1 | (Sym[-2 1; 1 2], Sym[0 0; 0 5]) |
第1成分を\(P\)（変換行列）、第2成分を\(D\)（対角行列）とすると、\(PDP^{-1}= A\)となるような関係です。実際にこれを計算してみれば、
1 | B.diagonalize()[1] * B.diagonalize()[2] * B.diagonalize()[1]^-1 |
\[ \begin{aligned}\left[ \begin{array}{rr}1&2\\2&4\end{array}\right]\end{aligned} \]
となって、確かに対角化できていることがわかります。
ここまで固有値を求めてきましたが、固有多項式を求めることもできます。
1 2 | B.charpoly() factor(B.charpoly().as_expr()) |
\[ \begin{aligned}\operatorname{PurePoly}{\left( \lambda^{2} – 5 \lambda, \lambda, domain=\mathbb{Z} \right)}\end{aligned} \]
\[ \begin{aligned}\lambda \left(\lambda – 5\right)\end{aligned} \]
回転行列
文字を含む行列でも、同様の方法が通用します。回転行列を例にしましょう。
1 2 3 4 5 6 | @vars x y z n α β λ R_α= Sym[cos(α) -sin(α); sin(α) cos(α)] R_α.eigenvects() R_α.diagonalize()[1] R_α.diagonalize()[2] sympy.simplify(R_α.diagonalize()[1] * R_α.diagonalize()[2]* R_α.diagonalize()[1]^-1) |
\[ \begin{aligned}\left[ \begin{array}{rr}\cos{\left(α \right)}&- \sin{\left(α \right)}\\\sin{\left(α \right)}&\cos{\left(α \right)}\end{array}\right]\end{aligned} \]
\[ \begin{aligned}\left[ \begin{array}{rr}\frac{\sin{\left(α \right)}}{\sqrt{- \sin^{2}{\left(α \right)}}}&- \frac{\sin{\left(α \right)}}{\sqrt{- \sin^{2}{\left(α \right)}}}\\1&1\end{array}\right]\end{aligned} \]
\[ \begin{aligned}\left[ \begin{array}{rr}- \sqrt{\left(\cos{\left(α \right)} – 1\right) \left(\cos{\left(α \right)} + 1\right)} + \cos{\left(α \right)}&0\\0&\sqrt{\left(\cos{\left(α \right)} – 1\right) \left(\cos{\left(α \right)} + 1\right)} + \cos{\left(α \right)}\end{array}\right]\end{aligned} \]
\[ \begin{aligned}\left[\begin{smallmatrix}\cos{\left(α \right)} & – \sin{\left(α \right)}\\\sin{\left(α \right)} & \cos{\left(α \right)}\end{smallmatrix}\right]\end{aligned} \]
回転行列の固有値、固有ベクトルは複素数です。\(\alpha = \frac{\pi}{2}\)のとき、きちんと計算できているか調べてみましょう。
1 2 3 | R_α.subs(α,PI/2) R_α.subs(α,PI/2).eigenvects() R_α.subs(α,PI/2).diagonalize() |
\[ \begin{aligned}\left[ \begin{array}{rr}0&-1\\1&0\end{array}\right]\end{aligned} \]
Iは虚数単位であり、虚部を持つ場合でも正しく計算できていることがわかります。
代数的重複度＝幾何学的重複度の例
固有値に重複がある例を考えてみましょう。
1 | C = Sym[2 1 1; 1 2 1; 1 1 2] |
\[ \begin{aligned}\left[ \begin{array}{rrr}2&1&1\\1&2&1\\1&1&2\end{array}\right]\end{aligned} \]
1 | C.eigenvects() |
固有値1の代数的重複度は2です。そしてその固有値に対して、固有ベクトル（固有空間の基底）が2つ見つかっています。
固有空間の次元（幾何学的重複度）は、固有空間の基底の個数を数えれば良いです。形式的に求めるなら、次のようにすれば良いでしょう。
1 2 | length(C.eigenvects()[1][3]) C.eigenvects()[1][2] == length(C.eigenvects()[1][3]) |
\[ \begin{aligned}2\end{aligned} \]
true
固有値1に対応する固有空間の次元は2であり、代数的重複度と一致しています。すべての固有値に対して、代数的重複度と幾何学的重複度が一致するとき、行列は対角化可能なのでした。
実際、この例は対角化可能です。
1 | C.diagonalize() |
代数的重複度≠幾何学的重複度の例
対角化できないような例を考えてみましょう。
1 2 | D = Rational(1,2)*Sym[13 -1; 9 7] D.eigenvects() |
\[ \begin{aligned}\left[ \begin{array}{rr}\frac{13}{2}&- \frac{1}{2}\\\frac{9}{2}&\frac{7}{2}\end{array}\right]\end{aligned} \]
1 2 | 1-element Vector{Tuple{Sym, Int64, Vector{Matrix{Sym}}}}: (5, 2, [[1/3; 1;;]]) |
固有値5の代数的重複度は2ですが、固有空間の次元は1になっています。これらが一致しないということは、対角化不可能です。
対角化を試してみると、
1 | D.diagonalize() |
1 2 | PyError ($(Expr(:escape, :(ccall(#= C:\Users\slime\.julia\packages\PyCall\3fwVL\src\pyfncall.jl:43 =# @pysym(:PyObject_Call), PyPtr, (PyPtr, PyPtr, PyPtr), o, pyargsptr, kw))))) <class 'sympy.matrices.common.MatrixError'> MatrixError('Matrix is not diagonalizable') |
エラーが返ってきます。「行列.is_diagonalizable()」で対角化可能かどうか判定できます。
1 | D.is_diagonalizable() |
対角化できなくても、それに近い形に変形することはできます。すなわち、ジョルダン標準形を求めてみましょう。「行列.jordan_form()」です。
1 | D.jordan_form() |
第1成分を\(P\)（変換行列）、第2成分を\(J\)（ジョルダン標準形）とするとき、\(PJP^{-1}=A\)となっています。実際に計算してみると、
1 | D.jordan_form()[1] * D.jordan_form()[2] * D.jordan_form()[1]^-1 |
\[ \begin{aligned}\left[ \begin{array}{rr}\frac{13}{2}&- \frac{1}{2}\\\frac{9}{2}&\frac{7}{2}\end{array}\right]\end{aligned} \]
と確かにジョルダン標準形が求められていることがわかりました。
以上、Julia（SymPy）で固有値、固有ベクトル、対角化、ジョルダン標準形を求める方法を紹介してきました。
今回は固有値が簡単になるような例を選びましたが、それでも対角化、ジョルダン標準形を求める計算は、手計算だと計算ミスすることもあります。コンピュータで簡単に求められると嬉しいですね。
ただし、3次の行列だとSymPyだとうまく求められない（対称行列の固有値が実数にならない）ケースがあります。その場合は、LinearAlgebra.jlなど、数値的に解く必要があるでしょう。
木村すらいむ（@kimu3_slime）でした。ではでは。
コロナ社 (2020-03-26T00:00:01Z)
￥7,353 (コレクター商品)
こちらもおすすめ
行列の対角化可能性の定義とメリット、例、同値条件について解説
なぜ行列式を学ぶのか？ 固有値・固有ベクトルの求め方：固有多項式の定義
固有空間の求め方、代数的・幾何学的重複度とは：部分空間となることの証明
対称行列の性質：内積による特徴づけ、逆行列、固有値、対角化について
