matlab 二阶导（海森矩阵）的数值计算（附代码和示例）

海森矩阵中就是单值函数对自变量（可以是向量，如$\mathbf{x}=[x_1,x_2,x_3,...]$）的二阶导数：

其中元素，如G的第一行第二列元素的定义如下：

可以看出是两个一阶导数的差再除以一个微小增量。如果$\mathbf{x}$是个二元自变量，那么：

Talk is cheap. Show me the code:

function [H]=hessian_numerical(f,x0,dx,dh)
    %计算数量场f在x0处的海森矩阵H（数值计算，差分距离dx）仅适用于实数
    n=length(x0);   
    H=zeros(n,n);
    for i=1:n
        for j=1:n
            x1=x0;
            x1(i)=x1(i)+dx;
            g2=(f(x1)-f(x0))/dx;%偏导定义

            x2=x1;
            x2(j)=x2(j)+dh;
            x3=x0;
            x3(j)=x3(j)+dh;
            g1=(f(x2)-f(x3))/dx;%偏导定义

            H(i,j)=(g1-g2)/dh;%偏导定义
        end
    end
end

如果自变量是复数，而单值函数是实数，那么可以把实部和虚部分开，当做二元函数考虑，分别求二阶偏导，相应代码修改如下：

function [H]=hessian_numerical_CtoR(f,x0,dx,dh)
    %计算数量场f在x0处的海森矩阵H（数值计算，差分距离dx）实值函数对复数的二阶导
    n=length(x0);   
    %% 实数
    Hr=zeros(n,n);
    for i=1:n
        for j=1:n
            x1=x0;
            x1(i)=x1(i)+dx;
            g2=(f(x1)-f(x0))/dx;%偏导定义

            x2=x1;
            x2(j)=x2(j)+dh;
            x3=x0;
            x3(j)=x3(j)+dh;
            g1=(f(x2)-f(x3))/dx;%偏导定义

            Hr(i,j)=(g1-g2)/dh;%偏导定义
        end
    end
     %% 复数
    Hi=zeros(n,n);
    for i=1:n
        for j=1:n
            x1=x0;
            x1(i)=x1(i)+1i*dx;
            g2=(f(x1)-f(x0))/dx;%偏导定义

            x2=x1;
            x2(j)=x2(j)+1i*dh;
            x3=x0;
            x3(j)=x3(j)+1i*dh;
            g1=(f(x2)-f(x3))/dx;%偏导定义

            Hi(i,j)=(g1-g2)/dh;%偏导定义
        end
    end
    H=Hr+1i.*Hi;%合成海森矩阵
end

注意此处并没有考虑目标函数是个复变函数，hessian_numerical_CtoR可能存在问题，请谨慎使用。

试一下hessian_numerical函数和已知解析解的对比：

fx =@(x) norm(x,2)^2; %目标函数（标量）
gx=@(x) 2.*x;%梯度解析解
Gx=@(x0) 2.*eye(3);%海森矩阵（二阶导数）解析解
rng(0)
x0=rand(3,1);
Gx_num=hessian_numerical(fx,x0,1e-6,1e-6)%数值计算海森矩阵
Gx_ana=Gx(x0)%海森矩阵解析式

结果表示虽然有些误差，但还是可以接受：

如果把微分量从1e-6调高到1e-3可以解决数值不稳定的问题：

Gx_num=hessian_numerical(fx,x0,1e-3,1e-3)%数值计算海森矩阵