inex_lsearch.m

ONE-DIMENSIONAL_OPTIMIZATION,主要用于最优化方法的设计以及最优化问题的求解

%  Program: inex_lsearch.m
%  Title: Inexact Line Search
%  Description: Implements Fletcher's inexact line search described in 
%  Algorithm 4.6.
%  Theory: See Practical Optimization Sec. 4.8
%  Input:
%    x:  initial point
%    s:  search direction
%    F:  objective function to be minimized along the direction of s  
%    G:  gradient of objective function F 
%   p1:  internal parameters that are required for the implementation of
%        the line search regardless of the application at hand.
%        It is a string (e.g. 'rho = 0.1') and can be a combination
%        several internal parameters (e.g., 'rho = 0.25; sigma=0.5').
%        Useful p1's include:                            default value
%        'rho=   ' defines right bracket                      0.1
%        'sigma= ' defines left bracket  (sigma >= rho)       0.1
%        'tau= '   defines minimum step for sectioning        0.1
%        'chi='                                               0.75
%   p2:  user-defined parameter vector. Note that p2 must be a vector
%        with all components numerically specified. The components 
%        of p2 must appear in the same order as they appear 
%        in function F and gradient G. For example, if p2=[a b], 
%        then F.m and G.m must be in the form of function z = F(x,p2)
%        and function z = G(x,p2)
%  Output:
%     z: acceptable value of alpha
%  Example 1: 
%  Perform inexact line search using the Himmelblau function
%     f(x1,x2) = (x1^2 + x2 - 11)^2 + (x1 + x2^2 - 7)^2
%  starting from point xk = [6 6]' along search direction s = [-1 -1]'
%  using the default parameter values. 
%  Solution: 
%  Execute the command
%     a1 = inex_lsearch([6 6]',[-1 -1]','f_himm','g_himm')
%  Example 2: 
%  Perform inexact line search using the Himmelblau funtion
%  starting from point xk = [6 6]' along search direction s = [-1 -1]'
%  with sigma = 0.5.
%  Solution: 
%  Execute the command
%     a2 = inex_lsearch([6 6]',[-1 -1]','f_himm','g_himm','sigma = 0.5')
%  Example 3: 
%  Perform inexact line search using the paramerized Himmelblau funtion
%     f(x1,x2,a,b) = (x1^2 + x2 - a^2)^2 + (x1 + x2^2 - b^2)^2
%  starting from point xk = [6 6]' along search direction s = [-1 -1]',
%  with parameters a = 3.2 and b = 2.6.
%  Solution: 
%  Execute the command
%     a3 = inex_lsearch([6 6]',[-1 -1]','f_himm_p','g_himm_p',[3.2 2.6])
%  Notes:
%  1. Command 
%       z = inex_lsearch(xk,s,F,G,p1,p2)
%     adds a new function inex_lsearch to MATLAB's vocabulary.
%  2. Do not use a semicolon in commands 
%      a1 = inex_lsearch([6 6]',[-1 -1]','f_himm','g_himm')
%      a2 = inex_lsearch([6 6]',[-1 -1]','f_himm','g_himm','sigma = 0.5')
%      a3 = inex_lsearch([6 6]',[-1 -1]','f_himm_p','g_himm_p',[3.2 2.6])
%    otherwise the acceptable value of alpha will not be displayed.
% 3. f_himm and g_himm are user defined functions implemented in m-files  
%    f_himm.m and g_himm.m, respectively, and are used to evaluate the 
%    Himmelblau function and its derivative. Similarly, f_himm_p and 
%    g_himm_p are user defined functions for the parameterized Himmelblau
%    function and its gradient.
%==========================================================================
function z = inex_lsearch(xk,s,F,G,p1,p2)
k = 0;
m = 0;
tau = 0.1;
chi = 0.75;
rho = 0.1;
sigma = 0.1;
mhat = 400;
epsilon = 1e-10;
xk = xk(:);
s = s(:);
parameterstring ='';
% evaluate given parameters:
  if nargin > 4,
    if isstr(p1),
      eval([p1 ';']);
    else
      parameterstring = ',p1';
    end
  end
  if nargin > 5,
    if isstr(p2),
      eval([p2 ';']);
    else
      parameterstring = ',p2';
    end
  end
% compute f0 and g0
  eval(['f0 = ' F '(xk' parameterstring ');']);
  eval(['gk = ' G '(xk' parameterstring ');']);
  m = m+2;
  deltaf0 = f0;
% step 2 Initialize line search
  dk = s;
  aL = 0;
  aU = 1e99;
  fL = f0;
  dfL = gk'*dk;
  if abs(dfL) > epsilon,
    a0 = -2*deltaf0/dfL;
  else
    a0 = 1;
 end
 if ((a0 <= 1e-9)|(a0 > 1)),
    a0 = 1;
 end
%step 3
 while 1,
    deltak = a0*dk;
    eval(['f0 = ' F '(xk+deltak' parameterstring ');']);
    m = m + 1;
%step 4
    if ((f0 > (fL + rho*(a0 - aL)*dfL)) & (abs(fL - f0) > epsilon) & (m < mhat))
      if (a0 < aU)
        aU = a0;
      end
      % compute a0hat using equation 7.65
      a0hat = aL + ((a0 - aL)^2*dfL)/(2*(fL - f0 + (a0 - aL)*dfL));
      a0Lhat = aL + tau*(aU - aL);
      if (a0hat < a0Lhat)
        a0hat = a0Lhat;
      end
      a0Uhat = aU - tau*(aU - aL);
      if (a0hat > a0Uhat)
        a0hat = a0Uhat;
      end
      a0 = a0hat;
    else
      eval(['gtemp =' G '(xk+a0*dk' parameterstring ');']);
      df0 = gtemp'*dk;
      m = m + 1;
      % step 6
      if (((df0 < sigma*dfL) & (abs(fL - f0) > epsilon) & (m < mhat) & (dfL ~= df0)))
        deltaa0 = (a0 - aL)*df0/(dfL - df0);
        if (deltaa0 <= 0)
          a0hat = 2*a0;
        else
          a0hat = a0 + deltaa0;
        end
        a0Uhat = a0 + chi*(aU - a0);
        if (a0hat > a0Uhat)
          a0hat = a0Uhat;
        end
        aL = a0;
        a0 = a0hat;
        fL = f0;
        dfL = df0;
      else
        break;
      end
    end
 end % while 1
 if a0 < 1e-5,
    z = 1e-5;
 else
    z = a0;
 end