秒懂神經網絡---真·模擬退火算法

著名的模擬退火算法，它是一種基於蒙特卡洛思想設計的近似求解最優化問題的方法。

 

一點歷史——如果你不感興趣，可以跳過

美國物理學家 N.Metropolis 和同仁在1953年發表研究復雜系統、計算其中能量分布的文章，他們使用蒙特卡羅模擬法計算多分子系統中分子的能量分布。這相當於是本文所探討之問題的開始，事實上，模擬退火中常常被提到的一個名詞就是Metropolis准則，后面我們還會介紹。

美國IBM公司物理學家 S.Kirkpatrick、C. D. Gelatt 和 M. P. Vecchi 於1983年在《Science》上發表了一篇頗具影響力的文章：《以模擬退火法進行最優化（Optimization by Simulated Annealing）》。他們借用了Metropolis等人的方法探討一種旋轉玻璃態系統（spin glass system）時，發覺其物理系統的能量和一些組合最優（combinatorial optimization）問題（著名的旅行推銷員問題TSP即是一個代表例子）的成本函數相當類似：尋求最低成本即似尋求最低能量。由此，他們發展出以 Metropolis 方法為本的一套算法，並用其來解決組合問題等的尋求最優解。

 

幾乎同時，歐洲物理學家 V.Carny 也發表了幾乎相同的成果，但兩者是各自獨立發現的；只是Carny“運氣不佳”，當時沒什么人注意到他的大作；或許可以說，《Science》雜志行銷全球，“曝光度”很高，素負盛名，而Carny卻在另外一本發行量很小的專門學術期刊《J.Opt.Theory Appl.》發表其成果因而並未引起應有的關注。

 

Kirkpatrick等人受到Metropolis等人用蒙特卡羅模擬的啟發而發明了“模擬退火”這個名詞，因為它和物體退火過程相類似。尋找問題的最優解（最值）即類似尋找系統的最低能量。因此系統降溫時，能量也逐漸下降，而同樣意義地，問題的解也“下降”到最值。

 

一、什么是退火——物理上的由來

 

在熱力學上，退火（annealing）現象指物體逐漸降溫的物理現象，溫度愈低，物體的能量狀態會低；夠低后，液體開始冷凝與結晶，在結晶狀態時，系統的能量狀態最低。大自然在緩慢降溫（亦即，退火）時，可“找到”最低能量狀態：結晶。但是，如果過程過急過快，快速降溫（亦稱「淬煉」，quenching）時，會導致不是最低能態的非晶形。

 

如下圖所示，首先（左圖）物體處於非晶體狀態。我們將固體加溫至充分高（中圖），再讓其徐徐冷卻，也就退火（右圖）。加溫時，固體內部粒子隨溫升變為無序狀，內能增大，而徐徐冷卻時粒子漸趨有序，在每個溫度都達到平衡態，最后在常溫時達到基態，內能減為最小（此時物體以晶體形態呈現）。

 

似乎，大自然知道慢工出細活：緩緩降溫，使得物體分子在每一溫度時，能夠有足夠時間找到安頓位置，則逐漸地，到最后可得到最低能態，系統最安穩。

 

二、模擬退火（Simulate Anneal）

 

如果你對退火的物理意義還是暈暈的，沒關系我們還有更為簡單的理解方式。想象一下如果我們現在有下面這樣一個函數，現在想求函數的（全局）最優解。如果采用Greedy策略，那么從A點開始試探，如果函數值繼續減少，那么試探過程就會繼續。而當到達點B時，顯然我們的探求過程就結束了（因為無論朝哪個方向努力，結果只會越來越大）。最終我們只能找打一個局部最后解B。

模擬退火其實也是一種Greedy算法，但是它的搜索過程引入了隨機因素。模擬退火算法以一定的概率來接受一個比當前解要差的解，因此有可能會跳出這個局部的最優解，達到全局的最優解。以上圖為例，模擬退火算法在搜索到局部最優解B后，會以一定的概率接受向右繼續移動。也許經過幾次這樣的不是局部最優的移動后會到達B 和C之間的峰點，於是就跳出了局部最小值B。

 

根據Metropolis准則，粒子在溫度T時趨於平衡的概率為exp(-ΔE/(kT))，其中E為溫度T時的內能，ΔE為其改變數,k為Boltzmann常數。Metropolis准則常表示為

Metropolis准則表明，在溫度為T時，出現能量差為dE的降溫的概率為P(dE)，表示為：P(dE) = exp( dE/(kT) )。其中k是一個常數，exp表示自然指數，且dE<0。所以P和T正相關。這條公式就表示：溫度越高，出現一次能量差為dE的降溫的概率就越大；溫度越低，則出現降溫的概率就越小。又由於dE總是小於0（因為退火的過程是溫度逐漸下降的過程），因此dE/kT < 0 ，所以P(dE)的函數取值范圍是(0,1) 。隨着溫度T的降低，P(dE)會逐漸降低。

