最速下降法和牛顿方法的Python实现和MATLAB实现

算法来源：《数值最优化方法~高立》

算法目的：实现函数的局部最优化寻找，以二元函数为例，展示了最速下降法和牛顿寻优的算法过程

主要Python模块：numpy,sympy

（1）Python实现

（2）MATLAB实现

（3）比较

（1）Python实现

A 最速下降法

[python] view plain copy

1. # -*- coding: utf-8 -*-  
2. """ 
3. Created on Sat Oct 01 15:01:54 2016 
4. @author: zhangweiguo 
5. """  
6. import sympy,numpy  
7. import math  
8. import matplotlib.pyplot as pl  
9. from mpl_toolkits.mplot3d import Axes3D as ax3  
10. #最速下降法，二维实验  
11. def SD(x0,G,b,c,N,E):  
12.     f = lambda x: 0.5 * (numpy.dot(numpy.dot(x.T, G), x)) + numpy.dot(b.T, x) + c  
13.     f_d = lambda x: numpy.dot(G, x)+b  
14.     X=x0;Y=[];Y_d=[];  
15.     xx=sympy.symarray('xx',(2,1))  
16.     n = 1  
17.     ee = f_d(x0)  
18.     e=(ee[0]**2+ee[1]**2)**0.5  
19.     Y.append(f(x0)[0,0]);Y_d.append(e)  
20.     a=sympy.Symbol('a',real=True)  
21.     print '第%d次迭代：e=%d' % (n, e)  
22.     while n<N and e>E:  
23.         n=n+1  
24.         yy=f(x0-a*f_d(x0))  
25.         yy_d=sympy.diff(yy[0,0],a,1)  
26.         a0=sympy.solve(yy_d)  
27.         x0=x0-a0*f_d(x0)  
28.         X=numpy.c_[X,x0]  
29.         Y.append(f(x0)[0,0])  
30.         ee = f_d(x0)  
31.         e = math.pow(math.pow(ee[0,0],2)+math.pow(ee[1,0],2),0.5)  
32.         Y_d.append(e)  
33.         print '第%d次迭代：e=%s'%(n,e)  
34.     return X,Y,Y_d  
35. if __name__=='__main__':  
36.     G=numpy.array([[21.0,4.0],[4.0,15.0]])  
37.     b=numpy.array([[2.0],[3.0]])  
38.     c=10.0  
39.     f = lambda x: 0.5 * (numpy.dot(numpy.dot(x.T, G), x)) + numpy.dot(b.T, x) + c  
40.     f_d = lambda x: numpy.dot(G, x) + b  
41.     x0=numpy.array([[-30.0],[100.0]])  
42.     N=40  
43.     E=10**(-6)  
44.     X, Y, Y_d=SD(x0,G,b,c,N,E)  
45.     X=numpy.array(X)  
46.     Y=numpy.array(Y)  
47.     Y_d=numpy.array(Y_d)  
48.     figure1=pl.figure('trend')  
49.     n=len(Y)  
50.     x=numpy.arange(1,n+1)  
51.     pl.subplot(2,1,1)  
52.     pl.semilogy(x,Y,'r*')  
53.     pl.xlabel('n')  
54.     pl.ylabel('f(x)')  
55.     pl.subplot(2,1,2)  
56.     pl.semilogy(x,Y_d,'g*')  
57.     pl.xlabel('n')  
58.     pl.ylabel("|f'(x)|")  
59.     figure2=pl.figure('behave')  
60.     x=numpy.arange(-110,110,1)  
61.     y=x  
62.     [xx,yy]=numpy.meshgrid(x,y)  
63.     zz=numpy.zeros(xx.shape)  
64.     n=xx.shape[0]  
65.     for i in xrange(n):  
66.         for j in xrange(n):  
67.             xxx=numpy.array([xx[i,j],yy[i,j]])  
68.             zz[i,j]=f(xxx.T)  
69.     ax=ax3(figure2)  
70.     ax.contour3D(xx,yy,zz)  
71.     ax.plot3D(X[0,:],X[1,:],Y,'ro--')  
72.     pl.show()  

结果：

第1次迭代：e=1401
第2次迭代：e=285.423850209
第3次迭代：e=36.4480186834
第4次迭代：e=7.42196747556
第5次迭代：e=0.947769462918
第6次迭代：e=0.192995789132
第7次迭代：e=0.0246451518433
第8次迭代：e=0.00501853110317
第9次迭代：e=0.000640855749362
第10次迭代：e=0.000130498466038
第11次迭代：e=1.66643765921e-05
第12次迭代：e=3.39339326369e-06
第13次迭代：e=4.33329103766e-07

