行列の累乗

【Pythonで学ぶ線形代数学講座(25)】行列の累乗計算

行列の累乗

同じ行列を $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