First version of the optimizer

2019-10-29 08:59:06 -07:00
parent f4aec2b35e
commit 4488cc110e
10 changed files with 521 additions and 5 deletions
@@ -0,0 +1,121 @@
+use std::default::Default;
+use crate::math::EPSILON;
+use crate::linalg::Vector;
+use crate::optimization::{F, DF};
+use crate::optimization::line_search::LineSearchMethod;
+
+pub trait FirstOrderOptimizer {
+    fn optimize<'a, X: Vector, LS: LineSearchMethod>(&self, f: &'a F<X>, df: &'a DF<X>, x0: &X, ls: &'a LS) -> OptimizerResult<X>;
+}
+
+#[derive(Debug, Clone)]
+pub struct OptimizerResult<X> 
+where X: Vector
+{
+    pub x: X,
+    pub f_x: f64
+}
+
+pub struct GradientDescent {
+    pub max_iter: usize,
+    pub g_rtol: f64,
+    pub g_atol: f64
+}
+
+impl Default for GradientDescent {
+    fn default() -> Self { 
+        GradientDescent {
+            max_iter: 10000,
+            g_rtol: EPSILON.sqrt(),
+            g_atol: EPSILON
+        }
+     }
+}
+
+impl FirstOrderOptimizer for GradientDescent
+{
+
+    fn optimize<'a, X: Vector, LS: LineSearchMethod>(&self, f: &'a F<X>, df: &'a DF<X>, x0: &X, ls: &'a LS) -> OptimizerResult<X> {        
+
+        let mut x = x0.clone();     
+        let mut fx = f(&x);
+
+        let mut gvec = x0.clone();   
+        let mut gnorm = gvec.norm2();        
+
+        let gtol = (gvec.norm2() * self.g_rtol).max(self.g_atol);        
+
+        let mut iter = 0;
+        let mut alpha = 1.0;        
+        df(&mut gvec, &x);         
+
+        while iter < self.max_iter && gnorm > gtol {
+            iter += 1;
+                        
+            let mut step = gvec.negative();
+
+            let f_alpha = |alpha: f64| -> f64 {
+                let mut dx = step.clone();
+                dx.mul_scalar_mut(alpha);
+                f(&dx.add_mut(&x)) // f(x) = f(x .+ gvec .* alpha)
+            };
+
+            let df_alpha = |alpha: f64| -> f64 {                
+                let mut dx = step.clone();
+                let mut dg = gvec.clone();
+                dx.mul_scalar_mut(alpha);
+                df(&mut dg, &dx.add_mut(&x)); //df(x) = df(x .+ gvec .* alpha)
+                gvec.dot(&dg)
+            };
+
+            let df0 = step.dot(&gvec);
+
+            let ls_r = ls.search(&f_alpha, &df_alpha, alpha, fx, df0);
+            alpha = ls_r.alpha;
+            fx = ls_r.f_x;
+            x.add_mut(&step.mul_scalar_mut(alpha));         
+            df(&mut gvec, &x);            
+            gnorm = gvec.norm2();            
+        }  
+
+        let f_x = f(&x);      
+
+        OptimizerResult{
+            x: x,
+            f_x: f_x
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {    
+    use super::*; 
+    use crate::linalg::naive::dense_vector::DenseVector;
+    use crate::optimization::line_search::Backtracking;
+    use crate::optimization::FunctionOrder;
+
+    #[test]
+    fn gradient_descent() { 
+
+            let x0 = DenseVector::from_array(&[-1., 1.]);
+            let f = |x: &DenseVector| {                
+                (1.0 - x.get(0)).powf(2.) + 100.0 * (x.get(1) - x.get(0).powf(2.)).powf(2.)
+            };
+
+            let df = |g: &mut DenseVector, x: &DenseVector| {                                         
+                g.set(0, -2. * (1. - x.get(0)) - 400. * (x.get(1) - x.get(0).powf(2.)) * x.get(0));
+                g.set(1, 200. * (x.get(1) - x.get(0).powf(2.)));                
+            };
+
+            let mut ls: Backtracking = Default::default();
+            ls.order = FunctionOrder::THIRD;
+            let optimizer: GradientDescent = Default::default();
+            
+            let result = optimizer.optimize(&f, &df, &x0, &ls);
+            
+            assert!((result.f_x - 0.0).abs() < EPSILON);
+            assert!((result.x.get(0) - 1.0).abs() < EPSILON);
+            assert!((result.x.get(1) - 1.0).abs() < EPSILON);
+
+    }
+}
@@ -0,0 +1,88 @@
+use crate::math::EPSILON;
+use crate::optimization::FunctionOrder;
+
+pub trait LineSearchMethod {
+    fn search<'a>(&self, f: &(dyn Fn(f64) -> f64), df: &(dyn Fn(f64) -> f64), alpha: f64, f0: f64, df0: f64) -> LineSearchResult;
+}
+
+#[derive(Debug, Clone)]
+pub struct LineSearchResult {
+    pub alpha: f64,
+    pub f_x: f64
+}
+
+pub struct Backtracking {
+    pub c1: f64,
+    pub max_iterations: usize,
+    pub phi: f64,
+    pub plo: f64,
+    pub order: FunctionOrder
+}
+
+impl Default for Backtracking {
+    fn default() -> Self { 
+        Backtracking {
+            c1: 1e-4,
+            max_iterations: 1000,
+            phi: 0.5,
+            plo: 0.1,
+            order: FunctionOrder::SECOND
+        }
+     }
+}
+
+impl LineSearchMethod for Backtracking {
+    
+    fn search<'a>(&self, f: &(dyn Fn(f64) -> f64), _: &(dyn Fn(f64) -> f64), alpha: f64, f0: f64, df0: f64) -> LineSearchResult {        
+
+        let (mut a1, mut a2) = (alpha, alpha);
+        let (mut fx0, mut fx1) = (f0, f(a1));        
+
+        let mut iteration = 0;        
+
+        while fx1 > f0 + self.c1 * a2 * df0 {
+            if iteration > self.max_iterations {
+                panic!("Linesearch failed to converge, reached maximum iterations.");
+            }
+
+            let a_tmp;
+
+            match self.order {
+
+                FunctionOrder::FIRST | FunctionOrder::SECOND => {
+                    a_tmp = - (df0 * a2.powf(2.)) / (2. * (fx1 - f0 - df0*a2))
+                },
+
+                FunctionOrder::THIRD => {
+
+                    let div = 1. / (a1.powf(2.) * a2.powf(2.) * (a2 - a1));
+                    let a = (a1.powf(2.) * (fx1 - f0 - df0*a2) - a2.powf(2.)*(fx0 - f0 - df0*a1))*div;
+                    let b = (-a1.powf(3.) * (fx1 - f0 - df0*a2) + a2.powf(3.)*(fx0 - f0 - df0*a1))*div;
+
+                    
+
+                    if (a - 0.).powf(2.).sqrt() <= EPSILON {
+                        a_tmp = df0 / (2. * b);
+                    } else {
+                        let d = f64::max(b.powf(2.) - 3. * a * df0, 0.);
+                        a_tmp = (-b + d.sqrt()) / (3.*a); //root of quadratic equation
+                    }
+                }
+            }
+
+            a1 = a2;            
+            a2 = f64::max(f64::min(a_tmp, a2*self.phi), a2*self.plo);
+
+            fx0 = fx1;
+            fx1 = f(a2);            
+
+            iteration += 1;
+        }
+
+        LineSearchResult {
+            alpha: a2,
+            f_x: fx1
+        }
+
+    }
+}
@@ -0,0 +1,14 @@
+pub mod first_order;
+pub mod line_search;
+
+use crate::linalg::Vector;
+
+type F<X: Vector> = dyn Fn(&X) -> f64;
+type DF<X: Vector> = dyn Fn(&mut X, &X);
+
+#[derive(Debug)]
+pub enum FunctionOrder {
+    FIRST,
+    SECOND,
+    THIRD
+}