行列を対角化してべき乗を求める。ついでに数列の一般項も求める。

14426 ワード
Python Python3 sympy 行列 Python テキストリンク

初めに

行列のべき乗を求めたい時が最近よくあって、その度に毎回調べていたので備忘録的に残す。
例題としてフィボナッチ数列の一般項の導出などする。

行列の対角化って何？(数式)

行列の対角化とは以下のように対角成分のみを残して他を0にすること。

A = 
\begin{pmatrix}
a & b \\
c & d
\end{pmatrix}
→
\begin{pmatrix}
e & 0 \\
0 & f
\end{pmatrix}
\\こんなの

行列のN乗計算するの面倒だなぁ～
スカラー倍になってくれたら楽なのになぁ～
というモチベーションで生み出されたと理解してるので、そんな方針で数式化していく。

A \vec{x} =\begin{pmatrix}
1 & 2 \\
4 & 3
\end{pmatrix}
\begin{pmatrix}
x \\
y
\end{pmatrix}

という数式がありAをxに左からかけるという操作がλをxにかけるという操作と等しいラムダが欲しいので以下のような式を考える。Eは単位行列

A \vec{x} = λE \vec{x}

これを移項して

\begin{align}
A \vec{x} - λE \vec{x} = 0 \\
= (A - λE) \vec{x} \\ \\
= \begin{pmatrix}
1-λ & 2 \\
4 & 3-λ
\end{pmatrix}
\begin{pmatrix}
x \\
y
\end{pmatrix}
= \begin{pmatrix}
0 \\
0
\end{pmatrix}
\end{align}

この時、xとyが0だと何も嬉しくないのでxとyが0ではない時を考える。
(A-λE)が線形従属であればxとyの取り方によって計算結果を0にすることが可能である。
例えば以下

\begin{align}
\begin{pmatrix}
2 & 2 \\
4 & 4
\end{pmatrix}
\begin{pmatrix}
x \\
y
\end{pmatrix}
= \begin{pmatrix}
0 \\
0
\end{pmatrix}\\
x=1\\
y=-1\\
とすれば成り立つ
\end{align}

で、
線形従属ならば行列式の結果が0になるので、det(A-λE)=0になるλを調べる

\begin{align}
\begin{vmatrix}
1-λ & 2 \\
4 & 3-λ
\end{vmatrix} 
= λ^2 - 4λ - 5 = (λ + 1)(λ - 5) = 0 \\
λ = -1, 5
\end{align}

次に各λの時のxとyを計算する。

\begin{align}
λ = -1のとき\\ \begin{pmatrix}
2 & 2\\
4 & 4
\end{pmatrix}
\begin{pmatrix}
x\\
y
\end{pmatrix} = 
\begin{pmatrix}
0\\
0
\end{pmatrix}
\\
x = 1, y = - 1\\\\

λ = 5のとき\\
\begin{pmatrix}
-4 & 2\\
4 & -2
\end{pmatrix} 
\begin{pmatrix}
x\\
y
\end{pmatrix} =
\begin{pmatrix}
0\\
0
\end{pmatrix}
\\
x = 1, y = 2

\end{align}

それぞれ掛けた時に対応するx,yとλがかけられるようにすると以下になる
ここでλを対角成分とした行列をD,
ラムダに対応するxとyの組をPとする。


\begin{align}
\begin{pmatrix}
1 & 1\\
- 1 & 2
\end{pmatrix}
\end{align}
\begin{pmatrix}
-1 & 0\\
0 & 5
\end{pmatrix}
 = P D = AP

※いま求めたDPはAxと等しく、Pはxの組であるのでAP = DPとなる

今欲しいのは対角化されたDなので

PD = AP \\
D = P^{-1}AP

となり対角化できる。
DをK乗すると

D = P^{-1}AP \\
D^k = (P^{-1}AP) * (P^{-1}AP) * .... * (P^{-1}AP) = P^{-1}A^kP\\
PD^kP^{-1} = A^k

と変形することで元のAをK乗したのと同じ結果になる(例えば100乗とか計算するとき100回の行列演算がスカラーのK乗と数回程度の行列演算になるのでとても楽になる)

