行列の累乗

行列の累乗

行列の累乗

 同じ行列を $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^k$ の計算方法を解説します。

行列の累乗の固有値

 行列 $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{4}\]
となり、一般に $A$ の $k$ 乗は
 
\[A^k=PD^kP^{-1}\tag{7}\]
で表すことができます。$D$ は対角行列なので、
 
\[D^k=\begin{bmatrix}\alpha&0\\0&\beta\end{bmatrix}^k
=\begin{bmatrix}\alpha^k&0\\0&\beta^k\end{bmatrix}
\tag{8}\]
のように簡単に計算することができます。

 前回と前々回の記事で扱った行列
 
\[A=\begin{bmatrix}3&1\\2&4\end{bmatrix}\tag{9}\]
の累乗を計算してみましょう。この行列の固有値は $\alpha=2,\ \beta=5$ であり、それぞれの固有値に対応する固有ベクトルは $\begin{bmatrix}1\\-1\end{bmatrix},\ \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{10}\]
となります。$P$ の逆行列を計算すると
 
\[P^{-1}=\frac{1}{3}\begin{bmatrix}2&-1\\1&1\end{bmatrix}\tag{11}\]
 対角化された行列 $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*}\]

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 乗を計算してみます。

# リストSLA035-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() を再掲します。

# リストSLA035-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 は対角行列なので、** 演算子を使って行列の累乗を計算できます (** は行列の要素ごとの累乗を計算する演算子ですが、対角行列については行列の累乗と同じ結果を返します)。

# リストSLA035-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.]]

 実行結果はリストSLA035-1 と一致しています。
 

SymPyによる行列の累乗計算

 行列の累乗の一般表式が必要なときは、計算機代数システム SymPy をインポートします (このパッケージは文字式を含む計算を実行できます)。SymPy の行列オブジェクトは ** 演算子で累乗を計算できます。

# リストSLA036

# Jupyter Notebook用コード

# 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)
\[\left[\begin{matrix}a^{3} + 2 a b c + b c d & a^{2} b + a b d + b^{2} c + b d^{2}\\a^{2} c + a c d + b c^{2} + c d^{2} & a b c + 2 b c d + d^{3}\end{matrix}\right]\]

** 演算子

 ** はオブジェクトの種類によって異なる結果を返す演算子です。** を NumPy配列 (ndarray) に作用させると各要素を累乗して返します。

# リストSLA037-1

import numpy as np

a = np.array([[1, 2],
              [3, 4]])

print(a**2)
[[ 1  4]
 [ 9 16]]

 これは行列の累乗とは異なる演算ですが、対角行列については
 
\[\begin{bmatrix}
\alpha&0&\\0&\beta&0\\0&0&\gamma
\end{bmatrix}^k
=\begin{bmatrix}
\alpha^k&0&\\0&\beta^k&0\\0&0&\gamma^k
\end{bmatrix}\]
が成り立つので、** を使って累乗を計算できます。

 行列オブジェクト (numpy.matrix) に対して ** 演算子を使用すると、行列の累乗を返します。

# リストSLA037-2

m = np.matrix([[1, 2],
               [3, 4]])

print(m**2)
[[ 7 10]
 [15 22]]

 SymPy の行列オブジェクト (sympy.Matrix) に対しても、** 演算子は行列の累乗を返します。

# リストSLA038

# Jupyter Notebook用コード

# SymPyをインポート
import sympy

# 数式をLaTeX形式で表示
sympy.init_printing()
sympy.var("a b c d")

# 行列オブジェクトを生成
m = sympy.Matrix([[a, b],
                  [c, d]])

# mの2乗を計算
m**2
\[\left[\begin{matrix}
a^{2} + b c & a b + b d\\a c + c d & b c + d^{2}
\end{matrix}\right]\]