浅谈随机化

· · 算法·理论

前言 - 我们为什么要学习随机化

随机化是一种极好的思想,当我们想不出正解的时候,我们就可以使用随机化。

并且,我们有的时候可以使用随机化过掉许多极难的题。

我们将会在这篇文章里,讲解随机化,由浅入深,一步步来。

从绿题到黑题,随机化所向披靡。

正文

Before we start - 了解我们要使用的工具

其实,随机化要使用的东西很少。

你看这下面的代码:

srand(k); // 其中 k 是任意一个正整数

就是我们的初始化操作。有的人会问了,为什么不使用

srand(time(0));

呢?因为 time(0) 具有一定的不确定性,在实践中可能会导致极其难调试。而且也不排除出题人会卡 time(0)

当我们想要调用一个随机数的时候,直接调用 rand() 就可以了。

而当我们想要打乱一个数组 a,其下标从 1n,我们可以使用这段代码:

random_shuffle(a + 1, a + n + 1);

简单的随机化

P2210 [USACO13OPEN] Haywire B

先读一下题再回来看。

这道题正解是状压 DP,但是我们可以用随机化 AC。

用一个数组 a 来表示奶牛的位置,每次随机打乱 a 数组并计算当前的值,答案为最小值。

#include <bits/stdc++.h>
using namespace std;
const int N = 12 + 5;
const int M = 5;
const int INF = 0x3f3f3f;
const double down = 0.996;
int n, a[N], fri[N][M];
signed main() {
    cin >> n;
    for (int i = 1; i <= n; i++) {
        a[i] = i;
    }
    for (int i = 1; i <= n; i++) {
        cin >> fri[i][1] >> fri[i][2] >> fri[i][3];
    }

    int temp = 2000000;
    int ans = INF;
    while (temp--) {
        random_shuffle(a + 1, a + n + 1);

        int tot = 0;
        for (int i = 1; i <= n; i++) {
            for(int j = 1; j <= 3; j++) {
                tot += abs(a[i] - a[fri[i][j]]);
            }
        }
        ans = min(ans, tot / 2);
    }

    cout << ans;
    return 0;
}

我们再给出另外一道题。

P5870 [SEERC 2018] Modern Djinn

此题正确题面。

使用随机化。

随机第 i 个人的心情为 mood_i,是 1 的时候心情好,0 的时候心情不好。代码就像下面这样:

for (int i = 1; i <= n; i++) {
    mood[i] = rand() & 1;
}

然后判断这时候可以实现多少个人的愿望:

int cnt = 0;
for (int i = 1; i <= m; i++) {
    if (mood[x[i]] && !mood[y[i]]) {
        ++cnt;
    }
}

cnt \geq \lfloor m/4 \rfloor +1 的时候输出答案:

if (cnt >= m / 4 + 1) {
    printf("%d\n" ,cnt);
    for (int i = 1; i <= m; i++) {
        if (mood[x[i]] && !mood[y[i]]) {
            printf("%d " , i);
        }
    }
    puts("");
    break;
}

但如果 cnt < \lfloor m/4 \rfloor +1,那么就重复之前的步骤,随机人的心情。

完整 AC 代码如下:

#include <bits/stdc++.h>
using namespace std;
const int N = 100000 + 5;
const int M = 200000 + 5;
int n, m;
int x[M], y[M];
bool mood[N];
int qcin()
{
    int t=0,f=1;
    char c=getchar(); 
    while(c<'0' || c>'9')
    {
        if(c=='-')
        {
            f=-1;
        }
        c=getchar();
    }
    while(c>='0' && c<='9')
    {
        t=t*10+c-'0';
        c=getchar();
    }
    return t*f;
}
signed main() {
    srand(72183);
    int T=qcin();
    while (T--) {
        n=qcin();
        m=qcin();
        for (int i = 1; i <= m; i++) {
            x[i]=qcin();
            y[i]=qcin();
        }
        while (true) {
            for (int i = 1; i <= n; i++) {
                mood[i] = rand() & 1;
            }
            int cnt = 0;
            for (int i = 1; i <= m; i++) {
                if (mood[x[i]] && !mood[y[i]]) {
                    ++cnt;
                }
            }
            if (cnt >= m / 4 + 1) {
                printf("%d\n" ,cnt);
                for (int i = 1; i <= m; i++) {
                    if (mood[x[i]] && !mood[y[i]]) {
                        printf("%d " , i);
                    }
                }
                puts("");
                break;
            }
        }
    }
    return 0;
}

