前回に引き続き、MicroPythonにおける「フルPythonにおける numpy 的な」演算ライブラリ ulab の機能を試用して行きたいと思います。今回のお題は「連立方程式」であります。別投稿にて「解答済」の問題があるので、MicroPythonでも解いてみようという趣向。しかし、今回は実装の違いによる問題に躓いてしまいました。トホホ。
※「MicroPython的午睡」投稿順 Indexはこちら
（今回使用のコード全文は末尾にあります）
ご存知のとおり「フルPython」の numpy は巨大なライブラリです。それに対応するMicroPythonの ulab もそれなりに大きなライブラリではあるのです。前回見た通り、ほとんど「同じ」コードで「ほぼ同等の」結果を得られることが多いと思われます。しかし、ケースによっては以下の２つの差異があからさまに見えることがあります。今回はまさにそんなケースでした。
- Python の numpy と MicroPythonの ulab の差異
- ulabの機種毎実装の差異
1は numpy と「似せて」作った ulab といえ、組み込み向けに「軽く」してあるための差異と言えます（まったく同じにするならば numpy そのものを載せるしかないし。）2は、標準的な MicroPython 実装(多分 pyboard向け）に対して、MicroPythonの機種毎移植に伴う実装のバラツキと思います。マイコンボード毎にハード構成などは大きく異なるし、移植者の考え方にも違いがある、と。本稿でターゲットといたしますのは、MicroPythonの実装の一つである MaixPy のさらにM5StickV向けの実装ということになります。
解こうとする課題—連立方程式
今回解こうとするのは、以下の投稿で取り上げさせていただいた回路の問題です。オリジナルはアナデバ社(ADI社）のStudentZoneの課題です。以下の別件投稿では解法の一つが連立方程式に帰着しました（アナデバ社の模範解答とは異なるアプローチですが。）それをそのままなぞって同じ結果が得られればOKという目論見です。
お手軽ツールで今更学ぶアナログ(31) スマホアプリで回路の宿題
Python 3.8.0上の numpy で連立方程式を解く方法
行列を中心とした数値演算ライブラリである numpy での連立方程式の解き方を簡単におさらいしておきます。２元連立方程式、中学生（小学生？）レベル。
x + y = 5 3*x + 4*y = 11
上の問題は、係数行列の逆行列を求めて右辺を掛けるという形で、Python上の numpy で解くことができます。このとき用いるのが numpy の「行列」表現である matrix です。なお、np.linalg.inv()というのが逆行列を求めるメソッドです。
import numpy as np
AM = np.matrix([
[1, 2],
[3, 4]
])
YM = np.matrix([
[5],
[11]
])
XM = np.linalg.inv(AM)*YM
print("Matrix:\n",XM)
しかし、これをそのまま MicroPython へもっていこうとするとまず引っかかるのが、
MicroPythonのulabには matrix が見当たらない
という点です。MicroPythonは、numpyのarray(ndarray、配列)と互換性のあるarrayで処理するのが基本のようです。行列と配列、似たようなもの？といいつつ違いがあります。そこでフルPythonのnumpyでも、以下のような書き方であれば MicroPythonの ulab と互換性をとることが可能だと思います。
import numpy as np
AA = np.array([
[1, 2],
[3, 4]
])
YA = np.array([
[5],
[11]
])
XA = np.dot(np.linalg.inv(AA), YA)
print("Array, dot : \n", XA)
ここで注意しなければならないのは、matrixであれば「逆行列 * Y」という*演算子で書けるのが arrayであると「np.dot(逆行列, Y)」という形で dot()を使わねばならないことです。arrayの場合 *演算子は要素毎掛け算のスカラー積になってしまいます。私はまずこれ（行列と配列の挙動の違い）につまづきました（お約束的。）
連立方程式は「解ける」のか？
「スマホアプリで回路の宿題」の投稿では、回路の結節点毎に電流で、小ループ毎に電圧で方程式を作って解きました。未知数にとった各ノードの電流よりも方程式の数の方が１個多い overdetermind な状態であります。しかし投稿では、手計算を諦め、連立方程式を解くのに「あろうことか」 Maxima 使わせていただきました（スマホへも移植されているのです。）Maxima走れば、過剰だろうが何だろうがお任せで結論出してくれます。１個式が多いという指摘のもと、ちゃんと解答が得られました。
しかし、numpy 使って逆行列を求める方法では、７個の未知数であれば７個の独立な方程式が無いと求まりませぬ。方程式は８個。それでエイヤー（何がエイヤか？）と最後の１個を取り除いて numpy で逆行列を求めると
numpy.linalg.LinAlgError: Singular matrix
Singular=非正則行列だと怒られました（singularを辞書引きました。）正則でなければ逆行列が求まりません。しかし、上の投稿で既に「答え」はMaximaから頂いているので解けないわけがない。そこで頭をまったく使わず総当たりです。末尾に掲げましたる「Python 3.8.0で動作確認済み、連立方程式を解くコード」を使い、方程式の７個セット毎に「行列式」をまず求めてみました。昔やりましたでせう。「行列式 が非0」ならば逆行列が求まるのであります。結果はこんな感じ
det(A0)= -20.99999999999999 det(A1)= -20.99999999999999 det(A2)= 21.0 det(A3)= -20.99999999999999 det(A4)= 20.99999999999999 det(A5)= 0.0 det(A6)= 0.0
つまり、最後とその前の方程式を抜いたセットでは求まらないけれど、それ以外の式を１個抜く分には求まるのでした。とりあえずキリのいい行列式の値が求まっているA2のセット（３つ目の方程式を抜いた）を使ってお答えを求めることにいたしました。結果はこんな感じ
Result : [[0.57142857] [0.42857143] [0.14285714] [0.28571429] [0.14285714] [0.42857143] [0.57142857]]
最初の 0.57142857（真の値は7分の4です。 Maximaによると）こそ知りたかった電流値mAであります。一応、Python3.8.0上では numpy 使ってお答えが求まることが分かりました。
MicroPythonでの連立方程式の求解
さてPython上で動作したコードをMicroPythonに持っていきます。ulabの使い方をまとめてくれているThe ulab book を参照したところでは、import numpy を import ulab と書き換えるだけで動きそうに思われました。しかし動作しません。調べると
- ulabの標準的な実装では ulab.linalg という階層があり、行列式det()も逆行列inv()もその中で定義されている
- M5StickV上のulab実装では ulab.linalg という階層はない。しかし、ulabの直下にdet()とinv()がある
そこに実装の大きな違いがあるとも知らず、「名前の違いね、簡単」と浅墓な理解で走らせてみると
input matrix is singular
とまた singular の指摘です。先ほど正則行列になる係数行列を選んで与えているので、singular の筈がありません。Python 3.8.0 上では答えが求まっているのだ！MicroPython上での動作を調べると行列式の値、0が返っています。おかしい。そこで、２ｘ２の行列、３ｘ３の行列、４ｘ４の行列の行列式を求めるプログラムを書いて Python 3.8.0 と M5StickV用上のMicroPythonで動作を比べてみました。なお、行列に詰めた係数の値は、以下の御本より引用させていただきました。（行列式の値を、御本の解答で確認するため）
末尾の「Python 3.8.0で動作確認済み、行列式を求めるコード」を走らせた時の結果がこちら
det(A)= -2.0000000000000004 det(B)= -56.00000000000002 det(C)= 27.999999999999986
小数点以下の奥深くを四捨五入するならば、石村先生の教科書の解答と一致しております。次に「M5StickV用のMicropython上で実行した行列式を求めるコード」を走らせた時の結果がこちら
det(A)= -2.0 det(B)= -56.0 det(C)= 0.0
およよ、３ｘ３行列までは求まっているけれど、４ｘ４行列は０だ。解こうとしている連立方程式は７ｘ７だから４ｘ４でコケているのじゃダメじゃん。３ｘ３行列ならば有名なサラスの公式があるので私でも計算できる（多分途中で計算間違いするに違いありませんが。）ちゃんと実装されてないんじゃないの？
この件について、MicroPython用の「フル」ulabのソースコードを見に行ったところ、その実装のままであれば大丈夫に見えます（動作確かめてないけど。）推測するに、
- ulab.linalg階層があれば多分「ちゃんと」実装されている
- M5StickVのMicroPython v0.5.0-98-g7ec09ea22-dirty on 2020-07-21実装上のdet()やinv()は「なんちゃって」実装（サラス公式レベル？）
であると。まあ、どうしても必要だったら自力で掃き出し法でも実装するか、ulabの本来コードを実装するしかないでしょう。多分、当面、そんな必要も無いでしょが。
MicroPython的午睡(8) M5StickV、ulab行列式、timeitデコレータ へ戻る
MicroPython的午睡(10) _thread、マルチスレッドは出来るけれども… へ進む
Python 3.8.0で動作確認済み、連立方程式を解くコード
import numpy as np
A = np.array([
#A, B, C, D, E, F, G
[1, 1, 0, 0, 0, 0, 0], #A+B=1
[0, 0, 0, 0, 0, 1, 1], #F+G=1
[1, 0, 1, 0, -1, 0, -1], #C+A=G+E
[0, 1, -1, -1, 0, 0, 0], #B=C+D
[0, 0, 0, 1, 1, -1, 0], #E+D=F
[1, -1, -1, 0, 0, 0, 0], #A-C-B=0
[0, 0, 1, -1, 1, 0, 0], #E-D+C=0
[0, 0, 0, 0, -1, -1, 1] #G-F-E=0
])
A0 = np.array([
#A, B, C, D, E, F, G
[0, 0, 0, 0, 0, 1, 1], #F+G=1
[1, 0, 1, 0, -1, 0, -1], #C+A=G+E
[0, 1, -1, -1, 0, 0, 0], #B=C+D
[0, 0, 0, 1, 1, -1, 0], #E+D=F
[1, -1, -1, 0, 0, 0, 0], #A-C-B=0
[0, 0, 1, -1, 1, 0, 0], #E-D+C=0
[0, 0, 0, 0, -1, -1, 1] #G-F-E=0
])
print("det(A0)=", np.linalg.det(A0))
A1 = np.array([
#A, B, C, D, E, F, G
[1, 1, 0, 0, 0, 0, 0], #A+B=1
[1, 0, 1, 0, -1, 0, -1], #C+A=G+E
[0, 1, -1, -1, 0, 0, 0], #B=C+D
[0, 0, 0, 1, 1, -1, 0], #E+D=F
[1, -1, -1, 0, 0, 0, 0], #A-C-B=0
[0, 0, 1, -1, 1, 0, 0], #E-D+C=0
[0, 0, 0, 0, -1, -1, 1] #G-F-E=0
])
print("det(A1)=", np.linalg.det(A1))
A2 = np.array([
#A, B, C, D, E, F, G
[1, 1, 0, 0, 0, 0, 0], #A+B=1
[0, 0, 0, 0, 0, 1, 1], #F+G=1
[0, 1, -1, -1, 0, 0, 0], #B=C+D
[0, 0, 0, 1, 1, -1, 0], #E+D=F
[1, -1, -1, 0, 0, 0, 0], #A-C-B=0
[0, 0, 1, -1, 1, 0, 0], #E-D+C=0
[0, 0, 0, 0, -1, -1, 1] #G-F-E=0
])
print("det(A2)=", np.linalg.det(A2))
A3 = np.array([
#A, B, C, D, E, F, G
[1, 1, 0, 0, 0, 0, 0], #A+B=1
[0, 0, 0, 0, 0, 1, 1], #F+G=1
[1, 0, 1, 0, -1, 0, -1], #C+A=G+E
[0, 0, 0, 1, 1, -1, 0], #E+D=F
[1, -1, -1, 0, 0, 0, 0], #A-C-B=0
[0, 0, 1, -1, 1, 0, 0], #E-D+C=0
[0, 0, 0, 0, -1, -1, 1] #G-F-E=0
])
print("det(A3)=", np.linalg.det(A3))
A4 = np.array([
#A, B, C, D, E, F, G
[1, 1, 0, 0, 0, 0, 0], #A+B=1
[0, 0, 0, 0, 0, 1, 1], #F+G=1
[1, 0, 1, 0, -1, 0, -1], #C+A=G+E
[0, 1, -1, -1, 0, 0, 0], #B=C+D
[1, -1, -1, 0, 0, 0, 0], #A-C-B=0
[0, 0, 1, -1, 1, 0, 0], #E-D+C=0
[0, 0, 0, 0, -1, -1, 1] #G-F-E=0
])
print("det(A4)=", np.linalg.det(A4))
A5 = np.array([
#A, B, C, D, E, F, G
[1, 1, 0, 0, 0, 0, 0], #A+B=1
[0, 0, 0, 0, 0, 1, 1], #F+G=1
[1, 0, 1, 0, -1, 0, -1], #C+A=G+E
[0, 1, -1, -1, 0, 0, 0], #B=C+D
[0, 0, 0, 1, 1, -1, 0], #E+D=F
[0, 0, 1, -1, 1, 0, 0], #E-D+C=0
[0, 0, 0, 0, -1, -1, 1] #G-F-E=0
])
print("det(A5)=", np.linalg.det(A5))
A6 = np.array([
#A, B, C, D, E, F, G
[1, 1, 0, 0, 0, 0, 0], #A+B=1
[0, 0, 0, 0, 0, 1, 1], #F+G=1
[1, 0, 1, 0, -1, 0, -1], #C+A=G+E
[0, 1, -1, -1, 0, 0, 0], #B=C+D
[0, 0, 0, 1, 1, -1, 0], #E+D=F
[1, -1, -1, 0, 0, 0, 0], #A-C-B=0
[0, 0, 0, 0, -1, -1, 1] #G-F-E=0
])
print("det(A6)=", np.linalg.det(A6))
A7 = np.array([
#A, B, C, D, E, F, G
[1, 1, 0, 0, 0, 0, 0], #A+B=1
[0, 0, 0, 0, 0, 1, 1], #F+G=1
[1, 0, 1, 0, -1, 0, -1], #C+A=G+E
[0, 1, -1, -1, 0, 0, 0], #B=C+D
[0, 0, 0, 1, 1, -1, 0], #E+D=F
[1, -1, -1, 0, 0, 0, 0], #A-C-B=0
[0, 0, 1, -1, 1, 0, 0] #E-D+C=0
])
Y = np.array([
[1],
[1],
[0],
[0],
[0],
[0],
[0]
])
X = np.dot(np.linalg.inv(A2), Y)
print("Result : \n", X)
Python 3.8.0で動作確認済み、行列式を求めるコード
import numpy as np
A = np.array([[1, 2],[3, 4]])
print("det(A)=", np.linalg.det(A))
B = np.array([[1, -5, 6],[-3, 0, 1],[2, 4, -2]])
print("det(B)=", np.linalg.det(B))
C = np.array([[0, 3, 2, -1],[-4, 6, 1, 5],[0, -2, 3, -1],[-1,0,0,2]])
print("det(C)=", np.linalg.det(C))
M5StickV用のMicropython上で実行した行列式を求めるコード
import ulab as np
A = np.array([[1, 2],[3, 4]])
print("det(A)=", np.det(A))
B = np.array([[1, -5, 6],[-3, 0, 1],[2, 4, -2]])
print("det(B)=", np.det(B))
C = np.array([[0, 3, 2, -1],[-4, 6, 1, 5],[0, -2, 3, -1],[-1,0,0,2]])
print("det(C)=", np.det(C))