一些求函数最值的算法

三分法

流程

三分法可以求区间上一个单峰函数的最值。顾名思义，三分法与二分法相似，但是会把区间分三段。

板子 P3382

在区间 \([l,r]\) 中任取两点 \(lmid,rmid\) （设 \(lmid\leq rmid\)），分类讨论。

当 \(f(lmid)>f(rmid)\)，说明 \([lmid,rmid]\) 不可能单调增。那么 \(rmid\) 一定在峰值之后，排除区间 \([rmid,r]\)。

反之，则排除区间 \([l,lmid]\)。

为了排除区间更大，通常取 \(lmid=mid-eps,rmid=mid+eps\)，一次排除几乎一半，复杂度 \(O(\log n)\)。

代码：

#include<bits/stdc++.h>
using namespace std;
int n;
double a[18],l,r;
double f(double x,double ans=0.0,double xi=1.0){
  for(int i=0;i<=n;i++)ans+=xi*a[i],xi*=x;
  return ans;
}
int main(){
  cin>>n>>l>>r;
  for(int i=n;i>=0;i--)cin>>a[i];
  while(r-l>0.000001){
    double lmid=(l+r)/2-0.0000001,rmid=(l+r)/2+0.0000001;
    if(f(lmid)>f(rmid))r=rmid;
    else l=lmid;
  }
  return cout<<l<<'\n',0;
}

三分套三分

当函数超过一维，就需要三分套三分。

例：P2571

最短路径应该有三段：点 \(A\) 传送到一点 \(E\)，然后沿平面走到 \(CD\) 上一点 \(F\)，然后传送到 \(D\)。设 \(x=\frac{AE}{AB},y=\frac{DF}{CD}\)，得到二元函数 \(f(x,y)=\frac{x\cdot AB}{P}+\frac{dis(E,F)}{R}+\frac{y\cdot CD}{Q}\)。显然这个函数固定其中一个变量都是单峰的。

于是在外层三分 \(x\)，在内层把 \(x\) 作为参数字母三分 \(y\)。

代码：

#include<bits/stdc++.h>
using namespace std;
long double ax,ay,bx,by,cx,cy,dx,dy,o,p,q,ab,cd;
long double dis(long double ax,long double ay,long double bx,long double by){
  return sqrtl((ax-bx)*(ax-bx)+(ay-by)*(ay-by));
}
long double f(long double x,long double y){
  long double ae=ab*x,df=cd*y,ex=ax+(bx-ax)*x,ey=ay+(by-ay)*x,fx=dx+(cx-dx)*y,fy=dy+(cy-dy)*y;
  return ae/o+dis(ex,ey,fx,fy)/q+df/p;
}
long double solve(long double x){
  long double l=0,r=1;
  while(r-l>0.000001){
    long double lmid=(l+r)/2-0.0000001,rmid=(l+r)/2+0.0000001;
    if(f(x,lmid)<f(x,rmid))r=rmid;
    else l=lmid;
  }
  return f(x,l);
}
int main(){
  cin>>ax>>ay>>bx>>by>>cx>>cy>>dx>>dy>>o>>p>>q,ab=dis(ax,ay,bx,by),cd=dis(cx,cy,dx,dy);
  long double l=0,r=1;
  while(r-l>0.000001){
    long double lmid=(l+r)/2-0.0000001,rmid=(l+r)/2+0.0000001;
    if(solve(lmid)<solve(rmid))r=rmid;
    else l=lmid;
  }
  return cout<<setprecision(2)<<fixed<<solve(l)<<'\n',0;
}

爬山算法

流程

爬山算法可以求区间上一个函数最值，一般用于单峰函数。

设当前位置为 \(x\)，爬山算法每次会在 \(x\) 左右找两个点 \(x_1,x_2\)，与原答案比较。为了让答案固定，需要让变化幅度逐渐减小。引入温度参数 \(T_0\) 初始温度，\(T\) 当前温度，\(T_k\) 结束温度，降温系数 \(\Delta\)。\(T\) 从 \(T_0\) 开始降温，每次循环就乘 \(\Delta\)，而变化的幅度应当随着 \(T\) 变小而变小。

代码：

