最优化 | 一维搜索与方程求根 | C++实现

4.2 黄金分割法（0.618）

4.2.1 原理

0.618法（黄金分割法）是一种常用的不需要导数信息的一维搜索方法。假设 $f(x)$ 在闭区间 $[a,b]$ 内是单谷的，存在最小值点 $x^{\ast }\in (a,b)$ ， 这时可以用 $0.618$ 法求最小值点。

在区间 $[a,b]$ 内对称地取两个点 $c<d$ 在所谓的黄金分割点上， 即
$$\frac{d-a}{b-a}=\rho \triangleq \frac{\sqrt{5}-1}{2}\approx 0.618,\frac{b-c}{b-a}=\rho ,$$
黄金分割比例 $\rho$ 满足
$$\frac{1-\rho }{\rho }=\rho ,$$

• 这样的比例的特点是，当区间缩小到 $[a,d]$ 时， $c$ 仍为缩小后的区间的右侧黄金分割点；当区间缩小到 $[c,b]$ 时， $d$ 仍为缩小后的区间的左侧黄金分割点。

• 这样选了 $c,d$ 两个点后，比较 $f(c)$ 和 $f(d)$ 的值，如果 $f(c)$ 较 小，则因为 $a<c<d,f(a)>f(c),f(c)<f(d)$, 最小值点 $x^{\ast }$ 一定在区间 $[a,d]$ 内，可以把搜索区间缩 小到 $[a,d]$, 以 $[a,d]$ 作为含有最小值点的区间再比较其中两个黄金分割点的函数值，而且这时 $c$ 是 $[a,d]$ 中右 侧的黄金分割点，只要再添加左侧的黄金分割点就可以了。

• 如果比较 $f(c)$ 和 $f(d)$ 时发现 $f(d)$ 较小，则把搜 索区间缩小到 $[c,b]$, 这时 $d$ 是 $[c,b]$ 内的左侧黄金分割点，只要再添加右侧黄金分割点。每次迭代使得区间长 度缩小到原来的0.618倍, 缩小后的两个黄金分割点中一个的坐标和函数值是已经算过的。

设规定的最小值点绝对误差限为 $ϵ$ ，算法如下:

0.618法具有线性收敛速度。

$$0.618^n(b-a)\le ϵ\Rightarrow n\ge \frac{\ln \frac{ϵ}{b-a}}{\ln 0.618}$$

4.2.2 c++实现

求解函数$f(x)=x^2-5$迭代$n$次后的最小值区间。每次迭代区间长度减小到原来的0.618倍。

//
// Created by CHH3213 on 2022/9/12.
//
#include 
#include
#include 

using namespace std;
#define EPS 1e-6
/*
 * 使用黄金搜索法求解迭代n次后的含最小值区间，：f(x)=x^2-5
 */
double func(double x){
    /**
     * 构造函数,单谷函数
     * x:自变量
     * return：函数值
     */
    return x*x-5;
}

vector<double> golden_section_search(double a, double b,int n){
    /**
     *  含最小值的区间[a,b]
     * n: 迭代次数
     */
    double rho = (sqrt(5)-1)/2.0;//0.618黄金分割比
    double c=rho*a+(1-rho)*b;//(a, b)中靠近a的黄金分割点
    double d = rho*b+(1-rho)*a;//(a,b)中靠近b的黄金分割点
    while(n--){
        if(func(c)>func(d)){
            //这时(c, b)是含最小值区间, d是(c, b)的中靠近c的黄金分割点, 仍用a, b表示含最小值区间，c表示靠近a的黄金分割点
            a=c;
            c=d;
            d= rho*b+(1-rho)*a;//(a,b)中靠近b的黄金分割点
        }else{
            //这时(a, d)是含最小值区间,c是区间(a,d)的靠近d的黄金分割点,仍用a, b表示含最小值区间，d表示靠近b的黄金分割点
            b=d;
            d=c;
            c=rho*a+(1-rho)*b;
        }
    }
    return {a,b};
}

