混沌系工学特論 配布資料 #5
担当 : 井上 純一 (情報科学研究科棟 8-13)
URL : http://chaosweb.complex.eng.hokudai.ac.jp/˜j inoue/

平成 16 年 11 月 22 日

目次
2

2.1

2.2

2

57
ノイズの効果を利用したアルゴリズム . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

最適化問題と統計力学

2.1.1
2.1.2

エネルギー 4 準位を持つ 2 体系の最適化 . . . . . . . . . . . . . . . . . . . . . . . . 59

2.1.3
2.1.4

最適化問題 : E(s1 , s2 ) = −J s1 s2 − h1 s1 − h2 s2 再考 . . . . . . . . . . . . . . . . . 67

マスター方程式とその定常解. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

ノイズレベルのスケジューリング . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

2.1.5 最適なノイズレベルのスケジューリング . . . . . . . . . . . . . . . . . . . . . . . . 70
平衡状態と物理量の期待値 : 内部エネルギー, 自由エネルギー, エントロピー etc. . . . . . . 73

最適化問題と統計力学
N 個の正の実定数 {a1 , a2 , · · · , aN } が与えられたとき, それらを A1 と A2 の 2 つの組に分け,


 




aj −
aj 



aj ∈A1(1)

aj ∈A2を最小にする問題を考えよう. これは個々のアイテム aj にラベル sj ∈ {+1, −1} を割り当て, アイテム ajがグループ A1 に属するのであれば, sj = +1, グループ A2 に属するのであれば, sj = −1 と約束すること
により次の量を最小化する問題に書き換えることができる.

 N




aj sj E(s) = 



(2)

j=1

これは 2 分割問題と呼ばれ, 典型的な組み合わせ最適化問題の一つである. この量 E(s) を最小にするアイ
テムの割り振り方 s = (s1 , s2 , · · · , sN ) とそのときの E(s) を求めるのがここでの課題となる.
さて, この種の問題を解かなければならない状況は実社会に数多く見受けられる. 例えば最近, プロ野球
では近鉄バッファローズとオリックスブルーウエーブが合併してオリックスバッファローズができ, さらに楽天ゴールデンイーグルス1 が誕生し, 支配下選手 70 名ほどをオリックス (A1 ) とゴールデンイーグル
ス (A2 ) で分配しなければならなくなった. 選手によって年棒額 (aj : j = 1, · · · , 70) もまちまちであるか

1 余談だが「楽天イーグルス」だけでいいものを, 何でまた「ゴールデン」という形容詞まで前にくっつけてしまったのであろう

か ? 感性は人ぞれぞれかも知れないが, この「ゴールデン」という語感が安っぽく響いてきてならない. 日本では動物系の名前が圧倒的に多いが決まりでもあるのだろうか. ちょっと野暮ったくは無いか ? それに比べてメジャーリーグでは「レッドソック
ス」などなかなか洒落た名前が多いと思う.

57

大学院 混沌系工学特論 2004

最適化問題の数理 part 1 (担当 : 井上 純一)

ら, 経営者側としては両球団で支払う総年棒額ができるだけ同じになるように s1 , s2 , · · · , s70 に ±1 を割り振り, 選手を両球団 A1 , A2 に分配したい. どうしたら良いだろう ? という具合に問題は設定されるわけである. アイテム数が 10 個程度の問題であれば, この解を探すことはさほど難しくはないが, アイテム数の増加とともに問題は困難になって, アイテム数が 70 であれば, 全件探索したとして 270 通り. これは非常に困難を極める. 従って, ここで我々が十分に時間をかけて取り組むにたる問題となるのである.
本節からこの種の問題に統計力学の立場から取り組むにはどうするのか, また, それから何がわかるのか,をできる限り平易に説明していく. 前節の [連想記憶の数理] でも述べたが, 統計力学ではミクロな量である
神経素子 Si の動きを追うのではなく, 重なり m = (1/N ) i ξi1 Si 等のマクロな量に着目し, この量を通じてシステムの振る舞いを調べた. これをここでの問題に当てはめるのであれば, エネルギー関数 (あるいはコスト関数) E(s) を通じてシステムの振る舞いを調べることになる. これは例えば「アイテム数の増加とともに, エネルギー関数の最小値がどのように振舞うか」という問いに対する答えを提供することになるであろう. それはつまるところ, 理論上どこまでエネルギーを小さくできるか, ということであるから, あるアルゴリズムを用いてエネルギー値が求まったときに, そのアルゴリズムはどの程度 (精度的に) 良いア
ルゴリズムなのかを, そのアルゴリズムは理論上の最小エネルギーにどの程度まで迫れているのか, という観点から評価できることになるであろう.
しかし, [連想記憶の数理] で扱った問題とやや異なるのは, 神経系での記憶の問題では重なりというマクロな量の振る舞いさえわかってしまえば連想記憶という機能について (我々の数理モデルの範囲で) 知ることができて, 個々の素子自体についての情報まで得られなくとも特に問題は無かったのではあるが, この最適化問題では実際にミクロなラベル自体の情報も極めて重要になってくるという点である. 「両球団が支払う総年棒の差の最小値」がわかっても, 実際にどの選手がどちらの球団に行ったのか s = (s1 , s2 , · · · , s70 )がわからなければ選手もファンも納得しないであろう. というわけで, ここで我々が学習する内容は大きくわけて次の 2 つである.

• アルゴリズムの精度評価の基準となるエネルギーの最小値等, マクロな量の解析的評価法• 具体的にエネルギーの最小値を与えるミクロな状態を求めるアルゴリズム
もちろん, この 2 つは互いに関連している. 1 番目の方法を説明するにはまだ準備が必要なので, まずはじめに 2 番目の方法について見て行こう.

2.1

ノイズの効果を利用したアルゴリズム

前節の連想記憶の最後に見たように, シナプス結合が対称である場合には神経回路網のエネルギー関数が定義できてE(S) =−1 
wij Si Sj
N ij

(3)

がそのエネルギー関数であった. そして, このエネルギー関数は素子のノイズが無い場合の非同期的状態
更新 :
Si (t + 1)

=

sgn 




wij Sj (t)

(4)

j=i

ここは 58 ページ目

大学院 混沌系工学特論 2004

