常见优化算法统一框架下的实现：最速下降法，partan加速的最速下降法，共轭梯度法，牛顿法，拟牛顿法，黄金分割法，二次插值法

常见优化算法实现

这里实现的主要算法有：

一维搜索方法：

黄金分割法

二次差值法

多维搜索算法

最速下降法

partan加速的最速下降法

共轭梯度法

牛顿法

拟牛顿法

使用函数表示一个用于优化的目标，包括其梯度函数和hessian矩阵函数

import numpy as np
import math

#用于测试的一个多元函数的例子
def f(x):
    return (x[0]-1)**2+5*(x[1]-5)**2+(x[2]-1)**2+5*(x[3]-5)**2

#f(x)函数的gradient向量计算函数
def g(x):
    return np.array([2*(x[0]-1),10*(x[1]-5),2*(x[2]-1),10*(x[3]-5)])

#f(x)函数的hessian矩阵的逆矩阵计算函数
def hi(x=None):
    h=[1/2,1/10,1/2,1/10]
    return np.diag(h)

拟牛顿法

def quasi_newton(f=f,x0=np.zeros(4),gradient=g,acc=0.001):
    k=0
    x=x0
    xp=None
    hpk=None
    gpk=None
    
    while True:
        gk=gradient(x)
        #print(gk)
        if np.sum(gk**2)<=acc:
            #print("迭代 %d 次"%(k+1))
            return x,np.round(f(x),5)
        if k==0:
            hik=np.eye(x0.shape[0])
        else:
            dx=x-xp
            dg=gk-gpk
            
            temp = (dx-np.dot(hpk,dg)).reshape((-1,1))
            hik=hpk + np.dot(temp,temp.transpose())/(np.dot(temp.transpose(),dg.reshape((-1,1))))
            #print(hik)
            
        pk=-1*np.dot(hik,gk)
        alpha,y=quadraticInterploation(lambda alpha:(f(alpha*pk+x)),0,10,0.001)
        #更新变量
        x=alpha*pk+x
        hpk=hik
        xp=x
        gpk=gk
        k+=1

共轭方向法

def conjugate_direction(f=f,x0=np.zeros(4),gradient=g,acc=0.001):
    k=0
    x=x0
    #设置初值
    gpk=x0
    ppk=x0
    
    while True:
        gk=gradient(x)
        #print(gk)
        if np.sum(gk**2)<=acc:
            #print("迭代 %d 次"%(k+1))
            return x,np.round(f(x),5)
        if k==0:
            pk=-1*gk
        else:
            betak=np.sum(gk*gk)/np.sum(gpk*gpk)
            pk=-1*gk+betak*ppk
        #lambda表达式可以使用上层函数中的变量，这样对于不同的上下文，就是不同的函数
        alpha,y=quadraticInterploation(lambda alpha:(f(alpha*pk+x)),0,10,0.001)
        x=alpha*pk+x
        ppk=pk
        gpk=gk
        k+=1

最速下降法

#最速下降法
def steepestDescent(f=f,x0=np.zeros(4),gradient=g,acc=0.001):
    k=0
    x=x0
    while True:
        gk=gradient(x)
        pk=-1*gk
        if np.sum(gk**2)<=acc:
            #print("迭代 %d 次"%(k+1))
            return x,f(x)
        #lambda表达式可以使用上层函数中的变量，这样对于不同的上下文，就是不同的函数
        alpha,y=quadraticInterploation(lambda alpha:(f(alpha*pk+x)),0,10,0.001)
        x=alpha*pk+x
        k+=1

牛顿法

def newton(f=f,x0=np.zeros(4),gradient=g,hessian=hi,acc=0.001):
    k=0
    x=x0
    while True:
        gk=gradient(x)
        hik=hessian(x)
        pk=-1*np.dot(gk,hik)
        if np.sum(gk**2)<=acc:
            #print("迭代 %d 次"%(k+1))
            return x,f(x)
        #lambda表达式可以使用上层函数中的变量，这样对于不同的上下文，就是不同的函数
        alpha,y=quadraticInterploation(lambda alpha:(f(alpha*pk+x)),0,10,0.001)
        x=alpha*pk+x
        k+=1

