srchbre_snn.m

来自「神经网络的工具箱, 神经网络的工具箱,」· M 代码 · 共 220 行

function [xmin,fmin] = srchbre_snn(func, x1, x2, varargin);
%SRCHBRE_SNN Brent's linear search algorithm (Brent's method).
%
% Syntax
%
%   [x_min,f_min] = srchbre_snn(func, x1, x2, varargin);
%
% Description
%
%   SRCHBRE_SNN finds the minimum of the function func(x, varargin) 
%               with respect to lambda on the line x = x1 + lambda*x2.
%
%   SRCBRE_SNN(func, x1, x2, varargin) takes
%     func     - name of the function
%     x1       - double array
%     x2       - double array
%     varargin - additional arguments voor func. 
%   and returns
%     x_min    - x which minimizes func(x, varargin)
%     f_min    - func(x_min, varargin)
%

tol = 1e-7;

[ax,bx,cx,dummy,fb] = mnbrak(func,x1,x2,0.01,1,varargin{:});

ITMAX = 100;
CGOLD = 0.3819660;

a = min(ax,cx);
b = max(ax,cx);
x = bx;
w = bx;
v = bx;
fx = fb;
fw = fb;
fv = fb;

e = 0;
iter = 0;
klaar = 0;
while (~klaar & iter < ITMAX),
      iter = iter+1;
      xm = (a+b)/2;
      tol1 = tol*abs(x)+1e-10;
      tol2 = 2*tol1;
      if (abs(x-xm) <= (tol2 - (b-a)/2)),
         xmin = x;
         fmin = fx;
         klaar = 1;
      end
      d = []; %onzin, to make mcc happy
      if (~klaar),
         if (abs(e) > tol1),
            r = (x-w)*(fx-fv);
            q = (x-v)*(fx-fw);
            p = (x-v)*q-(x-w)*r;
            q = 2*(q-r);
            if (q > 0), r = -r; end
            q = abs(q);
            etemp = e;
            e = d;
            if (abs(p) >= abs(0.5*q*etemp) | p <= q*(a-x) | p >= q*(b-x)),
               if (x >= xm), e = a-x; else e = b-x; end
               d = CGOLD*e;
            else
               d = p/q;
               u = x+d;
               if (u-a < tol2 | b-u < tol2),
                  d = sign(xm-x)*tol1;
               end
            end
         else
            if (x >= xm), e = a-x; else e = b-x; end
            d = CGOLD*e;
         end
         if (abs(d) >= tol1), u = x+d; else u = x + sign(d)*tol1; end
%#function wcf_snn
         fu = feval(func,x1+u*x2,varargin{:});
         if (fu <= fx),
            if (u >= x), a = x; else b = x; end
            v = w; w = x; x = u;
            fv = fw; fw = fx; fx = fu;
         else
            if (u < x), a = u; else b = u; end
            if (fu <= fw | w == x),
               v = w; w = u;
               fv = fw; fw = fu;
            elseif (fu <= fv | v == x | v == w),
               v = u;
               fv = fu;
            end
         end
     end
end
if (~klaar),
   disp('WATCH OUT: not sufficient number of iterations');
   xmin = x;
   fmin = fx;
end

xmin = x1+xmin*x2;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ax,bx,cx,fa,fb,fc] = mnbrak(func,x1,x2,ax,bx,varargin);

% MNBRAK searches values [ax,bx,cx] such that
% func(x1 + ax*x2) > func(x1 + bx*x2) and func(x1 + cx*x2) > func(x1 + bx*x2)
% where either ax < bx < cx or ax > bx > cx.

given = 0;
if (nargin >= 2),
   if (length(x1)),
      given = 1;
   end
end
if (~given),
   x1 = 0;
end
given = 0;
if (nargin >= 3),
   if (length(x2)),
      given = 1;
   end
end
if (~given),
   x2 = ones(size(x1));
end
given = 0;
if (nargin >= 4),
   if (length(ax)),
      given = 1;
   end
end
if (~given),
   ax = 0;
end
given = 0;
if (nargin >= 5),
   if (length(bx)),
      given = 1;
   end
end
if (~given),
   bx = 1;
   if (bx == ax),
      bx = 2;
   end
end

GOLD = 1.618034;
GLIM = 100;

fa = feval(func,x1+ax*x2,varargin{:});
fb = feval(func,x1+bx*x2,varargin{:});

if (fb > fa),    % wissel ax en bx
   dummy = [ax,bx];
   ax = dummy(2);
   bx = dummy(1);
   dummy = [fa,fb];
   fa = dummy(2);
   fb = dummy(1);
end

cx = bx + GOLD*(bx-ax);
fc = feval(func,x1+cx*x2,varargin{:});

klaar = 0;
while (fb > fc & ~klaar),
      r = (bx-ax)*(fb-fc);
      q = (bx-cx)*(fb-fa);
      u = (bx - (bx-cx)*q - (bx-ax)*r)/(2*sign(q-r)*max(abs(q-r),1e-20));
                 % parabolische schatting minimum
      ulim = bx + GLIM*(cx-bx);

      if ((bx-u)*(u-cx) > 0),   % schatting tussen bx en cx
         fu = feval(func,x1+u*x2,varargin{:});
         if (fu < fc),    % inderdaad een minimum tussen bx en cx
            ax = bx;
            bx = u;
            fa = fb;
            fb = fu;
            klaar = 1;
         elseif (fu > fb),   % minimum tussen ax en u
            cx = u;
            fc = fu;
            klaar = 1;
         else               % geen minimum gevonden
            u = cx + GOLD*(cx-bx);
            fu = feval(func,x1+u*x2,varargin{:});
         end
      elseif ((cx-u)*(u-ulim) > 0),    % schatting tussen cx en ulim
         fu = feval(func,x1+u*x2,varargin{:});
         if (fu < fc),
            bx = cx;
            cx = u;
            u = cx + GOLD*(cx-bx);
            fb = fc;
            fc = fu;
            fu = feval(func,x1+u*x2,varargin{:});
         end
      elseif ((u-ulim)*(ulim-cx) >= 0)  % schatting groter dan ulim
         u = ulim;
         fu = feval(func,x1+u*x2,varargin{:});
      else             % allemaal hopeloos
         u = cx + GOLD*(cx-bx);
         fu = feval(func,x1+u*x2,varargin{:});
      end
      if (~klaar),
         ax = bx;
         bx = cx;
         cx = u;
         fa = fb;
         fb = fc;
         fc = fu;
      end
end