【Coel.学习笔记】随机化算法：模拟退火与爬山法

简介

模拟退火 （\(\text{Simulate Anneal}\)） 和爬山法是随机化算法，二者的原理都在于通过随机生成答案并检查，把答案逐步缩小在一个可行的区间，尽可能地靠近正确答案。
在考场中，如果某道题目的正解比较难想（例如性质复杂的贪心）或者对应的函数取值可能数极大且难以二分/三分求解时，模拟退火和爬山法可以一定程度上得到正确答案并获得可观的分数。
顺带一提，现在的人工智能其实也是基于随机化算法的，例如遗传算法就是一种典型的随机化算法。
高情商：利用人工智能算法获取近似解
低情商：随机骗分

原理

下面主要讲模拟退火，爬山法只提一点。

爬山法

爬山法通常用来求解一个函数的最值，当然得是一个凸函数。

可能有人要问了，既然是凸函数，为什么不用三分？

一方面，某些时候看不出是凸函数，爬山法硬上有时能过；另一方面，如果函数有多维，就要写嵌套三分，三维套两个，四维套三个……如果有 \(k\) 维（比如下面这一题），即使耐得住性子写完，时间复杂度也是 \(O(n\log ^k n)\)，效率很低。

回到爬山法。爬山法的原理可以总结成“从一个胜利走向另一个胜利”（经典英语造句），即每次寻找一个最优解，按照当前给定的参数控制爬山距离，并根据当前局面选定爬山的方向。当然有可能跳过最优解或者陷在某一个局部最优解的位置，所以爬山法的局限性比较大。

模拟退火

顾名思义，模拟退火就是模拟物理学上的退火过程。类比分子在高温时能量更高、碰撞次数更大的原理（参见《化学反应原理》P28），我们定义几个量：

• 温度 \(T\)：相当于答案区间，控制每次随机区域。一开始的 \(T\) 为较大值，接下来逐步减小。
  \(T\) 有初始态 \(T_0\)、终止态 \(T_k\) 和降温系数 \(d\) 三个关联量。一般来说 \(T\) 取值较大，\(T_k\) 接近 \(0\)，\(d\) 略小于 \(1\)。这样每次让 \(T\) 乘以 \(d\) 就可以实现“降温”，让答案稳定下来。
• 能量差 \(\Delta E\)，表示新得到答案和原本答案的差值。如果 \(\Delta E<0\)，说明答案更优，应取该答案；反之，以一定概率取这个答案，从而减小落入局部最优解的可能性。
  通常取概率值为 \(e^{\frac{-\Delta E}{T}}\)，这样可以保证能量差越高取答案的可能越大。

下面这个图反映了落入局部最优解的坏处（来自 OI-Wiki）。

例题讲解

[JSOI2008]球形空间产生器 - 爬山法浅谈

洛谷传送门
给定 \(n+1\) 个 \(n\) 维空间上的点，求经过这些点的球的球心坐标。

解析：这道题正解是用高斯消元求解 \(n\) 元一次方程组（球的方程可以通过作差消掉高次项），效率确实很高，但使用爬山法求解写起来更直观。

具体操作为：把球心初始化成这些点的重心（即坐标均值），然后每次爬山时，通过尽可能让距离平均来寻找球心，远了就拉近一点，近了就推远一点。

当然这道题也可以用模拟退火做，这里只是简单介绍一下爬山法。事实上爬山法比较鸡肋，一般不会用quq

#include <iostream>
#include <iomanip>
#include <cmath>

using namespace std;

const int maxn = 20;

int n;
double dis[maxn][maxn], cans[maxn], ans[maxn], delta[maxn]; 

void hillClimbing() {
	double sum = 0;
	for (int i = 1; i <= n + 1; i++) {
		cans[i] = delta[i] = 0;
		for (int j = 1; j <= n; j++)
			cans[i] += (dis[i][j] - ans[j]) * (dis[i][j] - ans[j]);
		cans[i] = sqrt(cans[i]);
		sum += cans[i] / (n + 1); //求出距离均值
	}
	for (int i = 1; i <= n + 1; i++)
		for (int j = 1; j <= n; j++)
			delta[j] += (cans[i] - sum) * (dis[i][j] - ans[j]) / sum;
}

int main(void) {
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	cin >> n;
	for (int i = 1; i <= n + 1; i++)
		for (int j = 1; j <= n; j++) {
			cin >> dis[i][j];
			ans[j] += dis[i][j] / (n - 1); //预处理重心
		}
	for (double step = 1e4; step > 1e-6; step *= 0.99995) {
		hillClimbing();
		for (int i = 1; i <= n; i++)
			ans[i] += delta[i] * step;
	}
	for (int i = 1; i <= n; i++)
		cout << fixed << setprecision(3) << ans[i] << ' ';
	return 0;
}

