symmetricRankOne.m

Go to the documentation of this file.

%> \file
%>  Algorithm 13.2: SR1 method with trust region region. Implementation of algorithm 13.2 of \cite Bier15-book
%> 
%> @author Michel Bierlaire
%> @date Sat Mar 21 16:35:22 2015
%> @ingroup Algorithms
%> @ingroup chap13

%> Applies SR1 algorithm with trust region to solve \f$\min_x f(x)\f$ where \f$f:\mathbb{R}^n\to\mathbb{R}\f$.  The parameters of the method are taken from \cite ConnGoulToin00 (p. 117). 

%> @note Tested with \ref run0508sr1.m
%> @note Tested with \ref runRosenbrockSr1.m
%> @note Calls \ref dogleg
%> @note Calls \ref truncatedConjugateGradient

%> @param obj the name of the Octave function defining f(x) and its derivatives
%> @param x0 the starting point
%> @param delta0 radius of the initial trust region
%> @param eps  algorithm stops if \f$\|F(x)\| \leq \varepsilon \f$. 
%> @param tr method to solve the trust region subproblem. If 0, the dogleg method is used. If different from 0, the truncated conjugate gradient is used (default: 0). 
%> @param maxiter maximum number of iterations (Default: 100)
%> @return [solution,iteres,niter] 
%> @return solution: local minimum of the function
%> @return iteres: sequence of iterates generated by the algorithm. It contains n+2 columns. Columns 1:n contains the value of the current iterate. Column n+1 contains the value of the objective function. Column n+2 contains the value of the norm of the gradient. It contains maxiter rows, but only the first niter rows are meaningful.
%> @return niter: total number of iterations
function [solution,iteres,niter] = symmetricRankOne(obj,x0,delta0,eps, tr=0,maxiter=100)
  addpath("../chap12") ;
    STSR1 = "" ;
  iteres = zeros(1+ maxiter,4) ;

  eta1 = 0.01 ;
  eta2 = 0.9 ;
  k=0 ;
  xk = x0 ;
  n = size(x0,1) ;
  H = eye(n,n) ;
  [f,g] = feval(obj,xk)  ;
  iteres(1,:) = [xk'  f  norm(g) ] ;
  k = 0 ;
  delta = delta0 ;

# Store the radius of the trust region to generate figures like Fig 13.8.
  [fff,msg] = fopen("radius.dat","w") ;
  fprintf(fff,"%d  %e\n",k,delta) ;
  do
  k=k+1 ;

  if (tr == 0)
    [step,type] = dogleg(g,H,delta) ;
  else 
    [step,type] = truncatedConjugateGradient(g,H,delta) ;
  endif
  [fc,gc] = feval(obj,xk+step) ;
  y = gc - g ;
  num = f - fc;
  denom = -step'*g - 0.5 * step' * H * step ;
  rho = num / denom ;
  fprintf(fff,"%e %e %e ",num,denom,rho) ;
  if (rho < eta1)
    delta = norm(step) / 2.0 ;
    status = "- " ;
  else
    xk = xk + step ;
    f = fc ;
    g = gc ;
    if (rho >= eta2)
      delta = 2 * delta ;
      status = "++" ;
    else
      status = "+ " ;
    endif
    term = y - H * step ;
    sr1denom = step' * (term) ;
    STSR1 = "" ;
    if (abs(sr1denom) >= 1.0e-8 * norm(step) * norm(term)) 
      STSR1 = "SR1" ;
      H = H + (term * term') / sr1denom ;
    endif
  endif

  iteres(k+1,:) = [xk' f  norm(g) ] ;
  fprintf(fff,"%d  %e\n",k,delta) ;
  until (norm(g) <= eps || k >= maxiter)

  fclose(fff) ;
  solution = xk ;
  niter = k ;