scipy.optimizeのminimizeを使って多変数の最適化を行う例を紹介します。
曲線内の面積最大多角形を求める
【例題１】$y=x^2$と$y=1$に囲まれた領域に入る５角形の面積の最大値を求めよ
３点(-1,1),(0,0),(1,1)は固定し、残りの対称な２点の右側を$(x,y)$とすると５角形の面積$A(x)$は図の三角形の面積の和となるので、
\begin{align}
A(x)&= 2 \times (\frac{x}{2} + \frac{1-y}{2}) \\
&=-x^2+x+1 \\
最大値はA'(x) &= -2x+1 = 0 \ を解いて\\
x &= 0.5のとき
\end{align}
となります。これをscipy.optimizeのminimizeで同様に求めます。
scipy.optimizeのminimizeを使う
scipy.optimizeに関してはこのYoutube科学技術計算モジュール SciPy 解説 (1) 最適化が分かりやすいと思います。
基本的には以下の２つをminimizeに渡すと最小値とその時の変数の値とを返してくれます。
- 変数xの初期値の配列
- 変数xから評価する値を返す関数
1. 変数xの初期値の配列
今回の場合には２点(-1,1),(1,1)を除いた３点のx座標をx1,x2,x3の初期値$x0=[-0.5,0,0.5]$を配列で渡します。
一般のN角形では以下のコードで区間をN等分して両端を除いた配列を作ります
import numpy as np
N = 5
x0 = np.linspace(-1, 1, N)[1:-1]
print(x0)
# [-0.5 0. 0.5]
2. 変数xから評価する値を返す関数
ゴールは多角形の面積の最大値を求めることですが、minimizeは最小値を求めるので評価する関数は、多角形の下の台形の面積の総和を求めます。
\begin{align}
f(x)&= \frac{1}{2}\sum_{i=0}^{N-2}(x_{i+1}-x_{i})\times(y_{i+1}+y_{i}) \\
\end{align}
これをコードにすると、
def g(x): return x**2
def f(x):
x = np.concatenate(([-1], x, [1]))
y = g(x)
return np.dot(x[1:] - x[:-1], y[1:] + y[:-1])/2
scipyのminimizeを呼ぶ
このx0とfをminimizeに渡すだけで、答えが求まります。fの最小値はresult.funに返ってくるので2から引けば多角形の最大面積が得られます。
xの答えは初期値と変わらないのであまり面白くないですが、コードは関数やNを変えてもそのまま使えます。
result = scipy.optimize.minimize(f, x0)
print(result.x)
# [-0.5 0. 0.5]
print(2-result.fun)
# 1.25
関数をy=x^4にする
【例題１】$y=x^4$と$y=1$に囲まれた領域に入る５角形の面積の最大値を求めよ
まず同様にA(x)の最大値を求めると
\begin{align}
A(x)&= 2 \times (\frac{x}{2} + \frac{1-y}{2}) \\
&=-x^４+x+1 \\
A'(x) &= -4x^3+1 = 0 \ となるのは\\
x &= 0.63のとき
\end{align}
scipyのminimizeでも$g(x)$を以下に変えると同じ答えで 面積の最大値1.47247 が得られました。
def g(x): return x**4
:
print("["+", ".join(f"{x:.02f}" for x in result.x)+"]")
# [-0.63, -0.00, 0.63]
print(f"{2 - result.fun:.5f}")
# 1.47247
めでたしめでたしと思いましたが 実は驚くべき事実が待ってました。
1.47247は極大値だが最大値ではなかった！
scipyのminimizeを使うときに気をつけなければならないのが$x$の初期値です。この場合は答えの５角形が左右対称と決めつけて初期値に$0$を含めてしまいました。
そこでrandomのuniformを使って[-1,1]のxをランダムに選んだ値を初期値にしてみます。
x0 = [uniform(-1, 1) for _ in range(N-2)]
print(x0)
# [-0.37425418273992905, 0.09640569845950586, 0.47695576985059995]
:
print("["+", ".join(f"{x:.02f}" for x in result.x)+"]")
# [-0.58, 0.30, 0.71]
print(f"{2 - result.fun:.5f}")
# 1.47731
すると以下のような図形が答えとなり 前より大きい面積の最大値1.47731 が得られました。
もちろん初期値を乱数で選んでいるので必ず正しい答えが得られる保証はありませんので何度か実行する必要があります。
最大面積のN角形を求める
プログラムのNを変えれば良いだけなので$y=x^4$, $N=7$の場合を求めます。
N = 7
x0 = [uniform(-1, 1) for _ in range(N-2)]
print(x0)
# [0.3563244927134823, -0.20585104759327488, 0.24763729825802683, -0.22909622813792185, -0.10322665332406222]
:
print("["+", ".join(f"{x:.02f}" for x in result.x)+"]")
# [-0.75, -0.42, 0.28, 0.58, 0.81]
print(f"{2 - result.fun:.5f}")
# 1.54878
（開発環境：Google Colab）
この考え方はProject Euler, Problem 897: Maximal n-gon in a regionを解くのに役に立ちます