Simulated Annealing （模拟退火）
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)\)，还需要知道一阶导数（梯度）的信息\(\nabla f(x)\)；类推下去还有2阶优化方法，对应需要2阶导数，Hessian Matrix \(\nabla^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))}$$
其中\(T\ge 0\)为一个温度变量。这个分布对应物理上的Boltzmann分布（数学上又称为Gibbs分布），可以描述气体的统计性质，这里\(f(x)\)是能量函数。易见，当\(f(x)\)很高的时候\(p_T(x)\)会较小，而当\(f(x)\)较低的时候，\(p_T(x)\)会相应较大。这个分布有几个特殊的极端情况。第一个极端情况是当\(T\rightarrow \infty\)，也就是温度特别高的时候，对任意的\(x\) ，\(\frac{1}{T}f(x)\rightarrow 0\)，因而\(p_{T\rightarrow \infty}(x)\)趋于一个\(x\)的均匀分布。这对应于物理直观上温度高的时候分子杂乱无序的运动。另一个极端情况对应于\(T\rightarrow 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\rightarrow 0\)，对于非最优的\(x\notin X^*, \exp(-\frac{1}{T}(f(x)-f^*))\rightarrow 0\)，而若\(x\in X^*, \exp(-\frac{1}{T}(f(x)-f^*))=1\)。因此\(p_{T\rightarrow 0}(x)\)变为
$$p_{T\rightarrow 0}(x) = \left\{\begin{array}{cc} 0, & x\notin X^* \\ \frac{1}{|X^*|}, & x\in X^*\end{array}\right.$$
即为\(X^*\)上的均匀分布。物理直观上讲，当温度特别低的时候，低能量的分子会占绝大多数，甚至由气态变为液态或者固态，分子的运动的速度也会减弱。
基于这些性质，Simulated Annealing的基本想法就是将优化问题转化为一个采样（sampling）问题。在极限情况下，如果我们能成功地从\(p_{T\rightarrow 0}(x)\)中采到样本，那么根据上面的性质，所有采到的样本都是原问题的最优解，因为此时的\(p_{T\rightarrow 0}(x)\)是\(X^*\)上的均匀分布。另一方面，直接从\(p_{T\rightarrow 0}(x)\)中采样往往是困难的，因此Simulated Annealing常常从一个比较高的\(T\)开始，在温度高的时候，采出一堆样本，然后渐渐降低\(T\)（也即Annealing），同时相应调整这些样本，使得最终当\(T\rightarrow 0\)时，这些样本能够尽可能的接近\(p_{T\rightarrow 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\)非负：
$$\begin{align}\frac{dF(T)}{dT} =& \sum_x \frac{dp_T(x)}{dT} f(x) \\ =&\sum_x p_T(x) f(x) \frac{d\log 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{d\log 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)\right)$$
再带入上面的导数，可得
$$\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')\right) = \frac{1}{T^2}\left(\sum_x p_T(x) f(x)^2 - \left(\sum_x p_T(x) f(x)\right)^2\right)=\frac{1}{T^2}\mathrm{Var}_{p_T(x)}[f(x)]$$
这里\(\mathrm{Var}_{p_T(x)}[f(x)]\)为\(f(x)\)这个随机变量在\(x\sim 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_0\sim 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\rightarrow \infty\)时\(p(x_t)\rightarrow 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) = \min\left\{1, \frac{r(x|x')p(x')}{r(x'|x)p(x)}\right\}\)，
- 以\(A(x'|x)\)的概率接受新样本\(x'\)，以\(1-A(x'|x)\)的概率保留原样本\(x\)。
这个transition operator实际对应的transition分布为
$$q(x'|x) = \left\{\begin{array}{ll} A(x'|x) r(x'|x), & x'\ne x \\ \sum_{x'\ne x} (1-A(x'|x))r(x'|x) + r(x|x), & x'=x \end{array}\right.$$
易见\(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\ne x'\)时，
$$\begin{align} q(x'|x) p(x) =& A(x'|x)r(x'|x) p(x) \\=& \min\left\{1, \frac{r(x|x')p(x')}{r(x'|x)p(x)}\right\} r(x'|x)p(x) \\ =& \min\left\{r(x'|x)p(x), r(x|x')p(x')\right\} \\ =& \min\left\{\frac{r(x'|x)p(x)}{r(x|x')p(x')}, 1\right\}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)))} = \exp\left(-\frac{1}{T}(f(x') - f(x))\right)\)，在使用的\(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\rightarrow 0\)的时候，\(p(x)\)变得非常不平滑，在极端情况下，只有当\(x\in X^*\)时\(p(x)=\frac{1}{|X^*|}\)，\(x\notin 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能够发挥出效果，有时候还是一个凭经验判断的事情。