最適化問題の数理 part 1 (担当 : 井上 純一)で単調に減少することを学んだ. 従って, 問題が (3) 式の最小化として定式化される最適化問題に対する有効な∗アルゴリズムとして (4) 式が適用でき, 更新式 (4) 式を十分長い時間まわした挙句の収束点 S ∗ = (S1∗ , · · · , SN)をその答えとすることができるであろう. しかし, このアルゴリズム (4) 式でエネルギーは単調にしか減少
しないわけであるから, エネルギー関数 (3) が複数の谷を持つ場合にはその谷が最小値であるという保障はない. この意味で, こうした局所的な最小を回避するようなメカニズムをアルゴリズムの中に組み込む必要がありそうだ. そこで, 非常に簡単なシステムを考えることにより, この点を詳しく見て行こう.

2.1.1

エネルギー 4 準位を持つ 2 体系の最適化

いきなり, (3) 式のような多体系のエネルギー関数を扱うのではなく, 問題を恐ろしく簡単にしてしまおう. そこで次のような 2 体系のエネルギー関数の最小値を探す問題を考える.

E(s) =−J s1 s2 − h1 s1 − h2 s2(5)

ここで, s = (s1 , s2 ) であり, それぞれの成分は ±1 をとるものとする. この問題は簡単である. なぜならば,全ての状態 s の組み合わせは (s1 , s2 ) = (1, 1), (1, −1), (−1, 1), (−1, −1) の 4 つしかない. 従って, それぞれに対応するエネルギーの値が簡単に求まって

ε1≡

E(1, 1) = −J − h1 − h2

(6)

ε2

≡

E(1, −1) = J − h1 + h2

(7)

ε3

≡

E(−1, 1) = J + h1 − h2

(8)

ε4

≡

E(−1, −1) = −J + h1 + h2

(9)がそれらの値である. 従って, J > h1 > h2 , 例えば, J = 1, h1 = 0.5, h2 = 0.1 のように J1 , h1 , h2 を決めれば, ε1 , · · · , ε4 はその大きさで順序付けられて(10)

ε1 < ε4 < ε2 < ε3がそのエネルギー準位となる. これを模式的に描くと図 1 のようになる.

ここで, 我々が当初調べたかった問題 (2) と比べてかなりトーンダウンしてしまったではないか, と思われるかもしれない. たしかにこの 4 準位系は簡単ではあるが, 多数の要素からなるシステムに現れる現象である「相転移現象」を除けば我々がこれから直面する問題, そして, それをどのように回避するか等は全て
この簡単な問題 (5) を考えることから明らかになるので, まずはこれに集中しよう. 難しい問題 (2) はこれが十分理解できた後であらためて取り組むことにする. 何事も順序が大切なのである.
さて, ここでノイズを利用したアルゴリズムを次のように構成する.

✬✩

ノイズを利用したアルゴリズム :(i) 各時刻で任意に s1 , s2 のどちらか 1 つを選び, その符号を変える. 符号を変える前後の状態を s, sと書こう.
(ii) (i) の前後でのエネルギー差 : ∆E = E(s ) − E(s) を計算し, ∆E < 0, つまり, (i) の操作でエネ


ルギーが下がるのであれば, 新しい状態として s を採用し, ∆E > 0, つまり, エネルギーが増加してしまった場合でも e−∆E/T の確率で新しい状態 s を採用する.

(iii) 上記 (i)(ii) を繰り返す.
✫✪

ここは 59 ページ目

大学院 混沌系工学特論 2004

最適化問題の数理 part 1 (担当 : 井上 純一)

E

ε3

ε2
ε4
ε1

(1,1)

(1,-1)(-1,1)(-1,-1)S

図 1: エネルギー 4 準位.

このアルゴリズムは今後 [ノイズを利用したアルゴリズム] として本稿に現れ, 今回の内容の中心となるのでよく覚えておいてほしい. 上記の太字で書いた部分がノイズによる効果であり, そのノイズの大きさは T でコントロールされる. この効果により, エネルギーは単調には減少せず, ときどき増加することになる. このアルゴリズムを実際に計算機上で実行しょう. その際, 各時刻で状態が (s1 , s2 ) = (1, 1), (1, −1), (−1, 1), (−1, −1)の 4 つのそれぞれをとる確率 : P1 , P2 , P3 , P4 を各ステップで評価するために, 同じエネルギー関数 (5) 式を持つ複数の系 (アンサンブル) を用意し, (システムの中でその状態をとる数)/(アンサンブルのシステム総数)をもって, 各確率を評価することにしょう. この評価では例えば P1 は

P1=

状態 (1, 1) をとるシステムの数
アンサンブルを構成するシステムの総数(11)として見積もられることになる. 大して長いプログラムではないので, 下記に上記アルゴリズムのプログラ
ム例を載せておく. これは混沌系工学特論のホームページ :http://chaosweb.complex.eng.hokudai.ac.jp/˜ j inoue/KONTON2004/konton2004.html
からダウンロードできるので, 各自が必要とあらばこれを適時書き換えて実行してみると良いと思う.

/***************************************************//*
Algorithm for 2-spin problem*//*
J.Inoue*//***************************************************/
ここは 60 ページ目

大学院 混沌系工学特論 2004

最適化問題の数理 part 1 (担当 : 井上 純一)/***************************************************/#include<math.h>#include<stdio.h>#include<stdlib.h>/*************************************************//*--------Definitions of Parameters--------------*//*************************************************/#define N 10000/* アンサンブルのシステム総数 */#define p 0.7#define pp 0.7/* s1 の初期状態を選ぶための確率 *//* s2 の初期状態を選ぶための確率 */#define tmax 500 /* アルゴリズムをまわす回数の上限 */#define T 1.01/* ノイズレベル. アニーリングする場合にはこれをコメントアウト */#define Z (exp((J+h1+h2)/T)+exp((-J+h1-h2)/T)+exp((-J-h1+h2)/T)+exp((J-h1-h2)/T))#define uu0 (exp((J+h1+h2)/T)/Z) /* e1 をとる定常状態確率 */#define ud0 (exp((-J+h1-h2)/T)/Z) /* e2 をとる定常状態確率 */#define du0 (exp((-J+h2-h1)/T)/Z) /* e3 をとる定常状態確率 */#define dd0 (exp((J-h1-h2)/T)/Z)#define J 1.0/* e4 をとる定常状態確率 */#define h1 0.5#define h2 0.1/*****************************************************************//*------ "Seeds" of rundomnumbers appearing in this program -----*//*****************************************************************/#define SEEDNOISE 121#define SEEDNOISE2 142#define SEEDNOISE3 132#define SEEDNOISE4 45#define SEEDNOISE5 92/*****************************************************************/
int s1[N];int s2[N];double uu;
double ud;
double du;double dd;// double T; /* アニールする場合にはこのコメントをはずす *//*******************************************************//*
Function of generating randomnumber*//*
by Numerical recepies in C ran3*//*******************************************************/#define MBIG 1000000000#define MSEED 161803398#define MZ 0
#define FAC (1.0/MBIG)

