如何找单峰函数峰值

CZQ_King

2019-07-31 21:10:00

算法·理论

本文同步发布与我的博客,欢迎来踩!

找单峰函数的峰值有什么用? 怎么找单峰函数的峰值?

别急,让这篇文章来告诉你。

1.什么是单峰函数?

肯定要先听听度娘的解释啦!

单峰函数是在所考虑的区间中只有一个严格局部极大值(峰值)的实值函数。如果函数f(x)在区间[a, b]上只有唯一的最大值点C,而在最大值点C的左侧,函数单调增加;在点C的右侧,函数单调减少,则称这个函数为区间[a, b]上的单峰函数。——百度百科

引用百度的原图大概的说一下:

就是函数f(x)有一个最高点x_0(如图),然后它的右边单调下降,左边单调上升。

也可以这样解释:

如果f(x)最大:

寻找一个单峰函数的峰值有什么用?

2.怎么找单峰函数的峰值?

(0)注意

本节中以f(x)为代码中的单峰函数。左端点和右端点为010

这里使用f(x)=\sqrt[x]{x}f(x)峰值为\sqrt[e]{e},约等于1.444667861

double f(double x){
    return pow(x,1/x);
}

(1)暴力枚举

最普通的方法肯定是暴力啦!

这个不用说了吧,设置步长,从左到右扫记录最大值。

优点:相对稳定。

缺点:耗时长,精度低。

#include<bits/stdc++.h>
using namespace std;
const int pre=2;//精度可调
int main(){
    double p=pow(10,-pre);//计算步长
    double l=0,r=10,maxx=-1;//看情况来定maxx的初始值
    for(double i=l+pre;i<=r;i+=p)//由于使用x^(1/x)为f函数,所以i不能为0
        maxx=max(maxx,f(i));//寻找最大值
    cout<<fixed<<setprecision(pre)<<maxx;
}

时间复杂度:O(10^n),照着这样,1秒只能算到7位数左右。

(2)三分查找1

像这种带有单调性的函数,让人不由自主的想起——二分。

但是我们不能只取一个点,因为这是无法得出峰值在哪里的,所以我们需要取两个点,通过比较这两个点来得出峰值的位置。

俩学家阮大佬曾在他的博客中介绍了三分查找,主要是分一半再分一半的思想(midlr的中点,mmidmidr的中点)。如图:

比较f(mid)f(mmid)的值:

证明在链接里有,这里就照着我的思路再讲一遍吧(其实是差不多的):

f(mid)>f(mmid)的值时,若峰值在midd右边,会导致有2个峰(否则f(mid)不会大于f(mmid)),所以峰值必在midd左边。

同样可以证第二个。

既然这样,我们就可以进行三分,每次可以舍去原长度的\frac{1}{2}\frac{1}{4}

那么时间复杂度最优O(2\log_2n),最坏O(2\log_\frac{4}{3}n)

另外,在链接里有提及到写法,一种是控制迭代次数,一种是直接控制精度,接下来的本节的代码都是直接控制精度

优点:速度快。

缺点:不稳定。

#include<bits/stdc++.h>
using namespace std;
const int pre=10;//设置精度
int main(){
    double p=pow(10,-pre);//计算精度限制
    double l=0,r=10,mid,mmid;//初始化
    while(r-l>p){
        mid=(l+r)/2;//计算mid
        mmid=(mid+r)/2;//计算mmid
        if(f(mid)-f(mmid)>=0)r=mmid;//大于时舍去右边
        if(f(mid)-f(mmid)<=0)l=mid;//小于时舍去左边
    }
    cout<<fixed<<setprecision(pre)<<max(f(l),f(r));//输出更大的
}

这里有一个问题:

代码中1011行中,照着思路来说,应该是

if(f(mid)-f(mmid)>0)r=mmid;
else l=mid;

为什么要换成

if(f(mid)-f(mmid)>=0)r=mmid;
if(f(mid)-f(mmid)<=0)l=mid;

