【備忘録】pythonで逆行列を使いながらニュートン法を試してみる
最近sympyライブラリにはまっています。とても便利です。
ただ単に如何にsympyが便利か、ということのために練習の記録をアップしました。
今回のお題
ガーデン家具の製造業者が木枠とアルミ枠のガーデンチェアを製造しています。
木枠は1つ18ドルの製造経費、アルミ枠は10ドル、とし、製造台数を
木枠：x
アルミ枠：y
とする。
また、販売価格が
木枠：10+31*x**-0.5 + 1.3*y**-0.2
アルミ枠：5+15*y**-0.4 + 0.8*x**-0.08
となるとわかっているとき、
1日あたりの利益の最大値を考える
また例の本からです*1
先ず、グラフにしてみる
1日の収入zは以下になります。
（＊なんでわざわざこんな複雑なモデルにしたんでしょうか？著者は）
z = x*(10+31*x**-0.5 + 1.3*y**-0.2) -18*x + y*(5+15*y**-0.4 + 0.8*x**-0.08) -10*y
これをsympyを使ってプロットしてみます。
# -*- coding: utf-8 -*-
from sympy import *
z,x,y= symbols("z x y")
z = x*(10+31*x**-0.5 + 1.3*y**-0.2) -18*x + y*(5+15*y**-0.4 + 0.8*x**-0.08) -10*y
from sympy.plotting import plot3d
plot3d(z,(x,0.1,10),(y,0.1,10))
何故か「0」を入れると塁乗ができないよ！と怒られるので、
0.0 cannot be raised to a negative power
プロットの時は、xとyの開始の値を0.1にしておきました。
ランダム探索法を使ってみる
少し狡いですがグラフから少なくともxもyも10以下でしょう、と想定して、10までの中で先ずランダムに初期値を求め ｚの初期値zminを出します。
ループの中で、更にランダムにx,yを出し、その場合のzの値z1を求めます。
そのz1がzminより大きい間はループを続けるようにし、zminの値も新たにz1で更新します。
1000回ほど繰り返してとめた値をzを最大化するx,yの近似値だと考えます。
zmin,xmin,ymin = symbols("zmin ymin xmin")
import random
#seedは使わなくても問題ないですが何となく
random.seed(1)
#10以下だとふんでいるため0～10の間の乱数を作ります
x1 = random.uniform(0,10)
y1 = random.uniform(0,10)
#そして代入してzの初期値を求めます
zmin = z.subs([(x,x1),(y,y1)])
#以下、繰り返し部分です
for i in range(0,1000):
#新たな乱数でzを求める
x1 = random.uniform(0,6)
y1 = random.uniform(0,6)
z1 = z.subs([(x,x1),(y,y1)])
#新しいzの値z1の方が先ほどの値より大きいうちは繰り返し
if z1 > zmin:
xmin = x1
ymin = y1
#zminの値を更新するのを忘れないようにします
zmin = z1
print('xmin',xmin)
print('ymin',ymin)
print('zmin',zmin)
結果
('xmin', 4.657253897915613)
('ymin', 5.9093752268644355)
('zmin', 52.0719318635198)
木枠は4.6台、アルミ枠は5.9台（およそ5台と6台）で、利益は52ドル。
ということになりそうです
ニュートン法を使ってみる
ランダム探索法では、本当に当てずっぽうなので繰り返し回数を多くしないとなかなか近い値にはなりません。もう少し早く収束するべく、ニュートン方を使います。
先ほどは直接zの値の大小を比べてループさせていましたが、今回ニュートン法では、zの勾配（微分）を0とするようなx,yを求めます。
（＊先ほどのグラフで、zの最大値は極値を求めればよい、とわかっているからです）
xn,yn = symbols("yn xn")
xn = random.uniform(0,10)
yn = random.uniform(0,10)
for i in range(0,15):
N = Matrix([[xn],[yn]])
A = Matrix([[dzdx.diff(x),dzdx.diff(y)],
[dzdy.diff(x),dzdy.diff(y)]])
A = A.subs([(x,xn),(y,yn)])
F = Matrix([[dzdx.subs([(x,xn),(y,yn)])],
[dzdy.subs([(x,xn),(y,yn)])]])
#ニュートン法のキモの部分
#手作業で逆行列の要素を計算式にしなくてもsympyは大丈夫です
N = N - A.inv()*F
#xnとynの値を更新します
xn = N[0]
yn = N[1]
#確認で都度、出力してみます
print("N",i,N)
print("z",i,z.subs([(x,xn),(y,yn)]))#最終の値を出力
print("N",N)
print("z",z.subs([(x,N[0]),(y,N[1])]))
先ほど同様に初期値はランダムに求めて、行列Nとしています。
一般的にニュートン法では、偏導関数を成分とする行列をAとすると
F(x(0)) + A(x - x(0) )=0
を解いて
x(1) = x(0) - (Aの逆行列)*F(x(0))
と更新するのでした。
このとき Aの逆行列を求めるのが手計算ではしんどいのですが、sympyでは、A.inv()と記述するだけで済みます。ちょっと感動的ではないでしょうか？
試しに例をここに
from sympy import *
a,b,c,d = symbols("a b c d")
A = Matrix([[a,b],[c,d]])
print(A)
print(simplify(A.inv()))
結果は以下の通り、見慣れた公式の通りです。
Matrix([[a, b], [c, d]])
Matrix([
[ d/(a*d - b*c), -b/(a*d - b*c)],
[-c/(a*d - b*c), a/(a*d - b*c)]])
また、先ほどのコード結果は、
('N', 0, Matrix([[3.11021951910684],[4.39753563309424]]))
('z', 0, 50.4978633181626)
('N', 1, Matrix([[4.26447040602873],[5.59713059033704]]))
('z', 1, 51.9880863435632)
('N', 2, Matrix([[4.66015416182119],[5.84578821398569]]))
('z', 2, 52.0723541378789)
('N', 3, Matrix([[4.68944643469009],[5.85200277371353]]))
('z', 3, 52.0726917916350)
('N', 4, Matrix([[4.68958529813255],[5.85199380485953]]))
('z', 4, 52.0726917987936)
('N', 5, Matrix([[4.68958530121007],[ 5.8519938044504]]))
('z', 5, 52.0726917987936)
('N', 6, Matrix([[4.68958530121007],[ 5.8519938044504]]))
('z', 6, 52.0726917987936)
('N', 7, Matrix([[4.68958530121007],[ 5.8519938044504]]))
('z', 7, 52.0726917987936)
('N', 8, Matrix([[4.68958530121007],[ 5.8519938044504]]))
('z', 8, 52.0726917987936)
('N', 9, Matrix([[4.68958530121007],[ 5.8519938044504]]))
('z', 9, 52.0726917987936)
('N', 10, Matrix([[4.68958530121007],[ 5.8519938044504]]))
('z', 10, 52.0726917987936)
('N', 11, Matrix([[4.68958530121007],[ 5.8519938044504]]))
('z', 11, 52.0726917987936)
('N', 12, Matrix([[4.68958530121007],[ 5.8519938044504]]))
('z', 12, 52.0726917987936)
('N', 13, Matrix([[4.68958530121007],[ 5.8519938044504]]))
('z', 13, 52.0726917987936)
('N', 14, Matrix([[4.68958530121007],[ 5.8519938044504]]))
('z', 14, 52.0726917987936)
('N', Matrix([[4.68958530121007],[ 5.8519938044504]]))
('z', 52.0726917987936)
となります。
15回も繰り返しましたが、既に5回めで収束していますね。値もランダム探索法の結果とほぼ同じです。
---------------------------
sympyの逆行列の求め方について
ニュートン法について
