高级算法指北——从爬山算法到模拟退火
I 从爬山算法到模拟退火
i 定义
对于一个问题的解 \(x\)，我们常常能找到一个变量 \(y\)，使其仅与 \(x\) 相关联。我们称 \(y=f(x)\) 为解的对应函数。解取到 \(x_0\) 时对应的 \(y_0\)，称之为解 \(x_0\) 对应的函数值。一个问题的解的对应函数可能有很多，我们常常选择最显然、最容易解题的一个分析。某个问题的解可能不是一个数，也可能是一个状态或一种方案。
爬山算法/模拟退火是一种随机化算法。常用于问题的解的对应函数非单调函数时的求解或解的对应函数的最值。
“退火是一种金属热处理工艺，指的是将金属缓慢加热到一定温度，保持足够时间，然后以适宜速度冷却的过程”（节选自百度百科）。我们用程序模拟“以适宜速度冷却”这一过程，并在逐步冷却的过程中求解。
ii 单峰函数上的求解：爬山算法
对于单峰函数上的求解问题，我们可以选定一个初始答案，然后随机往左右两侧跳，比较新的答案与现在的答案对应函数值的大小关系。若新的答案更优，则跳到新的答案处。在更新的过程中，我们逐步减小随机跳的范围，当范围足够小时停止操作，则那个时候的答案就是最终的答案。而在整个过程中，随机范围的减小就对应着温度的降低。这个过程类似于往上爬山，故称之为“爬山算法”。
我们可以自己设出初始温度 \(temp\)，和每次操作降温的比例 \(d\)。一般而言，\(temp\) 是一个大整数，根据答案范围的不同而不同；而 \(d\) 是一个小于但接近 \(1\) 的数，取值范围大约在 \([0.985,0.997]\) 之间。由于题目的不同，你需要根据实际运行的结果调整各个变量的取值。每次随机下一个答案时，可以写作 ans+=(rand()*2-RAND_MAX)*temp，既保证了范围不过小，又保证了范围随温度的降低而减小。这里的 RAND_MAX 是 \(rand\) 函数出现的随机数最大值，是系统自带的一个常数。每次随机答案结束后，需要将 temp*=d，即降温。待降温到一个下限时，就停止随机。这个下限一般取 \(10^{-10}\) 到 \(10^{-15}\)，同样需要根据实际情况调整。
const double d=0.996;//降温参数
double temp=1000;//初始温度
while(temp>1e-15){
double nw=ans+(rand()+rand()-RAND_MAX)*temp,nwans=getans(nw);//getans 是求当前答案下的函数值
//这里假定是函数值越大越好，即单峰
if(nwans>maxn)maxn=nwans,ans=nw;//答案更优就更新答案
temp*=d;//降温
}
对于单峰函数上的求解，其正确性显然。因为我们每次更新到的答案只会往最终的那个峰值逼近。但对于多个峰值的函数，爬山算法很有可能求到某个值域内的局部最优解，而非整体最优解。优化方案可以是对于值域分块，每一块内爬山后，对所有爬山结果取最优值；也可以增大容错率，对于非最优的答案，以一定的概率去取，使我们能够跳出现目前的较劣的峰，而到达更优的峰。这就是模拟退火算法。
iii 模拟退火（Simulated Annealing）——对多峰下爬山的正确率优化
如同前面所说的，每随机到一个答案，我们就与现目前的最优答案进行比较。若新的答案更优，我们直接取新的答案；否则我们以一定的概率，跳到最新的答案上去，但最优答案仍然不变。随着温度的降低，我们逐步锁定了最优的峰，此时跳到较劣答案上的概率应当降低（注意：此时的答案范围不一定必要缩小）。此外，较劣答案与最优答案的函数值之差越大，跳到较劣答案上的概率也应当越低。
设跳到较劣答案上的概率函数为 \(p(x)\)，其中 \(x\) 为最新的较劣答案的函数值与目前最优答案的函数值之差。一般而言，我们选用函数 \(p(x)=e^{-\frac{x}{temp}}\)。不难发现，该函数值随 \(temp\) 的减小而减小，随 \(x\) 的减小而增大。我们可以用\(exp\)函数求得上述函数值，即 p=exp(-x/temp)。求得概率后，我们可以通过随机一个 \([0,1]\) 范围内的实数，并比较它与概率值的大小，确定是否跳到较劣答案上去。求随机实数可以写作 (double)rand()/1.0/RAND_MAX。
实现与爬山算法基本类似，只需要多加上跳到较劣答案上的判断即可。
void SA(){
double temp=1000;//温度
while(temp>1e-15){
double nw=ans+(rand()+rand()-RAND_MAX)*temp,nwans=getans(nw);//getans 是求当前答案下的函数值
//这里假定是函数值越大越好，即单峰
if(nwans>maxn)maxn=nwans,ans=nw,maxans=nw;//答案更优直接选，存下全局最优答案
else if(exp(-(maxn-nwans)/temp)>(double)rand()/1.0/RAND_MAX){//答案不优，有一定的概率跳过去.
ans=nw;//最优函数值不变，只更新下一步判断的答案
}
temp*=d;//降温
}
}
注意：我们一般把整个退火过程中的全局最优答案存下来，作为最终答案输出，而不是输出最后一次跳到的那个答案。
II 模拟退火的正确率优化
i 分块退火
当函数的峰很多时，模拟退火难以跑出最优解。可以将解的值域分块，对每一块内分别做模拟退火，最后对退火出的所有答案取最优即可。
ii 最优点偏移
我们可以在求出全局最优答案后，在全局最优答案附近执行若干次爬山的判断，试图寻找这附近有微小差异的更优解。这里的随机次数不宜过多，否则会挤占前面模拟退火算法的时间。
rep(i,1,100){//在最优解附近随机跳，使得答案更接近. 这里就不跳到稍微差一点的解上去了.不必要太多次，反而容易减少大SA的次数，导致WA.
double nw=maxans+(rand()+rand()-RAND_MAX)*temp,nwans=getans(nw);
if(nwans>maxn)maxn=nwans,ans=nw;
//不必降温了
}
iii 多次退火，循环卡时
我们可以多次进行退火算法，每一次从上一次找到的最优答案处开始退火，直到快要超时的时候再停止。
double s=clock();
while(1){
SA();
if((double)(clock()-s)/1.0/CLOCKS_PER_SEC>0.98)break;
}//clock()/CLOCKS_PER_SEC 表示程序运行的时间.
需要注意的是，有的操作在实际运行中可能导致负优化，需要结合实际情况选用。
III 经典例题
[JSOI2004] 平衡点
实际上要求的是与各个点的距离乘上权值的和最小的点，显然最终的解与各点的距离乘权值的和有函数关系，要求的是这个函数值最小的点。不难发现这个函数在 \(x\) 轴和 \(y\) 轴 方向上均具有单峰性（考虑正交分解），故爬山和模拟退火均可以使用。
事实上，只要各个参数选择合理，两种算法都可以达到极高的正确率。
- 爬山版本
#include<bits/stdc++.h>
#define rep(i,j,k) for(int i=j;i<=k;i++)
#define repp(i,j,k) for(int i=j;i>=k;i--)
using namespace std;
typedef long long ll;
const int N=1005,M=60,mo=1e9+7,inf=1e9+7;
const double d=0.996;//降温参数
void read(int &p){
int x=0,w=1;
char ch=0;
while(!isdigit(ch)){
if(ch=='-')w=-1;
ch=getchar();
}
while(isdigit(ch)){
x=(x<<1)+(x<<3)+ch-'0';
ch=getchar();
}
p=x*w;
}
//带权费马点问题，爬山也可以解决.
int n,x[N],y[N],w[N];
double ansx,ansy,minn;
double dis(double x,double y){
return sqrt(x*x+y*y);
}
double getans(double nx,double ny){//求权值和
double ret=0;
rep(i,1,n)
ret+=dis(nx-x[i],ny-y[i])*w[i];
return ret;
}
void SA(){
double temp=1000;//温度
while(temp>1e-15){
double nwx=ansx+(rand()+rand()-RAND_MAX)*temp,nwy=ansy+(rand()+rand()-RAND_MAX)*temp,nw=getans(nwx,nwy);
if(nw<minn)minn=nw,ansx=nwx,ansy=nwy;//答案更优直接选
temp*=d;//降温
}
}
int main(){
read(n);
double s=clock();
rep(i,1,n)
read(x[i]),read(y[i]),read(w[i]),ansx+=x[i],ansy+=y[i];
ansx=ansx/1.0/n,ansy=ansy/1.0/n,minn=getans(ansx,ansy);
while(((double)(clock()-s)/1.0/CLOCKS_PER_SEC)<0.97)
SA();
printf("%.3lf %.3lf",ansx,ansy);
return 0;
}
在总共的 \(3\) 次提交中，该程序的正确率达到了 \(100\%\)。
- 模拟退火版本
#include<bits/stdc++.h>
#define rep(i,j,k) for(int i=j;i<=k;i++)
#define repp(i,j,k) for(int i=j;i>=k;i--)
using namespace std;
typedef long long ll;
const int N=1005,M=60,mo=1e9+7,inf=1e9+7;
const double d=0.996;//降温参数
void read(int &p){
int x=0,w=1;
char ch=0;
while(!isdigit(ch)){
if(ch=='-')w=-1;
ch=getchar();
}
while(isdigit(ch)){
x=(x<<1)+(x<<3)+ch-'0';
ch=getchar();
}
p=x*w;
}
//带权费马点问题，模拟退火解决.
int n,x[N],y[N],w[N];
double ansx,ansy,minn,minx,miny;
double dis(double x,double y){
return sqrt(x*x+y*y);
}
double getans(double nx,double ny){//求权值和
double ret=0;
rep(i,1,n)
ret+=dis(nx-x[i],ny-y[i])*w[i];
return ret;
}
void SA(){
double temp=1000;//温度
while(temp>1e-15){
double nwx=ansx+(rand()+rand()-RAND_MAX)*temp,nwy=ansy+(rand()+rand()-RAND_MAX)*temp,nw=getans(nwx,nwy);
if(nw<minn)minn=nw,ansx=nwx,ansy=nwy,minx=nwx,miny=nwy;//答案更优直接选
else if(exp(-(nw-minn)/temp)>(double)rand()/1.0/RAND_MAX){
ansx=nwx,ansy=nwy;//最优答案不变.
}
temp*=d;//降温
}
rep(i,1,100){
double nwx=minx+(rand()+rand()-RAND_MAX)*temp,nwy=miny+(rand()+rand()-RAND_MAX)*temp,nw=getans(nwx,nwy);
if(nw<minn)minn=nw,minx=nwx,miny=nwy,ansx=nwx,ansy=nwy;
}
}
int main(){
read(n);
double s=clock();
rep(i,1,n)
read(x[i]),read(y[i]),read(w[i]),ansx+=x[i],ansy+=y[i];
ansx=ansx/1.0/n,ansy=ansy/1.0/n,minn=getans(ansx,ansy),minx=ansx,miny=ansy;
while(1){
SA();
if((double)(clock()-s)/1.0/CLOCKS_PER_SEC>0.98)break;
}
printf("%.3lf %.3lf",minx,miny);
return 0;
}
在总共的 \(12\) 次提交中，该程序的正确率达到了 \(100\%\)。
不难发现每一侧的金币数量是固定的，只需要用数组存两个组的金币种类，每次交换两组中的任意两个金币，作为新的答案即可。需要注意的是，这里的随机范围无需随温度的减小而减小。
int getans(){//1的x和2的y换
int sum1=0,sum2=0;
rep(i,1,len)
sum1+=g[i][0];
rep(i,1,len)
sum2+=g[i][1];
return abs(sum1-sum2);
}
void SA(){
int temp=3000000;
while(temp>1e-5){
int x=rand()%len+1,y=rand()%len+1;
swap(g[y][1],g[x][0]);//存的分别是两个集合里的金币价值.
int nw=getans();
if(nw<minn)minn=nw;
else if(exp(-(nw-minn)/temp)<(double)rand()/1.0/RAND_MAX)swap(g[y][1],g[x][0]);//没在概率范围内要调回来
temp*=d;
}
rep(i,1,50){
int x=rand()%len+1,y=rand()%len+1;
swap(g[y][1],g[x][0]);
int nw=getans();
if(nw<minn)minn=nw;
else swap(g[y][1],g[x][0]);
}
}
调参不困难，一遍就 AC。在 \(3\) 次提交中，该代码的通过率达到了 \(100\%\)。
完结撒花！
浙公网安备 33010602011771号