教えて!gooに以下のような問題が上がってた。
モンテカルロ法（乱数を用いる方法）によって，関数 f(x,y)= -x2 - y2 + 1 について， −1≤x≤1，−1≤y≤1 の範囲で， 乱数を 1,000,000 回繰り返し生成し，
関数 f(x,y) が最大となる (x, y) を 求めるプログラム max_monte.py を作成せよ．
この問題のプログラムと結果を教えてください。問題はこの文章から始まっていました。
# -*- coding: utf-8 -*
"""
max_monte.py プログラム
モンテカルロ法により関数の最大値を探索するプログラム
"""
この問題を見た時、僕が思ったのは、
「こんなバカな問題があるもんか!」
だった。
一体どこの学校がこんなバカな問題を出したのだろう。
一応解説しておく。
モンテカルロ法は乱数を用いたプログラムで、確かに強力な手法である。
ある図形の面積の近似値を調べるには至極便利だ。それは間違いない。
ただし、関数の最大値を求めるには適さないのだ。少なくとも、サイコロをいくら振ったトコロで「これが関数の最大値の近似値を得る独立変数です」とどうやって言うのか。言えるわけがないのである。
なお、関数の最大値、ないしは最小値を求める分野を最適化問題と言う。
実は、これはいまだに厄介な問題で、研究中の分野と言って良い。
と言うのも、我々は高校時代辺りから微分を習い始めて「未知の関数の概形を知る」技術を曲がりなりにも学び始めたが、一体例えば、四次元以上の関数の「概形」を知ったりするにはどうすれば良いのか。複雑な関数の場合は?
僕らは中学校で関数を習い始めたけど、この世には「簡単に図を描いて概形を把握可能な」関数ばっか、とは限らないのである。未知の空間で描くのが難しい関数を「想像する」のはすごく難しいのである。不可能、と言って良い。
そんなワケで、数学的にもコンピュータサイエンス的にも「関数の最大値/最小値を簡単に求めたい」と言う欲求は根強いのである。
しかし、「関数の最大値/最小値」をコンピュータサイエンス的に求める、ってのも遅々として進んでない、と言うのが現状だ。
そんな中、今回は1980年代に登場した「焼きなまし法」を紹介しようと思う。
およそ40年前に見つかった技術ではあるが、科学の長い歴史では殆ど「新参」の技術である。
さて、「焼きなまし法」も確率を使った手法ではある。
つまり、モンテカルロ法の近縁ではあるが、モンテカルロ法ではない。
と言うか、英語では確率を使った計算手法を、Randomized Algorithm(乱択アルゴリズム)と呼んだりする。
つまり、モンテカルロ法も焼きなまし法もRandomized Algorithmの一種だ、と言う事だ。そしてモンテカルロ法は乱択アルゴリズムの代表選手ではあるが、乱択アルゴリズム = モンテカルロ法、と言うわけでもない。
ところで、そもそも、一変数関数相手ならいざ知らず、二変数関数相手だと、単純なモンテカルロ法だと手に負えない範疇に入って来る。
と言うのも、モンテカルロ法は「乱数を使って」砂を撒くような行為をしていくわけだが、当然「線に砂を撒く」のと「面に砂を撒く」のでは後者の方が濃淡が激しくなっていくのは想像に難くないだろう。
つまり、乱数手法、ってのはある一定回数行うと「ある範疇にそこそこ綺麗に砂を撒ける」と言う前提があるからこそ信頼出来るのであって、もし濃淡が濃く現れるようになったら信頼度は減少するのである。これが3次元空間だったら・・・4次元空間だったら・・・と次元数を増やせば増やす程「砂を撒く」行為の信頼度は落ちていく。
これはあくまで比喩だが、貴方が1m四方の範疇で砂を撒け、って言われてやってみるのと1km四方で砂を撒け、って言われてやってみると、後者の方が均一度が悪くなる、ってのは「当たり前だよな」と納得出来るだろう。しかも両者共に「たった500gの砂で」と言われると「そりゃシャレにならんわ」って反応になると思う。
もちろん、後者の場合1tの砂を使っていいよ、って言われると・・・それを持てるか持てないかはさておき、可能性はグンと上がるだろう。
上がるが、コンピュータもそうなんだけど、コンピュータは繰り返しは得意ではあるが、どうせやるなら「なるたけ試行回数を少なく適切な結果を得たい」モンなのである。
それも当たり前だよな、って話である。
そしてこの問題は、「砂を撒く」のが二次元平面上なので、数直線上に「砂を撒く」より遥かに厄介なシチュエーションなのだ。
ホント、一体どこのバカがこんな問題出して「モンテカルロ法を使え」なんぞ言いだしたのか、そのツラを是非とも拝みたいモノである。
基本的な焼きなまし法のアルゴリズムはWikipediaに書いてある通りである。
そこに載ってる疑似コードをバカ正直に実装すれば、一変数関数相手の最小値を求める焼きなまし法はすぐ実装可能だろう。「一変数関数」で「最小値」を求めたい場合は、だ。
与題が厄介なのは
- 二変数関数相手である事
- 「最大値を求めよ」である事
だ。
なお、先にも書いた通り、このテのアルゴリズムは「未知のどんな概形になるか分からない」関数の最大値/最小値を求める為に開発されている。
しかし、さすがに与題の関数はそこまで意地が悪くはない。概形自体はすぐ分かる簡単なモノだ。
平たく言うと、x = 0、y = 0の時点で最大値は1になる。
つまり、答えはその近傍になれば良い。
それを得る為には次のようなコードを書く。
#!/usr/bin/env python3
from random import uniform, randint, normalvariate
from math import exp
# 与題の関数
def func(x, y):
return - x ** 2 - y ** 2 + 1
# 焼きなまし法
def simulated_annealing(func, bounds, n_iterations, step_size, temp):
# 初期値 x, yを一様乱数で適当に選ぶ
best = uniform(*bounds[0]), uniform(*bounds[1])
# 初期値での関数の値を求める
best_eval = func(*best)
# 現在の座標、関数の値を「現在の値」として退避させる
curr, curr_eval = best, best_eval
# 焼きなまし法心臓部
for i in range(n_iterations):
# x, yの近傍値を選ぶが、xとyのどちらを変化させるか、も乱数に依存させる
d = randint(0, 1)
candidate = ((normalvariate(curr[0], step_size)), curr[1]) if d == 1 else (curr[0], (normalvariate(curr[1], step_size)))
# 近傍値での関数の値を求める
candidate_eval = func(*candidate)
# 近傍値での関数の値が現在の関数の値より大きいかどうか調べる
if candidate_eval > best_eval:
# もし近傍値での関数の値が現在の関数の値より大きかったらbestを更新する
best, best_eval = candidate, candidate_eval
# 近傍値の関数の値を現在値の関数の値の差を求める
diff = curr_eval - candidate_eval
# 現時点での温度を計算する
t = temp / float(i + 1)
# 許諾か棄却かの確率を計算する
metropolis = exp(-diff / t)
# 現在値を保持するか更新するかを確率に従い決める
if diff < 0 or uniform(0, 1) < metropolis:
# 現在値を更新
curr, curr_eval = candidate, candidate_eval
return best, best_eval
if __name__ == '__main__':
# 探索範囲を決定
bounds = [[-1, 1], [-1, 1]]
# 繰り返しの最大値を指定
n_iterations = 1000000
# 移動距離を設定
step_size = 0.002
# 初期温度を設定
temp = 400
# 焼きなまし法の結果を出力
print("座標: {0}\n関数の最大値{1}".format(*simulated_annealing(func, bounds, n_iterations, step_size, temp)))
さて、Wikipediaの疑似コードの場合、そもそも現在値の近傍値を探る関数(文中ではNEIGHBOURになっている)は具体的には定義されていない。
その理由は、色んなシチュエーションが考えられて、「これだ」と言う関数が一般的に提案出来ないから、である。
実際、Web検索をしてみると、関数の最大値あるいは最小値を求めるより、巡回セールスマン問題、と言う問題を解く例を見かける方が多いだろう。つまり、「近傍」の定義は問題によって変わるわけだ。
ここではPythonのnormalvariate関数を利用して、性器正規乱数を用いて現在値の近傍値を選出してる。平均値に現在値を与え、標準偏差にステップ数を与える事によって、限りなく現在値に近い範疇を近傍値にしてるのだ。
また、関数のxとyを「同時に」変更するのは頂けない。xを変更するのか、yを変更するのか、と言うのも確率的に選出した方が良い。と言うのもxとyのいずれかが良い解を持ってても、「同時に」変化させるとせっかく求めてた「良い状態」が崩れる可能性が高いのである。結果、最終的な探索には良い結果を生まないだろう。
次は温度計算だ。WikipediaではTEMPERATURE関数として独立して定義されている。
これは具体的には何でも良い模様だ。重要なのは繰り返しが進むにつれて「温度が下がっていく」ような関数を設計すれば良い、と言うこと。それさえ守ればこのコードのような単純な計算で設定しようと何も問題はない。
さて、上のスクリプトを実行すれば、例えばこのような結果が得られるだろう。
座標: (1.4630158816264678e-06, -5.843012968705704e-06)
関数の最大値0.9999999999637188
>>>
関数の最大値は限りなく1に近い値が得られている。座標もそれぞれ0ではないが、かなり0に近くなった「小さい」値が得られてるのが分かるだろう。
と言うわけで実験は性交成功と言える。繰り返しの値、ステップ数、初期温度を色々と弄ってみるのもいいだろう。
なお、2変数関数に於いてはDual Annealing法と言うより優れた「焼きなまし法」の改良手法が提案されている模様である。
が、ネットを検索してみると、まだこなれてないようで、上手い具合にその手法を解説してるサイトは引っかかってこなかった。
そして、Scipyと言うライブラリにその名もズバリ、dual_anneling.pyと言うユーティリティが提供されている模様だが、これは「最小値探索」にだけ限定されてるようで、今回の与題には明らかに適さない。
また、NumpyもScipyもいまだPython3.9以降には対応してないようで、「どのPythonでも使用可能」ではないのが現状なので、今のトコ、一般的な意味では使用を見送らざるを得ないのである。
