線形代数学2(ジョルダン標準形から行列の指数関数)

行列と連立方程式

連立方程式を行列表記すると以下のようになる。
$$\left[\begin{matrix}1 & 2 & 3\\2 & 3 & 4\\1 & 1 & 2\end{matrix}\right]
\left[\begin{matrix}x\\y\\z\end{matrix}\right] =
\left[\begin{matrix}-1\\-2\\3\end{matrix}\right]$$
そもそもこの連立方程式は解を持つのだろうか。

rank (階数)

$$A=\left[\begin{matrix}1 & 2 & 3\\2 & 3 & 4\\1 & 1 & 2\end{matrix}\right],
B=\left[\begin{matrix}-1\\-2\\3\end{matrix}\right],
A_B=\left[\begin{matrix}1 & 2 & 3 & -1\\2 & 3 & 4 & -2\\1 & 1 & 2 & 3\end{matrix}\right]とする。$$

$A$の行ベクトルのうち線形独立なものの個数を$A$のrankといい${\rm rank}(A)$と書く。

${\rm rank}(A)$は$A$の列ベクトルのうち線形独立なものの個数にも等しい。

A_s = sym.Matrix([
    [1,2,3],
    [2,3,4],
    [1,1,2]
])
Ab_s = sym.Matrix([
    [1,2,3,-1],
    [2,3,4,-2],
    [1,1,2, 3]
])
A_s.rank()
Ab_s.rank()
A_n = np.array([
    [1,2,3],
    [2,3,4],
    [1,1,2]
])
Ab_n = np.array([
    [1,2,3,-1],
    [2,3,4,-2],
    [1,1,2, 3]
])
np.linalg.matrix_rank(A_n) 
np.linalg.matrix_rank(Ab_n)

$3$
$3$
${\rm rank}(A_B)\ge{\rm rank}(A)$は常に成り立つが、解を持つためには
${\rm rank}(A_B)={\rm rank}(A)$である必要がある。
${\rm rank}(X)$の最大値は$X$の縦横の次数の小さい方。
このとき$X$はフルランクであるという。フルランクでないときはどうなるか。
${\rm rank}(A)=2$のときは解に1次元分の自由度があり、${\rm dim(Ker}(A))=1$と書く。
このときは$AX$で移される$B$はKerの分を引いた2次元分の自由度しかなく、${\rm dim(Im}(A))=2$となる。

行列式

正方行列$A$の行列式を$|A|$と書く。
$(1, 0,\dots,0), (0, 1,\dots, 0), (0, 0,\dots, 1)$で張られる超立方体の体積が
$|A|$をかけると何倍になるかの意味合い。
負の数にもなるので$|A|$の絶対値が重要。
$|A|=0$なら次元が潰れているということ。

A_s = sym.Matrix([
    [1,2,3],
    [2,3,4],
    [1,1,2]
])
A_s.det()
A_n = np.array([
    [1,2,3],
    [2,3,4],
    [1,1,2]
])
np.linalg.det(A_n)

$-1$

1次正方行列の行列式は要素そのもの。
n+1次正方行列の行列式はn+1個のn次正方行列の和として帰納的に定義される。 $${\begin{vmatrix}1 & 2 & 3\\2 & 3 & 4\\1 & 1 & 2\end{vmatrix}} =-2\begin{vmatrix}2 & 3\\1 & 2\end{vmatrix} +3\begin{vmatrix}1 & 3\\1 & 2\end{vmatrix} -4\begin{vmatrix}1 & 2\\1 & 1\end{vmatrix} $$

J = 1
sum([(-1)**(loop_i+J)*A_s[J, loop_i]*A_s.minorMatrix(J, loop_i).det() for loop_i in range(3)])

$-1$
真面目に余因子展開すると項の数が$n!$になる。
手計算だと$0$が多い行や列で展開すると楽になる。

逆行列

$|A|!=0$なら$|A|$を分母にした式で逆行列$A^{-1}$が定義される。
$A^{-1}A=AA^{-1}=I$

A_s.inverse_ADJ()
np.linalg.inv(A_n)

$$A^{-1}=\left[\begin{matrix}-2 & 1 & 1\\0 & 1 & -2\\1 & -1 & 1\end{matrix}\right]$$

A_n.dot(np.linalg.inv(A_n))

array([[ 1., 0., 0.],
[ 0., 1., 0.],
[ 0., 0., 1.]])

逆行列が存在する行列は正則であるという。
$|(A^{-1})|=\frac{1}{|A|}$

inverse_GEとかinverse_LUとか他のアルゴリズムのもある。
先程のsin関数でsympyとnumpyの差があったので
理論上逆行列が存在しない行列に対してnumpyで逆行列を計算してみよう。

X = np.sin(np.array([
    [np.pi/2, 3*np.pi/4],
    [0,     np.pi]
]))
Y = np.linalg.inv(X)
print(X)
print(Y)
print(X.dot(Y))

