Skip to main content

conspire/math/optimize/line_search/
mod.rs

1#[cfg(test)]
2mod test;
3
4use crate::math::{
5    Erase, Jacobian, Quantity, Scalar, Solution, Style, StyledError, Tensor, styled_error,
6};
7use std::{
8    fmt::{self, Debug, Display, Formatter},
9    ops::Mul,
10};
11
12/// Available line search algorithms.
13#[derive(Clone, Debug)]
14pub enum LineSearch {
15    /// The Armijo condition.
16    Armijo {
17        control: Scalar,
18        cut_back: Scalar,
19        max_steps: usize,
20    },
21    /// Backtrack for errors.
22    Error { cut_back: Scalar, max_steps: usize },
23    /// The Goldstein conditions.
24    Goldstein {
25        control: Scalar,
26        cut_back: Scalar,
27        max_steps: usize,
28    },
29    /// The Wolfe conditions.
30    Wolfe {
31        control_1: Scalar,
32        control_2: Scalar,
33        cut_back: Scalar,
34        max_steps: usize,
35        strong: bool,
36    },
37    /// No line search.
38    None,
39}
40
41impl Default for LineSearch {
42    fn default() -> Self {
43        Self::Armijo {
44            control: 1e-3,
45            cut_back: 9e-1,
46            max_steps: 100,
47        }
48    }
49}
50
51impl Display for LineSearch {
52    fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
53        match self {
54            Self::Armijo { .. } => write!(f, "Armijo"),
55            Self::Error { .. } => write!(f, "Error"),
56            Self::Goldstein { .. } => write!(f, "Goldstein"),
57            Self::Wolfe { .. } => write!(f, "Wolfe"),
58            Self::None { .. } => write!(f, "None"),
59        }
60    }
61}
62
63impl LineSearch {
64    /// Backtrack on a merit function of the step size alone.
65    ///
66    /// The exact penalty function is not differentiable, its norm having a kink
67    /// wherever a constraint is satisfied, so its slope along the step is
68    /// supplied rather than recovered from a gradient.
69    pub fn backtrack_merit(
70        &self,
71        mut merit: impl FnMut(Scalar) -> Result<Scalar, String>,
72        value: Scalar,
73        slope: Scalar,
74        step_size: Scalar,
75    ) -> Result<Scalar, LineSearchError> {
76        if step_size <= 0.0 {
77            return Err(LineSearchError::NegativeStepSize(
78                format!("{self:?}"),
79                step_size,
80            ));
81        } else if slope <= 0.0 {
82            return Err(LineSearchError::NotDescentDirection(format!("{self:?}")));
83        }
84        let mut n = step_size;
85        match self {
86            Self::Armijo {
87                control,
88                cut_back,
89                max_steps,
90            } => {
91                let t = control * slope;
92                for _ in 0..*max_steps {
93                    if let Ok(trial) = merit(n)
94                        && value - trial >= n * t
95                    {
96                        return Ok(n);
97                    } else {
98                        n *= cut_back
99                    }
100                }
101                Err(LineSearchError::MaximumStepsReached(
102                    format!("{self:?}"),
103                    *max_steps,
104                ))
105            }
106            Self::Error {
107                cut_back,
108                max_steps,
109            } => {
110                for _ in 0..*max_steps {
111                    if merit(n).is_ok() {
112                        return Ok(n);
113                    } else {
114                        n *= cut_back
115                    }
116                }
117                Err(LineSearchError::MaximumStepsReached(
118                    format!("{self:?}"),
119                    *max_steps,
120                ))
121            }
122            Self::Goldstein {
123                control,
124                cut_back,
125                max_steps,
126            } => {
127                let t = control * slope;
128                let u = (1.0 - control) * slope;
129                let mut v;
130                for _ in 0..*max_steps {
131                    if let Ok(trial) = merit(n) {
132                        v = value - trial;
133                        if n * u < v || v < n * t {
134                            n *= cut_back
135                        } else {
136                            return Ok(n);
137                        }
138                    } else {
139                        n *= cut_back
140                    }
141                }
142                Err(LineSearchError::MaximumStepsReached(
143                    format!("{self:?}"),
144                    *max_steps,
145                ))
146            }
147            Self::Wolfe { .. } => panic!(
148                "The Wolfe conditions need the gradient of the merit function, which the exact penalty function does not have."
149            ),
150            Self::None => {
151                panic!("Cannot call backtracking line search when there is no algorithm.")
152            }
153        }
154    }
155    pub fn backtrack<X, J, D, W, E>(
156        &self,
157        mut function: impl FnMut(&X, Scalar) -> Result<Scalar, String>,
158        mut jacobian: impl FnMut(&X) -> Result<J, String>,
159        argument: &X,
160        jacobian0: &J,
161        decrement: &D,
162        step_size: Scalar,
163    ) -> Result<Scalar, LineSearchError>
164    where
165        J: Erase<Erased = E> + Jacobian,
166        D: Erase<Erased = E>,
167        E: Tensor,
168        X: Solution,
169        for<'a> &'a D: Mul<Quantity<W>, Output = X>,
170    {
171        if step_size <= 0.0 {
172            return Err(LineSearchError::NegativeStepSize(
173                format!("{self:?}"),
174                step_size,
175            ));
176        }
177        let mut n = step_size;
178        let f = if let Ok(value) = function(argument, 0.0) {
179            value
180        } else {
181            return Err(LineSearchError::InvalidStartingPoint(format!("{self:?}")));
182        };
183        let m = jacobian0.erase().full_contraction(decrement.erase());
184        if m <= 0.0 {
185            return Err(LineSearchError::NotDescentDirection(format!("{self:?}")));
186        }
187        let trial = |n: Scalar| decrement * Quantity::new(-n) + argument;
188        match self {
189            Self::Armijo {
190                control,
191                cut_back,
192                max_steps,
193            } => {
194                let mut f_n;
195                let t = control * m;
196                for _ in 0..*max_steps {
197                    f_n = function(&trial(n), n);
198                    if let Ok(value) = f_n
199                        && f - value >= n * t
200                    {
201                        return Ok(n);
202                    } else {
203                        n *= cut_back
204                    }
205                }
206                Err(LineSearchError::MaximumStepsReached(
207                    format!("{self:?}"),
208                    *max_steps,
209                ))
210            }
211            Self::Error {
212                cut_back,
213                max_steps,
214            } => {
215                for _ in 0..*max_steps {
216                    if function(&trial(n), n).is_ok() {
217                        return Ok(n);
218                    } else {
219                        n *= cut_back
220                    }
221                }
222                Err(LineSearchError::MaximumStepsReached(
223                    format!("{self:?}"),
224                    *max_steps,
225                ))
226            }
227            Self::Goldstein {
228                control,
229                cut_back,
230                max_steps,
231            } => {
232                let mut f_n;
233                let t = control * m;
234                let u = (1.0 - control) * m;
235                let mut v;
236                for _ in 0..*max_steps {
237                    f_n = function(&trial(n), n);
238                    if let Ok(value) = f_n {
239                        v = f - value;
240                        if n * u < v || v < n * t {
241                            n *= cut_back
242                        } else {
243                            return Ok(n);
244                        }
245                    } else {
246                        n *= cut_back
247                    }
248                }
249                Err(LineSearchError::MaximumStepsReached(
250                    format!("{self:?}"),
251                    *max_steps,
252                ))
253            }
254            Self::Wolfe {
255                control_1,
256                control_2,
257                cut_back,
258                max_steps,
259                strong,
260            } => {
261                let mut f_n;
262                let mut j_n;
263                let t_1 = control_1 * m;
264                let t_2 = control_2 * m;
265                let mut trial_argument = trial(n);
266                for _ in 0..*max_steps {
267                    f_n = function(&trial_argument, n);
268                    j_n = jacobian(&trial_argument);
269                    if let Ok(f_val) = f_n
270                        && let Ok(j_val) = j_n
271                        && f - f_val >= n * t_1
272                        && if *strong {
273                            j_val.erase().full_contraction(decrement.erase()).abs() < t_2.abs()
274                        } else {
275                            j_val.erase().full_contraction(decrement.erase()) < t_2
276                        }
277                    {
278                        return Ok(n);
279                    } else {
280                        n *= cut_back;
281                        trial_argument = trial(n)
282                    }
283                }
284                Err(LineSearchError::MaximumStepsReached(
285                    format!("{self:?}"),
286                    *max_steps,
287                ))
288            }
289            Self::None => {
290                panic!("Cannot call backtracking line search when there is no algorithm.")
291            }
292        }
293    }
294}
295
296/// Possible errors encountered during line search.
297pub enum LineSearchError {
298    InvalidStartingPoint(String),
299    MaximumStepsReached(String, usize),
300    NegativeStepSize(String, Scalar),
301    NotDescentDirection(String),
302}
303
304impl StyledError for LineSearchError {
305    fn message(&self, style: &Style) -> String {
306        let (h, c) = (style.headline, style.frame);
307        match self {
308            Self::InvalidStartingPoint(line_search) => format!(
309                "{h}Starting point is invalid.{c}\n\
310                In line search: {line_search}."
311            ),
312            Self::MaximumStepsReached(line_search, steps) => format!(
313                "{h}Maximum number of steps ({steps}) reached.{c}\n\
314                In line search: {line_search}."
315            ),
316            Self::NegativeStepSize(line_search, step_size) => format!(
317                "{h}Negative step size ({step_size}) encountered.{c}\n\
318                In line search: {line_search}."
319            ),
320            Self::NotDescentDirection(line_search) => format!(
321                "{h}Direction is not a descent direction.{c}\n\
322                In line search: {line_search}."
323            ),
324        }
325    }
326}
327
328styled_error!(LineSearchError);