Simulated Annealing，中文译作“模拟退火”，是一种很经典的随机搜索优化方法。1983年，影响深远的论文“Optimization by Simulated Annealing”由S. Kirkpatrick, C. D. Gelatt, Jr. 以及 M. P. Vecchi 联合发表在Science上，到现在已经有超过4万引用量了。Simulated Annealing是一种0阶优化算法，对于一个目标函数(f(x))，Simulated Annealing只需要能够计算(f(x))在任何一点(x)的值即可；相对应的，基于梯度的优化方法为1阶优化算法，因为我们除了(f(x))，还需要知道一阶导数（梯度）的信息( abla f(x))；类推下去还有2阶优化方法，对应需要2阶导数，Hessian Matrix ( abla^2 f(x))。神经网络模型的成功很大程度上可以归功于其模型end-to-end differentiable，因此可以很方便的使用1阶梯度优化方法。0阶优化方法需要的信息更少，因而会比梯度优化稍低效一些，但另一方面适用范围更广，例如可以用来优化不可导函数。
Simulated Annealing
假定我们需要优化的目标函数为(f(x))，我们的优化问题（这里为最简单的形式，忽略了可能有的关于(x)的限制条件等）为
$$min_x f(x).$$
定义一个概率分布
$$p_T(x) = frac{exp(-frac{1}{T}f(x))}{sum_x exp(-frac{1}{T}f(x))}$$
其中(Tge 0)为一个温度变量。这个分布对应物理上的Boltzmann分布（数学上又称为Gibbs分布），可以描述气体的统计性质，这里(f(x))是能量函数。易见，当(f(x))很高的时候(p_T(x))会较小，而当(f(x))较低的时候，(p_T(x))会相应较大。这个分布有几个特殊的极端情况。第一个极端情况是当(T ightarrow infty)，也就是温度特别高的时候，对任意的(x) ，(frac{1}{T}f(x) ightarrow 0)，因而(p_{T ightarrow infty}(x))趋于一个(x)的均匀分布。这对应于物理直观上温度高的时候分子杂乱无序的运动。另一个极端情况对应于(T ightarrow 0)，也是更有意思的一个情况。为了描述这个情况，我们定义优化问题的最优值(f^*=min_x f(x))，以及最优解的集合(X^*={x^*| f(x^*) = min_x f(x)})。当(|X^*|=1)时，这个优化问题有唯一最优解，而其他时候则有多个同样好的最优解。考虑如下的变换：
$$p_T(x)=frac{exp(-frac{1}{T}(f(x)-f^*))}{sum_x exp(-frac{1}{T}(f(x)-f^*))}$$
当(T ightarrow 0)，对于非最优的(x otin X^*, exp(-frac{1}{T}(f(x)-f^*)) ightarrow 0)，而若(xin X^*, exp(-frac{1}{T}(f(x)-f^*))=1)。因此(p_{T ightarrow 0}(x))变为
$$p_{T ightarrow 0}(x) = left{egin{array}{cc} 0, & x otin X^* \ frac{1}{|X^*|}, & xin X^*end{array} ight.$$
即为(X^*)上的均匀分布。物理直观上讲，当温度特别低的时候，低能量的分子会占绝大多数，甚至由气态变为液态或者固态，分子的运动的速度也会减弱。
基于这些性质，Simulated Annealing的基本想法就是将优化问题转化为一个采样（sampling）问题。在极限情况下，如果我们能成功地从(p_{T ightarrow 0}(x))中采到样本，那么根据上面的性质，所有采到的样本都是原问题的最优解，因为此时的(p_{T ightarrow 0}(x))是(X^*)上的均匀分布。另一方面，直接从(p_{T ightarrow 0}(x))中采样往往是困难的，因此Simulated Annealing常常从一个比较高的(T)开始，在温度高的时候，采出一堆样本，然后渐渐降低(T)（也即Annealing），同时相应调整这些样本，使得最终当(T ightarrow 0)时，这些样本能够尽可能的接近(p_{T ightarrow 0}(x))的样本。
如果对于任意的(T)我们都能够从(p_T(x))中采到样本，那么这些样本上目标函数的平均值为
$$F(T) = sum_x p_T(x) f(x)$$
易见，(F(0)=min_x f(x) = f^*)。我们可以证明(F(T))是一个单调函数，其值随着(T)降低也会降低。要证明这个结论只需要说明(F(T))的导数对于任意(T)非负：
$$egin{align}frac{dF(T)}{dT} =& sum_x frac{dp_T(x)}{dT} f(x) \ =&sum_x p_T(x) f(x) frac{dlog p_T(x)}{dT}end{align}$$
这其中
$$log p_T(x) = -frac{1}{T}f(x) - log sum_x exp(-frac{1}{T}f(x))$$
因而
$$frac{dlog p_T(x)}{dT} = frac{1}{T^2}f(x) -sum_x frac{exp(-frac{1}{T}f(x))}{sum_{x'} exp(-frac{1}{T}f(x'))} frac{1}{T^2}f(x) = frac{1}{T^2}left(f(x) - sum_x p_T(x) f(x) ight)$$
再带入上面的导数，可得
$$frac{dF(T)}{dT}=sum_x p_T(x) f(x) frac{1}{T^2} left(f(x) - sum_{x'} p_T(x') f(x') ight) = frac{1}{T^2}left(sum_x p_T(x) f(x)^2 - left(sum_x p_T(x) f(x) ight)^2 ight)=frac{1}{T^2}mathrm{Var}_{p_T(x)}[f(x)]$$
这里(mathrm{Var}_{p_T(x)}[f(x)])为(f(x))这个随机变量在(xsim p_T(x))时的方差（Variance）。显然方差总是一个非负值，而(frac{1}{T^2} > 0)，所以这个导数总是非负的。更进一步，只要(f(x))的值在(p_T(x) > 0)的地方稍微有一点点变化，(f(x))的方差就为正，因此这个导数在大多数情况下都是严格大于零的。这就证明了(F(T))的单调性。这个结论告诉我们，只要能够保证我们的样本始终跟着(T)变化，这些样本的平均目标函数的值会一直单调下降，在理想情况下，最终收敛于最优解(X^*)，得到最优值(f^*)。
值得注意的是，这里的结论是收敛到整个优化问题的最优解，而非局部最优解！我们都知道梯度优化只能保证收敛到局部最优解，这样一对比，就显出了Simulated Annealing的不一样。但是，Simulated Annealing也不是总是好使的，因为这个理论保证是在我们总能从(p_T(x))中顺利的采出样本的前提下才有的，而这个前提并不是总能实现的。
Sampling
现在是时候说一说到底怎样采样了。在统计学和机器学习里面，采样作为一个方向已经被研究了很多年。最流行的从复杂分布里采样的方法是基于Markov Chain理论的。对于一个目标分布(p(x))，我们可以任意地初始化一个样本(x_0sim p_0(x))，然后不断的将(x_t)按照某一种特定的概率分布进行转化，(x_{t+1}sim q(x|x_t))，这里(q(x'|x))是从一个值(x)变为另一个值(x')的概率。当这个分布(q)（又称为transition operator因为它将一个sample (x) transition成为另一个 (x')）满足特定条件的时候，可以证明当(t ightarrow infty)时(p(x_t) ightarrow p(x))。这个transition operator首先要满足的条件是需要使目标分布(p(x))在通过(q)进行变换之后不变，也即
$$p(x') = sum_x q(x'|x) p(x), forall x'$$
直观上理解，这个等式说的是如果(x)已经遵循(p(x))分布，经过(q)的变换之后得到的(x')仍然遵循(p(x'))的分布。所以，当不断用(q)进行变换，达到目标分布(p(x))之后，继续使用(q)进行变换，样本的分布就不会变了。实际使用中，更常用的是下面这个“Detailed balance”性质，它能够保证transition operator满足上面这个不变性
$$q(x|x')p(x') = q(x'|x) p(x)$$
这个很容易证明，因为
$$p(x') = sum_x q(x|x') p(x') = sum_x q(x'|x) p(x)$$
Markov Chain的正确性最重要的条件就是transition operator满足这个性质了，其他还有一些技术型的条件，在这里略过不谈。也就是说，对于任意目标分布(p(x))，只要我们能找到满足上面这个detailed balance的一个transition operator (q(x'|x))，通过不断的对样本进行变换，我们最终能保证样本的分布趋近于(p(x))，亦即我们一定能得到正确的样本。注意这个结论并不依赖于样本的初始分布(p_0(x))，从而不管最开始样本是什么样子的，只要不断用(q)进行变换，最后的样本还是正确的。这个性质使得用Markov Chain进行采样有了很大的普适性。
然而，直接根据(p(x))的形式来找一个transition operator在很多情况下还是很困难。所幸的是，我们有很多成熟的算法能够在很多情况下构造出一个transition operator。这些算法中最著名的应该是Metropolis-Hastings算法了。这个算法提供了一个构造transition operator的框架，它的想法是这样的，对于任意的目标分布(p(x))，以及任意的transition分布(r(x'|x))（从(x)变到(x')的概率分布），我们可以构造如下的transition operator：
- 根据(r(x'|x))采样得到(x')，
- 计算接受样本的概率(A(x'|x) = minleft{1, frac{r(x|x')p(x')}{r(x'|x)p(x)} ight})，
- 以(A(x'|x))的概率接受新样本(x')，以(1-A(x'|x))的概率保留原样本(x)。
这个transition operator实际对应的transition分布为
$$q(x'|x) = left{egin{array}{ll} A(x'|x) r(x'|x), & x' e x \ sum_{x' e x} (1-A(x'|x))r(x'|x) + r(x|x), & x'=x end{array} ight.$$
易见(q(x'|x))对任意(x)而言都是一个概率分布(可以验证(sum_{x'} q(x'|x)=1))。接下来可以证明这个(q(x'|x))是满足detailed balance的，当(x=x')时，(q(x|x) p(x) = q(x|x)p(x))显然成立；当(x e x')时，
$$egin{align} q(x'|x) p(x) =& A(x'|x)r(x'|x) p(x) \=& minleft{1, frac{r(x|x')p(x')}{r(x'|x)p(x)} ight} r(x'|x)p(x) \ =& minleft{r(x'|x)p(x), r(x|x')p(x') ight} \ =& minleft{frac{r(x'|x)p(x)}{r(x|x')p(x')}, 1 ight}r(x|x')p(x') \ =& A(x|x')r(x|x')p(x') = q(x|x')p(x')end{align}$$
因此，通过Metropolis-Hastings算法构造出来的transition operator满足detailed balance，从而保证是正确的transition operator。注意这里的这个构造适用于任意的(r(x'|x))（通常称为proposal分布，(r(x'|x))甚至不要求与目标分布(p(x))有任何的关系），这就使得Metropolis Hastings算法的适用性非常的广泛了。Metropolis Hastings算法的一个特例是当(r(x'|x) = r(x|x'))的时候，接受样本的概率可以简化为(A(x'|x)=min{1, frac{p(x')}{p(x)}})。在这种情况下，只要(p(x') ge p(x))，我们总是接受新样本(x')，反之则以(p(x')/p(x))的概率接受新样本。
如果套用到我们之前的Boltzmann分布中，(p(x')/p(x) = frac{exp(-frac{1}{T}f(x'))}{exp(-frac{1}{T}(f(x)))} = expleft(-frac{1}{T}(f(x') - f(x)) ight))，在使用的(r)分布满足(r(x'|x) = r(x|x'))的情况下，只要(f(x') < f(x))，我们总是接受新样本(x')，反之则以(exp(-frac{1}{T}(f(x') - f(x))))的概率接受新样本。满足这个条件的(r)的一个例子是，当(x)是一个binary vector的时候，(r)随机的把每一个bit按照固定的概率从0翻到1，或者从1翻到0，很容以看出对于任意的(x)和(x')，(r(x'|x) = r(x|x'))。另一个例子是当(r(x'|x))是以(x)为中心的高斯分布时，同样满足这个性质。
回到Simulated Annealing
现在让我们回到Simulated Annealing，实际应用中，我们常常从一个高温度的(p_T(x))开始，生成一堆样本（最初时都不需要严格从(p_T(x))中生成样本，按照Markov Chain的性质，只要应用transition operator足够多次就能保证样本的正确性）。当(T)很大的时候，(p_T(x))比较平滑，从(p_T(x))中采样会相对容易。这一点可以从上面的Metropolis-Hastings算法的样本接受概率看出来，当(p_T(x))比较平滑时，随机的变动(x)为(x')之后，其概率变化不会很大，因而接受的概率(min{1, p(x')/p(x)})不会特别小。而另一方面当(T ightarrow 0)的时候，(p(x))变得非常不平滑，在极端情况下，只有当(xin X^*)时(p(x)=frac{1}{|X^*|})，(x otin X^*)时(p(x)=0)，也就是说，只有当随机变化的(x')恰巧是一个最优解的时候，我们才有可能接受这个新样本，要碰巧找到这样一个(x')无疑是非常困难的。
Simulated Annealing的想法是，从一个高温度的(p_T(x))开始，不断的应用transition operator，并逐渐降低温度(T)，且保证每一次都使用对应温度下的transition operator。这样渐进地推进样本群体，尽量保证每一步的样本接受概率都不会太低。
本文的结尾，我们讲最后一个结论。假定一些样本(x)原本服从的分布为(p_0(x))，现按照温度为(T)的(p_T(x))对应的transition operator (q_T(x'|x))去变换样本，则新样本的分布为(p(x') = sum_x q_T(x'|x) p_0(x))，令
$$G(T) = sum_x p(x) f(x) = sum_x sum_{x'} q_T(x|x') p_0(x') f(x)$$
我们可以证明(G(T))为(T)的单调函数，(T)越低，(G(T))也越低。这个单调性结论也可以类似的通过导数非负来证明，过程在这里略去，因为涉及到Metropolis-Hastings算法构造的transition operator，所以会更复杂一些，但原理都差不多，最后也会得到一个方差。这个结论的一个引申情况是，当(p_0(x) = p_{T'}(x))的时候，(G(T') = sum_x sum_{x'} q_{T'}(x|x') p_{T'}(x') f(x) = sum_x p_{T'}(x) f(x) = F(T'))。因为(G(T))的单调性，若(T < T')，可得(G(T) < G(T')=F(T'))。因此只要使用低温度的transition operator，对应的样本平均目标函数的值总是会变的更优。
注意，这个结论不要求样本是(p_T(x))里的样本，相反，只要是使用了温度更低的transition operator，哪怕只使用一步，我们的样本平均目标函数值都会变得更优。这个结论从另一方面保证了simulated annealing的良好收敛性。
实际情况中，怎样降低温度，何时降低温度，在每一个温度上使用多少次transition operator都是会很大程度上影响simulated annealing效果的。虽然理论上simulated annealing保证收敛，但理论保证只在无限次应用transition operator的条件下成立，如果这些参数没有设置好，往往会在没有找到最优解之前温度就已经降得过低，样本无法再变化了。所以怎么样使得simulated annealing能够发挥出效果，有时候还是一个凭经验判断的事情。