优化算法之D.S.C法与D.S.C-Powell法及其C++代码实现

本文主要是介绍优化算法之D.S.C法与D.S.C-Powell法及其C++代码实现，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

---

注：本文内容大多借鉴于《最优计算方法》——陈开周
注：本文讨论的优化函数只限于二维。

---

一、凸函数

在介绍这两个算法之前，我们需要先了解一下凸函数的概念。因为大多数传统的优化算法所求的函数模型都是默认为凸函数，这两个算法也不例外，而对于其他的复杂函数，目前比较流行的解决办法是使用进化算法，本文将只讨论D.S.C法与D.S.C-Powell法。
所谓凸函数，用一句话来概括就是弦值大于弧值的函数，也就是在一个函数上任取两点，这两点连成一条弦，然后在这两点对应的横坐标之间再任取一点，该点对应在弦上的值永远比对应在函数上的值大，下面的图很好的展示了这一特征：

如上图所示，在凸函数中，对于函数f(x)，任意a,b∈R，c∈(a,b)或(b,a)，g(x)是过(a,f(a))与(b,f(b))的一条直线，则g( c ) > f( c )。下面给出一些不是凸函数的示例：

二、二次插值法

对于我们要谈到的两种算法除了要了解凸函数外，我们还需要知道二次插值法，其是D.S.C与D.S.C-Powell法的一部分。二次插值法的核心思想是在函数上取三个点，这三个点的函数值需要满足两头大，中间小，然后设一二次方程，将这三个点代入方程得到一个方程组，如下：
$$\{\begin{array}{c}a{x}_1^2+b{x}_1+c=f(x_1)\\ a{x}_2^2+b{x}_2+c=f(x_2)\\ a{x}_3^2+b{x}_3+c=f(x_3)\end{array}$$
化为矩阵的形式，也就是：
$$\begin{array}{r}\left[\begin{array}{ccc}{x}_1^2& x_1& 1\\ x_2^2& x_2& 1\\ x_3^2& x_3& 1\end{array}\right]\left[\begin{array}{c}a\\ b\\ c\end{array}\right]=\left[\begin{array}{c}f(x_1)\\ f(x_2)\\ f(x_3)\end{array}\right]\end{array}$$
我们把上式看作 A * X=f( X)的形式（有点奇怪），由于 A 是一个范德蒙行列式，所以 X 必有解。我们之后通过解这个线性方程组就可以得到二次函数的系数，从而确定二次函数的最低点（$x^{\ast }=-\frac{b}{2a}$）,$x^{\ast }$就是我们之前取得三点所在范围内可能得最低点，通过不断插值，我们会越来越接近真实值，下面给出了求解的一个大致过程：

如上图所示，黑线是原函数f(x)，我们在该函数上取a,b,c三点，满足两头大，中间小的特点，然后通过这三个点拟合出一条二次函数曲线（蓝线），通过$x^{\ast }=-\frac{b}{2a}$，我们可以得到该二次函数的最低点，同时该点又对应原函数上一点d，然后我们将a,d,b作为下一次取点，不断重复这一过程求d点，最后我们会无限接近原函数的最小值点。现在的问题是有没有一个公式可以直接通过我们取得三个点得信息得出d点的坐标，下面给出推导过程和结论：
首先， 根据克罗姆法则我们可以直接求解得到，
$$\begin{array}{r}a=\left|\begin{array}{ccc}{f}_1& x_1& 1\\ f_2& x_2& 1\\ f_3& x_3& 1\end{array}\right|/\left|\begin{array}{ccc}{x}_1^2& x_1& 1\\ x_2^2& x_2& 1\\ x_3^2& x_3& 1\end{array}\right|\end{array}$$
$$\begin{array}{r}b=\left|\begin{array}{ccc}{x}_1^2& f_1& 1\\ x_2^2& f_2& 1\\ x_3^2& f_3& 1\end{array}\right|/\left|\begin{array}{ccc}{x}_1^2& x_1& 1\\ x_2^2& x_2& 1\\ x_3^2& x_3& 1\end{array}\right|\end{array}$$
$$\begin{array}{r}{x}^{\ast }=-\frac{b}{2a}=-\frac{1}{2}\left|\begin{array}{ccc}{x}_1^2& f_1& 1\\ x_2^2& f_2& 1\\ x_3^2& f_3& 1\end{array}\right|/\left|\begin{array}{ccc}{f}_1& x_1& 1\\ f_2& x_2& 1\\ f_3& x_3& 1\end{array}\right|\end{array}$$
计算行列式可得
$$\begin{array}{r}{x}^{\ast }=\frac{1}{2}\frac{(x_2^2-x_3^2)f_1+(x_3^2-x_1^2)f_2+(x_1^2-x_2^2)f_3}{(x_2-x_3)f_1+(x_3-x_1)f_2+(x_1-x_2)f_3}\end{array}\text{\hspace{1em}(2-1)}$$
仔细观察这个式子，除了容易记忆之外几乎没有什么优点，计算量较大，抛开加减法不说，可以统计出乘除法进行了14次。我们的目的是计算$-\frac{b}{2a}$，有没有其他的办法来计算这一比值呢，答案是肯定的。
用方程组第一行减去第二行可得：
$$\begin{array}{r}a(x_1^2-x_3^2)+b(x_1-x_3)=f_1-f_3\end{array}$$
$$\begin{array}{r}b=\frac{(f_1-f_3)-a(x_1^2-x_3^2)}{(x_1-x_3)}=-a(x_1+x_3)+\frac{(f_1-f_3)}{(x_1-x_3)}\end{array}$$
记 $c_1=\frac{(f_1-f_3)}{(x_1-x_3)}$，$b=-a(x_1+x_3)+c_1$，则$\frac{b}{a}=-(x_1+x_3)+\frac{c_1}{a}$
$$\begin{gathered} \begin{array}{r}a=\left|\begin{array}{ccc}{f}_1& x_1& 1\\ f_2& x_2& 1\\ f_3& x_3& 1\end{array}\right|/\left|\begin{array}{ccc}{x}_1^2& x_1& 1\\ x_2^2& x_2& 1\\ x_3^2& x_3& 1\end{array}\right|\end{array}=-\frac{(x_2-x_3)f_1+(x_3-x_1)f_2+(x_1-x_2)f_3}{(x_1-x_2)(x_2-x_3)(x_3-x_1)} \\ =-\frac{(x_2-x_1+x_1-x_3)f_1+(x_3-x_1)f_2+(x_1-x_2)f_3}{(x_1-x_2)(x_2-x_3)(x_3-x_1)} \\ =-\frac{(x_1-x_3)(f_1-f_2)+(x_2-x_1)(f_1-f_3)}{(x_1-x_2)(x_2-x_3)(x_3-x_1)} \\ =\frac{\frac{(f_1-f_2)}{(x_1-x_2)}+\frac{(f_1-f_3)}{(x_3-x_1)}}{(x_2-x_3)} \\ =\frac{\frac{(f_1-f_2)}{(x_1-x_2)}-c_1}{(x_2-x_3)} \end{gathered}$$
设$c_2=\frac{\frac{(f_1-f_2)}{(x_1-x_2)}-c_1}{(x_2-x_3)}$，综上可得：
$$\{\begin{array}{c}{c}_1=\frac{(f_1-f_3)}{(x_1-x_3)}\\ c_2=\frac{\frac{(f_1-f_2)}{(x_1-x_2)}-c_1}{(x_2-x_3)}\\ x^{\ast }=\frac{1}{2}[(x_1+x_3)-\frac{c_1}{c_2}]\end{array}\text{\hspace{1em}(2-2)}$$
可以很明显的看出，由(2-2)式求最小点位置只需做5次乘法，大大提高了计算效率。另外，若选择的三个点满足$x_2-x_1=d，x_3-x_2=d$条件时，式(2-1)与式(2-2)还可以进一步简化：
$$\begin{array}{r}{x}^{\ast }=x_2+\frac{(f_1-f_3)d}{f_1-2f_2+f_3}\end{array}\text{\hspace{1em}(2-3)}$$
$$\{\begin{array}{c}{c}_1=-\frac{(f_1-f_3)}{2d}\\ c_2=\frac{f_1-2f_2+f_3}{2d^2}\\ x^{\ast }=\frac{1}{2}[(x_1+x_3)-\frac{c_1}{c_2}]\end{array}\text{\hspace{1em}(2-4)}$$
在下面的D.S.C与D.S.C-Powell法中将会用到式(2-1)与(2-3)。

三、D.S.C法

D.S.C法是1964年Davies,Swann,Campey一起提出一种插值法，用于求解凸函数的最优解。这种方法是基于二次插值法设计的，同时借鉴了0.618分割法，其可以很快的找到凸函数的全局最优解。其过程如下：

1. 给定初值$x_0$与搜索步长d
2. 令k = 0，若$f(x_k+d)<=f(x_k)$，转第3步，否则d = 0 - d，转第3步
3. $x_{k+1}=x_k+d$，计算$f(x_{k+1})$
4. 若$f(x_{k+1})<=f(x_k)$，令d = 2 * d，k = k + 1，转第3步，否则令$d=d/2，x_m=x_{k+1}，x_{m+1}=x_m-d，x_{m-1}=x_k，x_{m-2}=x_{k-1}$
5. 在{$x_{m-2},x_{m-1},x_m,x_{m+1}$}中，若$f(x_{m-1})>f(x_{m+1})$，则丢掉$x_m$，令$x_b=x_{m-1}$，否则丢掉$x_{m-2}$，令$x_b=x_{m+1}$
6. 令$x_a=x_b-d，x_c=x_b+d$，做二次插值：$x^{\ast }=x_b+(f(x_a)-f(x_c)\ast d)/(2\ast (f(x_a)-2\ast f(x_b)+f(x_c)))$
7. 若$|x_a-x_b|<\epsilon$，停止迭代，$x^{\ast }$即为所求最优解，否则以$x^{\ast }$为起点继续搜索，$x_0=x^{\ast }，d=d/4$，转第2步

上述过程大致思路就是先找一个方向搜索直到函数值不再下降为止，取这次的点与前两次的点做二次插值得到一个新的最小点，将此点作为初始点重复上述过程，这样就可以逼近最小值点。

四、D.S.C-Powell法

这个方法与D.S.C法不同之处在于找到两头大中间小的点后的处理方式，D.S.C法将二次插值后的点设为初始点继续搜索，而D.S.C-Powell法通过重复二次插值逼近最小值点。

1. 给定初值$x_0$与搜索步长d
2. 令k = 0，若$f(x_k+d)<=f(x_k)$，转第3步，否则d = 0 - d，转第3步
3. $x_{k+1}=x_k+d$，计算$f(x_{k+1})$
4. 若$f(x_{k+1})<=f(x_k)$，令d = 2 * d，k = k + 1，转第3步，否则令$d=d/2，x_m=x_{k+1}，x_{m+1}=x_m-d，x_{m-1}=x_k，x_{m-2}=x_{k-1}$
5. 在{$x_{m-2},x_{m-1},x_m,x_{m+1}$}中，若$f(x_{m-1})>f(x_{m+1})$，则丢掉$x_m$，令$x_2=x_{m-1}$，否则丢掉$x_{m-2}$，令$x_2=x_{m+1}，x_1=x_2-d，x_3=x_2+d$
6. 做二次插值： $x^{\ast }=\frac{1}{2}\frac{(x_2^2-x_3^2)f(x_1)+(x_3^2-x_1^2)f(x_2)+(x_1^2-x_2^2)f(x_3)}{(x_2-x_3)f(x_1)+(x_3-x_1)f(x_2)+(x_1-x_2)f(x_3)}$
7. 若$|x_1-x_2|<\epsilon$，停止迭代，$x^{\ast }$即为所求最优解，否则以$x^{\ast }$作为$x_2$，然后将从原来的三点中找到最靠近$x^{\ast }$的两个点分别作为$x_1$与$x_3$，转第6步

目前来讲，此法相对于D.S.C来讲要更好些，下面将给出程序设计源码与实验结果。

五、实验结果

测试用例：
$f(x)=(x+4{)}^4+3(x+4{)}^3$

测试结果：

六、C++源码

#include <iostream>
#include <cmath>
#include <algorithm>
#include <ctime>using namespace std;//待求解函数
struct func{double operator()(double x) const{return  pow((x + 4), 4) + 3 * pow((x + 4), 3);}
};//D.S.C算法
class D_S_C{public:D_S_C(double start, double step):x0{start},delta{step}{iter = 0;}double run(func f, double precision);int get_iter() const {return iter;}double abs(double x){ return (x > 0)?x:(0 - x); }
private:double x0; //起点double delta; //步长double xk; //保存当前x的位置double xk1; //保存上一次x的位置double xk2; //保存上上一次x的位置int iter; //迭代次数};//D.S.C-Powell 算法
class D_S_C_Powell{
public:D_S_C_Powell(double start, double step):x0{start},delta{step}{iter = 0;}double run(func f, double precision);int get_iter() const {return iter;}double abs(double x){ return (x > 0)?x:(0 - x); }
private:double x0; //起点double delta; //步长double xk; //保存当前x的位置double xk1; //保存上一次x的位置double xk2; //保存上上一次x的位置int iter; //迭代次数
};int main(){D_S_C dsc(-10,0.0001);D_S_C_Powell dscp(-10,0.0001);auto t1 = clock();cout << "D.S.C Method result:" << endl;cout << "iteration : " << dsc.get_iter() << "\nMin point: "<< dsc.run(func(), 0.000001) << endl;auto t2 = clock();cout << "time : " << t2 - t1 << endl;t1 = clock();cout << "D.S.C-Powell Method result:" << endl;cout << "iteration : " << dscp.get_iter() << "\nMin point: "<< dscp.run(func(), 0.000001) << endl;t2 = clock();cout << "time : " << t2 - t1 << endl;return 0;}double D_S_C::run(func f, double precision){double fval = 0;double lfval = 0;xk = x0;do{if( f(xk + delta) > f(xk) ){ //判断该朝哪边寻找最低点delta = 0 - delta;}xk2 = xk;xk1 = xk;xk = xk + delta;lfval = f(xk1);fval = f(xk);while(fval <= lfval){ //找最低点delta = delta + delta; //步数加倍xk2 = xk1;xk1 = xk;xk = xk + delta;lfval = fval;fval = f(xk);}//退出循环时说明已近越过最低点delta = delta / 2;//从xk,xk-delta,xk1,xk2中选取三个点，其满足两头大，中间小if(f(xk - delta) > lfval){ //丢掉xkxk = xk - delta;}else{ //丢掉xk2xk2 = xk1;xk1 = xk - delta;}//做二次插值xk = xk1 + (f(xk2) - f(xk)) * delta / (2 * (f(xk2) - 2 * f(xk1) + f(xk)));delta = delta / 4; //缩短步距iter++;}while(abs(xk - xk2) >= precision); // 符合精度停止迭代return xk;}double D_S_C_Powell::run(func f, double precision){double fval = 0;double lfval = 0;xk = x0;if( f(xk + delta) > f(xk) ){ //判断该朝哪边寻找最低点delta = 0 - delta;}xk2 = xk;xk1 = xk;xk = xk + delta;lfval = f(xk1);fval = f(xk);while(fval <= lfval){ //找最低点delta = delta + delta; //步数加倍xk2 = xk1;xk1 = xk;xk = xk + delta;lfval = fval;fval = f(xk);}//退出循环时说明已近越过最低点delta = delta / 2;//从xk,xk-delta,xk1,xk2中选取三个点，其满足两头大，中间小if(f(xk - delta) > lfval){ //丢掉xkxk = xk - delta;}else{ //丢掉xk2xk2 = xk1;xk1 = xk - delta;}do{//做二次插值double fxk = f(xk);double fxk1 = f(xk1);double fxk2 = f(xk2);double a = (xk1*xk1-xk2*xk2)*fxk + (xk2*xk2-xk*xk)*fxk1 + (xk*xk-xk1*xk1)*fxk2;double b = (xk1-xk2)*fxk + (xk2-xk)*fxk1 + (xk-xk1)*fxk2;double xmin = 0.5 * a / b;if(xmin < xk1){ //丢掉xkxk = xk1;xk1 = xmin;}else{ //丢掉xk2xk2 = xk1;xk1 = xmin;}iter++;}while(abs(xk - xk2) >= precision && abs(xk - xk1) >= precision); // 符合精度停止迭代return xk;}

这篇关于优化算法之D.S.C法与D.S.C-Powell法及其C++代码实现的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！

---