Skip to content
This repository was archived by the owner on Dec 27, 2019. It is now read-only.
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
76 changes: 42 additions & 34 deletions src/scalar/golden_ratio.rs
Original file line number Diff line number Diff line change
@@ -1,19 +1,19 @@
//! Golden ratio is to minimization what bisection is to root finding.
//! GoldenRatio searches within an interval for a local minimum. At
//! every iteration, the interval is decreased in size by a constant
//! GoldenRatio searches within an interval for a local minimum. At
//! every iteration, the interval is decreased in size by a constant
//! factor, until the desired precision is obtained.
//!
//! This algorithm is guaranteed to converge on a local minimum in a
//!
//! This algorithm is guaranteed to converge on a local minimum in a
//! finite amount of steps under very light smoothess criteria for the
//! target function.
//!
//!
//! In an iteration, the target function is calculated at 4 points:
//! +---------+----+---------+
//! iter 1 a b c d
//! The interval for the next iteration is chosen to be [a,c] if b<c,
//! and [b,d] if c<=b. The distances b-a, c-b, and d-c are chosen in such
//! a way that 3 out of 4 points can be reused, and only 1 new function
//! evaluation is required in the next iteration. If b<c this looks like:
//! a way that 3 out of 4 points can be reused, and only 1 new function
//! evaluation is required in the next iteration. If b<c this looks like:
//! +---------+----+---------+
//! iter 1 a b c d
//! +----+----+----+
Expand All @@ -23,13 +23,13 @@ use std::f64;

#[derive(Builder, Debug)]
pub struct GoldenRatio {
/// The width of the interval at which convergence is satisfactory.
/// The width of the interval at which convergence is satisfactory.
/// Smaller is more precise.
#[builder(default = "1e-8")]
pub xtol: f64,

/// The maximum number of iterations before the search terminates.
/// The number of function evaluations in a bracket search is
/// The number of function evaluations in a bracket search is
/// 2 + <number of iterations>.
/// Bigger is more precise.
#[builder(default = "1000")]
Expand All @@ -42,26 +42,30 @@ pub struct GoldenRatio {
const RATIO: f64 = 2.618033988749895; // 1.5 + 0.5*5f64.sqrt();

impl GoldenRatio {

/// Search for the minimum of `func` around `x0`.
/// The search region is expanded in both directions until an expansion
/// leads to an increase in the function value at the bound of the region.
///
/// Currently `minimize` makes, in some specific cases, at most around 2000
///
/// Currently `minimize` makes, in some specific cases, at most around 2000
/// additional calls to `func` to find a bracketing interval. For more details
/// see https://github.com/to266/optimize/issues/11
pub fn minimize<F>(&self, func: F, x0: f64) -> f64
where F: Fn(f64) -> f64 {
let left = x0 + self.explore(&func, x0, true);
let right = x0 + self.explore(&func, x0, false);
pub fn minimize<F>(&self, mut func: F, x0: f64) -> f64
where
F: FnMut(f64) -> f64,
{
let left = x0 + self.explore(&mut func, x0, true);
let right = x0 + self.explore(&mut func, x0, false);
self.minimize_bracket(func, left, right)
}

/// increase step until func stops decreasing, then return step
fn explore<F>(&self, func: F, x0: f64, explore_left: bool) -> f64
where F: Fn(f64) -> f64 {
let mut step = if explore_left {-1.} else {1.} *
f64::powi(2., f64::log2(f64::EPSILON + x0.abs()) as i32) * f64::EPSILON;
fn explore<F>(&self, mut func: F, x0: f64, explore_left: bool) -> f64
where
F: FnMut(f64) -> f64,
{
let mut step = if explore_left { -1. } else { 1. }
* f64::powi(2., f64::log2(f64::EPSILON + x0.abs()) as i32)
* f64::EPSILON;

let mut fprev = func(x0);
let mut fstepped = func(x0 + step);
Expand All @@ -70,13 +74,15 @@ impl GoldenRatio {
step *= 2.0;
fprev = fstepped;
fstepped = func(x0 + step);
}
}
step
}

/// Search for the minimum of `func` between `left` and `right`.
pub fn minimize_bracket<F>(&self, func: F, left: f64, right: f64) -> f64
where F: Fn(f64) -> f64 {
pub fn minimize_bracket<F>(&self, mut func: F, left: f64, right: f64) -> f64
where
F: FnMut(f64) -> f64,
{
let mut min = left.max(f64::MIN);
let mut max = right.min(f64::MAX);
let mut iter = 0;
Expand All @@ -86,7 +92,7 @@ impl GoldenRatio {
let mut f_b = func(x_b);
let mut f_c = func(x_c);

while (max-min).abs() > self.xtol && iter < self.max_iter {
while (max - min).abs() > self.xtol && iter < self.max_iter {
iter += 1;
if f_b < f_c {
max = x_c;
Expand All @@ -108,7 +114,6 @@ impl GoldenRatio {
}
}


#[cfg(test)]
mod tests {
use super::*;
Expand All @@ -118,37 +123,40 @@ mod tests {
let minimizer = GoldenRatioBuilder::default()
.xtol(1e-7)
.max_iter(1000)
.build().unwrap();
let f = |x: f64| (x-0.2).powi(2);
.build()
.unwrap();
let f = |x: f64| (x - 0.2).powi(2);
let res = minimizer.minimize_bracket(&f, -1.0, 1.0);

println!("res: {}", res);
assert!( (res - 0.2).abs() <= 1e-7);
assert!((res - 0.2).abs() <= 1e-7);
}

#[test]
fn golden_ratio_x0() {
let minimizer = GoldenRatioBuilder::default()
.xtol(1e-7)
.max_iter(1000)
.build().unwrap();
let f = |x: f64| (x-0.2).powi(2);
.build()
.unwrap();
let f = |x: f64| (x - 0.2).powi(2);
let res = minimizer.minimize(&f, 10.0);

println!("res: {}", res);
assert!( (res - 0.2).abs() <= 1e-7);
assert!((res - 0.2).abs() <= 1e-7);
}

#[test]
fn monotone_x0() {
let minimizer = GoldenRatioBuilder::default()
.xtol(1e-7)
.max_iter(1000)
.build().unwrap();
.build()
.unwrap();
let f = |x: f64| x;
let res = minimizer.minimize(&f, 10.0);

println!("res: {}", res);
assert!( res.is_infinite() );
assert!(res.is_infinite());
}
}
}
Loading