爬山算法

有的时候,简单的随机化无法满足我们的需求,就需要一些特别的算法。这里先介绍一个随机化算法:爬山算法。

摘自 OI wiki:

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

兔子逐渐变得清醒的过程就是降温过程,即温度参数在爬山的时候会不断减小。

爬山算法只能适用在单峰函数上面。

这里给出一道例题:

P3382 三分

这题可以使用三分解决,但我更想要使用爬山算法来解决。

根据上面的思想,我们可以写出下面的代码:

#include <bits/stdc++.h>
using namespace std;
const double down=0.996;
int n;
double a[10005];
double dx,now,l,r;
double energy(double x)
{
    double tot=0;
    for(int i=0; i<=n; i++)
    {
        tot+=a[i]*pow(x,i);
    }
    return tot;
}
void hillclimb()
{
    double t=4000;
    while(t>1e-17)
    {
        double tx=dx+(rand()*2-RAND_MAX)*t;
        if(tx<l) tx=l;
        if(tx>r) tx=r;

        double temp=energy(tx);
        double delte=now-temp;
        if(delte<0)
        {
            dx=tx;
            now=temp;
          }
        t*=down;
    }
}
int main()
{
    srand(114514);
    cin>>n>>l>>r;
    for(int i=0; i<=n; i++)
    {
        cin>>a[n-i];
    }
    dx=(l+r)/2;
    now=energy(dx);
    hillclimb();

    cout<<fixed<<setprecision(7)<<dx<<endl;
    return 0;
} 

我们可以背下来这个模板,毕竟挺有用的。

有的时候,

另外给出一道例题,也是正解三分。

P1883 【模板】三分 / 函数 / [ICPC-Chengdu 2010] Error Curves

容易判断这是一个单峰函数,使用爬山算法。

代码如下:

#include <bits/stdc++.h>
using namespace std;
const double down=0.985;
int n;
struct node{
    int  x,y,w;
}a[10005];
double dx,now;
double energy(double x)
{
    double tot=-1e20;
    for(int i=1; i<=n; i++)
    {
        tot=max(tot,a[i].x*x*x+a[i].y*x+a[i].w);
    }
    return tot;
}
void hillclimb()
{
    double t=3000;
    while(t>1e-17)
    {
        double tx=dx+(rand()*2-RAND_MAX)*t;
        if(tx<0) tx=0;
        if(tx>1000) tx=(1000);

        double temp=energy(tx);
        double delte=temp-now;
        if(delte<0)
        {
            dx=tx;
            now=temp;
          }
        else if(exp(-delte/t)*RAND_MAX>rand())
        {
            dx=tx;
        }
        t*=down;
    }
}
void run()
{
    cin>>n;
    for(int i=1; i<=n; i++)
    {
        cin>>a[i].x>>a[i].y>>a[i].w;
    }
    dx=500;
    now=energy(dx);

    hillclimb();
    cout<<fixed<<setprecision(4)<<energy(dx)<<endl;
}
int main()
{
    int t;
    cin>>t;
    while(t--) run();
    return 0;
} 

为了确保我们可以获得更加正确的答案,我们可以改一下 run 函数,一直让爬山运行到差不多到时限:

void run()
{
    cin>>n;
    for(int i=1; i<=n; i++)
    {
        cin>>a[i].x>>a[i].y>>a[i].w;
    }
    dx=500;
    now=energy(dx);
    while(((double)clock()/CLOCKS_PER_SEC)<0.9)
        hillclimb();
    cout<<fixed<<setprecision(4)<<energy(dx)<<endl;
}

先到这里吧,因为我们要先介绍另外一个算法。

模拟退火

模拟退火,就是对上面的爬山算法进行了一点改进,从而可以适用于多峰函数。

这里给出一个视频,有助于理解。

模拟退火就是如果新状态的解更优就直接接受这个状态,否则以一定概率接受新状态。

这么说讲不清楚,结合一道经典例题来说吧。

P1337 [JSOI2004] 平衡点 / 吊打XXX

这道题爬山算法是这么写:

