Pythonで数学の勉強：線形代数学1(行列を扱う)

# 線形代数学1(行列を扱う)
## ベクトル・行列とコンピュータでの表現
正直大学での線形代数がよくわからなかった。
巨大な行列式を手計算させられたり
固有値だとか固有ベクトルだの訳の分からないものを求めさせられて
ジョルダン標準形にしろだとか何が楽しいのかと。
連立方程式を高速に解くための手法が発展して他分野にも応用できるようになった。
みたいな紹介をされていたら違ったのかもしれないが。

線形代数を始めるにあたってまずベクトルとは何かでふと止まる。
ベクトルと行列をコンピュータ上で区別する必要性は？
ベクトルが行列の特別な場合であると実装して不都合はあるのだろうか。
ベクトル空間は和が可換な群でスカラー倍が定義できてとかあるけど、
結局コンピュータ上で実現するには配列でしょってなるし。
ヒルベルト空間とか言われてもﾊｧですし。
そもそもヒルベルト空間では完備性が求められるけど
コンピュータ上では実数が表現できないことによるギャップはどう埋めるんだろう。
$\mathcal{Z}$の完備化とかp進数あたりを理解すると議論できるようになるのかな。
誰か詳しい人教えてくれ。

そこでベクトルにせよ行列にせよ数が格子状に並んだものとして定義しよう。
ここでいう数とは実数か複素数だ。
そうすれば可換性やスカラー倍は当然成り立つし。
完備性についてはとりあえず考えないことにしよう。
自分でオリジナルのベクトル空間を定義して使うときには
そのときに改めて内積とかノルムの概念が必要になってくるだろう。

## ベクトル・行列の基本操作

SympyとNumpyを並列して行う。
数式をそのまま扱うときはSympy。
Scipyで数値計算やるときはNumpyだ。

### 行列の生成

**一般の行列**
$$\\left[\\begin{matrix}1 & 2 & 3\\\\4 & 5 & 6\\end{matrix}\\right]$$

```py3
A1_s = sym.Matrix([
    [1,2,3],
    [4,5,6]
])
A1_n = np.array([
    [1,2,3],
    [4,5,6]
])
```
ソースコードは行列の形が見えるように書くべし。

**横ベクトル**
$$\\left[\\begin{matrix}3 & 4 & 5\\end{matrix}\\right]$$

```py3
A2_s = sym.Matrix([[3, 4, 5]])
A2_n = np.array([[3, 4, 5]])
```
「横ベクトルは行列の行を部分的に取り出したものである」
という視点を持つならばソースコードはこうあるべきか。
行列と同じくカッコを二重につける。

**縦ベクトル**
$$\\left[\\begin{matrix}-1\\\\-2\\\\-3\\end{matrix}\\right]$$

```py3
A3_s = sym.Matrix([[-1, -2, -3]]).T
A3_n = np.array([[-1, -2, -3]]).T
```
縦ベクトルを律儀に書くと[[-1], [-2], [-3]]で一重カッコ[-1,-2,-3]でもいける。
が、それはいかにもバグを内在させそうで危なっかしい。
なので横ベクトルを作ってから.Tで縦横入れ替えて縦ベクトルにするという方法を使う。
常に二重カッコの方が一貫性保てていいのではないでしょうか。

**sympy行列→numpy行列の変換**

```py3
A1_n == sym.matrix2numpy(A1_s)
```
```text
array([[ True,  True,  True],
       [ True,  True,  True]], dtype=bool)
```
行けるとこまでsympyで行って途中からnumpyでの数値計算に移りたいときなど。

### 零行列
$$O:\\left[\\begin{matrix}0 & 0 & 0\\\\0 & 0 & 0\\end{matrix}\\right]
\\left[\\begin{matrix}0 & 0 & 0\\\\0 & 0 & 0\\\\0 & 0 & 0\\end{matrix}\\right]
\\left[\\begin{matrix}0 & 0 & 0\\end{matrix}\\right]
\\left[\\begin{matrix}0\\\\0\\\\0\\end{matrix}\\right]$$

```py3
sym.Matrix.zeros(2, 3)
sym.Matrix.zeros(3)
sym.Matrix.zeros(1, 3)
sym.Matrix.zeros(3, 1)

np.zeros([2, 3])
np.zeros([3, 3])
np.zeros([1, 3])
np.zeros([3, 1])
```

要素が全て0の行列を作る。
sympyの方は正方行列のときは引数1つでOK。

### 単位行列
$$I=\\left[\\begin{matrix}1 & 0 & 0\\\\0 & 1 & 0\\\\0 & 0 & 1\\end{matrix}\\right]$$

```py3
sym.eye(3)
np.eye(3)
```
正方行列で対角成分が全て1で他が0の行列を作る。

### One
$$1 = \\left[\\begin{matrix}1 & 1 & 1\\\\1 & 1 & 1\\\\1 & 1 & 1\\end{matrix}\\right]$$

```py3
sym.ones(3)
np.ones(3)
```
全要素が1の正方行列を生成する。

### 対角行列
$$\\left[\\begin{matrix}1 & 0 & 0\\\\0 & 2 & 0\\\\0 & 0 & 3\\end{matrix}\\right]$$

```py3
sym.diag(1, 2, 3)
np.diag([1, 2, 3])
```
sympyは要素を対角に据えた複雑な行列も作れる。

