ケーリー・ハミルトンの定理

≪【前の記事】ユニタリ行列

ケーリー・ハミルトンの定理

「固有値と固有ベクトル」で学んだように、行列 $A$ の固有値 $\lambda$ は固有方程式
 \[\det(A-\lambda I)=0\tag{1}\]
を解いて得られました。左辺は固有多項式とよばれる式で、$p(\lambda)$ のように表されます。行列 $A$ が二次正方行列であるとき、
 　\[A=\begin{bmatrix}a&c\\b&d\end{bmatrix}\tag{2}\]
とおいて、$p(\lambda)$ の具体的な表式を書き表してみると、
 \[\begin{align*}p(\lambda)&=\begin{vmatrix}a-\lambda & b\\c & d-\lambda\end{vmatrix}\\[6pt]&=(a-\lambda)(d-\lambda)-bc\\[6pt]&=\lambda^{2}-(a+d)\lambda^{1}+(ad-bc)\lambda^{0}\end{align*}\tag{3}\]
となります。$\lambda$ は固有値なので、もちろんスカラーですが、ここで敢えて $\lambda$ を行列 $A$ で置き換えて
 \[p(A)=A^2-(a+d)A^1+(ad-bc)I\tag{4}\]
という式を定義してみると、任意の二次正方行列 $A$ に対して $p(A)$ は必ず零行列になります。すなわち、
 \[A^2-(a+d)A+(ad-bc)I=O\tag{5}\]
が成り立ちます ($O$ は零行列)。これを ケーリー・ハミルトンの定理 とよびます（イギリスの数学者 アーサー・ケイリー とアイルランドの数学者 ウィリアム・ローワン・ハミルトン に由来します）。
 
実際に (5) の左辺を計算してみれば、零行列になることを簡単に確かめられるので証明は割愛します。$a+d$ は $A$ のトレース、$ad-bc$ は $A$ の行列式なので、(5) は
 \[A^2-(\mathrm{tr}A)A+(\det A)I=O\tag{6}\]
と書くこともできます。
 
NumPy で二次正方行列 $M$ の $p(M)$ を計算する式を定義して、ランダムに生成されたいくつかの行列について、$p(M)$ が零行列になることを確かめてみましょう。

# In[1]

import numpy as np

# ケーリー・ハミルトンの定理の左辺p(M)を計算する関数
def CayleyHamilton(M):
    I = np.eye(2)
    return np.linalg.matrix_power(M,2) - np.trace(M)*M + np.linalg.det(M)*I

np.random.seed(0)

for i in range(3):
    a, b, c, d = -10 + np.random.randint(20, size=4)
    M = np.array([[a, b], [c, d]])
    p_M = CayleyHamilton(M)
    print("M:\n{}\np(M):\n{}\n".format(M, p_M))

'''実行結果
M:
[[  2   5]
 [-10  -7]]
p(M):
[[0. 0.]
 [0. 0.]]

M:
[[-7 -3]
 [-1  9]]
p(M):
[[2.84217094e-14 0.00000000e+00]
 [0.00000000e+00 2.84217094e-14]]

M:
[[ 8 -6]
 [-4  2]]
p(M):
[[1.77635684e-15 0.00000000e+00]
 [0.00000000e+00 1.77635684e-15]]
'''

数値計算の性質上、若干の誤差は生じているものの、すべての行列 $M$ についてケーリー・ハミルトンの定理が成り立っていることがわかります。
 
ケーリー・ハミルトンの定理は行列 $A^{n}$ と $A^{n-1}$ (二次の場合は $A^2$ と $A$) の関係を示すので、これを漸化式のように用いて行列 $A$ の累乗 $A^n$ を計算できます。二次正方行列
 \[A=\begin{bmatrix}a&c\\b&d\end{bmatrix}\tag{8}\]
について、$A^n$ の一般表式を求めてみましょう。多項式
 \[p(\lambda)=\lambda^2-(\mathrm{tr}A)\lambda+\det A\tag{9}\]
について、
 
(ⅰ) $p(\lambda)=0$ が異なる 2 個の解 $\alpha,\ \beta$ をもつとき、すなわち $A$ の固有値が $\alpha$ と $\beta$ であるとき、解と係数の関係により、
 \[\lambda^2-(\alpha + \beta)\lambda + \alpha\beta=0\tag{11}\]
が成り立つので、この式が $A$ の固有方程式となります。ケーリー・ハミルトンの定理より、
 \[A^2-(\alpha + \beta)A + \alpha\beta I=0\tag{12}\]
が常に成り立ちます。この式を次のように 2 通りに書き表します。
 \[\begin{align*}A^2-\alpha A &= \beta A – \alpha \beta I\tag{13}\\[6pt]A^2-\beta A &= \alpha A – \alpha \beta I\tag{14}\end{align*}\]
少し整理すると、
 \[\begin{align*}A(A-\alpha I) &= \beta(A-\alpha I)\tag{15}\\[6pt]A(A-\beta I) &= \alpha(A-\beta I)\tag{16}\end{align*}\]
(15) の両辺に左から $A$ をかけてみると、
 \[A^2(A-\alpha I)=\beta A(A-\alpha I)\tag{16}\]
となりますが、右辺の $A(A-\alpha I)$ は (15) の左辺と同じなので、
 \[A^2(A-\alpha I)=\beta^2 (A-\alpha I)\tag{17}\]
となります。同様に繰り返し左から $A$ をかけると、
 \[A^{n-1}(A-\alpha I)=\beta^{n-1} (A-\alpha I)\tag{18}\]
