MicroPython的午睡(9) ulabで連立方程式を解く、機種固有実装の蹉跌

Joseph Halfmoon

前回に引き続き、MicroPythonにおける「フルPythonにおける numpy 的な」演算ライブラリ ulab の機能を試用して行きたいと思います。今回のお題は「連立方程式」であります。別投稿にて「解答済」の問題があるので、MicroPythonでも解いてみようという趣向。しかし、今回は実装の違いによる問題に躓いてしまいました。トホホ。

※「MicroPython的午睡」投稿順 Indexはこちら

(今回使用のコード全文は末尾にあります)

ご存知のとおり「フルPython」の numpy は巨大なライブラリです。それに対応するMicroPythonの ulab もそれなりに大きなライブラリではあるのです。前回見た通り、ほとんど「同じ」コードで「ほぼ同等の」結果を得られることが多いと思われます。しかし、ケースによっては以下の2つの差異があからさまに見えることがあります。今回はまさにそんなケースでした。

  1. Python の numpy と MicroPythonの ulab の差異
  2. ulabの機種毎実装の差異

1は numpy と「似せて」作った ulab といえ、組み込み向けに「軽く」してあるための差異と言えます(まったく同じにするならば numpy そのものを載せるしかないし。)2は、標準的な MicroPython 実装(多分 pyboard向け)に対して、MicroPythonの機種毎移植に伴う実装のバラツキと思います。マイコンボード毎にハード構成などは大きく異なるし、移植者の考え方にも違いがある、と。本稿でターゲットといたしますのは、MicroPythonの実装の一つである MaixPy のさらにM5StickV向けの実装ということになります。

解こうとする課題—連立方程式

今回解こうとするのは、以下の投稿で取り上げさせていただいた回路の問題です。オリジナルはアナデバ社(ADI社)のStudentZoneの課題です。以下の別件投稿では解法の一つが連立方程式に帰着しました(アナデバ社の模範解答とは異なるアプローチですが。)それをそのままなぞって同じ結果が得られればOKという目論見です。

お手軽ツールで今更学ぶアナログ(31) スマホアプリで回路の宿題

Python 3.8.0上の numpy で連立方程式を解く方法

行列を中心とした数値演算ライブラリである numpy での連立方程式の解き方を簡単におさらいしておきます。2元連立方程式、中学生(小学生?)レベル。

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の場合 *演算子は要素毎掛け算のスカラー積になってしまいます。私はまずこれ(行列と配列の挙動の違い)につまづきました(お約束的。)

連立方程式は「解ける」のか?

「スマホアプリで回路の宿題」の投稿では、回路の結節点毎に電流で、小ループ毎に電圧で方程式を作って解きました。未知数にとった各ノードの電流よりも方程式の数の方が1個多い overdetermind な状態であります。しかし投稿では、手計算を諦め、連立方程式を解くのに「あろうことか」 Maxima 使わせていただきました(スマホへも移植されているのです。)Maxima走れば、過剰だろうが何だろうがお任せで結論出してくれます。1個式が多いという指摘のもと、ちゃんと解答が得られました。

しかし、numpy 使って逆行列を求める方法では、7個の未知数であれば7個の独立な方程式が無いと求まりませぬ。方程式は8個。それでエイヤー(何がエイヤか?)と最後の1個を取り除いて numpy で逆行列を求めると

numpy.linalg.LinAlgError: Singular matrix

Singular=非正則行列だと怒られました(singularを辞書引きました。)正則でなければ逆行列が求まりません。しかし、上の投稿で既に「答え」はMaximaから頂いているので解けないわけがない。そこで頭をまったく使わず総当たりです。末尾に掲げましたる「Python 3.8.0で動作確認済み、連立方程式を解くコード」を使い、方程式の7個セット毎に「行列式」をまず求めてみました。昔やりましたでせう。「行列式 が非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

つまり、最後とその前の方程式を抜いたセットでは求まらないけれど、それ以外の式を1個抜く分には求まるのでした。とりあえずキリのいい行列式の値が求まっているA2のセット(3つ目の方程式を抜いた)を使ってお答えを求めることにいたしました。結果はこんな感じ

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が返っています。おかしい。そこで、2x2の行列、3x3の行列、4x4の行列の行列式を求めるプログラムを書いて Python 3.8.0 と M5StickV用上のMicroPythonで動作を比べてみました。なお、行列に詰めた係数の値は、以下の御本より引用させていただきました。(行列式の値を、御本の解答で確認するため)

石村園子先生著 「すぐわかる線形代数」東京図書 1994

末尾の「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

およよ、3x3行列までは求まっているけれど、4x4行列は0だ。解こうとしている連立方程式は7x7だから4x4でコケているのじゃダメじゃん。3x3行列ならば有名なサラスの公式があるので私でも計算できる(多分途中で計算間違いするに違いありませんが。)ちゃんと実装されてないんじゃないの?

この件について、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))