ここは 61 ページ目

大学院 混沌系工学特論 2004

最適化問題の数理 part 1 (担当 : 井上 純一)

float ran3(long *idum){
static int inext,inextp;
static long ma[56];static int iff=0;
long mj,mk;
int i,ii,k;if (*idum < 0 || iff == 0) {iff=1;
mj=MSEED-(*idum < 0 ? -*idum : *idum);
mj %= MBIG;
ma[55]=mj;
mk=1;
for (i=1;i<=54;i++) {ii=(21*i) % 55;
ma[ii]=mk;
mk=mj-mk;if (mk < MZ) mk += MBIG;
mj=ma[ii];}for (k=1;k<=4;k++)for (i=1;i<=55;i++) {ma[i] -= ma[1+(i+30) % 55];
if (ma[i] < MZ) ma[i] += MBIG;}
inext=0;
inextp=31;*idum= 1;}
if (++inext == 56) inext=1;
if (++inextp == 56) inextp=1;
mj=ma[inext]-ma[inextp];
if (mj < MZ) mj += MBIG;
ma[inext]=mj;
return mj*FAC;}#undef MBIG
#undef MSEED#undef MZ#undef FAC/* (C) Copr. 1986-92 Numerical Recipes Software 1+5-5i. */ここは 62 ページ目

大学院 混沌系工学特論 2004

最適化問題の数理 part 1 (担当 : 井上 純一)/* s1 の初期値を確率的に選ぶ */void set_initial(){
long idum=(-SEEDNOISE2);
int x;
for(x = 0; x < N-1;x++){if(p > ran3(&idum)){s1[x] = 1;
}else{
s1[x] = -1;}}}/* s2 の初期値を確率的に選ぶ */void set_initial2(){
long idum=(-SEEDNOISE3);
int x;
for(x = 0; x < N-1; x++){if(p > ran3(&idum)){s2[x] = 1;
}else{
s2[x] = -1;}}}/* 状態更新部分 */void metro(int k){long idum=(-SEEDNOISE+k);long idum2=(-SEEDNOISE5+k);
int x;
double r;for(x = 0; x < N-1; x++){r = ran3(&idum2);
if((r <= 0.5) && (ran3(&idum)< exp(-(2.0*J*s1[x]*s2[x]+2.0*h1*s1[x])/T))){s1[x] = -s1[x];}else{s1[x]= s1[x];}if((r > 0.5) && (ran3(&idum)< exp(-(2.0*J*s1[x]*s2[x]+2.0*h2*s2[x])/T))){s2[x] = -s2[x];}else{

ここは 63 ページ目

大学院 混沌系工学特論 2004

最適化問題の数理 part 1 (担当 : 井上 純一)

s2[x]= s2[x];}}}/*****************************************************//*
Main Program*//*****************************************************/
main(){FILE *pt;
int i,t;set_initial();set_initial2();if((pt = fopen("result.dat", "wt")) != NULL){
for(i=1,uu=0.0,ud=0.0,du=0.0,dd=0.0;i < N-1; i++){if((s1[i]==1) && (s2[i]==1)){uu = uu + 1.0/N;}
if((s1[i]==1) && (s2[i]==-1)){ud = ud + 1.0/N;}
if((s1[i]==-1) && (s2[i]==1)){du = du + 1.0/N;}
if((s1[i]==-1) && (s2[i]==-1)){dd = dd + 1.0/N;}}
fprintf(pt,"%d %lf %lf %lf %lf %lf %lf %lf %lf\n",0,uu,ud,du,dd,uu0,ud0,du0,dd0);for(t=0; t < tmax; t=t+1){/* アニールする場合には下記コメントをはずす */// T = 2.0/(log(1.0+2.0*t));// T = 2.0/(1.0+2.0*t);
metro(t);for(i=1,uu=0.0,ud=0.0,du=0.0,dd=0.0; i < N-1; i++){if((s1[i]==1) && (s2[i]==1)){uu = uu + 1.0/N;}
if((s1[i]==1) && (s2[i]==-1)){ud = ud + 1.0/N;}
if((s1[i]==-1) && (s2[i]==1)){du = du + 1.0/N;}
if((s1[i]==-1) && (s2[i]==-1)){dd = dd + 1.0/N;}}
fprintf(pt,"%d %lf %lf %lf %lf %lf %lf %lf %lf\n",t,uu,ud,du,dd,uu0,ud0,du0,dd0);}}
fclose(pt);}
結果の一例を図 2 に載せる. この図より, 十分に時間が経過すると, それぞれの状態をとる確率が定常値に収束していく様子がわかる. なお, この定常値への収束が初期値に依存するかどうかを見るために, s1 , s2の初期値を各アンサンブルで変えて (s1 , s2 ) = (1, 1) となる確率 P1 (グラフ中では P (uu) と表記) の時間変化をプロットしたものが図 3 である. この図から明らかに定常確率 (分布) は初期状態の選び方に依らないことがわかる. この性質, つまり, 定常分布が初期状態に依らない性質をエルゴード性と呼び, この性質が満たされていれば, そのシステムの物理量 At の長時間平均は定常分布を作るアンサンブルでの平均 A に

ここは 64 ページ目

大学院 混沌系工学特論 2004

最適化問題の数理 part 1 (担当 : 井上 純一)

0.7
P(uu)

0.6

0.5

0.4

0.3

P(dd)

0.2
P(ud)

0.1

P(du)

0
0

50

100

150

200

250

300

350

400

450

500

t