[[ 1.00000000e+00 7.07106781e-01]
[ 0.00000000e+00 1.22464680e-16]]
[[ 1.00000000e+00 -5.77396505e+15]
[ 0.00000000e+00 8.16561968e+15]]
[[ 1. 0.]
[ 0. 1.]]
Yの要素の絶対値が大きな数になっています。
sympyだと0除算エラーが出て理論値通りです。

固有値と固有ベクトル

$|A|$はベクトル$u$が$Au$でどう変化するかについての情報を与えているが、もう少し詳しく知りたい。
$A$をかけても方向が変わらないベクトルを固有ベクトルという。固有値を$\lambda$とおくと、
$$Au=\lambda u$$
とかける。
$$(A-\lambda I)u= 0$$
$|(A-\lambda I)| = 0$を特性多項式という。
特性多項式の解となる$\lambda$が固有値となる。
特性多項式は$\lambda$に関するn次方程式であるから、

n次正方行列の固有値は複素数の範囲でn個ある。
n次実正方行列の固有値に複素数が現れるときは共役複素数も固有値となる。

A_s = sym.Matrix([
    [1,2,3],
    [2,3,4],
    [1,1,2]
])
# A_s.eigenvals() 3次方程式がきれいにとけないとsympyは辛い
# A_s.eigenvects()

A_n = np.array([
    [1,2,3],
    [2,3,4],
    [1,1,2]
])
Vals, VectsT = np.linalg.eig(A_n)
Vects = np.transpose(VectsT)
print(Vals) # 固有値
print(A_n.dot(Vects[0])) # 固有ベクトルの1つ目
print(Vals[0] * Vects[0]) # 固有値の1つ目をかけたものと等しくなる

[ 5.97196077 -0.39542603 0.42346527]
[-3.08475786 -4.72087273 -1.9651832 ]
[-3.08475786 -4.72087273 -1.96518321]

B_s = sym.Matrix([
    [1, 1,-1],
    [0, 1, 0],
    [0, 0, 1]
])
Vals = B_s.eigenvals()
Vects = B_s.eigenvects()

$${1: 3}\\
\left [ \left ( 1, \quad 3, \quad \left [ \left[\begin{matrix}1\\0\\0\end{matrix}\right], \quad \left[\begin{matrix}0\\1\\1\end{matrix}\right]\right ]\right )\right ]$$
特性多項式が重解を持ち、かつ固有ベクトルが次元の数だけ求まらない場合もある。

対角化

正方行列$A$の固有ベクトル$u$を使って対角化ができる。
$P=(u_1, u_2,\dots, u_n)$とすると、
$$
P^{-1}AP=
\left[
\begin{array}{cccc}
\lambda_1 & 0 & \ldots & 0 \\
0 & \lambda_2 & \ldots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \ldots & \lambda_n
\end{array}\right]
$$

A_s = sym.Matrix([
    [-1, 1, 1],
    [ 1,-1, 1],
    [ 1, 1,-1]
])
P, M = A_s.diagonalize()
P * M * P.inverse_ADJ()

A_n = np.array([
    [-1, 1, 1],
    [ 1,-1, 1],
    [ 1, 1,-1]
])
L_n, P_n = np.linalg.eig(A_n) # 固有値と固有ベクトルを求める
M_n = np.diag(L_n) # 対角化
P_n.dot(M_n).dot(np.linalg.inv(P_n)) # 結果がもとに戻ることを確認

M
P

$$\left[\begin{matrix}-2 & 0 & 0\\0 & -2 & 0\\0 & 0 & 1\end{matrix}\right],
\left[\begin{matrix}-1 & -1 & 1\\1 & 0 & 1\\0 & 1 & 1\end{matrix}\right]$$

ジョルダン標準形

k次のジョルダン細胞を
$$J_{\lambda}(k)=\left[
\begin{array}{cccc}
\lambda & 1 & \ldots & 0 \\
0 & \lambda & \ddots & 0 \\
\vdots & \vdots & \ddots & 1\\
0 & 0 & \ldots & \lambda
\end{array}\right]$$
とする。
任意の正方行列は
$$
P^{-1}AP=
\left[
\begin{array}{cccc}
J_{\lambda 1} & 0 & \ldots & 0 \\
0 & J_{\lambda 2} & \ldots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \ldots & J_{\lambda n}
\end{array}\right]
$$
の形に書ける。
足りない固有ベクトルは
$$(A-\lambda I)^{l}u_{m+l}=u_m$$
で求める。
ここの部分がジョルダン細胞での2次以上の部分に相当する。

B_s = sym.Matrix([
    [1, 1,-1],
    [0, 1, 0],
    [0, 0, 1]
])
B_s.jordan_form()

