CGOptimizer.C 4.79 KB
 Francois Gygi committed Aug 30, 2011 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 //////////////////////////////////////////////////////////////////////////////// // // Copyright (c) 2011 The Regents of the University of California // // This file is part of Qbox // // Qbox is distributed under the terms of the GNU General Public License // as published by the Free Software Foundation, either version 2 of // the License, or (at your option) any later version. // See the file COPYING in the root directory of this distribution // or . // //////////////////////////////////////////////////////////////////////////////// // // CGOptimizer.C // //////////////////////////////////////////////////////////////////////////////// #include "CGOptimizer.h" #include #include  Francois Gygi committed Sep 07, 2011 21 #include  Francois Gygi committed Sep 28, 2011 22 #include "blas.h"  Francois Gygi committed Aug 30, 2011 23 24 using namespace std; ////////////////////////////////////////////////////////////////////////////////  Francois Gygi committed Sep 07, 2011 25 void CGOptimizer::compute_xp(const valarray& x, const double f,  Francois Gygi committed Sep 11, 2011 26  valarray& g, valarray& xp)  Francois Gygi committed Aug 30, 2011 27 28 { // Use the function value f and the gradient g at x to generate a new point xp  Francois Gygi committed Sep 29, 2011 29  // using the Fletcher-Reeves or Polak-Ribiere CG algorithm  Francois Gygi committed Aug 30, 2011 30 31  // return xp=x if the 2-norm of g is smaller than tol const double tol = 1.0e-18;  Francois Gygi committed Sep 28, 2011 32  const int one = 1;  Francois Gygi committed Aug 30, 2011 33 34 35 36 37 38 39 40  assert(x.size()==n_ && g.size()==n_ && xp.size()==n_); double fp; // define the descent direction if ( first_step_ ) { p_ = -g;  Francois Gygi committed Sep 29, 2011 41  gm_ = g;  Francois Gygi committed Aug 30, 2011 42 43 44 45  x0_ = x; f0_ = f;  Francois Gygi committed Mar 15, 2012 46  g0norm2_ = ddot(&n_,&g[0],&one,&g[0],&one);  Francois Gygi committed Aug 30, 2011 47 48 49 50 51 52 53 54 55 56  if ( g0norm2_ < tol ) { xp = x; return; } fp = -g0norm2_; fp0_ = fp; linmin_.reset(); alpha_ = linmin_.next_alpha(alpha_,f,fp); if ( debug_print )  Francois Gygi committed Sep 07, 2011 57  cout << " CGOptimizer: first_step: alpha=" << alpha_  Francois Gygi committed Aug 30, 2011 58  << " f=" << f << " fp=" << fp << endl;  Francois Gygi committed Sep 07, 2011 59   Francois Gygi committed Aug 30, 2011 60 61 62 63 64 65 66  xp = x0_ + alpha_ * p_; first_step_ = false; } else { // fp: derivative along the current descent direction p_ // fp = df(x0+alpha*p)/dalpha at x  Francois Gygi committed Mar 15, 2012 67  fp = ddot(&n_,&g[0],&one,&p_[0],&one);  Francois Gygi committed Sep 28, 2011 68  alpha_ = linmin_.next_alpha(alpha_,f,fp);  Francois Gygi committed Aug 30, 2011 69  if ( debug_print )  Francois Gygi committed Sep 07, 2011 70  cout << " CGOptimizer: alpha=" << alpha_  Francois Gygi committed Aug 30, 2011 71 72 73 74 75 76  << " f=" << f << " fp=" << fp << endl; if ( linmin_.fail() ) { // line minimization failed if ( debug_print )  Francois Gygi committed Sep 08, 2011 77  cout << " CGOptimizer: line minimization failed" << endl;  Francois Gygi committed Aug 30, 2011 78 79 80  // restart from current point p_ = -g;  Francois Gygi committed Sep 29, 2011 81  gm_ = g;  Francois Gygi committed Aug 30, 2011 82 83 84 85  x0_ = x; f0_ = f;  Francois Gygi committed Mar 15, 2012 86  g0norm2_ = ddot(&n_,&g[0],&one,&g[0],&one);  Francois Gygi committed Aug 30, 2011 87 88 89 90 91 92 93 94 95 96  if ( g0norm2_ < tol ) { xp = x; return; } fp = -g0norm2_; fp0_ = fp; linmin_.reset(); alpha_ = linmin_.next_alpha(alpha_,f,fp); if ( debug_print )  Francois Gygi committed Sep 07, 2011 97  cout << " CGOptimizer: restart after fail: alpha=" << alpha_  Francois Gygi committed Aug 30, 2011 98  << " f=" << f << " fp=" << fp << endl;  Francois Gygi committed Sep 07, 2011 99   Francois Gygi committed Aug 30, 2011 100 101  xp = x0_ + alpha_ * p_; first_step_ = false;  Francois Gygi committed Sep 08, 2011 102  return;  Francois Gygi committed Aug 30, 2011 103 104 105 106 107 108 109 110 111  } if ( linmin_.done() ) { // wolfe1_ && wolfe2_ are true at alpha_ if ( debug_print ) cout << " CGOptimizer: done with current descent direction" << endl; // define a new descent direction p_ using the Fletcher-Reeves formula assert(g0norm2_ > 0.0);  Francois Gygi committed Sep 28, 2011 112 113  #if 0  Francois Gygi committed Sep 29, 2011 114  // Fletcher-Reeves  Francois Gygi committed Mar 15, 2012 115  double beta = ddot(&n_,&g[0],&one,&g[0],&one) / g0norm2_;  Francois Gygi committed Sep 29, 2011 116 117 #else // Polak-Ribiere  Francois Gygi committed Mar 15, 2012 118 119  double beta = (ddot(&n_,&g[0],&one,&g[0],&one)- ddot(&n_,&gm_[0],&one,&g[0],&one)) / g0norm2_;  Francois Gygi committed Sep 29, 2011 120 121 #endif  Francois Gygi committed Sep 29, 2011 122  if ( beta_max_ > 0.0 && fabs(beta) > beta_max_ )  Francois Gygi committed Aug 30, 2011 123 124  { if ( debug_print )  Francois Gygi committed Sep 29, 2011 125 126  cout << " CGOptimizer: |beta| exceeds beta_max " << endl; beta = (beta > 0.0) ? beta_max_ : -beta_max_;  Francois Gygi committed Aug 30, 2011 127 128 129 130 131 132 133 134 135  } if ( debug_print ) cout << " CGOptimizer: beta = " << beta << endl; p_ = beta * p_ - g; x0_ = x; f0_ = f; // recalculate f0, fp0 // fp0 = d_e / d_alpha in direction pc_  Francois Gygi committed Mar 15, 2012 136 137  fp = ddot(&n_,&g[0],&one,&p_[0],&one); g0norm2_ = ddot(&n_,&g[0],&one,&g[0],&one);  Francois Gygi committed Sep 29, 2011 138 139  gm_ = g;  Francois Gygi committed Sep 28, 2011 140 141 142 143 144 145 146 147 148  fp0_ = fp; if ( fp > 0.0 ) { // p_ is not a descent direction // restart from current point if ( debug_print ) cout << " CGOptimizer: p_ not a descent direction" << endl; p_ = -g;  Francois Gygi committed Sep 29, 2011 149  gm_ = g;  Francois Gygi committed Sep 28, 2011 150 151 152 153  x0_ = x; f0_ = f;  Francois Gygi committed Mar 15, 2012 154  g0norm2_ = ddot(&n_,&g[0],&one,&g[0],&one);  Francois Gygi committed Sep 28, 2011 155 156 157  fp = -g0norm2_; fp0_ = fp; }  Francois Gygi committed Aug 30, 2011 158   Francois Gygi committed Jul 28, 2015 159 160 161 162 163 164  // set the starting alpha of the minimizer to be the current alpha_ if ( alpha_ < linmin_.alpha_max() ) linmin_.set_alpha_start(alpha_); else linmin_.set_alpha_start(0.5*linmin_.alpha_max());  Francois Gygi committed Aug 30, 2011 165 166 167 168  // reset the line minimizer linmin_.reset(); alpha_ = linmin_.next_alpha(alpha_,f,fp); if ( debug_print )  Francois Gygi committed Sep 07, 2011 169  cout << " CGOptimizer: restart: alpha=" << alpha_  Francois Gygi committed Aug 30, 2011 170 171 172 173 174  << " f=" << f << " fp=" << fp << endl; } xp = x0_ + alpha_ * p_; } }