模式识别：模拟退火算法的设计与实现

本节的目的是记录下面学习和掌握模拟退火(Simulated Annealing，简称SA算法)过程。

模拟退火算法是一种通用概率算法，用来在一个大的搜寻空间内寻找命题的最优解。这里分别使用随机模拟退火算法和确定性模拟退火算法。在MATLAB平台上进行编程。以寻找一个6-单元全连接网络的能量最小化模型。

參考书籍：Richard O. Duda, Peter E. Hart, David G. Stork 著《模式分类》

一、技术论述

1.随机方法

学习在构造模式分类器中起着中心的作用。一个通常的做法是先假设一个单參数或多參数的模型，然后根据训练样本来预计各參数的取值。

当模型相当简单而且低维时，能够採用解析的方法。如求函数导数，来显式求解方程以获得最优參数。

假设模型相对复杂一些，则能够通过计算局部的导数而採用梯度下降算法来解，如人工神经网络或其它一些最大似然方法。而对于更复杂的模型，经常会出现很多局部极值，上述方法的效果往往不尽人意。

假设一个问题越复杂，或者先验知识和训练样本越少，我们对能够自己主动搜索可行解的复杂搜索算法的依赖性就越强，如基于參数搜索的随即方法。

一个通常的做法是使搜索朝着预期最优解的区域前进，同一时候同意一定程度的随机扰动，以发现更好的解。

2.随机搜索

这里主要研究在多个候选解中搜索最优解的方法。假设给定多个变量si,i=1,2,…,N。当中每一个变量的数值都取两个离散值之中的一个（如-1和1）。

优化问题是这样描写叙述的：确定N个si的合适取值。时下述能量函数（又称为代价函数）最小：

当中w_ij是一个实对称矩阵，取值可正可负。当中令到自身的反馈权为零（即w_ii=0），这是由于非零w_ii仅仅是在能量函数E上添加一个与si无关的常数，不影响问题的本质。这个优化问题能够用网络和节点的方式表示。例如以下图所看到的。当中节点之间的链接相应每一个权值w_ij。

3.贪心算法的局限性

如上所述。对于求解有N个变量si的能量E最小化问题，除非N的取值非常小，否则往往无法直接求解。由于其构型数目高达N^2。假设使用贪心算法搜索最优的构型，须要先随机选取每一个节点的起始状态，然后顺序考查每一个节点从而计算与之相联系的si=1状态和si=-1状态的能量，最后选取能够减少能量的状态迁移。

这样的推断仅仅用到了直接与之相连的具有非零权值的相邻节点。

这样的贪心搜索算法成功找到最优解的可能性非常小，由于系统经常会陷入局部能量最小值，或根本就不收敛。因此须要选择其它搜索方法。

4.模拟退火

在热力学中。固体的退火过程主要由下面三部分组成：

1. 加温过程。其目的是增强粒子的热运动。使其偏离平衡位置。当温度足够高时，固体将熔解为液体。从而消除系统原先可能存在的非均匀态，使随后进行的冷却过程以某一平衡态为起点。

   熔解过程实际是系统的熵增过程，系统能量也随温度的升高而增大。

2. 等温过程。物理学的知识告诉我们，对于与周围环境交换热量而温度不变的封闭系统。系统状态的自发变化总是朝自由能减少的方向进行，当自由能达到最小时。系统达到平衡态。
3. 冷却过程。其目的是使粒子的热运动减弱并渐趋有序，系统能量逐渐下降，从而得到低能的晶体结构。

Metropolis等在1953年提出了重要性採样法，即以概率大小接受新状态。

详细而言，在温度T时, 由当前状态i产生新状态j，两者的能量分别为Ei和Ej，若Ej小于Ei。则接受新状态j为当前状态；否则。计算概率p(∆E)：

若p(∆E)大于[0,1]区间内的随机数。则仍旧接受新状态j为当前状态；若不成立则保留i为当前状态，当中k为玻尔兹曼常数。T为系统温度。上述重要性採样过程通常称为Metropolis准则：

1. 在高温下可接受与当前状态能量差较大的新状态；
2. 而在低温下基本仅仅接受与当前能量差较小的新状态。
3. 且当温度趋于零时，就不能接受比当前状态能量高的新状态。

1983年Kirkpatrick 等意识到组合优化与物理退火的类似性。并受到Metropolis 准则的启迪，提出了模拟退火算法。模拟退火算法是基于Monte Carlo 迭代求解策略的一种随机寻优算法，其出发点是基于物理退火过程与组合优化之间的类似性，模拟退火方法由某一较高初温開始。利用具有概率突跳特性的Metropolis抽样策略在解空间中进行随机搜索，伴随温度的不断下降，反复抽样过程，终于得到问题的全局最优。对照贪心算法。模拟退火算法基本的优势在于它使系统有可能从局部最小处跳出。

对于一个优化问题：

