常见优化算法实现

这里实现的主要算法有：

一维搜索方法：

黄金分割法

二次差值法

多维搜索算法

最速下降法

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加速方法相比最速下降方法并没有什么优势，主要原因是目标函数太简单，迭代次数太少

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