q_newton_3.c

an quasi-newton method used in non-linear programming and especially in optimization

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "math_util.h"
#define MAX_ITT_QN3 5000

void q_newton_3( fun, x, n, tol )
long double (*fun)( long double * );
long double *x, tol;
int n;
{
   long double *grad0, *grad1, *h, *x0, *x1, f0, f1, alpha, err1, gfh;
   long double *gamma, *d, dgamma, dgrad, gg, a, b, ggrad;
   int i, j;
   
   long double line_search_newt( long double(*)(), long double, long double *,
    						long double *, long double *,
                             long double *, int, long double );
   void fd_grad( long double *, long double(*)(), long double, long double *, int );

   if( tol <= 0.0 ) tol = 1.0e-4;
   if( n <= 0 ) return;
   
   gamma = a1d_allo_dbl( n );
   d = a1d_allo_dbl( n );
   grad0 = a1d_allo_dbl( n );
   grad1 = a1d_allo_dbl( n );
   h = a1d_allo_dbl( n );
   x0 = a1d_allo_dbl( n );
   x1 = a1d_allo_dbl( n );
   
   vec_copy( x0, x, n );
   
   f0 = (*fun)( x0 );
   fd_grad( grad0, fun, f0, x0, n );
   
   vec_copy( h, grad0, n );
   for( j = 0; j < n; j++ )
      h[j] = -h[j];

   gfh = 2.0*tol;
   
   for( i = 0; i < MAX_ITT_QN3; i++ ) {
      printf(" itt = %d, f = %Le\n",i,f0);
      err1 = enorm( h, n );
      if( err1 <= tol || gfh <= tol  ) {
         vec_copy( x, x0, n );
         a1d_free_dbl( x0 );
         a1d_free_dbl( x1 );
         a1d_free_dbl( h );
         a1d_free_dbl( d );
         a1d_free_dbl( gamma );
         a1d_free_dbl( grad0 );
         a1d_free_dbl( grad1 );
         printf("***  Norm( h ) <= eps\n");
         return;
      }
               
      alpha = line_search_newt( fun, f0, x0, &f1, x1, h, n, tol );
      
      for( j = 0; j < n; j++ ) {
         d[j] = alpha*h[j];
      }
 /*
      if( fabs(f1-f0) <= tol ) {
         vec_copy( x, x1, n );
         a1d_free_dbl( x0 );
         a1d_free_dbl( x1 );
         a1d_free_dbl( h );
         a1d_free_dbl( d );
         a1d_free_dbl( gamma );
         a1d_free_dbl( grad0 );
         a1d_free_dbl( grad1 );
         printf("***  |f(x1)-f(x0)| <= eps\n");
         return;
      }
*/     
      fd_grad( grad1, fun, f1, x1, n );
      
      for( j = 0; j < n; j++ ) {
         gamma[j] = grad1[j]-grad0[j];
      }
             
      f0 = f1;
      vec_copy( x0, x1, n );
      vec_copy( grad0, grad1, n );
/*
   Update h
*/
      dgamma = 0.0;
      gg = 0.0;
      dgrad = 0.0;
      ggrad = 0.0;
      for( j = 0; j < n; j++ ) {
         dgamma += d[j]*gamma[j];
         gg += gamma[j]*gamma[j];
         dgrad += d[j]*grad0[j];
         ggrad += gamma[j]*grad0[j];
      }
      a = (1.0+gg/dgamma)*dgrad/dgamma - ggrad/dgamma;
      b = dgrad/dgamma;
      for( j = 0; j < n; j++ )
         h[j] = -(grad0[j] + a*d[j] - b*gamma[j]);
     
      /*gfh = 0.0;for(j = 0;j<n;j++) gfh+=grad0[j]*h[j];gfh=fabs(gfh);*/
      gfh = 0.0;for(j = 0;j<n;j++) gfh+=d[j]*d[j]; gfh = sqrt(gfh);
         
   }
   printf("Error: Too many Quasi-Newton iterations\n");
   vec_copy( x, x0, n );
   a1d_free_dbl( x0 );
   a1d_free_dbl( x1 );
   a1d_free_dbl( h );
   a1d_free_dbl( d );
   a1d_free_dbl( gamma );
   a1d_free_dbl( grad0 );
   a1d_free_dbl( grad1 );
   return;   
}
/* 
Compute the gradient using finite difference 
*/
void fd_grad( grad, fun, f0, x0, n ) 
long double *grad, (*fun)( long double * ), *x0, f0;
int n;
{
   int i;
   long double f1, h, *x1;
   
   x1 = a1d_allo_dbl( n );
   
   for( i = 0; i < n; i++ ) {
      vec_copy( x1, x0, n );
      if( fabs((double)(FD_TOL*x0[i])) > MACH_EPS )
         h = fabs(FD_TOL*x0[i]);
      else
         h = FD_TOL;
      x1[i] += h;
      f1 = (*fun)( x1 );
      grad[i]  = (f1-f0)/h;
   }
   a1d_free_dbl( x1 );
}