Pythonで線形代数！〜行列・応用編（行列式・固有値）：数学×Pythonプログラミング入門（4/5 ページ）

[羽山博，著]

目標4： 行列を対角化して、行列のべき乗を簡単に求める

　4番目の目標は、行列の対角化です。行列の対角化を行うと行列のべき乗が簡単に求められるようになります。これはとても便利です。なお、これまで見てきたのは2次元の正方行列ばかりでしたが、n次元の正方行列でも行列式や固有値／固有ベクトルの考え方は同じです。そこで、以下のような3次元の行列を例としてやってみましょう。

　行列の対角化を行うには、やるべきことが幾つかあるので、3つのステップに分けて進めす。

• ステップ1： 固有値を対角行列Dにし、D²を求める
• ステップ2： 固有ベクトルを列ベクトルとして並べた行列Pの逆行列P⁻¹を求める
• ステップ3： A ＝ PDP⁻¹となることを確認する（この右辺への変形が対角化です）

　ステップ3までできたら、さらに、A² ＝ PD²P⁻¹となることも確認してみてください。一般に、Aⁿ ＝ PDⁿP⁻¹となり、行列Aのべき乗が簡単に求められるようになります。

　ステップ1は行列Aのn個の固有値λ₁,λ₂,...,λ_nを以下のような行列Dにするということです。大きな0は右上と左下が全て0であるということを表しています。

　これは、固有値を求め、NumPyのdiag関数の引数に指定すれば簡単にできます。このような行列Dを作った上で、D²を求めてみましょうということです。行列のべき乗はnumpy.linalgモジュールのmatrix_power関数に配列と指数を指定するだけで求められます。

　ステップ2は簡単です。逆行列はnp.linalgモジュールのinv関数に配列を指定するだけで求められます。

　ステップ3も簡単ですね。単にPDP⁻¹の値を求めてAと一致するか確認すればいいだけです。A²＝PD²P⁻¹となることも簡単に確認できますね。

4. 行列を対角化してべき乗を求めるためのコードを書く

　ここからはもうステップごとに淡々とコードを書くのみです。が、興味深い点も幾つかあるので、確認しつつ進めましょう。それらのポイントについては動画でも解説しています。ぜひご視聴ください。

動画3　対角化により行列のべき乗を求める

　まず、固有値を対角要素として持つ行列を作ります。

import numpy as np

A = np.array([[2, 3, 5],
              [3, -2, 4],
              [-2, 5, -2]])
lambda_value = np.linalg.eig(A)[0]  # 固有値だけを返す
D = np.diag(lambda_value)  # 固有値の対角行列を作る

print("固有値の対角行列：")
print(D)
print("固有値の対角行列を2乗したもの：")
print(np.linalg.matrix_power(D, 2))  # Dの2乗を求める
# 出力例：
# 固有値の対角行列：
# [[-5.57732961  0.          0.        ]
#  [ 0.         -0.71081023  0.        ]
#  [ 0.          0.          4.28813984]]
# 固有値の対角行列を2乗したもの：
# [[31.10660554  0.          0.        ]
#  [ 0.          0.50525119  0.        ]
#  [ 0.          0.         18.38814328]]

リスト8　固有値を対角行列にする
np.linalgモジュールのeig関数が返す値の0番目が固有値。NumPyのdiag関数を使って固有値を対角要素に持つ行列を作る。それをそのまま2乗すると、対角要素を2乗した行列が返される。

　リスト8で注目すべき点は、対角行列を2乗すると、各要素を2乗した行列になるということです。小数点以下のある値なので2乗した結果が分かりにくいですが、例えば、-5.57732961をそのまま2乗すると31.10660554になります。しかし、一般的な行列の場合、それぞれの要素をそのまま2乗しただけでは行列の2乗は求められません。以下の例を見ると一目瞭然です。

のとき、

　対角行列の場合は単に各要素を2乗するだけで済みます。これは2乗の場合だけでなく、何乗でも同じです。2乗ぐらいであれば大した手間ではありませんが、10乗、20乗となると（さすがに手作業でやろうという人はいないと思いますが）、必要な計算量がずいぶんと少なくなります。

　次は、固有ベクトルとその逆行列ですね。これは簡単です（リスト9を書く）。

import numpy as np

A = np.array([[2, 3, 5],
              [3, -2, 4],
              [-2, 5, -2]])
P = np.linalg.eig(A)[1]
Pinv = np.linalg.inv(P)

print("固有ベクトルを行列にしたもの(P)：")
print(P)
print("Pの逆行列：")
print(Pinv)  # Dの2乗を求める
# 出力例：
# 固有ベクトルを行列にしたもの(P)
# [[-0.24094928 -0.81389274  0.86890565]
#  [-0.62663324 -0.18340916  0.48310215]
#  [ 0.74113037  0.55130725  0.10777424]]
# Pの逆行列
# [[ 0.62263603 -1.23339331  0.50886926]
#  [-0.92616272  1.4579606   0.93161521]
#  [ 0.45600647  1.02363048  1.01374245]]

リスト9　固有ベクトルを行列にしたものと、その逆行列を求める
np.linalgモジュールのeig関数が返す値の1番目が固有ベクトルを配列にしたもの。逆行列が求められない場合はエラーとなるが、この例では逆行列が求められる。返される個々の値については気にしなくてもよい。

　さて、いよいよ対角化です。PDP⁻¹を求めてみましょう（リスト8とリスト9に続けてリスト10を書く）。

print(P @ D @ Pinv)  # リスト8とリスト9の続きに書く
# 出力例：
# [[ 2.  3.  5.]
#  [ 3. -2.  4.]
#  [-2.  5. -2.]]

リスト10　A＝PDP⁻¹を確認する
P、D、Pinvはリスト6で求めたもの。結果はちゃんとAになっている！

　対角化ができると、行列Aのべき乗が簡単に求められます。リスト11の通りです。

print("A^2：")
print(np.linalg.matrix_power(A, 2))
print("PD^2Pinv：")
print(P @ np.linalg.matrix_power(D, 2) @ Pinv)
# 出力例：
# Aの2乗
# [[  3  25  12]
#  [ -8  33  -1]
#  [ 15 -26  14]]
# PD^2Pinv
# [[  3.  25.  12.]
#  [ -8.  33.  -1.]
#  [ 15. -26.  14.]]

リスト11　A²とPD²P⁻¹が等しいことを確認する
A²とPD²P⁻¹は確かに等しくなっている。やはり2乗の場合だけでなく、何乗の場合でも同じ。単純に行列Aのべき乗を求めるには時間がかかるが、対角化すればDのべき乗を求め、左から固有ベクトルを行列にしたもの、右からその逆行列を掛ければ答えが求められる。Dのべき乗はリスト8で見たように簡単に求められる。

　行列式や固有値／固有ベクトルは、計算の手順が複雑なので、手計算でやろうとするとそれだけで嫌気が差してしまいます。「線形代数の難所」と言われるのもそのためかもしれません。しかし、NumPyなどを使って計算すれば簡単に答えが求められます。その分、本質的な意味を理解したり、応用したりすることにリソースが割けるのではないかと思います。

　以下の練習問題では、行列式と固有値／固有ベクトルの求め方をしっかり身に付けられるように、これまでに見てきた関数を使ってできる応用例を取り上げてみました。