1 #ifndef STAN_OPTIMIZATION_BFGS_LINESEARCH_HPP
2 #define STAN_OPTIMIZATION_BFGS_LINESEARCH_HPP
4 #include <boost/math/special_functions/fpclassify.hpp>
12 namespace optimization {
34 template<
typename Scalar>
36 const Scalar &x1,
const Scalar &f1,
const Scalar &df1,
37 const Scalar &loX,
const Scalar &hiX) {
38 const Scalar c3((-12*f1 + 6*x1*(df0 + df1))/(x1*x1*x1));
39 const Scalar c2(-(4*df0 + 2*df1)/x1 + 6*f1/(x1*x1));
40 const Scalar &c1(df0);
42 const Scalar t_s = std::sqrt(c2*c2 - 2.0*c1*c3);
43 const Scalar s1 = - (c2 + t_s)/c3;
44 const Scalar s2 = - (c2 - t_s)/c3;
50 minF = loX*(loX*(loX*c3/3.0 + c2)/2.0 + c1);
54 tmpF = hiX*(hiX*(hiX*c3/3.0 + c2)/2.0 + c1);
61 if (loX < s1 && s1 < hiX) {
62 tmpF = s1*(s1*(s1*c3/3.0 + c2)/2.0 + c1);
70 if (loX < s2 && s2 < hiX) {
71 tmpF = s2*(s2*(s2*c3/3.0 + c2)/2.0 + c1);
102 template<
typename Scalar>
103 Scalar
CubicInterp(
const Scalar &x0,
const Scalar &f0,
const Scalar &df0,
104 const Scalar &x1,
const Scalar &f1,
const Scalar &df1,
105 const Scalar &loX,
const Scalar &hiX) {
106 return x0 +
CubicInterp(df0, x1-x0, f1-f0, df1, loX-x0, hiX-x0);
113 template<
typename FunctorType,
typename Scalar,
typename XType>
114 int WolfLSZoom(Scalar &alpha, XType &newX, Scalar &newF, XType &newDF,
116 const XType &x,
const Scalar &f,
const Scalar &dfp,
117 const Scalar &c1dfp,
const Scalar &c2dfp,
const XType &p,
118 Scalar alo, Scalar aloF, Scalar aloDFp,
119 Scalar ahi, Scalar ahiF, Scalar ahiDFp,
120 const Scalar &min_range) {
121 Scalar d1, d2, newDFp;
127 if (std::fabs(alo-ahi) < min_range)
131 alpha = 0.5*(alo+ahi);
134 d1 = aloDFp + ahiDFp - 3*(aloF-ahiF)/(alo-ahi);
135 d2 = std::sqrt(d1*d1 - aloDFp*ahiDFp);
139 - (ahi - alo) * (ahiDFp + d2 - d1) / (ahiDFp - aloDFp + 2*d2);
140 if (!boost::math::isfinite(alpha) ||
141 alpha < std::min(alo, ahi) + 0.01 * std::fabs(alo - ahi) ||
142 alpha > std::max(alo, ahi) - 0.01 * std::fabs(alo - ahi))
143 alpha = 0.5 * (alo + ahi);
146 newX = x + alpha * p;
147 while (func(newX, newF, newDF)) {
148 alpha = 0.5 * (alpha + std::min(alo, ahi));
149 if (std::fabs(std::min(alo, ahi) - alpha) < min_range)
151 newX = x + alpha * p;
153 newDFp = newDF.dot(p);
154 if (newF > (f + alpha * c1dfp) || newF >= aloF) {
159 if (std::fabs(newDFp) <= -c2dfp)
161 if (newDFp*(ahi-alo) >= 0) {
221 template<
typename FunctorType,
typename Scalar,
typename XType>
224 XType &x1, Scalar &f1, XType &gradx1,
226 const XType &x0,
const Scalar &f0,
const XType &gradx0,
227 const Scalar &c1,
const Scalar &c2,
228 const Scalar &minAlpha) {
229 const Scalar dfp(gradx0.dot(p));
230 const Scalar c1dfp(c1*dfp);
231 const Scalar c2dfp(c2*dfp);
233 Scalar alpha0(minAlpha);
234 Scalar alpha1(alpha);
237 XType prevDF(gradx0);
241 int retCode = 0, nits = 0, ret;
244 x1.noalias() = x0 + alpha1 * p;
245 ret = func(x1, f1, gradx1);
247 alpha1 = 0.5 * (alpha0 + alpha1);
250 newDFp = gradx1.dot(p);
251 if ((f1 > f0 + alpha * c1dfp) || (f1 >= prevF && nits > 0)) {
252 retCode = WolfLSZoom(alpha, x1, f1, gradx1,
256 alpha0, prevF, prevDFp,
261 if (std::fabs(newDFp) <= -c2dfp) {
266 retCode = WolfLSZoom(alpha, x1, f1, gradx1,
271 alpha0, prevF, prevDFp,
278 std::swap(prevDF, gradx1);
Probability, optimization and sampling library.
int WolfeLineSearch(FunctorType &func, Scalar &alpha, XType &x1, Scalar &f1, XType &gradx1, const XType &p, const XType &x0, const Scalar &f0, const XType &gradx0, const Scalar &c1, const Scalar &c2, const Scalar &minAlpha)
Perform a line search which finds an approximate solution to: satisfying the strong Wolfe conditions...
Scalar CubicInterp(const Scalar &df0, const Scalar &x1, const Scalar &f1, const Scalar &df1, const Scalar &loX, const Scalar &hiX)
Find the minima in an interval [loX, hiX] of a cubic function which interpolates the points...