行列の累乗
同じ行列を $k$ 回乗算することを 行列の累乗 とよび、$A^k$ で表します:
\[A^k=\begin{bmatrix}a&b\\c&d\end{bmatrix}\begin{bmatrix}a&b\\c&d\end{bmatrix}\cdots\begin{bmatrix}a&b\\c&d\end{bmatrix}\]
行列 $A$ の固有値 $\lambda$ に対応する固有ベクトルが $\boldsymbol{x}$ であるとします:
\[A\boldsymbol{x}=\lambda\boldsymbol{x}\tag{1}\]
両辺の左側からさらに $A$ を乗じてみます。
\[\begin{align*}A(A\boldsymbol{x})&=\lambda(A\boldsymbol{x})\\[6pt]A^2\boldsymbol{x}&=\lambda(\lambda\boldsymbol{x})\\[6pt]A^2\boldsymbol{x}&=\lambda^2\boldsymbol{x}\end{align*}\]
$A^2$ の固有値は $\lambda^2$ であり、その固有ベクトルは $A$ の固有ベクトルと同じです。
同様にして、$A^k$ の固有値は $\lambda^k$ であり、固有ベクトルは $\boldsymbol{x}$ となります:
\[A^k\boldsymbol{x}=\lambda^k\boldsymbol{x}\tag{2}\]
すなわち、$A$ の固有値と固有ベクトルがわかっていれば、$A^k$ の具体的な表式を知らなくても、その固有値と固有ベクトルは簡単に求めることができます。
対角化による行列の累乗計算
行列 $A$ を対角化できれば、行列の累乗 $A^k$ を簡単に求めることができます。対角化の式を再掲します。
\[P^{-1}AP=D\tag{3}\]
この式を変形すると、行列 $A$ は $P$ と $D$ を使って
\[A=PDP^{-1}\tag{4}\]
と表すことができます。$A$ の $2$ 乗は
\[A^2=(PDP^{-1})(PDP^{-1})\tag{5}\]
と書けますが、$P^{-1}P$ は単位行列となるので、
\[A^2=PD^2P^{-1}\tag{6}\]
となります。同様にして $A$ の $3$ 乗は
\[A^3=(PDP^{-1})(PDP^{-1})(PDP^{-1})=PD^3P^{-1}\tag{7}\]
となり、一般に $A$ の $k$ 乗は
\[A^k=PD^kP^{-1}\tag{8}\]
で表すことができます。$D$ は対角行列なので、
\[D^k=\begin{bmatrix}\alpha&0\\0&\beta\end{bmatrix}^k=\begin{bmatrix}\alpha^k&0\\0&\beta^k\end{bmatrix}\tag{9}\]
のように簡単に計算することができます。
前回と前々回の記事で扱った行列
\[A=\begin{bmatrix}3&1\\2&4\end{bmatrix}\tag{10}\]
の累乗を計算してみましょう。この行列の固有値は $\alpha=2,\ \beta=5$ であり、それぞれの固有値に対応する固有ベクトルは
\[\begin{bmatrix}1\\-1\end{bmatrix},\quad \begin{bmatrix}1\\2\end{bmatrix}\]
なので、
\[P=\begin{bmatrix}1&1\\-1&2\end{bmatrix},\quad D=\begin{bmatrix}2&0\\0&5\end{bmatrix}\tag{11}\]
となります。$P$ の逆行列を計算すると
\[P^{-1}=\frac{1}{3}\begin{bmatrix}2&-1\\1&1\end{bmatrix}\tag{12}\]
対角化された行列 $D$ を使って $A^3$ を計算してみます。
\[\begin{align*}A^3&=PD^3P^{-1}\\[6pt]&=\frac{1}{3}\begin{bmatrix}1&1\\-1&2\end{bmatrix}\begin{bmatrix}8&0\\0&125\end{bmatrix}\begin{bmatrix}2&-1\\1&1\end{bmatrix}\\[6pt]&=\frac{1}{3}\begin{bmatrix}141&117\\234&258\end{bmatrix}=\begin{bmatrix}47&39\\78&86\end{bmatrix}\end{align*}\tag{13}\]
numpy.linalg.matrix_power()
numpy.linalg.matrix_power(a, n) を使えば、累乗を高速計算できます。
例として、3 次正方行列
\[A=\begin{bmatrix}1&2&3\\4&5&6\\7&8&9\end{bmatrix}\]
の 5 乗を計算してみます。
# numpy_matrix_power
# In[1]
import numpy as np
# 3次正方行列
# [[1 2 3]
# [4 5 6]
# [7 8 9]]
a = np.arange(1, 10).reshape(3, 3)
a_5 = np.linalg.matrix_power(a, 5)
print(a_5)
# [[121824 149688 177552]
# [275886 338985 402084]
# [429948 528282 626616]]
実務上は numpy.linalg.matrix_power() を使えばよいのですが、次は敢えて対角化を使って累乗計算してみます。前回記事で作成した対角化関数 eig_diag() を再掲します。
# In[2]
# 対角化関数
def eig_diag(x):
eig = np.linalg.eig(x)
e = np.diag(eig[0])
p = eig[1]
return e, p
a を eig_diag() に渡して対角化された行列と、固有ベクトルを並べた行列を取得して a の 5 乗を計算します。d は対角行列なので、** 演算子を使って行列の累乗を計算できます (** は行列の要素ごとの累乗を計算する演算子ですが、対角行列については行列の累乗と同じ結果を返します)。
# In[3]
# aを対角化
e = eig_diag(a)
# 対角化された行列
d = e[0]
# 固有ベクトルを並べた行列
p = e[1]
# pの逆行列
ip = np.linalg.inv(p)
a_5 = p @ (d**5) @ ip
print(a_5)
# [[121824. 149688. 177552.]
# [275886. 338985. 402084.]
# [429948. 528282. 626616.]]
実行結果は In[1] と一致しています。
SymPyによる行列の累乗計算
行列の累乗の一般表式が必要なときは、計算機代数システム SymPy をインポートします (このパッケージは文字式を含む計算を実行できます)。SymPy の行列オブジェクトは ** 演算子で累乗を計算できます。
# sympy_matrix_power
# In[1]
# SymPyをインポート
import sympy
# 数式をLaTeXで表示
sympy.init_printing()
# 記号の定義
sympy.var("a b c d")
# 行列オブジェクトを生成
m = sympy.Matrix([[a, b],
[c, d]])
# m**3を計算して結果を展開
sympy.expand(m**3)
** 演算子
** はオブジェクトの種類によって異なる結果を返す演算子です。** を NumPy配列 (ndarray) に作用させると各要素を累乗して返します。
# python_elements_power
# In[1]
import numpy as np
a = np.array([[1, 2],
[3, 4]])
print(a**2)
# [[ 1 4]
# [ 9 16]]
これは行列の累乗とは異なる演算ですが、対角行列については
\[\begin{bmatrix}\alpha&0&0\\0&\beta&0\\0&0&\gamma\end{bmatrix}^k=\begin{bmatrix}
\alpha^k&0&0\\0&\beta^k&0\\0&0&\gamma^k\end{bmatrix}\]
が成り立つので、** を使って累乗を計算できます。
行列オブジェクト (numpy.matrix) に対して ** 演算子を使用すると、行列の累乗を返します。
# python_matrix_power
# In[1]
m = np.matrix([[1, 2],
[3, 4]])
print(m**2)
# [[ 7 10]
# [15 22]]
SymPy の行列オブジェクト (sympy.Matrix) に対しても、** 演算子は行列の累乗を返します。
# sympy_matrix_power
# In[1]
import sympy
# 数式をLaTeX形式で表示
sympy.init_printing()
sympy.var("a b c d")
# 行列オブジェクトを生成
m = sympy.Matrix([[a, b],
[c, d]])
# mの2乗を計算
m**2
a^{2} + b c & a b + b d\\a c + c d & b c + d^{2}
\end{matrix}\right]\]
下記は誤植と思われるので、ご確認ください。
(6)式の次の、A^3=… の式番号で、(4) → (7)
以降の式番号が1ずつ増えることになりますが…
あちこちに誤植あって、申し訳ないです。(_ _*)
直しておきました。本当に感謝です。