再谈模拟退火

源于现实的启发性算法：模拟退火与混合策略

前言

模拟退火（Simulated Annealing, SA）在算法竞赛圈素来以“玄学”著称，广泛地被用于骗分。这类方法看似不需要过多思考，参数一设，成败全看天命（和脸黑不黑）。

但在我上大学接触机器学习后，发现这个被戏称为“骗分大法”的算法，其实有着严谨的理论。更有意思的是，如果将它与梯度下降结合，就能搞出一种强力混合算法。

一、模拟退火算法：源于热力学的全局优化方法

1.1 物理溯源：从固体退火到算法抽象

模拟退火的灵感来自于固体退火这一工艺流程，主要有两个阶段：

1. 升温阶段：将固体加热至高温，粒子获得足够能量，处于无序的高能量状态（对应算法中的随机乱跑）；
2. 降温阶段：缓慢且稳定降低温度，粒子热运动减弱，最终形成规则晶体（对应算法收敛到最优解）。

算法层面将这一过程抽象为：

• 能量：目标函数值（Function Loss）。我们的目标是让能量越低越好（最小化问题）；
• 温度：控制探索（Exploration）与利用（Exploitation）权衡的核心参数。高温时“瞎逛”寻找新大陆，低温时“内卷”精细打磨；
• 状态转移：从当前解生成邻域新解，并决定是否接受它。

值得注意的是，粒子热运动本质上都是通过统计学来预测的，所以存在一种可能即粒子能量可能反常升高，而这是模拟退火能够不同于其他算法的核心所在，即在温度下降过程中有一定概率接受劣解（比当前最优解差的解）。

1.2 算法核心流程拆解

流程严格对应物理退火，分为五步：

1. 初始化：设定初始解、初始温度 \(T_0\)、衰减系数 \(\alpha\)、终止温度 \(T_{end}\)；
2. 迭代降温：在当前温度下，多次生成邻域新解；
3. Metropolis判定：关键一步！决定是否接受新解；
4. 温度衰减：\(T = T \times \alpha\)，逐步收缩探索范围；
5. 终止输出：温度降至冰点，输出历史最优解。

1.3 核心公式：Metropolis接受准则

这是模拟退火能够跳出局部最优解的关键。设当前能量 \(E_{now}\)，新解能量 \(E_{new}\)，能量差 \(\Delta E = E_{new}-E_{now}\)：

1. 若 \(\Delta E < 0\)：新解更优（能量降低），无条件接受；
2. 若 \(\Delta E \geq 0\)：新解更劣（能量升高），以概率 \(P\) 概率性接受：
   \[P = \exp\left(-\frac{\Delta E}{T}\right) \]

直观理解：

• 温度 \(T\) 很高时，\(\exp\) 的结果接近 1，算法接受大部分劣解
• 温度 \(T\) 接近 0 时，\(\exp\) 的结果接近 0，算法基本只接受更好的解（类似于普通贪心）

1.4 实操案例：模拟退火求解物理平衡点

以洛谷 P1337《平衡点/吊打XXX》为例。
注：虽然这个特定的物理问题本质上是一个凸优化问题（能量函数是一个单峰的函数，只有一个全局最优解，不存在局部陷阱），用梯度下降甚至三分法就能解决，但它非常适合用来演示模拟退火的标准流程。

题目大意：\(n\) 个重物系在绳子上，求绳结最终静止的 \((x, y)\) 坐标。本质是求系统总重力势能最小的点。

目标函数：

\[E(x,y) = \sum_{i=1}^n w_i \cdot \sqrt{(x-x_i)^2+(y-y_i)^2} \]

代码实现：