図 2: (s1 , s2 ) = (1, 1), (1, −1), (−1, 1), (−1, −1) をとる確率 P (uu) = P1 (E = ε1 ), P (ud) = P2 (E = ε2 ), P (du) = P3 (E =ε3 ), P (dd) = P4 (E = ε4 ) の時間変化. アンサンブルのシステム総数 N = 10000, ノイズレベル T = 1. 横に引かれた破線はマスター方程式 : (14)(15)(16)(17) の定常解である平衡分布 (確率). これらの定常確率はエネルギー 4 準位を反映している. つまり,
ε1 < ε4 < ε2 < ε3 の順に P1 > P4 > P2 > P3 となっている.

等しくなる. つまり, 式で書けば
T

A

=

1 
At
T →∞ T
t=0
lim

(12)

が成り立つ. 従って, 定常分布さえ書き下すことができれば, 後は学部で習った確率・統計の知識を用いて,その分布での期待値を計算することにより, 注目する一つのアンサンブルの物理量に関する長時間平均が計算できるということになる2 . 具体例は後に見ることにして, 次に, ここでの計算機シミュレーションで得られた定常分布が解析的に得られないかどうかを検討することにしょう.

2.1.2

マスター方程式とその定常解

図 2 で得られた計算機シミュレーションでの定常分布が解析的に得られれば何かと都合が良いであろう.
そこで, これを調べるために「確率の保存式」であるマスター方程式を用いる. これは一般的にdP (s)dt=s =sw(s → s)P (s ) −w(s → s )P (s)

(13)

s=s



と書ける. ここで, w(s → s) は状態 s から状態 s への遷移確率を表す. 従って, 上式 (13) の右辺第 1 項


は s 以外の全ての状態 s から状態 s に流入する確率の流れを, また, 第 2 項は状態 s から, それ以外の状態

2 式 (12) で長時間平均を計算する際にはサンプリング ((12) 式中の和) の開始時刻をアルゴリズムのスタート時刻 t = 0 ではな

く, システムが平衡状態に落ち込んだ時刻 t = t0 とする場合がある. T の大きさが十分に大きければ平衡状態へ落ち込むまでのサンプリング点の効果はほとんど結果に効かないであろうと予想できるが, 現実の計算機シミュレーションでは T を有限値で置き換えるので, この t0 の設定は重要になってくる場合がある.

ここは 65 ページ目

大学院 混沌系工学特論 2004

最適化問題の数理 part 1 (担当 : 井上 純一)

0.9
0.8
0.7
P(uu)

0.6
0.5
0.4
0.3
0.2
0.1
0
0

50

100

150

200

250

300

350

400

450

500

t

図 3: (s1 , s2 ) = (1, 1) をとる確率 P (uu) の時間変化を初期状態を変えてプロット. アンサンブルのシステム数総 N = 10000, ノイズレベル T = 1.
s へ出て行く確率の流れを表しており, この両者の差がシステムが状態 s をとる確率の変化率を決める, という意味を (13) 式は持っている. そこで, このマスター方程式を我々が議論してきたエネルギー関数が (5)式で与えられ, 状態が [ノイズを利用したアルゴリズム] で更新されるシステムに適用すれば, 4 状態のそれぞれをとる確率は次のマスター方程式に従う.
dP1= P2 + P3 − e−2(J+h2 )/T + e−2(J+h1 )/T P1
dt

dP2= P3 − 2 + e−2(h1 −h2 )/T P2 + e−2(J+h2 )/T P1 + e−2(J−h1 )/T P4
dt
dP3= e−2(h1 −h2 )/T P2 − 2P3 + e−2(J−h2 )/T P4
dt

dP4= P2 + P3 − e−2(J−h1 )/T + e−2(J−h2 )/T P4
dt

(14)
(15)
(16)
(17)

ここで, 当然 P1 + P2 + P3 + P4 = 1 が成り立っていることに注意しなければならない.
さて, システムが定常状態になっているのであれば, (dP1 /dt) = (dP2 /dt) = (dP3 /dt) = (dP4 /dt) = 0 が成り立つはずであるから, これから得られる方程式 (14),(15),(16),(17) を P1 , P2 , P3 , P4 に関して解くこと
により, 多少面倒な計算の結果

P1

=

P2

=

P3

=

P4

=

P (1, 1) = 

s1 =±1

P (1, −1) = P (−1, 1) = e(J+h1 +h2 )/Te(J+h1 +h2 )/T=(Js1 s2 +h1 s1 +h2 s2 )/TZ
s2 =±1 ee(−J+h1 −h2 )/Te(−J+h1 −h2 )/T
=
(Js
s
+h
s
+h
s
)/T
1 2
1 1
2 2
Z
s2 =±1 e(19)e(−J−h1 +h2 )/Te(−J−h1 +h2 )/T=(Js1 s2 +h1 s1 +h2 s2 )/TZ
s2 =±1 e(20)e(J−h1 −h2 )/Te(J−h1 −h2 )/T=(Js1 s2 +h1 s1 +h2 s2 )/TZ
s2 =±1 e(21)

s1 =±1s1 =±1P (−1, −1) = (18)

s1 =±1



ここは 66 ページ目

大学院 混沌系工学特論 2004

最適化問題の数理 part 1 (担当 : 井上 純一)が得られるので, 一般に状態 s = (s1 , s2 ) をとる確率は

P (s)e(Js1 s2 +h1 s1 +h2 s2 )/Te−E(s)/T=(Js1 s2 +h1 s1 +h2 s2 )/TZ
s1 =±1s2 =±1 e= P (s1 , s2 ) = と書ける. ここで, Z は確率の規格化因子に相当し,

Z =



e−E(s)/T

(22)

(23)

s1 =±1 s2 =±1と書き直すことができ, 統計力学ではこの Z のことを分配関数と呼んでいる. また, 分布 (22) はボルツマ
ン分布, あるいは, ギブス分布と呼ばれる. 従って, [ノイズを利用したアルゴリズム] の定常分布はボルツ
マン分布になることがわかった. なお, 図 2 にはこの定常分布から見積もられる各確率を横線で書き込んであるが, この図から [ノイズを利用したアルゴリズム] から求められる定常確率はいずれもこのボルツマン分布に収束している様子が見て取れる.

2.1.3

最適化問題 : E(s1 , s2 ) = −J s1 s2 − h1 s1 − h2 s2 再考

