P1337 [JSOI2004] 平衡点 / 吊打XXX

saadzin · 2026-10-11 15:23:18 · 个人记录

提供一种 CEM 的解法，在此题中仅仅需要 129ms 并且如果通过合理的调参可以做到更小。

https://www.luogu.com.cn/record/302448513

什么是 CEM

CEM（Cross-Entropy Method，交叉熵方法）是一种基于概率分布的随机优化算法，其核心思想是：通过不断采样、筛选优质解，并利用优质样本更新概率分布，使后续采样逐渐集中于最优解附近。

为什么选择 CEM

我们可以考虑最小化物品到绳结的最小加权距离，证明如下：

设绳结位置为 (x,y)，第 i 个固定点坐标为 (x_i,y_i)，对应重物的重力为 w_i。

由于绳子长度不变，绳结到固定点的距离增加时，对应重物的悬挂长度减少，重力势能随之增加。因此，整个系统的重力势能与各固定点到绳结的加权距离之和仅相差一个常数。

根据最小势能原理，系统的稳定平衡位置对应势能的最小值，所以可以将原问题转化为：

\min_{x,y}\sum_{i=1}^{n}w_i\sqrt{(x-x_i)^2+(y-y_i)^2}

此外，欧氏距离是凸函数，非负加权求和仍然具有凸性，因此该目标函数不存在非全局的局部极小值。

综上，最小化加权距离之和能够正确求得绳结的稳定平衡位置。

由此可见，本题目标函数为加权欧氏距离之和，具有凸性，不存在非全局的局部极小值。

因此，我们可以采用 Gaussian CEM（高斯交叉熵方法）进行求解。该算法通过正态分布生成候选解，筛选目标函数值较小的精英样本，并利用这些样本更新分布的均值和标准差，使搜索范围逐渐向最优区域集中。

由于本题只有两个决策变量，采用独立的二维正态分布即可完成采样与参数更新，实现简单，且无需计算梯度。

需要注意的是，Gaussian CEM 属于随机优化算法，不保证有限次迭代必然得到全局最优解，但在本题的低维凸优化结构下，通常能够取得较好的收敛效果。

Gaussian CEM 的大致思路如下：

初始化分布：分别为横、纵坐标设置正态分布的均值 μ 和标准差 σ。
随机采样：每轮生成 N 个候选位置，计算各位置对应的目标函数值。
精英筛选（Elite Selection）：将候选解按照目标函数值升序排序，按照 \rho 的精英解比例保留前 $ K=\lfloor \rho N\rfloor
分布更新：根据精英样本的均值和标准差，重新估计下一轮的高斯分布参数。
迭代收敛：重复上述过程，直到达到最大迭代次数或满足停止条件，同时维护历史最优解。

换句话说，Gaussian CEM 的核心思想就是通过不断迭代更新正态分布的参数，使采样逐渐集中于优质解所在的区域。

具体而言，每轮按照当前正态分布随机生成一批候选位置，计算对应的目标函数值，从中筛选出较优的部分作为精英样本，再重新拟合正态分布的均值与标准差，使下一轮采样更倾向于当前发现的优质解区域。

由于本题的目标函数为凸函数，不存在非全局的局部极小值，因此我们不需要专门引入均匀分布采样等全局探索机制来逃离局部极小值，可以直接使用 Gaussian CEM 进行迭代优化。

虽然目标函数具有凸性，但有限样本 CEM 仍可能因方差过早收缩而停滞，因此应合理设置采样数量、精英比例和终止条件，并维护历史最优解。

为使下一轮采样分布更倾向于产生优质解，我们采用最大似然估计拟合当前的精英样本。由于采样分布服从高斯分布，其最大似然估计具有解析解，即精英样本的均值与方差，因此可以直接据此更新采样分布参数。

最大似然估计的拟合公式如下：

\mu_x=\frac{1}{K}\sum_{i=1}^{K}x_i, \mu_y=\frac{1}{K}\sum_{i=1}^{K}y_i

\sigma_x^2=\frac{1}{K}\sum_{i=1}^{K}(x_i-\mu_x)^2,\sigma_y^2=\frac{1}{K}\sum_{i=1}^{K}(y_i-\mu_y)^2

#include<bits/stdc++.h>
using namespace std;
#define int long long
int n,st=0;
double X[1005],Y[1005],w[1005],minx=10001.1,maxx=-10001.1,miny=10001.1,maxy=-10001.1;
double ans=1e10,ansx,ansy;
struct ret{
    double ans,x,y;
    ret(double ans,double x,double y): ans(ans),x(x),y(y) {}
    bool operator < (const ret &a) const{
        return ans<a.ans;
    }
};
vector<ret> d;
double f(double x,double y){
    double sum=0;
    for(int i=1;i<=n;i++){
        sum+=sqrt((X[i]-x)*(X[i]-x)+(Y[i]-y)*(Y[i]-y))*w[i];
    }
    return sum;
}
signed main(){
    cin>>n;
    for(int i=1;i<=n;i++){
        cin>>X[i]>>Y[i]>>w[i];
        minx=min(minx,X[i]),maxx=max(maxx,X[i]);
        miny=min(miny,Y[i]),maxy=max(maxy,Y[i]);
    }
    mt19937 gen(time(0));
    double nx=(minx+maxx)/2,sx=(maxx-minx)/2,ny=(miny+maxy)/2,sy=(maxy-miny)/2;
    for(int i=1;i<=1000;i++){
        normal_distribution<> disx(nx,sx);
        normal_distribution<> disy(ny,sy);
        d.clear();
        for(int j=1;j<=200;j++){
            double x=disx(gen),y=disy(gen);
            d.push_back(ret(f(x,y),x,y));
        }
        sort(d.begin(),d.end());
        int lim=d.size()*0.2;
        if(d[0].ans<ans){
            ans=d[0].ans;
            ansx=d[0].x;
            ansy=d[0].y;
            st=0;
        }
        else st++;
        double totx=0,toty=0;
        for(int i=0;i<lim;i++){
            totx+=d[i].x,toty+=d[i].y;
        }
        nx=totx/lim,ny=toty/lim;
        totx=0,toty=0;
        for(int i=0;i<lim;i++){
            totx+=(nx-d[i].x)*(nx-d[i].x);
            toty+=(ny-d[i].y)*(ny-d[i].y);
        }
        sx=totx/lim,sy=toty/lim;
        sx=sqrt(sx),sy=sqrt(sy);
        if(st>=20) break;
    }
    printf("%.3f %.3f",ansx,ansy);
    return 0;
}