#include <bits/stdc++.h>
using namespace std;
double deltaT = 0.996;
double initT = 3000;
double goalT = 1e-15;
double a[10001],b[10001],c[10010];
int n;
inline double cal(double x,double y){
	double sum = 0;
	for (int i=1;i<=n;i++){
		sum+=sqrt(pow(x-a[i],2)+pow(y-b[i],2))*c[i];
	}
	return sum;
}
inline pair<double,double> SA(double x,double y){
	double nowT = initT;
	double nowx = x,nowy = y;
	double nowans = cal(nowx,nowy);
	while (nowT>=goalT){
		double newx = nowx+(rand()*2-RAND_MAX)*nowT;
		double newy = nowy+(rand()*2-RAND_MAX)*nowT;
		double newans = cal(newx,newy);
		if (newans<nowans){
			nowx = newx;
			nowy = newy;
			nowans = cal(nowx,nowy);
		}
		else {
			double tmp = exp((nowans-newans)/nowT);
			double l = (double)rand()/RAND_MAX*1.0;
			if (tmp >= l){
				nowx = newx;
				nowy = newy;
				nowans = cal(nowx,nowy);
			}
		}
		nowT*=deltaT;
	}
	return {nowx,nowy};
}
signed main(){
	srand(time(0));
	cin>>n;
	double x1 = 0,y1 = 0;
	for (int i=1;i<=n;i++){
		cin>>a[i]>>b[i]>>c[i];
		x1+=a[i];
		y1+=b[i];
	}
	pair<double,double> ans;
	double minn = 10000000000.0;
	for (int i=1;i<=10;i++){
		pair<double,double> tmp = SA(x1/n,y1/n);
		if (minn>cal(tmp.first,tmp.second)){
			ans = tmp;
			minn = cal(tmp.first,tmp.second);
		}
	}
	
	cout<<fixed<<setprecision(3)<<ans.first<<" "<<ans.second<<endl;
	return 0;
}

---

二、梯度下降

2.1 核心原理

如果说模拟退火是“到处乱跑碰运气”，那梯度下降（Gradient Descent, GD）就是“看着地形找路”。
其数学含义是：沿梯度的反方向更新参数，能最快逼近极小值。

这句话的意思是根据函数的导数的陡峭程度，判断下一步该往哪里走，走多大步。这种能较快收敛得到最佳结果的方法常被用在机器学习中求解最优参数。

参数更新公式：

\[\boldsymbol{\theta}_{new} = \boldsymbol{\theta}_{old} - \eta \cdot \nabla E(\boldsymbol{\theta}_{old}) \]

其中 \(\eta\) 是学习率（Learning Rate），决定了步子迈多大；\(\nabla E\) 是梯度向量，决定了往哪个方向迈。

2.2 原生梯度下降的“硬伤”

梯度下降是贪心算法，这导致它有几个致命弱点：

1. 易陷局部最优：如果是凹凸不平的复杂函数，它掉进第一个小坑（局部最优）就出不来了，因为它拒绝往高处走，这也是贪心的通病。
2. 对初始值极度敏感：寻优结果高度依赖初始解的选择，如果初始解在最优解的谷里则容易得到最优解，否则不容易得到。
3. 鞍点停滞：在梯度为0的平坦区域，它会直接“躺平”不再更新。

---

三、SA-GD：模拟退火与梯度下降的结合

3.1 融合动机

算法	特点	缺陷
纯模拟退火	乱跑，不怕坑	盲目，收敛慢，效率低
纯梯度下降	定向，跑得快	贪心，掉坑里出不来

我们希望结合两者：在大体方向上顺着梯度快跑（GD），但在每一步行走时加入随机的“抖动”（SA），防止掉入小坑出不来。

3.2 注意点

许多人可能会恍然大悟：原来如此，那么我先算梯度下降得到新解，再用SA判断接不接受不就万事大吉了。
这是错的！
因为梯度下降的方向必然是能量减少的方向（\(E_{new} < E_{now}\)），这会导致Metropolis准则中 \(\Delta E\) 永远小于0，算法永远接受新解。这样SA模块就失效了，算法退化回了纯梯度下降。

正确的结合方式：在梯度更新中引入噪声。

\[\boldsymbol{\theta}_{new} = \boldsymbol{\theta}_{old} - \eta \nabla E + \mathbf{Noise}(T) \]

这个 \(\mathbf{Noise}(T)\) 是一个随温度衰减的随机扰动项。

