九州大学 海洋システム工学専攻
システム設計特論 （木村）6) 数値最適化（２）ランダムサーチ
ダシミュレーテッドアニーリング（ＳＡ）
滑降シンプレックス法
プ
授業の資料等はhttp://sysplan.nams.kyushu-u.ac.jp/gen/index.html
p y p
yjp g1

ベンチマーク問題： Rastrigin関数
コスト関数
スト関数

n



f ( x1 , x2 ,  xn )  10n   xi  10 cos( 2 xi )
i 1

ただし2

 5.12  xi  5.12

n=２の場合のn２の場合のコスト関数の景観・局所解が多数存在する多峰性関数
・探索域中心（全て０）に最適解が存在し、その周辺に格子状に複数の局所解を持つ。

ベンチマーク問題： Rosenbrock関数
コスト関数
スト関数

n



f ( x1 , x2 ,  xn )   100 ( x1  xi )  ( xi  1)i2

ただし2 22 2.048  xi  2.048

n=２の場合のコスト関数の景観
極小点は１つしか無いが極小点は１つしか無いが、そこへたどり着くには長くて曲がりくねった谷を辿らねばならない=変数間の依存関係が強いベンチマーク問題： Schwefel関数
コスト関数
スト関数

nf ( x1 , x2 ,  xn )    xi sini 1 xi

ただし 512  xi  512

n=２の場合のコスト関数の景観・局所解が多数存在する多峰性関数・最適解が探索領域の境界付近に存在

4

勾配法は多峰性関数には使えない！

勾配法＝勾配を下って極小値を探す方法

多峰性関数に勾配法は使えない！

【解決策】

１．ランダムに動く
２．多点で探す5

勾配法＝勾配を下って極小値を探す方法

多峰性関数に勾配法は使えない！

【解決策】

１．ランダムに動く
２．多点で探す6

【ランダムサーチ・網羅探索による最適化】
モンテカルロ法とも呼ばれる
も ばれ
全数探索と同義な場合と区別する場合がある
・最も単純な方法で、どんな問題に対しても適用可能
最も単純な方法で どんな問題に対しても適用可能
・時間さえかければ、理屈の上ではいつかは最適解を発見できる・今まで生成した解候補と評価値を参考にしないので、探索効率は悪い
探索性能の下限としての比較対象

ランダムサーチによる最適化のプロセス：
新しい解候補をランダムに生成

解候補の評価

それまで見つけた最良の解候補よりも高い評価値だった場合、最良解候補を置き換える7

【ランダムサーチの探索性能は？】10個のパラメータから構成されるある評価関数について、最適解を含むある領域を一様分布を利用したランダムサーチによって探索する。
このとき、ランダムサーチの１回の解候補生成において全パラメータが最適なパラメータの±5%以内となる確率と、ランダムサ チによ てそのような解を生成するまでにかかる解候補生成回数の期待値
ランダムサーチによってそのような解を生成するまでにかかる解候補生成回数の期待値を計算せよ。 ただし最適解は探索領域の端部 5%以内には無いものとする。
最適パラメータの±5%の領域
最適パラメ
タの±5%の領域＝10%ある1次元の領域（100%）
ヒント１： ランダムサーチの探索1回あたりで最適パラメータの±5%以内の値になる確率は？

8

【ランダムサーチの探索性能は？】10個のパラメータから構成されるある評価関数について、最適解を含むある領域を一様分布を利用したランダムサーチによって探索する。
このとき、ランダムサーチの１回の解候補生成において全パラメータが最適なパラメータの±5%以内となる確率と、ランダムサ チによ てそのような解を生成するまでにかかる解候補生成回数の期待値
ランダムサーチによってそのような解を生成するまでにかかる解候補生成回数の期待値を計算せよ。 ただし最適解は探索領域の端部 5%以内には無いものとする。
最適パラメータの±5%の領域
最適パラメ
タの±5%の領域＝10%ある1次元の領域（100%）
ヒント１： ランダムサーチの探索1回あたりで最適パラメータの±5%以内の値になる確率は？
101p  10 

よって初めて上記のような解を生成するまでの解候補生成回数の期待値は、確率 の事象が発生するまでの回数（幾何分布）の期待値に等しいから確率pの事象が発生するまでの回数（幾何分布）の期待値に等しいから
19
 1010p

