海森矩陣中就是單值函數(shù)對自變量（可以是向量，如$\mathbf{x}=[x_1,x_2,x_3,...]$）的二階導數(shù)：

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

可以看出是兩個一階導數(shù)的差再除以一個微小增量。如果$\mathbf{x}$是個二元自變量，那么：

Talk is cheap. Show me the code:

function [H]=hessian_numerical(f,x0,dx,dh)
    %計算數(shù)量場f在x0處的海森矩陣H（數(shù)值計算，差分距離dx）僅適用于實數(shù)
    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

如果自變量是復數(shù)，而單值函數(shù)是實數(shù)，那么可以把實部和虛部分開，當做二元函數(shù)考慮，分別求二階偏導，相應代碼修改如下：

function [H]=hessian_numerical_CtoR(f,x0,dx,dh)
    %計算數(shù)量場f在x0處的海森矩陣H（數(shù)值計算，差分距離dx）實值函數(shù)對復數(shù)的二階導
    n=length(x0);   
    %% 實數(shù)
    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
     %% 復數(shù)
    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

注意此處并沒有考慮目標函數(shù)是個復變函數(shù)，hessian_numerical_CtoR可能存在問題，請謹慎使用。

試一下hessian_numerical函數(shù)和已知解析解的對比：

fx =@(x) norm(x,2)^2; %目標函數(shù)（標量）
gx=@(x) 2.*x;%梯度解析解
Gx=@(x0) 2.*eye(3);%海森矩陣（二階導數(shù)）解析解
rng(0)
x0=rand(3,1);
Gx_num=hessian_numerical(fx,x0,1e-6,1e-6)%數(shù)值計算海森矩陣
Gx_ana=Gx(x0)%海森矩陣解析式

結(jié)果表示雖然有些誤差，但還是可以接受：

如果把微分量從1e-6調(diào)高到1e-3可以解決數(shù)值不穩(wěn)定的問題：

Gx_num=hessian_numerical(fx,x0,1e-3,1e-3)%數(shù)值計算海森矩陣

文章來源地址http://www.zghlxwxcb.cn/news/detail-775940.html

到了這里，關于matlab 二階導（海森矩陣）的數(shù)值計算（附代碼和示例）的文章就介紹完了。如果您還想了解更多內(nèi)容，請在右上角搜索TOY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關文章，希望大家以后多多支持TOY模板網(wǎng)！