大好きな対称行列をsympyでスペクトル分解する

# はじめに

『秋葉原ロボット部 理論グループ Advent Calendar 2023』20日目の投稿です
線形代数が大好きで、スキが高じて量子力学の勉強会にも参加しています
線形代数のなかでも個人的に一番好きな『対称行列』について書いてみました
実数から複素数に世界を広げると、エルミート行列になります、エルミート行列は量子力学の中心的な役割を演じます

コメントなどございましたらありがたいです（お手柔らかに）


# 大好きな対称行列

この記事は、[大好きな対称行列をsympyで固有値分解する](https://qiita.com/osawat/items/154e185953cf1afcabf5)の続編です
固有値分解の先のスペクトル分解に話を展開します

## 対称行列の性質

対称行列は、以下のようなとても良い性質を持っています
対称行列のこれらの性質をsympyを使って検証していきます

話題の中心は行列の固有値分解で、対称行列がどのように振る舞うかが見どころです

1. 固有値はすべて実数
1. 正規直交する固有ベクトルを選ぶことができる
    - $S = [x_1, \ x_2, \ x_3]$
    - $x_1 \perp x_2, \ x_2 \perp x_3, \ x_3 \perp x_1$
1. その固有ベクトル行列$S$は、直交行列$Q$とすることができる
    - $Q = [q_1, \ q_2, \ q_3]$
    - $Q^{-1} = Q^T$
    - $Q^TQ=QQ^T = I \ (単位行列)$

4. スペクトル分解

```math
\begin{aligned}
A &= Q \Lambda Q^{-1} \\
&= Q \Lambda Q^T \\
&= \lambda_1 Q_1 + \lambda_2 Q_2 + \lambda_3 Q_3 \ \ \ (Q_i = q_i q_i^T)\\
&= \lambda_1 q_1 q_1^T + \lambda_2 q_2 q_2^T + \lambda_3 q_3 q3^T
 \end{aligned}
```

## とにかく、sympyで試してみる

sympyはPythonで使える数式を扱えるライブラリです
行列の計算はとても手間がかかります、人生を無駄使いしないようツールを使いましょう
jupyter notebookベースで使うことを前提としてsympyでコードを記述します

ここから、以下の$3 \times 3$の対称行列について考察します
(※ちなみに、この行列は『量子力学10講(谷村省吾)』の(3.101) P51に登場する対称行列で、3つの固有値のうち、2つが重解となるものです)

```math
A = \begin{bmatrix} 3 & -1 & -1 \\
-1 & 3 & -1 \\
-1 & -1 & 3 \end{bmatrix}
```
# 初期設定
jupyter notebookでsympyを使うときのお約束です
```py
# 初期設定
from sympy import *
init_printing()
```

# 対称行列の定義

```py
A = Matrix([
    [3, -1, -1],
    [-1,  3, -1],
    [-1, -1,  3],
])
A
```

```math
\begin{bmatrix} 3 & -1 & -1 \\
-1 & 3 & -1 \\
-1 & -1 & 3 \end{bmatrix}
```


## 対称行列の固有値分解


3x3行列の対角化についての概要

```math
\begin{array}{l}
A\vec{x_1} = \lambda \vec{x_1} & (1) \\
A\vec{x_2} = \lambda \vec{x_2} & (2) \\
A\vec{x_3} = \lambda \vec{x_3} & (3) \\
\\ 

(1),(2),(3)は以下のようにまとめられる\\
A \underbrace{[ \vec{x_1} \ \vec{x_2} \ \vec{x_3}}_{S}] =  \underbrace{[ \vec{x_1} \ \vec{x_2} \ \vec{x_3} ]}_{S}
\underbrace{ \begin{bmatrix}
\lambda_1  & 0 & 0 \\
0 & \lambda_2 & 0 \\
0 & 0 & \lambda_3
\end{bmatrix}}_{\Lambda} \\
AS = S \Lambda \\
\\
これから \\
A = S \Lambda S^{-1} \ \ (固有値分解) \\
\Lambda = S^{-1} AS \ \ (対角化)\\
\\
\end{array}
```

# sympyで固有値分解してみる

行列`A`に対して`.diagonalize()`メソッドを使って固有値分解します
- 戻り値
    - 行列`S`：固有ベクトルが列ベクトルとなっている行列
        - 特に引数として`normalize=True`とすると固有ベクトルが正規化される
    - 行列`Lambdas`： 対角成分が固有値となっている対角行列

```py
# 固有値分解
S, Lambdas = A.diagonalize(normalize=True)
S, Lambdas # 固有ベクトル、固有値ベクトル
```

```math
\left(
\begin{bmatrix}
\frac{\sqrt{3}}{3} & -\frac{\sqrt{2}}{2} & -\frac{\sqrt{2}}{2} \\
\frac{\sqrt{3}}{3} & \frac{\sqrt{2}}{2}, &0 \\
\frac{\sqrt{3}}{3} &  0 & \frac{\sqrt{2}}{2} \\
\end{bmatrix}
 \ , \
\begin{bmatrix}
1 & 0 & 0 \\
0 & 4 & 0 \\
0 & 0 & 4 \\
\end{bmatrix}
\right)
```

## ちゃんと固有値分解されているか復元して確認

$A=S\Lambda S^{-1}$として元の$A$になるかを確認 ($AS=S\Lambda$)

```py
# 検算 : 元通りになるか確認
S * Lambdas * S.inv()
```
```math
\begin{bmatrix} 3 & -1 & -1 \\
-1 & 3 & -1 \\
-1 & -1 & 3 \end{bmatrix}
```
- 復元できました

## 結果確認

- 固有値は`Lambdas`の対角成分で、$1, 4, 4$ であり$4$が重解となっています
- 固有ベクトルは固有ベクトル行列`S`のそれぞれの列ベクトルです
    - ここでは正規化するよう`normalize=True`としたのでノルムは$1$です

```py
# 固有ベクトルと固有値を取り出す
x1, x2, x3 = S[:,0] , S[:,1] , S[:,2]  # 固有ベクトル
lambda1, lambda2, lambda3 = Lambdas.diagonal()  # 固有値
x1, x2, x3, lambda1, lambda2, lambda3
```

```math
\left(
\begin{bmatrix}
\frac{\sqrt{3}}{3} \\
\frac{\sqrt{3}}{3} \\
\frac{\sqrt{3}}{3} \\
\end{bmatrix}
 \ , \
\begin{bmatrix}
-\frac{\sqrt{2}}{2} \\
\frac{\sqrt{2}}{2} \\
0  \\
\end{bmatrix}
 \ , \
\begin{bmatrix}
-\frac{\sqrt{2}}{2} \\
0 \\
\frac{\sqrt{2}}{2} \\
\end{bmatrix}
\ , \  1 \ , \ 4  \ , \ 4 
\right)
```


## ちゃんと固有ベクトルが直交しているか確認、アレッ！？

固有ベクトルが直交しているかどうかは、内積をとればわかります

```py
# 内積
# (x1,x2), (x1,x3), (x2,x3)
x1.T*x2,  x1.T*x3, x2.T*x3
```

```math
\left([0],[0],[\frac{1}{2}]\right)
```

上記から、以下のことがわかりました

- $x1 \perp x2 \ , \ x1 \perp x3$
- $x2 \not \perp x3$

つまり、$x2$ と $x3$は直交していないってことです
理由は、これらの固有値は重解で同じ$4$であったからですが・・・
- 対称行列の固有ベクトルは
    - **固有値が異なれば必ず直交します**
    - 固有値が同じものでも直交するように固有ベクトルを選ぶことができます


<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/71601/8b912b71-a746-295a-22a5-284f812cc285.jpeg" width="60%">

- 固有値が同じものの固有ベクトルは、ひと手間いるようです
    - つまり、上図のように$x_1$は、**$x_2$と$x_3$で張られる面**に直交していますが
    - $x_2$と$x_3$どうしは直交していません
    - ただし、固有値が$\lambda=4$の固有ベクトルは**この面上**の2つのベクトルを選べばよいので
        - グラム・シュミットの直交化により直交するようにします

<hr>
また、$\lambda=4$を固有値を持つような固有ベクトル、以下のように計算されます

```math
\begin{aligned}
(A-\lambda I)z &= 0 \\
\left(
\begin{bmatrix} 3 & -1 & -1 \\
-1 & 3 & -1 \\
-1 & -1 & 3 \end{bmatrix}
-
\begin{bmatrix} 4 & 0 & 0 \\
0 & 4 & 0 \\
0 & 0 & 4 \end{bmatrix}
\right)
\begin{bmatrix} z_1\\z_2\\z_3\end{bmatrix}
&=0\\
\begin{bmatrix} -1 & -1 & -1 \\
-1 & -1 & -1 \\
-1 & -1 & -1 \end{bmatrix}
\begin{bmatrix} z_1\\z_2\\z_3\end{bmatrix}
&=0 \\
\\
\therefore z_1 + z_2 + z_3 = 0 \ \ (平面の方程式)
\end{aligned}
```
上記方程式を満たすように2つの正規直交ベクトルを選べばよいのです


# 同じ固有値の固有ベクトルを直交化する

## グラム・シュミットの直交化によりx2, x3を直交化

sympyには、そのものズバリの関数がありますので、ありがたく使わせてもらいます

```py
# x2, x3を正規直交化 (グラム・シュミットの正規直交化)
x2, x3 = GramSchmidt([x2, x3], True)
x2, x3
```

```math
\left(
\begin{bmatrix}
\frac{\sqrt{2}}{2} \\
\frac{\sqrt{2}}{2} \\
0 \\
\end{bmatrix}
 \ , \
\begin{bmatrix}
-\frac{\sqrt{6}}{6} \\
-\frac{\sqrt{6}}{6} \\
\frac{\sqrt{6}}{3} \\
\end{bmatrix}
\right)
```

## ちゃんと直交するか確認
```py
# 内積
# (x1,x2), (x1,x3), (x2,x3)
x1.T*x2,  x1.T*x3, x2.T*x3
```

```math
\left( [0], [0], [0]\right)
```

- 直交してますね

# 正規直交化された固有ベクトル

正規直交する列ベクトルをまとめて固有ベクトル行列を作成する

```py
## 正規直交する列ベクトルを行列にまとめる
S = Matrix.hstack(v1,v2,v3)
S
```

```math
\begin{bmatrix}
\frac{\sqrt{3}}{3} & -\frac{\sqrt{2}}{2} & -\frac{\sqrt{6}}{6} \\
\frac{\sqrt{3}}{3} & \frac{\sqrt{2}}{2}, & -\frac{\sqrt{6}}{6} \\
\frac{\sqrt{3}}{3} &  0 & \frac{\sqrt{6}}{6} \\
\end{bmatrix}
```

## ちゃんと固有値分解されているか復元して確認

$A=S\Lambda S^{-1}$からちゃん元の$A$になるかを確認 ($AS=S\Lambda$)

```py
# 検算 : 元通りになるか確認
S * Lambdas * S.inv()
```


```math
\begin{bmatrix} 3 & -1 & -1 \\
-1 & 3 & -1 \\
-1 & -1 & 3 \end{bmatrix}
```

- 復元しました

# 直交行列

対称行列の固有ベクトルは直交行列とすることができます
- 上記の固有ベクトル行列$S$は直交行列となっています
- 以降、直交行列らしく$Q$として扱います（$Q=S$）

## 直交行列とは

直交行列は以下の性質を持ちます

- $Q^{-1} = Q^T$
- $Q^TQ = QQ^T = I$

## 固有ベクトル行列が直交行列となのか確認する

 $Q^{-1} = Q^T$となるか

```py
Q = S # SをQとおく
Q.inv() == Q.T
```

- なりました

```math
True
```

$Q^TQ = Q^TQ = I$となるか

```py
Q*Q.T, Q*Q.T
```

```math
\left(
\begin{bmatrix} 1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1 \end{bmatrix} \ , \
\begin{bmatrix} 1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1 \end{bmatrix}
\right)
```

- なりました

# スペクトル分解

## スペクトル分解とは

スペクトル分解のひとつの見方は、独立な射影の線形結合とも言えます（と思う）
つまり、もとものの行列の振る舞いを射影の線形結合として表現する
(※後で、見直す)

```math
\begin{aligned}
A &= Q \Lambda Q^{T} \\
&=
Q
\begin{bmatrix} 
\lambda_1 & 0 & 0\\0 & \lambda_2 & 0 \\0 & 0 & \lambda_3\end{bmatrix}
Q^{T} \\
&=
Q
\left(
\begin{bmatrix} \lambda_1 & 0 & 0\\0 & 0 & 0\\0 & 0 & 0\end{bmatrix}
+
\begin{bmatrix} 0 & 0 & 0\\0 & \lambda_2 & 0\\0 & 0 & 0\end{bmatrix}
+
\begin{bmatrix} 0 & 0 & 0\\0 & 0 & 0\\0 & 0 & \lambda_3\end{bmatrix}
\right)
Q^{T} \\
&=
\lambda_1 Q\begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}Q^{T}
+
\lambda_2 Q\begin{bmatrix} 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{bmatrix}Q^{T}
+
\lambda_3 Q\begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix}Q^{T} \\
&=
\lambda_1  \begin{bmatrix} q_1 & q_2 & q_3 \end{bmatrix}
\begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}
\begin{bmatrix} q_1^T \\ q_2^T \\ q_3^T \end{bmatrix}
+
\lambda_2  \begin{bmatrix} q_1 & q_2 & q_3 \end{bmatrix}
\begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}
\begin{bmatrix} q_1^T \\ q_2^T \\ q_3^T \end{bmatrix}
+
\lambda_3  \begin{bmatrix} q_1 & q_2 & q_3 \end{bmatrix}
\begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}
\begin{bmatrix} q_1^T \\ q_2^T \\ q_3^T \end{bmatrix}\\
&= \lambda_1 q_1 q_1^T + \lambda_2 q_2 q_2^T + \lambda_3 q_3 q_3^T \\
&= \lambda_1 Q_1 + \lambda_2 Q_2 + \lambda_3 Q_3 \ \ (\because Q_i = q_iq_i^T)
\end{aligned}
```

## sympyでスペクトル分解してみる

```py
q1, q2, q3 = Q[:,0], Q[:,1], Q[:,2]
q1, q2, q
```

```math
\left(
\begin{bmatrix}
\frac{\sqrt{3}}{3} \\
\frac{\sqrt{3}}{3} \\
\frac{\sqrt{3}}{3} \\
\end{bmatrix}
 \ , \
\begin{bmatrix}
\frac{\sqrt{2}}{2} \\
\frac{\sqrt{2}}{2} \\
0 \\
\end{bmatrix}
 \ , \
\begin{bmatrix}
-\frac{\sqrt{6}}{6} \\
-\frac{\sqrt{6}}{6} \\
\frac{\sqrt{6}}{3} \\
\end{bmatrix}
\right)
```

```py
Q1, Q2, Q3 = q1*q1.T, q2*q2.T, q3*q3.T
Q1, Q2, Q3
```

```math
\left(
\begin{bmatrix} \frac{1}{3} & \frac{1}{3} & \frac{1}{3} \\
\frac{1}{3} & \frac{1}{3} & \frac{1}{3} \\
\frac{1}{3} & \frac{1}{3} & \frac{1}{3} \end{bmatrix} 
\ , \
\begin{bmatrix} \frac{1}{2} & -\frac{1}{2} & 0 \\
-\frac{1}{2} & \frac{1}{2} & 0 \\
0 & 0 & 0 \end{bmatrix} 
\ , \
\begin{bmatrix} \frac{1}{6} & \frac{1}{6} & -\frac{1}{3} \\
\frac{1}{6} & \frac{1}{6} & -\frac{1}{3} \\
-\frac{1}{3} & -\frac{1}{3} & \frac{2}{3} \end{bmatrix} 
\right)
```

## スペクトル分解したものがもとに復元できるか

$A = \lambda_1 Q_1 + \lambda_2 Q2 + \lambda_3 Q3$

```py
lambda1 * Q1 + lambda2 * Q2 + lambda3 * Q3
```
```math
\begin{bmatrix} 3 & -1 & -1 \\
-1 & 3 & -1 \\
-1 & -1 & 3 \end{bmatrix}
```

- 復元できた

## 射影行列

射影行列の特徴

- $P=P^2$
- $P=P^T$

```py
Q1*Q1==Q1, Q2*Q2==Q2, Q3*Q3==Q3 
```

```math
(True, True, True)
```

```py
Q1==Q1.T, Q2==Q2.T, Q3==Q3.T 
```

```math
(True, True, True)
```

## 直和

```py
Q1 + Q2 + Q3
```

```math
\begin{bmatrix} 1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1 \end{bmatrix}
```

```py
Q1*Q2, Q1*Q3, Q2*Q3
```

```math
\left(
\begin{bmatrix} 0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0 \end{bmatrix}
\,\
\begin{bmatrix} 0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0 \end{bmatrix}
\,\
\begin{bmatrix} 0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0 \end{bmatrix}
\right)
```

# まとめ

対称行列は本当にいい仕事をしてくれます
この投稿をアップデートしながら伝えていけたらと思います

<hr>

参考文献

-『ストラング線形代数イントロダクション』
-『量子力学10講』