CGOptimizer.C 4.18 KB
Newer Older
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 <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// CGOptimizer.C
//
////////////////////////////////////////////////////////////////////////////////
#include "CGOptimizer.h"
#include <iostream>
#include <cassert>
21
#include <numeric>
22 23
using namespace std;
////////////////////////////////////////////////////////////////////////////////
24
double CGOptimizer::norm2(std::valarray<double>& v)
25
{
26
  return inner_product(&v[0],&v[v.size()],&v[0],0.0);
27 28
}
////////////////////////////////////////////////////////////////////////////////
29
void CGOptimizer::compute_xp(const valarray<double>& x, const double f,
30
                             valarray<double>& g, valarray<double>& xp)
31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
{
  // Use the function value f and the gradient g at x to generate a new point xp
  // using the Fletcher-Powell CG algorithm
  // return xp=x if the 2-norm of g is smaller than tol
  const double tol = 1.0e-18;

  assert(x.size()==n_ && g.size()==n_ && xp.size()==n_);

  double fp;
  // define the descent direction
  if ( first_step_ )
  {
    p_ = -g;

    x0_ = x;
    f0_ = f;

    g0norm2_ = norm2(g);
    if ( g0norm2_ < tol )
    {
      xp = x;
      return;
    }
    fp = -g0norm2_;
    fp0_ = fp;
    linmin_.reset();
    alpha_ = linmin_.next_alpha(alpha_,f,fp);
    if ( debug_print )
59
      cout << "  CGOptimizer: first_step: alpha=" << alpha_
60
           << " f=" << f << " fp=" << fp << endl;
61

62 63 64 65 66 67 68
    xp = x0_ + alpha_ * p_;
    first_step_ = false;
  }
  else
  {
    // fp: derivative along the current descent direction p_
    // fp = df(x0+alpha*p)/dalpha at x
69
    fp = inner_product(&g[0],&g[g.size()],&p_[0],0.0);
70 71
    double new_alpha = linmin_.next_alpha(alpha_,f,fp);
    if ( debug_print )
72
      cout << "  CGOptimizer: alpha=" << alpha_
73 74 75 76 77 78
           << " f=" << f << " fp=" << fp << endl;

    if ( linmin_.fail() )
    {
      // line minimization failed
      if ( debug_print )
79
        cout << "  CGOptimizer: line minimization failed" << endl;
80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97

      // restart from current point
      p_ = -g;

      x0_ = x;
      f0_ = f;

      g0norm2_ = norm2(g);
      if ( g0norm2_ < tol )
      {
        xp = x;
        return;
      }
      fp = -g0norm2_;
      fp0_ = fp;
      linmin_.reset();
      alpha_ = linmin_.next_alpha(alpha_,f,fp);
      if ( debug_print )
98
        cout << "  CGOptimizer: restart after fail: alpha=" << alpha_
99
           << " f=" << f << " fp=" << fp << endl;
100

101 102
      xp = x0_ + alpha_ * p_;
      first_step_ = false;
103
      return;
104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127
    }

    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);
      double beta = norm2(g) / g0norm2_;
      if ( beta_max_ > 0.0 && beta > beta_max_ )
      {
        if ( debug_print )
          cout << "  CGOptimizer: beta exceeds beta_max " << endl;
        beta = min(beta, beta_max_);
      }
      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_
128
      fp0_ = inner_product(&g[0],&g[g.size()],&p_[0],0.0);
129 130 131 132 133 134 135 136
      g0norm2_ = norm2(g);
      fp = fp0_;
      assert(fp<0.0 && "CGOptimizer: p_ not a descent direction");

      // reset the line minimizer
      linmin_.reset();
      alpha_ = linmin_.next_alpha(alpha_,f,fp);
      if ( debug_print )
137
        cout << "  CGOptimizer: restart: alpha=" << alpha_
138 139 140 141 142 143 144 145 146 147
             << " f=" << f << " fp=" << fp << endl;
    }
    else
    {
      // not done with the current descent direction
      alpha_ = new_alpha;
    }
    xp = x0_ + alpha_ * p_;
  }
}