#include <bits/stdc++.h>
using namespace std;
const double down=0.996;
int n;
struct node{
    double x,y,w;
}a[1005];
double dx,dy,now;
double energy(double x,double y)
{
    double tot=0;
    for(int i=1; i<=n; i++)
    {
        tot+=sqrt((a[i].x-x)*(a[i].x-x)+(a[i].y-y)*(a[i].y-y))*a[i].w;
    }
    return tot;
}
void hillclimb()
{
    double t=4000;
    while(t>1e-15)
    {
        double tx=dx+(rand()*2-RAND_MAX)*t;
        double ty=dy+(rand()*2-RAND_MAX)*t;
        double temp=energy(tx,ty);
        double delte=temp-now;
        if(delte<0)
        {
            dx=tx;
            dy=ty;
            now=temp;

          }

        t*=down;
    }
}
int main()
{
    cin>>n;
    for(int i=1; i<=n; i++)
    {
        cin>>a[i].x>>a[i].y>>a[i].w;
        dx+=a[i].x;
        dy+=a[i].y;
    }
    dx/=n,dy/=n;
    now=energy(dx,dy);
    while(((double)clock()/CLOCKS_PER_SEC)<0.9) hillclimb();
    cout<<fixed<<setprecision(3)<<dx<<" "<<dy;
    return 0;
} 

这个代码只能拿 67 分,我们要使用模拟退火优化。其实模拟退火没有什么改变,就是在 hillclimb 做出了一点改变,但是为了区分算法,我们就把 hillclimb 改成 SA 吧(Simulated annealing 的英语缩写)。

SA 的代码如下:

void SA()
{
    double t=4000;
    while(t>1e-15)
    {
        double tx=dx+(rand()*2-RAND_MAX)*t;
        double ty=dy+(rand()*2-RAND_MAX)*t;
        double temp=energy(tx,ty);
        double delte=temp-now;
        if(delte<0)
        {
            dx=tx;
            dy=ty;
            now=temp;

          }
        else if(exp(-delte/t)*RAND_MAX>rand())
        {
            dx=tx;
            dy=ty;
        }

        t*=down;
    }
}

其中,能够体现模拟退火“以一定概率接受新状态”的部分就是:

else if(exp(-delte/t)*RAND_MAX>rand())
{
    dx=tx;
    dy=ty;
}

选择使用哪个算法

有的时候,我们需要使用模拟退火,但有的时候使用爬山算法,更有的时候简单的随机化就可以搞定。这需要我们的考虑。下面给出几个例题。

P2503 [HAOI2006] 均分数据

很明显,这道题数据小的可怕,一般来说简单的随机化就可以搞定。思路确定了,直接给出代码:

#include <bits/stdc++.h>
using namespace std;
const int N = 20 + 5;
const int INF = 0x3f3f3f;
const double down = 0.996;
int n, a[N], b[N], m;
signed main() {
    cin >> n >> m;
    double x_ = 0;
    for (int i = 1; i <= n; i++) {
        cin >> a[i];
        x_ += a[i];
    }
    x_ /= m;

    int temp = 500000;
    double ans = INF;
    while (temp--) {
        random_shuffle(a + 1, a + n + 1);
        memset(b, 0, sizeof(b));
        for (int i = 1; i <= n; i++) {
            int k = 1;
            for (int j = 2; j <= m; j++) {
                if (b[k] > b[j]) {
                    k = j;
                }
            }
            b[k] += a[i];
        }
        double tot = 0;
        for (int i = 1; i <= m; i++) {
            tot += (b[i] - x_) * (b[i] - x_);
        }
        tot /= m;
        if (tot < ans) {
            ans = tot;
        }
    }

    cout << fixed << setprecision(2) << sqrt(ans);
    return 0;
}

P3878 [TJOI2010] 分金币

这道题我们注意到这题没有任何单峰函数的迹象可言,又因为数据过于庞大,直接使用模拟退火。

又由于这题又多组数据不好卡时间,我们直接只让每一个数据的模拟退火运行 50 遍。

代码如下:

