Implement different rules for step size in derivatives approximation

src/algo/trust_region.rs

@@ -37,7 +37,7 @@ use thiserror::Error;
 
 use crate::{
     core::{Domain, Function, Optimizer, Problem, RealField as _, Solver, System},
-    derivatives::{Gradient, Hessian, Jacobian},
+    derivatives::{Gradient, Hessian, Jacobian, StepRule},
 };
 
 /// Specification for initial value of trust region size.
@@ -77,6 +77,14 @@ pub struct TrustRegionOptions<P: Problem> {
     /// Determines whether steps that increase the error can be accepted.
     /// Default: `true`.
     allow_ascent: bool,
+    /// Relative epsilon used in Jacobian and gradient computations.
+    /// Default: `sqrt(EPSILON)`.
+    eps_sqrt: P::Field,
+    /// Relative epsilon used in Hessian computations.
+    /// Default: `cbrt(EPSILON)`.
+    eps_cbrt: P::Field,
+    /// Step rule used in Jacobian, gradient and Hessian computations.
+    step_rule: StepRule,
     #[getset(skip)]
     prefer_greater_magnitude_in_cauchy: bool,
 }
@@ -111,6 +119,9 @@ impl<P: Problem> Default for TrustRegionOptions<P> {
             rejections_thresh: 10,
             allow_ascent: true,
             prefer_greater_magnitude_in_cauchy: false,
+            eps_sqrt: P::Field::EPSILON_SQRT,
+            eps_cbrt: P::Field::EPSILON_CBRT,
+            step_rule: StepRule::default(),
         }
     }
 }
@@ -225,6 +236,8 @@ impl<R: System> Solver<R> for TrustRegion<R> {
             rejections_thresh,
             allow_ascent,
             prefer_greater_magnitude_in_cauchy,
+            eps_sqrt,
+            step_rule,
             ..
         } = self.options;
 
@@ -264,7 +277,7 @@ impl<R: System> Solver<R> for TrustRegion<R> {
 
         // Compute r(x) and r'(x).
         r.eval(x, rx);
-        jac.compute(r, x, scale, rx);
+        jac.compute(r, x, scale, rx, eps_sqrt, step_rule);
 
         let rx_norm = rx.norm();
 
@@ -706,6 +719,9 @@ impl<F: Function> Optimizer<F> for TrustRegion<F> {
             accept_thresh,
             rejections_thresh,
             allow_ascent,
+            eps_sqrt,
+            eps_cbrt,
+            step_rule,
             ..
         } = self.options;
 
@@ -742,8 +758,8 @@ impl<F: Function> Optimizer<F> for TrustRegion<F> {
 
         // Compute f(x), grad f(x) and H(x).
         let mut fx = f.apply(x);
-        grad.compute(f, x, scale, fx);
-        hes.compute(f, x, scale, fx);
+        grad.compute(f, x, scale, fx, eps_sqrt, step_rule);
+        hes.compute(f, x, scale, fx, eps_cbrt, step_rule);
 
         let estimate_delta = *delta == zero;
         if estimate_delta {

src/analysis.rs

@@ -1,13 +1,12 @@
 //! Various supporting analyses.
 
 use nalgebra::{
-    convert, storage::StorageMut, ComplexField, DimName, Dyn, IsContiguous, OVector, RealField,
-    Vector, U1,
+    convert, storage::StorageMut, ComplexField, DimName, Dyn, IsContiguous, OVector, Vector, U1,
 };
 
 use crate::{
-    core::{Domain, RealField as _, System},
-    derivatives::Jacobian,
+    core::{Domain, RealField, System},
+    derivatives::{Jacobian, StepRule},
 };
 
 /// Estimates magnitude of the variable given lower and upper bounds.
@@ -55,7 +54,14 @@ where
 
     // Compute r'(x) in the initial point.
     r.eval(x, rx);
-    let jac1 = Jacobian::new(r, x, &scale, rx);
+    let jac1 = Jacobian::new(
+        r,
+        x,
+        &scale,
+        rx,
+        R::Field::EPSILON_SQRT,
+        StepRule::default(),
+    );
 
     // Compute Newton step.
     let mut p = rx.clone_owned();
@@ -70,7 +76,14 @@ where
 
     // Compute r'(x) after one Newton step.
     r.eval(x, rx);
-    let jac2 = Jacobian::new(r, x, &scale, rx);
+    let jac2 = Jacobian::new(
+        r,
+        x,
+        &scale,
+        rx,
+        R::Field::EPSILON_SQRT,
+        StepRule::default(),
+    );
 
     // Linear variables have no effect on the Jacobian matrix. They can be
     // recognized by observing no change in corresponding columns (i.e.,

src/core/domain.rs

@@ -5,10 +5,10 @@ use std::iter::FromIterator;
 use fastrand::Rng;
 use na::{Dim, DimName};
 use nalgebra as na;
-use nalgebra::{storage::StorageMut, OVector, RealField, Vector};
+use nalgebra::{storage::StorageMut, OVector, Vector};
 
 use crate::analysis::estimate_magnitude_from_bounds;
-use crate::core::Sample;
+use crate::core::{RealField, Sample};
 
 /// Domain for a problem.
 #[derive(Clone)]