[POJ2420]A Star not a Tree? - 模拟退火入门

洛谷传送门
给定平面直角坐标系上的 \(n\) 个点，找到一个点使得到这些点的距离最小（即费马点），输出距离之和（保留整数）。

解析：这道题的答案是一个凸函数，可以用三分法求解，这里演示一下模拟退火。
模拟退火写起来其实非常无脑xwx

#include <iostream>
#include <algorithm>
#include <cmath>
#include <random>
#include <iomanip>

using namespace std;

const int maxn = 150;
const double maxTime = 0.95;

random_device rd;
mt19937 Rand(rd());
uniform_real_distribution<double> Dis(0, 1); //随机三件套

struct Point {
    double x, y;
} q[maxn];

int n;
double ans;

inline double rand(double l, double r) { //得到 [l, r) 的数字
    return (double) Dis(Rand) * (r - l) + l;
}

inline double dis(Point a, Point b) {
    double dx = a.x - b.x, dy = a.y - b.y;
    return sqrt(dx * dx + dy * dy);
}

double calcDistance(Point p) { //计算距离和，顺便更新答案
    double res = 0;
    for (int i = 1; i <= n; i++) res += dis(p, q[i]);
    ans = min(ans, res);
    return res;
}

void simulateAnneal() {
    Point cur = {rand(0, 1e4), rand(0, 1e4)};
    for (double T = 1e4; T > 1e-4; T *= 0.99) {
        Point p = {rand(cur.x - T, cur.x + T), rand(cur.y - T, cur.y + T)}; //在范围内随机取一个新点
        double delta = calcDistance(p) - calcDistance(cur);
        if (exp(-delta / T) > rand(0, 1)) cur = p; //概率取解
    }
}

int main(void) {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    int T;
    cin >> T;
    for (int j = 1; j <= T; j++) {
        ans = 1e8;
        cin >> n;
        for (int i = 1; i <= n; i++) cin >> q[i].x >> q[i].y;
        for (int i = 1; i <= 100; i++) simulateAnneal();
        /*如不是多组数据则可以利用“卡时”技巧，下一题会详细讲解*/
        cout << fixed << setprecision(0) << ans << '\n';
        if (j != T) cout << '\n';
    }
    return 0;
}

[JSOI/AHOI2014]保龄球 - 非数学的模拟退火

洛谷传送门
进行若干轮保龄球，对于每一轮，全倒时下一轮有双倍奖励分，补中时下一轮第一次有双倍奖励分，失误无奖励分。特别地，如果第 \(n\) 轮全中则可以额外进行一个奖励轮且有双倍奖励分，奖励轮只能进行一次。

某人可以随意调整这 \(n\) 轮的得分情况，但必须保证调整前后奖励轮的情况一致。求出在该操作下可以得到的最高分数。

解析：这道题看起来不像上一题那么“数学”，不过本质上都是求最优解；如果用动态规划或者搜索，会发现状态太多不好储存，且速度很慢，考虑使用模拟退火。事实上，之前提到的很多最优解问题（动态规划、二分、网络流……）都可以用模拟退火，只不过答案不够精确罢了。

不难发现，这道题的调整为保龄球情况的排列，对于排列上的模拟退火，我们可以随机选取两个数据交换一下，把它作为新的状态。如果状态合法，就按照解的优劣，以一定概率更新答案；否则把位置交换回去，还原现场。此时温度 \(T\) 的“区间长度”意义被淡化，仅仅用来控制退火次数以及取答案的概率。

由于模拟退火次数比较难估计，这里采用卡时技巧。C++11 的chrono 头文件提供了一个叫 chrono::system_clock::now().time_since_epoch().count() 的函数（好长的名字……），可以得到以纳秒为单位的时间，除以 \(10^9\) 就是秒为单位了。

时间上限可以设置为一个略小于时间限制的值（不能等于，否则还没输出就 TLE 了），实际发现数据很水， \(0.05\) 秒也可以过。

#include <iostream>
#include <cmath>
#include <chrono>
#include <random>
#include <vector>

using namespace std;

const int maxn = 80;

int n, ans;
bool bonus;

random_device rd; 
mt19937 Rand(rd());
uniform_real_distribution<double> check(0, 1);

pair<int, int> a[maxn];

inline auto getTime() { //得到秒为单位的时间（函数名有点长……）
    return chrono::system_clock::now().time_since_epoch().count() / 1e9;
}

