CGOptimizer.C 4.57 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 Sep 28, 2011 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 Sep 28, 2011 67 68  fp = ddot_(&n_,&g[0],&one,&p_[0],&one); 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 Sep 28, 2011 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 115 116 117 118 119 120 121  // Fletcher-Reeves double beta = ddot_(&n_,&g[0],&one,&g[0],&one) / g0norm2_; #else // Polak-Ribiere double beta = (ddot_(&n_,&g[0],&one,&g[0],&one)- ddot_(&n_,&gm_[0],&one,&g[0],&one)) / g0norm2_; #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 Sep 28, 2011 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 154 155 156 157  x0_ = x; f0_ = f; g0norm2_ = ddot_(&n_,&g[0],&one,&g[0],&one); fp = -g0norm2_; fp0_ = fp; }  Francois Gygi committed Aug 30, 2011 158 159 160 161 162  // reset the line minimizer linmin_.reset(); alpha_ = linmin_.next_alpha(alpha_,f,fp); if ( debug_print )  Francois Gygi committed Sep 07, 2011 163  cout << " CGOptimizer: restart: alpha=" << alpha_  Francois Gygi committed Aug 30, 2011 164 165 166 167 168  << " f=" << f << " fp=" << fp << endl; } xp = x0_ + alpha_ * p_; } }