一些求函数最值的算法
三分法
流程
三分法可以求区间上一个单峰函数的最值。顾名思义,三分法与二分法相似,但是会把区间分三段。
板子 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。
-
玄学。
[[杂项]]