さて, 我々の当初の問題に戻ろう. 我々の問題はエネルギー関数 (5) 式の最小値とその最小値を与える状
態 (s1 , s2 ) を求めることであった. そこで, [ノイズを利用したアルゴリズム] から得られる定常分布 (22) 式の形を見てみると, この式でノイズレベル T を十分に小さく設定すると, 低いエネルギーを与える状態が比較的高頻度で出現し, T = 0 では最小エネルギーをとる状態 (s1 , s2 ) = (1, 1) のみが確率 1 で現れることがわかるであろう. これは我々の期待していた結果である. つまり, P1 = 1, P2 = P3 = P4 = 0 という結果
が [ノイズを利用したアルゴリズム] から得られるのであれば, アンサンブルの中の全てのシステムの状態
は (s1 , s2 ) = (1, 1) という「正解」となっているはずである. そこで, ノイズレベルを T = 0.8, 0.6, 0.4, 0.2で固定した場合の計算機シミュレーション結果を見てみよう. その結果を図 4 に載せた. この図からまず,

T = 0.8 ではエネルギーの 1 番高い P (du), 及び, 2 番目に高い P (ud) が比較的小さな値をとり, 一番エネ
ルギーの低い P (uu) が最も大きく, 次にエネルギーの低い P (dd) が 2 番目に高い値に収束している. これは図 1 に示したエネルギー準位 ε1 < ε4 < ε2 < ε3 がきちんと反映されている. そして T = 0.4 ではエネル
ギーの高い P (du), P (ud) が完全にゼロになり, 2 番目にエネルギーの低い P (dd) が減少に転じ, P (uu) が増加し始める. そして, T = 0.2 に固定して [ノイズを利用したアルゴリズム] を走らせると t = 500 程度ではまだ定常分布に収束していないとは言え, この「好ましい傾向」は依然として見られる. そこで, 調子に乗って一挙に T = 0.01 までノイズレベルを下げて, このノイズレベルに固定してアルゴリズムを走らせて見よう. 結果を図 5 に載せる.
残念ながらこの図は我々の期待 (希望的観測) を裏切っている. P (uu), P (dd) とも完全に一定値に収束しているが, 我々が望む結果 : P (uu) = 1, P (dd) = 0 ではない. 何が問題だったのであろうか ?
そこで, 図 1 を見ながらもう一度考えてみると, [ノイズを利用したアルゴリズム] では 1 ステップで s1
か s2 のいずれか一方のみの状態変化しか許されていないので, 状態 (−1, −1) から状態 (1, 1) へ確率が流入するためには, 必ず, 状態 (1, −1) か (−1, 1) を経由しなければならないことがわかる. 従って, 例えば,(−1, −1) (E = ε4 ) から, (−1, 1) (E = ε3 ) に移るためには ∆E = ε3 − ε4 = 2(J − h2 ) の「エネルギー障壁」を飛び越えなければならない. この計算機シミュレーションでは J = 1, h2 = 0.1 と選んでいるので,[ノイズを利用したアルゴリズム] より, このエネルギー障壁を飛び越える確率は T = 0.01 のときP ((−1, −1) → (−1, 1)) =e−∆E/T = e−1.8/0.01 = e−180(24)であり, この逆数でもって状態 (−1, −1) に留まる滞在時間を定めれば, アルゴリズムをスタートさせてから, おおよそ t = e180 までの間, この状態 (−1, −1) に居座ることになる. これは事実上, (−1, −1) に状態が落ち込めば, そこからは永遠に脱出できないことを意味する. 従って, 初期状態の選び方が悪く, 初期状
ここは 67 ページ目

大学院 混沌系工学特論 2004

最適化問題の数理 part 1 (担当 : 井上 純一)

0.8

0.9
P(uu)

0.7

P(uu)

0.8
0.7

0.6

T=0.8

0.5

0.6

T=0.6

0.5
0.4

0.4
0.3

0.3
P(dd)

0.2

0.2
P(dd)

0.1

0.1

0
0

50

100

150

200

250

300

350

400

450

500

0
0

50

100

150

200

250

t

300

350

400

450

500

300

350

400

450

500

t

1

0.9

0.9

P(uu)

0.8
P(uu)

0.8

0.7

0.7
0.6

T=0.4

T=0.2

0.6
0.5
0.5
0.4

0.4

0.3

0.3

P(dd)

0.2

0.2
P(dd)

0.1

0.1

0
0

50

100

150

200

250

300

350

400

t

450

500

0
0

50

100

150

200

250

t

図 4: 左上から右下へ, T = 0.8, 0.6, 0.4, 0.2 とノイズレベルを固定した場合の各確率の時間変化の様子.
態 (−1, −1) を比較的高頻度でアンサンブルの中の各システムに与えてしまえば P (uu) は低くなってしま
い, 逆にアンサンブルの中の全てのシステムの初期状態を (1, 1) に選べば, P (uu) = 1 となることが予想される. 定常分布が初期状態に依存するという意味において, 実質的にはエルゴード性が破れてしまっているわけである3 . 実際にこれを見るために, T = 0.01 のときに複数の初期状態を選んだ場合の計算機シミュ
レーションによる P (uu) の時間発展を図 6 に載せよう. この図から予想通り, P (uu) の値は初期値の選び方に依存している.

2.1.4

ノイズレベルのスケジューリング

以上の結果から, いきなり T = 0 と置いてしまうと [ノイズを利用したアルゴリズム] は「ランダムシューティング」と何ら変わりなく, また, ゼロではないにしても T を低く設定しすぎれば (−1, −1) から脱出するのは我慢できないほどの時間がかかってしまうことがわかった. では T を時間的に変化させ,アルゴリズムのスタート段階では T を大きめにとって全ての状態を等確率で生成させるようにし, ある程度のアンサンブルが状態 (1, 1) に落ち込んだのを見はからって T を徐々に下げていくという戦略はどうであろうか ? これはうまく行きそうである. なぜならば, T = ∞ では完全に等確率 1/4 で 4 つの状態

(1, 1), (1, −1), (−1, 1), (−1, −1) が実現するが, 図 2 で見たようにノイズレベル T を T = 1 程度まで下げて

3 もちろん, ノイズレベル T が完全にゼロでは無いわけであるから, 本当はエルゴード性が破れているわけではなく, 無限の時間が経過しさえすれば, いずれは (熱) 平衡状態に達する. しかし, 我々の我慢できる時間の範囲内ではあたかもエルゴード性が破れているように見え, 平衡状態になることはとても期待できない.

ここは 68 ページ目

大学院 混沌系工学特論 2004

