pianofisica

Mathematics & Physics, Maxima, a bit Python & Wolfram, and Arts

Python(SymPy)で学ぶ行列の固有値とJordan標準形

数学の具体的な計算にPython(SymPy)を使って、数学もPython(SymPy)も同時に学んでしまいましょう。今回はPython(SymPy)を使って行列の計算をしてみたいと思います。行列の固有値固有ベクトルを求め、対角化する方法、Jordan標準形を求める方法をみてみます。

本記事はPython(SymPy)で行列の扱い方を解説した過去記事

pianofisica.hatenablog.com

の続編という位置付けになりますので、あわせて読んでみてください。



固有値固有ベクトル

行列  A固有値固有ベクトルとは、行列に関する方程式

 \qquad Au=\lambda u

をみたす数  \lambda と(非自明な)ベクトル  u のことをいいます。

固有値  \lambda に対しては、以下の固有方程式

 \qquad\det(A-\lambda\,I)=0

を満たすことが要請されます。この方程式の解として、行列  A固有値が求まります。

行列の対角化

さて、 n 次正方行列  A固有値  \lambda_i \ (i=1,\dots,n) に対する固有ベクトル  u_i が線形独立である場合には、行列  P

 \qquad\displaystyle{P=\left(\begin{array}{cccc} u_1 & u_2 & \dots & u_n \end{array}\right)}

によって導入すると、この行列  P は正則であることがわかり、また  u_i が行列  A固有ベクトルであることから

 \qquad\displaystyle{AP=\left(\begin{array}{cccc} \lambda_1u_1 & \lambda_2u_2 & \dots &\lambda_n u_n \end{array}\right)}

より、これに左から行列  P^{-1} を作用させて

\qquad \displaystyle{\Lambda=P^{-1}AP=\left(\begin{array}{cccc} \lambda_1 & 0 & \dots & 0 \\  0 & \lambda_2 & \dots & 0 \\ \vdots & & & \vdots \\  0 & 0 & \dots & \lambda_n \end{array}\right)}

という対角行列を作ることができます。このような対角行列を作ることのご利益ですが、行列  A の累乗が

\qquad A^k=(P\Lambda P^{-1})^k=P\Lambda^kP^{-1}\qquad\quad \Lambda^k=\left(\begin{array}{cccc} \lambda_1^k & 0 & \dots & 0 \\  0 & \lambda_2^k & \dots & 0 \\ \vdots & & & \vdots \\  0 & 0 & \dots & \lambda_n^k \end{array}\right)

によって簡単に計算できるということがあげられます。



具体例:2次正方行列の場合

固有値固有ベクトルを具体的にみるために行列

\qquad \displaystyle{A=\left(\begin{array}{cc} 1 & 2 \\ 3 & 2 \end{array}\right)}

を考えてみます。行列をPythonに入力するには

import sympy as sp
sp.init_printing()

A = sp.Matrix([[1,2],[3,2]])

とします。

固有方程式  \det(A-r\,I)=0行列式

\qquad \displaystyle{\det\left(\begin{array}{c} 1-r & 2 \\ 3 & 2-r \end{array}\right)}

sp.var('r') 
L = r*sp.eye(2)
D = (A-L).det()

をゼロとして r について方程式を解きます:

sp.solve (D, r)

により、行列の固有値

 \qquad\displaystyle{\lambda_1=-1 \qquad \lambda_2=4}

と求まります。実際、Pythonの組み込み関数を使って

A.eigenvals()

により、固有値とその縮重度を知ることができます(いまの場合、固有値は−1と4でその縮重度はどちらも1です)。さらに

A.eigenvects()

では、固有値とその縮重度に加え、それらに対応する固有ベクトルを出力します:

 \qquad\displaystyle{u_1=\left(\begin{array}{cc} -1 \\ 1 \end{array}\right) \qquad u_2=\left(\begin{array}{cc} 2/3 \\ 1 \end{array}\right)}

そこで行列

 \qquad\displaystyle{P=\left(\begin{array}{cc} u_1 & u_2  \end{array}\right)=\left(\begin{array}{cc} -1 & 2/3 \\ 1 & 1 \end{array}\right)}

を導入して  P^{-1}AP を計算すると

P = sp.Matrix([[-1,sp.Rational(2,3)],[1,1]])
R = P.inv()*A*P

により

 \qquad\displaystyle{R=P^{-1}AP=\left(\begin{array}{cc} -1 & 0 \\ 0 & 4 \end{array}\right)}

となることがわかります。以上により行列が対角化できました。

 

行列が対角化できない場合

以上みたように  n 次正方行列に対し、線型独立な固有ベクトル n 個ある場合には行列を対角化することができますが、一般には線形独立な固有ベクトルの個数は  n より小さくなります。

例えば行列