把优化问题的求解过程与统计热力学的热平衡问题进行类比，通过模拟高温物体退火过程的方法。试图找到优化问题的全局最优或近似全局最优解。

同意随着參数的调整，目标函数能够偶尔向能量添加的方向发展（相应于能量有时上升）。以利于跳出局部极小的区域，随着假想温度的下降（相应于物体的退火），系统活动性减少，终于稳定在全局最小所在的区域。

5.两种模拟退火算法

两种模拟退火算法，即随机模拟退火和和确定性模拟退火算法的实现过程例如以下所看到的：

随机模拟退火算法收敛非常慢，部分原因在于当中搜索的所有的构型空间的离散本质，即构型空间是一个N维超立方体。

每一次搜索轨迹都仅仅能沿着超立方体的一条边。状态仅仅能落在超立方体的顶点上，因此失去了完整的梯度信息。

而梯度信息是能够用超立方体内部的连续状态值提供的。一种更快的方法就是下面的确定性模拟退火算法：

二、实验结果讨论

构造一个6-单元全连接网络，能量函数使用公式：

当中网络的连接权值矩阵例如以下：

设计步骤主要包含下面几个部分：
编敲代码[E, s_out] = RandomSimulatedAnnealing(T_max, time_max, c, s, w)。实现以上算法1所述的随机模拟退火算法。这里须要设定下面參数： T_max=10。T(m+1) =c*T(m)，c=0.9，进行实验。能量随温度下降次数的变化曲线如图2所看到的（由于模拟退火算法所得到的结果有一定的随机性，因此下面步骤均执行四次算法进行观察）。四次所得到的终于构型s如图3所看到的。

改变參数：初始温度：T_max=5，T(m+1) =c*T(m)，c=0.5，进行实验，能量随温度下降次数的变化曲线如图4所看到的。四次所得到的终于构型s如图5所看到的。

编敲代码[E, s_out] = DeterministicAnnealing(T_max, time_max, c, s, w)，实现以上算法2所述的确定性模拟退火算法。

这里须要设定下面參数： T_max=10。T(m+1) =c*T(m)，c=0.9，进行实验，能量随温度下降次数的变化曲线如图6所看到的，四次所得到的终于构型s如图7所看到的。

改变參数：初始温度：T_max=5，T(m+1) =c*T(m)，c=0.5，进行实验，能量随温度下降次数的变化曲线如图8所看到的。四次所得到的终于构型s如图9所看到的。

结论：图2、3给出了多次随机模拟退火算法的执行结果，能够看到构型s不一定全然一样；能量函数E的波形在经过若干次逐渐递减的震荡后基本收敛到全局最小值-19。当改变T(1)=5,c=0.5时，从图4中可观察到能量函数E的波形极速下降，并达到较小的值，中间少了一个温度渐变和震荡调整的过程。

图6、7给出多次确定性模拟退火算法的执行结果，且每次得到的终于构型s均一致；能量函数E的波形在经过平缓递减后收敛到全局最小值-19。不出现随机模拟退火中剧烈震荡的情况。当改变T(1)=5,c=0.5时。从图8中可观察到能量函数E的波形相同呈现出极速下降的态势。并达到较小的值，中间少了一个温度渐变和调整的过程。

综合以上的实验结果。我们发现随机退火和确定性退火均能给出类似的终于解。但对于一些大规模的现实问题，随机模拟退火的执行速度非常慢，而相比之下确定性退火算法要快非常多，有时能够快2~3个数量级。

另外，初温T(1)和温度下降系数c的选择对算法性能也有非常大影响。初温的确定应折衷考虑优化质量和优化效率。经常用法包含下面几种：

• 均匀抽样一组状态，以各状态目标值的方差为初温。
• 随机产生一组状态，确定两两状态间的最大目标值差|∆max|，然后根据差值。利用一定的函数确定初温。

  譬如T(1)=-∆/ln
  p_r。当中p_r为初始接受概率。

• 利用经验公式给出。

模拟退火算法设计中包含三个重要的函数：状态产生函数、状态接受函数、温度更新函数。同一时候在程序设计时，需遵循内循环终止准则、外循环终止准则。这些环节的设计将决定模拟退火算法的优化性能。

三、实验结果

四、简单代码实现

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 随机模拟退火函数
% 输入參数：
%   T_max：初始温度
%   time_max：最大迭代次数
%   c：温度下降比率
%   s：初始构型
%   w：权值矩阵
% 输出參数：
%   E：能量变化矩阵
%   s_out：经过算法计算后的构型
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [E, s_out] = RandomSimulatedAnnealing(T_max, time_max, c, s, w)

[x, y] = size(s);

time = 1;        % 迭代次数
T(time) = T_max; % 初始温度设置