int newRecord() { // 计算当前得分
    int res = 0;
    for (int i = 1; i <= n + bonus; i++) {
        res += a[i].first + a[i].second;
        if (i <= n) {
            if (a[i].first == 10) res += a[i + 1].second;
            if (a[i].first + a[i].second == 10) res += a[i + 1].first;
        }
    }
    ans = max(ans, res);
    return res;
}

void simulateAnneal() {
    uniform_int_distribution<int> dis(1, n + bonus + 1);
   for (double T = 1e4; T > 1e-4; T *= 0.99) {
       int pos1 = dis(Rand), pos2 = dis(Rand), now = newRecord();
       swap(a[pos1], a[pos2]);
       if (!((a[n].first == 10) ^ bonus)) {
           int delta = newRecord() - now;
           if (exp(delta / T) < check(Rand)) swap(a[pos1], a[pos2]); //答案较劣，一定概率还原
		/*这里 delta 不带负号是因为答案劣才还原，而不是答案优还原*/
       } else swap(a[pos1], a[pos2]); // 答案不合法，直接还原
   }
}

int main(void) {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    auto start = getTime(); //起始时间
    cin >> n;
    for (int i = 1; i <= n; i++) cin >> a[i].first >> a[i].second;
    if (a[n].first == 10) { //判断是否有额外轮
        bonus = true;
        cin >> a[n + 1].first >> a[n + 1].second;
    }
    while (getTime() - start < 0.05) simulateAnneal();
    cout << ans;
    return 0;
}

[HAOI2006] 均分数据 - 模拟退火配合贪心

洛谷传送门
已知 \(n\) 个正整数 \(a_1,a_2 ... a_n\) 。今要将它们分成 \(m\) 组，使得各组数据的数值和最平均，即各组的均方差最小。均方差公式如下：

\[\sigma = \sqrt{\frac 1n \sum\limits_{i=1}^n(\overline x - x_i)^2},\overline x = \frac 1n \sum\limits_{i=1}^n x_i \]

其中 \(\sigma\) 为均方差，\(\bar{x}\) 为各组数据和的平均值，\(x_i\) 为第 \(i\) 组数据的数值和。

解析：暴力做法为枚举每个数据所在的位置，方案可以（？）用 \(n\) 位 \(m\) 进制储存。当然时间复杂度为 \(O(m^n)\)，跑得过就有鬼了ovo

由于更换某几个数据后答案差别不会太大，所以对应的“函数”具有一定的连续性，可以使用模拟退火。事实上，这是模拟退火可用的一个重要条件。

但是如果直接用模拟退火做的话答案区间过大，很难得到正确答案，考虑进行优化。使用贪心思想，从前往后找要选择的数据，每次把数据插入到数值和最小的数组中，可以证明这样做可以让均方差尽可能小。当然数据插入顺序要用模拟退火随机修改，从而提高贪心的正确率。

当然本题也可以不用模拟退火，而是 random_shuffle 随机生成数据顺序，不过正确率没有模拟退火高。

#include <iostream>
#include <cmath>
#include <chrono>
#include <iomanip>
#include <algorithm>
#include <random>
#include <cstring>

const int maxn = 50;

using namespace std;

int n, m;
int a[maxn], s[maxn];
double ans = 1e8;

random_device rd;
mt19937 Rand(rd());
uniform_real_distribution<double> check(0, 1);

inline auto getTime() {
	return chrono::system_clock::now().time_since_epoch().count() / 1e9;
}

double calcVariance() {
	memset(s, 0, sizeof(s));
	double sum = 0, res = 0;
	for (int i = 1; i <= n; i++) {
		int k = 1;
		for (int j = 2; j <= m; j++)
			if (s[j] < s[k]) k = j; // 贪心寻找数据和最小的数组
		s[k] += a[i];
	}
	for (int i = 1; i <= m; i++) sum += (double) s[i] / m;
	for (int i = 1; i <= m; i++)
		res += (s[i] - sum) * (s[i] - sum);
	res = sqrt(res / m);
	ans = min(ans, res);
	return res;
}

void simulateAnneal() {
	static uniform_int_distribution<int> dis(1, n);
	random_shuffle(a + 1, a + n + 1); //初始化一个排列（其实也可以不初始化）
	for (double T = 1e4; T >= 1e-4; T *= 0.99) {
		int x = dis(Rand), y = dis(Rand);
		double now = calcVariance(), delta;
		swap(a[x], a[y]);
		delta = calcVariance() - now;
		if (exp(-delta / T) < (double) check(Rand)) swap(a[x], a[y]);
	}
}

int main(void) {
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	auto start = getTime();
	cin >> n >> m;
	for (int i = 1; i <= n; i++) cin >> a[i];
	simulateAnneal(); //数据很水，所以跑一遍就够了
	cout << fixed << setprecision(2) << ans;
	return 0;
}