最適化問題の数理 part 1 (担当 : 井上 純一)

1

0.8
P(uu)

0.6

T=0.01
0.4
P(dd)

0.2

0
0

50

100

150

200

250

300

350

400

450

500

t

図 5: T = 0.01 にした場合の各確率時間発展. 横に走る破線はマスター方程式の定常解である平衡分布.
いくとエネルギーの低い順 : ε1 < ε4 < ε2 < ε3 に出現確率の大小関係 : P (uu) > P (dd) > P (ud) > P (du)が決まった. 従って, システムがあるノイズレベル T での (熱) 平衡状態に到達する限り, 得られる定常確率はエネルギー準位をそのまま反映するのである. さらに, 各ノイズレベルではアルゴリズムの初期値を再設定するわけではなく, 直前のノイズレベルでのアンサンブルの状態をそのまま使うわけであるから, いき
なりノイズレベルを T = 0 に設定したときのようなランダムシューティングにはならず, 必ずアルゴリズ
ムの全プロセスを通じて最小エネルギー状態が最も高頻度に現れるという「好ましい性質」「好ましい状態」を引きずりながら新しい状態が生成されていくことになる. よって, ノイズレベルを徐々に下げつつ,つまり, 各ノイズレベルではシステムが平衡状態になるのを待ってから次の (より小さな) ノイズレベルに設定する, というように [ノイズを利用したアルゴリズム] を修正すれば, 十分に長い時間をかけさえすれば確率 1 で最小エネルギー状態が求められる (つまり, P (uu) = 1, P (dd) = P (ud) = P (du) = 0) と推察でき
る4 . 前節 : [連想記憶の数理] ではノイズと言えば神経素子の誤動作の原因としてネガティブな要素であると扱われていたわけであるが, これをうまく用いてアルゴリズムを構成すれば好ましい情報処理が実現できる, というこの講義で見る最初の例がこれである.
さて, 早速これを確かめるため, T に時間依存性を持たせT (t) =
T (t) =

2
1 + 2t2log(1 + 2t)

(25)
(26)

のような 2 通りで下げた場合の結果を図 7 に載せる. この図より, T ∼ 1/t : (25) でノイズレベルを下げてしまうと, 下げ方が早すぎてうまくいっていない. 一方, T ∼ 1/ log t : (26) で下げると, P (uu) = 1 へ
4 確率 1 で最小状態を得るためにはアニーリングプロセスの最後のつめの段階 (T 0) で熱平衡状態を作る必要があり, 前述のように, これには極めて長い時間ががかる. 従って, 実用的にはどこか満足の行くところでアニーリングプロセスを打ち切ってしまう場合が多い. 近似解, それも精度の良い近似解なら OK とするわけである.

ここは 69 ページ目

大学院 混沌系工学特論 2004

最適化問題の数理 part 1 (担当 : 井上 純一)

P(uu)
1
0.9
0.8
0.7

T=0.01

0.6
0.5
0.4
0.3
0.2
0.1
0
0

50

100

150

200

250

300

350

400

450

500

t

図 6: T = 0.01 の場合に初期状態を様々変えてみたときの P (uu) の時間変化. 見かけ上, エルゴード性が破れている.
順調に向かっていることがわかる. このようにノイズレベル T に時間的な変化をもたせて徐々に下げていく [ノイズを利用したアルゴリズム] のことをシミュレーテッドアニーリングと呼んでおり, 様々な工学的場面で用いられている. 一般的に, エネルギー関数が極めて特殊な状態依存性を持たなければノイズレベル
を 1/ log t より十分に遅く下げていけば, 十分な時間の経過後には確率 1 で最小エネルギーを与える状態が得られることが知られている.

2.1.5

最適なノイズレベルのスケジューリング

ここではシミュレーテッドアニーリングにおける最適なノイズレベルのスケジューリングがT (t) ∼

1
log t

(27)

であることを議論しよう. 前にみたように, 状態 (1, 1) の方からも, (−1, −1) の方からもある有限のエネル
ギー障壁がある場合を我々は議論してきたわけであるから, 図 8 のような 2 準位系を考えよう. つまり, 図

8 に記された準位 1 が E = −J s1 s2 − h1 s1 − h2 s2 の場合の (−1, −1) に, 準位 2 が (1, 1) と見るわけである.
B はエネルギー準位 1 から 2 へ飛び移るのに越えなければならないエネルギー障壁の高さであり, ∆ は準
位 1, 2 間のエネルギー差である. B, ∆ とエネルギー E = −J s1 s2 − h1 s1 − h2 s2 との間にある対応関係は
B = ε3 − ε4 = 2(J − h2 ), B + ∆ = ε3 − ε1 = 2(J + h1 ) と考えればよい. 従って, 我々が計算機シミュレー
ションで確認した (26) 式で与えられるノイズレベルのスケジューリングがなぜうまく行ったのかを解析的にチェックするには図 8 のような 2 つの状態のいずれから他方へ行くのに有限のエネルギー障壁があるような 2 準位系を調べれば十分なのである.
さて, このシステムに対して [ノイズを利用したアルゴリズム] を適用すると, アンサンブルを構成する各システムの状態は 1,2 の間を行ったり来たりする.
ここは 70 ページ目

大学院 混沌系工学特論 2004

最適化問題の数理 part 1 (担当 : 井上 純一)

0.7

1
P(uu)

p(uu)

0.9

0.6
0.8

T=2/(1+2t)
0.5

0.7

T=2/log(1+2t)
0.6

0.4
P(dd)

0.5

0.3

0.4
0.3

0.2

0.2
0.1

P(dd)

0.1
0

0
0

50

100

150

200

250

300

350

400

450

500

0

50

100

150

200

250

t

300

350

400

450

500

t

図 7: ノイズレベルを T (t) = 2/(1 + 2t) (左) , T (t) = 2/ log(1 + 2t) (右) で下げた場合の各確率の時間変化の様子.
時刻 t に状態 1 をアンサンブル中の各システムがとる確率を P1 (t), 状態 2 をとる確率を P2 (t) としよう.
このとき, 確率 P1 の微小時間 ∆t の間の増加分 [P1 (t + ∆t) − P (t)] は次のように表すことができる.

P1 (t + ∆t) − P (t)= [P2 ×w(2→1) − P1 ×w(1→2)]∆t(28)