#include <bits/stdc++.h>
#define int long long
using namespace std;
const int N = 30 + 5;
const double down = 0.996;
const int INF = 2e10;
int n, a[N], now;
int energy() {
    double tot1 = 0, tot2 = 0;
    for (int i = 1; i <= (n + 1) / 2; i++) {
        tot1 += a[i];
    }
    for (int i = (n + 1) / 2 + 1; i <= n; i++) {
        tot2 += a[i];
    }

    return abs(tot1 - tot2);
}
void SA() {
    double t = 4000;
    while (t > 1e-15)
    {
        int tx=(rand()) % ((n + 1) / 2) + 1;
        int ty=(rand()) % ((n + 1) / 2) + ((n + 1) / 2);
        swap(a[tx], a[ty]);
        int temp = energy();
        int delte = temp - now;
        if (delte < 0) {
            now = temp;
        }
        else if (exp(-delte / t) * RAND_MAX < rand()) {
            swap(a[tx], a[ty]);
        }

        t *= down;
    }
}
void solve() {
    cin >> n;
    for (int i = 1; i <= n; i++) {
        cin >> a[i];
    }

    now = INF;

    int cs = 50;
    while (cs--) {
        SA();
    }

    cout << now << endl;
}
signed main() {
    int t;
    cin >> t;
    while(t--) {
        solve();
    }
    return 0;
}

P5544 [JSOI2016] 炸弹攻击1

这题很明显是一道模拟退火,直接按照模板写出代码来。

但是这样过不了 HACK,我们开始乱搞:

首先,我们可以先以 \left \{ \frac{\sum_{i=1}^{m} p_i}{m} ,\frac{\sum_{i=1}^{m} q_i}{m} \right \} 为初始点位进行一次模拟退火,然后分别将每一个敌人设为初始点位进行一次模拟退火,这样会让答案越来越接近最优解。但是如果只这样做了话,会 TLE 一堆点,得不偿失。

我们随便动一下参数,再进行去重,然后让设置每一个敌人为初始点位进行一次模拟退火的概率为 \frac{3}{4}。这样就可以过 HACK 了。

AC 代码:

#include <bits/stdc++.h>
#define int long long
using namespace std;
const int N = 10 + 5;
const int M = 1e3 + 5;
const double down = 0.999;
int n, m, R, ans;
struct peo{
    double x, y;
}b[M];
struct bui{
    double x, y, r;
}a[N];
int bst;
double tx, ty;
double dist(double x1, double y1, double x2, double y2) {
    return sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
}
int energy(double x, double y) {
    double r = R;
    for (int i = 1; i <= n; i++) {
        r = min(r, dist(x, y, a[i].x, a[i].y) - a[i].r);
    }
    int tot = 0;
    for (int i = 1; i <= m; i++) {
        if(dist(x, y, b[i].x, b[i].y) <= r) {
            tot++;
        }
    }
    return tot;
}
void SA() {
    double t = 2000;
    while (t > 1e-10) {
        double dx = tx + (rand() * 2 - RAND_MAX) * t;
        double dy = ty + (rand() * 2 - RAND_MAX) * t;
        int now = energy(dx, dy);
        double delta = bst - now;
        if (delta < 0) {
            tx = dx;
            ty = dy;
            bst = now;
            ans =max(ans, now);
        }
        else if (exp(delta / t) * RAND_MAX < rand()) {
            tx = dx;
            ty = dy;
            bst = now;
        }
        t *= down;
    }
}
map<pair<int, int>, bool> vis;
signed main() {
    srand(132052);
    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;
        tx += b[i].x;
        ty += b[i].y;
    }
    tx /= m;
    ty /= m;
    bst = max(bst, energy(tx, ty));
    SA();
    for (int i = 1; i <= m; i++) {
        tx = b[i].x;
        ty = b[i].y;
        if (((double)clock() / CLOCKS_PER_SEC) > 0.9) break;
        if(!vis[{tx, ty}] && rand() % 4 != 0) {
            SA();
            vis[{tx, ty}] = true;
        }

    }

    cout << ans;
    return 0;
}

P6505 Run Away

这道题很简单,既可以用退火写,又可以用爬山。虽然是一道黑题,但是如果我们看了上面的代码,改一改吊打 XXX 的代码就可以过。因为是黑题,留作习题。

习题

P7812 [JRKSJ R2] Dark Forest

P4044 [AHOI2014/JSOI2014] 保龄球

P6505 Run Away

结束语

我本来想说什么富有哲学的话,但是我没有能力。

随机化是一种很好的解题技巧,能够在我们想不出正解的时候骗部分分。

写这个真的是累死我了,看在这个份上,点个赞吧?

本文章总字数:12000 字