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

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

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

主要Python模块：numpy,sympy

（1）Python实现

（2）MATLAB实现

（3）比较

（1）Python实现

A 最速下降法

# -*- coding: utf-8 -*-
"""
Created on Sat Oct 01 15:01:54 2016
@author: zhangweiguo
"""
import sympy,numpy
import math
import matplotlib.pyplot as pl
from mpl_toolkits.mplot3d import Axes3D as ax3
#最速下降法，二维实验
def SD(x0,G,b,c,N,E):
f = lambda x: 0.5 * (numpy.dot(numpy.dot(x.T, G), x)) + numpy.dot(b.T, x) + c
f_d = lambda x: numpy.dot(G, x)+b
X=x0;Y=[];Y_d=[];
xx=sympy.symarray('xx',(2,1))
n = 1
ee = f_d(x0)
e=(ee[0]**2+ee[1]**2)**0.5
Y.append(f(x0)[0,0]);Y_d.append(e)
a=sympy.Symbol('a',real=True)
print '第%d次迭代：e=%d' % (n, e)
while n<N and e>E:
n=n+1
yy=f(x0-a*f_d(x0))
yy_d=sympy.diff(yy[0,0],a,1)
a0=sympy.solve(yy_d)
x0=x0-a0*f_d(x0)
X=numpy.c_[X,x0]
Y.append(f(x0)[0,0])
ee = f_d(x0)
e = math.pow(math.pow(ee[0,0],2)+math.pow(ee[1,0],2),0.5)
Y_d.append(e)
print '第%d次迭代：e=%s'%(n,e)
return X,Y,Y_d
if __name__=='__main__':
G=numpy.array([[21.0,4.0],[4.0,15.0]])
b=numpy.array([[2.0],[3.0]])
c=10.0
f = lambda x: 0.5 * (numpy.dot(numpy.dot(x.T, G), x)) + numpy.dot(b.T, x) + c
f_d = lambda x: numpy.dot(G, x) + b
x0=numpy.array([[-30.0],[100.0]])
N=40
E=10**(-6)
X, Y, Y_d=SD(x0,G,b,c,N,E)
X=numpy.array(X)
Y=numpy.array(Y)
Y_d=numpy.array(Y_d)
figure1=pl.figure('trend')
n=len(Y)
x=numpy.arange(1,n+1)
pl.subplot(2,1,1)
pl.semilogy(x,Y,'r*')
pl.xlabel('n')
pl.ylabel('f(x)')
pl.subplot(2,1,2)
pl.semilogy(x,Y_d,'g*')
pl.xlabel('n')
pl.ylabel("|f'(x)|")
figure2=pl.figure('behave')
x=numpy.arange(-110,110,1)
y=x
[xx,yy]=numpy.meshgrid(x,y)
zz=numpy.zeros(xx.shape)
n=xx.shape[0]
for i in xrange(n):
for j in xrange(n):
xxx=numpy.array([xx[i,j],yy[i,j]])
zz[i,j]=f(xxx.T)
ax=ax3(figure2)
ax.contour3D(xx,yy,zz)
ax.plot3D(X[0,:],X[1,:],Y,'ro--')
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 牛顿方法

# -*- coding: utf-8 -*-
"""
Created on Sat Oct 01 15:01:54 2016
@author: zhangweiguo
"""
import numpy
import math
import matplotlib.pyplot as pl
from mpl_toolkits.mplot3d import Axes3D as ax3
#Newton寻优，二维实验
def NE(x0,y0,N,E):
X1=[];X2=[];Y=[];Y_d=[]
n = 1
ee = g(x0,y0)
e=(ee[0,0]**2+ee[1,0]**2)**0.5
X1.append(x0)
X2.append(y0)
Y.append(f(x0,y0))
Y_d.append(e)
print '第%d次迭代：e=%s' % (n, e)
while n<N and e>E:
n=n+1
d=-numpy.dot(numpy.linalg.inv(G(x0,y0)),g(x0,y0))
#d=-numpy.dot(numpy.linalg.pinv(G(x0,y0)),g(x0,y0))
x0=x0+d[0,0]
y0=y0+d[1,0]
ee = g(x0, y0)
e = (ee[0,0] ** 2 + ee[1,0] ** 2) ** 0.5
X1.append(x0)
X2.append(y0)
Y.append(f(x0, y0))
Y_d.append(e)
print '第%d次迭代：e=%s'%(n,e)
return X1,X2,Y,Y_d
if __name__=='__main__':
f = lambda x,y: 3*x**2+3*y**2-x**2*y #原函数
g = lambda x,y: numpy.array([[6*x-2*x*y],[6*y-x**2]]) #一阶导函数向量
G = lambda x,y: numpy.array([[6-2*y,-2*x],[-2*x,6]]) #二阶导函数矩阵
x0=-2;y0=4
N=10;E=10**(-6)
X1,X2,Y,Y_d=NE(x0,y0,N,E)
X1=numpy.array(X1)
X2=numpy.array(X2)
Y=numpy.array(Y)
Y_d=numpy.array(Y_d)
figure1=pl.figure('trend')
n=len(Y)
x=numpy.arange(1,n+1)
pl.subplot(2,1,1)
pl.semilogy(x,Y,'r*')
pl.xlabel('n')
pl.ylabel('f(x)')
pl.subplot(2,1,2)
pl.semilogy(x,Y_d,'g*')
pl.xlabel('n')
pl.ylabel("|f'(x)|")
figure2=pl.figure('behave')
x=numpy.arange(-6,6,0.1)
y=x
[xx,yy]=numpy.meshgrid(x,y)
zz=numpy.zeros(xx.shape)
n=xx.shape[0]
for i in xrange(n):
for j in xrange(n):
zz[i,j]=f(xx[i,j],yy[i,j])
ax=ax3(figure2)
ax.contour3D(xx,yy,zz,15)
ax.plot3D(X1,X2,Y,'ro--')
pl.show()




