conspire/math/optimize/
mod.rs

1#[cfg(test)]
2mod test;
3
4mod constraint;
5mod gradient_descent;
6mod line_search;
7mod newton_raphson;
8mod strategy;
9mod tolerance;
10mod trust_region;
11
12pub use constraint::EqualityConstraint;
13pub use gradient_descent::GradientDescent;
14pub use line_search::{LineSearch, LineSearchError};
15pub use newton_raphson::NewtonRaphson;
16pub use strategy::SolveStrategy;
17pub use tolerance::Tolerances;
18pub use trust_region::TrustRegion;
19
20use crate::{
21    math::{
22        Erase, Jacobian, Quantity, Scalar, Solution, Style, StyledError, Tensor, Vector,
23        assert::AssertionError,
24        matrix::square::SquareMatrixError,
25        sparse::{CscMatrix, SparseError, SparseSolver},
26        styled_error,
27    },
28    units::UnitDiv,
29};
30use std::{
31    fmt::{Debug, Display},
32    ops::Mul,
33};
34
35/// The step size taking a decrement of type `D` to an increment of `X`.
36pub type StepSize<D, X> = Quantity<<<X as Tensor>::Unit as UnitDiv<<D as Tensor>::Unit>>::Output>;
37
38/// Zeroth-order root-finding algorithms.
39pub trait ZerothOrderRootFinding<F, X> {
40    fn root(
41        &self,
42        function: impl FnMut(&X) -> Result<F, String>,
43        initial_guess: X,
44        equality_constraint: EqualityConstraint,
45    ) -> Result<X, OptimizationError>;
46}
47
48/// First-order root-finding algorithms.
49pub trait FirstOrderRootFinding<F, J, X> {
50    fn root(
51        &self,
52        function: impl FnMut(&X) -> Result<F, String>,
53        jacobian: impl FnMut(&X) -> Result<J, String>,
54        initial_guess: X,
55        equality_constraint: EqualityConstraint,
56        sparse: Option<SparseSolver>,
57    ) -> Result<X, OptimizationError>;
58}
59
60/// First-order root-finding algorithms that hand out each increment before
61/// applying it.
62///
63/// The solver keeps the iteration; the increment is only lent to the caller so
64/// that whatever was eliminated from the system can be carried along with it.
65///
66/// The increment is lent whole, with the step it is about to be scaled by
67/// alongside. Elimination solves one direction for the eliminated variables and
68/// the retained ones together, so shortening the step has to shorten both by
69/// the same amount, exactly as it would if nothing had been eliminated. Handing
70/// over the shortened increment instead would invite a fresh solve against it,
71/// which is a different direction rather than less of the same one.
72///
73/// A step is offered before it is taken. The caller is asked to report whether
74/// the state it arrives at is admissible, and only later told to keep it.
75pub trait FirstOrderRootFindingIncremental<F, J, X> {
76    fn root_incremental(
77        &self,
78        function: impl FnMut(&X) -> Result<F, String>,
79        jacobian: impl FnMut(&X) -> Result<J, String>,
80        update: impl FnMut(&X, &Vector, Scalar, bool) -> Result<(), String>,
81        initial_guess: X,
82        equality_constraint: EqualityConstraint,
83        sparse: Option<SparseSolver>,
84    ) -> Result<X, OptimizationError>;
85}
86
87/// First-order optimization algorithms.
88pub trait FirstOrderOptimization<F, J, X> {
89    fn minimize(
90        &self,
91        function: impl FnMut(&X) -> Result<F, String>,
92        jacobian: impl FnMut(&X) -> Result<J, String>,
93        initial_guess: X,
94        equality_constraint: EqualityConstraint,
95    ) -> Result<X, OptimizationError>;
96}
97
98/// Second-order optimization algorithms.
99pub trait SecondOrderOptimization<F, J, H, X> {
100    fn minimize(
101        &self,
102        function: impl FnMut(&X) -> Result<F, String>,
103        jacobian: impl FnMut(&X) -> Result<J, String>,
104        hessian: impl FnMut(&X) -> Result<H, String>,
105        initial_guess: X,
106        equality_constraint: EqualityConstraint,
107        sparse: Option<SparseSolver>,
108    ) -> Result<X, OptimizationError>;
109}
110
111/// Second-order optimization algorithms that hand out each increment before
112/// applying it.
113///
114/// The counterpart of [`FirstOrderRootFindingIncremental`] for problems with an
115/// energy to descend, and the increment is lent on the same terms.
116///
117/// What the line search measures is the energy of the whole state, eliminated
118/// variables included. Each trial is offered through the same update, so the
119/// eliminated variables are already standing where the trial puts them by the
120/// time the energy there is asked for.
121pub trait SecondOrderOptimizationIncremental<F, J, H, X> {
122    #[expect(clippy::too_many_arguments)]
123    fn minimize_incremental(
124        &self,
125        function: impl FnMut(&X) -> Result<F, String>,
126        jacobian: impl FnMut(&X) -> Result<J, String>,
127        hessian: impl FnMut(&X) -> Result<H, String>,
128        update: impl FnMut(&X, &Vector, Scalar, bool) -> Result<(), String>,
129        initial_guess: X,
130        equality_constraint: EqualityConstraint,
131        sparse: Option<SparseSolver>,
132    ) -> Result<X, OptimizationError>;
133}
134
135/// First-order root-finding algorithms for problems split into global and local variables.
136#[expect(clippy::too_many_arguments)]
137pub trait FirstOrderRootFindingBlock<U, V, Ru, Rv, Kuu, Kvu, Kuv, Kvv> {
138    fn root_block(
139        &self,
140        residual_global: impl FnMut(&U, &V) -> Result<Ru, String>,
141        residual_local: impl FnMut(&U, &V) -> Result<Rv, String>,
142        tangents: impl FnMut(&U, &V) -> Result<(Kuu, Kvu, Kuv, Kvv), String>,
143        initial_guess: (U, V),
144        constraint_global: (CscMatrix, Vector),
145        constraint_local: (CscMatrix, Vector),
146        sparse: Option<SparseSolver>,
147        strategy: SolveStrategy,
148    ) -> Result<(U, V), OptimizationError>;
149}
150
151/// Second-order optimization algorithms for problems split into global and local variables.
152#[expect(clippy::too_many_arguments)]
153pub trait SecondOrderOptimizationBlock<F, U, V, Ru, Rv, Kuu, Kvu, Kuv, Kvv> {
154    fn minimize_block(
155        &self,
156        function: impl FnMut(&U, &V) -> Result<F, String>,
157        residual_global: impl FnMut(&U, &V) -> Result<Ru, String>,
158        residual_local: impl FnMut(&U, &V) -> Result<Rv, String>,
159        tangents: impl FnMut(&U, &V) -> Result<(Kuu, Kvu, Kuv, Kvv), String>,
160        initial_guess: (U, V),
161        constraint_global: (CscMatrix, Vector),
162        constraint_local: (CscMatrix, Vector),
163        sparse: Option<SparseSolver>,
164        strategy: SolveStrategy,
165    ) -> Result<(U, V), OptimizationError>;
166}
167
168trait BacktrackingLineSearch<J, X>
169where
170    Self: Debug,
171{
172    fn backtracking_line_search<D, E>(
173        &self,
174        mut function: impl FnMut(&X, Scalar) -> Result<Scalar, String>,
175        mut jacobian: impl FnMut(&X) -> Result<J, String>,
176        argument: &X,
177        jacobian0: &J,
178        decrement: &D,
179        step_size: Scalar,
180    ) -> Result<Scalar, OptimizationError>
181    where
182        J: Erase<Erased = E> + Jacobian,
183        D: Erase<Erased = E> + Tensor,
184        E: Tensor,
185        X: Solution,
186        <X as Tensor>::Unit: UnitDiv<<D as Tensor>::Unit>,
187        for<'a> &'a D: Mul<StepSize<D, X>, Output = X>,
188    {
189        if matches!(self.get_line_search(), LineSearch::None) {
190            Ok(step_size)
191        } else {
192            self.get_line_search()
193                .backtrack(
194                    &mut function,
195                    &mut jacobian,
196                    argument,
197                    jacobian0,
198                    decrement,
199                    step_size,
200                )
201                .map_err(|error| OptimizationError::upstream(error, self))
202        }
203    }
204    fn get_line_search(&self) -> &LineSearch;
205}
206
207/// Possible errors encountered during optimization.
208pub enum OptimizationError {
209    Intermediate(String),
210    MaximumStepsReached(usize, String),
211    NotMinimum(String, String),
212    Upstream(String, String),
213    SingularMatrix,
214    UnsymmetricMatrix,
215}
216
217impl OptimizationError {
218    pub fn upstream(error: impl Display, context: &(impl Debug + ?Sized)) -> Self {
219        Self::Upstream(format!("{error}"), format!("{context:?}"))
220    }
221}
222
223impl From<String> for OptimizationError {
224    fn from(error: String) -> Self {
225        Self::Intermediate(error)
226    }
227}
228
229impl StyledError for OptimizationError {
230    fn message(&self, style: &Style) -> String {
231        let (h, c) = (style.headline, style.frame);
232        match self {
233            Self::Intermediate(message) => message.to_string(),
234            Self::MaximumStepsReached(steps, solver) => format!(
235                "{h}Maximum number of steps ({steps}) reached.{c}\n\
236                In solver: {solver}."
237            ),
238            Self::NotMinimum(solution, solver) => format!(
239                "{h}The obtained solution is not a minimum.{c}\n\
240                For solution: {solution}.\n\
241                In solver: {solver}."
242            ),
243            Self::SingularMatrix => format!("{h}Matrix is singular."),
244            Self::UnsymmetricMatrix => format!("{h}Matrix is not symmetric."),
245            Self::Upstream(error, solver) => format!(
246                "{error}{c}\n\
247                In solver: {solver}."
248            ),
249        }
250    }
251}
252
253styled_error!(OptimizationError);
254
255impl From<OptimizationError> for String {
256    fn from(error: OptimizationError) -> Self {
257        error.to_string()
258    }
259}
260
261impl From<OptimizationError> for AssertionError {
262    fn from(error: OptimizationError) -> Self {
263        Self {
264            message: error.to_string(),
265        }
266    }
267}
268
269impl From<SquareMatrixError> for OptimizationError {
270    fn from(_error: SquareMatrixError) -> Self {
271        Self::SingularMatrix
272    }
273}
274
275impl From<SparseError> for OptimizationError {
276    fn from(error: SparseError) -> Self {
277        match error {
278            SparseError::Singular => Self::SingularMatrix,
279            SparseError::Unsymmetric => Self::UnsymmetricMatrix,
280        }
281    }
282}
conspire/math/optimize/mod.rs

conspire/math/optimize/
mod.rs