```py3
A = sym.Matrix([
    [1,2,3],
    [4,5,6]
])
X = sym.diag(sym.eye(2), sym.ones(2), 3, A)
```
$$
\left[\begin{matrix}1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\
0 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\\\
0 & 0 & 1 & 1 & 0 & 0 & 0 & 0\\\\
0 & 0 & 1 & 1 & 0 & 0 & 0 & 0\\\\
0 & 0 & 0 & 0 & 3 & 0 & 0 & 0\\\\
0 & 0 & 0 & 0 & 0 & 1 & 2 & 3\\\\
0 & 0 & 0 & 0 & 0 & 4 & 5 & 6\end{matrix}\right]
$$
これならジョルダン標準形での初期化も一発ですね。

### 行列の要素へのアクセスと変更
$$\left[\begin{matrix}1 & 2 & 3\\\\4 & 5 & 6\end{matrix}\right]$$
4の部分を10に変更したいとする。

```py3
A1_s = sym.Matrix([
    [1,2,3],
    [4,5,6]
])
A1_n = np.array([
    [1,2,3],
    [4,5,6]
])
A1_s[1,0] #1, 0要素
A1_s[3] # 全体としてインデックスは3つ目
A1_s[1,0] = 10
A1_n[1,0]
A1_n[1,0] = 10
```

$$\left[\begin{matrix}1 & 2 & 3\\\\10 & 5 & 6\end{matrix}\right]$$
縦と横の2値を与えるのがどうみても無難。

部分行列の抽出や変更もできる。

```py3
A1_s[:,0:2] = sym.Matrix([
    [ 7, 8],
    [11,12]
])
A1_n[:,0:2] = np.array([
    [ 7, 8],
    [11,12]
])
```
:は全要素の抽出を表す。a:bはインデックスaから始まりb個の抽出。

$$\left[\begin{matrix}7 & 8 & 3\\\\11 & 12 & 6\end{matrix}\right]$$

### 行列の結合

```py3
A2_s = A1_s.row_join(sym.eye(2))
A3_s = A1_s.col_join(sym.eye(3))
A2_n= np.concatenate((A1_n, np.eye(2)), axis=1)
A3_n= np.concatenate((A1_n, np.eye(3)), axis=0)
```
$$A_2=\left[\begin{matrix}7 & 8 & 3 & 1 & 0\\\\11 & 12 & 6 & 0 & 1\end{matrix}\right],\ 
A_3=\left[\begin{matrix}7 & 8 & 3\\\\11 & 12 & 6\\\\1 & 0 & 0\\\\0 & 1 & 0\\\\0 & 0 & 1\end{matrix}\right]$$

### 行列の行または列の削除

```py3
A2_s.col_del(1)
A3_s.row_del(1)
A2_nr = np.delete(A2_n, 1, axis=1)
A3_nr = np.delete(A3_n, 1, axis=0)
```
$$\left[\begin{matrix}7 & 3 & 1 & 0\\\\11 & 6 & 0 & 1\end{matrix}\right], \ 
\left[\begin{matrix}7 & 8 & 3\\\\1 & 0 & 0\\\\0 & 1 & 0\\\\0 & 0 & 1\end{matrix}\right]$$
sympyの削除は破壊的。

### スカラーとの演算

```py3
sym.diag(2, 4, 6) / 2 + 10 * sym.ones(3)
A_n = np.array([
    [1,2,3],
    [4,5,6]
])
B_n = A1_n = np.array([
    [0,1,0],
    [1,0,1]
])
A_n * (-1) + 10 * B_n + 100
```
$$\left[\begin{matrix}11 & 10 & 10\\\\10 & 12 & 10\\\\10 & 10 & 13\end{matrix}\right]$$
array([[ 99, 108,  97],
       [106,  95, 104]])

スカラー倍は両者ともできる。スカラーとの加算はnumpyのみ可で全要素に加算される。
sympyで同じことやるにはonesにスカラー倍したのを足せばいいんじゃないかな。

### 関数の適用
要素ごとに関数を適用する。

```py3
A_s = sym.Matrix([
    [Pi/2, 3*Pi/4],
    [Pi/3,     Pi]
])
A_sr = A_s.applyfunc(lambda x: sym.sin(x))
A_n = np.array([
    [np.pi/2, 3*np.pi/4],
    [np.pi/3,     np.pi]
])
np.sin(A_n)
```
$$A_{sr} = \left[\begin{matrix}1 & \frac{\sqrt{2}}{2}\\\\
\frac{\sqrt{3}}{2} & 0\end{matrix}\right]$$
array([[  1.00000000e+00,   7.07106781e-01],
       [  8.66025404e-01,   1.22464680e-16]])
sympyの関数適用は破壊的ではないので受け皿が必要。
numpyは普通の変数と同じく関数に入れるだけっていう楽さだけど
ただのsin関数で早速浮動小数の誤差が出てる。

### 行列の積

```py3
A_s = sym.Matrix([
    [1,2,3],
    [4,5,6]
])
B_s = sym.Matrix([
    [1,2],
    [0,1],
    [0,1]
])
A_s * B_s
A_n = np.array([
    [1,2,3],
    [4,5,6]
])
B_n = np.array([
    [1,2],
    [0,1],
    [0,1]
])
A_n.dot(B_n)
```
$$AB=\left[\begin{matrix}1 & 7\\\\4 & 19\end{matrix}\right]$$

numpyはdotを使わないといけないのがちょっとねぇ。
あっちはすべての要素に同じ演算を適用させるのが基本。
つまりSIMDを意識してるってことで。