という関係式が得られます。(16) についても同じようにして、
 \[A^{n-1}(A-\beta I)=\alpha^{n-1} (A-\beta I)\tag{19}\]
となります。(19) × $\alpha$ – (18) × $\beta$ を計算すると、
 \[(\alpha – \beta)A^n=\alpha^{n}(A-\beta I)-\beta^{n}(A-\alpha I)\tag{20}\]
となります。これを整理して、
 \[A^n=\frac{\alpha^n-\beta^n}{\alpha-\beta}A-\frac{\alpha^n\beta-\alpha\beta^n}{\alpha-\beta}\tag{21}\]
を得ます。

(ⅱ) $p(\lambda)=0$ が重解をもつとき、すなわち $\alpha=\beta$ のとき、(15) は
 \[A(A-\alpha I)=\alpha(A-\alpha I)\tag{22}\]
となるので、両辺に左から $A$ を繰り返し乗じると、
 \[A^n(A-\alpha I)=\alpha^n(A-\alpha I)\tag{23}\]
を得ます。両辺を $\alpha^{n+1}$ で割ると、
 \[\frac{A^{n+1}}{\alpha^{n+1}}-\frac{A^n}{\alpha^n}=\frac{A-\alpha I}{\alpha}\tag{24}\]
となります。ここで数列でよく使われるテクニックを用います。(24) の $A$ の次数を 1 つずつ下げた式を書き並べてみます。ただし、$n\geq 2$ とします。
 \[\begin{align*}\frac{A^{n}}{\alpha^{n}}-\frac{A^{n-1}}{\alpha^{n-1}}&=\frac{A-\alpha I}{\alpha}\\[6pt]\frac{A^{n-1}}{\alpha^{n-1}}-\frac{A^{n-2}}{\alpha^{n-2}}&=\frac{A-\alpha I}{\alpha}\\[6pt]\frac{A^{n-2}}{\alpha^{n-2}}-\frac{A^{n-3}}{\alpha^{n-3}}&=\frac{A-\alpha I}{\alpha}\\[6pt]&\vdots\\[6pt]\frac{A^{3}}{\alpha^{3}}-\frac{A^{2}}{\alpha^{2}}&=\frac{A-\alpha I}{\alpha}\\[6pt]\frac{A^{2}}{\alpha^{2}}-\frac{A^{1}}{\alpha^{1}}&=\frac{A-\alpha I}{\alpha}\end{align*}\tag{25}\]
これらを全部足し合わせると、隣り合う項が打ち消し合って、
 \[\frac{A^{n}}{\alpha^{n}}-\frac{A}{\alpha}=\frac{n(A-\alpha I)}{\alpha}\tag{26}\]
となります。この式を整理すると、
 \[A^n=n\alpha^{n-1}A-\alpha^{n}(n-1)I\tag{27}\]
を得ます。これは $n=1$ のときも成り立ちます。
 
(21) と (27) を Python で実装してみます。

# In[2]

# ケーリー・ハミルトンの定理を用いて行列の累乗を計算する関数
# M：二次正方行列　n：指数
# eps：固有値が縮退しているとみなす差分閾値
def matrix_power_2d(M, n, eps):

    # Mの固有値
    a, b = np.linalg.eigvals(M)

    # 単位行列
    I = np.eye(2)

    # 固有値の差の絶対値がeps未満であれば固有値は縮退(重複)しているとして処理
    if abs(a-b) < eps:
        Mn = n * (a**(n-1)) * M -(a**n) * (n - 1) * I

    # 異なる2つの固有値が存在する場合の処理
    else:
        Mn = (a**n-b**n) * A/(a-b) - (a**n*b-a*b**n) * I/(a-b)

    # 固有値とM^nを返す
    return (a, b), Mn

上のコードにもコメントしてありますが、数値計算では誤差を避けられないので、本来は固有値が縮退 (重複) している場合であっても、linalg.eigvals() の戻り値が僅かに異なることがあります。たとえば、行列
 \[A=\begin{bmatrix}-2 & -1\\4 & -6\end{bmatrix}\tag{28}\]
は重複する負整数の固有値 $\lambda=-4$ をもちますが、linalg.eigvals() に渡してみると次のような実行結果となります。

# In[3]

A = np.array([[-2, -1],
              [ 4, -6]])

# Aの固有値を求める
eig = np.linalg.eigvals(A)

print(eig)

# [-3.99999997 -4.00000003]

　得られた固有値の差分の絶対値は本当に僅かです。一応確かめてみましょう。

# In[4]

# 得られた固有値の差分の絶対値
e = abs(eig[0]-eig[1])

print(e)

# 6.664001839240541e-08

ですから、matrix_power_2d() 関数では、固有値の差分の絶対値が eps で指定した値以下ならば、縮退しているものとみなして、つまり式 (27) に基づいて処理するようにしてあります。それでは、行列
 \[A=\begin{bmatrix}3 & -9\\0 & 7\end{bmatrix}\tag{29}\]
を matrix_power_2d() に渡して、$A^5$ を計算してみます。

# In[3]

A = np.array([[3, -9],
              [0,  7]])

matrix_power_2d(A, 5, 1e-4)

# ((3.0, 7.0),
#  array([[   243., -37269.],
#         [     0.,  16807.]]))

　numpy.linalg.matrix_power() にも同じ行列を渡して検算します。

# In[4]

np.linalg.matrix_power(A, 5)

# array([[   243, -37269],
#        [     0,  16807]])

値がぴたり一致しました。数値計算では特に有用とはいえませんが、手計算で二次正方行列の累乗を計算するときは便利です（三次以上は固有方程式が複雑になるので事実上使えません）。