OI学习笔记3：最优性算法
最优性算法
一、二分法
1、二分法。
（1）优势。
- 二分法可以将部分线性处理转化为对数处理。
- 二分法的时间复杂度通常为 \(O(\log N)\)。
（2）使用条件。
- 二分法仅能在问题状态空间严格单调的情况下使用。
（3）二分算法（详见OI学习笔记11：二分）。
2、三分法。
（1）函数的凸性。
- 上凸函数：对于函数 \(f(x)\) 上任意两点 \(x_1，x_2\)，令 \(x=q_1x_1+q_2x_2\)，其中 \(q_1+q_2=1\)，若有 \(f(x)\geq q_1f(x_1) +q_2f(x_2)\) 恒成立，则称 \(f(x)\) 为上凸函数。
- 下凸函数：对于函数 \(f(x)\) 上任意两点 \(x_1，x_2\)，令 \(x=q_1x_1+q_2x_2\)，其中 \(q_1+q_2=1\)，若有 \(f(x)\leq q_1f(x_1) +q_2f(x_2)\) 恒成立，则称 \(f(x)\) 为下凸函数。
- 琴生不等式：\(f(\sum q_ix_i)\leq \sum q_if(x_i)\)
（2）整数三分。
- 三分法与二分法的步骤相近，参照代码可简单理解。
while(l < r){
int lmid = l + (r - l) / 3;
int rmid = r - (r - l) / 3;
if(f(lmid) > f(rmid)) r = rmid;
else l = lmid;
}
（3）实数三分。
while(r - l > eps){//eps为精确度。
double lmid = l + (r - l) / 3;
double rmid = r - (r - l) / 3;
if(f(lmid) > f(rmid)) r = rmid;
else l = lmid;
}
二、启发式最优算法
1、爬山算法。
- 一种简单的贪心算法。
- 容易陷入局部最优解，无法搜出全局最优解。
2、模拟退火算法（SA）。
（1）原理。
- 来源于固体退火原理，是一种概率算法。
（2）算法流程。
- 初始化，给定初温 \(T\)，最低温度 \(T_{min}\)，同一温度下迭代次数 \(N\)，初始状态 \(x\)，初始化结果 \(result\)。
- 产生新状态 \(x'\)，产生新结果 \(result'\)。
- 设计合适的状态产生函数。函数要求具有全空间分散性或局部区域性。
- 避免状态的迂回搜索。
- 计算函数差值 \(\Delta f=f(x') - f(x)\)。
- 若 \(\Delta f \geq 0\)，则接受此新解作为当前解。
- 若 \(\Delta f < 0\)，以概率 \(P=\exp(\frac{-\Delta f}{k\cdot T})\) 大于 \([0,1)\) 区间某个随机数时接受此新解作为当前解。\(k\) 通常取 \(1\)。
- 对当前温度迭代 \(N\) 次。
- 若满足终止条件，输出当前解，结束算法。
- 否则降温，继续迭代。
（3）伪代码。
while(T > Tmin){
dE = J(Y(i + 1)) - J(Y(i));
if(dE >= 0)
Y(i+1) = Y(i);
else if(exp(dE / T) > random(0, 1)) Y(i + 1) = Y(i);
T *= r;//0 < r < 1。r过大，拉长搜索过程，增大搜到全局最优解可能。r过小，反之。
i++;
}//实际使用时根据时间限制和局部性调整参数，主要考虑T,r。
（4）例题。
POJ2420。
求平面上一个点，使其到给定的n个点的距离和最小，即费马点。
//AC代码
#include<cstdio>
#include<cmath>
#include<cstdlib>
using namespace std;
const double eps = 0.00000001;
const double PI = acos(-1.0);
struct Point{
double x, y;
Point(const double &x, const double &y){
this->x = x;
this->y = y;
}
Point(){}
inline void read(){
scanf("%lf%lf", &x, &y);
}
inline double length(){
return sqrt(x * x + y * y);
}
}a[105], p;
double ans;
int n;
typedef Point Vector;
Vector operator - (const Point &a, const Point &b){return Vector(a.x - b.x, a.y - b.y);}
Vector operator + (const Point &a, const Point &b){return Vector(a.x + b.x, a.y + b.y);}
Vector operator * (const double &k, const Point &b){return Vector(k * b.x, k * b.y);}
double calc(Point p){
double res = 0.0;
for(int i = 1; i <= n; i++){
res += (a[i]-p).length();
}
return res;
}
int main(){
srand(233);
scanf("%d", &n);
for(int i = 1; i <= n; i++){
a[i].read();
p.x += a[i].x;
p.y += a[i].y;
}
p.x /= (double)n;
p.y /= (double)n;
ans = calc(p);
double T = 100000.0;
while(T > eps){
double bestnow = 10000000.0;
Point bp;
for(int i = 1; i <= 10; i++){
double rad = (double)(rand()%10000 + 1) / 10000.0 * 2.0 * PI;
Point tp = p + T * Point(cos(rad), sin(rad));
double now = calc(tp);
if(now < bestnow){
bestnow = now;
bp = tp;
}
}
if(bestnow < ans || exp((ans - bestnow) / T) * 10000.0 > (double)(rand() % 10000)){
ans = bestnow;
p = bp;
}
T *= 0.99;
}
printf("%.0f\n", ans);
return 0;
}
3、遗传算法（GA）。
（1）原理：模拟生物种群进化过程，近似最优。
（2）编码:需将问题的解编码成字符串形式才能使用遗传算法，最简单的一种编码形式是二进制编码。
（3）适应度函数。
- 用于评价某个染色体的适应度，用 \(f(x)\) 表示。
- 适应度函数与目标函数是正相关的，目标函数变形后即为适应度函数。
- 有时需区分染色体的适应度函数与问题的目标函数。
（4）基本操作。
1.选择。
- 选择一些染色体产生下一代。
- 常用策略是比例选择。算法实现：
int RWS(){
m = 0; r = random(0, 1);
for(int i = 1; i <= N; i++){
m = m + P[i];
if(r <= m) return i;
}
}
- 交叉。
- 两条染色体交换部分基因，构造下一代的两条新染色体。
- 此行为以一定概率 \(P_c\) 发生。
- 变异。
- 新染色体中的基因会以一定概率出错。
- 变异发生概率记为 \(P_m\)。
（5）伪代码。
do{
计算种群Pop中每一个个体的适应度F(i);
初始化空种群newPop;
do{
根据适应度以比例选择算法从种群Pop中选出两个个体。
if(random(0, 1) < Pc)
对两个个体按Pc执行交叉操作。
if(random(0, 1) < Pm)
对两个个体按Pm执行变异操作。
将两个新个体加入种群newPop中。
}while(M个子代被创建);
用newPop取代Pop;
}while(任何染色体得分超过Tf,或繁殖代数超过G);
- 优化：精英主义选择：将每一代中的最优解原封不动地复制到下一代中去。
4、蚁群算法。
- 该算法理论尚不成熟，不做说明。
5、粒子群优化算法（PSO）。
（1）原理：鸟群觅食的仿生学模拟。
（2）算法流程：
- 初始化粒子群，给每个粒子分配随机位置与随机速度。
- 计算每个粒子的适应度。
- 根据适应度更新 \(pbest，gbest\)，更新粒子位置与速度。
- 迭代公式：
- \(v_{i+1}=v_i\cdot w+c\cdot rand()\cdot(pbest+gbest-2\cdot x_i)\)。
- \(w\) 为惯性因子，位于 \(0\) ~ \(1\) 之间。\(c\) 为学习因子，常取 \(2\)。
- \(x_{i+1}=x_i+v_i\)。如果出界，就将速度反转。
- 迭代公式：
- 判断是否达到迭代次数或解超过最优界限。
- 如果没有，返回第二步。
（3）例题。
给定一函数的各项系数，求其在区间 \([l,r]\) 的最值。
标程如下：
#include<bits/stdc++.h>
using namespace std;
const int cnt = 100;
int n;
double xs[15], l, r, by = -1e233, bx;
double f(double x){
double y = 0;
for(int i = n + 1; i >= 1; i --){
y += xs[n - i + 2] * pow(x, i - 1);
}
return y;
}
double Rand(){
return (double)rand()/RAND_MAX;
}
struct node{
double xv, x, y, besty, bestx;
}b[105];
inline void update(int a){
b[a].xv = b[a].xv * 0.5 + Rand() * 2 * (bx + b[a].bestx - b[a].x * 2);
b[a].x += b[a].xv;
if(b[a].x < l) b[a].x = l, b[a].xv *= -1;
if(b[a].x > r) b[a].x = r, b[a].xv *= -1;
b[a].y = f(b[a].x);
if(b[a].y > b[a].besty){
b[a].bestx = b[a].x;
b[a].besty = b[a].y;
}
}
int main(){
scanf("%d%lf%lf", &n, &l, &r);
for(int i = 1; i <= n + 1; i++){
scanf("%lf", &xs[i]);
}
srand(xs[1] + xs[n]);
for(int i = 1; i <= cnt; i++){
b[i].x = b[i].bestx = l + Rand() * (r - l);
b[i].xv = 0;
b[i].y = b[i].besty = f(b[i].x);
if(by < b[i].y){
bx = b[i].bestx;
by = b[i].besty;
}
}
for(int k = 1; k <= 100; k++){
for(int i = 1; i <= cnt; i++){
update(i);
if(by < b[i].besty){
bx = b[i].bestx;
by = b[i].besty;
}
}
}
printf("%.5lf\n", bx);
return 0;
}
三、牛顿迭代法
1、牛顿迭代法。
（1）几何定义。
- 选择曲线某点切线作为曲线的线性逼近（在切点附近，切线与曲线近似）。
- 像这样不断迭代，切线的根式解会不断收敛于曲线的某一个根式解。
（2）代数定义。
- 设 \(f(x)\) 为一曲线方程，由 \(f(x)\) 在 \(x_n\) 处的切线可以得到切线的根式解 \(x_{n+1}\)。
- \(\because f(x)\) 在 \(x_n\) 处的切线方程为 \(f'(x_n)(x-x_n)+f(x_n)=0\)，
\(\therefore x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}······(1)\)。
（3）从泰勒展开看。
- \(f(x)=\displaystyle \lim_{n \to \infty}{\sum_{i=0}^{n}{\frac{f^i(x_0)}{i!}}(x-x_0)^i+R_n(x)}\)。
- 可以看出 \(f(x)\) 的线性近似就是其一阶展开 \(\phi(x)=f(x_0)+f'(x_0)(x-x_0)\)。
（4）收敛性与时间复杂度。
- 收敛性：其收敛的充分条件为：\(f\) 二阶可导，只要初始值位于零点 \(x\) 的某一个去心邻域 \(\mathring{U}(x, \delta)\) 内，牛顿迭代法必定收敛。
- 时间复杂度为 \(O(\log n)\)。
（5）特例。
- 驻点：一阶导数为 \(0\) 的点，该点上 \((1)\) 式无效。
- 远离：如 \(y = x^{\frac{1}{3}}，x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}=-2x_n\)，切点对应的根点会远离。
- 循环震荡：如 \(y=|x|^{\frac{1}{2}}，x_{n+1}=-x_n\)，形成循环。
（6）求出所有根。
- 根据代数分解定理，复多项式 \(p(x)=a_0+a_1x+a_2x^2+...+a_nx^n\) 可被唯一分解为多个复单项式的乘积 \(p(x)=a(x-x_1)(x-x_2)...(x-x_n)\)。
- 根据该定理，每次逼近一个解 \(x_i\) 之后用多项式除法除去 \((x-x_i)\)，循环往复直至结果多项式为 \(0\) 时，就求出了所有解。
（7）使用牛顿迭代法的条件。
- 整个定义域内 \(f(x)\) 二阶可导。
- 要求多个根，需枚举不同初始值，避免陷入局部邻域。
（8）割线法。
- 用弦代替切线，用割线与横轴交点的横坐标作为近似。
- 由 \(x_n\) 和 \(x_{n-1}\) 确定的割线确定 \(x_{n+1}\)。
\(\: x_{n+1}=x_n-\frac{x_n-x_{n-1}}{f(x_n) - f(x_{n-1})}\)。 - 要两个初始值才能开始迭代。
2、应用
（1）给定 \(w\)，求 \(\sqrt{w}\)。
- 令 \(g(x)=\sqrt{x}\)，有 \(g^{-1}(\sqrt{x})=x\)，即 \(g^{-1}(x)=x^2\)。
- 设 \(f(x)=g^{-1}(x)-w=x^2-w\)，由于此函数二阶可导，所以迭代 \(f(x)\) 可得 \(\sqrt{w}\)。
- 该思路将开根运算转化为多次迭代的四则运算，适用于高精度。
- 代码：
#include<bits/stdc++.h>
using namespace std;
int w;
inline double f(double x){ return x * x - w; }
inline double _f(double x){ return 2 * x; }
inline double csqrt(){
double x = 2;
while(abs(f(x)) > 0.00001){
x = x - f(x)/_f(x);
}
return x;
}
int main(){
scanf("%d", &w);
printf("%.10lf", csqrt());
return 0;
}
（2）求 \(x^2\leq n\) 的最大整数解。
- \(x\)为整型，若迭代使近似解越界，则停止迭代，输出解。
（3）\(\exp\) 运算：给定 \(w\)，求 \(e^w\)。
- 令 \(g(x)=e^x\)，则 \(g^{-1}(e^x)=x\)，即 \(g^{-1}(x)=\ln x\)。
- 设 \(f(x)=g^{-1}(x)-w=\ln x-w\)，迭代 \(f(x)\) 即可。
（4）给定 \(g(x)=y\)，求 \(g(w)\)。
- 可化为求 \(f(x)=g^{-1}-w\)的零点。
- 反函数的导数=原函数导数的反函数。
（5）求函数极值。
- 可转化为对导函数迭代，因为导函数的零点就是原函数的极点。
- \(x_{n+1}=x_n-\frac{f'(x_n)}{f''(x_n)}\)。
四、分数规划
1、定义。
-
分数规划，又称 \(01\) 分数规划，用于在所有可行解中求出给定分式的极值。
-
\((S=\frac{\sum_{i=1}^{n}{a_i\cdot w_i}}{\sum_{i=1}^{n}{b_i\cdot w_i}})_{max/min}\)
-
具体而言，对于给定的 \(a_i，b_i\)，求一组 \(w_i\in [0,1]\)，最大/最小化分式 \(S\) 的值。
- 即，若将 \(a_i，b_i\) 看作一个物品的两个信息，则要求选出若干个物品使得 \(\frac{\sum{a}}{\sum{b}}\) 最大。
- 本质是一类最优比率问题，可以视作背包问题的变种。
- 有时还会给出一些额外限制，如：
- 分母大小不小于 \(W\)。
- \(w_i\) 为 \(1\) 的个数为 \(k\) 个。