（2）MATLAB实现

A 最速下降法

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

再调用一下即可

function [X Y Y_d]=SD(x0)
%%%最速下降法
%%%采用精确搜索
tic
G=[21 4;4 15];
%G=[21 4;4 1];
b=[2 3]';
c=10;
%x=[x1;x2]
f=@(x)0.5*(x'*G*x)+b'*x+c;%目标函数
f_d=@(x)G*x+b;%一阶导数
%设定终止条件
N=100;E=10^(-6);

n=1;
e=norm(abs(f_d(x0)));
X=x0;Y=f(x0);Y_d=e;
syms a real
while n<N && e>E
n=n+1;
f1=f(x0-a*f_d(x0));
a0=solve(diff(f1,'a',1));
a0=double(a0);
x0=x0-a0*f_d(x0);
X(:,n)=x0;Y(n)=f(x0);
e=norm(abs(f_d(x0)));
Y_d(n)=e;
end
toc
figure(1)
subplot(2,1,1)
semilogy(Y,'r*')
title(['函数最优值：' num2str(Y(n))])
ylabel('f(x)')
xlabel('迭代次数')
subplot(2,1,2)
semilogy(Y_d,'g<')
ylabel('f''(x)')
xlabel('迭代次数')
title(['导函数最优值：' num2str(Y_d(n))])
figure(2)
x1=-110:1:110;y1=x1;
[X1 Y1]=meshgrid(x1,y1);
nn=length(x1);
Z1=zeros(nn,nn);
for i=1:nn
for j=1:nn
Z1(i,j)=f([X1(i,j);Y1(i,j)]);
end
end
hold on
contour(X1,Y1,Z1)
plot(X(1,:),X(2,:),'o-','linewidth',1)
end

调用：
x0为初始点，可替换

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

结果：





B 牛顿方法

function [X Y Y_d]=NT(x0)
%%%固定步长为1，Newton方法
%%%设定终止条件
tic
N=100;E=10^(-6);
%%%f为原函数，g为导函数，G为二阶导数
f=@(x)3*x(1).^2+3*x(2).^2-x(1).^2*x(2);
g=@(x)[6*x(1)-2*x(1)*x(2);6*x(2)-x(1)^2];
G=@(x)[6-2*x(2),-2*x(1);-2*x(1),6];
e=norm(g(x0));
X=x0;Y=f(x0);Y_d=e;
n=1;
while n<100 && e>E
n=n+1;
d=-inv(G(x0))*g(x0);
%d=-pinv(G(x0))*g(x0);
x0=x0+d;
X(:,n)=x0;Y(n)=f(x0);
e=norm(g(x0));
Y_d(n)=e;
end
toc
figure(1)
subplot(2,1,1)
semilogy(Y,'r*')
title(['函数最优值：' num2str(Y(n))])
ylabel('f(x)')
xlabel('迭代次数')
subplot(2,1,2)
semilogy(Y_d,'g<')
ylabel('f''(x)')
xlabel('迭代次数')
title(['导函数最优值：' num2str(Y_d(n))])
figure(2)
x1=-6:0.05:6;y1=x1;
[X1 Y1]=meshgrid(x1,y1);
nn=length(x1);
Z1=zeros(nn,nn);
for i=1:nn
for j=1:nn
Z1(i,j)=f([X1(i,j);Y1(i,j)]);
end
end
hold on
contour(X1,Y1,Z1,20)
plot(X(1,:),X(2,:),'go-','linewidth',1)
end

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

结果：





（3）比较


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