Eternity

Until The End

随机化算法

发表于 分类于 ，
本文字数： 14k 阅读时长 ≈ 13 分钟

“人，在不断地碰壁中，黑暗里摸索出光明。”

这是一个大概率不会使你 $\textcolor{Green}{\text{Accepted}}$ 的算法，但其适用范围就和暴力一样，非常广阔。而且，得分比暴力高。

适用范围：记答案函数为 $f(x)$ 在一个 $K$ 维空间里，值域为 $R$ ，则暴力的时间复杂度为 $R^K$ ，而我们将整个答案区间覆盖为线，面，或更高维的连续性问题，那么，答案如果是这个连续性区域的某一个极值（最极值）的话，模拟退火会在一个优秀的时空中将其找出，但是，会出现不定的误差。满足这类性质的题有 $\text{Dp}$ 和计算几何。

模拟退火

退火是一种金属热处理工艺，指的是将金属缓慢加热到一定温度，保持足够时间，然后以适宜速度冷却。目的是降低硬度，改善切削加工性；消除残余应力，稳定尺寸，减少变形与裂纹倾向；细化晶粒，调整组织，消除组织缺陷。准确的说，退火是一种对材料的热处理工艺，包括金属材料、非金属材料。而且新材料的退火目的也与传统金属退火存在异同。

摘自《百度百科》。

如果新状态的解更优则修改答案，否则以一定概率接受新状态。

温度函数

设立一个温度（步长）$delta$ ，表示当前的随机点记为 $x$ ，则我们当前搜索的区间就是 $[x-delta,x+delta]$ ，每一次减小温度，使区间变窄并确定结点，此操作称为退火。

温度分为初始温度，终止温度，衰减系数。可以表示为 $T_0,T_E,k$ ，其中 $k$ 是一个 $(0,1)$ 的小数。

一般而言，衰减系数越大，找到最优解的概率越大，但超时的概率也更大。

随机迭代

每一次迭代过程中，在当前区间内随机一个点，并求出当前点的函数值，令目标函数是一种能量函数，设新点为 $x’$ ，原来的点为 $x$ ，函数值为 $f(x)$ ，当前温度为 $T$ ，假设我们现在要求这个函数的最小值，有 $\Delta E=f(x’)-f(x)$ ，则有两种情况：

  • $\Delta E<0$ ，则跳到新点；
  • $\Delta E>0$ ，则以一定概率跳到新点。

跳过去的概率可以这样计算：

这个 $e$ 表示当前函数值。

那么，模拟退火就是重复这个过程，并逐渐降低温度以减小区间和到达终止温度以停止程序。但一般而言，在整个过程之后，也可以在得到最终解之后以最终解为原始点再次退火，以保证最终解的周围不存在更优解。

所以，一般都不会只做一次模拟退火。只要在时限之内，尽可能多随机几次。

一张来自 $\text{OI-Wiki}$ 的经典图：

可见，随着温度降低，取值变得稳定。

例题

P1337

据说这道题是求 $n$ 个带权类费马点，可以用三分套三分求解。

先给出模拟退火的伪代码板子：

很好，非常不清晰的伪代码，果然还是写不成 $\text{Sukwants}$ 那样呢。

实测，一次模拟退火能够得到的分数在 $66pt\sim 89pt$ ，三次模拟退火可能会超时，两次模拟退火可能合适。分数在 $89pts$ 之间波动，$9$ 次 $\textcolor{Green}{\text{AC}}$ 。