B 牛顿方法

[python] view plain copy

1. # -*- coding: utf-8 -*-  
2. """ 
3. Created on Sat Oct 01 15:01:54 2016 
4. @author: zhangweiguo 
5. """  
6. import numpy  
7. import math  
8. import matplotlib.pyplot as pl  
9. from mpl_toolkits.mplot3d import Axes3D as ax3  
10. #Newton寻优，二维实验  
11. def NE(x0,y0,N,E):  
12.     X1=[];X2=[];Y=[];Y_d=[]  
13.     n = 1  
14.     ee = g(x0,y0)  
15.     e=(ee[0,0]**2+ee[1,0]**2)**0.5  
16.     X1.append(x0)  
17.     X2.append(y0)  
18.     Y.append(f(x0,y0))  
19.     Y_d.append(e)  
20.     print '第%d次迭代：e=%s' % (n, e)  
21.     while n<N and e>E:  
22.         n=n+1  
23.         d=-numpy.dot(numpy.linalg.inv(G(x0,y0)),g(x0,y0))  
24.         #d=-numpy.dot(numpy.linalg.pinv(G(x0,y0)),g(x0,y0))  
25.         x0=x0+d[0,0]  
26.         y0=y0+d[1,0]  
27.         ee = g(x0, y0)  
28.         e = (ee[0,0] ** 2 + ee[1,0] ** 2) ** 0.5  
29.         X1.append(x0)  
30.         X2.append(y0)  
31.         Y.append(f(x0, y0))  
32.         Y_d.append(e)  
33.         print '第%d次迭代：e=%s'%(n,e)  
34.     return X1,X2,Y,Y_d  
35. if __name__=='__main__':  
36.     f = lambda x,y: 3*x**2+3*y**2-x**2*y                    #原函数  
37.     g = lambda x,y: numpy.array([[6*x-2*x*y],[6*y-x**2]])   #一阶导函数向量  
38.     G = lambda x,y: numpy.array([[6-2*y,-2*x],[-2*x,6]])    #二阶导函数矩阵  
39.     x0=-2;y0=4  
40.     N=10;E=10**(-6)  
41.     X1,X2,Y,Y_d=NE(x0,y0,N,E)  
42.     X1=numpy.array(X1)  
43.     X2=numpy.array(X2)  
44.     Y=numpy.array(Y)  
45.     Y_d=numpy.array(Y_d)  
46.     figure1=pl.figure('trend')  
47.     n=len(Y)  
48.     x=numpy.arange(1,n+1)  
49.     pl.subplot(2,1,1)  
50.     pl.semilogy(x,Y,'r*')  
51.     pl.xlabel('n')  
52.     pl.ylabel('f(x)')  
53.     pl.subplot(2,1,2)  
54.     pl.semilogy(x,Y_d,'g*')  
55.     pl.xlabel('n')  
56.     pl.ylabel("|f'(x)|")  
57.     figure2=pl.figure('behave')  
58.     x=numpy.arange(-6,6,0.1)  
59.     y=x  
60.     [xx,yy]=numpy.meshgrid(x,y)  
61.     zz=numpy.zeros(xx.shape)  
62.     n=xx.shape[0]  
63.     for i in xrange(n):  
64.         for j in xrange(n):  
65.             zz[i,j]=f(xx[i,j],yy[i,j])  
66.     ax=ax3(figure2)  
67.     ax.contour3D(xx,yy,zz,15)  
68.     ax.plot3D(X1,X2,Y,'ro--')  
69.     pl.show()  

（2）MATLAB实现

A 最速下降法

先建立调用函数，输入不同的系数矩阵G和不同的初始迭代点会有不同的表现

再调用一下即可

[plain] view plain copy

