行列のべき乗を求めたい。

行列の対角化を学習すると、必ず出てくるのが行列のべき乗を求める演習問題です。
今回はこの定番の演習問題をPythonを使って解いてみたいと思います。

まずは`np.linalg`を使って直接求めてみる。

今回はAとして以下の様な行列をとることにします。

import numpy as np
A = np.array([[0.,1.,0.],[0.,0.,1.],[-2.,1.,2.]])

np.linalgを使って、例えばAの5乗を求めたければ以下の様に書けば求まります。

print(np.linalg.matrix_power(A,5))

実際に上のコードを実行すると以下の行列が表示されます。これがAの5乗ということになります。
[[-10. 1. 10.]
[-20. 0. 21.]
[-42. 1. 42.]]

余談ですが、以下の様に書くとnp.linalgの機能一覧が表示されます。

import numpy as np
help(np.linalg)

対角化してみる。

では次に一度対角してからAのべき乗を求めてみたいと思います。
以下の様に書くと、Aの固有値と、固有ベクトルを並べた行列が一気に求まります。

print(np.linalg.eig(A))

今回の場合は、固有値は2,1,-1で。固有ベクトルを並べた行列は
array([[-0.21821789, 0.57735027, -0.57735027],
[-0.43643578, 0.57735027, 0.57735027],
[-0.87287156, 0.57735027, -0.57735027]])
となります。
これをPとおくことにします。

P = np.array([[-0.21821789,  0.57735027, -0.57735027],
              [-0.43643578,  0.57735027,  0.57735027],
              [-0.87287156,  0.57735027, -0.57735027]])

次にPの逆行列を求めます。

P_inv = np.linalg.inv(P)

あとはP_inv,A,Pを順に掛けていきます。

X = np.dot(P_inv,A)
B = np.dot(X,P)

すると、BはAを対角化した行列、すなわち対角成分が固有値の2,1,-1である対角行列となることが確認できます。

行列Aのべき乗を求める。

Bのn乗は、対角成分が(2のn乗),(1のn乗(つまり1)),((-1)のn乗)である対角行列です。今はAの5乗を求めたいのでn=5とすればBの5乗は対角成分が32,1,-1の対角行列です。
よってあとはこの対角行列にPとPの逆行列を掛けていけばOKということになります。

B_5 = np.array([[32.,0.,0.],[0.,1.,0.],[0.,0.,-1.]])
Y = np.dot(P,B_5)
A_5 = np.dot(Y,P_inv)
# Aの5乗を表示。
print(A_5)

実際にこのコードを実行すると、下のような結果が表示されます。
[[-1.00000000e+01 1.00000000e+00 1.00000000e+01]
[-2.00000000e+01 4.70962477e-16 2.10000000e+01]
[-4.20000000e+01 1.00000000e+00 4.20000000e+01]]

まとめ

まとめ

import numpy as np
A = np.array([[0.,1.,0.],[0.,0.,1.],[-2.,1.,2.]])

# Aの5乗を計算。
print(np.linalg.matrix_power(A,5))

# Aの固有値と固有ベクトルを計算。

B_5 = np.array([[32.,0.,0.],[0.,1.,0.],[0.,0.,-1.]])
Y = np.dot(P,B_5)
A_5 = np.dot(Y,P_inv)
# Aの5乗を表示。
print(A_5)

因みに以下の様にすると計算時間を測れますが、3次正方行列くらいだとほとんど変わらないですね。(´-ω-`;)

start=time.time()
print(np.linalg.matrix_power(A,1000))
float(time.time())-start

start=time.time()
P = np.array([[-0.21821789,  0.57735027, -0.57735027],
              [-0.43643578,  0.57735027,  0.57735027],
              [-0.87287156,  0.57735027, -0.57735027]])
P_inv = np.linalg.inv(P)
B_5 = np.array([[np.power(2,1000),0.,0.],[0.,1.,0.],[0.,0.,np.power(-1,1000)]])
Y = np.dot(P,B_5)
A_5 = np.dot(Y,P_inv)
print(A_5)
float(time.time())-start