int main(){
    vector<double> res = golden_section_search(-2,2,10);
    cout<<res[0]<<' '<<res[1]<<endl;
    return 0;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50

迭代10次的结果为

-0.0263112 0.00621124
1

用黄金分割法求解最小值也很简单，只要迭代次数够多，或者设置$b-a$的精度达到预设的阈值即可。

进一步地，如果是求$f(x)=x^2-5=0$的根，那么可以构造损失函数，求解损失函数的最小值也就是求解方程的根。代码相比之前的改动很小：

//
// Created by CHH3213 on 2022/9/12.
//
#include 
#include
#include 

using namespace std;
#define EPS 1e-6
/*
 * 使用黄金搜索法求解迭代n次后的方程的根所在区间：f(x)=x^2-5=0
 */
double func(double x){
		//构造损失函数(x^2-5)^2
    return pow(x*x-5,2);
}

vector<double> golden_section_search(double a, double b,int n){
    /**
     *  含最小值的区间[a,b]
     * n: 迭代次数
     */
    double rho = (sqrt(5)-1)/2.0;//0.618黄金分割比
    double c=rho*a+(1-rho)*b;//(a, b)中靠近a的黄金分割点
    double d = rho*b+(1-rho)*a;//(a,b)中靠近b的黄金分割点
    while(n--){
        if(func(c)>func(d)){
            //这时(c, b)是含最小值区间, d是(c, b)的中靠近c的黄金分割点, 仍用a, b表示含最小值区间，c表示靠近a的黄金分割点
            a=c;
            c=d;
            d= rho*b+(1-rho)*a;//(a,b)中靠近b的黄金分割点
        }else{
            //这时(a, d)是含最小值区间,c是区间(a,d)的靠近d的黄金分割点,仍用a, b表示含最小值区间，d表示靠近b的黄金分割点
            b=d;
            d=c;
            c=rho*a+(1-rho)*b;
        }
    }
    return {a,b};
}

int main(){
    vector<double> res = golden_section_search(-3,3,20);
    cout<<res[0]<<' '<<res[1]<<endl;

    return 0;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

4.3 抛物线法

4.3.1 原理

为了用二次多项式逼近 $f(x)$ ，只 要有 $f(x)$ 函数图像上的三个点就可以。仍假设 $f(x)$ 在闭区间 $[a,b]$ 内是单谷的，存在最小值点 $x^{\ast }\in (a,b)$ 。 设 $c\in (a,b)$, 利用 $(a,f(a)),(c,f(c)),(b,f(b))$ 三个点求得二次插值多项式 $\phi (x)$ ，用 $\phi (x)$ 的最小值点 作为下一个近似值。

二次插值多项式知识点可以参考这篇博客。

设插值多项式为
$$\phi (x)={\beta }_0+{\beta }_{1}x+{\beta }_2{x}^2,$$
其中 ${\beta }_2>0$, 在给定了 $a<c<b$ 和相应函数值后，可以解出
β 2 = 1 b − a [ f ( b ) − f ( c ) b − c − f ( c ) − f ( a ) c − a ] , β 1 = f ( c ) − f ( a ) c − a − β 2 ( c + a ) ,

• 如果 $\tilde{x}<c$, 则以 $a,\tilde{x},c$ 为新的插值节点继续计算插值多项式并求插值多项式的最小值点；
• 如果 $\tilde{x}>c$, 则以 $c,\tilde{x},b$ 为新的插值节点继续计算插值多项式并求插值多项式的最小值点，如此迭代直至三个插值点足够接 近。

适当条件下抛物线法达到超线性收敛速度。

4.3.2 c++实现

求解函数$f(x)=x^2-5$的最小值点。

//
// Created by CHH3213 on 2022/9/12.
//
#include
#include 
/*
 * 使用抛物线法求解 x*x-5的最小值点
 */
using namespace std;

#define DELTA 1e-4
#define EPS 1e-6
double func(double x){
    /**
     * 构造函数
     * x:自变量
     * return：函数值
     */
    return x*x-5;
}

double quadratic_fit_search(double a,double b,double c,int n){
    /**
     * 抛物线法；
     * 区间(a,b)；
     * a
     while(n--){
         double beta_1 = 1/(b-a)*((func(b)-func(c))/(b-c)-(func(c)- func(a))/(c-a));
         double beta_2 = (func(c)-func(a))/(c-a)-beta_2*(c+a);
         double x = -beta_1/(2*beta_2);
         if(x<c){
             b=c;
             c=x;
         }else{
             a=c;
             c=x;
         }
     }
     return c;
}

int main(){
    cout<<quadratic_fit_search(-3,3,1,10)<<endl;
    return 0;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

4.3.3 改进

节 $\mathrm{4.3.1}$ 的算法 迭代不保证最小值点在拟合抛物线所用的三个点构成的区间内。下面是改进的做法。

对一元单谷函数 $f(x)$ ， 取三个点 $a<b<c$ 使得 $f(a)>f(b)$, $f(c)>f(b)$ 。

在迭代过程中要求每一步产生的新的 $a<b<c$ 都 满足 $f(a)>f(b),f(c)>f(b)$ ， 即最小值一定在 $(a,c)$ 内。我们可以直接利用Lagrange插值公式写出二次多项 式，并用导数等于零解出多项式的最小值点，可得从 $a,b,c$ 计算下一个点的迭代公式:
$$x^{\ast }=\frac{1}{2}\frac{y_a\left(b^2-c^2\right)+y_b\left(c^2-a^2\right)+y_c\left(a^2-b^2\right)}{y_a(b-c)+y_b(c-a)+y_c(a-b)}$$
产生 $x^{\ast }$ 后，利用 $a,b,c,x^{\ast }$ 这四个点中找到含最小值的缩小的区间并仍用 $a,b,c$ 三个点表示。

4.3.4 c++实现

求解函数$f(x)=x^2-5$的最小值点。

//
// Created by CHH3213 on 2022/9/12.
//
#include
#include 
/*
 * 使用抛物线法求解 x*x-5的最小值点
 */
using namespace std;

#define DELTA 1e-4
#define EPS 1e-6
double func(double x){
    /**
     * 构造函数
     * x:自变量
     * return：函数值
     */
    return x*x-5;
}

double quadratic_fit_search_2(double a,double b,double c,int n){
    /**
     * 抛物线法；
     * 区间(a,c)；
     * a
    while(n--){
        double y_a=func(a),y_b=func(b),y_c=func(c);
        double x = 0.5*(y_a*(b*b-c*c)+y_b*(c*c-a*a)+y_c*(a*a-b*b))/(y_a*(b-c)+y_b*(c-a)+y_c*(a-b));
        double y_x = func(x);
        if(x>b){
             if(y_x>y_b){   //用 a < b < x 作为新的起点
                 c=x;
             }else{     //用 b < x < c作为新的起点
                 a=b;
                 b=x;
             }
        }else if(x<b){
             if(y_x>y_b){   //用 x < b < c作为新的起点
                 a=x;
             }else{     //用 a < x < b作为新的起点
                 c=b;
                 b=x;
             }
        }

     }
    return b;

}

int main(){
    cout<<quadratic_fit_search_2(-2,2,3,50)<<endl;
    return 0;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58

同样地，如果要求方程的根，构造损失函数，求损失函数的最小值点就是在求方程的根。