Main Content

使用梯度和黑塞矩阵的最小化

此示例说明如何使用显式三对角黑塞矩阵 H(x) 求解非线性最小化问题。问题可以表示为找出使以下方程最小化的 x

f(x)=i=1n-1((xi2)(xi+12+1)+(xi+12)(xi2+1)),

其中 n = 1000。

此示例末尾的辅助函数 brownfgh 计算 f(x),其梯度 g(x),和黑塞矩阵 H(x)。要指定 fminunc 求解器使用导数信息,请使用 optimoptions 设置 SpecifyObjectiveGradientHessianFcn 选项。要将黑塞矩阵与 fminunc 结合使用,必须使用 'trust-region' 算法。

options = optimoptions(@fminunc,'Algorithm','trust-region',...
    'SpecifyObjectiveGradient',true,'HessianFcn','objective');

将参数 n 设置为 1000,并将初始点 xstart 设置为 -1(对于奇数分量)和 +1(对于偶数分量)。

n = 1000;
xstart = -ones(n,1);
xstart(2:2:n) = 1;

f 的最小值。

[x,fval,exitflag,output] = fminunc(@brownfgh,xstart,options);
Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

检查解和求解过程。

disp(fval)
   2.8709e-17
disp(exitflag)
     1
disp(output)
         iterations: 7
          funcCount: 8
           stepsize: 0.0039
       cgiterations: 7
      firstorderopt: 4.7948e-10
          algorithm: 'trust-region'
            message: 'Local minimum found....'
    constrviolation: []

函数 f(x) 是平方值乘幂的和,因此为非负值。解 fval 几乎为零,显然是最小值。退出标志 1 也表示 fminunc 已求得解。output 结构体表明 fminunc 只需 7 次迭代就找到了解。

显示该解的最大和最小元素。

disp(max(x))
   1.1987e-10
disp(min(x))
  -1.1987e-10

该解非常接近于当 x 的所有元素均满足 x = 0 时的点。

辅助函数

以下代码创建 brownfgh 辅助函数。

function [f,g,H] = brownfgh(x)
%BROWNFGH  Nonlinear minimization problem (function, its gradients
% and Hessian)
% Documentation example        

%   Copyright 1990-2008 The MathWorks, Inc.

% Evaluate the function.
  n=length(x); y=zeros(n,1);
  i=1:(n-1);
  y(i)=(x(i).^2).^(x(i+1).^2+1)+(x(i+1).^2).^(x(i).^2+1);
  f=sum(y);
%
% Evaluate the gradient.
  if nargout > 1
     i=1:(n-1); g = zeros(n,1);
     g(i)= 2*(x(i+1).^2+1).*x(i).*((x(i).^2).^(x(i+1).^2))+...
             2*x(i).*((x(i+1).^2).^(x(i).^2+1)).*log(x(i+1).^2);
     g(i+1)=g(i+1)+...
              2*x(i+1).*((x(i).^2).^(x(i+1).^2+1)).*log(x(i).^2)+...
              2*(x(i).^2+1).*x(i+1).*((x(i+1).^2).^(x(i).^2));
  end
%
% Evaluate the (sparse, symmetric) Hessian matrix
  if nargout > 2
     v=zeros(n,1);
     i=1:(n-1);
     v(i)=2*(x(i+1).^2+1).*((x(i).^2).^(x(i+1).^2))+...
            4*(x(i+1).^2+1).*(x(i+1).^2).*(x(i).^2).*((x(i).^2).^((x(i+1).^2)-1))+...
            2*((x(i+1).^2).^(x(i).^2+1)).*(log(x(i+1).^2));
     v(i)=v(i)+4*(x(i).^2).*((x(i+1).^2).^(x(i).^2+1)).*((log(x(i+1).^2)).^2);
     v(i+1)=v(i+1)+...
              2*(x(i).^2).^(x(i+1).^2+1).*(log(x(i).^2))+...
              4*(x(i+1).^2).*((x(i).^2).^(x(i+1).^2+1)).*((log(x(i).^2)).^2)+...
              2*(x(i).^2+1).*((x(i+1).^2).^(x(i).^2));
     v(i+1)=v(i+1)+4*(x(i).^2+1).*(x(i+1).^2).*(x(i).^2).*((x(i+1).^2).^(x(i).^2-1));
     v0=v;
     v=zeros(n-1,1);
     v(i)=4*x(i+1).*x(i).*((x(i).^2).^(x(i+1).^2))+...
            4*x(i+1).*(x(i+1).^2+1).*x(i).*((x(i).^2).^(x(i+1).^2)).*log(x(i).^2);
     v(i)=v(i)+ 4*x(i+1).*x(i).*((x(i+1).^2).^(x(i).^2)).*log(x(i+1).^2);
     v(i)=v(i)+4*x(i).*((x(i+1).^2).^(x(i).^2)).*x(i+1);
     v1=v;
     i=[(1:n)';(1:(n-1))'];
     j=[(1:n)';(2:n)'];
     s=[v0;2*v1];
     H=sparse(i,j,s,n,n);
     H=(H+H')/2;
  end
end

相关主题