ここで, w(2→1), w(1→2) 等はそれぞれ状態 2 から 1 へ飛び移る遷移確率, 状態 1 から 2 へ飛び移る遷移確率である. ここで, [ノイズを利用したアルゴリズム] に従えば

w(2→1) =e−(B+∆)/T(29)

w(1→2) =

−B/T

(30)

e

でそれぞれの遷移確率が与えられることになる. また, ∆t → 0 の極限のもとで (28) 式の左辺はP1 (t + ∆t) − P1 (t) =dP1
dt∆t + O((∆t)2 )(31)となる. また, 確率の保存の式から P1 = 1 − P2 ≡ P とすると, (31) 式からは結局

dP= −[e−(B+∆)/T + e−B/T ]P + e−(B+∆)/T
dt

(32)

が得られる. 我々は既にみたマスター方程式を導出することができたわけである.
さて, 今の場合には状態 2 をアンサンブル中のほとんど全てのシステムがとるように確率 P = P1 ができるだけ小さくあってくれれば「最適化問題」の立場からは好ましいわけであるから, (32) 式の右辺をノ
イズレベル T の関数

L(T ) ≡−[e−(B+∆)/T + e−B/T ]P + e−(B+∆)/T(33)とみなし, 各時刻 t で L(T ) を最小化する T の条件を求めることにしよう (そうすれば結果的に dP/dtを最大化, つまり P の増加率を最大化することになり, 我々の目標にかなう). そのような条件は直ちに(∂L/∂T ) = 0 であることがわかるので, 実際この条件を書き下してみると−+(B + ∆) −(B+∆)/TBe+ 2 e−B/T P − [e−(B+∆)/T + e−B/T ]T2
T(B + ∆) −(B+∆)/T
e
=0
T2

∂P
∂T
(34)
ここは 71 ページ目

大学院 混沌系工学特論 2004

最適化問題の数理 part 1 (担当 : 井上 純一)

B

∆

1

2

図 8: 2 準位系.
ここで, 状態 1,2 間のエネルギー差 ∆ は状態 1,2 間のエネルギー障壁の高さと比べて小さいと仮定すると,(34) 式は次のように書き直すことができる.
dP
dT

=

−

e−∆/T
B
B
P+ 2
2
T
T
1 + e−∆/T− [e−(B+∆)/T + e−B/T ]P + e−(B+∆)/T(35)

上の微分方程式を解く前に, まずは次の斉次方程式 :

dP
dt

= −

B
P
T2

(36)

をまず解くことにしよう. これは変数分離形だから難なく解けて P (T ) = C eB/T が得られる. 次にこの積分定数を T の関数とみなして (定数変化法), P (T ) = C(T ) eB/T とし, これを (35) 式に代入してみると,

dP
dT

=

B
B
e−∆/T
dC(T ) B
B
B
B
e T − 2 C(T ) e T = − 2 C(T ) e T + 2
dT
TTT
1 + e−∆/T(37)となり, 整理すると

B 2 e−(B+∆)/T
dC(T )= 2dT
T
1 + e−∆/T(38)となる. 次に 1/T = x と変数変換すると,

dC
dx

= −

e−(B+∆)x
1 + e−∆x(39)

e−(B+∆)xdx
1 + e−∆x(40)と簡略化され, 形式的に (39) 式の解は

C(x) =−

ここは 72 ページ目

大学院 混沌系工学特論 2004

最適化問題の数理 part 1 (担当 : 井上 純一)と求めることができる. そこで e−∆x = y と置けば e−∆x dx = −dy/∆ であり, e−Bx = (e−∆x )B/∆ = y B/∆となるので (40) 式は

1∆

C(x) =



B

1
y∆
dy 
1+y
∆



B

(41)

y ∆ dyと変形できる. これを積分し, y = e−∆x を用いて元に戻せばC(T ) =

(B+∆)
1
e− T
(B + ∆)

(42)

∆
1
e− T
(B + ∆)(43)

P (B + ∆)(44)が得られ, P (T ) は

P (T ) =となる. 従って T = Topt の場合

e

− T∆

=

opt

が成り立つ. また,

e
e

− TB

− T∆

=

opt

−(B+∆)/Topt

(e

=

opt

B

B

) ∆ = {P (B + ∆)} ∆{P (B + ∆)}

(45)

B
∆ +1

(46)

なので, これら (44), (45), (46) 式を (39) 式に代入して整理すれば

dP
dt

B

B

= −[P ∆ +1 + P ]P + (B + ∆)P ∆ +1(47)が得られる. P ≤ 1 であることに注意して主要項のみ残すと

dP
dT

B

−P ∆ +2=(48)が得られ, これを積分すれば

∆
B



P

−∆
B

∆

t− B

(49)

が得られる. この (49) 式を (44) 式に代入して,log(B + ∆) −

∆
log
B

∆
B

−

∆
logt
B

= −

∆
Topt

(50)

となるので, t→∞ の極限で, 最適なノイズレベルのスケジューリングは

Topt

=

1
B
∼
log t
log t

(51)

である. 従って, 我々が前に計算機シミュレーションを行った際に T = 2/ log(1 + 2t) と選んでうまくいったという結果 (図 7 の右側) が解析的に裏付けられたことになる.

2.2

平衡状態と物理量の期待値 : 内部エネルギー, 自由エネルギー, エントロピー etc.

ここまでの議論で [ノイズを利用したアルゴリズム] を動作させることにより, アンサンブルの中の各システムが状態 s = (s1 , s2 ) を割合 :P (s1 , s2 ) =e(Js1 s2 +h1 s1 +h2 s2 )/Te−E(s)/T=(Js1 s2 +h1 s1 +h2 s2 )/T−E(s)/Ts∈±1 es1 =±1s2 =±1 e(52)

ここは 73 ページ目

大学院 混沌系工学特論 2004

最適化問題の数理 part 1 (担当 : 井上 純一)でとることがわかった. 時間が十分に経過し, システムが平衡状態に達すれば, 注目する 1 つのシステムの物理量の長時間平均は分布 P (s1 , s2 ) に関する平均で置き換えることができる.
ところで, 我々はエネルギー :E(s) =−J s1 s2 − h1 s1 − h2 s2(53)がいったい何に関する物理的対象のモデルを表すのか, についてはここまで明確にしてこなかったのではあるが, これは例えば, 2 つの磁石が互いに大きさ J の相互作用を及ぼしあっており, 個々の磁石にはそれぞれ, 上向きに (N 極方向に) h1 , h2 の強さの磁場がかかっている系として見ることができる. あるいは, 2 つの神経素子からなる回路が 1 つのパターンを Hebb 則により記憶しており, 個々の神経素子にはそのパター
ン方向へ素子を向けさせるようなバイアスが働いている神経系, つまり, 神経系のエネルギー関数が−J ξ11 ξ21 s1 s2 − h1 ξ11 s1 − h2 ξ21 s2