#include<bits/stdc++.h>
using namespace std;
int n;
double a[18],l,r;
double f(double x,double ans=0.0,double xi=1.0){
  for(int i=0;i<=n;i++)ans+=xi*a[i],xi*=x;
  return ans;
}
const double d=0.999,tk=1e-10;
double solve(double x){
  double t=10,ans=-1145141919;
  while(t>tk){
    if(x+t<=r&&f(x+t)>ans)ans=f(x+t),x+=t;
    else if(x-t>=l&&f(x-t)>ans)ans=f(x-t),x-=t;
    t*=d;
  }
  return x;
}
int main(){
  cin>>n>>l>>r;
  for(int i=n;i>=0;i--)cin>>a[i];
  return cout<<max(max(solve(l),solve(r)),solve((l+r)/2))<<'\n',0;
}

爬山算法可以容易的增加维度，对于每一维都尝试移动即可。

缺点显而易见，答案可能陷在局部最优解。

优化

为了水过多峰函数，可以挑不同的位置多爬几次取最优解，比如上面的代码取了三个初始位置。

另外作为玄学算法的一种，模拟退火的近亲，调参也是不可缺少的。不过可能没有模拟退火难调。

模拟退火

流程

当爬山算法到达一个局部最优解时，如果爬不出来，那么就会固定在那里，因此当函数不是单峰时，爬山算法是很不靠谱的。

模拟退火是一个著名的随机化骗分算法，常用于求解某个问题的最优解。

几个参数：\(T_0\) 初始温度，\(T\) 当前温度，\(T_k\) 结束温度，降温系数 \(\Delta\)，当前最优解 \(E_0\)，上一个被接受的解 \(E_1\)，新的解 \(E\)，新的解与原解的差 \(\Delta E\)。

首先将 \(T\) 设为 \(T_0\)，每次降温将 \(T\) 乘 \(\Delta\)，直到 \(T<T_k\)。同时，每次生成一个新的解 \(E\)。这个解应该由 \(E_1\) 随机变化而成，并且随着 \(T\) 变小而变小。这样解就会随着降温而稳定。

当新的解更优，即 \(\Delta E<0\)，则接受解，将 \(E_0,E_1\) 设为 \(E\)。否则根据 Metropolis 准则接受，有 \(e^{\frac{-\Delta E}{T}}\) 的概率接受，此时仅将 \(E_1\) 设为 \(E\)。其中，\(\Delta E\) 可以取 \(E-E_0\) 或 \(E-E_1\)。

以 UVA10228 为例：

求 \(n\) 边形的费马点，直接上模拟退火板子。

代码：

#include<bits/stdc++.h>
using namespace std;
mt19937 engine(chrono::steady_clock::now().time_since_epoch().count());
template<typename T>T randint(T l,T r){
  return uniform_int_distribution<T>(l,r)(engine);
}
template<typename T>T random(T l,T r){
  return uniform_real_distribution<T>(l,r)(engine);
}
chrono::steady_clock::time_point s=chrono::steady_clock::now();
long double get_time(){
  return chrono::duration<long double,milli>(chrono::steady_clock::now()-s).count();
}
int t,n;
long double x[1005],y[1005],ansx,ansy,ans=1e18;
long double calc(long double px,long double py,long double ans=0){
  for(int i=1;i<=n;i++)ans+=sqrtl((px-x[i])*(px-x[i])+(py-y[i])*(py-y[i]));
  return ans;
}
const long double t0=10000000,d=0.996,tk=1e-14;
void sa(){
  long double t=t0,tempx=ansx,tempy=ansy;
  while(t>tk){
    long double nowx=tempx+random(-32768.0,32768.0)*t,nowy=tempy+random(-32768.0,32768.0)*t,now=calc(nowx,nowy),delta=now-ans;
    if(delta<0)ansx=tempx=nowx,ansy=tempy=nowy,ans=now;
    else if(random(0.0,1.0)<exp(-delta/t))tempx=nowx,tempy=nowy;
    t*=d;
  }
}
int main(){
  cin>>t;
  while(t--){
    cin>>n;
    for(int i=1;i<=n;i++)cin>>x[i]>>y[i],ansx+=x[i],ansy+=y[i];
    ansx/=n,ansy/=n,ans=calc(ansx,ansy),sa(),cout<<(int)(ans+0.5)<<'\n';
    if(t)putchar('\n'); 
    ansx=ansy=0,ans=1e18,memset(x,0,sizeof(x)),memset(y,0,sizeof(y)); 
  }
  return 0;
}

小寄巧

模拟退火非常玄学，有时 WA 需要一些玄学的修改。

• 调参，修改 \(T_0,T_k,\Delta\)。这几个参数变更对结果影响很大，而且要结合具体题目。

• 在时间允许时多跑几遍模拟退火。

• 卡时，一直跑模拟退火直到即将 TLE。

• 玄学。

[[杂项]]