Probably AC Code 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
#include<bits/stdc++.h>
#define re register
typedef long long ll;
template<class T>
inline void read(T &x)
{
    x=0;
    char ch=getchar(),t=0;
    while(ch<'0'||ch>'9') t|=ch=='-',ch=getchar();
    while(ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
    if(t) x=-x;
}
template<class T,class ...T1>
inline void read(T &x,T1 &...x1)
{
    read(x),read(x1...);
}
template<class T>
inline void write(T x)
{
    if(x<0) putchar('-'),x=-x;
    if(x>9) write(x/10);
    putchar(x%10+'0');
}
template<>
inline void write(char c)
{
    putchar(c);
}
template<>
inline void write(char *s)
{
    while(*s) putchar(*s++);
}
template<class T,class ...T1>
inline void write(T x,T1 ...x1)
{
    write(x),write(x1...);
}
const int MAXN=1e4+10;
const double T0=1e6,TE=0.001,k=0.999;
int N,X[MAXN],Y[MAXN],Val[MAXN];
double ansX,ansY,Dist;
inline double randSeed()
{
    return (double) rand()/RAND_MAX;
}
inline double calc(double x,double y)
{
    double res=0;
    for(int i=1;i<=N;++i)
    {
        double dx=X[i]-x,dy=Y[i]-y;
        res+=std::sqrt(dx*dx+dy*dy)*Val[i];
    }
    if(res<Dist) Dist=res,ansX=x,ansY=y;
    return res;
}
inline void simuilateAnneal()
{
    double T=T0;
    double nowX=ansX,nowY=ansY;
    while(T>TE)
    {
        double nxtX=nowX+T*(randSeed()*2-1);
        double nxtY=nowY+T*(randSeed()*2-1);
        double delta=calc(nxtX,nxtY)-calc(nowX,nowY);
        if(exp(-delta/T)>randSeed()) nowX=nxtX,nowY=nxtY;
        T*=k;
    }
    for(int i=1;i<=1000;++i)
    {
        double nxtX=ansX+T*(randSeed()*2-1);
        double nxtY=ansY+T*(randSeed()*2-1);
        calc(nxtX,nxtY);
    }
}
int main()
{
    // freopen(".in","r",stdin);
    // freopen(".out","w",stdout);
    srand(time(NULL));				//一般而言，随机种子不能固定
    read(N);
    for(int i=1;i<=N;++i)
    {
        read(X[i],Y[i],Val[i]);
        ansX+=X[i],ansY+=Y[i];
    }
    ansX/=N,ansY/=N,Dist=calc(ansX,ansY);
    simuilateAnneal();
    simuilateAnneal();			//次数看玄学
    printf("%.3lf %.3lf",ansX,ansY);
    return 0;
}
/*
3
0 0 1
0 2 1
1 1 1
*/


不对的话，本人不负任何责任。

P2210

求出一个序列，使得题目所要求的代价最小。

那么，对于构造一个序列的题，就直接随机化就好了，每一次随机两个数并交换。

如果是纯暴力的话，时间复杂度是 $\mathcal O(n!)$ ，还是能过不少的，不过我们需要的时间复杂度是 $\mathcal O(rand())$ （bushi）

跑一遍模拟退火的模板，随机交换即可。随机种子套了四重玄学（本尊，神秘的东方力量，我保佑我自己，玄学玄学），然后居然一遍就过了。

Probably AC Code 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
#include<bits/stdc++.h>
#define gh() getchar()
#define re register
typedef long long ll;
template<class T>
inline void read(T &x)
{
    x=0;
    char ch=gh(),t=0;
    while(ch<'0'||ch>'9') t|=ch=='-',ch=gh();
    while(ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+(ch^48),ch=gh();
    if(t) x=-x;
}
template<class T,class ...T1>
inline void read(T &x,T1 &...x1)
{
    read(x),read(x1...);
}
template<class T>
inline void write(T x)
{
    if(x<0) putchar('-'),x=-x;
    if(x>9) write(x/10);
    putchar(x%10+'0');
}
template<>
inline void write(char c)
{
    putchar(c);
}
template<>
inline void write(char *s)
{
    while(*s) putchar(*s++);
}
template<class T,class ...T1>
inline void write(T x,T1 ...x1)
{
    write(x),write(x1...);
}
template<class T>
inline bool checkMax(T &x,T &y)
{
    return x<y?x=y,1:0;
}
template<class T>
inline bool checkMin(T &x,T &y)
{
    return x>y?x=y,1:0;
}
const int MAXN=15;
const int INF=0x3f3f3f3f;
const double T0=1e5,Te=0.0001,Tk=0.9999;
const double MAX_TIME=0.8;
int N,Pos[MAXN],ret=INF;
int Fri[MAXN][4];
inline int randSeed()
{
    return rand()%N+1;
}
inline int calc()
{
    int res=0;
    for(int i=1;i<=N;++i)
        for(int j=1;j<=3;++j)
            res+=abs(Pos[i]-Pos[Fri[i][j]]);
    return res;
}
inline void simulate()
{
    for(double T=T0;T>Te;T*=Tk)
    {
        int x=randSeed(),y=randSeed();
        if(x==y) continue;
        std::swap(Pos[x],Pos[y]);
        int nxtres=calc();
        if(nxtres<=ret) ret=nxtres;
        else if(exp((ret-nxtres)/T)>(double)rand()/RAND_MAX) std::swap(Pos[x],Pos[y]);
    }
}
int main()
{
    // freopen("simulate.in","r",stdin);
    // freopen("simulate.out","w",stdout);
    srand(time(NULL));
    srand(rand()+19260817);
    srand(rand()+20061112);
    srand(rand());
    read(N);
    for(int i=1;i<=N;++i)
    {
        for(int j=1;j<=3;++j)
            read(Fri[i][j]);
        Pos[i]=i;
    }
    ret=calc();
    while((double)clock()/CLOCKS_PER_SEC<MAX_TIME) simulate();
    write(ret>>1);
    return 0;
}
/*

*/

P2503

求均方差，随机改组即可，然后暴力统计答案。

Probably AC Code 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
#include<bits/stdc++.h>
#define gh() getchar()
#define re register
typedef long long ll;
template<class T>
inline void read(T &x)
{
    x=0;
    char ch=gh(),t=0;
    while(ch<'0'||ch>'9') t|=ch=='-',ch=gh();
    while(ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+(ch^48),ch=gh();
    if(t) x=-x;
}
template<class T,class ...T1>
inline void read(T &x,T1 &...x1)
{
    read(x),read(x1...);
}
template<class T>
inline void write(T x)
{
    if(x<0) putchar('-'),x=-x;
    if(x>9) write(x/10);
    putchar(x%10+'0');
}
template<>
inline void write(char c)
{
    putchar(c);
}
template<>
inline void write(char *s)
{
    while(*s) putchar(*s++);
}
template<class T,class ...T1>
inline void write(T x,T1 ...x1)
{
    write(x),write(x1...);
}
template<class T>
inline bool checkMax(T &x,T &y)
{
    return x<y?x=y,1:0;
}
template<class T>
inline bool checkMin(T &x,T &y)
{
    return x>y?x=y,1:0;
}
const int MAXN=22;
const int INF=0x3f3f3f3f;
const double T0=1e6,TE=1e-6,Tk=0.999;
const double MAX_TIME=0.5;
int N,M,Val[MAXN],Bel[MAXN];
double ans=INF,avg;
inline int randN()
{
    return rand()%N+1;
}
inline int randM()
{
    return rand()%M+1;
}
inline double calc()
{
    static double sum[MAXN];
    memset(sum,0,sizeof(sum));
    for(int i=1;i<=N;++i) sum[Bel[i]]+=Val[i];
    double res=0;
    for(int i=1;i<=M;++i) res+=(avg-sum[i])*(avg-sum[i]);
    res/=M;
    return res;
}
inline void simulate()
{
    std::random_shuffle(Val+1,Val+N+1);
    for(double T=T0;T>TE;T*=Tk)
    {
        int nowX=randN();
        int bel=Bel[nowX];
        int nxtX=randM();
        if(bel==nxtX) continue;
        Bel[nowX]=nxtX;
        double nowres=calc();
        if(nowres<=ans) ans=nowres;
        else if(exp((ans-nowres)/T)*RAND_MAX<(double)rand()) Bel[nowX]=bel;
    }
}
int main()
{
    // freopen("simulate.in","r",stdin);
    // freopen("simulate.out","w",stdout);
    srand(19260817);
    srand(rand()+114514);
    srand(rand()+20070722);
    srand(rand());
    read(N,M);
    for(int i=1;i<=N;++i) read(Val[i]),Bel[i]=randM(),avg+=Val[i];
    avg/=M;
    ans=calc();
    while((double)clock()/CLOCKS_PER_SEC<MAX_TIME) simulate();
    printf("%.2lf",std::sqrt(ans));
    return 0;
}
/*
6 3
1 2 3 4 5 6
*/

技巧

分块模拟退火

暴力怎么少得了分块呢（恼

单峰函数还好说，一下子就可以固定最优解，但如果是多峰函数且分布极不含规律的话，可以考虑值域分块，然后在每一段都跑一次模拟退火，然后择取最优解。

卡时

暴力算法，最忌讳的就是超时。

在 C++ 的某一个库（应该是 ctime）里，有一个函数 clock() 可以返回程序的运行时间。那么，我们就可以使在时间限制内不断执行模拟退火，使超时的概率和误差的概率都相应减少。

1
while((double)clock()/CLOCKS_PER_SEC<MAX_TIME) simulateAnneal();

其中的 MAX_TIME 是略小于该题目时间限制的数（单位为秒，准确而言 clock() 返回的是毫秒值，但我们加了一个 CLOCKS_PER_SEC 就变成秒了）。一般如果时限为 $1.00s$ 的话，可以定义 MAX_TIME=0.8 。

（实测这个优化比模拟退火本身更玄学）

玄学调参

如果时间很充裕，但最优解概率不大，可以选择调大 $T_0,k$ 并调小 $T_E$ 。

如果最优解稳定，但时间卡的有些紧，可以选择调小 $T_0,k$ 并调大 $T_E$ 。

一般而言，设定：

1
const double T0=1e6,TE=1e-6,Tk=0.993

对于一般非酋而言，也可以选择调整随机种子。

1
2
3
srand(time(NULL));
srand(rand()+19260817);
srand(rand());

多重祝福。

如果还是不行，考虑打正解吧。

打乱数列

对一些题目，数列的顺序与最优答案无影响的话，可以考虑使用 std::random_shuffle() 函数对整个数组进行打乱重组。需要注意的是，似乎这个函数已在 C++14 被废弃，C++17 被移除。但使用 C++17 提交又不会报错。有一定的风险。

AC Code 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
#include<bits/stdc++.h>
#define gh() getchar()
#define re register
typedef long long ll;
template<class T>
inline void read(T &x)
{
    x=0;
    char ch=gh(),t=0;
    while(ch<'0'||ch>'9') t|=ch=='-',ch=gh();
    while(ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+(ch^48),ch=gh();
    if(t) x=-x;
}
template<class T,class ...T1>
inline void read(T &x,T1 &...x1)
{
    read(x),read(x1...);
}
template<class T>
inline void write(T x)
{
    if(x<0) putchar('-'),x=-x;
    if(x>9) write(x/10);
    putchar(x%10+'0');
}
template<>
inline void write(char c)
{
    putchar(c);
}
template<>
inline void write(char *s)
{
    while(*s) putchar(*s++);
}
template<class T,class ...T1>
inline void write(T x,T1 ...x1)
{
    write(x),write(x1...);
}
template<class T>
inline bool checkMax(T &x,T &y)
{
    return x<y?x=y,1:0;
}
template<class T>
inline bool checkMin(T &x,T &y)
{
    return x>y?x=y,1:0;
}
const int MAXN=22;
const int INF=0x3f3f3f3f;
const double T0=1e6,TE=1e-6,Tk=0.999;
const double MAX_TIME=0.5;
int N,M,Val[MAXN],Bel[MAXN];
double ans=INF,avg;
inline int randN()
{
    return rand()%N+1;
}
inline int randM()
{
    return rand()%M+1;
}
inline double calc()
{
    static double sum[MAXN];
    memset(sum,0,sizeof(sum));
    for(int i=1;i<=N;++i) sum[Bel[i]]+=Val[i];
    double res=0;
    for(int i=1;i<=M;++i) res+=(avg-sum[i])*(avg-sum[i]);
    res/=M;
    return res;
}
inline void simulate()
{
    std::random_shuffle(Val+1,Val+N+1);
    for(double T=T0;T>TE;T*=Tk)
    {
        int nowX=randN();
        int bel=Bel[nowX];
        int nxtX=randM();
        if(bel==nxtX) continue;
        Bel[nowX]=nxtX;
        double nowres=calc();
        if(nowres<=ans) ans=nowres;
        else if(exp((ans-nowres)/T)*RAND_MAX<(double)rand()) Bel[nowX]=bel;
    }
}
int main()
{
    // freopen("simulate.in","r",stdin);
    // freopen("simulate.out","w",stdout);
    srand(19260817);
    srand(rand()+114514);
    srand(rand()+20070722);
    srand(rand());
    read(N,M);
    for(int i=1;i<=N;++i) read(Val[i]),Bel[i]=randM(),avg+=Val[i];
    avg/=M;
    ans=calc();
    while((double)clock()/CLOCKS_PER_SEC<MAX_TIME) simulate();
    printf("%.2lf",std::sqrt(ans));
    return 0;
}
/*
6 3
1 2 3 4 5 6
*/

爬山算法

解决单峰函数的随机化算法，可以部分替代三分。

爬山算法是一种局部择优的方法，采用启发式方法，是对深度优先搜索的一种改进，它利用反馈信息帮助生成解的决策。

因为三分只能解决二维的单峰问题，而如果是三维的单峰问题的话，就需要使用三分套三分。而如果维数极高的话，就需要用很多个三分连环套，套到吐。而且三分容易出错（就和二分一样），所以，不如爬山。

爬山算法的实现其实类似于模拟退火，也会引入一个温度函数，并进行降温。

类比地说，爬山算法就像是一只兔子喝醉了在山上跳，它每次都会朝着它所认为的更高的地方（这往往只是个不准确的趋势）跳，显然它有可能一次跳到山顶，也可能跳过头翻到对面去。不过没关系，兔子翻过去之后还会跳回来。显然这个过程很没有用，兔子永远都找不到出路，所以在这个过程中兔子冷静下来并在每次跳的时候更加谨慎，少跳一点，以到达合适的最优点。

截自 $\text{OI-Wiki}$ 。

况且，爬山算法只适用于解决单峰函数，并没有模拟退火优秀。

例题

P4035

这道题的正解应该是高斯消元，但如果对高斯消元不熟练的话，就可以考虑使用爬山算法，能够骗到比暴力高一点的分。

爬山算法其实并不随机，甚至连 rand() 都没用到。而且如果你交同一份代码，也不会像模拟退火那样随机 $\textcolor{red}{\text{Wrong Answer}}$ 几个点，只会在同样的地方犯同样的错误。

但是，你会发现，你调整 $T_0,T_E,k$ 的任意一个的任意一个值，都会导致结果大变。（这不是仅靠样例就能调出来的）所以，其实还是可以考虑模拟退火。

AC Code 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
#include<bits/stdc++.h>
#define gh() getchar()
#define re register
typedef long long ll;
template<class T>
inline void read(T &x)
{
    x=0;
    char ch=gh(),t=0;
    while(ch<'0'||ch>'9') t|=ch=='-',ch=gh();
    while(ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+(ch^48),ch=gh();
    if(t) x=-x;
}
template<class T,class ...T1>
inline void read(T &x,T1 &...x1)
{
    read(x),read(x1...);
}
template<class T>
inline void write(T x)
{
    if(x<0) putchar('-'),x=-x;
    if(x>9) write(x/10);
    putchar(x%10+'0');
}
template<>
inline void write(char c)
{
    putchar(c);
}
template<>
inline void write(char *s)
{
    while(*s) putchar(*s++);
}
template<class T,class ...T1>
inline void write(T x,T1 ...x1)
{
    write(x),write(x1...);
}
template<class T>
inline bool checkMax(T &x,T &y)
{
    return x<y?x=y,1:0;
}
template<class T>
inline bool checkMin(T &x,T &y)
{
    return x>y?x=y,1:0;
}
const int MAXN=15;
int N;
double Dir[MAXN][MAXN],ans[MAXN],Dist[MAXN],Delta[MAXN];
inline void calc()
{
    double avg=0;
    for(int i=0;i<=N;++i)
    {
        Dist[i]=Delta[i]=0;
        for(int j=0;j<N;++j)
            Dist[i]+=(Dir[i][j]-ans[j])*(Dir[i][j]-ans[j]);
        Dist[i]=std::sqrt(Dist[i]);
        avg+=Dist[i]/(N+1);
    }
    for(int i=0;i<=N;++i)
        for(int j=0;j<N;++j)
            Delta[j]+=(Dist[i]-avg)*(Dir[i][j]-ans[j])/avg;
}
int main()
{
    // freopen("simulate.in","r",stdin);
    // freopen("simulate.out","w",stdout);
    read(N);
    for(int i=0;i<=N;++i)
        for(int j=0;j<N;++j)
        {
            scanf("%lf",&Dir[i][j]);
            ans[j]+=Dir[i][j]/(N+1);
        }   
    for(double T=1e4;T>1e-6;T*=0.99995)
    {
        calc();
        for(int i=0;i<N;++i) ans[i]+=Delta[i]*T;
    }
    for(int i=0;i<N;++i) printf("%.3lf ",ans[i]);
    return 0;
}
/*
2
0.0 0.0
-1.0 1.0
1.0 0.0
*/
相关文章