焼きなまし法のコツ Ver. 1.2

↓改訂版 Ver1.3はこちらから↓
shindannin.hatenadiary.com

2016年8月18日 Ver. 1.2 内容を更新しました。
4.遷移確率と温度の部分を全般的に修正

この記事は、Competitive Programming Advent Calendar Div2012の24日目の記事です。
http://partake.in/events/3fcea6d7-0bab-4597-82db-86439aadb1b9
素晴らしい企画をしていただいたtanzakuさん(https://twitter.com/_tanzaku_)、今年もどうもありがとうございます!

焼きなまし法そのものの解説は少しだけにして、焼きなまし法の使い方のコツをメインに書こうと思います。

焼きなまし法の概要

2015年5月26日 内容を更新しました。
最適化(=最大値か最小値を求める)の手法の1つです。以下は、てきとーな擬似コードです。

あらかじめ、計算する時間をTを決める。
以下のスコープ内のコードを時間T経過するまで繰り返す。
{
    適当に、次の状態(近傍)を選ぶ。
    (1) 次の状態のスコアが、今の状態のスコアより良くなった場合は、常に、状態を次の状態にする(遷移)。
    (2) 次の状態のスコアが、今の状態のスコアより悪くなった場合でも、ある確率で、状態を次の状態にする。
        確率は以下のように決める。
       (2.A) 序盤(時刻t=0に近いとき)ほど高確率で、終盤(t=Tに近いとき)ほど低確率。
       (2.B) スコアが悪くなった量が大きければ大きいほど、低確率。
}
最後に、今まで一番スコアのよかった状態を選ぶ。

次の状態が悪くてもそれを選ぶ可能性がある(2)の部分が、焼きなまし法の最大の特徴です。

「状態」「スコア」というのがピンとこないかもしれませんが、

  • もし、巡回セールスマン問題であれば、
    • 状態 →ルート
    • 次の状態(近傍)→元のルートをちょっとだけかえてみたルート
    • スコア →総距離
  • もし、y=f(x)の最大値を求める問題であれば、
    • 状態 →x
    • 次の状態(近傍)→ x' = x + a (aは小さめの値)
    • スコア →y

といったかんじです。

一応、参考までに、Marathon Match 60で使ったコードの断片を貼っておきます。全然難しくないです。

2015年4月20日 追加
以下のコード例は、上記の(2.B)の部分を無くした、スコアが悪くなった量は考慮しない焼きなまし法です。もっと効果的な方法「4. 遷移確率と温度」は必ず見てください。

vector <int> PolymerPacking::mirrorSequenceSI(ll msec)
{
    
    vector <int> currentSeq; // 現在の状態。この問題では、vector <int>が状態だけど、問題によって大きく異なる。

    // 現在の状態の初期化。この部分は問題によっていろいろなので、あまり気にしないでください。
    const int L = 5;
    const int MAXMIR = 5;
    for (int i = 0; i < MAXMIR; i++)
    {
        currentSeq.push_back(0);
    }

    vector <int>    bestSeq         = currentSeq;           // ベストの状態
    double          bestScore       = getScore(currentSeq); // ベストスコア。getScoreは、状態からスコアを求める関数
    double          currentScore    = bestScore;            // 現在のスコア
    {
        const ll saTimeStart    = getTime();            // 焼きなまし開始時刻。getTimeは、時間を返す関数
        const ll saTimeEnd      = m_startTime+msec;     // 焼きなまし終了時刻。m_startTimeはこのプログラム自体の開始時間
        ll saTimeCurrent        = saTimeStart;          // 現在の時刻

        while ((saTimeCurrent = getTime()) < saTimeEnd) // 焼きなまし終了時刻までループ
        {
            vector <int> nextSeq(currentSeq);       // 次の状態
            nextSeq[rand()%MAXMIR] = rand()%(L-1);  // 次の状態を求める。ランダムで1要素選んで、変えてみる。

            const double nextScore = getScore(nextSeq);

            // 最初t=0のときは、スコアが良くなろうが悪くなろうが、常にnextを使用
            // 最後t=Tのときは、スコアが改善したときのみ、nextを使用
            const ll T = saTimeEnd-saTimeStart;         // 焼きなましにかける時間
            const ll t = saTimeCurrent-saTimeStart;     // 焼きなまし法開始からの時間
            const ll R = 10000;
            // スコアが悪くなったときでも、次の状態に移動する遷移する場合はtrue。1次関数を使用。
            const bool FORCE_NEXT = R*(T-t)>T*(rand()%R);

            // スコアが良くなった or 悪くなっても強制遷移
            if (nextScore>currentScore || (FORCE_NEXT && currentScore>0.0f) )
            {
                currentScore = nextScore;
                currentSeq = nextSeq;
                printf("current Score=%.8f time=%lld\n",currentScore,t);
            }

            // ベストスコアは保存しておく。
            if(currentScore>bestScore)
            {
                bestScore = currentScore;
                printf("New best Score=%.8f time=%lld\n",bestScore,t);
                bestSeq = currentSeq;
            }
        }
    }

    // ベストスコアをしたときの状態を返す。
    return bestSeq;
}

焼きなまし法のコツ

が多いところほど、自分では重要だと思ってます。自信がないので、「こういうやり方もあるよ」「それはない…」「嘘乙」など、アドバイス大歓迎です!

1. 高速化

たくさんwhileループ内を回せば回すほど、スコアが伸びる可能性があります。基本的にはボトルネックを調べて、そこを高速化しましょう。だいたい以下のところがボトルネックになるでしょう。

★★★★ 変更がある変数だけを保存して、戻す

先ほどのコードだと、

        vector <int> nextSeq(currentSeq);       // 次の状態
        nextSeq[rand()%MAXMIR] = rand()%(L-1);	// 次の状態を求める。ランダムで1要素選んで、変えてみる。


とやっていますが、これでは時間がかかります。。次の状態を全てコピーしたのち、1つの要素を変えています(しかもvectorで動的にメモリ確保もしてます…)。高速化するなら、

        const int index = rand()%MAXMIR; // 変更する要素
        const int backup_value = currentSeq[index]; // 変更前の値
        currentSeq[index] = rand()%(L-1);
        // 以下currentSeqをnextSeqとみなす。
        const double nextScore = getScore(currentSeq);

        // 最初t=0のときは、スコアが良くなろうが悪くなろうが、常にnextを使用
        // 最後t=Tのときは、スコアが改善したときのみ、nextを使用
        const ll T = saTimeEnd-saTimeStart;         // 焼きなましにかける時間
        const ll t = saTimeCurrent-saTimeStart;     // 焼きなまし法開始からの時間
        const ll R = 10000;
        // スコアが悪くなったときでも、次の状態に移動する遷移する場合はtrue確率。この場合は単に1次関数を使用。
        const bool FORCE_NEXT = R*(T-t)>T*(rand()%R);

        // スコアが良くなった or 悪くなっても強制遷移
        if (nextScore>currentScore || (FORCE_NEXT && currentScore>0.0f) )
        {
            currentScore = nextScore;
            // currentSeq = nextSeq; ←すでにnextSeqになってる。いらない。
            printf("current Score=%.8f time=%lld\n",currentScore,t);
        }
        else
        {
            // アンドゥ。次の状態に移動しないので、バックアップしてた値に戻す
            currentSeq[index] = backup_value;
        }

とやりましょう。

★★ 逆の手で、戻す

逆の手をやれば、前の状態に高速に戻せるのであれば、そうしましょう。例えば、上下左右に駒をうごかせるというパズル問題であれば、左に動かした後、右に動かせば、戻せるかもしれません。

★★★★ スコア関数の高速化

スコア関数の求め方は問題によって大きく異なるので、計算時間も異なってきます。ボトルネックになってたら高速化しましょう。
また、スコア関数の正確さを犠牲にしてでも、高速化したほうが良い場合もあります。1つのコツとして覚えて置きましょう。

★★ C言語のrand()を使わない。

C言語のrand()は遅いので、高速なxor128を使いましょう。

unsigned int randxor()
{
    static unsigned int x=123456789,y=362436069,z=521288629,w=88675123;
    unsigned int t;
    t=(x^(x<<11));x=y;y=z;z=w; return( w=(w^(w>>19))^(t^(t>>8)) );
}
★★ 時間計測の関数を呼びすぎない。

上にあげたコードは、時間計測の関数を呼びすぎです。。getTime関数は高速ではなく、それがボトルネックになる可能性もあります。例えば、最初のコードなら、以下のように100回に1回だけ時間を測るというようにできます。

        while ((saTimeCurrent = getTime()) < saTimeEnd) // 焼きなまし終了時刻までループ
        {
            for(int loop=0; loop<100; ++loop)
            {

ただ、TopCoderだとループ回数を増やしすぎると、ループしている間に、制限時間オーバーしてしまうという可能性もあります。そういった点も考えて、適切な回数だけ時間をはかるようにしましょう。残り時間が少なくなったら、ループ回数を減らすといった安全策をとる方法もあります。



あと、高速化に夢中になる前に、

★★★★ 高速化がどれだけ効果的か検討する。

のも重要です。ためしに10秒だけやってる計算を、あえて100秒とかに伸ばしてみて、結果がどうなるか見てみましょう。

  • 時間をかければかけるほどスコアが伸びる
    • 高速化しよう&スコアが素早く伸びるように次の状態を決めよう。山登り法に変えるのも良い。
  • 時間をかけてもスコアが伸びない & スコアが高い
    • 最適解がでてる可能性が高い。計算を打ち切って、余った時間を他に使う
  • 時間をかけてもスコアが伸びない & スコアが低い
    • 次の状態の決め方、スコア関数の見直しを再検討する。焼きなまし法を使わない法がよいときもあるのも忘れずに

2. 次の状態(状態近傍)をどう決めるか?

ここは鍵になる部分ですが、コツとしてまとめるのは難しいところですね…。

★★★★ スコアが改善しやすくなるように、次の状態を選ぶ。


これは、問題ごとによく考えて対処しなければいけないので、あまり「コツ」というかんじではないのですが…。

[問題]
タスク数aの仕事Aと、タスク数bの仕事Bがある。作業員はx人いて、それぞれ、1時間あたりにこなせる仕事Aのタスク数と、仕事Bのタスク数は決まっている。ただし、仕事Aの仕事Bのどちらかしか行うことができない。仕事Aと仕事Bの両方が終わるのに時間が最小になるように、作業員を仕事Aか仕事Bのどちらかに割り当てよ。

焼きなまし法で解くなら、まず最初にてきとーに作業員をA,Bに割り振り、ここからスコア=作業時間として、改善させていきます。こういった問題を解くときに、次の状態として以下の2つを考えてみましょう。

  1. 移動:ある作業員をA→B または B→A に移動させる。
  2. 交換:Aの作業員の誰かと、Bの作業員の誰かを交換する。

移動2回は交換になることもあるので、とりあえず移動をさせておけば、交換もたまにしてくれるし良さそうに見えますが、そうとは限りません。この場合は、交換を主に行い、たまに移動させるぐらいがちょうど良いと思います。移動では、人数が偏り、片方の仕事の作業時間が増えやすくなります。2回移動すれば、交換になって、作業時間が減る可能性がありますが、

  • 移動:良スコア → 悪いスコア → とても良いスコア
  • 交換:良スコア    → とても良いスコア

と考えると、移動だけだと、スコアが悪い状態を経由しないとダメなので、運がよくないと、とても良いスコアまで辿りつけません。その点、交換ならば、人員が偏った悪いスコアになりやすい状態をスキップできるので、スムーズにとても良いスコアまで改善できます。

もちろん、正しい人数割り当てもしなければいけないので、移動も少しは試したほうがいいですが、こういった点も考えて次の状態を決めると、スコアが改善しやすくなります。

★★ Greedyに、次の状態を選ぶ。

先ほどの問題であれば、「仕事Aと仕事Bそれぞれから、常に作業時間が一番多くかかっている人を1人ずつ選び、交換」とGreedyに、次の状態を選ぶことができます。こういう選び方をすると、ベストスコアは序盤はどんどん伸びていきますが、局所最適に陥りやすいです。また、ありがちなのが、Greedyに選ぶアルゴリズムが、計算時間のボトルネックになってしまうパターンです。複雑なアルゴリズムより、ランダムでとにかくいろいろ試すのが焼きなまし法の強みです。Greedyにやるにせよ、多少ランダムを入れたり、計算時間はかからない単純な選び方をしたりする工夫が必要です。

★★ 序盤は大きく変化、終盤は小さく変化させたほうがいい。

山登り法と同じ考えです。
単純にy=f(x)のyを最大化する問題なら、

  • 序盤(t=0に近いとき) 次のx'は x'=x+rand()%1001 - 500
  • 終盤(t=Tに近いとき) 次のx'は x'=x+rand()% 11 - 5

とか、
作業員割り当ての問題なら、

  • 序盤(t=0に近いとき) 交換を100回やる
  • 終盤(t=Tに近いとき) 交換を 1回やる

といったかんじです。スコアが局所最適でハマっている場合は有効です。ただ、焼きなまし法自体が、局所最適にはまりにくいアルゴリズムなので、あまり効果がない場合も多いです。

★★ 良い初期状態からスタートする

2015年5月26日 追加しました。
これは問題に大きく依存する項目です。初期状態はてきとーでも良い場合と、すこし時間をかけて計算して良い状態から始めたほうが良い場合があります。特に、初期状態が悪いと状態を改善する近傍の選択肢が少なすぎる、逆に言えば、状態を改善すればするほどさらに改善できる可能性が増えるタイプの問題は、初期状態に少し時間を書けて求めて、それから焼きなましたほうが良い結果がでる場合があります。

3. スコア関数

スコア関数(scoring function)とも、評価関数(evaluation function)とも呼びます。

★★★★★ 良い状態を正しく評価できるような、スコア関数を新たに用意する

問題のスコア関数をそのまま使うと、焼きなまし法がうまく行かないことがあります。良い状態を正しく評価できるような、新たなスコア関数を用意しましょう。

[問題]
 (大きく略)…  i番目の人の満足値をs(i)とする。「満足度の最小値 min s(i) 」を、最大化しなさい。

この問題は「満足度の最小値 min s(i) 」を最大化する問題なので、素直に考えると、スコア=「満足度の最小値 min s(i) 」としたくなりそうですが、これは良くありません。仮にそうした場合、

  状態A s(0)=  10, s(1)=10, s(2)=  10 のとき、スコアはmin{s(0),s(1),s(2)} = 10
  状態B s(0)=1000, s(1)=10, s(2)=1000 のとき、スコアはmin{s(0),s(1),s(2)} = 10

と状態Aと状態Bのスコアは同じになってしまいます。これはスコアは同じといえ、状態Bのほうが明らかに良い状態です。なぜなら、s(1)の点が仮に300点に増えたとしたら

  状態A' s(0)=  10, s(1)=300, s(2)=  10 のとき、スコアはmin{s(0),s(1),s(2)} =  10
  状態B' s(0)=1000, s(1)=300, s(2)=1000 のとき、スコアはmin{s(0),s(1),s(2)} = 300

と状態Bのほうは一気にスコアが増えます。しかし元のスコア関数のままだと状態Aと状態Bは同じスコアで、状態A→状態Bはスコア改善になりません。つまり焼きなまし法を行っても、状態A→状態Bには運が良くないと移動しません。
この問題であれば、例えば、
スコア=1番小さいs(i) + 0.0001 * 2番目に小さいs(i)
とすれば

  状態A s(0)=  10, s(1)=10, s(2)=  10 のとき、スコアは 10 + 0.0001 *   10  = 10.001
  状態B s(0)=1000, s(1)=10, s(2)=1000 のとき、スコアは 10 + 0.0001 * 1000  = 10.1

と状態Bの方が好スコアになり、状態A→状態Bに移動しやすくなります。
ここは問題によって大きく有効かどうか大きく異なってきますが、必ず考慮に入れたほうがいいでしょう。結果に大差がつきやすいポイントです。

★★★スコアが改善しやすくなるようなスコア関数を決める。

2015年5月26日 追加しました
上記の言い換えともいえるのですが、わりと焼きなまし法のスコア関数をつくるときに、問題のスコア関数をそのまま使ったり、違う状態を正しく評価しなかったりすると、同じスコアになりがちで、青い線のような飛び飛びの平坦な関数になることが多いです。近傍に同じ値が続くと、焼きなまし法をやっても状態が行ったりきたりするだけで、解がなかなか改善しません。赤いように、山登りしやすいスコア関数を作れないか考えましょう。必ずしも、問題のスコア関数の最大値と、焼きなまし法で使うスコア関数の最大値が一致する必要はありません。最大値をとりうるような状態(argmax)を一致させることが重要です。
f:id:shindannin:20150526014955p:plain

4. 遷移確率と温度

たいていの焼きなまし法の解説書には、expを使った方法が載っています。
焼きなまし法 - Wikipedia

遷移確率 P(e, e' , T) は e' < e の時に 1 を返すよう定義されている(すなわち下りの移動は必ず実行される)。
そうでない場合の確率はexp((e−e' )/T) と定義 

注:wikipediaの引用なので、本文の設定と違います。
    ここでのTは温度(=時間とともに減っていく値)です。
    最小化問題なのでe−e'となってます(最大化問題ならe'-eです)。
時刻tに置ける温度Tを調整する

最初の例の以下のコードの部分です。

     const bool FORCE_NEXT = R*(T-t)>T*(rand()%R);

上の例だと、

  • t=0のとき、常にtrue
  • t=Tのとき、常にfalse

になるような線形関数にしています。だいたいの場合は、これで問題ありませんが、バリエーションはいろいろ考えられます。実際には、ベストスコアの更新状況をみて、更新があまりに序盤と終盤に固まっているときは、2乗してみるとか、ロジスティック曲線をつかうとか、その程度で良いと思います。

★★★★exp(e-e')/Tを使う。(e-e':スコアの悪化量)

2015年5月26日 重要度を★★★から上げました。全般的に内容を更新しました。
wikipediaでは以下のように書いていますが、exp(e-e')/Tの形にしたほうが良いでしょう。(理由は、次回の改定時に書きます)

焼きなまし法 - Wikipedia

関数 P は e' −e の値が増大する際には確率を減らす値を返すように設定される。
つまり、ちょっとしたエネルギー上昇の向こうに極小がある可能性の方が、
どんどん上昇している場合よりも高いという考え方である。しかし、この条件は
必ずしも必須ではなく、上記の必須条件が満たされていればよい。

ただ、このexpの式、分子はスコア差なので分かりますが、分母の温度をどう決めれば良いかは分かりづらいです。目安としては、たまには許容しても良いスコアの悪化量を温度とすると良いです。 最大の山を越えられるようにと考えると、考えられるスコアの最大値-最小値を初期温度を設定したほうが良いと思います。(理由は、次回の改定時に書きます)

double startTemp   =  100;
double endTemp     =   10;
double temp        = startTemp + (endTemp - startTemp) * t / T;
double probability = exp((nextScore - currentScore) / temp);  // 最大化なので、分子はnextScore - currentScoreで良い。
bool FORCE_NEXT    = probability > (double)(mXor.random() % R) / R;

時刻t=0(最初)でtemp=startTemp、時刻t=T(最後)でtemp=endTempとなります。

例えばスコアが100悪化したとき、nextScore - currentScore = -100なので、startTempを100にしてみます。endTempはとりあえず10にしてみましょう。
時間t=0であれば、exp(-100/100) = exp(- 1) = 0.36787944117と、3回に1回ぐらいは遷移します。
時間t=Tであれば、exp(-100/ 10) = exp(-10) = 0.00004539992と、まず遷移はしません。
このような感じでやると、良いと思います。

5. アレンジ

★★★★ 山登り

最初の例をFORCE_NEXT=0にすれば、山登り法になります。そもそも関数が凸であれば、山登り法も焼きなまし法も結果は一緒で、しかもrand()を使わない分、山登り法のほうが高速で、スコア改善のチャンスも多くなります。簡単に試せるので、試しましょう。

多点スタート

2015年5月26日 重要度を★★から下げました。内容を更新しました。
焼きなまし法をつかっても、局所解にはまりやすい場合は、初期値を毎回かえた、多点スタートにする方法もあります。山登りと同時に使う方法もあります。ただ、多点スタートで結果がよくなるケースは、そもそもスコア関数や近傍が良くないので、そこを見直したほうがよいかもしれません。

★★★ 1変数ごとに焼きなまし

焼きなまし法で、たとえば、y = f(x1,x2,x3,x4,x5)と、状態変数が多いとき、次の状態のパターンが多くなり、結果的にあまり多くの状態を試せず、良い解が求められないときがあります。このような場合は、x2=x3=x4=x5を固定して、x1を焼きなまし、x1=x3=x4=x5を固定して、x2を焼きなまし、x1=x2=x4=x5を固定して、x3を焼きなましというふうに1変数ずつ焼きなまししていく方法があります。xの値がそれぞれ独立性が高い(?)場合は、この方法は有効です。 ←要検討…