随机化算法
UKE_Automation · · 算法·理论
1 随机化算法简介
随机化算法，是一种十分玄学的做法。
百度百科对其的定义是：
随机化算法（randomized algorithm），是这样一种算法，在算法中使用了随机函数，且随机函数的返回值直接或者间接的影响了算法的执行流程或执行结果。就是将算法的某一步或某几步置于运气的控制之下，即该算法在运行的过程中的某一步或某几步涉及一个随机决策，或者说其中的一个决策依赖于某种随机事件。
听上去十分高级，实际上就是，程序的行为及结果会由随机数来控制的算法，叫随机化算法。
在 OI 中，我们最常见的随机化算法是输出一个随机数作为结果。这常被用于走投无路时的最后的骗分之举，但这确实也可以称之为随机化算法。
而一些高级的算法，譬如此文要提到的爬山法和模拟退火，本质上与输出随机数是没有区别的。然而，他们终究是比随机数更加准确的，对于有些题目能够拿到
本文将介绍两种基本随机化算法：爬山法、模拟退火。
2 爬山法
2.1 概念
爬山法是一种局部择优的算法，它采用启发式方法。
什么意思呢？我们举一个例子，比如我们将题目转化为一个
贪心地想，我们对于一个
持续下去，我们就能找到最优解……吗？
考虑下图：
我们在找的时候大概率会停在红色处，然而绿色才是最优解。但是，这种方法对于单峰函数一定有效。
而这种策略，我们称之为爬山法，因为他确实像一个
2.2 实现
我们考虑一个问题：什么叫轻微扰动？
我们引入一个关键：温度参数。每次按照温度参数的一定比例进行调整，然后每次结束后进行降温。
我们先看一个例子：
double t = 10000, p = 0.997;
while(t >= 1e-15) {
double nx = ansx + t * (rand() * 2 - RAND_MAX);
int nw = calc(nx);
int delta = nw - answ;
if(delta > 0) {
ansx = nx;
ansy = ny;
answ = nw;
}
t *= p;
}此处，t 为温度参数，p 为降温速率，1e-15 为结束标志。我们对当前最优解 ansx 进行扰动得到 nx，计算函数值并比较，如果更有就转移。最后一轮迭代完成，降温。
在实际中，
这是爬山法的基本代码，但是这很明显有一个问题：
它只对单峰函数保证正确率。
怎么解决呢？继续看
3 模拟退火
3.1 概念
模拟退火的名称，来源于金属退火。但是你并不需要了解金属退火是什么
我们看爬山法为什么不对，爬山法只选取了较优解，而完全抛弃掉了较次的解。其实与贪心的缺点是一样的：目光短浅。
那么，模拟退火的任务就很简单了：利用随机化，保留较次解。
3.2 过程
先概括：如果更优就接受，否则以一定概率接受。
我们继续沿用爬山法的温度系统（其实准确的说爬山法的温度系统是参考模拟退火），设当前温度为
这样我们就以一个概率去接受一个当前较次的解，以换取更优的解。
还有一个小技巧：在模拟退火后，在当前答案周边在进行几次小幅度的扰动，尝试得到更优解。
同时，我们在模拟退火过程中取答案最大值，可以进一步提高精确度。
最后放一张图，红线就是当前的
3.3 实现
3.3.1 主要函数
void SA() {
double t = 10000, p = 0.997;
while(t >= eps) {
double nx = ansx + t * (rand() * 2 - RAND_MAX);
int nw = calc(nx);
int delta = nw - answ;
if(delta > 0) {
ansx = nx;
answ = nw;
}
else if(exp(-delta / t) * RAND_MAX < rand()) { //
ansx = nx; //
} //
t *= p;
}
for(int i = 1; i <= 1000; i++) {
double nx = ansx + t * (rand() * 2 - RAND_MAX);
int nw = calc(nx);
int delta = nw - answ;
if(delta > 0) {
ansx = nx;
answ = nw;
}
}
}仔细观察后发现，只有带注释的三行与爬山法不同，这就是接受较次解的部分。
下面就是刚刚提到的再次小范围迭代，其实就是爬山法（因为此时已经几乎是一个单峰函数了）。
3.3.2 小技巧
在实际应用中，我们往往需要跑多个模拟退火以得到较为准确的方法。一种粗暴的方式是多写几遍。但是我们有文明的做法，是一个暴力算法的常用技巧：卡时。
如下：
while((double)clock()/CLOCKS_PER_SEC < 0.9)
SA();3.3.3 完整代码
在模拟退火中，影响答案的往往就是温度相关的参数以及 RP。
以 P1337 [JSOI2004] 平衡点 / 吊打XXX 为例，给出模板：
#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int Maxn = 2e5 + 5;
int n;
struct node{
double x, y, w;
}a[Maxn];
double ansx, ansy, answ;
double calc(double x, double y) { //物理知识
double res = 0;
for(int i = 1; i <= n; i++) {
double dx = x - a[i].x;
double dy = y - a[i].y;
res += sqrt(dx * dx + dy * dy) * a[i].w;
}
return res;
}
void SA() {
double t = 100000, p = 0.997;
while(t >= 1e-15) {
double nx = ansx + t * (rand() * 2 - RAND_MAX);
double ny = ansy + t * (rand() * 2 - RAND_MAX);
double nw = calc(nx, ny);
double delta = nw - answ;
if(delta < 0) { //注意这个退火求最小值
answ = nw;
ansx = nx;
ansy = ny;
}
else {
if(exp(-delta / t) * RAND_MAX > rand()) { //注意这个退火求最小值
ansx = nx;
ansy = ny;
}
}
t *= p;
}
for(int i = 1; i <= 1000; i++) {
double nx = ansx + t * (rand() * 2 - RAND_MAX);
double ny = ansy + t * (rand() * 2 - RAND_MAX);
double nw = calc(nx, ny);
double delta = nw - answ;
if(delta < 0) {
answ = nw;
ansx = nx;
ansy = ny;
}
}
}
int main() {
srand(time(NULL));
cin >> n;
for(int i = 1; i <= n; i++) {
cin >> a[i].x >> a[i].y >> a[i].w;
ansx += a[i].x, ansy += a[i].y;
}
ansx /= n, ansy /= n; //以平均值作为迭代起点
answ = calc(ansx, ansy);
while((double)clock()/CLOCKS_PER_SEC < 0.9) //卡时
SA();
printf("%.3lf %.3lf", ansx, ansy);
return 0;
}我们发现，求最大值和最小值的代码是不同的。事实上，对于求最大值：
if(delta < 0)
if(exp(-delta / t) * RAND_MAX > rand())而最小值：
if(delta > 0)
if(exp(-delta / t) * RAND_MAX < rand())只有两行不同。
3.4 习题
P5544 [JSOI2016] 炸弹攻击1
代码如下：
#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int Maxn = 2e5 + 5;
const double eps = 1e-18;
int n, m, r;
struct building {
int x, y, r;
}a[Maxn];
struct node {
int x, y;
}b[Maxn];
double ansx, ansy;
int answ;
int calc(double x, double y) {
double maxr = r;
for(int i = 1; i <= n; i++) {
double dx = x - a[i].x;
double dy = y - a[i].y;
double dr = sqrt(dx * dx + dy * dy);
maxr = fmin(maxr, dr - a[i].r);
} //最大半径
int cnt = 0;
for(int i = 1; i <= m; i++) {
double dx = x - b[i].x;
double dy = y - b[i].y;
double dr = sqrt(dx * dx + dy * dy);
if(dr <= maxr) {
cnt++;
}
} //统计杀伤数
return cnt;
}
void SA() {
double t = 10000, p = 0.997;
while(t >= eps) {
double nx = ansx + t * (rand() * 2 - RAND_MAX);
double ny = ansy + t * (rand() * 2 - RAND_MAX);
int nw = calc(nx, ny);
int delta = nw - answ;
if(delta > 0) {
ansx = nx;
ansy = ny;
answ = nw;
}
else if(exp(-delta / t) * RAND_MAX < rand()) {
ansx = nx;
ansy = ny;
}
t *= p;
}
for(int i = 1; i <= 1000; i++) {
double nx = ansx + t * (rand() * 2 - RAND_MAX);
double ny = ansy + t * (rand() * 2 - RAND_MAX);
int nw = calc(nx, ny);
int delta = nw - answ;
if(delta > 0) {
ansx = nx;
ansy = ny;
answ = nw;
}
}
}
void solve() {
while((double)clock()/CLOCKS_PER_SEC < 0.9) //卡时
SA();
}
int main() {
srand(time(NULL));
ios::sync_with_stdio(0);
cin >> n >> m >> r;
for(int i = 1; i <= n; i++) {
cin >> a[i].x >> a[i].y >> a[i].r;
}
for(int i = 1; i <= m; i++) {
cin >> b[i].x >> b[i].y;
ansx += b[i].x, ansy += b[i].y;
}
ansx /= m, ansy /= m; //以平均值作为迭代起点
answ = calc(ansx, ansy);
solve();
cout << answ;
return 0;
}（其实这两道题用爬山法也可以通过）
4 总结
这是常见的两个随机化算法。但毕竟是随机化，有时可能需要十几次提交才有一个 AC。但是，对于蒟蒻不会的题，能够拿下七八十分，很不错了。
（第二题就是江苏 2016 的省选，一道紫题，模拟退火可以玄学 AC。对于加强后的题目，是一道黑题，模拟退火仍然可以拿到不少分数。）