データの転がる音がする

気が向いたら取り組んだことをまとめて記事にしたいと思います。機械学習やソフトウェアの話、はたまた読んだ本についてなど、内容は気が向くままにです。

関数の最適化をしたり可視化をしたり

今日は凝縮物理の期末レポートを提出してきました!これで無事卒業できそうです。笑

ベイズ統計や機械学習を扱っていると、尤度の最大化など、関数の大域的な最適化にしばしば出会います。きっと実装には手間がかかるんだろうなと思っていたのですが、Scipyに用意されているoptimizeクラスを用いたら驚くほど簡単でした。

大域的最適化の手法

そもそも大域的最適化ってどうやるのか?大域的最適化のアルゴリズムにはいくつか種類があります。

(1) Basin hopping法

    初期値->局所最適を繰り返す

(2) Brute法

     格子状に力任せで探索 

(3)  Differential_evolution法

     進化的アルゴリズムを使い、上手に探索

今回は最適化にはscipy.optimizeに用意された、differential_evolution関数を用います。実際のコードでは、differential_evolution( 最適化したい関数名 , 変数の定義域 ) と打ち込むだけで結果が帰ってきます。非常に簡単です。

関数の例:自発的対称性を破るワイン型ポテンシャルV(x)

今日は期末レポートでひたすら波動関数微分積分していたので、懐かしの物理の分野から。宇宙論素粒子論には「自発的対称性の破れ」という概念があり、この議論から神の粒子と呼ばれるHiggs粒子の存在が導かれます。この理論で出てくる有名なワインボトル型のポテンシャル関数は、局所最小な値がいくつもある形をしているので、今回は可視化・局所最小探索にこの関数を使ってみます。ただ、実際には局所最小値は1つになって欲しいので、このポテンシャルに、傾きの項を加えてあげます。

Ex1: 1次元ポテンシャル

下のコードを実行するとポテンシャルが可視化されます。左が良く見るワイン型で、右が摂動項を加えて傾かせたポテンシャル。コード後半では、傾きを加えたポテンシャルの最小値探索を実行しています。変数の定義域 = bounds を定めて、それを differential_evolution関数に代入すれば結果が返ってくるので本当に簡単。図の赤い点が大域的最小値です。

# 1次元ポテンシャルの可視化
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import rosen, differential_evolution, minimize, basinhopping
from mpl_toolkits.mplot3d import Axes3D

# パラメータ
mu = 8
lamda = 1

# ポテンシャルの定義
def V1(x):
    return -(mu ** 2)* (x**2) + lamda*(x **4)

def V1pert(x):
    return -(mu ** 2)* (x**2) + lamda*(x **4) - 100*x

# ポテンシャルの可視化
x = np.arange(-10, 10, 0.1)
y = V1(x)
plt.plot(x,y)
plt.show()

y = V1pert(x)
plt.plot(x, y)
plt.show()

# 最適化の実行
bounds = [(-10.0, 10.0)]
res = differential_evolution(V1pert, bounds)
print res

# 最適化の結果の可視化
x = np.arange(-10, 10, 0.1)
y = V1pert(x)
plt.plot(x, y)
plt.scatter(res.x, V1pert(res.x), c="r")
plt.show()

f:id:idkphys-030zzz:20160813021630p:plain f:id:idkphys-030zzz:20160813021829p:plain

Ex2: 2次元ポテンシャル

次に2次元ポテンシャルを解析してみます。こっちの方が見た目に楽しいかも。やっていることはEx1と一緒です。 グリッド上の3次元プロットが新しかったので、その実装の仕方は記憶に留めておこうと思います。

# 2次元ポテンシャルの可視化と最適化
#パラメータ
mu = 12
lamda = 1
alpha = 100

# ポテンシャルの定義
def V2(x):
    return  -(mu ** 2)* (x[0]**2+x[1]**2) + lamda*(x[0]**2+x[1]**2)**2

def V2pert(x):
    return  -(mu ** 2)* (x[0]**2+x[1]**2) + lamda*(x[0]**2+x[1]**2)**2 +alpha*(x[0]+x[1])

# 3次元プロット
x = np.arange(-10.0, 10.0, 1)
y = np.arange(-10.0, 10.0, 1)

xx, yy = np.meshgrid(x, y, sparse=True)
z = -(mu ** 2)* (xx**2+yy**2) + lamda*(xx**2+yy**2)**2

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_wireframe(xx, yy, z)

xx, yy = np.meshgrid(x, y, sparse=True)
z = -(mu ** 2)* (xx**2+yy**2) + lamda*(xx**2+yy**2)**2 + alpha*(xx+yy)

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_wireframe(xx, yy, z)

f:id:idkphys-030zzz:20160813021633p:plain f:id:idkphys-030zzz:20160813023218p:plain

疲れたので今日はここまで!

広告を非表示にする