Commit 8551a20a authored by Jason Rhinelander's avatar Jason Rhinelander

Improvements to single_peak_search

- Added tol_abs parameter which defaults to 1e-20.
- renamed existing tolerance parameter to tol_rel
- Unrolled the recursion into a loop
- Added check and fix for midleft > midright, which can happen when the
  optimum is close to 0.
- Added check for numerical precision limit: if the calculated midpoint
  is not distinguishable from the endpoint, we're as good as it's going
  to get.  (This, combined with the previous problem, could cause an
  infinite loop for an optimum at 0).
parent 5cf27380
......@@ -260,13 +260,21 @@ class Stepper final {
* function value.
* \param left the left edge of the domain to consider
* \param right the right edge of the domain to consider
* \param tolerance the relative size of the domain at which the algorithm stops. In particular,
* the algorithm stops once \f$\frac{right - left}{max\{\|left\|, \|right\|\}} \leq tolerance\f$. The default, if the
* argument is omitted, is \f$10^{-12}\f$.
* \param tol_rel the relative size of the domain at which the algorithm stops. In particular, the
* algorithm stops if \f$\frac{right - left}{max\{\|left\|, \|right\|\}} \leq tol_{rel}\f$. The
* default, if the argument is omitted, is \f$10^{-10}\f$. Note that the algorithm might
* alternatively stop because of the `tol_abs` value. Note also that there is also a numerical
* lower bound on this value: if the midpoint calculation results in a midpoint exactly numerically
* equal to an end point, the algorithm stops.
* \param tol_abs the absolute size of the domain at which the algorithm stops. In particular, the
* algorithm stops if \f$right - left \leq tol_{abs}\f$. Note that the algorithm might
* alternatively stop because of the `tol_rel` value.
*/
double single_peak_search(
const std::function<double(const double &)> &f,
const double &left, const double &right,
const double &tolerance = 1e-12);
double left,
double right,
double tol_rel = 1e-10,
double tol_abs = 1e-20);
}
......@@ -58,47 +58,59 @@ const double sps_left_mid_ = 1.5 - std::sqrt(1.25);
// Relative right midpoint position. Equal to \f$\varphi - 1\f$.
const double sps_right_mid_ = 1 - sps_left_mid_;
double sps_(
const std::function<double(const double &)> &f,
const double &left, const double &midleft, const double &midright, const double &right,
const double &fl, const double &fml, const double &fmr, const double &fr,
const double &tolerance) {
if ((right - left) / std::max(std::fabs(left), std::fabs(right)) <= tolerance) {
// Prefer the end-points for ties (the max might legitimate be an end-point), and prefer
// left over right (for no particularly good reason).
return fl >= fml && fl >= fmr && fl >= fr ? left :
fr >= fmr && fr >= fml ? right :
fml >= fmr ? midleft : midright;
}
if (fml >= fmr) {
// midleft is the higher point, so we can exclude everything right of midright.
// midright -> new right
// midleft -> new midright
// and calculate the new midleft position and evaluate it for the next iteration
const double newmidleft = left + (midright - left) * sps_left_mid_;
return sps_(f, left, newmidleft, midleft, midright, fl, f(newmidleft), fml, fmr, tolerance);
}
else {
// midright is higher, so exclude everything left of midleft.
// midleft -> new left
// midright -> new midleft
// and calculated the new midright position and evaluate for the next iteration
const double newmidright = midleft + (right - midleft) * sps_right_mid_;
return sps_(f, midleft, midright, newmidright, right, fml, fmr, f(newmidright), fr, tolerance);
}
}
double single_peak_search(
const std::function<double(const double &)> &f,
const double &left, const double &right,
const double &tolerance) {
double left,
double right,
double tol_rel,
double tol_abs) {
if (tol_rel < 0) tol_rel = 0;
if (tol_abs < 0) tol_abs = 0;
double midleft = left + sps_left_mid_*(right - left);
double midright = left + sps_right_mid_*(right - left);
double fl = f(left), fml = f(midleft), fmr = f(midright), fr = f(right);
while (
tol_abs < right - left
and
tol_rel < (right - left) / std::max(std::fabs(left), std::fabs(right))
) {
if (fml >= fmr) {
// midleft is the higher point, so we can exclude everything right of midright.
right = midright; fr = fmr;
midright = midleft; fmr = fml;
midleft = left + (right - left) * sps_left_mid_;
fml = f(midleft);
// If the midpoint is closer to the endpoint than is numerically distinguishable, finish:
if (midleft == left) break;
}
else {
// midright is higher, so exclude everything left of midleft.
left = midleft; fl = fml;
midleft = midright; fml = fmr;
midright = left + (right - left) * sps_right_mid_;
fmr = f(midright);
// If the midpoint is closer to the endpoint than is numerically distinguishable, finish:
if (midright == right) break;
}
const double size = right - left;
const double midleft = left + sps_left_mid_*size;
const double midright = left + sps_right_mid_*size;
// Sometimes, we can run into numerical instability that results in midleft > midright,
// particular when the optimum is close to 0. If we encounter that, swap midleft and
// midright.
if (midleft > midright) {
std::swap(midleft, midright);
std::swap(fml, fmr);
}
}
return sps_(f, left, midleft, midright, right, f(left), f(midleft), f(midright), f(right), tolerance);
// Prefer the end-points for ties (the max might legitimate be an end-point), and prefer
// left over right (for no particularly good reason).
return fl >= fml && fl >= fmr && fl >= fr ? left :
fr >= fmr && fr >= fml ? right :
fml >= fmr ? midleft : midright;
}
}
......@@ -48,3 +48,16 @@ TEST(Maximize, PosQuadraticEnd) {
EXPECT_EQ(-11, single_peak_search(f, -11, -6));
EXPECT_EQ(-2.875, single_peak_search(f, -11, -2.875));
}
TEST(Maximize, Zero) {
auto f = [] (const double &x) -> double { return -(x*x); };
EXPECT_EQ(0, single_peak_search(f, 0, 100));
EXPECT_NEAR(0, single_peak_search(f, -10, 100, 1e-10, 0), 1e-150);
EXPECT_NEAR(0, single_peak_search(f, -10, 100, 0, 1e-10), 1e-10);
EXPECT_NEAR(0, single_peak_search(f, -10, 100, 0, 1e-20), 1e-20);
EXPECT_NEAR(0, single_peak_search(f, -10, 100, 0, 1e-100), 1e-100);
EXPECT_NEAR(0, single_peak_search(f, -10, 100, 0, 0), 1e-150);
EXPECT_EQ(0, single_peak_search(f, -1, 0));
EXPECT_EQ(0, single_peak_search(f, -1e30, 0));
}
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment