大好きな射影行列をsympyで記述してみる

# はじめに

『秋葉原ロボット部 理論グループ Advent Calendar 2023』6日目の投稿です
線形代数が大好きで、スキが高じて量子力学の勉強会にも参加しています
線形代数のなかでも個人的に一番好きな『射影行列』について書いてみました
コメントなどございましたらありがたいです（お手柔らかに）

# ベクトル **b** を **a** ベクトルに正射影

ベクトル$\vec{b}$をベクトル$\vec{a}$に直交するように射影します
射影後のベクトルを$\vec{p}$とします
(ベクトル$\vec{a}$に直角に光線をあてて$\vec{b}$の影を投影するのに例えられると思う)


![スクリーンショット 2023-12-06 000748.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/71601/3e640404-f32b-8305-db27-9fb5af570830.png)


## 射影行列を考える

- ベクトル$\vec{p}$は、図からからわかるように、
    - ベクトル$\vec{a}$と並行で
    - その大きさは、$\vec{a}$の$x$倍(スカラー)だと考えられます


$$
\begin{array}{l}
\vec{p} = x \vec{a}
\end{array}
$$


- ここでポイントは、ベクトル$a$とベクトル$\vec{p} -\vec{b}$が直交していることです
    - 2つのベクトルが直交するということは、その内積は0になりますので
       - これから$x$を求めることができます

```math
\begin{array}{l}
\vec{a} \perp (\vec{p}-\vec{b}) \ (直交)\Rightarrow \ \vec{a}^T(\vec{p}-\vec{b}) = \vec{0} \ (内積=0)\\
\vec{a}^T(\vec{p}-\vec{b}) = \vec{a}^T(x\vec{a}-\vec{b}) = x \vec{a}^T \vec{a} -\vec{a}^T \vec{b} = \vec{0} \\
x \ \vec{a}^T\vec{a} = \vec{a}^T \vec{b} \\
x = \frac{ \vec{a}^T \vec{b}}{\vec{a}^T\vec{a}} \\
\therefore \vec{p} = x \vec{a} =  \frac{ \vec{a}^T \vec{b}}{\vec{a}^T\vec{a}} \vec{a}
\end{array}
```

ここから本番です！
ここで、少しだけ上記を変形します

```math
\begin{array}{l}
\vec{p} = \left( \frac{ \vec{a}^T \vec{b}}{\vec{a}^T\vec{a}} \right) \vec{a}
= \left(\frac{ \vec{a}\vec{a}^T}{\vec{a}^T\vec{a}} \right) \vec{b} = P \vec{x} \\
\because P = \frac{ \vec{a}\vec{a}^T}{\vec{a}^T\vec{a}} \ \ (射影行列)
\end{array}
```
上記の$P$は射影行列となります

どのようなベクトル$b$を持ってきても、射影行列$P$を作用させれば、たちどころにベクトル$a$に射影されたベクトルが得られます (何でも射影できる$P$はカッコイイと思います)

## sympyで試してみる

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

```py
# 初期設定
from sympy import *
init_printing()

# ベクトルの定義
a = Matrix([1,1])
b = Matrix([1,0])
a, b
```

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

```py
# 射影行列
P = a*a.T/(a.T*a)[0]
P
```

```math
\begin{bmatrix} \frac{1}{5} & \frac{2}{5} \\ \frac{2}{5} & \frac{4}{5} \end{bmatrix}
```

```py
# 射影
p = P*b
p
```


```math
\begin{bmatrix} \frac{1}{5}  \\ \frac{2}{5} \end{bmatrix}
```
# ベクトル **b** を **a1, a2** ベクトルが張る平面に正射影

上記は直線への射影でしたが、ここでは平面への射影を考えます
ベクトル$\vec{b}$をベクトル$\vec{a_1}, \vec{a_2}$によって張られる平面に直交するように射影します
射影後のベクトルを$\vec{p}$とします
($\vec{a_1},\vec{a_2}$からなる平面に直角に光線をあてて$\vec{b}$の影を投影するのに例えられると思う)

