いつ役立つか
自己回帰モデルが暴走するかどうか調べたい!
自己回帰モデルとは
例えば、今日の$ξ(t)$は昨日の$ξ(t-1)$、おとといの$ξ(t-2)$、さらにその前の日$ξ(t-3)$と、今日の$u(t)$から、次のように決まるとします。
$ξ(t)=-0.5ξ(t-1)+0.34ξ(t-2)+0.08ξ(t-3)+2u(t)$
これを行列で表現すると、
x(t)=
\begin{pmatrix}
-0.5&0.34&0.08\\
1&0&0\\
0&1&0
\end{pmatrix} \times
x(t-1)
となります。
以降、これを一般化した
x(t)=Ax(t-1)
を扱います。
1次元の場合
x(t)=7x(t-1)
これは$t$が大きくなる度に7倍されるので、
x(t)=7^tx(0)
となり、$t→∞$の時、$7^t→∞$なので暴走します。
一方で係数が$0.2$の場合は暴走しない。
これは簡単ですね。
対角行列の場合
x(t)=
\begin{pmatrix}
5&0&0\\
0&-3&0\\
0&0&0.8
\end{pmatrix} \times
x(t-1)
の場合、対角行列のべき乗は成分をそのままべき乗するだけでよいので、
x(t)=
\begin{pmatrix}
5&0&0\\
0&-3&0\\
0&0&0.8
\end{pmatrix}^t
x(0)
と書ける。
一般化
対角行列でない行列の場合は、無理やり対角行列に変換してしまいます。
元の変数$x(t)$に対し、何か正則行列$P$を持ってきて
x(t)=Py(t)
で別の変数に変換する方法を考えてみましょう。
この時、$x(t)=Ax(t-1)$はどう変換されるのでしょうか？
まず、$x(t)=Py(t)$を
y(t)=P^{-1}x(t)=P^{-1}Ax(t-1)\\
=P^{-1}A(Py(t-1))=(P^{-1}AP)y(t-1)
$Λ=P^{-1}AP$と置くと、
y(t)=Λy(t-1)
と置ける。
この$Λ$がもし対角行列なら、
y(t)=Λ^ty(0)
で簡単に$y(t)$が求まり、同時に
x(t)=Py(t)=PΛty(0)=PΛ^tP^{-1}x(0)
でxも求まってめでたしめでたし。
この手順「都合の良い$P$を持ってきて$P^{-1}AP$を対角行列にする」ことを対角化と呼びます。
やりたいことは
P^{-1}AP≡Λ=diag(λ_1,...,λ_n)
のような$P$を見つけることです。
一般に、正方行列$A$に対して
Ap=λp\\
p≠o
を満たす数$λ$とベクトル$p$をそれぞれ固有値、固有ベクトルと呼びます。
Pythonで固有値、固有ベクトルを計算
Numpyの機能を駆使すれば簡単に計算が可能です。
次の行列で計算してみましょう。
A=
\begin{pmatrix}
5&3&4\\
6&8&-8\\
6&9&-9
\end{pmatrix}
import numpy as np
def get_eigenpairs(arr):
#np.linalg.eigは第一引数に固有値、第二引数に固有ベクトルを出力する
w, v = np.linalg.eig(arr)
eigenpairs = []
#vは1に規格化されているため。このままでは良く分からない。
#そのため、各列の0を除く最小値で割り、数値を元に戻している。
for i, val in enumerate(w):
vec = v[:, i] / np.min(np.abs(v[:, i][v[:, i] != 0]))
eigenpairs.append((val, vec))
eigenpairs.sort(key=lambda x:x[0])
return eigenpairs
A = np.array([[5,3,-4],[6,8,-8],[6,9,-9]])
get_eigenpairs(A)
# array([[-1.0, array([1., 2., 3.])],
# [2.0, array([-1., -3., -3.])],
# [3.0, array([1., 2., 2.])]], dtype=object)
固有ベクトルが3つ得られました。これらは固有値$-1, 2,3$にそれぞれ対応しています。
そこで固有ベクトルを3つ並べて
P=
\begin{pmatrix}
1&-1&1\\
2&-3&2\\
3&-3&2
\end{pmatrix}
とし、$P^{-1}AP$が対角行列になるか確認してみましょう。
P = np.empty((3,3))
# 3つの固有ベクトルを3×3の行列に格納
for i,val in enumerate(get_eigenpairs(a)):
P[i] = val[1]
np.rint(np.dot(np.dot(np.linalg.inv(P.T),a),P.T))
# array([[-1., 0., 0.],
# [ 0., 2., -0.],
# [ 0., 0., 3.]])
無事、$P^{-1}AP=diag(-1,2,3)$になることが確認できました。
これで簡単に自己回帰モデルの暴走について調べることができますね。