• 高温时：噪声很大，梯度指东，你可能往西，主打一个“由于能量过高而胡乱冲撞”；
• 低温时：噪声较小，不再有较大的跳跃，算法乖乖顺着梯度滑向坑底。

3.3 SA-GD 完整代码实现

注：为了代码演示的清晰性，这里我们使用线性回归作为示例。虽然线性回归本质上也是凸优化问题（单纯 GD 即可解决），但这套‘梯度+噪声’的代码框架完全适用于神经网络等复杂的非凸优化场景。

#include <bits/stdc++.h>
using namespace std;

// === 全局配置与随机数生成 ===
mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());
uniform_real_distribution<double> dist(0.0, 1.0);

inline double read(){ 
    double x = 0.0;
    int f = 1;
    char c = getchar();
    while(!isdigit(c) && c != '.'){
        if (c == '-') f = -1;
        c = getchar();
    }
    while(isdigit(c)) {
        x = x * 10 + (c - '0');
        c = getchar();
    }
    if(c == '.'){
        c = getchar();
        double decimal = 0.1;
        while(isdigit(c)){
            x += (c-'0')*decimal;
            decimal *= 0.1;
            c = getchar();
        }
    }
    return x*f;
}

int n;
vector<double> true_w;      // 实际权重 
vector<vector<double>> simple_x; // 示例x 
vector<double> simple_y;    // 示例y 

// === 核心计算函数 ===

// 计算当前的预测值 (y_pred)
// 参数: w-权重向量, b-偏差
inline vector<double> predict(const vector<double>& w, double b) {
    vector<double> preds(100);
    for (int i = 0; i < 100; i++) {
        preds[i] = b;
        for (int j = 0; j < n; j++) {
            preds[i] += w[j] * simple_x[j][i];
        }
    }
    return preds;
}

// 计算能量（即损失函数 MSE）
// 这是 SA 算法判断解优劣的标准
inline double calc_energy(const vector<double>& preds) {
    double loss = 0;
    for (int i = 0; i < 100; i++) {
        double diff = preds[i] - simple_y[i];
        loss += diff * diff;
    }
    return loss / 100.0; // Mean Squared Error
}

// 计算梯度 (Gradient)
// 返回: pair<w的梯度向量, b的梯度>
pair<vector<double>, double> calc_gradient(const vector<double>& preds) {
    vector<double> w_grad(n, 0.0);
    double b_grad = 0.0;

    for (int i = 0; i < 100; i++) {
        double diff = preds[i] - simple_y[i]; // (y_pred - y_true)
        b_grad += diff; 
        for (int j = 0; j < n; j++) {
            w_grad[j] += diff * simple_x[j][i];
        }
    }
    
    // 平均化梯度
    b_grad /= 100.0;
    for(int j=0; j<n; j++) w_grad[j] /= 100.0;
    
    return {w_grad, b_grad};
}