\qquad \displaystyle{B=\left(\begin{array}{cc} 1 & 2 \\ -2 & 5 \end{array}\right)}

を考えると

B = sp.Matrix([[1,2],[-2,5]])
B.eigenvects()

より、この行列は固有値3(縮重度は2)を持ちますが、固有ベクトルが(定数倍の自由度を除いて)ただ一つしかありません:

 \qquad\displaystyle{v=\left(\begin{array}{cc} 1 \\ 1 \end{array}\right)}

このような場合には、上で議論したような正則行列が作れず、行列を対角化することができません。


 


Jordan標準形

上の例でみている行列  B は、対角化はできませんが、Jordan標準形と呼ばれる次の形

 \qquad\displaystyle{J=\left(\begin{array}{cc} 3 & 1 \\ 0 & 3 \end{array}\right)}

に、適当な正則行列  Q を使って変形することができます:

 \qquad\displaystyle{Q^{-1}BQ=J}

このようなJordan標準形と呼ばれる形をした行列を考えることのメリットとして、その累乗が

 \qquad\displaystyle{J^n=\left(\begin{array}{cc} 3^n & n\times 3^{n-1} \\ 0 & 3^n \end{array}\right)}

と計算できることがあげられます(これは数学的帰納法で簡単に示せますね)。すると、対角行列の場合と同様に、行列  B の累乗を

 \qquad\displaystyle{B^n=(QJQ^{-1})^n=QJ^nQ^{-1}}

として簡単に計算することができます。

また、Jodran標準形の対角部分と非対角部分を

 \qquad\displaystyle{J_{\rm S}=\left(\begin{array}{cc} 3 & 0 \\ 0 & 3 \end{array}\right)\qquad\qquad J_{\rm N}=\left(\begin{array}{cc} 0 & 1 \\ 0 & 0 \end{array}\right)}

と分けてみましょう( J=J_{\rm S}+J_{\rm N})。このとき対角部分の累乗が

 \qquad\displaystyle{J_{\rm S}^n=\left(\begin{array}{cc} 3^n & 0 \\ 0 & 3^n \end{array}\right)}

である一方、非対角部分は2乗するとゼロになります:

 \qquad\displaystyle{J_{\rm N}^n=0\qquad (n\geq2)}

実際

JN = sp.Matrix([[0,1],[0,0]])
JN**2

また、対角部分と非対角部分は交換します:

 \qquad\displaystyle{J_{\rm S}J_{\rm N}-J_{\rm N}J_{\rm S}=0}

具体的に確かめてみましょう。

JS = sp.Matrix([[3,0],[0,3]])
JS*JN-JN*JS

以上の性質を合わせると、行列  J の累乗が

 \qquad\displaystyle{J^n=\left(J_{\rm S}+J_{\rm N}\right)^n=J_{\rm S}^n+nJ_{\rm S}^{n-1}J_{\rm N}}

で与えられることがわかります。これは当然ながら上で与えた一般式と一致します。このような行列の分解をJordan分解といい、 J_{\rm S} を行列  J の半単純成分、 J_{\rm N} を冪零成分といいます。行列の作用の線形性から、これまでの議論を逆にたどって、行列  B の半単純成分、冪零成分がそれぞれ

 \qquad\displaystyle{B_{\rm S}=QJ_{\rm S}Q^{-1}\qquad\qquad B_{\rm N}=QJ_{\rm N}Q^{-1}}

で与えられることがわかります。

Pythonには与えられた行列のJordan標準形と変換行列を求める組み込みがあって

B.jordan_form()

とすると、行列  B をJordan標準形  J に変形する変換行列  Q

 \qquad\displaystyle{Q=\left(\begin{array}{cc} -2 & 1 \\ -2 & 0 \end{array}\right)}

と求まります。実際に  Q^{-1}BQ を計算すると

Q = sp.Matrix([[-2,1],[-2,0]])
Q.inv()*B*Q

より

 \qquad\displaystyle{Q^{-1}BQ=\left(\begin{array}{cc} 3 & 1 \\ 0 & 3 \end{array}\right)=J}

となることがわかります。以上により行列をJordan標準形に変形できました。そして行列  B の半単純成分、冪零成分がそれぞれ

BS = Q*JS*Q.inv()
BN = Q*JN*Q.inv()
[BS,BN]

より

 \qquad\displaystyle{B_{\rm S}=QJ_{\rm S}Q^{-1}=\left(\begin{array}{cc} 3 & 0 \\ 0 & 3 \end{array}\right)\qquad\qquad B_{\rm N}=QJ_{\rm N}Q^{-1}=\left(\begin{array}{cc} -2 & 2 \\ -2 & 2 \end{array}\right)}

であることがわかります。

 


キーワードPython、行列、固有値、対角化、Jordan標準形、Jordan分解、半単純成分、冪零成分