常见优化算法实现 这里实现的主要算法有： 一维搜索方法：
黄金分割法
二次差值法
多维搜索算法最速下降法
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加速方法相比最速下降方法并没有什么优势，主要原因是目标函数太简单，迭代次数太少
拟牛顿法相比最速下降法也没有什么优势，我想也是基于同样的原因