signed main(){
    srand(114514); // 保持数据生成的随机种子一致，方便对比
    n = read(); 
    
    // 读取真实权重
    for (int i=0; i<n; i++) true_w.push_back(read());
    double true_b = read(); 

    // === 生成模拟数据===
    for (int i=0; i<n; i++){ 
        simple_x.push_back(vector<double>());
        double maxn = -200.0, minn = 200.0;
        for (int j=0; j<100; j++){
            double tmp = (1 + rand()%100) * 0.1;
            maxn = max(maxn, tmp);
            minn = min(minn, tmp);
            simple_x[i].push_back(tmp);
        }
        // 归一化
        for (int j=0; j<100; j++){
            simple_x[i][j] = (simple_x[i][j]-minn)/(maxn-minn);
        }
    }
    
    // 生成带噪声的 y
    vector<double> tmp_y(100);
    for (int i=0; i<100; i++){
        tmp_y[i] = true_b + (1 + rand()%100)*0.01;
    }
    for (int i=0; i<n; i++){ 
        for (int j=0; j<100; j++){
            tmp_y[j] += true_w[i] * simple_x[i][j];
        }
    }
    simple_y = tmp_y;

    // === SA-GD 参数设置 ===
    int train_times = read(); 
    double init_lr = read(); 
    
    // 模拟退火参数
    double T = 100.0;       // 初始温度
    double endT = 1e-6;     // 终止温度
    double decayT = 0.99;   // 温度衰减系数
    // 根据训练次数动态调整衰减速度，确保能跑完 train_times
    if (train_times > 0) decayT = pow(endT/T, 1.0/train_times);

    // 初始化当前解
    vector<double> now_w(n, 0);
    double now_b = 0;
    
    // 记录全局最优解（防止最后停在劣解上）
    vector<double> best_w = now_w;
    double best_b = now_b;
    double min_loss = 1e18;

    // 当前的学习率
    double lr = init_lr;

    // === SA-GD 主循环 ===
    for (int iter = 1; iter <= train_times; iter++){
        // 1. 计算当前状态的预测和能量
        vector<double> current_preds = predict(now_w, now_b);
        double current_loss = calc_energy(current_preds);

        // 更新全局最优
        if (current_loss < min_loss) {
            min_loss = current_loss;
            best_w = now_w;
            best_b = now_b;
        }

        // 2. 计算梯度 (理性的方向)
        auto grads = calc_gradient(current_preds);
        
        // 3. 生成新解 (梯度下降 + 随机噪声)
        vector<double> next_w = now_w;
        double next_b = now_b;
        
        // 噪声系数随温度降低而减小
        double noise_scale = T * lr * 5.0; 

        for(int j=0; j<n; j++) {
            // 梯度项 + 噪声项
            double noise = (dist(rng) * 2 - 1) * noise_scale; 
            next_w[j] = now_w[j] - lr * grads.first[j] + noise;
        }
        double noise_b = (dist(rng) * 2 - 1) * noise_scale;
        next_b = now_b - lr * grads.second + noise_b;

        // 4. 计算新解的能量
        vector<double> next_preds = predict(next_w, next_b);
        double next_loss = calc_energy(next_preds);
        double delta_E = next_loss - current_loss;

        // 5. Metropolis 准则 (SA的核心)
        if (delta_E < 0) {
            // 新解更好，接受
            now_w = next_w;
            now_b = next_b;
        } else {
            // 新解更差，概率接受 (跳出局部最优)
            if (exp(-delta_E / T) > dist(rng)) {
                now_w = next_w;
                now_b = next_b;
            }
        }

        // 6. 协同衰减
        T *= decayT;
        // 学习率也可以随温度衰减，以便后期精细收敛
        // lr *= 0.999; 
    } 

    // 输出全局最优解
    for (int i=0; i<n; i++){
        printf("%.2lf ", best_w[i]); 
    }
    printf("%.2lf\n", best_b); 

    return 0;
}

---

四、实用调参指南

SA-GD 虽然强大，但参数也更多了。以下是“炼丹”心得：

1. 关于噪声系数：
   代码中的 noise_scale 极其重要。

   • 如果太小，算法退化成纯梯度下降，跳不出坑。
   • 如果太大，噪声淹没了梯度信息，算法退化成纯随机乱跑。
   • 经验法则：噪声的量级在初期应该能让解产生明显的坐标偏移，但不至于飞出地图边界。

2. 关于学习率 \(\eta\)：

   • 这里的学习率不仅控制步长，还配合梯度的数值大小。如果梯度值很大（如物理引力问题在近距离时），学习率要设极小（如 \(10^{-3}\)），否则就算找对地方也有可能因为步长太大跳过去。

3. 多点启动（Restart）：
   对于非凸的极难问题，不要指望一次SA-GD就能成。随机生成10个不同的初始坐标，跑10次算法，最后取最好的那个。这是工程上最稳的策略。

五、总结

1. 模拟退火（SA）：靠随机游走，概率性接受劣解，拥有全局视野，但收敛慢。
2. 梯度下降（GD）：靠学习率和导数值，贪心逼近极值，速度快但容易掉坑。
3. SA-GD 混合：在梯度下降的理性步伐中，混入随温度衰减的随机噪声。用梯度保证效率，用噪声保证可逆性。

下次再遇到让人头秃的黑箱优化问题，不妨试试给你的梯度下降加个模拟退火，说不定有奇效！