使用partan加速的最速下降法

def partan(f=f,x0=np.zeros(4),gradient=g,acc=0.001,N=3):
    k=0
    x=x0
    xp1=x0
    xp2=x0
    while True:
        if k>=N and k%3==0:
            pk=x-xp2
        else:
            gk=gradient(x)
            pk=-1*gk
        if np.sum(pk**2)<=acc:
            #print("迭代 %d 次"%(k+1))
            return x,f(x)
        #lambda表达式可以使用上层函数中的变量，这样对于不同的上下文，就是不同的函数
        alpha,y=quadraticInterploation(lambda alpha:(f(alpha*pk+x)),0,10,0.001)
        xp2=xp1
        xp1=x
        x=alpha*pk+x
        k+=1

一维搜索的黄金分割方法

def goldenSegmantation(f,a,b,acc):
    x1=a+0.382*(b-a)
    x2=b-(x1-a)
    R=f(x1);G=f(x2)
    #因为浮点数的舍入误差，可能导致a,b的大小逆转
    while abs(b-a)>acc and a<=x1<x2<=b:
       #print(abs(b-a))
        if R>G:
            a=x1
            x1=x2
            R=G
            x2=b-(x1-a)
            G=f(x2)
        else:
            b=x2
            x2=x1
            G=R
            x1=a+(b-x2)
            R=f(x1)
    return (a+b)/2.0,f(((a+b)/2.0))

一维搜索的二次差值方法

def quadraticInterploation(f,a,b,acc):
    assert(a<b)
    x1=a;x2=(a+b)/2;x3=b
    f1=f(x1);f2=f(x2);f3=f(x3)
    while True:
        c1=(f3-f1)/(x3-x1);c2=((f2-f1)/(x2-x1)-c1)/(x2-x3)
        xp=0.5*(x1+x3-c1/c2)
        fp=f(xp)
        if abs(xp-x2)<acc or not a<=x1<x2<x3<=b:
            if fp<f2:
                return xp,fp
            else:
                return x2,f2
        if x2<xp:
            if f2<fp:
                x3=xp;f3=fp
            else:
                x1=x2;f1=f2
                x2=xp;f2=fp
        else:
            if f2<fp:
                x1=xp;f1=fp
            else:
                x3=x2;f3=f2
                x2=xp;f2=fp

测试一维搜索方法

%timeit(goldenSegmantation(lambda x:(x**4-5),-1,1,0.0001))
%timeit(quadraticInterploation(lambda x:(x**4-5),-1,1,0.00001))
%timeit(goldenSegmantation(lambda x:(x**2-5*x+6),-10,10,0.00000005))
%timeit(quadraticInterploation(lambda x:(x**2-5*x+6),-10,10,0.000001))
%timeit(goldenSegmantation(math.sin,-1*math.pi,0,0.000001))
%timeit(quadraticInterploation(math.sin,-1*math.pi,0,0.0000001))                

11.3 µs ± 58.8 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
3.09 µs ± 18.5 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
12.5 µs ± 47.9 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
5.44 µs ± 27.1 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
7.97 µs ± 33.1 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
2.05 µs ± 19.5 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

结果分析

对于不同的目标函数，二次插值的速度均大于黄金分割方法

测试高维搜索方法

%timeit steepestDescent()
%timeit partan()
%timeit conjugate_direction()
%timeit newton()
%timeit quasi_newton()

236 µs ± 2.39 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
297 µs ± 2.39 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
197 µs ± 1.49 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
125 µs ± 276 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
224 µs ± 1.28 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

结果分析

从结果看出来，partan加速方法相比最速下降方法并没有什么优势，主要原因是目标函数太简单，迭代次数太少

拟牛顿法相比最速下降法也没有什么优势，我想也是基于同样的原因