1
$\begingroup$

I am trying to find the minimum of the function using GNU scientific library, package Multidimensional Minimization. The method I am using is Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm which is implemented in function 

gsl_multimin_fdfminimizer_vector_bfgs2

Unfortunately, after the first iteration with

gsl_multimin_fdfminimizer_iterate

I receive an error with code 27: iteration is not making progress towards solution.

I tried several initial points and different tolerance/step combinations, but still the same error pops up. Could you please tell me what could be the culprit here?

The function I am trying to minimize is Gupta-Potential for Aluminum 13 atoms cluster.

I also noticed that forces calculated on the first step are the SAME as gradients after first interation. Should it be like that? Shouldn't step come into the picture somehow? I use

gsl_multimin_fdfminimizer_set (s, &my_func, x, 0.001, 1e-10);

And function (Gupta potential) does NOT change after iteration (I print it before and after) though gradients are big.

// HERE IS MY GUPTA POTENTIAL//

double
      my_f (const gsl_vector *v, void *params)
     {
       double *p = (double *)params;
        int ii,jj;
        double vv = 0;
        double sum1 = 0;
        double sum2 = 0;

          for(ii=1;ii<=NN;ii++){
          X[ii]=gsl_vector_get (v, 3*ii-3);
          Y[ii]=gsl_vector_get (v, 3*ii-2);
          Z[ii]=gsl_vector_get (v, 3*ii-1 );
         }

        for(ii = 1; ii <= NN; ii++){
                sum1 = 0;
                for(jj = 1; jj <= NN; jj++){

                        if(jj!=ii){
                                sum1 += uA*exp(0 - uP*((rr(ii,jj)/uR0) - 1.));
                        }
                }
                sum2 = 0;
                for(jj = 1; jj <= NN; jj++){
                        if(jj!=ii){
                                sum2 += uXi*uXi*exp(0 - 2*uQ*((rr(ii,jj)/uR0) - 1.));
                        }
                }
                vv += sum1 - sqrt(sum2);
        }
        return vv;
 }

// HERE ARE MY GRADIENT FUNCTIONS//

double dVdr(int ii, int jj){
        double vr;
        double sum1 = 0;
        double sum2 = 0;
        int kk;

        vr = 0. - 2.*uA*(uP/uR0)*exp(0. - uP*((rr(ii,jj)/uR0) - 1.));

        for(kk = 1; kk <= NN; kk++){
                if(kk != ii){
                        sum1 += uXi*uXi*exp(0 - 2*uQ*((rr(ii,kk)/uR0) - 1.));
                }
                if(kk != jj){
                        sum2 += uXi*uXi*exp(0 - 2*uQ*((rr(kk,jj)/uR0) - 1.));
                }
        }

        vr -= 0.5*uXi*uXi*(0. - 2.*uQ/uR0)*exp(0. - 2.*uQ*((rr(ii,jj)/uR0) - 1.))*( (1./sqrt(sum1)) + (1./sqrt(sum2)) );

        return vr;
}

//-----------------------------------------------
 void Force(int jj, double & fx, double & fy, double & fz){
        int mm;
        double dxx, dyy, dzz, r;
        fx = 0;
        fy = 0;
        fz = 0;
   for(mm = 1; mm <= NN; mm++){

    if(mm != jj){
                        dxx = X[mm] - X[jj];
                        dyy = Y[mm] - Y[jj];
                        dzz = Z[mm] - Z[jj];

     r = sqrt(dxx*dxx + dyy*dyy + dzz*dzz);
     fx += dVdr(mm,jj)*(dxx)/r;
     fy += dVdr(mm,jj)*(dyy)/r;
     fz += dVdr(mm,jj)*(dzz)/r;
   }
   }
}
//--------------------------------------------
 void
     my_df (const gsl_vector *v, void *params,
            gsl_vector *df)
     {
       int ii;
       double fx,fy,fz;
       double *p = (double *)params;

          for(ii=1;ii<=NN;ii++){

          X[ii]=gsl_vector_get (v, 3*ii-3);
          Y[ii]=gsl_vector_get (v, 3*ii-2);
          Z[ii]=gsl_vector_get (v, 3*ii-1 );
         }

          for(ii=1;ii<=NN;ii++){
          Force(ii,fx,fy,fz);
          gsl_vector_set(df, 3*ii-3, fx );
          gsl_vector_set(df, 3*ii-2, fy );
          gsl_vector_set(df, 3*ii-1, fz );
         }
     }
//----------------------------------
| cite | improve this question | |
$\endgroup$
  • $\begingroup$ what is the function you are trying to minimize? $\endgroup$ – GertVdE Mar 23 '12 at 15:51
  • $\begingroup$ @GertVdE Thanks for your attention. I updated the question. $\endgroup$ – Vitaliy Kaurov Mar 23 '12 at 16:05
  • $\begingroup$ What values of uR0 are you searching over? If uR0 is very small the problem will be ill-conditioned. Also, you need to post your corresponding gradient evaluation function. $\endgroup$ – Aron Ahmadia Mar 23 '12 at 16:40
  • $\begingroup$ @AronAhmadia uR0 is around 2.6 - so it's not small. Thanks for the comment! $\endgroup$ – Vitaliy Kaurov Mar 23 '12 at 18:16
  • $\begingroup$ @VitaliyKaurov Would you be so kind to tell us whether you solved the issue? $\endgroup$ – GertVdE Mar 26 '12 at 14:34
3
$\begingroup$

Most likely, the programmed gradient does not match the programmed objective.

You can check this by choosing a few random $x, p$ (including your starting point as $x$, if you provided one) and printing for $s_i=10^{20-i}$, $i=0:40$ the quotients

$q_i = (f(x+ s_ip)-f(x))/(s_ig(x)^Tp)$,

where $g$ is your gradient. The sequence of $q$'s must converge to 1, well, almost, as sooner or later numerical instability sets in and forces the $q_i$ to be zero (and somewhat random one or two steps earlier).

If this is not the case,either the routine for $f$ or that for $g$ is in error. You can check the individual components of $g$ by picking $p$ as the $(0,1)$-unit vectors in turn.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ @Vitaliy: My original formula for $q_i$ was wrong - now corrected! $\endgroup$ – Arnold Neumaier Mar 25 '12 at 12:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.