【焼きなまし法(simulated annealing)による最適化】
統計物理学からのアナロジー：
統計物理学からのアナロジ
・物質の分子結晶はいろいろな状態を遷移・エネルギーの低い状態へ遷移しやすい傾向
・温度を下げるとエネルギー極小状態へ到達
温度を下げると ネ ギ 極小状態 到達・温度を一気に下げると局所的極小状態へ・温度をゆっくり下げると大域的極小状態へ

「状態」
→ 解候補
「エネルギー」 → 評価値
エネルギー最小状態＝最適解とした最適化のプロセスシステムがある状態 xi になっている確率
自由エ
エネルギー

exp( f ( xi ) T )P( xi | T )  exp( f ( x) T )xS

ギブス分布（ボルツマン分布）
Tは温度パラメータ

f (x)

状態空間

x

エネルギーの低い状態ほど確率が高い

特徴：
１）時間をかければかけるほど解候補が改善され、十分時間をかければ最適解の発見が保障される
２）連続関数最適化，および組合せ問題最適化の両方に適用できる10

【焼きなまし法(simulated annealing)による最適化】
処理手順
１） 初期解候補 x を生成
２） 何らかの確率分布によって解候補 x から次の解候補 x ' を生成（ランダムに x の近傍を選ぶことが多い）
３） 次の確率に従って乱数Zを生成：
ただしZは0または1の値をとる
ただし温度関数４）

T (t ) kln(t  2)

1
P( Z  1)  f ( x' )  f ( x) 1  exp
T (t )
Z  1 のとき x を x' に置換える、

５）時刻 t を1進めて２）から繰返す

これでは収束が遅すぎて非実用的k/t の関数にすることもある
11

「近傍」の解候補の生成方法
１） 数値最適化問題の場合： 解候補 x は連続値パラメータなので、中心値 x , 分散 2 の正規分布によって生成： x '  N ( x,  2 )分散のパラメータは設計者が問題に応じて適当に与える
時間とともに小さくしていく場合もある

x' x２） 組合せ最適化の場合：
まず解候補の「近傍」を定義する必要がある→ 問題毎に個別に設定： 実行可能な表現になるよう工夫を要する
例１） ＴＳＰ
順列中の都市を入れ替える
順列中の部分順列を移動するなど

ＡＢＣＤＥＦ → ＡＥＣＤＢＦＡＢＣＤＥＦ → ＡＣＤＥＢＦ

例２） ＳＡＴ
論理変数の論理を反転させる

１０１１１０１ → １００１１０１

12

【参考】

大腸菌 走性大腸菌の走性2～3秒間の直線的運動と、0.1秒ほどの方向転換の繰り返し

化学的濃度が低い（要するにエサが少ない）場合には活発に動き回り、濃度の高い領域ではあまり動かないようにすることでエサのある場所へ集まる

焼きなまし法によく似た探索をしている

13

勾配法＝勾配を下って極小値を探す方法

多峰性関数に勾配法は使えない！

【解決策】

１．ランダムに動く
２．多点で探す14

【制約のない関数最適化】 滑降シンプレックス法
関数の勾配を計算せず、多点探索で最小値を得る
関数の勾配を計算せず多点探索で最小値を得る
N次元変数関数に対してN+1コ以上のシンプレックス（探索点） x j について以下を定義
1) 目的関数値が最大の点 xh← 関数最小化なので、最悪の点
2) 目的関数値が2番目に大きい点 xs
3) 目的関数値が最小の点 xl← 関数最小化なので、最良点
4) $i = h$なる点を除いた全ての x j の重心 x g

xh を以下の xr で置き換える：
xr  (1   ) x g  xh ，ただし   0 は反射係数(操作2：拡大) x g  xr 方向に沿って xr を以下の xe に置き換える：  0 は拡大係数，ただしxe  xr  (1   ) x g(操作3：縮小) xh を以下の xc で置き換える：
xc   xh  (1   ) x g ，，ただし 0    1 は縮小係数
縮 係数(操作4：収縮) シンプレックス全体を xl の方向へ半分に縮小する，ただし i  1,  n  1xi  0.5( xl  xi )(操作1：反射)

これらの操作を、次に述べる手順で組合せ，シンプレックスを更新する。
基本的に，最悪点をシンプレックス重心の逆側へ移動して関数値を小さくしていく。
が  1,   0.5,   2 が経験的に良い

15

滑降シンプレックス法の動作

最悪点

xh〔スタート〕

・シンプレックスを構成する探索点の中で、プ評価の最悪点、２番目に悪い点、最良点
シンプレックス重心座標を求める

xl

最良点2番目に悪い点

xg
xs

最悪点を除いた全シンプレックスの重心

