拟牛顿法matlab

以下是一个基本的拟牛顿法的 MATLAB 代码示例：

function [x, fval, exitflag, output] = myquasinewton(fun,x0,options,varargin)
% MYQUASINEWTON finds the minimum of a function of several variables.
%   X = MYQUASINEWTON(FUN,X0) starts at X0 and finds a minimum of the function
%   described in FUN. FUN is a function handle. The syntax for FUN is
%   Y = FUN(X) where X is a vector and Y is a scalar. FUN can also accept a 
%   second input, the value returned by OPTIMSET. 
%
%   X = MYQUASINEWTON(FUN,X0,OPTIONS) minimizes with the default optimization
%   parameters replaced by values in the structure OPTIONS, created with the
%   OPTIMSET function. See OPTIMSET for details. 
%
%   [X, FVAL] = MYQUASINEWTON(FUN,X0,OPTIONS) returns the value of the objective
%   function at X.
%
%   [X, FVAL, EXITFLAG] = MYQUASINEWTON(FUN,X0,OPTIONS) returns an EXITFLAG that
%   describes the exit condition. Possible values of EXITFLAG and the 
%   corresponding exit conditions are
%
%     1  First-order optimality measure falls below tolerance specified 
%        in OPTIONS.TolFun.
%     2  Change in x is smaller than the tolerance specified in 
%        OPTIONS.TolX.
%     0  Maximum number of function evaluations or iterations reached.
%    -1  Algorithm was terminated by the output function.
%
%   [X, FVAL, EXITFLAG, OUTPUT] = MYQUASINEWTON(FUN,X0,OPTIONS) returns a
%   structure OUTPUT with the following fields:
%
%     iterations  Number of iterations.
%     funccount   Number of function evaluations.
%     message     Exit message.
%
%   EXAMPLES:
%
%   FUN can be specified using @:
%      X = myquasinewton(@myfun,x0)
%   where myfun is a MATLAB function such as:
%      function y = myfun(x)
%      y = x(1)^2 + x(2)^2;
%
%   FUN can also be an anonymous function:
%      X = myquasinewton(@(x) sin(x(1))+cos(x(2)),x0)
%
%   FUN can include additional arguments through the VARARGIN parameter:
%      X = myquasinewton(@(x) myfun2(x,c),x0,options,c,d,e)
%   where c,d,e, etc. are additional arguments.
%
%   OPTIONS can be created with the OPTIMSET function:
%      options = optimset('Param1',value1,'Param2',value2,...)
%   See OPTIMSET for details. Commonly used options are:
%      Display     Level of display [ off | iter | notify | final ]
%      MaxIter     Maximum number of iterations allowed.
%      TolFun      Termination tolerance on the function value.
%      TolX        Termination tolerance on X.
%
%   See also FMINSEARCH, FZERO, FUNCTION_HANDLE, OPTIMSET.

% Parse inputs
if nargin < 2
    error('MYQUASINEWTON requires at least two input arguments.');
end
if nargin < 3 || isempty(options)
    options = optimset('GradObj','on');
end
if nargin < 4
    varargin = {};
end

% Set algorithm parameters
DisplayIter = strcmp(options.Display,'iter');
TolFun = options.TolFun;
TolX = options.TolX;
MaxIter = options.MaxIter;
GradObj = strcmp(options.GradObj,'on');
DerivStep = options.DerivativeCheck;

% Initialize variables
x = x0(:); % Make sure x is a column vector
n = length(x);
fval = feval(fun,x,varargin{:});
funccount = 1;
iterations = 0;
delta = 1e-8; % Initial step size
H = eye(n); % Initial Hessian approximation

% Main loop
while iterations < MaxIter

    % Compute gradient and Hessian approximation
    if GradObj
        [fval,grad] = feval(fun,x,varargin{:});
        funccount = funccount + 1;
    else
        grad = zeros(n,1);
        for i = 1:n
            xplus = x;
            xplus(i) = xplus(i) + DerivStep;
            xminus = x;
            xminus(i) = xminus(i) - DerivStep;
            grad(i) = (feval(fun,xplus,varargin{:}) - ...
                       feval(fun,xminus,varargin{:})) / (2*DerivStep);
            funccount = funccount + 2;
        end
    end
    p = -H * grad;
    pk = p' * grad;

    % Check for convergence
    if norm(p) < TolX
        exitflag = 2;
        message = sprintf('Change in x too small: %e',norm(p));
        break;
    elseif norm(grad) < TolFun
        exitflag = 1;
        message = sprintf('First-order optimality measure falls below tolerance: %e',norm(grad));
        break;
    elseif funccount >= options.MaxFunEvals
        exitflag = 0;
        message = 'Maximum number of function evaluations reached.';
        break;
    elseif iterations >= MaxIter
        exitflag = 0;
        message = 'Maximum number of iterations reached.';
        break;
    end

    % Perform line search
    t = 1;
    while feval(fun,x+t*p,varargin{:}) > fval + t*delta*pk
        t = t/2;
    end

    % Update x and Hessian approximation
    xnew = x + t*p;
    s = xnew - x;
    y = (feval(fun,xnew,varargin{:}) - fval) / t - grad' * s;
    if s' * y >= eps * norm(s) * norm(y) % Positive curvature test
        H = H + (s-H*y)*(s-H*y)'/(s'*y) + ...
            (y-H*s)*(y-H*s)'/(y'*s);
    end
    x = xnew;
    iterations = iterations + 1;

    % Display status
    if DisplayIter
        fprintf('Iteration %d: fval=%g\n',iterations,fval);
    end

end

% Output results
output.iterations = iterations;
output.funccount = funccount;
output.message = message;

该算法使用了 BFGS 方法来近似 Hessian 矩阵的逆。要使用该代码，您需要提供一个函数句柄来计算目标函数及其梯度。您还可以通过选项结构来指定优化参数，例如最大迭代次数、函数评估次数等。在运行算法后，该代码将返回一个近似最优解、目标函数的值、退出标志和一个输出结构，其中包含有关迭代次数、函数评估次数和退出消息的信息。