1. function [X Y Y_d]=SD(x0)  
2. %%%最速下降法  
3. %%%采用精确搜索  
4. tic  
5. G=[21 4;4 15];  
6. %G=[21 4;4 1];  
7. b=[2 3]';  
8. c=10;  
9. %x=[x1;x2]  
10. f=@(x)0.5*(x'*G*x)+b'*x+c;%目标函数  
11. f_d=@(x)G*x+b;%一阶导数  
12. %设定终止条件  
13. N=100;E=10^(-6);  
14.   
15. n=1;  
16. e=norm(abs(f_d(x0)));  
17. X=x0;Y=f(x0);Y_d=e;  
18. syms a real  
19. while n<N && e>E  
20.     n=n+1;  
21.     f1=f(x0-a*f_d(x0));  
22.     a0=solve(diff(f1,'a',1));  
23.     a0=double(a0);  
24.     x0=x0-a0*f_d(x0);  
25.     X(:,n)=x0;Y(n)=f(x0);  
26.     e=norm(abs(f_d(x0)));  
27.     Y_d(n)=e;  
28. end  
29. toc  
30. figure(1)  
31. subplot(2,1,1)  
32. semilogy(Y,'r*')  
33. title(['函数最优值：' num2str(Y(n))])  
34. ylabel('f(x)')  
35. xlabel('迭代次数')  
36. subplot(2,1,2)  
37. semilogy(Y_d,'g<')  
38. ylabel('f''(x)')  
39. xlabel('迭代次数')  
40. title(['导函数最优值：' num2str(Y_d(n))])  
41. figure(2)  
42. x1=-110:1:110;y1=x1;  
43. [X1 Y1]=meshgrid(x1,y1);  
44. nn=length(x1);  
45. Z1=zeros(nn,nn);  
46. for i=1:nn  
47.     for j=1:nn  
48.         Z1(i,j)=f([X1(i,j);Y1(i,j)]);  
49.     end  
50. end  
51. hold on  
52. contour(X1,Y1,Z1)  
53. plot(X(1,:),X(2,:),'o-','linewidth',1)  
54. end  

调用：

x0为初始点，可替换

[html] view plain copy

1. x0=[-30;100];  
2. [X Y Y_d]=SD(x0);  

结果：

B 牛顿方法

[plain] view plain copy

1. function [X Y Y_d]=NT(x0)  
2. %%%固定步长为1，Newton方法  
3. %%%设定终止条件  
4. tic  
5. N=100;E=10^(-6);  
6. %%%f为原函数，g为导函数，G为二阶导数  
7. f=@(x)3*x(1).^2+3*x(2).^2-x(1).^2*x(2);  
8. g=@(x)[6*x(1)-2*x(1)*x(2);6*x(2)-x(1)^2];  
9. G=@(x)[6-2*x(2),-2*x(1);-2*x(1),6];  
10. e=norm(g(x0));  
11. X=x0;Y=f(x0);Y_d=e;  
12. n=1;  
13. while n<100 && e>E  
14.     n=n+1;  
15.     d=-inv(G(x0))*g(x0);  
16.     %d=-pinv(G(x0))*g(x0);  
17.     x0=x0+d;  
18.     X(:,n)=x0;Y(n)=f(x0);  
19.     e=norm(g(x0));  
20.     Y_d(n)=e;  
21. end  
22. toc  
23. figure(1)  
24. subplot(2,1,1)  
25. semilogy(Y,'r*')  
26. title(['函数最优值：' num2str(Y(n))])  
27. ylabel('f(x)')  
28. xlabel('迭代次数')  
29. subplot(2,1,2)  
30. semilogy(Y_d,'g<')  
31. ylabel('f''(x)')  
32. xlabel('迭代次数')  
33. title(['导函数最优值：' num2str(Y_d(n))])  
34. figure(2)  
35. x1=-6:0.05:6;y1=x1;  
36. [X1 Y1]=meshgrid(x1,y1);  
37. nn=length(x1);  
38. Z1=zeros(nn,nn);  
39. for i=1:nn  
40.     for j=1:nn  
41.         Z1(i,j)=f([X1(i,j);Y1(i,j)]);  
42.     end  
43. end  
44. hold on  
45. contour(X1,Y1,Z1,20)  
46. plot(X(1,:),X(2,:),'go-','linewidth',1)  
47. end  

[html] view plain copy

1. %x0=[1.5;1.5];  
2. %x0=[-2;4];  
3. x0=[0;3];  
4. %x0=[0;3]时矩阵为奇异的  
5. %矩阵G为病态或奇异的时候，pinv能适当改善结果，更差的时候需采用正则化的方法  
6. [X Y Y_d]=NT(x0);  
7. disp('[1.5;1.5]的最优值：')  
8. disp(min(Y))  

结果：

（3）比较

相同的算法，在不同的环境下，由于截断参数不同，会有些许的差异，MATLAB更易上手，Python却更严谨，但Python在符号运算时却有诸多不便，两者之间的差异，只能说各有所长各有所短，能结合起来用也许效果会更好。