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