我們將一次向較差解的移動看做一次溫度跳變過程，我們以概率P(dE)來接受這樣的移動。也就是說，在用固體退火模擬組合優化問題，將內能E模擬為目標函數值 f，溫度T演化成控制參數 t，即得到解組合優化問題的模擬退火演算法：由初始解 i 和控制參數初值 t 開始，對當前解重復“產生新解→計算目標函數差→接受或丟棄”的迭代，並逐步衰減 t 值，算法終止時的當前解即為所得近似最優解，這是基於蒙特卡羅迭代求解法的一種啟發式隨機搜索過程。退火過程由冷卻進度表(Cooling Schedule)控制，包括控制參數的初值 t 及其衰減因子Δt 、每個 t 值時的迭代次數L和停止條件S。

 

總結起來就是：

• 若f( Y(i+1) ) <= f( Y(i) ) (即移動后得到更優解)，則總是接受該移動；
• 若f( Y(i+1) ) > f( Y(i) ) (即移動后的解比當前解要差)，則以一定的概率接受移動，而且這個概率隨着時間推移逐漸降低（逐漸降低才能趨向穩定）相當於上圖中，從B移向BC之間的小波峰時，每次右移（即接受一個更糟糕值）的概率在逐漸降低。如果這個坡特別長，那么很有可能最終我們並不會翻過這個坡。如果它不太長，這很有可能會翻過它，這取決於衰減 t 值的設定。

 

關於普通Greedy算法與模擬退火，有一個有趣的比喻：

• 普通Greedy算法：兔子朝着比現在低的地方跳去。它找到了不遠處的最低的山谷。但是這座山谷不一定最低的。這就是普通Greedy算法，它不能保證局部最優值就是全局最優值。
• 模擬退火：兔子喝醉了。它隨機地跳了很長時間。這期間，它可能走向低處，也可能踏入平地。但是，它漸漸清醒了並朝最低的方向跳去。這就是模擬退火。

 

通過一個實例來編程演示模擬退火的執行。特別地，我們這里所采用的實例是著名的“旅行商問題”（TSP，Traveling Salesman Problem），它是哈密爾頓回路的一個實例化問題，也是最早被提出的NP問題之一。

TSP是一個最常被用來解釋模擬退火用法的問題，因為這個問題比較有名，我們這里不贅言重述，下面直接給出C++實現的代碼：

#include <iostream>
#include <string.h>
#include <stdlib.h>
#include <algorithm>
#include <stdio.h>
#include <time.h>
#include <math.h>

#define N     30      //城市數量
#define T     3000    //初始溫度
#define EPS   1e-8    //終止溫度
#define DELTA 0.98    //溫度衰減率

#define LIMIT 1000   //概率選擇上限
#define OLOOP 20    //外循環次數
#define ILOOP 100   //內循環次數

using namespace std;

//定義路線結構體
struct Path
{
    int citys[N];
    double len;
};

//定義城市點坐標
struct Point
{
    double x, y;
};

Path bestPath;        //記錄最優路徑
Point p[N];       //每個城市的坐標
double w[N][N];   //兩兩城市之間路徑長度
int nCase;        //測試次數

double dist(Point A, Point B)
{
    return sqrt((A.x - B.x) * (A.x - B.x) + (A.y - B.y) * (A.y - B.y));
}

void GetDist(Point p[], int n)
{
    for(int i = 0; i < n; i++)
        for(int j = i + 1; j < n; j++)
            w[i][j] = w[j][i] = dist(p[i], p[j]);
}

void Input(Point p[], int &n)
{
    scanf("%d", &n);
    for(int i = 0; i < n; i++)
        scanf("%lf %lf", &p[i].x, &p[i].y);
}

void Init(int n)
{
    nCase = 0;
    bestPath.len = 0;
    for(int i = 0; i < n; i++)
    {
        bestPath.citys[i] = i;
        if(i != n - 1)
        {
            printf("%d--->", i);
            bestPath.len += w[i][i + 1];
        }
        else
            printf("%d\n", i);
    }
    printf("\nInit path length is : %.3lf\n", bestPath.len);
    printf("-----------------------------------\n\n");
}

void Print(Path t, int n)
{
    printf("Path is : ");
    for(int i = 0; i < n; i++)
    {
        if(i != n - 1)
            printf("%d-->", t.citys[i]);
        else
            printf("%d\n", t.citys[i]);
    }
    printf("\nThe path length is : %.3lf\n", t.len);
    printf("-----------------------------------\n\n");
}