ここまでが前座

利用するシチュエーション(python使用)

こんな数列が与えられて、一般項を求めてくださいとか言われた時

S = 0, 2, 2, 6, 10, 22, 42, 86, 170, ...

まず項数が増えるときのルールを観察する。
どうやら以下のルールで増えてるっぽいことがわかった

項数=k\\
S_{k+1} = \left\{
\begin{array}{ll}
2S_k + 2 & k(mod2)=0 \\
2S_k - 2 & k(mod2)=1
\end{array}
\right.

かけたら次の項に遷移するSkがSk+1に遷移する行列Aを考える

\begin{align}
S_{k+1} = AS_k
\\
\begin{pmatrix}
S_{k+1} \\
2(-1)^{k+1}
\end{pmatrix}=
\begin{pmatrix}
2 & -1\\
0 & -1
\end{pmatrix}
\begin{pmatrix}
A_k \\
2(-1)^{k}
\end{pmatrix}
\end{align}

ここでAを対角化する。
python君では嬉しいことに行列を渡すと「行列の対角化って何？」の項のPとDを取得できる。
ここではPとA、P^(-1)からDを求めている

import sympy as sp

sp.init_printing()
A = sp.Matrix( [[2,-1],[0,-1]])

# 「行列の対角化って何？」の項のPと同じものを取得
P= A.eigenvects()[0][2][0].col_insert(1, A.eigenvects()[1][2][0])
P_inv = P.inv()

# 対角行列取得
D = P_inv * A * P

# 変数kを設定
k = sp.symbols("k")

# 一般項を出力
# 上の式がSkになっている
print(P * D**k * P.inv() * sp.Matrix([[0],[-2]])[0])

出力結果
一般項は下記式

for k in range(10):
    ans = (P * D**k * P.inv() * sp.Matrix([[0],[-2]]))[0]
    print(f"{k}項目:{ans}")

一般項の出力結果確認

0項目:0
1項目:2
2項目:2
3項目:6
4項目:10
5項目:22
6項目:42
7項目:86
8項目:170
9項目:342

このように行列の掛け算で表現できる数列の一般項は行列の対角化で簡単に作成できることがある。

例題 フィボナッチ数列の一般項(python使用)

3項間漸化式の場合も見てみよう。

S = 0, 1, 1, 2, 3, 5, 8, 13, 21, ...

このように表現できる

項数=k\\
S_{k+1} = S_{k} + S_{k-1}

かけたら次の項に遷移するSkがSk+1に遷移する行列Aを考える

\begin{align}
S_{k+1} = AS_k
\\
\begin{pmatrix}
S_{k+1} \\
S_{k}
\end{pmatrix}=
\begin{pmatrix}
1 & 1\\
1 & 0
\end{pmatrix}
\begin{pmatrix}
S_k \\
S_{k-1}
\end{pmatrix}
\end{align}

ここでSを対角化する

import sympy as sp

sp.init_printing()
A = sp.Matrix( [[1,1],[1,0]])

# 「行列の対角化って何？」の項のPと同じものを取得
P= A.eigenvects()[0][2][0].col_insert(1, A.eigenvects()[1][2][0])
P_inv = P.inv()

# 対角行列取得
# D = P_inv * A * P
D = sp.Matrix(np.diag(list([A.eigenvects()[0][0], A.eigenvects()[1][0]])))

# 変数kを設定
k = sp.symbols("k")
# 一般項を出力

# sp.Matrix([[1],[0]]))はSk1とSk0
# [1]でSkの一般項が取得できる
sp.simplify(P * D**k * P_inv *  sp.Matrix([[1],[0]]))[1]

出力結果
k番目の一般項が以下になる

for k in range(10):
    ans = sp.simplify(P * D**k * P_inv *  sp.Matrix([[1],[0]]))[1]
    print(f"{k}項目:{ans}")

一般項の出力結果確認

0項目:0
1項目:1
2項目:1
3項目:2
4項目:3
5項目:5
6項目:8
7項目:13
8項目:21
9項目:34