行列の対角化

当サイトではアフィリエイトプログラムを利用して商品を紹介しています。

【Pythonで学ぶ線形代数学講座(24)】行列の対角化

行列の対角化

2 次の正方行列 $A=\begin{bmatrix}a&b\\c&d\end{bmatrix}$ について、
 \[\boldsymbol{p}=\begin{bmatrix}p_1\\p_2\end{bmatrix},\quad \boldsymbol{q}=\begin{bmatrix}q_1\\q_2\end{bmatrix}\tag{1}\]
を固有値 $\alpha,\ \beta$ に対する $A$ の固有ベクトルとします。すなわち、
 \[A\begin{bmatrix}p_1\\p_2\end{bmatrix}=\alpha\begin{bmatrix}p_1\\p_2\end{bmatrix},\quad A\begin{bmatrix}q_1\\q_2\end{bmatrix}=\beta\begin{bmatrix}q_1\\q_2\end{bmatrix}\tag{2}\]
が成り立っているとします。この式は
 \[A\begin{bmatrix}p_1&q_1\\p_2&q_2\end{bmatrix}=\begin{bmatrix}p_1&q_1\\p_2&q_2\end{bmatrix}\begin{bmatrix}\alpha&0\\0&\beta\end{bmatrix}\tag{3}\]
という形にまとめることができます。
 \[P=\begin{bmatrix}p_1&q_1\\p_2&q_2\end{bmatrix},\quad D=\begin{bmatrix}\alpha&0\\0&\beta\end{bmatrix}\tag{4}\]
とおくと、
 \[AP=PD\tag{5}\]
と表せます。ベクトル $\boldsymbol{p}$ と $\boldsymbol{q}$ が互いに線形独立であるときには、$P$ に逆行列 $P^{-1}$ が存在します。(5) の両辺に左側から逆行列 $P^{-1}$ をかけると
 \[P^{-1}AP=D\tag{6}\]
を得ます。すなわち、行列 $A$ を $P$ と $P^{-1}$ で挟むと固有値を対角成分にもつ行列 $D$ が得られます。この操作を行列の対角化、行列 $P$ を対角化行列といいます。一般の $n$ 次正方行列 $A$ についても、$n$ 個の互いに線形独立な固有ベクトルが存在すれば、$A$ は必ず対角化できます。

一例として前回記事で扱った行列
 \[A=\begin{bmatrix}3&1\\2&4\end{bmatrix}\]
を考えます。この行列の固有値は $\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{13}\]
となります。$P$ の逆行列を計算すると
 \[P^{-1}=\frac{1}{3}\begin{bmatrix}2&-1\\1&1\end{bmatrix}\tag{14}\]
となります。$P^{-1}AP$ を計算してみると、
 \[\begin{align*}P^{-1}AP&=\frac{1}{3}\begin{bmatrix}2&-1\\1&1\end{bmatrix}\begin{bmatrix}3&1\\2&4\end{bmatrix}\begin{bmatrix}1&1\\-1&2\end{bmatrix}\\[6pt]&=\frac{1}{3}\begin{bmatrix}6&0\\0&15\end{bmatrix}=\begin{bmatrix}2&0\\0&5\end{bmatrix}\end{align*}\]
となって確かに
 \[P^{-1}AP=D\tag{6}\]
が成り立っています。

eig_diag()

NumPy を使って行列を対角化してみます。
前回記事で解説した numpy.linalg.eig() は固有値と固有ベクトルを (1 次元配列, 2 次元配列) のタプルで返す関数でした。これを少し改造して、固有値を対角行列の形で返すようにします。

# numpy_diagonalization

# In[1]

import numpy as np

# 対角化関数
def eig_diag(x):
    eig = np.linalg.eig(x)
    e = np.diag(eig[0])
    p = eig[1]
    return e, p

eig_diag() に行列 a を渡すと、対角化された行列 d と固有ベクトルを格納した 2 次元配列が返ります。先の例と同じ 2 次正方行列を定義して、行列式が 0 ではないことを確認しておきます。

# In[2]

# 2次正方行列を定義
a = np.array([[3, 1],
              [2, 4]])

det_a = np.linalg.det(a)

print(det_a)
# 10.000000000000002

a を eig_diag() に渡してみます。

# In[3]

# aの固有値と固有ベクトルを取得
e = eig_diag(a)

# 対角化された行列と固有ベクトルを別々の変数に格納
d = e[0]
p = e[1]

print(d)
print(p)

# [[2. 0.]
#  [0. 5.]]
# [[-0.70710678 -0.4472136 ]
#  [ 0.70710678 -0.89442719]]