因为我在测试的过程中,发现若换成带else的代码,那么当pre>15的时候,由于C++的精度问题会将\frac{x+r}{2}的值算成lx很接近l),那么就会算出mid=l,然后又算出mmid=l。那么mid=mmid就算出了f(mid)-f(mmid)等于0,那么就会执行else后的语句l=mid,然而mid本来就等于l,也就是说什么都没有变,最终导致了死循环。

另外,我们并没有提到过当f(mid)=f(midd)时的情况,因为在此时峰值三种情况都有可能(如图),所以只能靠“蒙”。

然而度娘告诉我们:我们应该“蒙”中间,也就是说把左右都去掉(如图):

于是就有了这样一段代码。

(3)三分查找2

为了让算法变得更稳定,我们可以把这两个点设在三等分处,这样每次都一定会舍弃原长度的\frac{1}{3}。时间复杂度O(2log_3n)

优点:速度快,相对更稳定。

#include<bits/stdc++.h>
using namespace std;
const int pre=10;
int main(){
    double p=pow(10,-pre);
    double l=0,r=10,m1,m2;
    while(r-l>p){
        m1=(l+r+l)/3;
        m2=(l+r+r)/3;//这两句需要注意!!!
        if(f(m1)-f(m2)>=0)r=m2;
        if(f(m1)-f(m2)<=0)l=m1;
    }
    cout<<fixed<<setprecision(pre)<<max(f(l),f(r));
}

在代码中有一处是需要注意的。

千万不要写成m1=(l+r)3和m2=(l+r)2/3!!

很多人认为二分时m=\frac{1}{2}(l+r),然而自然就想到三分是应该是\frac{1}{3}(l+r)\frac{2}{3}(l+r)

这样写是错误的,例如l=3,r=4的时候m1=\frac{7}{3}=2.333\dots,m2=\frac{14}{3}=4.666\dots。根本不在lr的范围内!

因为我们截取的是一段,所以原长度的\frac{1}{3}应该是\frac{r-l}{3},也就是说m1应该是l+\frac{r-l}{3}=\frac{2l+r}{3}m2=r-\frac{r-l}{3}=\frac{2r+l}{3}

那为什么二分是m=\frac{l+r}{2}

其实也是一样的道理,m=l+\frac{r-l}{2}=\frac{l+r}{2}

(4)“二分查找”

我们发现只要m1l的距离和m2的到r距离相等,算法就会稳定。如果我们想需要速度更快,应该怎么处理呢?

可以想到我们如果把m1m2取得接近于中点,那么一次比较就可以舍去几乎\frac{1}{2}的长度了!

如图中(C是一个较小的常数):

优点:速度比普通三分快

缺点:语言浮点运算误差必须小

#include<bits/stdc++.h>
using namespace std;
const int pre=10;
int main(){
    double p=pow(10,-pre);
    double l=0,r=10,m;
    while(r-l>=p){
        m=(l+r)/2;//取中点 
        if(f(m-p)-f(m+p)>=0)r=m;
        if(f(m-p)-f(m+p)<=0)l=m;//比较中点两端 
    }
    cout<<fixed<<setprecision(pre)<<max(f(l),f(r));
}

注意:如果常数C太小,那么C++就会将mid-Cmid+C看作相等,然后l=mid,r=mid,然后结束循环,最终你得到的答案只是f(mid)的值。

时间复杂度大约是O(\log_2n)

(5)求导+二分

为什么可以这样做?

如果你学过高中数学,你会知道导数的零点就是原函数的峰值。

因为是单峰,所以只需要找lr之间的使f(x)0的数。

那么问题就转化成了二分

做法:

首先,你需要先对函数求导。

个人不太支持这种,毕竟复杂的函数很难求导,如果是多项式求导就比较简单一些,直接用幂函数求导公式。

a_0+a_1x+a_2x^2+\dots+a_{n-1}x^{n-1}+a_nx^n

求导后是

a_1+2a_2x+3a_3x^2+\dots+na_nx^{n-1}

优点:用着舒服,基本不会出bug

缺点:复杂函数难求导

