hessian.m

来自「基于matlab的约束非线性规划算法库」· M 代码 · 共 42 行

function H = hessian(func,X,v,params)
%Call: H = hessian(func,X,v,params)
%Approximate the Hessian at X of the function defined in func
%H=sigma(v(k)*Gk) where k=1,2, ......, length(v) and Gk is an
%approximation of the Hessian corresponding to the k_th function
%
m=length(v);
n=length(X);
eps2=eps^(1/3);
eps1=eps2;
for k=1:n
    xk=X(k);
    epsk=max(abs(xk),1)*eps2;
    for j=1:k
        xj=X(j);
        epsj=max(abs(xj),1)*eps1;
        X(k)=xk+epsk;
        X(j)=X(j)+epsj;
        [f1 dummy]=feval(func,X,0,params{:});
        X(k)=xk;
        X(j)=xj;
        X(k)=xk+epsk;
        X(j)=X(j)-epsj;
        [f2 dummy]=feval(func,X,0,params{:});
        X(k)=xk;
        X(j)=xj;
        X(k)=xk-epsk;
        X(j)=X(j)-epsj;
        [f4 dummy]=feval(func,X,0,params{:});
        X(k)=xk;
        X(j)=xj;
        X(k)=xk-epsk;
        X(j)=X(j)+epsj;
        [f3 dummy]=feval(func,X,0,params{:});
        X(j)=xj;
        ww=(f1-f2-f3+f4)/(4*epsk*epsj);
        S=sum(ww.*v);
        H(j,k)=S;
        H(k,j)=S;
    end % of for j
    X(k)=xk;
end % of for k