Hokkaido University of Education
Title
多点局所探索に基づく大域的最適化アルゴリズム（１）
Author(s)
金光, 秀雄; 今野, 英明; 高橋, 伸幸
Citation
北海道教育大学紀要, 自然科学編, 58(2): 1-10
Issue Date
2008/02
URL
http://s-ir.sap.hokkyodai.ac.jp/dspace/handle/123456789/66
Rights
ABSTRACT
We will report on the development of a multi-start local search algorithm. At first, we developed a local search algorithm and tested this algorithm. Next, we investigated and developed a simple multi-start local search algorithm. Finally, we developed and demonstrated multi level single linkage algorithm as an effective multi-start local search algorithms.
1 はじめに
与えられた関数をある制約条件のもとでの最小（大） 化する解を求める最適化問題は，工学，理学，経済学 などの多くの分野で登場する．このような問題の中で 唯一の解が存在する問題については，降下法の枠組み の中で理論やアルゴリズム面から十分深い研究がな され，非常に大規模な問題まで解けるようになってい る．一方，複数の局所最適解を有する多峰性関数を最 適化する問題も数多く存在しており，このような問題 の大域的最適解を求める解法は，決定論的な方法 [8] と確率論的な方法 [4, 7, 20] に大別される．決定論的 な方法は，解への収束が確実に保証され，多くの研究 がおこなわれているが，問題の先見情報を要請し，問 題構造によっては計算量の指数的な増大を招く恐れが ある．確率論的な手法の中で複数の初期点から局所 最適化を行う多点局所探索法（多スタート局所探索 法）は，探索領域に複数の点を生成・選択し，選択さ れた点に局所探索法を適用するもので，古くから研 究され，比較的少ない関数評価回数で高精度に大域的 最適解を求められ [14]，収束性についても研究されて いる [5] ．さらに，確率論的な枠組みの中で，国内で は SA (Simulated Anealing) 法 [13] や GA (Genetic Algorithm)法 [1, 2, 15] などの新しい方法の研究が盛 んである．しかし，これらの方法は十分な収束性が示 されておらず，問題構造との関わりが不明なパラメー タ設定が多く，この設定でのユーザ負担が大きい．一 方，多点局所探索法は欧州を中心に研究されてきた が，SA 法や GA 法の研究が盛んな国内ではこの方法 は殆ど研究されていない．本稿では，多点局所探索法 について，収束性やアルゴリズムの観点から検討し， それを実装して動作の確認をすることを目的とする． 第 2 章では，準備として，本稿で扱う問題や仮定 および記法について述べ，第 3 章では，局所探索法の アルゴリズムと簡単な実行例を示す．第 4 章では，多 点局所探索の概要を示し，単純な多点局所探索法にお いて大域的最適解に収束する確率を与えたときに必要 な初期点の数について検討し，簡単な実行例を示す． 第 5 章では，多点局所探索法の中で有効といわれる 多段階単連結（multi-level single linkage; MLSL) 法 を示し，幾つかのテスト関数に対して簡単な数値実験 をおこない，その結果を示す．2 準備
2.1 問題と仮定 ここでは，次のような大域的最適化問題 (PG) を 扱う． (PG) Minimize f (x) ≡ f (x1, x2, · · · , xn) Subject to Li≤ xi≤ Ui, (i = 0, 1, . . . , n) この問題において，最小化する関数 f(x) は目的関 数（objective function）と呼ばれ，制約条件は各変数 の上下限制約である．この問題を上下限制約条件付き 大域的最適化問題と呼ぶことがある．この制約条件を 満足する領域を S とすると，S は次式で表される． S = { x ≡ (x1, x2, · · · , xn)| Li≤ xi≤ Ui} この S は実行可能領域（feasible region）と呼ばれる． この問題における孤立極小点は，点 x∗を中心とし た半径 δ > 0 の開球 B(x∗, δ) が存在し， ∀x ∈ (B(x∗, δ) ∩ S) \ {x∗}, f(x∗) < f (x) (1) を満足するする点 x∗である．この問題における目的 関数 f(x) は十分滑らかで，全ての極小点は有限個の 孤立極小点 x∗1, x∗2, . . . , x∗M からなると仮定し，この 解の集合を X∗ とし，この集合の添字集合を IM と すると，これらは次のように表される． X∗={ x∗1, x∗2, . . . , x∗M} (2) IM ={ 1, 2, . . . , M } (3) このとき，各極小点 x∗i への収束領域を D(x∗i)と し，その収束領域の相対的大きさ（relative size）を η(x∗i) = µ(D(x∗i)) / µ(S), (i = 1, 2, . . . , M ) とする (µ(·) はルベーク測度)．ここで収束領域とは，局所探 索法 Lopt が初期点 x(0)i から出発すれば，ある極小 点 x∗i ∈ X∗ へ収束するような初期点すべてからなる S の部分集合を意味する．すなわち，D(x∗i)は次式 で表される． D(x∗i) ={ x | x∗i ← Lopt(f, S, x(0)i ), x(0)i ∈ S } (4) また，任意の初期点 x(0)i に対し Lopt を適用したと き X∗ の中の極小点のいずれかに必ず収束すること を仮定すると，x∗i への収束領域 D(x∗i)やその収束領 域の相対的大きさ η(x∗i)に対して，次式が成立する． D(x∗i)∩ D(x∗j) =∅, ∀i = ∀j ∈ IM D(x∗1)∪ D(x∗2)∪ · · · ∪ D(x∗M) = S η(x∗1) + η(x∗2) +· · · + η(x∗M) = 1 (5) ここで，大域的最適解（最小点）が唯一存在するこ とを仮定し，これを x∗∗ とすると，この点 x∗∗ は孤 立局所最適解の中で関数の値が最も小さい点であるか ら，次のように表すことができる． f∗∗= min{ f(x∗1), f (x∗2), . . . , f (x∗M)} x∗∗= argmin{ f(x∗1), f (x∗2), . . . , f (x∗M)} (6) 2.2 記法 本稿では，以降の記述において，次の記法を用いる． 1. 右辺の式または変数 Eval の値を左辺の変数 v に代入することを，“v ← Eval; ” のように表す． 2. p 個の入力パラメータ a1, a2, . . . , apを与えて，q 個の出力 b1, b2, . . . , bq が得られる名前 A name のアルゴリズムを次のように表す． (b1, b2, . . . , bq)← A name (a1, a2, . . . , ap) なお，出力が一つの場合は，出力パラメータの括 弧 (·) を省略する． 3. アルゴリズムの選択，繰り返しを以下のように 記す． • 選択： if 条件 then 文 T1; . . . , 文 Tr; else 文 F1; . . . , 文 Fs; fi ; これは，“条件” が成立すれば，“ then ” 以 降の部分：“ 文 T1; . . . , 文 Tr; ” を実行し， 条件が成立しなければ，“ else ” 以降の部分： “ 文 F1; . . . , 文 Fs; ”を実行することを示 している．“条件” が成立しない場合がなけれ ば，“ else 文 F1; . . . , 文 Fs; ”は記述しな い．また，実行する 文 が 1 つだけであれば， “ fi ; ” を省略する． • 繰り返し： for 変数 ← 初期値 to 変数の終値 do 文 L1; . . . 文 Lt; od ; これは，変数を初期値から終値まで 1 ずつ増 やしながら，“ do . . . od ; ” で囲まれた部分： “文 L1; . . . 文 Lt; ”を繰り返すことを示し ている．これも，実行する文が 1 つであれば， “od ; ” を省略する． 4. 有限の要素からなる集合 A に対して，A の j 番 目の要素を Aj で，その集合 A の個数を |A| で 表す．3 局所探索法
3.1 問題と最適性条件 局所探索法は，与えられた問題の極小点（局所最適 解）を求める方法である．この方法は最適化法の中で は局所最適化法と呼ばれ，本稿で対象とする目的関数 が複数の極小点を有する多峰性関数の場合は非線形な 関数となるので，非線形最適化法とも呼ばれる [9, 21]． ここでは，まず制約無し局所最適化問題について述 べる．この問題 (PL)は，次のように定式化される． (PL) Local Minimize f(x) (7) この問題において，点 x∗が制約無し最適化問題の 極小点で，f が x∗ の近傍で微分可能であるならば ∇f(x∗) = 0 (8) という 1 次の必要条件が成立する． 3.2 降下法 制約無し局所最適化問題を解く方法の枠組みとして 降下法（descent method）がある．これは，目的関数 f と初期点 x0∈ Rn を与え，降下する方向への探索 を繰り返して，極小点 x∗ とその関数値 f∗ を求める 方法で，アルゴリズムの概要は以下のようになる． アルゴリズム 1（降下法） (f∗, x∗)← Loptd(f, x(0)) Ld1) 初期点 x(0) の関数値を評価する．繰り返し 回数 k を k ← 0 とする． Ld2) 探索方向 d(k) を降下方向に定める． Ld3) 現在の点 x(k) から降下方向へ関数値が減少 する点を探索する． Ld4) 探索して見つけられた点を新たな点 x(k+1) とし，収束条件を満足したら終了．そうでなけれ ば，繰り返し回数を k ← k + 1 として，ステップ Ld2) へ戻る． このような降下法において，ステップ Ld2) の降 下方向 d(k) の定め方によって，最急降下法，ニュー トン法，共役勾配法，準ニュートン法などがあるが， ここでは，これらの中で最も有効とされる BFGS 更 新公式による準ニュートン法について述べる． 準ニュートン法では，ステップ Ld2) の探索方向 d(k)を定めるために，現在の点 x(k)の関数値 f(x(k)) とこの点の勾配ベクトル∇f(x(k))∈ Rn および次式 で表される正定値対称行列 B(k)∈ Rn×n の情報を用 いる．具体的な d(k) の定め方は，3.4 で述べる． for ∀y = 0, yTB(k)y > 0. (9) 3.3 直線探索 降下法のステップ Ld3) の探索は直線探索（line search）と呼ばれ，これはある点 x(k) から与えられ た方向 d(k)で定まる直線上で，関数 f(p(k)+ αd(k)) を α に関して最小化して，そのステップ幅 α(k)を求 める手法である．直線探索には正確な最小点を求める 正確な直線探索法 (exact line search) と，ある基準 も満足する範囲内でのステップ幅を求める正確でない 直線探索 (inexact line search) がある．ここでは，正 確でない直線探索の一つとして，より少ない関数評価 で局所探索の収束の保証が得られる Armijo の基準を 用いた線形探索法を用いる．現在の点 x(k)と探索方 向 d(k)およびパラメータ 0 < ξ < 1/2, 0 < τ < 1 を 与えて，Armijo の基準を満たすステップ α(k)を求め るアルゴリズムを以下に示す． アルゴリズム 2（Armijo の基準による直線探索法） α(k)← LSrchA(x(k), d(k), ξ, τ ) LS1) 初期ステップ幅 α0← 1 とし，繰り返し回数 を i ← 0 とする． LS2) ξ に対して，次のような Armijo の基準 f (x(k)+ αid(k))≤ f(x(k)) + ξαi∇f(x(k))Td(k) を満足するならば，α(k)← αi として直線探索を 終了．そうでなければ，ステップ LS3) へ進む． LS3) ステップ幅を αi+1← ταi とし，繰り返し回 数を i ← i + 1 として，ステップ LS2) へ戻る． 3.4 準ニュートン法 直線探索を用いた準ニュートン法（BFGS 公式）の アルゴリズムについて述べる．これは，目的関数 f， 初期点 x(0)，正数 η および 直線探索に与えるパラ メータ (ξ, τ ) を与え，解 x∗ とその関数値 f∗求める もので，アルゴリズムは次のようになる． アルゴリズム 3（準ニュートン法） (f∗, x∗)← Loptq(f, x(0), η, ξ, τ ) Lq1) 正定値対称行列 H(0)∈ Rn×n を定め，繰り 返し回数 k を k ← 0 とする．Lq2) もし ∇f (x(k)) = 0ならば終了．そうでなけ れば，次式により探索方向 d(k) を求める． d(k)=−H(k)∇f(x(k)) Lq3) 次の直線探索によりステップ幅 α(k)を求め， 次の点 x(k+1)← x(k)+ α(k)d(k)を定める． α(k)← LSrchA(x(k), d(k), ξ, τ ) Lq4) 行列 H(k) の正定値性を数値的に保つた め に ，(y(k))Ts(k) ≥ η y(k) s(k) ならば H(k+1)← H(k) とし，そうでなければ，次の BFGS公式により行列 H(k+1) を定める． H(k+1)← H(k)+ 1 +(y (k))TH(k)y(k) (s(k))Ty(k) · s(k)(s(k))T (s(k))Ty(k) −H(k)y(k)(s(k))T + s(k)(H(k)y(k))T (s(k))Ty(k) 但し，s(k)= x(k+1)− x(k), y(k)=∇f(x(k+1))− ∇f(x(k)). k ← k + 1 とおいてステップ Lq2) へ戻る. このアルゴリズムで，H(0) は通常単位行列に取られ， H(k)は (9) 式の B(k)の逆行列である． 3.5 プログラムの実装と数値例 プログラムは C 言語で作成され，以下のファイル に分けて実装した． 1) defc.h: 定数の定義 2) defp.h: 関数プログラムのプロトタイプ宣言 3) util.h: マクロの定義 4) bfgs.c: BFGS 公式による準ニュートン法 5) gfdif.c: 勾配∇f(x) の差分近似（後述） 6) lsearch.c: 直線探索法 7) prnt.c: 結果の出力 これらの中で，4)∼7) については，茨木・福島の書物 [9]に書かれている FORTRAN のプログラムをもと に，C 言語で一部修正を加えて作成した． また，勾配ベクトル∇f(x(k))の計算は，多峰性関 数のように複雑な関数を求める場合は計算が困難にな る場合があるので，解析的な計算の代わりに前進差分 による近似で行った．具体的には，∇f(x(k)) の第 i 成分 ∂f (∂xx(k)i )は，次式によって計算した． ∂f (x(k)) ∂xi ≈ f (x(k)+ δei)− f(x(k)) δ (10) ここで，eiは，i 番目の成分だけが 1 で残りの成分が 0の単位ベクトルであり，この例では，δ = 1 × 10−7 とした．また，アルゴリズムに与えるパラメータは， それぞれ，η = 10−4, ξ = 10−3, τ = 0.5 とした． 以上のようにして実装したプログラムの動作の一例 として，次の Rosenbrock 関数を取り上げる． f (x) = f (x1, x2) = 100(x2− x21)2+ (1− x1)2 なお，この問題において，初期点は x0= (−1.2, 1) で その関数値は f0= 24.2 である．この関数は，唯一の 極小点 x∗= (1, 1) を有し，その関数値は f∗ = 0で ある． この問題における解への収束過程を表 1 と 図 1 に 示す． 表 1：各繰り返し回数 k における点の座標 (x1, x2)， 関数値 f(x1, x2)および関数評価回数 Nf e k x1 x2 f (x1, x2) Nf e 0 −1.20000 1.00000 24.20000000 3 5 −0.48526 0.20704 2.28688083 24 10 −0.12268 −0.00361 1.29523422 40 15 0.25398 0.05557 0.56453009 58 20 0.55483 0.29243 0.22191510 75 25 0.82242 0.67364 0.03228395 90 30 0.98146 0.96262 0.00038514 106 35 1.00001 1.00001 5.06×10−11 125 36 1.00000 1.00000 3.07×10−14 130 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 x2 x1 "rosenbrk.xy" using 1:2 "../prog/lopt/ort.d" using 2:3 図 1: Rosenbrock 関数における解への収束過程
4 多点局所探索法
4.1 アルゴリズム
本稿では，局所最適化手法 Lopt を初期点を変えな がら適当な回数だけ適用する多点局所最適化法（mul-tistart local search method; MSLS）について考える．
この手法では，試行回数 N を与えて，最小値 ˜f∗∗ と 最小点 ˜x∗∗ を求めるもので，そのアルゴリズムの概 要は，次の通りである． アルゴリズム 4（多点局所探索法 1） ( ˜f∗∗, ˜x∗∗)← MSLS1(f, S, N) M1) 探索領域内に一様に N 個の点をサンプルし， これらの点での関数値を評価する． M2) サンプルした点から一定の基準で点を選択 する． M3) 選択されたそれぞれの点を開始点として局所 探索をおこなう． M4) 局所探索法で得られた点の中で関数値が最も 小さいものを ˜f∗∗ とし，それを与える点を ˜x∗∗ として，出力する． 上記のステップ M2) で，特別な基準を設けずに 「すべての点を」選択し，それに局所探索を適用した 場合のアルゴリズムを示し，解を見出す確率について 以降で検討する． 関数 f ，実行可能領域 S，最小点を見出す確率 P ， 最小点の収束領域の相対的大きさ η(x∗∗)を与えたと き，最小値 ˜f∗∗と 大域的最小点 ˜x∗∗ を求める多点局 所最適化手法 M SLS2 のアルゴリズムは，次のよう になる． アルゴリズム 5（多点局所探索法 2） ( ˜f∗∗, ˜x∗∗)← MSLS2(f, S, P, η(x∗∗)) ML1) 試行回数 N を，確率 P と最小点の収束領域 の相対的大きさ η(x∗∗)から求める（詳細は 4.2 を参照）．最小値の初期値を ˜f∗∗← + ∞ とする． ML2) 次のように N 回繰り返す． for i = 1 to N do ML3) S 内に一様乱数により点を生成する手続き Rand を用いて点 x(0)i を x(0)i ← Rand(S) ; のように発生させる． ML4) 点 x(0)i を初期点として局所最適化手法 Lopt を次のように適用し，局所最小解 x∗i とその 点での最適値 fi∗ を求める． (fi∗, x∗i)← Lopt(f, S, x(0)i ) ML5) Lopt で得られた fi∗, x∗i から，これまで の大域的最小値 ˜f∗∗ と大域的最小解 ˜x∗∗ を 次のように更新する． if fi∗< ˜f∗∗ then ˜f∗∗← fi∗; ˜x∗∗← x∗i; fi ; od ; ML6) ˜f∗∗ と ˜x∗∗ を結果として出力する． return ( ˜f∗∗, ˜x∗∗) . このアルゴリズムでは，局所最適化法が N 回繰り 返されるので，局所最適化法の平均計算量を L とす ると，全体の計算量は，O(N L) となる． 4.2 大域的最小点が確率 P 以上で求められるため に必要なサンプル数 アルゴリズム 5 のステップ ML2) において，大域 的最小点を確率 P 以上で求めるために必要な試行回 数 N を求めてみよう．このアルゴリズムにおいて， 1つのサンプル点が大域的最小点の収束領域に入る確 率 Pinと入らない確率 Pout は， Pin= η(x∗∗) Pout= 1− η(x∗∗) (11) で与えられる．したがって N 回の試行においてすべ てのサンプル点が大域的最小点の収束領域に入らない 確率 PN out は， PoutN = (1− η(x∗∗))N (12) となる．よって N 回の試行において少なくとも 1 つ のサンプル点が大域的最小点の収束領域に入る確率 PN in は 1− PoutN となるから，(12) 式より PinN = 1− (1 − η(x∗∗))N (13) 上式が P となるために必要な試行回数 N は，1 − (1− η(x∗∗))N = P を解くことで N = log[1log(1− η(x− P )∗∗)] (14) となる [6]．上式の試行回数の右辺値は一般に小数点 付きの値になるが，もともと N は確率 P 以上での 試行回数なので，小数点以下の部分は切り上げること にする．したがって， N = log(1− P ) log[1− η(x∗∗)] (15) となる．4.1 の アルゴリズム 5 のステップ ML1) で，試行回数 N を求める部分は，(15) 式によって求 めればよい．
最小点の収束領域の相対的大きさ η(x∗∗)と 最小点 を見出す確率 P における試行回数 N を表 2 に示す． この表から，η(x∗∗) = 0.3 で 99.9% の確率で最小点 を求めたいときには，試行回数 N は 20 回必要とな り，η(x∗∗) = 0.01 で 99.5% の確率で最小点を求めた いときには，試行回数 N は 528 回必要となることが わかる． 通常，多峰性関数の大域的最適化で解く問題の困難 さが議論されることがある．困難さの要因として，問 題の次元が高いときや峰の数が多いとき，問題が困難 になるといわれているが，表 2 からわかるように，多 点局所探索法においては，最小点の収束領域の相対的 大きさ η(x∗∗)が小さいときが，本質的に困難な問題 となることがわかる． 表2：最小点の収束領域の相対的大きさη(x∗∗) と最小点を見出す確率P との関係 XXXXX確率P XXXXη(x∗∗) 0.500 0.400 0.300 0.200 0.900 4 5 7 11 0.950 5 6 9 14 0.990 7 10 13 21 0.995 8 11 15 24 0.999 10 14 20 31 XXXXXXXXX 確率P η(x∗∗) 0.100 0.050 0.020 0.010 0.900 22 45 114 230 0.950 29 59 149 299 0.990 44 90 228 459 0.995 51 104 263 528 0.999 66 135 342 688 XXXXXXXXX 確率P η(x∗∗) 0.005 0.002 0.001 0.0005 0.900 460 1151 2302 4605 0.950 598 1497 2995 5990 0.990 919 2301 4603 9209 0.995 1058 2647 5296 10594 0.999 1379 3451 6905 13813 4.3 プログラムの実装と数値実験 M SLS2 のアルゴリズムは，以下の C 言語のプロ グラムで実装した． 1) msls2.c: M SLS2 の関数プログラム 2) randnv.c: ステップ ML3 の Rand(S) を用 いて S 内に xi を求める． 3) mtrand.c: Mersenne Twister による [0, 1)
区間の乱数を 1 個発生させる． これらの中で，msls2.c の中で randnv.c の関数を 呼び出しており，さらに randnv.c の中で mtrand.c の関数を呼び出している．局所探索法は，前節の BFGS 公式による準ニュートン法のプログラム（bfgs.c）を 用いている． 3)の Mersenne Twister は，松本等 [12] によって 提案された乱数生成法で，219937− 1 の長い周期と 623次元で均等分布することが知られており，乱数の 発生も比較的高速である．また，この乱数で生成され た乱数列は，diehard test，load test，ultimate load
testなどの統計検定の全てをパスしており，信頼性も
高い．
次の 2 変数の Six Hump Camel Back 関数でテス トした． f (x1, x2) =4x21− 2.1x41+ x1x2+13x61− 4x22+ 4x42 − 2.3 ≥ x1≥ 2.3, −1.5 ≥ x2≥ 1.5 このテスト関数は 6 つの局所最適解を有し, 次のよう な 2 つの最小点 x∗∗1 , x∗∗2 と最小値 f∗∗ をもつ． f∗∗=−1.031628 x∗∗1 = (0.089842, −0.712656), x∗∗2 = (−0.089842, 0.712656) この関数の概形を 図 2 に示す． f(x,y) -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 x1 -1.5 -1 -0.5 0 0.5 1 1.5 x2 -5 0 5 10 15 20 25
図 2：Six Hump Camel Back 関数
この実験で用いたパラメータは，P = 0.99，η = 0.3 で，残りの局所最適化のパラメータについては，3.5 の実験と同様である．以上のように実装したプログ ラムとパラメータを与えて，上記で述べたテスト関 数に適用したところ，13 回の局所探索（N = 13）を 実行し，最小値 f∗∗ = −1.0316285 と大域的最適解 x∗∗= (0.0898420, −0.7126564) を得ることができた． なお，これに要した関数評価回数は Nf e = 492回で ある．また，13 回の試行中 2 つの大域的最適解に収 束したのは 5 回であった． 関数の概形を 図 2 に示す． f(x,y) -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 x1 -1.5 -1 -0.5 0 0.5 1 1.5 x2 -5 0 5 10 15 20 25 図 S H C l B k関数
この関数に，M SLS2 のアルゴリズムを適用した 結果の挙動を図 3 に示す． -1.5 -1 -0.5 0 0.5 1 1.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 x2 x1 "6humpc.xy" using 1:2 "../prog/msls/or.d" using 2:3 "../prog/msls/ors.d" using 2:3
図 3：Six Hump Camel Back 関数に M SLS2 アルゴ リズムを適用した結果の挙動 図 3 で，実行可能領域の境界に近い点を初期点 x(0) としたとき，その次の点 x(1) が大きく離れて，近接 の極小点を飛び越えて収束する場合がある．これは， 局所探索 Loptq において，実行可能領域の境界でこの 関数の勾配が非常に大きいため，ステップ Lq2) で， d(0)=−H(0)∇f(x(0)) (H(0):単位行列) で求まる方 向ベクトルの大きさ d(0) が非常に大きくなり，ス テップ Lq3) の x(1)← x(0)+ α(0)d(0)(通常 α(0)= 1) で求まる x(1)が x(0) と遠く離れてしまうことが原因 となっていると考えられる．これを避けるためには， 一般に d(k) の大きさを制限する必要がある．例え ば，0 < φ < 1 として， d(k)← min{ d(k) , φ · min 1≤i≤n{Ui− Li} } とするのが良いであろう．このとき φ は，問題の構 造にも依存するが，今回の問題の場合には φ = 0.1 程 度で良いと思われる．一般には，解く問題の各極小点 での単峰領域1に内接する円の直径を φ の目安にすれ ばよいであろう． また，図 3 の挙動において始めの 6 つの初期点と 局所探索法で最終的に得られた結果を表 3 に示す． この表において，各局所探索法における関数評価回 数は 23∼44 回の範囲にあり，関数値については，概 ね 8 桁の精度で合致している．13 回全部の各局所探 索法における関数評価回数は 23∼50 回の範囲にあり， 評価回数の平均は約 32.8 回である．
表 3: Six Hump Camel Back 関数に M SLS2 アル ゴリズムを適用した結果の挙動（一部） k x1 x2 f(x1, x2) Nfe 0 −1.156583 −0.721309 2.22678238 3 5 −0.089837 0.712653 −1.03162845 24 0 −1.788813 1.186345 2.38999690 3 7 −1.607118 −0.568656 2.10425031 29 0 2.229322 0.724524 9.54600241 3 9 −1.703606 0.796084 −0.21546382 38 0 1.293896 −1.162203 2.76594998 3 11 1.607142 0.568635 2.10425032 44 0 −1.432259 0.268663 1.59322922 3 7 −1.703606 0.796085 −0.21546382 30 0 0.442093 0.663022 0.01176557 3 5 −0.089846 0.712675 −1.03162845 23
5 多段階単連結法
5.1 アルゴリズム 前節の図 3 でわかるように，M SLS2 のアルゴリ ズムでは，同じ局所最適解に重複して収束するという 問題がある．この問題を解消する方法としてクラスタ リング法 [16] がある．このアルゴリズムの概要は次 のようになる． アルゴリズム 6 (多点局所探索法 3 (クラスタ法)) ( ˜f∗∗, ˜x∗∗)← MSLSc(S, f, N ) MC1) 探索領域内に一様に N 個の点をサンプル し，これらの点での関数を評価する． MC2) 選択された点から一定の基準で点を選択 する． MC3) 選択された点に対してクラスタリングし， 各クラスタを代表する点を開始点として局所探 索をおこなう． MC4) 局所探索法で得られた点の中で関数値が最 も小さいものを ˜f∗∗とし，それを与える点を ˜x∗∗ として，出力する． クラスタリング法では，点の集まりであるクラスタ から代表する点一つを局所最適化の開始点とするた め，同じ局所最適解に重複して収束するという状況 は，大幅に緩和されると考えられる．しかしながら， 点の集まりを把握するクラスタリング法には多くの方 法が存在するのが現状であり，利用する方法により構 成されるクラスタは異なり，クラスタリング法自身の アルゴリズムも比較的複雑なものが多い．多段階単連結（Multi Level Single Linckage; 以降
MLSLと呼ぶ）法は，Timmer によって提案され [19]，
Rinnooy Kanと Timmer の共著論文 [17] により広く 知られるようになり，多点局所探索に基づく手法で最 も効率的な手法として知られている [14]． このアルゴリズムにおいて，k 回目の繰り返しにお ける局所探索手続きは，他のサンプル点が注目してい るサンプル点を中心とした閾値 rk 内に入る場合を除 いたあらゆる点に対して，適用される．この閾値 rk は次式で表される． rk= 1 √ π σ · µ(S) · Γ 1 +n 2 log N N 1/n (16) ここで, σ は正定数, µ(S) は実行可能領域 S のルベー ク測度, Γ(·) はガンマ関数および N はサンプル点の 数である． さらに, このアルゴリズムは，予想される極小点の 個数 [5]: |X∗|(γN−1) / (γN−|X∗|−2) が，すでに見 出された異なる極小点の個数|X∗| に 0.5 を加えた値 を超えないところまで，サンプリングを繰り返す. し たがって, このアルゴリズムは，もし次式が成立すれ ば停止する． |X∗| (γN − 1) γN − |X∗| − 2 ≤ |X ∗| + 0.5 (17) ここで，γ (0 < γ < 1) は，サンプルした点を減少させ る割合である．すなわち，γN 個の小さい関数値を持 つ点だけが，(16) 式に基づく局所探索に用いられる． 最初の繰り返しで, N(0)（= N ）が与えられ，その 後の繰り返しでは，N(k+1)(k = 0, 1, . . .) と N は以 下のように更新される． N(k+1)←1 γ(2|X ∗|2+ 3|X∗| + 2) , N ← N + N(k+1) この MLSL アルゴリズムは，探索領域 S 上の関数 f ，初期サンプル点 N ，サンプルの減少率 γ (0 < γ < 1)および正数 σ を与えて，局所最適解の集合 X∗ と 極小関数値の集合 F∗を求める．このアルゴリズムの ステップは次のようになる． アルゴリズム 7（多点局所探索法 4（MLSL 法）） (F∗, X∗)← MLSL (S, f, N(0), γ, σ) ; ML1) 初期化として，次のように設定する． X∗← ∅ ; F∗← ∅ ; k ← 1 ; N ← N(0). ML2) N(k) 個の点 xi(i = 1, 2, . . . , N(k))を一様 乱数 Rand(S) を用いて発生させ，これらの関数 値を評価する． ML3) サンプルされた点の中で，関数値の小さな 点を γN(k)個選択し，Nr← γN(k) とする． ML4) rkを (16) 式で決定し， f(xj) < f (xi)かつ xj − xi < rk を満足する点 xj が存在する ような点を除いたそれぞれの新しい点 xi (i = 1, 2, . . . , Nr)に局所探索を適用し，得られた局所 最適解 x∗ をその集合 X∗ に加える． ML5) もし停止条件 (17) が満足されるならば，停 止する．そうでなければ，新しいサンプルサイ ズを N(k+1)←1γ(2|X∗|2+ 3|X∗| + 2) に設定し, N ← N + N(k+1) とし, k ← k + 1 として, ステッ プ ML2) に戻る． このアルゴリズムは大変単純で, クラスタリング法 も必要ない．このアルゴリズムの困難さは，新しい点 に対する近傍点を求めることに計算量を要することで ある． また，このアルゴリズムはその単純さにもかかわら ず, 理論的性質は非常に強い．もし σ ≥ 4 ならば, た とえサンプリングが永遠に続いても，局所探索の総数 は確率 1 で有限である [17]． 5.2 プログラムの実装と数値実験 M LSL のアルゴリズムは，mlsl.c という C 言語 のプログラムで実装されている．ステップ ML2 の一 様乱数 Rand(S) を発生させて点 xi を得るプログラ ムは，4.3 と同様に，randnv.c と mtrand.c を用い ている．局所探索法は，3.5 で示した BFGS 公式に よる準ニュートン法のプログラム（bfgs.c）を用い ている． 以上のプログラムを実装し，数値実験した結果を示 す．実験に用いたテスト関数は以下の (1)∼(3) に示す 2変数関数で，実験で用いたパラメータは，N(1)= 40， γ = 0.25，σ = 4 である．
(1) Six Hump Camel Back関数
この関数に要した最終的なサンプル数は N = 69 となり，5 回の局所探索を実行し，最小値 f∗∗ = −1.0316285 と大域的最適解 x∗∗ = (−0.0898421, −0.7126564) を得ることができた． なお，これに要した関数評価回数は Nf e = 229 回である．この関数評価回数は，前節の M SLS2 のアルゴリズムの Nf e= 492回と比較して，半 分以下の関数評価回数で解を見出していること から，M SLS2 のアルゴリズムより効率がよい ことがわかる．
(2) Goldstein and Price関数
表される． f (x) ={(1 + (x1+ x2+ 1)2(19− 14x1+ 3x21− 14x2+ 6x1x2+ 3x22)}{30 + (2x1− 3x2)2 (18− 32x1+ 12x1x1+ 48x2− 36x1x2 + 27x22)} − 2 ≥ x1≥ 2, −2 ≥ x2≥ 2 この関数は，4 つの局所最適解を有し，最小点 x∗∗ = (0, −1) で，その関数値が f∗∗ = 3 と なる． この関数の探索に要した最終的なサンプル数は N = 56 となり，5 回の局所探索を実行し，最 小値 f∗∗ = 3.0000000 と大域的最適解 x∗∗ = (0.0000000, −1.0000000) が得られた．この結果 に要した関数評価回数は Nf e= 334回である． (3) Branin’s RCOS関数 この関数の形とその実行可能領域は，次のように 表される． f (x) = x2− 5.1 4π2x 2 1+ 5 πx1− 6 2 + 10 1− 1 8π cos(x1) + 10 − 5 ≥ x1≥ 10, 0 ≥ x2≥ 15, この関数は，3 つの局所最適解を有し，これら 3 つ がいずれも最小点となり，これらの最小点での関数値 は f∗∗= 5/(4π) ≈ 0.397887 となる． この関数の探索に要した最終的なサンプル数は N = 136 となり，この探索の中で 5 回の局所探索をお こない，最小値 f∗∗ = 0.3978874 と大域的最適解 x∗∗ = (9.4247780, 2.4750000) を得た．この結果に 要した関数評価回数は Nf e= 204回である．
6 おわりに
多点局所探索に基づく大域的最適化法の収束性やア ルゴリズムの検討について，若干解説的な内容を含ん でいるが，その大枠を述べた．このような解説的な記 述を含めたのは，この方法が長い伝統を有し理論的な 背景も充実しており，海外でも評価が高い [14, 20] に もかかわらず，国内では殆ど研究・開発・利用されて いない現状で，この方法をもう少し国内でも広めてみ たいと考えたからである．多点局所探索法は，基本ア ルゴリズム自身は非常に単純なので，多くの理論的な 考察も得られている [5]．単純な多点局所探索法にお いて，通常は点のサンプル数を入力パラメータとして 利用者が与えるが，本稿では，最小点への収束領域の 相対的大きさと最小点を見出す確率を入力パラメー タとして与え，それらをもとにサンプル数を求めるア ルゴリズムを提案した．パラメータの数は一つ増える が，両者のパラメータは，問題の構造に直接関わるパ ラメータである．また，両者のパラメータからサンプ ル数を具体的な数値の表として求めた．これから，多 点局所探索法を用いる場合，問題の困難さが最小点へ の収束領域の相対的大きさが本質的であることが明ら かになった．また，多くの多点局所探索法では，局所 探索法をブラックボックスとして利用しているので， それが大域的最適化法の中でどのような挙動をしてい るのかについて，明らかにされていなかった．今回の 例で，多峰性関数を有する問題において局所探索法が かなり複雑な挙動を示しており，大域的最適化法と組 み合わせるときは，初期のステップ幅を制限する必要 があることが明らかになった． 多点局所探索の中で特に有効とされる [14] MLSL 法は，特に西欧でその後の研究が広くおこなわれた [11, 18]が，特に原論文 [17, 19] でのアルゴリズムの 実装や数値実験の結果に不明確なところも多いとい う問題をいくつか抱えている．本稿では，このような 点も少しでも明らかにしたいという意図もあり，アル ゴリズムとして原論文より少し明確な記述をした．但 し，原論文のアルゴリズムで実装されている近傍点探 索のステップについては，原論文では Bentley のア ルゴリズム [3] を実装しているが，このアルゴリズム を実装すると高次元での空間計算量が膨大になること が T¨orn 等により指摘されている [20] ので，このス テップは実装しなかった． 今後の課題として，本稿で示した多点局所探索法 は，解が実行可能領域の境界に存在するときに，その 解を見出すことができない．この問題を解消するため に，解が境界上に存在するときにも，その解を見出せ るように局所探索法を改善する必要がある．また，今 回は小規模な問題でのテストしかおこなわなかったた め，今後はより大規模な問題に対するテストを行い， 他の方法と比較していく必要がある．7 謝辞
本稿での研究成果の一部は，学長裁量経費の個人研 究支援の経費によった．この研究機会と支援を頂いた ことについて，謝意を表します．参考文献
[1] 足立進，澤井秀文：遺伝子重複に基づく進化的計 算法：関数最適化への適用, 情報処理学会論文誌， Vol.49, No.11, 2663–2671, (2001). [2] 相澤彰子: スキーマ貪欲な遺伝的探索アルゴリズ ムの構成, 電子情報通信学会論文誌，Vol. J 78-D-II, No.1, 94–104, (1995).[3] Bentley, J.L., Weide, A.C. and Yao, A.C.: Opti-mal expected-time algorithms for closest point problems, ACM Transactions on Mathematical Software, Vol.6, 563–580, (1980).
[4] Boender, C.G.E., Rinnooy Kan, A.H.G. and Timmer, G.T.: A Stochastic Method for Glo-bal Optimization, Mathematical Programming, Vol.22, 125–140, (1982).
[5] Boender, C.G.E., Rinnooy Kan, A.H.G. and Timmer, G.T.: Beyesian stopping rules for mul-tistart global optimization methods, Mathe-matical Programming, Vol.37, 59–80, (1987).
[6] Brooks, S.H.: A Discussion of Random Meth-ods for Seeking Maxima, Operations Research, Vol.6, 244–251, (1958).
[7] Dzemyda, G., ˇSaltenis, V. and ˇZilinskas, A.: Stochastic and Global Optimization, Kluwer Acad. Publ., Dordrecht, pp.233, (2002).
[8] Floudas, C.A.: Deterministic Global Optimi-zation —Theory, Methods and Applications—, Kluwer A.P. , pp.739, (2000). [9] 茨木俊秀，福島雅夫: FORTRAN 77 による最適 化プログラミング, 岩波書店, pp.480, (1991). [10] 金光秀雄, 宮腰政明, 新保勝: 極小値集合による 単峰性および多峰性関数の定義とその性質, 電 子情報通信学会論文誌, Vol.J79-A, No.11, 1844– 1851, (1996).
[11] Locatelli, M. and Schoen, F. : Random Linkage: a family of acceptance/rejection algorithms for global optimisation, Mathematical Program-ming, Vol.85, 379–396, (1999).
[12] Matsumoto, M. and Nishimura, T. : Mersenne Twister: A 623-Dimensionally Equidistributed Uniform Pseudo-Random Number Generator, ACM Transactions on Modeling and Computer Simulation, Vol.8, No.1, 3–30, (1998).
[13] 三木光範, 廣安知之, 笠井誠之, 小野景子: 適応的
近傍を持つ温度並列シミュレーテッドアニーリン グ, 情報処理学会論文誌，Vol.42, No.4, 745–753, (2001).
[14] Nemhouser, G.L, Rinnooy Kan, A.H.G. and Todd, M.J. (Eds.): Optimization, Handbooks in Operations Research and Management Sci-ence Vol.1, North-Holland, pp.709, (1989).
[15] 小野功, 山村雅幸，喜多一: 実数値 GA とその
応用, 人工知能学会誌，Vol.15, No.2, 259–266, (2000).
[16] Rinnooy Kan, A.H.G. and Timmer, G.T. : Stochastic Global Optimization Methods. Part I: Clustering Methods, Mathematical Program-ming, Vol.39, 27–56, (1987).
[17] Rinnooy Kan, A.H.G. and Timmer, G.T. : Stochastic Global Optimization Methods. Part II: Multi Level Methods, Mathematical Pro-gramming, Vol.39, 57–78, (1987).
[18] Sch¨oen, F. : Random and Quasi-Random Link-age Methods in Global Optimization, Journal of Global Optimization, Vol.13, No.4, 445-454, (1998).
[19] Timmer, G.T.: Global optimization: A stochas-tic approach, Ph.D. Dissertation, Econometric Institute, Erasmus University, Rotterdam, Hol-land, pp.175, (1984).
[20] T¨orn, A. and ˇZilinskas, A.：Global Optimiza-tion, Lecture Notes in Compter. Science. 350, Springer-Verlag, p.255, (1989). [21] 矢部博: 最適化とその応用, 数理工学社, pp.263, (2006). （金光 秀雄 函館校教授） （今野 英明 函館校准教授） （高橋 伸幸 函館校教授）
