conspire/math/optimize/
mod.rs

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

conspire/math/optimize/
mod.rs