16

滑降シンプレックス法の動作〔まずは反射操作〕

最悪点

xh

xl

最良点2番目に悪い点

xg
xs
〔反射〕

xr
17

滑降シンプレックス法の動作

最悪点

xh

xg
2番目に悪い点〔反射の結果が最良点より良ければ 拡大操作〕
良ければ、拡大操作〕

xl

最良点

xs
xr

xe

〔拡大〕

18

滑降シンプレックス法の動作

最悪点・拡大の結果が反射の結果拡大の結果が反射の結果より悪ければ、最悪点を反射点と置き換えて最初へxhxg
2番目に悪い点・拡大の結果が反射の結果より良ければ、最悪点を拡大点と置き換えて最初へ

xl

最良点

xs
xr

xe

〔拡大〕

19

滑降シンプレックス法の動作

最悪点〔反射の結果が最悪点より悪ければ 縮小操作〕
悪ければ、縮小操作〕

xh〔縮小〕

xcxg
2番目に悪い点

xl

最良点

xsxr

xe20

滑降シンプレックス法の動作

最悪点

・さもなければ〔収縮〕操作
さもなければ〔収縮〕操作
新たに生成したシンプレックス〔縮小〕 の関数値を評価して最初へ

xh

xcxg
2番目に悪い点・縮小の結果が最悪点より良ければ、最悪点を縮小点と置き換えて最初へ

xl

最良点

xs
xr

xe21

滑降シンプレックス法の動作

最悪点

xh

〔収縮〕

2番目に悪い点

xl

最良点

xs22

シンプレックス全体を最良点の方向へ移動 （最良点は動かない）

滑降シンプレックス法の動作

最悪点

xh

〔反射の結果が最悪点よりましだが２番目に悪い点より悪い
場合、最悪点を反射点に置換
えて縮小操作〕

〔縮小〕

xl

最良点

xg
2番目に悪い点

xsxc
xr

xe23

滑降シンプレックス法の動作

最悪点

xh・縮小の結果が最悪点より良ければ、最悪点を縮小点と置き換えて最初へ・さもなければ〔収縮〕操作
さもなければ〔収縮〕操作
新たに生成したシンプレックス〔縮小〕 の関数値を評価して最初へ

xl

最良点

xg
2番目に悪い点

xsxc
xr

xe

24

滑降シンプレックス法の動作

最悪点〔反射の結果が２番目に悪い点よりましだが最良点ほどではない場合、最悪点を反射点に置き換えて最初へもどる〕

xh

xl

最良点

2番目に悪い点

xg
xs
〔反射〕

xr
25

【例：冗長自由度アームのリーチング問題】●各関節の角度をパラメ タとして、ア●各関節の角度をパラメータとしてアーム先端がターゲット座標に一致するパラメータ
ム先端がタ ゲット座標に 致するパラメ タを求める問題。ただし、各関節はリンク同士を±90度以内で繋ぐ制約があるので、なるべく関節の可動角度中央付近の角度にすることが好ましい。
コスト関数＝（アーム先端座標とターゲット座標の偏差の2乗）＋（各関節の可動中央角度からの偏差の2乗）

局所解が高々数個の場合、滑降シンプレックス法
による探索が最も適する

26

３Ｄモデルと画像認識による特徴点データとのマッチング目的• モデルとの比較による建造進行状況把握
• 拡張現実における正確な位置での表示
評価対象• 特徴点データとモデルデータの距離を最小外れ値

誤差範囲内
マ チング
マッチング
PTAMからの特徴点データ

3Dモデルデータ

３Ｄモデルと画像認識による特徴点データとのマッチング

立体物の画像特徴点の抽出

元データ
デ タ

特徴点座標

滑降シンプレックス法
マッチング• コストは特徴点とモデルとの距離の二乗• モデルを平行移動・回転・拡大縮小して特徴点集合へ重ねる (x,y,z,α,β,γ,size)の7パラメータ探索• マトリックス計算とクォータニオンを併用
トリ クス計算とクタ オンを併用

探索後

３Ｄモデルと画像認識による特徴点データとのマッチング

コスト関数の景観：多峰性

xyx-yxz
x-z

y-z

α-β

このくらいの多峰性でも、滑降シンプレックス法による探索が十分通用する

α-γ

β-γ

29

滑降シンプレックス法と最小メジアン法（ロバスト推定法の 種）の組合せ
最小メジアン法（ロバスト推定法の一種）の組合せ

外れ値を含む特徴点データとのマッチング

模型2b

