Merge pull request #28 from denehoffman/formatting_and_fixes

Formatting and fixes
denehoffman · Sep 9, 2024 · 2b3cbc7 · 2b3cbc7
2 parents 7bb0f01 + 833ae74
commit 2b3cbc7
Show file tree

Hide file tree

Showing 8 changed files with 289 additions and 65 deletions.
diff --git a/Cargo.lock b/Cargo.lock
diff --git a/Cargo.toml b/Cargo.toml
@@ -18,6 +18,7 @@ num = "0.4.3"
 
 [dev-dependencies]
 criterion = "0.5.1"
+float-cmp = "0.9.0"
 
 [[bench]]
 name = "rosenbrock_benchmark"

diff --git a/README.md b/README.md
@@ -65,21 +65,31 @@ To minimize this function, we could consider using the Nelder-Mead algorithm:
 use ganesh::prelude::*;
 use ganesh::algorithms::NelderMead;
 
-let problem = Rosenbrock { n: 2 };
-let nm = NelderMead::default();
-let mut m = Minimizer::new(nm, 2);
-let x0 = &[2.0, 2.0];
-m.minimize(&problem, x0, &mut ()).unwrap();
-println!("{}", m.status);
+fn main() -> Result<(), Infallible> {
+    let problem = Rosenbrock { n: 2 };
+    let nm = NelderMead::default();
+    let mut m = Minimizer::new(nm, 2);
+    let x0 = &[2.0, 2.0];
+    m.minimize(&problem, x0, &mut ())?;
+    println!("{}", m.status);
+    Ok(())
+}
 ```
 
 This should output
 ```shell
 MSG:       term_f = STDDEV
-X:         [0.9999999946231828, 0.9999999884539057]
-F(X):      0.00000000000000009170942877687133
-N_EVALS:   160
+X:         +1.000 ± 0.707
+           +1.000 ± 1.416
+F(X):      +0.000
+N_F_EVALS: 159
+N_G_EVALS: 0
 CONVERGED: true
