一、算法原理
无约束优化问题的目标函数为:
多维极值的牛顿法,不同于之前介绍过的梯度下降法,该方法引入了二阶导数的信息,假设当前迭代到第 kk 次,将目标函数在自变量 xk 处展开为二阶泰勒级数: (1)
f(x)两端同时对 x求导,导数为零求得的极值点就是下一次迭代的取值 xk+1。
(2)
式中为梯度,为海塞矩阵。
将记做g ,记做 H,g 与 H 的形式分别如下:
则式(2)可化简为 (3)
综上,
式中-Hk^(-1)gk即为搜索方向。
二、matlab代码
f=@(x1,x2) (x1-4)^2+(x2+2)^2+1;
[X,result]=Min_Newton(f,[-0.6 0.2],1e-6,100)
function [X,result]=Min_Newton(f,x0,eps,n)
TiDu=gradient(sym(f),symvar(sym(f)));% 计算出梯度表达式
Haisai=jacobian(TiDu,symvar(TiDu));%计算出海塞矩阵表达式
Var_Tidu=symvar(TiDu); %梯度表达式中变量的个数
Var_Haisai=symvar(Haisai); %海塞矩阵中变量的个数
Var_Num_Tidu=length(Var_Tidu); %梯度的维数
Var_Num_Haisai=length(Var_Haisai); %海塞矩阵的维数
TiDu=matlabFunction(TiDu);%将梯度表达式转换为匿名函数
flag = 0;
if Var_Num_Haisai == 0 %海塞矩阵变量的个数为零,也就是说海塞矩阵是常数
Haisai=double((Haisai));
flag=1; %海塞矩阵为常量的标志
end
%求当前点梯度与海赛矩阵的逆
f_cal='f(';
TiDu_cal='TiDu(';
Haisai_cal='Haisai(';
for k=1:length(x0) %求得初始变量的x0的元素个数
f_cal=[f_cal,'x0(',num2str(k),'),'];%组装f_cal=f(x0(k))求得该点函数值
for j=1: Var_Num_Tidu %求得梯度中的元素个数
if char(Var_Tidu(j)) == ['x',num2str(k)]
TiDu_cal=[TiDu_cal,'x0(',num2str(k),'),'];%组装TiDu_cal=TiDu_cal(x0(k)求得该点梯度值
end
end
for j=1:Var_Num_Haisai
if char(Var_Haisai(j)) == ['x',num2str(k)]
Haisai_cal=[Haisai_cal,'x0(',num2str(k),'),'];%组装Haisai_cal=Haisai_cal(x0(k)求得该点海塞矩阵的值
end
end
end
Haisai_cal(end)=')'; %完成海塞矩阵封装
TiDu_cal(end)=')';%完成梯度封装
f_cal(end)=')';%完成函数封装
switch flag %根据标志位判断海塞矩阵表达式中是否有参数
case 0 %有参数
Haisai=matlabFunction(Haisai);
dk='-eval(Haisai_cal)^(-1)*eval(TiDu_cal)';
case 1 %无参数
dk='-Haisai^(-1)*eval(TiDu_cal)';
Haisai_cal='Haisai';
end
i=1;
while i < n %设置最大迭代次数n
if rcond(eval(Haisai_cal)) < 1e-6 %计算海塞矩阵的条件数 条件数越大,逆阵结果越不稳定
disp('海赛矩阵病态!'); %条件数超出范围,示为病态矩阵
break;
end
x0=x0(:)+eval(dk); %eval函数将字符串转换为指令
if norm(eval(TiDu_cal)) < eps
X=x0;
result=eval(f_cal);
return;
end
i=i+1;
end
disp('无法收敛!'); %超出迭代范围
X=[];
result=[];
end