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

# 線形代数学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$の列ベクトルのうち線形独立なものの個数にも等しい。

```py3
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$なら次元が潰れているということ。

``` py3
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}
$$

```py3
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$

```py3
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]$$

```py3
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で逆行列を計算してみよう。

```py3
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次実正方行列の固有値に複素数が現れるときは共役複素数も固有値となる。

```py3
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]

```py3
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\]
$$

```py3
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次以上の部分に相当する。

```py3
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 )$$

## 種々の行列の性質

### トレース

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

```py3
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$

### 転置行列

```py3
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}$となる行儀の良い行列。

### 複素共役

```py3
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]$$

### 随伴行列

```py3
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。

```py3
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]$$

```py3
A_n = np.array([
    [1,2, 0],
    [0,3, 1],
    [0,0,-1]
])
np.linalg.matrix_power(A_n, 10)
```
```text
array([[    1, 59048, 14762],
       [    0, 59049, 14762],
       [    0,     0,     1]])
```

## 行列の指数関数

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

```py3
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$の固有値の実部が全て負である必要がある。
固有値を求めなくてもラウス・フルビッツの判定とかリアプノフ安定とかで判定できる。