+COV:       
+  ┌             ┐
+  │ 0.500 1.000 │
+  │ 1.000 2.005 │
+  └             ┘
 ```
 
 # Bounds

diff --git a/src/algorithms/bfgs.rs b/src/algorithms/bfgs.rs
@@ -1,7 +1,7 @@
 use nalgebra::{DMatrix, DVector, RealField, Scalar};
 use num::Float;
 
-use crate::{Algorithm, Bound, Function, Status};
+use crate::{convert, Algorithm, Bound, Function, Status};
 
 use super::line_search::{LineSearch, StrongWolfeLineSearch};
 
@@ -68,12 +68,12 @@ pub struct BFGS<T: Scalar, U, E> {
     status: Status<T>,
     x: DVector<T>,
     g: DVector<T>,
-    p: DVector<T>,
     h_inv: DMatrix<T>,
     f_previous: T,
     terminator_f: BFGSFTerminator<T>,
     terminator_g: BFGSGTerminator<T>,
     line_search: Box<dyn LineSearch<T, U, E>>,
+    max_step: T,
     error_mode: BFGSErrorMode,
 }
 
@@ -117,7 +117,6 @@ where
             status: Default::default(),
             x: Default::default(),
             g: Default::default(),
-            p: Default::default(),
             h_inv: Default::default(),
             f_previous: T::infinity(),
             terminator_f: BFGSFTerminator {
@@ -127,6 +126,7 @@ where
                 tol_g_abs: Float::cbrt(T::epsilon()),
             },
             line_search: Box::new(StrongWolfeLineSearch::default()),
+            max_step: convert!(1e8, T),
             error_mode: Default::default(),
         }
     }
@@ -150,7 +150,7 @@ where
 
 impl<T, U, E> Algorithm<T, U, E> for BFGS<T, U, E>
 where
-    T: RealField + Float,
+    T: RealField + Float + Default,
 {
     fn initialize(
         &mut self,
@@ -159,6 +159,8 @@ where
         bounds: Option<&Vec<Bound<T>>>,
         user_data: &mut U,
     ) -> Result<(), E> {
+        self.status = Status::default();
+        self.f_previous = T::infinity();
         self.h_inv = DMatrix::identity(x0.len(), x0.len());
         self.x = Bound::to_unbounded(x0, bounds);
         self.g = func.gradient_bounded(self.x.as_slice(), bounds, user_data)?;
@@ -178,23 +180,23 @@ where
         bounds: Option<&Vec<Bound<T>>>,
         user_data: &mut U,
     ) -> Result<(), E> {
-        self.p = -&self.h_inv * &self.g;
+        let d = -&self.h_inv * &self.g;
         let (valid, alpha, f_kp1, g_kp1) = self.line_search.search(
             &self.x,
-            &self.p,
-            None,
+            &d,
+            Some(self.max_step),
             func,
             bounds,
             user_data,
             &mut self.status,
         )?;
         if valid {
-            let s = self.p.scale(alpha);
+            let dx = d.scale(alpha);
             let grad_kp1_vec = g_kp1;
-            let y = &grad_kp1_vec - &self.g;
+            let dg = &grad_kp1_vec - &self.g;
             let n = self.x.len();
-            self.update_h_inv(i_step, n, &s, &y);
-            self.x += s;
+            self.update_h_inv(i_step, n, &dx, &dg);
+            self.x += dx;
             self.g = grad_kp1_vec;
             self.status
                 .update_position((Bound::to_bounded(self.x.as_slice(), bounds), f_kp1));
@@ -236,12 +238,49 @@ where
                 if covariance.is_none() {
                     covariance = hessian.pseudo_inverse(Float::cbrt(T::epsilon())).ok();
                 }
-                self.status.cov = covariance
+                self.status.set_cov(covariance);
             }
             BFGSErrorMode::ApproximateHessian => {
-                self.status.cov = Some(self.h_inv.clone());
+                self.status.set_cov(Some(self.h_inv.clone()));
             }
         }
         Ok(())
     }
 }
+
+#[cfg(test)]
+mod tests {
+    use std::convert::Infallible;
+
+    use float_cmp::approx_eq;
+
+    use crate::{prelude::*, test_functions::Rosenbrock};
+
+    use super::BFGS;
+
+    #[test]
+    fn test_bfgs() -> Result<(), Infallible> {
+        let algo = BFGS::default();
+        let mut m = Minimizer::new(algo, 2).with_max_steps(10000);
+        let problem = Rosenbrock { n: 2 };
+        m.minimize(&problem, &[-2.0, 2.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[2.0, 2.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[2.0, -2.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[-2.0, -2.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[0.0, 0.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[1.0, 1.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        Ok(())
+    }
+}
diff --git a/src/algorithms/lbfgs.rs b/src/algorithms/lbfgs.rs
@@ -3,7 +3,7 @@ use std::collections::VecDeque;
 use nalgebra::{DVector, RealField, Scalar};
 use num::Float;
 
-use crate::{Algorithm, Bound, Function, Status};
+use crate::{convert, Algorithm, Bound, Function, Status};
 
 use super::line_search::{LineSearch, StrongWolfeLineSearch};
 
@@ -68,14 +68,14 @@ pub struct LBFGS<T: Scalar, U, E> {
     status: Status<T>,
     x: DVector<T>,
     g: DVector<T>,
-    p: DVector<T>,
     f_previous: T,
     terminator_f: LBFGSFTerminator<T>,
     terminator_g: LBFGSGTerminator<T>,
     line_search: Box<dyn LineSearch<T, U, E>>,
     m: usize,
     y_store: VecDeque<DVector<T>>,
     s_store: VecDeque<DVector<T>>,
+    max_step: T,
     error_mode: LBFGSErrorMode,
 }
 
@@ -124,7 +124,6 @@ where
             status: Default::default(),
             x: Default::default(),
             g: Default::default(),
-            p: Default::default(),
             f_previous: T::infinity(),
             terminator_f: LBFGSFTerminator {
                 tol_f_abs: T::epsilon(),
@@ -136,6 +135,7 @@ where
             m: 10,
             y_store: VecDeque::default(),
             s_store: VecDeque::default(),
+            max_step: convert!(1e8, T),
             error_mode: Default::default(),
         }
     }
@@ -148,7 +148,7 @@ where
     fn g_approx(&self) -> DVector<T> {
         let m = self.s_store.len();
         let mut q = self.g.clone();
-        if m <= 2 {
+        if m < 1 {
             return q;
         }
         let mut a = vec![T::zero(); m];
@@ -171,7 +171,7 @@ where
 
 impl<T, U, E> Algorithm<T, U, E> for LBFGS<T, U, E>
 where
-    T: RealField + Float,
+    T: RealField + Float + Default,
 {
     fn initialize(
         &mut self,
@@ -180,6 +180,8 @@ where
         bounds: Option<&Vec<Bound<T>>>,
         user_data: &mut U,
     ) -> Result<(), E> {
+        self.status = Status::default();
+        self.f_previous = T::infinity();
         self.x = Bound::to_unbounded(x0, bounds);
         self.g = func.gradient_bounded(self.x.as_slice(), bounds, user_data)?;
         self.status.inc_n_g_evals();
@@ -198,34 +200,38 @@ where
         bounds: Option<&Vec<Bound<T>>>,
         user_data: &mut U,
     ) -> Result<(), E> {
-        self.p = -self.g_approx();
+        let d = -self.g_approx();
         let (valid, alpha, f_kp1, g_kp1) = self.line_search.search(
             &self.x,
-            &self.p,
-            None,
+            &d,
+            Some(self.max_step),
             func,
             bounds,
             user_data,
             &mut self.status,
         )?;
         if valid {
-            let dx = self.p.scale(alpha);
-            self.s_store.push_back(dx.clone());
-            if self.s_store.len() > self.m {
-                self.s_store.pop_front();
-            }
+            let dx = d.scale(alpha);
             let grad_kp1_vec = g_kp1;
-            self.y_store.push_back(&grad_kp1_vec - &self.g);
-            if self.y_store.len() > self.m {
-                self.y_store.pop_front();
+            let dg = &grad_kp1_vec - &self.g;
+            let sy = dx.dot(&dg);
+            let yy = dg.dot(&dg);
+            if sy > T::epsilon() * yy {
+                self.s_store.push_back(dx.clone());
+                self.y_store.push_back(dg);
+                if self.s_store.len() > self.m {
+                    self.s_store.pop_front();
+                    self.y_store.pop_front();
+                }
             }
             self.x += dx;
             self.g = grad_kp1_vec;
             self.status
                 .update_position((Bound::to_bounded(self.x.as_slice(), bounds), f_kp1));
         } else {
-            self.status.update_message("LINE SEARCH FAILED");
-            self.status.set_converged();
+            // reboot
+            self.s_store.clear();
+            self.y_store.clear();
         }
         Ok(())
     }
@@ -262,9 +268,46 @@ where
                 if covariance.is_none() {
                     covariance = hessian.pseudo_inverse(Float::cbrt(T::epsilon())).ok();
                 }
-                self.status.cov = covariance
+                self.status.set_cov(covariance);
             }
         }
         Ok(())
     }
 }
+
+#[cfg(test)]
+mod tests {
+    use std::convert::Infallible;
+
+    use float_cmp::approx_eq;
+
+    use crate::{prelude::*, test_functions::Rosenbrock};
+
+    use super::LBFGS;
+
+    #[test]
+    fn test_lbfgs() -> Result<(), Infallible> {
+        let algo = LBFGS::default();
+        let mut m = Minimizer::new(algo, 2);
+        let problem = Rosenbrock { n: 2 };
+        m.minimize(&problem, &[-2.0, 2.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[2.0, 2.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[2.0, -2.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[-2.0, -2.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[0.0, 0.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[1.0, 1.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        Ok(())
+    }
+}