凸優化之無約束優化（一維搜索方法：二分法、牛頓法、割線法）

 

1、二分法（一階導）

二分法是利用目標函數的一階導數來連續壓縮區間的方法，因此這里除了要求 f 在 [a₀,b₀] 為單峰函數外，還要去 f(x) 連續可微。

（1）確定初始區間的中點 x⁽⁰⁾=(a₀+b₀)/2 。然后計算 f(x) 在 x⁽⁰⁾ 處的一階導數 f'(x⁽⁰⁾)， 如果 f'(x⁽⁰⁾) >0 , 說明極小點位於 x⁽⁰⁾ 的左側，也就是所，極小點所在的區間壓縮為[a₀,x⁽⁰⁾]；反之，如果 f'(x⁽⁰⁾) <0，說明極小點位於x⁽⁰⁾的右側，極小點所在的區間壓縮為[x⁽⁰⁾，b₀]；如果f'(x⁽⁰⁾) = 0，說明就是函數 f(x) 的極小點。

（2）根據新的區間構造x⁽¹⁾，以此來推，直到f'(x^(k)) = 0，停止。

可見經過N步迭代之后，整個區間的總壓縮比為（1/2）^N，這比黃金分割法和斐波那契數列法的總壓縮比要小。

#ifndef _BINARYSECTION_H_
#define _BINARYSECTION_H_

typedef float (* PtrOneVarFunc)(float x);
void BinarySectionMethod(float a, float b, PtrOneVarFunc fi, float epsilon);

#endif

#include<iostream>
#include<cmath>
#include "BinarySection.h"

using namespace std;

void BinarySectionMethod(float a, float b, PtrOneVarFunc tangent, float epsilon)
{
    float a0,b0,middle;
    int k;
    k = 1;
    a0 = a;
    b0 = b;
    middle = ( a0 + b0 )/2;
    
    while( abs(tangent(middle)) - epsilon > 0 )
    {
#ifdef _DEBUG
        cout<<k++<<"th iteration:x="<<middle<<",f'("<<middle<<")="<<tangent(middle)<<endl;
#endif
 
        if( tangent(middle) > 0)
        {
            b0 = middle;
        }
        else
        {
            a0 = middle;
        }
        middle =( a0+b0)/2;
    }
    
    cout<<k<<"th iteration:x="<<middle<<",f'("<<middle<<")="<<tangent(middle)<<endl;
}

#include<iostream>
#include "BinarySection.h"


float TangentFunctionofOneVariable(float x)
{
    return 14*x-5;//7*x*x-5*x+2;
}

int main()
{
    BinarySectionMethod(-50, 50, TangentFunctionofOneVariable, 0.001);
    return 0;
}

1th iteration:x=0,f'(0)=-5
2th iteration:x=25,f'(25)=345
3th iteration:x=12.5,f'(12.5)=170
4th iteration:x=6.25,f'(6.25)=82.5
5th iteration:x=3.125,f'(3.125)=38.75
6th iteration:x=1.5625,f'(1.5625)=16.875
7th iteration:x=0.78125,f'(0.78125)=5.9375
8th iteration:x=0.390625,f'(0.390625)=0.46875
9th iteration:x=0.195312,f'(0.195312)=-2.26562
10th iteration:x=0.292969,f'(0.292969)=-0.898438
11th iteration:x=0.341797,f'(0.341797)=-0.214844
12th iteration:x=0.366211,f'(0.366211)=0.126953
13th iteration:x=0.354004,f'(0.354004)=-0.0439453
14th iteration:x=0.360107,f'(0.360107)=0.0415039
15th iteration:x=0.357056,f'(0.357056)=-0.0012207
16th iteration:x=0.358582,f'(0.358582)=0.0201416
17th iteration:x=0.357819,f'(0.357819)=0.00946045
18th iteration:x=0.357437,f'(0.357437)=0.00411987
19th iteration:x=0.357246,f'(0.357246)=0.00144958
20th iteration:x=0.357151,f'(0.357151)=0.000114441

 

2、牛頓法（二階導）

前提：f 在 [a₀,b₀] 為單峰函數， 且[a₀,b₀] 在極小點附近，不能離的太遠否則可能無法收斂。

這里進一步要求f(x)連續二階可微。對於函數 f(x) 上一點 x^(k)，我們可以使用泰勒公式構造一個多項式函數 q(x)=f(x^(k))+f'(x^(k))(x-x^(k))+1/2 * f''(x^(k))(x-x^(k))²,對 f（x）在 x^(k) 附近進行局部二次擬合，q(x)可以看為f(x)的（在x^(k)附近的局部）近似（當f(x)為二次函數時，q(x)=f(x)），因此求f(x)的極小點可以轉化為求q(x)的極小點。

0=q'(x)=f'(x^(k))+f''(x^(k))(x-x^(k))  =>   x = x^(k)  -  f'(x^(k))/f''(x^(k))  因此可以選擇 x^(k+1) = x^(k)  -  f'(x^(k))/f''(x^(k))

牛頓法能夠不斷地迫使目標函數f(x)的一階導數趨於0。x^(k+1)由g（x^(k)）的切線產生：

x^(k+1)為g（x^(k)）的切線與x軸交點。隨着x^(k) -> x^*,  g(x)->0。（g(x)為f(x)的近似）

 

注意1：牛頓法不需要計算函數值，但需要計算一階導數和二階導數值，且要求 f''(x)>0，如果f''(x)<0，則可能存在從 x^(k) 到 x^(k+1)反向調整，收斂到極大點，即越來越偏離極小點。

f''(x^(k))>0, x^(k)->x^*