while (time < (time_max + 1)) % (T(time) > T_min) 
    for i = 1:1000
        num = ceil(rand(1) * y); % 选择产生一个1到y之间的随机数
        for j = 1:y
            e(j) = w(num, j) * s(num) * s(j);
        end
        Ea(time) = -1 / 2 * sum(e);
        Eb(time) = -Ea(time);
        if Eb(time) < Ea(time)
            s(num) = -s(num);
        elseif (exp(-(Eb(time) - Ea(time)) / T(time)) > rand())
            s(num) = -s(num);
        else
            s(num) = s(num);
        end
    end
    % 计算能量E
    E(time) = 0;
    for it = 1:6
        for jt = 1:6
            E(time) = E(time) + w(it, jt) * s(it) * s(jt);
        end
    end
    E(time) = E(time) * (-0.5);
    s_out(time,:) = s;  % 每次形成的构型
    time = time + 1;
    T(time) = T(time - 1) * c;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 确定性模拟退火函数
% 输入參数：
%   T_max：初始温度
%   time_max：最大迭代次数
%   c：温度下降比率
%   s：初始构型
%   w：权值矩阵
% 中间函数：
%   tanh(l / T)：响应函数，该函数有一个隐含的又一次规格化的作用
% 输出參数：
%   E：能量变化矩阵
%   s_out：经过算法计算后的构型
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [E, s_out] = DeterministicAnnealing(T_max, time_max, c, s, w)

[x, y] = size(s);

time = 1;        % 迭代次数
T(time) = T_max; % 初始温度设置

while (time < (time_max + 1))
    num = ceil(rand(1) * y); % 选择产生一个1到y之间的随机数
        for j = 1:y
            e(j) = w(num, j) * s(j);
        end
    l(time) = sum(e);
    s(num) = tanh(l(time) / T(time));
    % 计算能量E
    E(time) = 0;
    for it = 1:6
        for jt = 1:6
            E(time) = E(time) + w(it, jt) * s(it) * s(jt);
        end
    end
    E(time) = E(time) * (-0.5);
    s_out(time,:) = s; % 每次形成的构型
    time = time + 1;
    T(time) = T(time - 1) * c;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%模拟退火算法实验
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
close all;

% 网络的连接权值矩阵
w = [ 0 5 -3 4 4 1;...
      5 0 -1 2 -3 1;...
     -3 -1 0 2 2 0;...
      4 2 2 0 3 -3;...
      4 -3 2 3 0 5;...
      1 1 0 -3 5 0];

num = 6; % 总共产生6个数 
s_in = rand(1,num); % 生成1和-1的随机序列 
s_in(s_in > 0.5) = 1; 
s_in(s_in < 0.5) = -1;
disp(['初始构型S为：',num2str(s_in)]);

% 下面是随机模拟退火算法
T_max = 10;        % 初始温度设置
time_max = 100;    % 最大迭代次数
c = 0.9;           % 温度变化比率
[E1, s_out1] = RandomSimulatedAnnealing(T_max, time_max, c, s_in, w);
subplot(221),plot(E1);grid on;
title(['T(1) = ',num2str(T_max),',c = ',num2str(c),',随机模拟退火算法能量变化曲线']);
disp(['T(1) = 10,c = 0.9,随机模拟退火算法终于构型S为：',num2str(s_out1(time_max,:))]);

T_max = 5;         % 初始温度设置
time_max = 100;    % 最大迭代次数
c = 0.5;           % 温度变化比率
[E2, s_out2] = RandomSimulatedAnnealing(T_max, time_max, c, s_in, w);
subplot(222),plot(E2);grid on;
title(['T(1) = ',num2str(T_max),',c = ',num2str(c),',随机模拟退火算法能量变化曲线']);
disp(['T(1) = 5,c = 0.5,随机模拟退火算法终于构型S为：',num2str(s_out2(time_max,:))]);

% 下面是确定性模拟退火算法
T_max = 10;        % 初始温度设置
time_max = 100;    % 最大迭代次数
c = 0.9;           % 温度变化比率
[E3, s_out3] = DeterministicAnnealing(T_max, time_max, c, s_in, w);
subplot(223),plot(E3);grid on;
title(['T(1) = ',num2str(T_max),',c = ',num2str(c),',确定性模拟退火算法能量变化曲线']);
disp(['T(1) = 10,c = 0.9,确定性模拟退火算法终于构型S为：',num2str(s_out3(time_max,:))]);

T_max = 5;         % 初始温度设置
time_max = 100;    % 最大迭代次数
c = 0.5;           % 温度变化比率
[E4, s_out4] = DeterministicAnnealing(T_max, time_max, c, s_in, w);
subplot(224),plot(E4);grid on;
title(['T(1) = ',num2str(T_max),',c = ',num2str(c),',确定性模拟退火算法能量变化曲线']);
disp(['T(1) = 5,c = 0.5,确定性模拟退火算法终于构型S为：',num2str(s_out4(time_max,:))]);

參考链接：http://www.cnblogs.com/growing/archive/2010/12/16/1908255.html