【Sympy演算処理】線形代数再入門①「複式簿記」から「単位行列=単位ベクトル」へ。

**連立一次方程式**も**行列演算**も、それ自体はコンピューターに丸投げすれば一瞬で解いてくれます。むしろその時何が起こってるかちゃんとイメージするのが難しいという…

# 1.行列演算の出発点としての「複式簿記」。
[【Python演算処理】行列演算の基本④大源流における記述統計学との密接な関連性?](https://qiita.com/ochimusha01/items/d36169ad13872a34ebc6)
「**とある喫茶店におけるコーヒーと紅茶の1週間の売り上げ**」から出発。

<table border="1" class="dataframe">  <thead>    <tr style="text-align: right;">      <th></th>      <th>Coffe</th>      <th>Tea</th>      <th>Total</th>    </tr>  </thead>  <tbody>    <tr>      <th>0</th>      <td>18</td>      <td>10</td>      <td>8900</td>    </tr>    <tr>      <th>1</th>      <td>4</td>      <td>5</td>      <td>2950</td>    </tr>    <tr>      <th>2</th>      <td>41</td>      <td>6</td>      <td>14400</td>    </tr>    <tr>      <th>3</th>      <td>49</td>      <td>23</td>      <td>22750</td>    </tr>    <tr>      <th>4</th>      <td>4</td>      <td>9</td>      <td>4350</td>    </tr>    <tr>      <th>5</th>      <td>41</td>      <td>6</td>      <td>14400</td>    </tr>    <tr>      <th>6</th>      <td>33</td>      <td>35</td>      <td>22150</td>    </tr>  </tbody></table>

①これをそのまま行列表現したのが**拡大係数行列**(Augmented Coefficient Matrix)」となる。

```python
import sympy as sp
import numpy as np
import pandas as pd
a,b,x,y = sp.symbols('a,b,x,y')
A = sp.Matrix([[300,350]])
X = sp.Matrix([[18,4,41,49,4,41,33],[10,5,6,23,9,6,35]])
X1=X.row_insert(2, A*X)
x=X1.transpose()
df=pd.DataFrame(np.matrix(x),columns=['Coffe', 'Tea', 'Total'])
sp.init_printing()
display(A)
display(X) 
display(A*X)
print(sp.latex(Aa)+sp.latex(X)+"="+sp.latex(Aa*X))
display(x)
org=df.to_html()
print(org.replace('\n', ''))
```

```math
\left[\begin{matrix}18 & 10 & 8900\\4 & 5 & 2950\\41 & 6 & 14400\\49 & 23 & 22750\\4 & 9 & 4350\\41 & 6 & 14400\\33 & 35 & 22150\end{matrix}\right]
```
②ここでコーヒーの単価は300円、紅茶は350円。毎日の売り上げを行列演算結果に対応させると以下となる。

```python
import sympy as sp
a,b,x,y = sp.symbols('a,b,x,y')
A = sp.Matrix([[a,b]])
X = sp.Matrix([[x],[y]])
AX=A*X
expr00 =AX[0]
Aa=A.subs([(a,300),(b,350)])
Xr=sp.randMatrix(2,7, min=0, max=50)
sp.init_printing()
display(Aa)
display(X) 
display(Aa*X)
print(sp.latex(Aa)+sp.latex(X)+"="+sp.latex(Aa*X))
display(Xr) 
display(Aa*Xr)
print(sp.latex(Aa)+sp.latex(Xr)+"="+sp.latex(Aa*Xr))
```

```math
\left[\begin{matrix}300 & 350\end{matrix}\right]\left[\begin{matrix}x\\y\end{matrix}\right]=\left[\begin{matrix}300 x + 350 y\end{matrix}\right] \\
\left[\begin{matrix}300 & 350\end{matrix}\right]\left[\begin{matrix}18 & 4 & 41 & 49 & 4 & 41 & 33\\10 & 5 & 6 & 23 & 9 & 6 & 35\end{matrix}\right]=\left[\begin{matrix}8900 & 2950 & 14400 & 22750 & 4350 & 14400 & 22150\end{matrix}\right]
```
③ここから**複式簿記**(Double-entry bookkeeping)と**一次方程式**(Linear Equation)と**行列演算**(Matrix Operation)の関係が浮かび上がってくる。

```python
import sympy as sp
a,b,c,x,y = sp.symbols('a,b,c,x,y')
A = sp.Matrix([[a,b,c]])
X = sp.Matrix([[x],[y],[-1]])
AX=A*X
expr00 =AX[0]
expr01 = sp.solve(expr00,y)
expr02 = sp.solve(expr00,x)
expr03 = sp.solve(expr00,c)
Aa=A.subs([(a,300),(b,350)])
expr00a = expr00.subs([(a,300),(b,350)])
expr01a = expr01[0].subs([(a,300),(b,350)])
expr02a = expr02[0].subs([(a,300),(b,350)])
expr03a = expr03[0].subs([(a,300),(b,350)])
sp.init_printing()
display(Aa)
display(X) 
display(Aa*X)
print(sp.latex(Aa)+sp.latex(X)+"="+sp.latex(Aa*X))
display(expr00a)
print(sp.latex(expr00a)+"=0")
display(expr01a)
print("y="+sp.latex(expr01a))
display(expr02a)
print("x="+sp.latex(expr02a))
display(expr03a)
print(sp.latex(expr03a)+"=c")
```

```math
\left[\begin{matrix}300 & 350 & c\end{matrix}\right]\left[\begin{matrix}x\\y\\-1\end{matrix}\right]=\left[\begin{matrix}- c + 300 x + 350 y\end{matrix}\right] \\
- c + 300 x + 350 y=0 \\
y=\frac{c}{350} - \frac{6 x}{7} \\
x=\frac{c}{300} - \frac{7 y}{6} \\
300 x + 350 y=c
```
(拡大係数行列の部分集合として現れる)連立一次方程式の係数を並べた**係数行列**(Coefficient Matrix)を**小行列**(Minor Determinant)と呼びます。それは(自分自身との直積としての)**正方行列**(Square Matrix=列数と行数が等しい行列)を構成し、**単位行列**(Identity Matrix)」と呼ばれたりもする訳ですね。それが**複式簿記**でいうと適切な個数の売り上げデータからの逆算で**単価表**を作成する事に対応する訳です。

```math
\left[\begin{matrix}1 & 0 & 300\\0 & 1 & 350\\300 & 350 & -\end{matrix}\right]
```
この様に**連立一次方程式**や**行列演算**を解くとは、まさに**対角成分**(Diagonal Element)のみを抽出しそれを求める過程に他なりません。ちなみに残りを**非対角成分**(Nondiagonal Element)と呼びます。**成分**(Elemennt)も**元**(Elemennt)も英語では同じ言葉。要するに数学ではこれらを同じ集合論上の概念として扱うのです。

# 2.連立一次方程式と行列演算の関係
[行列の定義と演算規則](https://python.atelierkobato.com/matrix-vector/)

①Sympyで連立一次方程式を解く

```python
import sympy as sp
a,b,c,d,e,f,x,y = sp.symbols('a,b,c,d,e,f,x,y')
eq1=sp.Eq(a*x+b*y, e)
eq2=sp.Eq(c*x+d*y, f)
expr01 =sp.solve ([eq1, eq2], [x, y]) 
#表示
sp.init_printing()
display(eq1)
display(eq2) 
display(expr01)
print(sp.latex(expr01))
```

```math
a x + b y = e\\
c x + d y = f\\
\left\{ x : \frac{- b f + d e}{a d - b c}, \  y : \frac{a f - c e}{a d - b c}\right\}
```
ここで分母に現れるad-bcを**行列式**(Determinant)といい、その解が0になる時2式は単なる比例関係を表すだけとなって解が無数に存在する事になり(不定)、マイナスになる時は解が「反転」します。そして計算上矛盾が生じる時、解は存在しません(不能)。

「**不定**」の例(比例式が答えとなる)

```python
import sympy as sp
a,b,c,d,e,f,x,y = sp.symbols('a,b,c,d,e,f,x,y')
eq10=sp.Eq(a*x+b*y, e)
eq20=sp.Eq(c*x+d*y, f)
#代入
eq11=eq10.subs(([(a, 1), (b, -2),(e,3)]))
eq21=eq20.subs(([(c, 3), (d, -6),(f,9)]))
#解法
expr01 =sp.solve ([eq11, eq21], [x, y]) 
#表示
sp.init_printing()
display(eq11)
display(eq21) 
display(expr01)
print(sp.latex(expr01))
```

```math
x - 2 y = 3\\
3 x - 6 y = 9\\
\left\{ x : 2 y + 3\right\}
```
「**不能**」の例(空集合[]が答えとなる)

```python
import sympy as sp
a,b,c,d,e,f,x,y = sp.symbols('a,b,c,d,e,f,x,y')
eq10=sp.Eq(a*x+b*y, e)
eq20=sp.Eq(c*x+d*y, f)
#代入
eq11=eq10.subs(([(a, 1), (b, -2),(e,3)]))
eq21=eq20.subs(([(c, 3), (d, -6),(f,10)]))
#解法
expr01 =sp.solve ([eq11, eq21], [x, y]) 
#表示
sp.init_printing()
display(eq11)
display(eq21) 
display(expr01)
print(sp.latex(expr01))
```

```math
x - 2 y = 3\\
3 x - 6 y = 10\\
\left[ \right]
```
②行列演算との対応

```python
import sympy as sp
a,b,c,d,x,y = sp.symbols('a,b,c,d,x,y')
mt1= sp.Matrix([[a,b],[c,d]])
mt2= sp.Matrix([[x],[y]])
expr01 =mt1*mt2 
#表示
sp.init_printing()
display(mt1)
display(mt2) 
display(expr01)
print(sp.latex(mt1)+sp.latex(mt2)+"="+sp.latex(expr01))
```

```math
\left[\begin{matrix}a & b\\c & d\end{matrix}\right]\left[\begin{matrix}x\\y\end{matrix}\right]=\left[\begin{matrix}a x + b y\\c x + d y\end{matrix}\right]
```
ここから**単位行列**に到達する為に「**現在の係数行列にそれを掛けると単位行列となる**」**逆行列**(Inverse Matrix)の概念を導入。多くの線形代数入門書ではここで最初の難所を迎えるが「**面倒な計算はコンピューターに丸投げする**」式の考え方ではどうという事もない。

```python
import sympy as sp
a,b,c,d,x,y = sp.symbols('a,b,c,d,x,y')
mt1= sp.Matrix([[a,b],[c,d]])
mt2= sp.Matrix([[x],[y]])
expr01 =mt1**(-1)
#表示
sp.init_printing()
display(mt1)
display(expr01)
print(sp.latex(mt1)+"^{-1}="+sp.latex(expr01))
```

```math
\left[\begin{matrix}a & b\\c & d\end{matrix}\right]^{-1}=\left[\begin{matrix}\frac{d}{a d - b c} & - \frac{b}{a d - b c}\\- \frac{c}{a d - b c} & \frac{a}{a d - b c}\end{matrix}\right]
```
上掲の2つの式を合算すると同じ答えに到達する。

```python
import sympy as sp
a,b,c,d,x,y = sp.symbols('a,b,c,d,x,y')
mt1= sp.Matrix([[a,b],[c,d]])
mt2= sp.Matrix([[x],[y]])
expr01 =mt1**(-1)*mt2
#表示
sp.init_printing()
display(mt1)
display(expr01)
print(sp.latex(mt1)+"^{-1}"+sp.latex(mt2)+"="+sp.latex(expr01))
```

```math
\left[\begin{matrix}a & b\\c & d\end{matrix}\right]^{-1}\left[\begin{matrix}x\\y\end{matrix}\right]=\left[\begin{matrix}- \frac{b y}{a d - b c} + \frac{d x}{a d - b c}\\\frac{a y}{a d - b c} - \frac{c x}{a d - b c}\end{matrix}\right]
```
③「**面倒な計算はコンピューターに丸投げする**」式の考え方では計算が合うかだけが重要。そして実際、上述の式に具体的に同じ数値を代入すると、確かに計算結果が一致する。

* 連立一次方程式の場合

```python
import sympy as sp
a,b,c,d,e,f,x,y = sp.symbols('a,b,c,d,e,f,x,y')
eq10=sp.Eq(a*x+b*y, e)
eq20=sp.Eq(c*x+d*y, f)
#代入
eq11=eq10.subs(([(a, 5), (b, -1),(e,9)]))
eq21=eq20.subs(([(c, 2), (d, 3),(f,7)]))
#解法
expr01 =sp.solve ([eq11, eq21], [x, y]) 
#表示
sp.init_printing()
display(eq11)
display(eq21) 
display(expr01)
print(sp.latex(expr01))
```

```math
5 x - y = 9\\
2 x + 3 y = 7\\
\left\{ x : 2, \  y : 1\right\}
```

* 行列演算の場合

```python
import sympy as sp
a,b,c,d,x,y = sp.symbols('a,b,c,d,x,y')
mt10= sp.Matrix([[a,b],[c,d]])
mt20= sp.Matrix([[x],[y]])
mt11=mt10.subs(([(a, 5), (b, -1),(c, 2), (d, 3)]))
mt21=mt20.subs(([(x,9),(y,7)]))
expr00=mt11**(-1)
expr01 =mt11**(-1)*mt21
#表示
sp.init_printing()
display(mt11)
display(mt21)
display(expr00)
print(sp.latex(mt11)+"^{-1}"+"="+sp.latex(expr00))
display(expr01)
print(sp.latex(mt11)+"^{-1}"+sp.latex(mt21)+"="+sp.latex(expr01))
```

```math
\left[\begin{matrix}5 & -1\\2 & 3\end{matrix}\right]^{-1}=\left[\begin{matrix}\frac{3}{17} & \frac{1}{17}\\- \frac{2}{17} & \frac{5}{17}\end{matrix}\right]\\
\left[\begin{matrix}5 & -1\\2 & 3\end{matrix}\right]^{-1}\left[\begin{matrix}9\\7\end{matrix}\right]=\left[\begin{matrix}2\\1\end{matrix}\right]
```
ところで上掲の状況は以下の様に**基底**の**線形結合**に変形する事も出来て、そしてこの形だと**行列/ベクトル演算における基底のスカラー倍**(この例ではx成分が2倍、y成分が1倍)を説明するのです。

```math
\left[\begin{matrix}5 & -1\\2 & 3\end{matrix}\right]\left[\begin{matrix}x\\y\end{matrix}\right]=\left[\begin{matrix}5 x - y\\2 x + 3 y\end{matrix}\right]=x\left[\begin{matrix}5\\2\end{matrix}\right]+y\left[\begin{matrix}-1\\3\end{matrix}\right]
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/597772/00cc3a03-ae8b-0643-835a-fd7cbb3f8d16.png)
**数学的構造**(Mathematical Structure)の概念って、こんな思わぬ形で人間が普段使いしている多種多様な諸概念を貫いていたりするんですね。
[【用語集】ユージン・ウィグナー『自然科学における数学の理不尽なまでの有効性(1959年)』](https://ochimusha02.hatenadiary.com/entry/2022/02/20/192053)
そんな感じで以下続報…