CGOptimizer.C 4.57 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
#include "blas.h"
23 24
using namespace std;
////////////////////////////////////////////////////////////////////////////////
25
void CGOptimizer::compute_xp(const valarray<double>& x, const double f,
26
                             valarray<double>& g, valarray<double>& xp)
27 28
{
  // Use the function value f and the gradient g at x to generate a new point xp
29
  // using the Fletcher-Reeves or Polak-Ribiere CG algorithm
30 31
  // return xp=x if the 2-norm of g is smaller than tol
  const double tol = 1.0e-18;
32
  const int one = 1;
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;
41
    gm_ = g;
42 43 44 45

    x0_ = x;
    f0_ = f;

46
    g0norm2_ = ddot_(&n_,&g[0],&one,&g[0],&one);
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 )
57
      cout << "  CGOptimizer: first_step: alpha=" << alpha_
58
           << " f=" << f << " fp=" << fp << endl;
59

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
67 68
    fp = ddot_(&n_,&g[0],&one,&p_[0],&one);
    alpha_ = linmin_.next_alpha(alpha_,f,fp);
69
    if ( debug_print )
70
      cout << "  CGOptimizer: alpha=" << alpha_
71 72 73 74 75 76
           << " f=" << f << " fp=" << fp << endl;

    if ( linmin_.fail() )
    {
      // line minimization failed
      if ( debug_print )
77
        cout << "  CGOptimizer: line minimization failed" << endl;
78 79 80

      // restart from current point
      p_ = -g;
81
      gm_ = g;
82 83 84 85

      x0_ = x;
      f0_ = f;

86
      g0norm2_ = ddot_(&n_,&g[0],&one,&g[0],&one);
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 )
97
        cout << "  CGOptimizer: restart after fail: alpha=" << alpha_
98
           << " f=" << f << " fp=" << fp << endl;
99

100 101
      xp = x0_ + alpha_ * p_;
      first_step_ = false;
102
      return;
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);
112 113

#if 0
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

122
      if ( beta_max_ > 0.0 && fabs(beta) > beta_max_ )
123 124
      {
        if ( debug_print )
125 126
          cout << "  CGOptimizer: |beta| exceeds beta_max " << endl;
        beta = (beta > 0.0) ? beta_max_ : -beta_max_;
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_
136 137
      fp = ddot_(&n_,&g[0],&one,&p_[0],&one);
      g0norm2_ = ddot_(&n_,&g[0],&one,&g[0],&one);
138 139

      gm_ = g;
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;
149
        gm_ = g;
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;
      }
158 159 160 161 162

      // reset the line minimizer
      linmin_.reset();
      alpha_ = linmin_.next_alpha(alpha_,f,fp);
      if ( debug_print )
163
        cout << "  CGOptimizer: restart: alpha=" << alpha_
164 165 166 167 168
             << " f=" << f << " fp=" << fp << endl;
    }
    xp = x0_ + alpha_ * p_;
  }
}