rustyqlib/core/optimization/
mod.rs

1//! Multi-dimensional optimization for model calibration, one algorithm
2//! per file. This is the fitting layer for every parametric model —
3//! Heston today ([`equity::heston::calibrate`](crate::equity::heston)),
4//! SABR / Nelson-Siegel or any other least-squares fit tomorrow — so the
5//! machinery lives in one place, like [`solvers`](crate::core::solvers)
6//! does for 1-D root finding.
7//!
8//! **Gradient-based** (analytic gradient optional — central finite
9//! differences fill in when absent):
10//!
11//! - [`steepest_descent`]: robust baseline, linear convergence;
12//! - [`conjugate_gradient`]: nonlinear CG (Polak-Ribiere with restarts);
13//! - [`bfgs`]: quasi-Newton with the inverse-Hessian update — the
14//!   default choice for smooth unconstrained problems;
15//! - [`levenberg_marquardt`]: for least squares `min sum r_i(x)^2`
16//!   specifically — the standard calibration workhorse.
17//!
18//! **Gradient-free**:
19//!
20//! - [`nelder_mead`]: the simplex method, for noisy or non-smooth
21//!   objectives;
22//! - [`differential_evolution`]: seeded, bounded global search for
23//!   multimodal landscapes (e.g. a cold-start calibration before a
24//!   gradient polish).
25//!
26//! Every algorithm can be called directly, or through the pluggable
27//! [`minimize`] with a [`Method`] enum and a [`Problem`] description:
28//!
29//! ```
30//! use rustyqlib::core::optimization::{minimize, Method, OptimConfig, Problem};
31//! let rosenbrock = |x: &[f64]| {
32//!     (1.0 - x[0]).powi(2) + 100.0 * (x[1] - x[0] * x[0]).powi(2)
33//! };
34//! let problem = Problem::scalar(&rosenbrock, vec![-1.2, 1.0]);
35//! let fit = minimize(&OptimConfig::default(), Method::Bfgs, &problem).unwrap();
36//! assert!((fit.x[0] - 1.0).abs() < 1e-5 && (fit.x[1] - 1.0).abs() < 1e-5);
37//! ```
38use crate::core::errors::RustyQLibError;
39
40pub mod bfgs;
41pub mod conjugate_gradient;
42pub mod differential_evolution;
43pub mod levenberg_marquardt;
44pub(crate) mod line_search;
45pub mod nelder_mead;
46pub(crate) mod numerics;
47pub mod steepest_descent;
48
49pub use bfgs::bfgs;
50pub use conjugate_gradient::conjugate_gradient;
51pub use differential_evolution::differential_evolution;
52pub use levenberg_marquardt::levenberg_marquardt;
53pub use nelder_mead::nelder_mead;
54pub use steepest_descent::steepest_descent;
55
56/// Result of an optimization run.
57#[derive(Debug, Clone)]
58pub struct OptimResult {
59    /// The best parameter vector found.
60    pub x: Vec<f64>,
61    /// Objective value at `x` (for least squares: the sum of squared
62    /// residuals).
63    pub value: f64,
64    /// Iterations (gradient methods), simplex steps, or generations.
65    pub iterations: usize,
66    /// True when the method's convergence criterion was met.
67    pub converged: bool,
68}
69
70/// Optimizer configuration: convergence tolerance and iteration cap.
71///
72/// The tolerance applies to each method's natural criterion: the
73/// gradient infinity norm (gradient-based), the simplex value spread
74/// (Nelder-Mead), the population value spread (differential evolution),
75/// or the cost decrease (Levenberg-Marquardt).
76#[derive(Debug, Clone, Copy)]
77pub struct OptimConfig {
78    pub tol: f64,
79    pub max_iter: usize,
80}
81
82impl Default for OptimConfig {
83    fn default() -> Self {
84        Self { tol: 1e-8, max_iter: 500 }
85    }
86}
87
88impl OptimConfig {
89    pub fn new(tol: f64, max_iter: usize) -> Self {
90        Self { tol, max_iter }
91    }
92}
93
94/// The pluggable algorithm choice for [`minimize`].
95#[derive(Debug, Clone, Copy, PartialEq, Eq)]
96pub enum Method {
97    SteepestDescent,
98    ConjugateGradient,
99    Bfgs,
100    LevenbergMarquardt,
101    NelderMead,
102    DifferentialEvolution,
103}
104
105/// An optimization problem: a scalar objective or a residual vector
106/// (least squares), whatever derivatives are available, a start, and
107/// optional bounds.
108pub struct Problem<'a> {
109    pub f: Option<&'a dyn Fn(&[f64]) -> f64>,
110    pub residuals: Option<&'a dyn Fn(&[f64]) -> Vec<f64>>,
111    pub gradient: Option<&'a dyn Fn(&[f64]) -> Vec<f64>>,
112    pub jacobian: Option<&'a dyn Fn(&[f64]) -> Vec<Vec<f64>>>,
113    pub x0: Vec<f64>,
114    /// Per-parameter `(lo, hi)` box, required by differential evolution.
115    pub bounds: Option<Vec<(f64, f64)>>,
116    /// RNG seed for stochastic methods (differential evolution).
117    pub seed: u64,
118}
119
120impl<'a> Problem<'a> {
121    /// A scalar-objective problem.
122    pub fn scalar(f: &'a dyn Fn(&[f64]) -> f64, x0: Vec<f64>) -> Self {
123        Self { f: Some(f), residuals: None, gradient: None, jacobian: None, x0, bounds: None, seed: 42 }
124    }
125
126    /// A least-squares problem `min sum r_i(x)^2`.
127    pub fn least_squares(residuals: &'a dyn Fn(&[f64]) -> Vec<f64>, x0: Vec<f64>) -> Self {
128        Self { f: None, residuals: Some(residuals), gradient: None, jacobian: None, x0, bounds: None, seed: 42 }
129    }
130
131    pub fn with_gradient(mut self, gradient: &'a dyn Fn(&[f64]) -> Vec<f64>) -> Self {
132        self.gradient = Some(gradient);
133        self
134    }
135
136    pub fn with_jacobian(mut self, jacobian: &'a dyn Fn(&[f64]) -> Vec<Vec<f64>>) -> Self {
137        self.jacobian = Some(jacobian);
138        self
139    }
140
141    pub fn with_bounds(mut self, bounds: Vec<(f64, f64)>) -> Self {
142        self.bounds = Some(bounds);
143        self
144    }
145
146    pub fn with_seed(mut self, seed: u64) -> Self {
147        self.seed = seed;
148        self
149    }
150}
151
152/// Minimize `problem` with the chosen [`Method`].
153///
154/// Scalar methods accept either problem form (a residual problem is
155/// minimized as its sum of squares); `LevenbergMarquardt` requires
156/// residuals and `DifferentialEvolution` requires bounds.
157pub fn minimize(
158    cfg: &OptimConfig,
159    method: Method,
160    problem: &Problem,
161) -> Result<OptimResult, RustyQLibError> {
162    // the scalar view of the problem, however it was posed
163    let sum_sq;
164    let f: &dyn Fn(&[f64]) -> f64 = match (problem.f, problem.residuals) {
165        (Some(f), _) => f,
166        (None, Some(r)) => {
167            sum_sq = move |x: &[f64]| r(x).iter().map(|e| e * e).sum::<f64>();
168            &sum_sq
169        }
170        (None, None) => return Err(RustyQLibError::invalid_input("optimization problem", "Problem has neither a scalar objective nor residuals")),
171    };
172    match method {
173        Method::SteepestDescent => Ok(steepest_descent(cfg, f, problem.gradient, &problem.x0)),
174        Method::ConjugateGradient => Ok(conjugate_gradient(cfg, f, problem.gradient, &problem.x0)),
175        Method::Bfgs => Ok(bfgs(cfg, f, problem.gradient, &problem.x0)),
176        Method::LevenbergMarquardt => {
177            let r = problem
178                .residuals
179                .ok_or(RustyQLibError::invalid_input("optimization problem", "Levenberg-Marquardt needs residuals: use Problem::least_squares"))?;
180            Ok(levenberg_marquardt(cfg, r, problem.jacobian, &problem.x0))
181        }
182        Method::NelderMead => Ok(nelder_mead(cfg, f, &problem.x0)),
183        Method::DifferentialEvolution => {
184            let bounds = problem
185                .bounds
186                .as_deref()
187                .ok_or(RustyQLibError::invalid_input("optimization problem", "differential evolution needs bounds: use Problem::with_bounds"))?;
188            Ok(differential_evolution(cfg, f, bounds, problem.seed))
189        }
190    }
191}
192
193#[cfg(test)]
194mod tests {
195    use super::*;
196
197    fn sphere(x: &[f64]) -> f64 {
198        (x[0] - 1.0).powi(2) + (x[1] + 2.0).powi(2)
199    }
200
201    #[test]
202    fn every_scalar_method_is_pluggable_on_one_problem() {
203        let f = |x: &[f64]| sphere(x);
204        let cfg = OptimConfig::new(1e-10, 2000);
205        for method in [
206            Method::SteepestDescent,
207            Method::ConjugateGradient,
208            Method::Bfgs,
209            Method::NelderMead,
210            Method::DifferentialEvolution,
211        ] {
212            let problem = Problem::scalar(&f, vec![4.0, 4.0])
213                .with_bounds(vec![(-10.0, 10.0), (-10.0, 10.0)]);
214            let r = minimize(&cfg, method, &problem).unwrap();
215            assert!(
216                (r.x[0] - 1.0).abs() < 1e-3 && (r.x[1] + 2.0).abs() < 1e-3,
217                "{method:?}: {:?}",
218                r.x
219            );
220        }
221    }
222
223    #[test]
224    fn least_squares_problems_feed_scalar_methods_too() {
225        // residuals of the sphere: r = (x0 - 1, x1 + 2)
226        let r = |x: &[f64]| vec![x[0] - 1.0, x[1] + 2.0];
227        let problem = Problem::least_squares(&r, vec![5.0, 5.0]);
228        let fit = minimize(&OptimConfig::default(), Method::Bfgs, &problem).unwrap();
229        assert!(fit.value < 1e-10, "{fit:?}");
230    }
231
232    #[test]
233    fn missing_requirements_error_clearly() {
234        let f = |x: &[f64]| sphere(x);
235        let no_bounds = Problem::scalar(&f, vec![0.0, 0.0]);
236        assert!(minimize(&OptimConfig::default(), Method::DifferentialEvolution, &no_bounds).is_err());
237        assert!(minimize(&OptimConfig::default(), Method::LevenbergMarquardt, &no_bounds).is_err());
238    }
239}
rustyqlib/core/optimization/mod.rs

rustyqlib/core/optimization/
mod.rs
