Hessian矩阵H hessian矩阵的正定负定

3.2 无约束问题的MATLAB解法

3.2.1 知识准备

1、Hessian阵、正定阵与负定阵

黑塞矩阵（Hessian矩阵）：是一个多元函数的二阶偏导数构成的方阵，描述了函数的局部曲率。黑塞矩阵常用于牛顿法解决优化问题，利用黑塞矩阵可判定多元函数的极值问题。在工程实际问题的优化设计中，为了使问题简化，常常将目标函数在某点邻域展开成泰勒多项式来逼近原函数，此时函数在某点泰勒展开式的矩阵形式中会涉及到黑塞矩阵。

正定矩阵：

（1）广义定义：设是阶方阵，如果对任何非零向量，都有> 0，其中 表示的转置，就称为正定矩阵。

（2）狭义定义：一个阶的实对称矩阵是正定的的条件是当且仅当对于所有的非零实系数向量，都有。其中表示的转置。

负定矩阵：

狭义定义：设是实对称矩阵。如果对任意的实非零列矩阵有，就称A为负定矩阵。

2、黑塞矩阵用于判断极值点的性质

驻点：多元函数偏导数为零的点。
如果在驻点处黑塞矩阵为正定阵，则在该点取极小值；如果在驻点处黑塞矩阵为负定阵，则在该点处取极大值；如果在驻点处黑塞矩阵为不定阵，则该驻点不是极值点。

3、MATLAB中的subs函数

MATLAB中是符号计算函数,表示将符号表达式中的某些符号变量替换为指定的新的变量,常用调用方式为: 表示将符号表达式中的符号变量替换为新的值。

4、MATLAB中的符号计算

MATLAB最强大的莫过于它的数值计算功能，但是它的符号计算功能也不能被忽视。值得玩味的是，符号计算对于初学MATLAB编程的人来说，反而更易于从数学的角度进行理解。或者形象地说，符号计算就是让MATLAB像人一样做数学题，其结果形式也是我们平时用笔和演草纸做数学题时写的结果。

3.2.2 无约束极值问题的符号法求解

clc,clear
syms x y//定义符号变量
f=x^3-y^3+3*x^2+3*y^2-9*x;
df=jacobian(f);//求一阶偏导数
d2f=jacobian(df);//求黑塞矩阵
[xx,yy]=solve(df);//求驻点
xx=double(xx);yy=double(yy);//转化成双精度浮点型数据，下面判断特征值的正负，必须是数值型数据
for i=1:length(xx)
     a=subs(d2f,{x,y},{xx(i),yy(i)});
     b=eig(a);//求矩阵的特征值
     f=subs(f,{x,y},{xx(i),yy(i)});
     if all(b>0)
        fprintf('(%f,%f)是极小值点，对应的极小值为%f\n',xx(i),yy(i),f);
     elseif all(b<0)
        fprintf('(%f,%f)是极大值点，对应的极大值为%f\n',xx(i),yy(i),f);
     elseif any(b>0)&any(b<0)
        fprintf('(%f,%f)不是极值点\n',xx(i),yy(i));
     else
        fprintf('无法判断(%f,%f)是否是极值点\n',xx(i),yy(i));
     end
end

3.2.3 无约束问题的数值解

[x,fval]=fminunc(fun,x0,options);

其中，是一个.m文件定义的函数（或者是一个匿名函数），当只有一个返回值时，它的返回值是函数；当有两个返回值时，第二个返回值是第一个返回值的梯度向量；当有三个返回值时，它的第三个返回值是第一个返回值的黑塞矩阵。

[x,fval]=fminsearch(fun,x0,options);

具体用法与上面相同。
例题：

//在求解极小值时，可以使用函数的梯度
function [f,g]=fun3(x)
f=100*(x(2)-x(1)^2)^2+(1-x(1))^2;
g=[-400*x(1)*(x(2)-x(1)^2)-2*(1-x(1));200*(x(2)-x(1)^2)];//是原函数的梯度向量
end
options=optimset('GradObj','on');
[x,y]=fminunc('fun3',rand(1,2),options);
//也可以利用函数的二阶导数
function [f,df,d2f]=fun4(x)
f=100*(x(2)-x(1)^2)^2+(1-x(1))^2;
df=[-400*x(1)*(x(2)-x(1)^2)-2*(1-x(1));200*(x(2)-x(1)^2)];//是原函数的梯度向量
d2f=[-400*x(2)+1200*x(1)^2+2,-400*x(1);-400*x(1),200];
end
options=optimset('GradObj','on','Hessian','on);
[x,y]=fminunc('fun4',rand(1,2),options);
//一般来说，提供的信息量越大，计算越快，精度越高。
//上面的一二阶导数，都可以先用符号法计算得到

//匿名函数用法，这样就不用新建一个函数文件，但这样给默认函数能够提供的信息也相应的减少
f=@(x) 100*(x(2)-x(1)^2)^2+(1-x(1))^2;

3.2.3 函数零点和方程组的解

求解方程：

clc,clear
//第一种解法
xishu=[1,-1,2,-3];//系数从高次幂到低次幂排列
x0=roots(xishu);
//第二种解法
syms x
x0=solve(x^3-x^2+2*x-3);
x0=vpa(x0,6);//化成小数格式的数据
//第三种解法
y=@(x) x^3-x^2+2*x-3;
x=fsolve(y,rand);//只能求给定初始值附件的一个零点

求解方程组：

//符号法
syms x y
[x,y]=solve(x^2+y-6,y^2+x-6);
//数值解
f=@(x) [x(1)^2+x(2)-6;x(2)^2+x(1)-6];
xy=fsolve(f,rand(2,1));//只能求给定初始值附近的一组解

PS.optimset函数的用法可以参考MATLAB自带的帮助文档。

参考文献

[1]司守奎，孙玺菁. 数学建模算法与应用. 北京：国防工业出版社，2011.