模型2a3Dモデル

【復習】 ロバスト推定法 ・データに例外値（ノイズ）を多く含む場合、これらを自動的に除去した回帰を行う
最小メジアン法 (LMedS：Least Median of Squares)• ランダムに幾つかのサンプルを抽出し、最小二乗法に当てはめることを繰り返す
LMedS  min med  i2
• LMedS基準• 全データの半分が外れ値でも大きな影響を及ぼさない回帰された各直線と全データとの間の誤差の２乗を計算し 昇の間の誤差の２乗を計算し、昇順に並べ、メディアンになるデータの２乗誤差値を読取る

この値が最小のモデルが最良
ランダムに取り出したサンプルに対し多重回帰を行う

最小メジアン法におけるノイズ処理
ランダムに取り出したサンプルに対し多重回帰を行う

YX

パラメータ（（ここでは1次元で表現）
では 次元で表現）

最小メジアン法におけるノイズ処理
１）ランダムに取り出したサンプルに対し多重回帰を行う

Y

２）全データとモデルの誤差を計算
X

パラメータ（（ここでは1次元で表現）
では 次元で表現）

最小メジアン法におけるノイズ処理
１）ランダムに取り出したサンプルに対し多重回帰を行う

Y

２）全データとモデルの誤差を計算
X

パラメータ
（
（ここでは1次元で表現）
では 次元で表現）

最小メジアン法におけるノイズ処理
１）ランダムに取り出したサンプルに対し多重回帰を行う

Y

２）全データとモデルの誤差を計算
X

パラメータ
（
（ここでは1次元で表現）
では 次元で表現）

最小メジアン法におけるノイズ処理
１）ランダムに取り出したサンプルに対し多重回帰を行う

２）各回帰モデルにおいて全データとモデルの誤差を計算
３）各回帰モデルにおける全データのモデルとの誤差を大きさ順に並べる
４）誤差の中間値（メディアン）が最小のモデルを採択

緑色のモデルを採択

この部分の誤差は必ずメディアン値以下

ノイズデータや例外値を除外
イズデ タや例外値を除外
中間値
メディアン

最小メジアン法におけるノイズ処理
１）ランダムに取り出したサンプルに対し多重回帰

滑降シンプレックス法を用いてサンプルと３Ｄモデル表面の距離が最小になる角度等を行うのパラメータを探索

２）各回帰モデルにおいて全データとモデルの誤差を計算
３）各回帰モデルにおける全データのモデルとの誤差を大きさ順に並べる
４）誤差の中間値（メディアン）が最小のモデルを採択

緑色のモデルを採択

この部分の誤差は必ずメディアン値以下

ノイズデータや例外値を除外
イズデ タや例外値を除外
中間値
メディアン

外れ値を含む特徴点データとのマッチング
模型2b(特徴点数：1393)
外れ値除外
探索成功

F=15
q=1000

外れ値除外
探索成功

結果
• 外れ値を良く除外• 全データ数に対し外れ値の割合が少ない

F=30
q=1000今回の講義のまとめ
勾配法 勾配を下 て極小値を探す方法
勾配法＝勾配を下って極小値を探す方法

多峰性関数に勾配法は使えない！

【解決策】

探索性能の下限

１．ランダムに動く → ランダムサーチ・ＳＡ
２．多点で探す

次回→ 滑降シンプレックス法
パラメータ数が5～6個程度または局所解が十数個程度ならお奨め

両方の特徴を併せ持つ最適化手法：遺伝的アルゴリズム
39参考文献
１）長尾 智晴： 最適化アルゴリズム，昭晃堂（２０００）．
２）W.H.Press, B.P.Flannery, S.A.Teukolsky and W.T.Vetterling:
ニューメリカルレシピ・イン・シー
カ
ピ
Ｃ言語による数値計算のレシピ，技術評論社（１９９３）．
言語
る数値計算
ピ 技術評論社
３）小野 功，佐藤 浩，小林 重信：
単峰性正規分布交叉ＵＮＤＸを用いた実数値ＧＡによる関数最適化，
単峰性正規分布交叉ＵＮＤＸを用いた実数値ＧＡによる関数最適化
人工知能学会誌 Vol.14, No.6, pp.1146—1155 (1999)．
４）樋口隆英，筒井茂義，山村雅幸
４）樋口隆英筒井茂義 山村雅幸 : 実数値GAにおけるシンプレックス交叉の提案 ,人工知能学会誌 , Vol.16 , No.1 , pp146-155 , 2001

40