行列の対角化

行列の対角化

行列の対角化

 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}$ が互いに1次独立であるときには、$P$ に逆行列 $P^{-1}$ が存在します。(5) の両辺に左側から逆行列 $P^{-1}$ をかけると
 
\[P^{-1}AP=D\tag{6}\]
を得ます。すなわち、行列 $A$ を $P$ と $P^{-1}$ で挟むと固有値を対角成分にもつ行列 $D$ が得られます。この操作を 行列の対角化、行列 $P$ を 対角化行列 といいます。一般の $n$ 次正方行列 $A$ についても、$n$ 個の互いに1次独立な固有ベクトルが存在すれば、$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}\]
が成り立っています。

diag_eigval()

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

# SLA_034-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 ではないことを確認しておきます。

# SLA_034-2

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

det_a = np.linalg.det(a)

print(det_a)
10.000000000000002

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

# SLA_034-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 以降で導入された比較的新しい演算子です)。

# SLA_034-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 次正方行列であれば手計算で対角化できますが (それでも固有値を求めるのは面倒な作業です)、大きなサイズの行列では現実的に手計算では不可能なので、コンピュータの力を借りることになります。
 

退屈なことはPythonにやらせよう ―ノンプログラマーにもできる自動化処理プログラミング

新品価格
¥3,996から
(2019/8/6 11:42時点)

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

 $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 で大きな行列を対角化してトレースを計算してみます。

# SLA_034-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*}\]
 コードSLA_034-5 の対角化された行列 d の行列式を計算して、a の行列式と比較してみます (a の行列式はすでに計算されています)。

# SLA_034-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