Path GetNext(Path p, int n)
{
    Path ans = p;
    int x = (int)(n * (rand() / (RAND_MAX + 1.0)));
    int y = (int)(n * (rand() / (RAND_MAX + 1.0)));
    while(x == y)
    {
        x = (int)(n * (rand() / (RAND_MAX + 1.0)));
        y = (int)(n * (rand() / (RAND_MAX + 1.0)));
    }
    swap(ans.citys[x], ans.citys[y]);
    ans.len = 0;
    for(int i = 0; i < n - 1; i++)
        ans.len += w[ans.citys[i]][ans.citys[i + 1]];
    cout << "nCase = " << nCase << endl;
    Print(ans, n);
    nCase++;
    return ans;
}

void SA(int n)
{
    double t = T;
    srand((unsigned)(time(NULL)));
    Path curPath = bestPath;
    Path newPath = bestPath;
    int P_L = 0;
    int P_F = 0;
    while(1)       //外循環，主要更新參數t，模擬退火過程
    {
        for(int i = 0; i < ILOOP; i++)    //內循環，尋找在一定溫度下的最優值
        {
            newPath = GetNext(curPath, n);
            double dE = newPath.len - curPath.len;
            if(dE < 0)   //如果找到更優值，直接更新
            {
                curPath = newPath;
                P_L = 0;
                P_F = 0;
            }
            else
            {
                double rd = rand() / (RAND_MAX + 1.0);
                //如果找到比當前更差的解，以一定概率接受該解，並且這個概率會越來越小
                if(exp(dE / t) > rd && exp(dE / t) < 1)
                    curPath = newPath;
                P_L++;
            }
            if(P_L > LIMIT)
            {
                P_F++;
                break;
            }
        }
        if(curPath.len < bestPath.len)
            bestPath = curPath;
        if(P_F > OLOOP || t < EPS)
            break;
        t *= DELTA;
    }
}

int main(int argc, const char * argv[]) {

    freopen("TSP.data", "r", stdin);
    int n;
    Input(p, n);
    GetDist(p, n);
    Init(n);
    SA(n);
    Print(bestPath, n);
    printf("Total test times is : %d\n", nCase);
    return 0;
}

注意由於是基於蒙特卡洛的方法，所以上面代碼每次得出的結果並不完全一致。你可以通過增加迭代的次數來獲得一個更優的結果。

我們這里需要說明的是，在之前的文章里，我們用求最小值的例子來解釋模擬退火的執行：如果新一輪的計算結果更前一輪之結果更小，那么我們就接受它，否則就以一個概率來拒絕或接受它，而這個拒絕的概率會隨着溫度的降低（也即是迭代次數的增加）而變大（也就是接受的概率會越來越小）。

但現在我們面對一個TSP問題，我們如何定義或者說如何獲取下一輪將要被考察的哈密爾頓路徑呢？在一元函數最小值的例子中，下一輪就是指向左或者向右移動一小段距離。而在TSP問題中，我們可以采用的方式其實是很多的。上面代碼中GetNext()函數所采用的方式是隨機交換兩個城市在路徑中的順序。例如當前路徑為 A->B->C->D->A，那么下一次路徑就可能是A->D->C->B->A，即交換B和D。而在文獻【3】中，作者采樣的代碼如下（我們截取一個片段，完整代碼請參考原文）：

public class Tour{ ... ... // Creates a random individual public void generateIndividual() { // Loop through all our destination cities and add them to our tour for (int cityIndex = 0; cityIndex < TourManager.numberOfCities(); cityIndex++) { setCity(cityIndex, TourManager.getCity(cityIndex)); } // Randomly reorder the tour Collections.shuffle(tour); } ... ... }

可見作者的方法是把上一輪路徑做一個隨機的重排（這顯然也是一種策略）。

TSP.data的數據格式如下，第一行的數字表示一個有多少座城市，第2至最后一行，每行有兩個數字表示，城市的坐標（平面直角坐標系）。例如：
6
20 80
16 84
23 66
62 90
11 9
35 28
最后請讀者自行編譯執行程序並觀察分析輸出的結果。

 

 

參考文獻與推薦閱讀材料

【1】關於哈密爾頓問題和TSP問題請參考下面兩個資料以了解更多：

• 從哈密爾頓路徑談NP問題
• William J. Cook，迷茫的旅行商：一個無處不在的計算機算法問題，人民郵電出版社，2013

【2】上面的C++代碼在下面這個兩個帖子中都有給出，原作者無法考證

• http://www.henufz.cn/bencandy.php?fid=151&id=1894 嚇尿了，這居然是高中的信息學科競賽！！！原來我們在高中就輸了....

• http://blog.csdn.net/acdreamers/article/details/10019849

【3】關於TSP問題的一個Java語言實現的源碼，請參考

• http://www.theprojectspot.com/tutorial-post/simulated-annealing-algorithm-for-beginners/6