E(s) =(54)

で与えられ, 記憶したパターンを ξ 1 = (ξ11 , ξ21 ) = (1, 1) とした場合であると考えることもできる.
さて, このシステムのエネルギー E のアンサンブルでの期待値は β ≡ T −1 と置けば
 E =E(s1 , s2 )P (s1 , s2 )s1 =±1 s2 =±1



= −

Z

s1 =±1β(Js1 s2 +h1 s1 +h2 s2 )
s2 =±1 (Js1 s2 + h1 s1 + h2 s2 )e
β(Js1 s2 +h1 s1 +h2 s2 )s1 =±1s2 =±1 e

∂
= −
log Z
∂β
 
≡
eβ(Js1 s2 +h1 s1 +h2 s2 )(55)

s1 =±1 s2 =±1= eβ(J+h1 +h2 ) + eβ(−J+h1 −h2 ) + eβ(−J−h1 +h2 ) + eβ(J−h1 −h2 )(56)と表すことができる. 従って, エネルギーの期待値を求めるためには, 分配関数 Z, あるいは自由エネルギーと呼ばれる次の関数 :F=−T log Z

(57)

E

=

∂
(βF )
∂β

(58)

から

という操作を通じて求めることができる.
この自由エネルギーはとても便利な関数である. なぜならば, 任意の物理量 A(s1 , s2 ) の期待値 A(s1 , s2 )が算出したければ, A(s1 , s2 ) に共役な外場 (パラメータ) a を用いて aA(s1 , s2 ) を分配関数の指数部分に加えた自由エネルギー :F= −T log






e

β(Js1 s2 +h1 s1 +h2 s2 )+aA(s1 ,s2 )(59)

s1 =±1 s2 =±1を作り, これからA(s1 , s2 )1
= − lim
T a→0

∂F
∂a




=

β(Js1 s2 +h1 s1 +h2 s2 )
s2 =±1 A(s1 , s2 ) eβ(Js1 s2 +h1 s1 +h2 s2 )s1 =±1s2 =±1 es1 =±1(60)という操作を実行すればよい. この自由エネルギーはここで示したような計算技法上だけではなく, 前節[連想記憶の数理] で触れたような相転移現象に関する重要な情報を与えてくれる. それに関しては追々見ていくことになるなるであろう.
ここは 74 ページ目

大学院 混沌系工学特論 2004

最適化問題の数理 part 1 (担当 : 井上 純一)

また, このシステムのエントロピー S は熱力学的関係式 :=F

E − TS(61)

から, 自由エネルギー, 及び, 内部エネルギー E を用いてS=

1
( E − F)T

(62)

からみ積もることができる. 図 9 に内部エネルギー, 自由エネルギー, 及び, エントロピーのノイズレベル
依存性をプロットする. この図より, 例えば, T = 0 のノイズゼロの場合には内部エネルギー, つまり, エネ
1.5
1

S

0.5
0
-0.5

<E>
-1
-1.5

F

-2
-2.5
-3
0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

T

図 9: 我々が考えてきたシステム : E = −Js1 s2 − h1 s1 − h2 s2 の内部エネルギー E, 自由エネルギー F , 及び, エントロピー Sのノイズレベル依存性.

ルギーの期待値は ε1 = −J − h1 − h2 = −1.0 − 0.5 − 0.1 = −1.6 であるエネルギーの最小状態に一致する.
ノイズレベルが有限であるのであれば, ε2 , ε3 , ε4 のエネルギーをとるシステムも有限の割合で存在するので, その期待値は ε1 からずれた値をとるが, 完全にノイズレベルがゼロであるのであれば, アンサンブルを構成する全てのシステムが (s1 , s2 ) = (1, 1) の状態をとり, 従って, アンサンブルの「期待値」といっても,それは完全に ε1 に一致するわけである.
一方, エントロピーは全ての可能な状態が (1, 1) の一つだけであることから, ノイズ T = 0 では状態数の対数として log 1 = 0 となっていることがわかる. また, 十分にノイズが大きければ, 4 つの状態の全てが等確率で現れることになるので, T → ∞ で log 4 に漸近していくことが見て取れる.
ここまでの準備ができた段階であらためて我々が調べたかった問題 :


N




aj sj E(s) = 



(63)

j=1

を見ていくことにしょう (次回).
ここは 75 ページ目

大学院 混沌系工学特論 2004

最適化問題の数理 part 1 (担当 : 井上 純一)✬✩

問 8 : [ノイズを利用したアルゴリズム] を用いて以下の問いに答えよ.

(1) [ノイズを利用したアルゴリズム] をノイズレベル T で動作させることにより, 平衡状態を作り,得られるアンサンブルを用いてエネルギーの期待値 :アンサンブルのシステム総数（N )EiiE =アンサンブルのシステム総数 (N )を求めよ. また, これを図 9 の E (つまり, 手計算からの結果) と比較せよ.

(2) s1 , s2 の平衡分布での期待値 s1 , s2 を自由エネルギー経由で手計算により求めよ.
(3) [ノイズを利用したアルゴリズム] でノイズレベル T の平衡状態を作り, 得られるアンサンブルよりアンサンブルのシステム総数 (N )

s1

=

s2

=

[s1 ]ii
アンサンブルのシステム総数 (N )
アンサンブルのシステム総数 (N )[s2 ]ii
アンサンブルのシステム総数 (N )を計算し, (2) で求めた s1 , s2 と比較せよ.

(4) [ノイズを利用したアルゴリズム] を用いて, アンサンブル中の 1 つのシステムに注目した場合のエネルギー E の長時間平均 :
T

E

=

∞
1 
[Et (s1 , s2 )]T →T∞ Tt=0

lim

(T∞ は十分に大きな数) を求め, (1) で得られた E の値と比較せよ. E = E が成り立っているであろうか ?
✫✪ここは 76 ページ目