$$\left ( \left[\begin{matrix}-2 & 0 & 0\\0 & -1 & 1\\0 & 1 & 1\end{matrix}\right], \quad \left[\begin{matrix}1 & 1 & 0\\0 & 1 & 0\\0 & 0 & 1\end{matrix}\right]\right )$$

種々の行列の性質

トレース

正方行列の対角成分の和をトレースという。

A_s = sym.Matrix([
    [1-I, 2,3+I],
    [  0, 1,  0],
    [  0, 0, -I]
])
A_s.trace()
A_n = np.array([
    [1-1j, 2, 3+1j],
    [   0, 1,    0],
    [   0, 0,  -1j]
])
A_n.trace()

$2−2i$

転置行列

A_s = sym.Matrix([
    [1, 2, 3],
    [0, 1, 0],
    [0, 0,-1]
])
A_s.T
A_n = np.array([
    [1, 2, 3],
    [0, 1, 0],
    [0, 0,-1]
])
A_n.T

$$A=\left[\begin{matrix}1 & 2 & 3\\0 & 1 & 0\\0 & 0 & -1\end{matrix}\right],{}^TA=\left[\begin{matrix}1 & 0 & 0\\2 & 1 & 0\\3 & 0 & -1\end{matrix}\right]$$

$A={}^TA$で要素がすべて実数となる行列を実対称行列という。

実対称行列の固有値は必ず実数で実直交行列$P$で対角化できる。
直交行列とは${}^TP=P^{-1}$となる行儀の良い行列。

複素共役

A_s = sym.Matrix([
    [1-I, 2,3+I],
    [  0, 1,  0],
    [  0, 0, -I]
])
A_sr = A_s.adjoint()
A_n = np.array([
    [1-1j, 2, 3+1j],
    [   0, 1,    0],
    [   0, 0,  -1j]
])
A_nr = A_n.conjugate()

$$\left[\begin{matrix}1 + i & 0 & 0\\2 & 1 & 0\\3 - i & 0 & i\end{matrix}\right]$$

随伴行列

A_s = sym.Matrix([
    [1-I, 2,3+I],
    [  0, 1,  0],
    [  0, 0, -I]
])
AH_s = A_s.H # 随伴行列

$A$を転置して複素共役をとった行列を随伴行列$A^*$とかく。
もとの行列と随伴行列が等しいものをエルミート行列という。

エルミート行列の固有値は必ず実数でユニタリ行列$U$で対角化できる。
ユニタリ行列とは$U^* U=UU^*=I$となる行列。

行列のn乗

$$J_{\lambda}(k)^n=
\left[
\begin{array}{cccc}
\lambda^n & { }_nC_1 \lambda^{n-1} & \ldots & { }_nC_k \lambda^{n-k} \\
0 & \lambda^n & \ldots & {}_n C_k-1 \lambda^{n-k+1} \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \ldots & \lambda^n
\end{array}\right]$$

$A=P^{-1}MPと書けるので$
$A^n=P^{-1}M^nP$

$A^n$が収束するには全ての固有値の絶対値が1より小さいか、
1の固有値がある場合はジョルダン細胞の最大次数が1。

A_s = sym.Matrix([
    [1,2, 0],
    [0,3, 1],
    [0,0,-1]
])
A_s ** n

$$\left[\begin{matrix}1 & 3^{n} - 1 & \frac{\left(-1\right)^{n}}{4} + \frac{3^{n}}{4} - \frac{1}{2}\\0 & 3^{n} & - \frac{\left(-1\right)^{n}}{4} + \frac{3^{n}}{4}\\0 & 0 & \left(-1\right)^{n}\end{matrix}\right]$$

A_n = np.array([
    [1,2, 0],
    [0,3, 1],
    [0,0,-1]
])
np.linalg.matrix_power(A_n, 10)

array([[    1, 59048, 14762],
       [    0, 59049, 14762],
       [    0,     0,     1]])

行列の指数関数

$$e^{A}=\sum_{k=0}^{\infty}\frac{A^k}{k!}$$

A_s = sym.Matrix([
    [1,2, 0],
    [0,3, 1],
    [0,0,-1]
])
A_s.exp()

$$\left[\begin{matrix}e & - e + e^{3} & - \frac{e}{2} + \frac{1}{4 e} + \frac{e^{3}}{4}\\0 & e^{3} & - \frac{1}{4 e} + \frac{e^{3}}{4}\\0 & 0 & e^{-1}\end{matrix}\right]$$

ある種の線形システムでは解が$x=e^{At}$の形で与えられる。
この系が$t \to \infty$で安定であるには$A$の固有値の実部が全て負である必要がある。
固有値を求めなくてもラウス・フルビッツの判定とかリアプノフ安定とかで判定できる。

Pythonで数学の勉強：線形代数学2(ジョルダン標準形、行列の指数関数)