p の逆行列 ip を求めて、ip, a, p の行列積が d に等しいことを確認してみます。このように複数の行列の積を計算するときは @ 演算子を使うと便利です (Python 3.5 以降で導入された比較的新しい演算子です)。

# In[4]

# pの逆行列
ip = np.linalg.inv(p)

# ip,a,pの行列積を計算
d2 = ip @ a @ p

print(d2)
# [[2.00000000e+00 5.55111512e-16]
#  [0.00000000e+00 5.00000000e+00]]

若干の誤差はありますが、d と等しい値が出力されています。
2 次正方行列であれば手計算で対角化できますが (それでも固有値を求めるのは面倒な作業です)、大きなサイズの行列では現実的に手計算では不可能なので、コンピュータの力を借りることになります。

対角化された行列のトレース

$A$ と $D$ のトレースは一致します。
 
\[\mathrm{tr}A=\mathrm{tr}D\tag{7}\]
この定理は $\mathrm{tr}XY=\mathrm{tr}YX$ を使って簡単に証明できます:
 
\[\mathrm{tr}D=\mathrm{tr}PAP^{-1}=\mathrm{tr}P^{-1}PA=\mathrm{tr}A\]
NumPy で大きな行列を対角化してトレースを計算してみます。

# In[5]

# 10×10正方行列を作成
a = np.arange(1, 101).reshape(10, 10)

# aの行列式を計算
det_a = np.linalg.det(a)

# 対角化された行列と固有ベクトルを取得
e = eig_diag(a)

# 対角化された行列
d = e[0]

# aのトレース
tr_a = np.trace(a)

# dのトレース
tr_d = np.trace(d)

print("det_a : {}".format(det_a))
print("tr_a : {}".format(tr_a))
print("tr_d : {}".format(tr_d))

# det_a : 3.496590348060813e-113
# tr_a : 505
# tr_d : (505.00000000000017+0j)

対角化された行列の行列式

$A$ と $D$ の行列式は一致します:
 \[\det A=\det D\tag{8}\]
この定理は $\det(XY)=\det X\det Y$ を使って証明できます:
 \[\begin{align*}\det D&=\det(PP^{-1})\det D\\[6pt]&=\det P\det P^{-1}\det D\\[6pt]&=\det P\det D\det P^{-1}\\[6pt]&=\det(PD)\det P^{-1}\\[6pt]&=\det(PDP^{-1})=\det A\end{align*}\]
In[3] の対角化された行列 d の行列式を計算して、a の行列式と比較してみます (a の行列式はすでに計算されています)。

# In[6]

# dの行列式
det_d = np.linalg.det(a)

print("det_a : {}".format(det_a))
print("det_d : {}".format(det_d))

# det_a : 3.496590348060813e-113
# det_d : 3.496590348060813e-113

 

コメント

  1. HNaito より:

    下記は誤植と思われますので、ご確認ください。
    diag_eigval( ) → eig_diag( )
    プログラムの間数名に合わせた方がいいと思いました。

  2. あとりえこばと より:

    ChatGPT が行列の対角化を先生と生徒の対話形式で解説してくれました。

    生徒: 先生、行列の対角化って何ですか?聞いたことはありますが、実際に何をするものなのか分かりません。

    先生: それは素晴らしい質問です!行列の対角化は、行列を対角行列に変換することです。対角行列とは、非対角要素がすべてゼロで、対角要素だけが値を持つ行列のことです。

    生徒: なるほど、でもなぜ行列を対角行列に変換したいのでしょうか?

    先生: 行列を対角行列に変換すると、いくつかの計算が非常に簡単になります。たとえば、行列のべき乗や指数関数の計算が簡単になるんです。

    生徒: ほんとうですか?具体的な例を教えてもらえますか?

    先生: もちろんです!例えば、行列Aが対角行列Dで対角化されるとします。すると、Aのべき乗は簡単に計算できます。Aのn乗は、Dの対角要素をそれぞれn乗した対角行列になるんです。

    生徒: なるほど、それは便利ですね!でも、どうやって行列を対角化するのでしょうか?

    先生: 行列を対角化するためには、固有値と固有ベクトルを使います。まず、行列Aの固有値を求め、それぞれの固有値に対応する固有ベクトルを求めます。それから、それらの固有ベクトルを並べた行列をPとすると、Pの逆行列P^-1と対角行列Dを使って次のように表すことができます:A = PDP^-1。

    生徒: なるほど、固有値と固有ベクトルを使ってPとDを求めて、それを使ってAを表すんですね。でも、それはどうやって計算するのでしょう?

    先生: 固有値と固有ベクトルを求める方法はいくつかありますが、一般的には行列の特性方程式を解くことで求めることができます。特性方程式は、AからλIを引いた行列の行列式がゼロになるようなλを求める方程式です。