这里给出P3382 【模板】三分法的求导+二分的代码(\sqrt[x]{x}求导太难了其实是不会

因为懒,我这里直接放出小黑AWM大佬的题解。

#include<cstdio>
using namespace std;
int n;
double a[20],L,R,k;
double f(double x){
    double ans=0;
    for(int i=n;i>=1;i--) ans=ans*x+a[i];//常数项没了
    return ans;
}
int main(){
    scanf("%d%lf%lf",&n,&L,&R);
    for(int i=n;i>=0;i--)scanf("%lf",&a[i]),a[i]*=i;//a[i]*=i求导
    while(R-L>=1e-6){//二分
        double mid=(R+L)/2;
        f(mid)>0?L=mid:R=mid;//等价于if(f(mid)>0)L=mid;else R=mid;
    }
    printf("%.5lf\n",L);//输出答案
    return 0;
}

(6)0.618

然而这些做法都不能很好地优化。(除了求导二分)

最好的做法是循环利用上一次比较时剩下的点。使一个点被去掉后另一个点正好是下一次比较的其中一个点。这就是0.618法,是优选法的其中一种。

那么这两个点应该取在哪里呢?

![](https://s2.ax1x.com/2019/07/26/enWnF1.png) $\frac{a}{b}$约等于$0.618$。 那么$a=\frac{r-l}{1+\frac{1+\sqrt{5}}{2}}$,也就是是$\frac{r-l}{\phi+1}$。 这样,我们就确定了两个点的位置了。 **优点:重复利用测试点,速度极快。** ```cpp #include<bits/stdc++.h> #define F 0.38196601125010515179541316563436//1-φ using namespace std; const int pre=10; int main(){ double p=pow(10,-pre); double l=0,r=10,m1=l+F*(r-l),m2=r-F*(r-l); while(r-l>p){ if(f(m1)>=f(m2)){ r=m2;//去掉右边 m2=m1;//替换 m1=l+F*(r-l);//计算新的m1 if(f(m1)>f(m2))continue;//如果等于先不要结束 } if(f(m1)<=f(m2)){ l=m1;//去掉左边 m1=m2;//替换 m2=r-F*(r-l);//计算新的m2 } } cout<<fixed<<setprecision(pre)<<f(l); } ``` ## $(7)$本节小结 ### $I.$哪种算法更好? 应该是$0.618$法了。 如果说函数比较方便求导,那肯定推荐求导二分,毕竟只需要$log_2n$的时间。 附带一张图关于以上的算法的速度(没有写求导二分): $\color{orange}\text{暴力枚举}\color{black}\text{、}\color{purple}\text{三分查找1}\color{black}\text{、}\color{green}\text{三分查找2}\color{black}\text{、}\color{blue}\text{“二分查找”}\color{black}\text{、}\color{red}\text{0.618法}

Draw\ by\ GeoGebra

II.其他玄学的方法

直接看P3382的题解你就会发现有很多玄学大法:什么粒子群优化(by yuy_)啊、模拟退火(by Soulist)啊、梯度下降法(by 平衡树森林)啊、甚至还出现了四分(by Youngsc_AFO)。(顺便膜拜这些大佬)

III.小优化

若题目要求多项式的峰值,可以通过秦九韶算法简化多项式,从而使效率更高。

多项式a_0+a_1x+a_2x^2+a_3x^3+\dots+a_nx^n,若直接计算,则需要进行\frac{n(n+1)}{2}乘法和n次加法。

但是我们学过乘法分配律,你会发现除a_0项外每项都带有x,那么分配起来:a_0+x(a_1+a_2x+a_3x^2+\dots+a_nx^{n-1})

在括号内除a_1项外都含有x,再次分配:a_0+x(a_1+x(a_2+a_3x+\dots+a_nx^{n-2}))

以此类推,最后就会变成:a_0+x(a_1+\dots+x(a_{n-2}+x(a_{n-1}+a_nx)))。这个式子最终只需要进行n次乘法和n次加法。

int f(int x){
    int ans=0;
    for(int i=0;i<=n;i++)
        ans+=a[i]*pow(x,i);
    return ans;
}//普通算法
int f(int x){
    int ans=a[n];
    for(int i=n-1;i>=0;i--)ans=ans*x+a[i];
    return ans;
}//秦九韶算法

3.例题

(1)原创题目

题意:给你一张边长为n的正方形纸,在四个角落各剪一个边长为x的小正方形(0<x<\frac n2),拼成一个无盖长方体,请输出这个长方体体积最大时的体积、高和底面正方形的边长(保留5位小数)。

我们知道,它的体积是x(n-2x)^2,可以知道当x=0x=\frac n2时体积为0(还有一个x=\frac n2是重根),那么0\frac n2之间一定只有一个峰。类似下图:

所以可以用三分,三分x的长度(即三分盒子的高)。

代码:

#include<bits/stdc++.h>
using namespace std;
double n,l,r,m1,m2;
double f(double x){
    return x*(n-2*x)*(n-2*x);//返回体积
}
int main(){
    cin>>n;
    r=n/2.0;//题目说了x<n/2
    while(r-l>1e-7){
        m1=(l+r+l)/3;
        m2=(l+r+r)/3;
        if(f(m1)-f(m2)>=0)r=m2;
        if(f(m1)-f(m2)<=0)l=m1;
    }
    cout<<fixed<<setprecision(5)<<f(l)<<" "<<l<<" "<<sqrt(f(l)/l);//a*a*h=f(h)
}

当你测试几个测试点的时候,你就会发现答案似乎是一些循环小数。

确实是这样,这里是公式V=\frac{2n^3}{27},h=\frac n6,a=\frac{2n}{3},所以:

正解就有另一种解法:

#include<bits/stdc++.h>
using namespace std;
double n;
int main(){
    cin>>n;
    cout<<fixed<<setprecision(5)<<2*n*n*n/27<<" "<<n/6<<" "<<2*n/3;//直接套用公式
}

(2)POJ3737-UmBasketella

题意:输入圆锥表面积S,输出其最大体积,以及此时圆锥的高和底面半径。

先证明它是单峰函数(只看正数部分):

这只是其中的$S=1$的情况~~我也不会证只是给你们看一下知道是就行~~ 圆锥表面积公式:$S=\pi r(r+\sqrt{r^2+h^2})

圆锥体积公式:V=\frac{1}{3}\pi r^2h

知道这些后,就可以开始三分了。

做法:三分r的值。通过表面积Sr计算出h的值,再通过rh计算出圆锥的体积。

#include<iostream>
#include<algorithm>
#include<iomanip>
#include<math.h>//poj不给用万能头QAQ
#define M_PI acos(-1)
/*上面这句其实可以不加,因为math.h里有,但如果不加在poj就会CE*/
using namespace std;
double S,V,H;//表面积,体积,高
double f(double R){
    H=sqrt((S/(M_PI*R)-R)*(S/(M_PI*R)-R)-R*R);//通过表面积和底面半径算出高
    V=M_PI*R*R*H/3;//先存到V里
    return V;//返回体积
}
int main(){
    while(cin>>S){//有多个数据
        double l=0,r=sqrt(S/M_PI),m1,m2,t=50;//初始化
        /*因为表面积是πr(r+l)那么r大概就是sqrt(表面积/π)*/
        while(t--){//控制迭代次数,50次
            m1=(l+l+r)/3;
            m2=(l+r+r)/3;
            if(f(m1)<f(m2))l=m1;
            else r=m2;
            /*这部分基本一样*/
        }
        cout<<fixed<<setprecision(2)<<V<<endl<<H<<endl<<l<<endl;//照题意输出
    }
    return 0;
}

当然,你也可以三分h的值,只是通过Sh的值比较难求出r的值。

(3)本章小结

其实一般的三分题,只需要:

  1. 记住模板
  2. 选择要三分的数
  3. 修改f(x)函数
  4. 设置好左端点和右端点
  5. 调整精度

就可以轻松AC。如果你想做更多三分题,戳我。

4.结束

蒟蒻第一次写文章,可能还不够熟练,如果文章有什么问题欢迎各位大佬指正,如果有什么漏掉的也欢迎各位大佬补充!

参考: