1 // Tony : lambda=0.25 ou lambda=0.025 ???
2 // cf. CG.hpp
3 //
4 // The minimizer along the search direction is returned by the function.
5 //
6 /*
7 CubicLineSearch() class implements the efficient line search procedure
8 described in Dennis and Schbabel's book entitled "Numerical
9 Methods for Unconstrained and Nonlinear Equations. The
10 objective is to perform a unidimensional search for a minimum
11 point along a specified direction in a multidimensional
12 space.
13 */
14 // This procedure seems to be fairly robust. It has worked for a fairly broad class of
15 // problems from optimization of standard test functons in optimization theory and
16 // to hard geophysical problems as stacking power optimization and amplitude
17 // seismogram inversion
18
19 #ifndef CUBIC_LINE_SEARCH_HH
20 #define CUBIC_LINE_SEARCH_HH
21
22 #include "LineSearch.hpp"
23 #include <limits.h>
24 #include <float.h>
25
26 template< class LS >
27 class CubicLineSearch : public LS {
28 typedef typename LS::Real Real;
29 typedef typename LS::Param Param;
30 typedef typename LS::Vect Vect;
31 typedef typename LS::VMat VMat;
32 typedef LS LineSearch;
33 typedef tNRJ< Param, Vect, VMat, Real > NRJ;
34
35 public:
36 // a constructor with the default delta
37 CubicLineSearch(NRJ* f, int iter);
38 // a constructor with the specified delta
39 CubicLineSearch(NRJ* f, int iter, Vect* delta);
40 /* The parameter $delta$ is not used by the line search
41 itself. Rather it is used in the numerical computation
42 of the derivatives using centered differences. For
43 example the derivative of f(x) at the point x0 would be
44 given by
45 \[f(x0 - delta) - f(x0 + delta) / 2 * delta\]
46 */
47 ~CubicLineSearch( );
48
49 // search for the minimum model along a 1-D direction
50 Param search( /// initial model for the line search
51 const Param& m0,
52 /// search direction
53 Vect& direction,
54 /// product of search direction and gradient
55 Real descent,
56 /// a parameter
57 double lambda);
58 /* The parameter $lambda$ controls the accuraccy of the
59 line search. $lambda = .25$ is a good choice.
60 The minimizer along the search direction is returned
61 by the function. */
62
63 // lambda=0.25 ou lambda=0.025 ???
64 // cf. CG.hpp
65 };
66
67 template< class LS >
CubicLineSearch(NRJ * f,int it)68 CubicLineSearch< LS >::CubicLineSearch(NRJ* f, int it) : LS(f) {
69 this->iterMax = it;
70 }
71
72 template< class LS >
CubicLineSearch(NRJ * f,int it,Vect * interval)73 CubicLineSearch< LS >::CubicLineSearch(NRJ* f, int it, Vect* interval) : LS(f, interval) {
74 this->iterMax = it;
75 }
76
77 template< class LS >
~CubicLineSearch()78 CubicLineSearch< LS >::~CubicLineSearch( ) {}
79
80 // Code for the Cubic Line Search
81 template< class LS >
search(const Param & current_solution,Vect & p,Real slope,double lambda)82 typename CubicLineSearch< LS >::Param CubicLineSearch< LS >::search(const Param& current_solution,
83 Vect& p, Real slope,
84 double lambda) {
85
86 int tst = 0; // useful vars
87 Real alpha = 0., alpha_prev = 0;
88 Real fprev = 0;
89 Real new_m = 0, old_m = 0;
90 Param new_solution(current_solution);
91 cout << " search " << p.max( ) << endl;
92 assert(p.max( ) < 1e100);
93 old_m = this->nrj->getVal(current_solution); // Evaluation at the current solution
94 this->iterNum = 0;
95 this->iterNum++; // iteration counter
96 alpha = 1.; // updating step
97
98 new_solution = this->update(current_solution, 1, alpha, p); // updating
99 new_m = this->nrj->getVal(new_solution); // Evaluation at the
100 // new solution
101 this->iterNum++;
102
103 // Implementing Goldstein's test for alpha too small
104 while (new_m < old_m + (1. - lambda) * alpha * slope && this->iterNum < this->iterMax) {
105 alpha *= 3;
106 new_solution = this->update(current_solution, 1, alpha, p);
107 new_m = this->nrj->getVal(new_solution);
108 this->iterNum++;
109 }
110 if (this->iterNum == this->iterMax) cerr << "Alpha overflowed! \n";
111
112 // Armijo's test for alpha too large
113 alpha_prev = alpha; // H.L. Deng, 6/13/95
114 while (new_m > old_m + lambda * alpha * slope && this->iterNum < this->iterMax) {
115 Real alpha2 = alpha * alpha, alpha_tmp = 0.;
116 Real f1 = new_m - old_m - slope * alpha;
117 Real a = 0, b = 0;
118 Real c = 0, cm11 = 0, cm12 = 0, cm21 = 0, cm22 = 0;
119 Real disc = 0;
120
121 if (tst == 0) {
122 alpha_tmp = -slope * alpha2 / (f1 * 2.);
123 // tentative alpha
124
125 tst = 1;
126 } else {
127 Real alpha_prev2 = alpha_prev * alpha_prev;
128 Real f2 = fprev - old_m - alpha_prev * slope;
129
130 c = 1. / (alpha - alpha_prev);
131 cm11 = 1. / alpha2;
132 cm12 = -1. / alpha_prev2;
133 cm21 = -alpha_prev / alpha2;
134 cm22 = alpha / alpha_prev2;
135
136 a = c * (cm11 * f1 + cm12 * f2);
137 b = c * (cm21 * f1 + cm22 * f2);
138 disc = b * b - 3. * a * slope;
139
140 if ((Abs(a) > FLT_MIN) && (disc > FLT_MIN))
141 alpha_tmp = (-b + sqrt(disc)) / (3. * a);
142 else
143 alpha_tmp = slope * alpha2 / (2. * f1);
144
145 if (alpha_tmp >= .5 * alpha) alpha_tmp = .5 * alpha;
146 }
147 alpha_prev = alpha;
148 fprev = new_m;
149
150 if (alpha_tmp < .1 * alpha)
151 alpha *= .1;
152 else
153 alpha = alpha_tmp;
154
155 new_solution = this->update(current_solution, 1, alpha, p);
156 new_m = this->nrj->getVal(new_solution);
157 this->iterNum++;
158 }
159 if (this->iterNum == this->iterMax) {
160 cerr << "Alpha underflowed! \n";
161 cerr << this->iterMax;
162 }
163
164 this->value = new_m;
165 this->appendSearchNumber( );
166 return (new_solution); // # of iterations
167 }
168
169 #endif
170