Aが対称行列であった場合、直交行列Uを用いて、
\mathbf{A} = \mathbf{U} \mathbf{D} \mathbf{U}^T
のように対角化することができます(Dは対角行列)。このとき、
\mathbf{A}^n = \mathbf{U} \mathbf{D}^n \mathbf{U}^T
が成立するので、Aのn乗の計算がとても容易になります。
この度、その対角化の計算(固有値分解)を容易に行うことができるpythonコードを作成しました。
func.py
def Eigenvalue_decomp(A):
# 対称行列Aを、直交行列Uを用いて対角化する。
import numpy as np
if A.shape[0]!=A.shape[1]: #Aが正方行列でない場合エラーを返す
raise Exception("error! A is not square matrix!")
if (A==A.T).all()==False: #Aが対称行列でない場合エラーを返す
raise Exception("error! A is not symmetric matrix!")
la, U = np.linalg.eig(A) # eigen value and eigen vector is resolved.
U,_=np.linalg.qr(U) # Gram-Schmidt orthonormalization
D=np.diag(la)
return U,D
これで、
In[1]: import func
In[2]: U,D=func.Eigenvalue_decomp(A)
と実行すると、numpyで作成された対称行列Aを対角化するための直交行列U、及び対角行列Dを得ることができます。
テスト
2次行列
\mathbf{A}=
\begin{pmatrix}
5.0 & 2.0 \\
2.0 & 7.0
\end{pmatrix}
In[3]: A=np.array([[5.0,2.0],[2.0,7.0]])
In[4]: U,D=func.Eigenvalue_decomp(A)
In[5]: print U.dot(D.dot(U.T)) #A=U*D*Utを計算
[[ 5. 2.]
[ 2. 7.]]
In[6]: print A.dot(A) #A*Aを計算
[[ 29. 24.]
[ 24. 53.]]
In[7]: print U.dot(D.dot(D.dot(U.T))) #A*A=U*D*D*Utを計算
[[ 29. 24.]
[ 24. 53.]]
3次行列
\mathbf{A}=
\begin{pmatrix}
1.0 & \sqrt{2.0} & 0 \\
\sqrt{2.0} & 1.0 & \sqrt{2.0} \\
0 & \sqrt{2.0} & 1.0
\end{pmatrix}
In[8]: A=np.array([[1.0,np.sqrt(2.0),0],[np.sqrt(2.0),1.0,np.sqrt(2.0)],[0,np.sqrt(2.0),1.0]])
In[9]: U,D=func.Eigenvalue_decomp(A)
In[10]: print U.dot(D.dot(U.T)) #A=U*D*Utを計算
[[ 1.00000000e+00 1.41421356e+00 -3.33066907e-16]
[ 1.41421356e+00 1.00000000e+00 1.41421356e+00]
[ -3.33066907e-16 1.41421356e+00 1.00000000e+00]]
In[11]: print A.dot(A) #A*Aを計算
[[ 3. 2.82842712 2. ]
[ 2.82842712 5. 2.82842712]
[ 2. 2.82842712 3. ]]
In[12]: print U.dot(D.dot(D.dot(U.T))) #A*A=U*D*D*Utを計算
[[ 3. 2.82842712 2. ]
[ 2.82842712 5. 2.82842712]
[ 2. 2.82842712 3. ]]
3次行列(固有値に重解が含まれる場合)
\mathbf{A}=
\begin{pmatrix}
0 & -\sqrt{2.0} & \sqrt{2.0} \\
-\sqrt{2.0} & 1.0 & 1.0 \\
\sqrt{2.0} & 1.0 & 1.0
\end{pmatrix}
この行列Aは、解析的に解くと固有値に重解(λ=2)が含まれます。
In[13]: A=np.array([[0,-np.sqrt(2.0),np.sqrt(2.0)],[-np.sqrt(2.0),1.0,1.0],[np.sqrt(2.0),1.0,1.0]])
In[14]: U,D=func.Eigenvalue_decomp(A)
In[15]: print U.dot(D.dot(U.T)) #A=U*D*Utを計算
[[ -2.85835307e-16 -1.41421356e+00 1.41421356e+00]
[ -1.41421356e+00 1.00000000e+00 1.00000000e+00]
[ 1.41421356e+00 1.00000000e+00 1.00000000e+00]]
In[16]: print A.dot(A) #A*Aを計算
[[ 4.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 4.00000000e+00 -4.44089210e-16]
[ 0.00000000e+00 -4.44089210e-16 4.00000000e+00]]
In[17]: print U.dot(D.dot(D.dot(U.T))) #A*A=U*D*D*Utを計算
[[ 4.00000000e+00 3.16592776e-16 -5.71585644e-16]
[ 3.16592776e-16 4.00000000e+00 -4.44089210e-16]
[ -5.71585644e-16 -4.44089210e-16 4.00000000e+00]]
いずれのケースでも、UとDの計算を正しく行うことができているようです。
Pythonのnumpyには、もともとnumpy.linalg.matrix_power関数という行列のn乗を計算することができる関数があるので、「対称行列のn乗を計算する」という目的で上記の関数を利用することはないでしょうが、対称行列の固有値分解を行いたい時にはとても有効であると思います。
参考URL：
http://mathtrain.jp/symmetriceigen
http://ksmzn.hatenablog.com/entry/2014/02/16/004921
http://takegue.hatenablog.com/entry/2014/07/26/222141
http://docs.scipy.org/doc/numpy-1.10.0/reference/generated/numpy.linalg.matrix_power.html