LineMinimizer.C 8.17 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
////////////////////////////////////////////////////////////////////////////////
//
// 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// LineMinimizer.C
//
////////////////////////////////////////////////////////////////////////////////

#include "LineMinimizer.h"
#include <iostream>
#include <cmath>
#include <cassert>
using namespace std;
////////////////////////////////////////////////////////////////////////////////
double LineMinimizer::interpolate(void)
{
27 28 29 30 31
  const double dalpha = alpha_high - alpha_low;
  // use psi'(alpha_low), psi(alpha_low_), psi(alpha_high)

  double new_alpha;
  if ( use_psi )
32
  {
33 34 35 36 37 38
    // use psi
    double psip_low = psip(fp_low);
    double psi_low = psi(alpha_low,f_low);
    double psi_high = psi(alpha_high,f_high);
    new_alpha = alpha_low - 0.5 * ( psip_low * dalpha * dalpha ) /
                ( psi_high - psi_low - psip_low * dalpha );
39 40 41
  }
  else
  {
42 43 44 45 46 47 48 49 50 51 52 53 54 55 56
    // use f
    // quadratic interpolation using f_low, fp_low, f_high
    // new_alpha = alpha_low - 0.5 * ( fp_low * dalpha * dalpha ) /
    //                  ( f_high - f_low - fp_low * dalpha );
    if ( fp_low*fp_high < 0 )
    {
      // secant
      new_alpha = alpha_low - fp_low *
        ( (alpha_high-alpha_low)/(fp_high-fp_low) );
    }
    else
    {
      // midpoint
      new_alpha = 0.5 * (alpha_low+alpha_high);
    }
57
  }
58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77

  if ( debug_print )
  {
    cout << "LineMinimizer: interpolate: [alpha_low,alpha_high] = ["
         << alpha_low << "," << alpha_high << "]" << endl;
    cout << "LineMinimizer: interpolate: f_low, f_high: "
         << f_low << " " << f_high << endl;
    cout << "LineMinimizer: interpolate: fp_low, fp_high: "
         << fp_low << " " << fp_high << endl;
    cout << "LineMinimizer: interpolate: new_alpha: " << new_alpha << endl;
  }
  if ( alpha_low < alpha_high )
  {
    assert ( new_alpha >= alpha_low && new_alpha <= alpha_high );
  }
  else
  {
    assert ( new_alpha <= alpha_low && new_alpha >= alpha_high );
  }
  return new_alpha;
78 79 80 81 82 83 84 85 86 87 88 89
}
////////////////////////////////////////////////////////////////////////////////
double LineMinimizer::next_alpha(double alpha, double f, double fp)
{
  if ( done_ || fail_ )
    return alpha;

  if ( first_use )
  {
    first_use = false;
    f0 = f;
    fp0 = fp;
90 91 92 93 94 95 96
    alpha_low = 0;
    f_low = f0;
    fp_low = fp0;
    alpha_high = alpha_max_;
    assert(alpha_max_ > alpha_start_);
    if ( debug_print )
      cout << "LineMinimizer: first use: f0, fp0: " << f0 << " " << fp0 << endl;
97 98
    return alpha_start_;
  }
99 100 101
  if ( debug_print )
    cout << "LineMinimizer: next_alpha(" << alpha << ","
         << f << "," << fp << ")" << endl;
102 103 104 105

  bool wolfe1 = f < f0 + sigma1_ * alpha * fp0;
  bool wolfe2 = fabs(fp) < sigma2_ * fabs(fp0);

106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
  if ( debug_print )
  {
    cout << "LineMinimizer: wolfe1: f = " << f << endl;
    cout << "LineMinimizer: wolfe1: f0 + sigma1_ * alpha * fp0 = "
         << f0 + sigma1_ * alpha * fp0 << endl;
    cout << "LineMinimizer: wolfe1/wolfe2: " << wolfe1 << "/" << wolfe2 << endl;
  }

  // check if alpha satisfies both wolfe1 and wolfe2 and return
  if ( wolfe1 && wolfe2 )
  {
    done_ = true;
    return alpha;
  }

121 122
  if ( !bracketing )
  {
123 124 125 126
    // we have not entered the bracketing phase yet
    // Enter bracketing mode if condition U1 holds: psi(alpha) > psi(alpha_low)
    // Note: U1 is equivalent to wolfe1(alpha) == false
    if ( psi(alpha,f) > psi(alpha_low,f_low) )
127
    {
128 129 130 131 132 133 134 135 136 137 138 139 140 141 142
      // wolfe1(alpha) == false
      // we can enter the bracketing phase
      // enter the bracketing phase with alpha_low = 0, alpha_high = alpha

      if ( debug_print )
      {
        cout << "LineMinimizer: entering bracketing: wolfe1==false" << endl;
        cout << "LineMinimizer: psi(alpha), psip(alpha):"
             << psi(alpha,f) << " " << psip(fp) << endl;
        cout << "LineMinimizer: psi(alpha_low, psi(alpha_high):"
             << psi(alpha_low,f_low) << " " << psi(alpha_high,f_high) << endl;
        cout << "LineMinimizer: psip(alpha_low, psip(alpha_high):"
           << psip(fp_low) << " " << psip(fp_high) << endl;
      }

143
      bracketing = true;
144
      use_psi = true;
145 146 147
      alpha_high = alpha;
      f_high = f;
      fp_high = fp;
148 149 150
      assert(psi(alpha_low,f_low)<=0);
      assert(psip(fp_low)*(alpha_high-alpha_low)<=0);
      assert(psi(alpha_low,f_low)<=psi(alpha_high,f_high));
151 152 153
      return interpolate();
    }

154 155 156
    // check if wolfe1 is ok and f'(alpha)(alpha-alpha_low)  > 0
    // Need to test only the second inequality at this point
    if ( fp*(alpha-alpha_low) > 0 )
157
    {
158 159 160
      // enter bracketing mode with alpha_high = alpha_low, alpha_low = alpha
      if ( debug_print )
        cout << "LineMinimizer: entering bracketing: case U3" << endl;
161 162

      bracketing = true;
163 164 165 166
      use_psi = false;
      alpha_high = alpha_low;
      f_high = f_low;
      fp_high = fp_low;
167 168 169
      alpha_low = alpha;
      f_low = f;
      fp_low = fp;
170 171 172
      assert(psi(alpha_low,f_low)<=0);
      assert(psip(fp_low)*(alpha_high-alpha_low)<=0);
      assert(psi(alpha_low,f_low)<=psi(alpha_high,f_high));
173 174 175
      return interpolate();
    }

176 177
    // Condition U2 holds
    // increase alpha
178

179 180 181 182 183 184 185
    if ( debug_print )
      cout << "LineMinimizer: U2, increase alpha" << endl;

    const double delta = 1.1;
    double new_alpha = std::min(alpha+delta*(alpha-alpha_low), alpha_max_);
    if ( new_alpha == alpha_max_ )
      done_ = true;
186 187 188 189
    return new_alpha;
  }
  else
  {
190 191
    // we are already in bracketing mode
    nstep_++;
192
    if ( nstep_max_ > 0 && nstep_ > nstep_max_ )
193
    {
194 195
      if ( debug_print )
        cout << "LineMinimizer: fail, nstep_max" << endl;
196

197
      fail_ = true;
198 199 200
      return alpha;
    }

201 202 203 204 205 206 207 208 209 210 211 212 213 214 215
    if ( debug_print )
    {
      cout << "LineMinimizer: bracketing mode: [alpha_low,alpha_high] = ["
           << alpha_low << "," << alpha_high << "]" << endl;
      cout << "LineMinimizer: bracketing mode: f_low, f_high: "
           << f_low << " " << f_high << endl;
      cout << "LineMinimizer: bracketing mode: fp_low, fp_high: "
           << fp_low << " " << fp_high << endl;
      cout << "LineMinimizer: bracketing mode: psi(alpha), psip(alpha):"
           << psi(alpha,f) << " " << psip(fp) << endl;
      cout << "LineMinimizer: bracketing mode: psi(alpha_low, psi(alpha_high):"
           << psi(alpha_low,f_low) << " " << psi(alpha_high,f_high) << endl;
      cout << "LineMinimizer: bracketing mode: psip(alpha_low, psip(alpha_high):"
           << psip(fp_low) << " " << psip(fp_high) << endl;
    }
216

217 218
    // check U1: psi(alpha) > psi(alpha_low)
    if ( psi(alpha,f) > psi(alpha_low,f_low) )
219
    {
220 221 222 223 224 225 226 227
      if ( debug_print )
        cout << "LineMinimizer: bracketing, U1" << endl;
      alpha_high = alpha;
      f_high = f;
      fp_high = fp;
      assert(psi(alpha_low,f_low)<=0);
      assert(psip(fp_low)*(alpha_high-alpha_low)<=0);
      assert(psi(alpha_low,f_low)<=psi(alpha_high,f_high));
228 229 230 231
      return interpolate();
    }
    else
    {
232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261
      // at this point psi(alpha) <= psi(alpha_low)
      // test condition U2: psi'(alpha)*(alpha_low-alpha) > 0
      if ( psip(fp)*(alpha_low-alpha) > 0 )
      {
        if ( debug_print )
          cout << "LineMinimizer: bracketing, U2" << endl;
        alpha_low = alpha;
        f_low = f;
        fp_low = fp;
        assert(psi(alpha_low,f_low)<=0);
        assert(psip(fp_low)*(alpha_high-alpha_low)<=0);
        assert(psi(alpha_low,f_low)<=psi(alpha_high,f_high));
        return interpolate();
      }
      else
      {
        if ( debug_print )
          cout << "LineMinimizer: bracketing, U3" << endl;
        alpha_high = alpha_low;
        f_high = f_low;
        fp_high = fp_low;
        alpha_low = alpha;
        f_low = f;
        fp_low = fp;
        assert(psi(alpha_low,f_low)<=0);
        assert(psip(fp_low)*(alpha_high-alpha_low)<=0);
        assert(psi(alpha_low,f_low)<=psi(alpha_high,f_high));
        use_psi = false;
        return interpolate();
      }
262 263 264
    }
  }
}