简介

本算法用来求解凸函数极值点的问题，由我在写ACM习题时想到，在网上并未找到这样的算法，拿出来给大家分享一下，如果网上没有的话，我决定给它起名叫做 CSF迭代法，如果这个算法早已经存在，那就当看个笑话吧。

by 蔡少斐 西安交大 软件53

算法步骤

1. 确定一个常数$n$，表示区间的划分粒度，以及一个常数$eps$表示精度。

2. 先从定义域$[l,r]$中等差地取出n个点,$x_1,x_2,...,x_n$。

   $x_1 = l + \frac{(r-l)}{n+1},x_i = x_{i-1} +\frac{r-l}{n+1}$

3. 计算$f(x_1),...,f(x_n)$，挑选出使得$f(x_t)$最大的点$x_t$。

4. 判断定义域区间长度是否小于$eps$，若是执行$6$，否则执行$5$。

5. 把定义域缩小到$[x_{t-1},x_{t+1}]$，执行2。

   $x_0 = l,x_{n+1} = r$

6. 输出$x_t$

算法证明

用反证法：假设存在一次迭代，使得极值点不在$[x_{t-1},x_{t+1}]$里面，那么设这个极值点为$x_p$，并且有$x_p > x_{t+1}$，那么可以知道区间$[x_{t-1},x_{t+1}]$是单调递增的。也就是说$f(x_{t+1}) >f(x_t)$，因此我们在算法的第2个步骤时将选择$x_{t+1}$而不是$x_t$，矛盾。
因而极值一定在我们所迭代的区间内。

算法复杂度分析

设区间长度为$len$,划分粒度为$n$，精度要求为$eps$，递归深度为$dep$
可以列出递归方程如下：
$T(len) = n+T(\frac{len}{(n+1)/2})$
其中满足方程：
$(\frac{2}{n+1})^{dep}len =eps$
解得：
$O(len) =n log_{\frac{2}{n+1}}^{\frac{eps}{len}}$

特殊性

如果令划分粒度$n = 2$的话，等价于三分算法。

算法实现

double csf(double l,double r,int n = 2,double eps = 0.001){  static const double INF = 1e9;double x;while(r - l > eps){double step = (r-l)/(n+1);double mx = -INF;for(double i = l+step;i < r;i += step) if(mx < f(i)) mx = f(i),x = i;l = x-step;r = x+step;}return x;
} 

参数分析

函数总共有4个参数$l,r,n,eps$,其中有3个参数$l,r,eps$都是给定的。
总是选取$l = -10000000，r= 10000000$

当选取$eps = 0.001$时候

当选取$eps = 0.000001$的时候

结论

$n$选取$4$或6<script type="math/tex" id="MathJax-Element-37">6</script>的时候函数效果较好。

实验代码

#include <bits/stdc++.h>
using namespace std;
double f(double x){return -(x-5)*(x-786);
}
int cnt;
double csf(double l,double r,int n = 2,double eps = 0.001){  static const double INF = 1e9;double x;while(r - l > eps){double step = (r-l)/(n+1);double mx = -INF;for(double i = l+step;i < r;i += step) if(mx < f(i)) mx = f(i),x = i,cnt++;l = x-step;r = x+step;}return x;
} 
int main(){freopen("1.csv","w",stdout);for(int i = 2;i <= 100;i+=1) {cnt = 0;double x = csf(-10000000,10000000,i,0.000001);printf("%d,%d\n",i,cnt);}return 0;
}