本文是关于数值优化算法里最基础的最速下降法(Steepest Descent,SD)还有Newton法的Matlab实现

本文参考了高立《数值最优化方法》北京大学出版社

马昌凤 《最优化方法及其 Matlab程序设计》电子版

似乎应该在写代码之前先把理论复习一下，但是由于懒，还是直接复制粘贴写好的代码来得爽

在讲具体算法之前本应该介绍一下线搜索算法的框架和一些准则的，但是还没有写好，先发这个吧

Armijo准则+回溯的线搜索方法

function[alpha,flag]=Armijo(fun,gfun,x,d)% abs:使用Armijo准则回溯进行非精确线搜索

% input:待搜索函数fun，梯度gfun,初始点x,方向d

% output:搜索结果alpha

%设置初始参数

rho=10e-3;

alpha=1;beta=0.5;

m=0;m_max=20;flag=0;

while (m < m_max )

flag=(feval(fun,x+alpha*d) <= feval(fun,x)+ rho*alpha*feval(gfun,x)'*d);

if flag == 1

break;

end

m=m+1;

alpha=alpha*beta;

end

SD法就是沿着下降方向最快的方向搜索的方法

function [x,val,iter,flag] =SD(fun,gfun,x0)

%abs:使用最速下降法(steepest descent)求解无约束优化问题

%input:目标函数f,梯度函数g,初始点x0

%output:最优解x,最优值y,迭代次数iter,flag表示是否正常终止

%设置初始参数

k_max=10000; %最大迭代次数

rho=0.5;sigma=0.4;flag=0;

k=0; epsilon=1e-5;xk=x0

while(k < k_max)

gk=gfun(xk); %计算梯度

dk=-gk; %计算搜索方向

if(norm(dk) < epsilon) flag=1;break; end

m=0; mk=0;

[alpha,] = Armijo(fun,gfun,xk,dk)

xk=xk+alpha*dk;

k=k+1;

end

x=xk;

val=fun(xk);

iter=k;

end

基本的Newton法

function [x,val,iter,flag] = Newton(fun,gfun,Hess,x0)

%abs:使用Newton法求解无约束优化问题

%input:目标函数fun,梯度函数gfun,黑塞矩阵Hess，初始点x0

%output:最优解x,最优值val,迭代次数iter

%设置初始参数

k_max=5000; %最大迭代次数

k=0; epsilon=1e-5;

xk=x0;

while(k < k_max)

gk=gfun(xk); %计算梯度

Gk=Hess(xk);

dk=-Gk\gk; %计算搜索方向

if(norm(dk) < epsilon) flag=1; break; end

xk=xk+dk;

k=k+1;

end

x=xk;

val=fun(xk);

iter=k;

end

阻尼牛顿法：使用Newton法的方向加上线搜索技术

function [x,val,iter,flag] = dampNewton(fun,gfun,Hess,x0)

%abs:使用阻尼Newton法求解无约束优化问题

%input:目标函数fun,梯度函数gfun,黑塞矩阵Hess，初始点x0

%output:最优解x,最优值val,迭代次数iter，flag表示是否在最大迭代次数内终止

%设置初始参数

k_max=5000; %最大迭代次数

k=0; epsilon=1e-5;

xk=x0;

while(k < k_max)

gk=gfun(xk); %计算梯度

Gk=Hess(xk);

dk=-Gk\gk; %计算搜索方向

if(norm(dk) < epsilon) flag=1; break; end

[alpha,] = Armijo(fun,gfun,xk,dk)

xk=xk+alpha*dk;

k=k+1;

end

x=xk;

val=fun(xk);

iter=k;

end

下面两种方法都是在上面的算法实现上稍加改进，就不贴代码了(其实我也没写，因为感觉差不多)

由于Hessian矩阵不一定是正定的，因此可以考虑混合牛顿法：在Hessian矩阵可逆时使用Newton法的方向，在Hessian矩阵不可逆时采取最速下降方向

另外，若Hessian矩阵不可逆可以考虑加上一个单位矩阵的倍数来修正得到正定矩阵，修正系数可以考虑以几何级数增加来搜索.这就是Levenberg-Marquardt方法.这种思想在设计求解最小二乘问题的同名算法时也用到了.

下面是一些测试的例子，里面的Rosenbrock是数值最优化里比较常用的检验函数

%测试Armijio准则

%[alpha,flag]=Armijo(@Rosenbrock,@DRosenbrock,[-1 1]',[1 -2]')

%测试SD

%[x,val,iter,flag] =SD(@Rosenbrock,@DRosenbrock,[-1.2,1]')

%测试Newton法

%[x,val,iter,flag] = Newton(@Rosenbrock,@DRosenbrock,@HRosenbrock,[-1.2,1]')

%测试阻尼Newton法

% [x,val,iter,flag] = dampNewton(@Rosenbrock,@DRosenbrock,@HRosenbrock,[-1.2,1]')

按照我上面的参数运行了一下，SD用了大概1w次，Newton法5次，阻尼Newton法大概20次.

函数里有一些参数会影响到算法的循环次数，可以自己调整