タイトル書いてから、そんな用語あるのかな、と思って web記事を漁ってみたら、あったので続行。
よくあるシミュレーションなんかでは、 非常にたくさんのパラメタで初期条件が決まって、 出て来る値は一つの実数、 なんていうことがあります。 単純に考えると引数をたくさん持つ関数という事です。
シミュレーションを行う目的というと、最適な条件を探すためだったりするわけですが、 これら多数の変数の定義域をそれぞれ探っていると 大変なことになってしいます。 だから、最適条件を効率的に探すには、何かうまいやり方を考えねば。
きっと誰かがとっくに考えるはずだと思い、 まえの職場の同僚である藤原さんに訊いてみました。 そこで教わったのが数値演算の黄金時代(1960年代)に発明された、 Nelder-Mead法、通称 downhill simplex method です。 もしこれが、 irb なんかでその場でこしらえた Proc に対して実行できたら 面白い、じゃなくて便利じゃないですか。 Ruby はかなり動的な言語なので、そういう使い方もうまく書けるし、 きっとできるはず。 そこでさっそく Proc のメソッド dhsmplx として実装してみましたよ。 使い方はこんな感じです。
irb(main):2397:0> himmelblau=lambda{|x, y|(x**2 + y -11)**2 + (x + y**2 - 7)**2}
#<Proc:0xb29514e4@(irb):2397>
irb(main):2398:0> simplex
[[0, 0], [1, 0], [0, 1]]
irb(main):2399:0> result= himmelblau.dhsmplx(simplex); ""
""
irb(main):030:0> result[0]
[[0, 0], [1, 0], [0, 1]]
irb(main):031:0> result[-1]
[[2.99999997498766, 1.99999998913814], [2.9999999732276, 1.99999998837382], [2.99999998080073, 1.99999999166252]]
一行目で関数を作ってます。 たかが2変数の初等関数ですけど。 2行めは2次元単体を作っています。これは座標3つからなるarrayです。 これを引数として dhsmplx を実行すると、戻り値は 収束までの simplex の履歴からなる array です。 つまりこの例では [3, 2] で収束しているぽいですね。 その時の値を計算すると
himmelblau.call(*result[-1][0]) 3.05871121083329e-14
およそゼロということでよろしかったでしょうか。
元の関数をみるとゼロが lower bound になるのは自明ですが、 本当にゼロになるような x,y があるかどうかは自明ではないし どこで最小値をとるかとなると、 少なくとも私には皆目見当もつきません。 ノイマンとかラマヌジャンだったら判るかもだけど。 だから、けっこう便利です。
ソースを掲載して説明するまえに、 数理的な背景をちょっと書いておきたいと思います。 そもそも simplex とは何か？ simplex は、自転車ヲタク以外は「単体」といいます。 ちなみに自転車ヲタクはこれを「サンプレックス」と読みます。 単体には次元があります。 1次元単体は線分です。 2次元単体は三角、3次元単体は4面体。 n次元単体は有界な凸多面体であって、その面は n-1次元単体という関係です。 再帰でスマソ。
n次元ベクトル空間はn個の基底の線形和でもれなく調べ尽くせます。 これを幾何学的対象におきかえると、 基底の数はn個ですから、それに原点を追加しまして、 線形独立なn+1個の点になります。 これらn+1個の点が張る convex hull が n次元単体です。
だから、n次元単体を使えば、n次元ベクトル空間を洩れなく 捜索できるわけですが、では、これを使うと要領よく捜索できる、 というのはどのような仕組みによるのでしょうか。
アルゴリズムは単体の各点で問題になっている関数の値を計算し、 一番成績の悪い点をマシなものと置き換える、という操作を繰り返すというものです。 置き換え方は、一番成績の悪い点を、単体の重心を軸として対称移動したものの 成績で場合わけします。
反射点の成績がそれなりなら、最下位と反射点をいれかえたものを 新たな単体として使えばよい。 反射点の成績が一番良い時は、反射方向に正解があるわけで、 この場合も最下位と反射点を入れ換えてもいいけれど、 ここは少し最適化して、 反射方向にもう少し伸ばした点のほうが成績が良い可能性があるわけで、 そっちを使う場合も考慮する。
反射点がかえって成績が悪い場合とは、単体の近所に正解がある場合で、 単体を小さくして行くと事情が判って来る。 まず、最悪点を重心方向に少し詰めた点の成績をみる、 これがそれなりだったら最悪点と詰めた点を入れ換えたものを新たな単体にする。 これもあいかわらず悪い成績だったら、一番成績の良い点の方に、 他の点を全部ちょっと詰めたものを新たな単体として使う。
こうして、単体が足をのばしたりひっこめたりしながら正解に じわじわと近付いていくのがこのアルゴリズムのだいたいのイメージです。
では、以下に dhsmplx のソースを掲載します。
class Proc
def dhsmplx(simplex, hist=[])
hist.push(simplex)
ss=simplex[0]; size= simplex.map{|v|ss.mplus(v*-1).abs}.inject(0){|r, v|r+=v}
if size < 0.00000001
hist
else
ref = simplex.map{|v|
[self.call(*v), v]}.sort
best = ref[0]; worse = ref[-2]; worst = ref[-1]
nsimplex = ref[0..-2].map{|e|e[1]} ## remove worst
cp = simplex.transpose.map{|ax|ax.ave} ## centroid
cr = (cp*2.0).mplus(worst[1]*(-1)) ## reflection of the worst
rv = self.call(*cr)
if best[0] <= rv and rv < worse[0] ## case intermediate
dhsmplx(nsimplex.push(cr), hist)
else
if rv < best[0] ## best
ep = (cp*3.0).mplus(worst[1]*(-2))
if self.call(*ep) < best[0]
nsimplex.push(ep)
else
nsimplex.push(cr)
end
dhsmplx(nsimplex, hist)
else ## case worst
ctrp = (cp.mplus(worst[1]))*0.5 ## contract point
if self.call(*ctrp) < worse[0]
dhsmplx(nsimplex.push(ctrp), hist)
else
nsimplex = ref[1..-1].map{|e|e[1]}.map{|v|(best[1].mplus(v))*0.5}.push(best[1]) ## conctact others to the best point
dhsmplx(nsimplex, hist)
end
end
end
end
end
end
案外長いですな。
あとは、 珍しく適宜コメントが入っているので読めばだいたい判ると思いますが、 若干、標準で装備されていない謎なメソッドが入っているのを説明しておきます。 array のメソッド mplus というのは array をベクトルもしくは行列とみなしてその加法演算をするメソッドです。 かけざん(*)はベクトルもしくは行列の乗算です。 最初の4行が終了条件、あとが本体です。 任意次元で動くように、Ruby の可変長引数機能を使っています。 実装は、新しい単体を作って再帰呼出する構成です。
アルゴリズムの構成から自明ですが、 global solution にたどりつかないおそれがあります。 また、こちらはあまり自明ではありませんが、 初期値の選び方と対象とする関数によっては、変なところに収束します。 numerical recipes をみると、収束先をスタート地点にしてもう一回やってみろ、 みたいな事が書いてあります。
さて、すこしこのアルゴリズムについて考察してみます。 どこが有利かというと、 一回に動かす点は多くの場合、一つ、という事が有利です。 計算プロセスのうち、一番資源を必要とするのは 多くの場合、 self.call の部分ですからこれが少なくて済むのが賢いですね。 この実装ではそうなっていませんが、動かなかった点における 結果は使い廻せるからです。
演習: 次元をdとする。self.call は定数と仮定して、 時間計算量をd の関数として見積もろう。
何変数関数でも動くし、irb で Proc を作ってすぐに計算できますから、 それなりに便利ですよ。
ソースの利用と再配布は Ruby 配布条件に準じます。