![スクリーンショット 2023-12-06 001113.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/71601/bf9496de-8397-ba9f-b571-429ce4c25a56.png)


## 射影行列を考える

- ベクトル$\vec{p}$は
    - ベクトル$\vec{a_1}$とベクトル$\vec{a_2}$の線形結合で表せる

```math
\begin{align}
\vec{p} &= x_1 \vec{a_1} + x_2 \vec{a_2}　\\
&= \begin{bmatrix}\vec{a_1} & \vec{a_2}\end{bmatrix} \begin{bmatrix}x_1 \\ x_2\end{bmatrix} \\
&= A\vec{x}
\end{align}
```


ここでのポイントは、ベクトル$\vec{a_1}$とベクトル$\vec{a_2}$は、それぞれベクトル$\vec{p} -\vec{b}$と直交していることです
ベクトルが直交するということは、その内積は0になりますので

```math
\begin{array}{l}
\begin{cases}
\ \vec{a_1} \perp (\vec{p}-\vec{b}) \ (直交)\Rightarrow \ \vec{a_1}^T(\vec{p}-\vec{b}) = \vec{0} \ (内積=0)\\
\ \vec{a_2} \perp (\vec{p}-\vec{b}) \ (直交)\Rightarrow \ \vec{a_2}^T(\vec{p}-\vec{b}) = \vec{0} \ (内積=0)\\
\end{cases} \\
\begin{cases}
\ \vec{a_1}^T(\vec{p}-\vec{b}) = \vec{a_1}^T(A\vec{x}-\vec{b}) =  \vec{0} \\
\ \vec{a_2}^T(\vec{p}-\vec{b}) = \vec{a_2}^T(A\vec{x}-\vec{b}) =  \vec{0} \\
\end{cases} \\
2つの式を一つにまとめて\\
\begin{bmatrix} \vec{a_1}^T \\ \vec{a}^T \end{bmatrix} (A\vec{x}-\vec{b}) = A^T  (A\vec{x}-\vec{b})= A^TA \vec{x} - A^Tb = \vec{0} \\
A^TA \vec{x} = A^Tb\\
\vec{x} = (A^TA)^{-1}A^Tb\\
\therefore \vec{p} = A\vec{x} = A(A^TA)^{-1}A^Tb = P\vec{b}\\
\because P =  A(A^TA)^{-1}A^T \left( 射影行列\right) 
\end{array}
```

上記の$P$は射影行列となります

どのようなベクトル$b$を持ってきても、射影行列$P$を作用させれば、たちどころにベクトル$a_1, a_2$の張る平面に射影されたベクトルが得られます

## sympyで試してみる

超簡単な例ですが

```py
# ベクトル定義
a1 = Matrix([1,0,0])
a2 = Matrix([1,1,0])
b = Matrix([2,4,6])
a1, a2, b
```


```math
\left(
\begin{bmatrix}1 \\ 0 \\ 0 \end{bmatrix}
\begin{bmatrix}0 \\ 1 \\ 0 \end{bmatrix}
\begin{bmatrix}2 \\ 4 \\ 6 \end{bmatrix}
\right)
```


```py
# 列ベクトルa1, a2をAとして並べる
A = Matrix.hstack(a1, a2)
A
```

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

```py
# 射影行列
P = A*(A.T*A).inv()*A.T
P
```

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

```py
# 射影
p = P*b
p
```

```math
\begin{bmatrix} 2 \\ 4 \\  0 \end{bmatrix}
```

# まとめ

直線への射影行列と平面への射影行列はとてもよく似ています


- 直線への射影 $P = \frac{ \vec{a}\vec{a}^T}{\vec{a}^T\vec{a}}$
- 平面への射影 $P =  A(A^TA)^{-1}A^T$

直線の射影行列の式から平面への射影行列の式が類推できますね

- 行列に分数はありませんから、分母は$-1$乗されています(逆行列)
- 分数の部分はちょうど$AA^T$の真ん中に割り込むように配置されています

<hr>

参考文献

『ストラング線形代数イントロダクション』