f''(x^(k))<0, x^(k)遠離x^*

注意2：如果初始x^(k)的調整幅度太大可能會導致x的調整序列在極小點左右波動。

初始點g'(x⁽⁰⁾)/g''(x⁽⁰⁾)太大，造成算法失效

#ifndef _NEWTON_H_
#define _NEWTON_H_

#define NEWTON_ZERO  0.0001

typedef float (* PtrOneVarFunc)(float x);



void NewtonMethod(float x0, PtrOneVarFunc oneorderderivative, PtrOneVarFunc twoorderderivative, float epsilon);


#endif

#include<iostream>
#include<cmath>
#include "Newton.h"

using namespace std;

void NewtonMethod(float x0, PtrOneVarFunc oneorderderivative, PtrOneVarFunc twoorderderivative, float epsilon)
{
    float xk,newtonstep;
    int k;
    k = 0;
    xk = x0;

    while( abs(oneorderderivative(xk)) - epsilon > 0)
    {

#ifdef _DEBUG
        cout<<k<<"th iteration xk="<<xk<<",f'("<<xk<<")="<<oneorderderivative(xk)<<endl;
#endif
        k++;

        newtonstep = twoorderderivative(xk);
        if(abs(newtonstep) - NEWTON_ZERO < 0)
        {
            cout<<"f''("<<xk<<")=0, algorithms terminate!"<<endl;
            return ;
        }

        newtonstep = oneorderderivative(xk) / newtonstep;
        xk = xk - newtonstep;
    }

    cout<<k<<"th iteration xk="<<xk<<",f'("<<xk<<")="<<oneorderderivative(xk)<<endl;
}

#include<iostream>
#include "Newton.h"


float OneOrderDerivative(float x)
{
    return (14*x-5);//f(x)=7*x*x-5*x+2;
}


float TwoOrderDerivative(float x)
{
    return 14;//f''(x)=14;
}

int main()
{
    NewtonMethod(-109.23, OneOrderDerivative, TwoOrderDerivative, 0.001);
    return 0;
}

0th iteration xk=-109.23,f'(-109.23)=-1534.22
1th iteration xk=0.35714,f'(0.35714)=-4.57764e-05

（3）割線法

前提：f 在 [a₀,b₀] 為單峰函數， 且[a₀,b₀] 在極小點附近，不能離的太遠否則可能無法收斂。

 牛頓法需要f(x)的二階導數，如果二階導數不存在，可以采用不同點的一階導數對其近似，如 f''(x^(k))=(f'(x^(k))-f'(x^(k-1)))/(x^(k)-x^(k-1)),將其帶入牛頓迭代公式，可以得到新的迭代公式：

 x^(k+1) = x^(k)  -  f'(x^(k))  *  (x^(k)-x^(k-1)) /(f'(x^(k))-f'(x^(k-1)))   <=>   x^(k+1) = ( f'(x^(k))x^(k-1)  - f'(x^(k-1)) x^(k)  )    /   (    f'(x^(k))  -  f'(x^(k-1))  )

這個方法需要兩個初始點 x^(-1)和x⁽⁰⁾,可以看出割線法不需要計算函數值f(x^(k))，二階導數值f''(x^(k))，它使用x^(k)和x^(k-1)之間的割線產生x^(k+1)：

x^(k+1)由過x^(k)和x^(k-1)的割線產生

#ifndef _SECANT_H_
#define _SECANT_H_

#define SECANT_ZERO  0.1

typedef float (* PtrOneVarFunc)(float x);



void SecantMethod(float x0, float x_1, PtrOneVarFunc oneorderderivative, float epsilon);


#endif

#include<iostream>
#include<cmath>
#include "Secant.h"

using namespace std;

void SecantMethod(float x_1, float x0, PtrOneVarFunc oneorderderivative, float epsilon)
{
    float xk,xk_1,t1,t_1,tk;
    int k;
    k = 0;
    xk = x0;
    xk_1 = x_1;

    if(xk - xk_1 < SECANT_ZERO)
    {
        cout<<"x0 == x_1,algorithm terminate."<<endl;
        return;
    }

    while( abs(oneorderderivative(xk)) - epsilon > 0)
    {
#ifdef _DEBUG
        cout<<k<<"th iteration xk_1="<<xk_1<<",f'("<<xk_1<<")="<<oneorderderivative(xk_1)<<",xk="<<xk<<",f'("<<xk<<")="<<oneorderderivative(xk)<<endl;
#endif

        k++;

        t1 = oneorderderivative(xk);
        t_1 = oneorderderivative(xk_1);

        tk = xk;
        xk = (t1*xk_1 - t_1*xk)/(t1-t_1);
        xk_1 = tk;
    }

        cout<<k<<"th iteration xk_1="<<xk_1<<",f'("<<xk_1<<")="<<oneorderderivative(xk_1)<<",xk="<<xk<<",f'("<<xk<<")="<<oneorderderivative(xk)<<endl;

}

#include<iostream>
#include "Secant.h"


float OneOrderDerivative(float x)
{
    return (14*x-5);//f(x)=7*x*x-5*x+2;
}


int main()
{
    SecantMethod(-109.23, -20, OneOrderDerivative, 0.001);
    return 0;
}

0th iteration xk_1=-109.23,f'(-109.23)=-1534.22,xk=-20,f'(-20)=-285
1th iteration xk_1=-20,f'(-20)=-285,xk=0.357142,f'(0.357142)=-1.03116e-05

（4）擬拋物線插值法

迭代算法跟割線法類似，不同點為使用三個點的函數值近似公式中的兩個一階導數。