1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) 2011 The Regents of the University of California
4 //
5 // This file is part of Qbox
6 //
7 // Qbox is distributed under the terms of the GNU General Public License
8 // as published by the Free Software Foundation, either version 2 of
9 // the License, or (at your option) any later version.
10 // See the file COPYING in the root directory of this distribution
11 // or <http://www.gnu.org/licenses/>.
12 //
13 ////////////////////////////////////////////////////////////////////////////////
14 //
15 // LineMinimizer.cpp
16 //
17 ////////////////////////////////////////////////////////////////////////////////
18 
19 #include "LineMinimizer.h"
20 #include <iostream>
21 #include <cmath>
22 #include <cassert>
23 using namespace std;
24 ////////////////////////////////////////////////////////////////////////////////
interpolate(void)25 double LineMinimizer::interpolate(void)
26 {
27   const double dalpha = alpha_high - alpha_low;
28   // use psi'(alpha_low), psi(alpha_low_), psi(alpha_high)
29 
30   double new_alpha;
31   if ( use_psi )
32   {
33     // use psi
34     double psip_low = psip(fp_low);
35     double psi_low = psi(alpha_low,f_low);
36     double psi_high = psi(alpha_high,f_high);
37     new_alpha = alpha_low - 0.5 * ( psip_low * dalpha * dalpha ) /
38                 ( psi_high - psi_low - psip_low * dalpha );
39   }
40   else
41   {
42     // use f
43     // quadratic interpolation using f_low, fp_low, f_high
44     // new_alpha = alpha_low - 0.5 * ( fp_low * dalpha * dalpha ) /
45     //                  ( f_high - f_low - fp_low * dalpha );
46     if ( fp_low*fp_high < 0 )
47     {
48       // secant
49       new_alpha = alpha_low - fp_low *
50         ( (alpha_high-alpha_low)/(fp_high-fp_low) );
51     }
52     else
53     {
54       // midpoint
55       new_alpha = 0.5 * (alpha_low+alpha_high);
56     }
57   }
58 
59   if ( debug_print )
60   {
61     cout << "LineMinimizer: interpolate: [alpha_low,alpha_high] = ["
62          << alpha_low << "," << alpha_high << "]" << endl;
63     cout << "LineMinimizer: interpolate: f_low, f_high: "
64          << f_low << " " << f_high << endl;
65     cout << "LineMinimizer: interpolate: fp_low, fp_high: "
66          << fp_low << " " << fp_high << endl;
67     cout << "LineMinimizer: interpolate: new_alpha: " << new_alpha << endl;
68   }
69   if ( alpha_low < alpha_high )
70   {
71     assert ( new_alpha >= alpha_low && new_alpha <= alpha_high );
72   }
73   else
74   {
75     assert ( new_alpha <= alpha_low && new_alpha >= alpha_high );
76   }
77   return new_alpha;
78 }
79 ////////////////////////////////////////////////////////////////////////////////
next_alpha(double alpha,double f,double fp)80 double LineMinimizer::next_alpha(double alpha, double f, double fp)
81 {
82   if ( done_ || fail_ )
83     return alpha;
84 
85   if ( first_use )
86   {
87     first_use = false;
88     f0 = f;
89     fp0 = fp;
90     alpha_low = 0;
91     f_low = f0;
92     fp_low = fp0;
93     alpha_high = alpha_max_;
94     assert(alpha_max_ > alpha_start_);
95     if ( debug_print )
96       cout << "LineMinimizer: first use: f0, fp0: " << f0 << " " << fp0 << endl;
97     return alpha_start_;
98   }
99   if ( debug_print )
100     cout << "LineMinimizer: next_alpha(" << alpha << ","
101          << f << "," << fp << ")" << endl;
102 
103   bool wolfe1 = f < f0 + sigma1_ * alpha * fp0;
104   bool wolfe2 = fabs(fp) < sigma2_ * fabs(fp0);
105 
106   if ( debug_print )
107   {
108     cout << "LineMinimizer: wolfe1: f = " << f << endl;
109     cout << "LineMinimizer: wolfe1: f0 + sigma1_ * alpha * fp0 = "
110          << f0 + sigma1_ * alpha * fp0 << endl;
111     cout << "LineMinimizer: wolfe1/wolfe2: " << wolfe1 << "/" << wolfe2 << endl;
112   }
113 
114   // check if alpha satisfies both wolfe1 and wolfe2 and return
115   if ( wolfe1 && wolfe2 )
116   {
117     done_ = true;
118     return alpha;
119   }
120 
121   if ( !bracketing )
122   {
123     // we have not entered the bracketing phase yet
124     // Enter bracketing mode if condition U1 holds: psi(alpha) > psi(alpha_low)
125     // Note: U1 is equivalent to wolfe1(alpha) == false
126     if ( psi(alpha,f) > psi(alpha_low,f_low) )
127     {
128       // wolfe1(alpha) == false
129       // we can enter the bracketing phase
130       // enter the bracketing phase with alpha_low = 0, alpha_high = alpha
131 
132       if ( debug_print )
133       {
134         cout << "LineMinimizer: entering bracketing: wolfe1==false" << endl;
135         cout << "LineMinimizer: psi(alpha), psip(alpha):"
136              << psi(alpha,f) << " " << psip(fp) << endl;
137         cout << "LineMinimizer: psi(alpha_low, psi(alpha_high):"
138              << psi(alpha_low,f_low) << " " << psi(alpha_high,f_high) << endl;
139         cout << "LineMinimizer: psip(alpha_low, psip(alpha_high):"
140            << psip(fp_low) << " " << psip(fp_high) << endl;
141       }
142 
143       bracketing = true;
144       use_psi = true;
145       alpha_high = alpha;
146       f_high = f;
147       fp_high = fp;
148       assert(psi(alpha_low,f_low)<=0);
149       assert(psip(fp_low)*(alpha_high-alpha_low)<=0);
150       assert(psi(alpha_low,f_low)<=psi(alpha_high,f_high));
151       return interpolate();
152     }
153 
154     // check if wolfe1 is ok and f'(alpha)(alpha-alpha_low)  > 0
155     // Need to test only the second inequality at this point
156     if ( fp*(alpha-alpha_low) > 0 )
157     {
158       // enter bracketing mode with alpha_high = alpha_low, alpha_low = alpha
159       if ( debug_print )
160         cout << "LineMinimizer: entering bracketing: case U3" << endl;
161 
162       bracketing = true;
163       use_psi = false;
164       alpha_high = alpha_low;
165       f_high = f_low;
166       fp_high = fp_low;
167       alpha_low = alpha;
168       f_low = f;
169       fp_low = fp;
170       assert(psi(alpha_low,f_low)<=0);
171       assert(psip(fp_low)*(alpha_high-alpha_low)<=0);
172       assert(psi(alpha_low,f_low)<=psi(alpha_high,f_high));
173       return interpolate();
174     }
175 
176     // Condition U2 holds
177     // increase alpha
178 
179     if ( debug_print )
180       cout << "LineMinimizer: U2, increase alpha" << endl;
181 
182     double new_alpha = std::min(alpha+delta_*(alpha-alpha_low), alpha_max_);
183     if ( new_alpha == alpha_max_ )
184       done_ = true;
185     return new_alpha;
186   }
187   else
188   {
189     // we are already in bracketing mode
190     nstep_++;
191     if ( nstep_max_ > 0 && nstep_ > nstep_max_ )
192     {
193       if ( debug_print )
194         cout << "LineMinimizer: fail, nstep_max" << endl;
195 
196       fail_ = true;
197       return alpha;
198     }
199 
200     if ( debug_print )
201     {
202       cout << "LineMinimizer: bracketing mode: [alpha_low,alpha_high] = ["
203            << alpha_low << "," << alpha_high << "]" << endl;
204       cout << "LineMinimizer: bracketing mode: f_low, f_high: "
205            << f_low << " " << f_high << endl;
206       cout << "LineMinimizer: bracketing mode: fp_low, fp_high: "
207            << fp_low << " " << fp_high << endl;
208       cout << "LineMinimizer: bracketing mode: psi(alpha), psip(alpha):"
209            << psi(alpha,f) << " " << psip(fp) << endl;
210       cout << "LineMinimizer: bracketing mode: psi(alpha_low, psi(alpha_high):"
211            << psi(alpha_low,f_low) << " " << psi(alpha_high,f_high) << endl;
212       cout << "LineMinimizer: bracketing mode: psip(alpha_low, psip(alpha_high):"
213            << psip(fp_low) << " " << psip(fp_high) << endl;
214     }
215 
216     // check U1: psi(alpha) > psi(alpha_low)
217     if ( psi(alpha,f) > psi(alpha_low,f_low) )
218     {
219       if ( debug_print )
220         cout << "LineMinimizer: bracketing, U1" << endl;
221       alpha_high = alpha;
222       f_high = f;
223       fp_high = fp;
224       assert(psi(alpha_low,f_low)<=0);
225       assert(psip(fp_low)*(alpha_high-alpha_low)<=0);
226       assert(psi(alpha_low,f_low)<=psi(alpha_high,f_high));
227       return interpolate();
228     }
229     else
230     {
231       // at this point psi(alpha) <= psi(alpha_low)
232       // test condition U2: psi'(alpha)*(alpha_low-alpha) > 0
233       if ( psip(fp)*(alpha_low-alpha) > 0 )
234       {
235         if ( debug_print )
236           cout << "LineMinimizer: bracketing, U2" << endl;
237         alpha_low = alpha;
238         f_low = f;
239         fp_low = fp;
240         assert(psi(alpha_low,f_low)<=0);
241         assert(psip(fp_low)*(alpha_high-alpha_low)<=0);
242         assert(psi(alpha_low,f_low)<=psi(alpha_high,f_high));
243         return interpolate();
244       }
245       else
246       {
247         if ( debug_print )
248           cout << "LineMinimizer: bracketing, U3" << endl;
249         alpha_high = alpha_low;
250         f_high = f_low;
251         fp_high = fp_low;
252         alpha_low = alpha;
253         f_low = f;
254         fp_low = fp;
255         assert(psi(alpha_low,f_low)<=0);
256         assert(psip(fp_low)*(alpha_high-alpha_low)<=0);
257         assert(psi(alpha_low,f_low)<=psi(alpha_high,f_high));
258         use_psi = false;
259         return interpolate();
260       }
261     }
262   }
263 }
264