optimization_solvers/newton/
mod.rs

1use super::*;
2
3pub mod projected_newton;
4pub use projected_newton::*;
5pub mod spn;
6pub use spn::*;
7#[derive(derive_getters::Getters)]
8pub struct Newton {
9    tol: Floating,
10    decrement_squared: Option<Floating>,
11    x: DVector<Floating>,
12    k: usize,
13}
14
15impl Newton {
16    pub fn new(tol: Floating, x0: DVector<Floating>) -> Self {
17        Newton {
18            tol,
19            decrement_squared: None,
20            x: x0,
21            k: 0,
22        }
23    }
24}
25
26impl ComputeDirection for Newton {
27    fn compute_direction(
28        &mut self,
29        eval: &FuncEvalMultivariate,
30    ) -> Result<DVector<Floating>, SolverError> {
31        let hessian = eval
32            .hessian()
33            .clone()
34            .expect("Hessian not available in the oracle");
35        //[TODO]: Boyd recommends several alternatives to the solution of Newton system which take advantage of prior information about sparsity/banded bandwidth of the hessian.
36        match hessian.try_inverse() {
37            Some(hessian_inv) => {
38                let direction = -&hessian_inv * eval.g();
39                // we compute also the squared newton decrement
40                self.decrement_squared = Some((hessian_inv * &direction).dot(&direction));
41                Ok(direction)
42            }
43            None => {
44                warn!(target:"newton","Hessian is singular. Using gradient descent direction.");
45                Ok(-eval.g())
46            }
47        }
48    }
49}
50
51impl LineSearchSolver for Newton {
52    fn xk(&self) -> &DVector<Floating> {
53        &self.x
54    }
55    fn k(&self) -> &usize {
56        &self.k
57    }
58    fn xk_mut(&mut self) -> &mut DVector<Floating> {
59        &mut self.x
60    }
61    fn k_mut(&mut self) -> &mut usize {
62        &mut self.k
63    }
64    fn has_converged(&self, _: &FuncEvalMultivariate) -> bool {
65        match self.decrement_squared {
66            Some(decrement_squared) => decrement_squared * 0.5 < self.tol,
67            None => false,
68        }
69    }
70}
71
72#[cfg(test)]
73mod newton_test {
74    use super::*;
75
76    #[test]
77    pub fn newton_morethuente() {
78        std::env::set_var("RUST_LOG", "info");
79
80        let _ = Tracer::default()
81            .with_stdout_layer(Some(LogFormat::Normal))
82            .build();
83        let gamma = 1222.0;
84        let oracle = |x: &DVector<Floating>| -> FuncEvalMultivariate {
85            let f: f64 = 0.5 * (x[0].powi(2) + gamma * x[1].powi(2));
86            let g = DVector::from(vec![x[0], gamma * x[1]]);
87            let hessian = DMatrix::from_iterator(2, 2, vec![1.0, 0.0, 0.0, gamma]);
88            FuncEvalMultivariate::new(f, g).with_hessian(hessian)
89        };
90
91        // Linesearch builder
92
93        let mut ls = MoreThuente::default();
94
95        // newton builder
96        let tol = 1e-8;
97        let x_0 = DVector::from(vec![1.0, 1.0]);
98        let mut nt = Newton::new(tol, x_0);
99
100        // Minimization
101        let max_iter_solver = 1000;
102        let max_iter_line_search = 100;
103
104        nt.minimize(&mut ls, oracle, max_iter_solver, max_iter_line_search, None)
105            .unwrap();
106
107        println!("Iterate: {:?}", nt.xk());
108
109        let eval = oracle(nt.xk());
110        println!("Function eval: {:?}", eval);
111        println!("Gradient norm: {:?}", eval.g().norm());
112        println!("tol: {:?}", tol);
113
114        let convergence = nt.has_converged(&eval);
115        println!("Convergence: {:?}", convergence);
116
117        assert!((eval.f() - 0.0).abs() < 1e-6);
118    }
119
120    #[test]
121    pub fn newton_backtracking() {
122        std::env::set_var("RUST_LOG", "info");
123
124        let _ = Tracer::default()
125            .with_stdout_layer(Some(LogFormat::Normal))
126            .build();
127        let gamma = 1222.0;
128        let oracle = |x: &DVector<Floating>| -> FuncEvalMultivariate {
129            let f: f64 = 0.5 * (x[0].powi(2) + gamma * x[1].powi(2));
130            let g = DVector::from(vec![x[0], gamma * x[1]]);
131            let hessian = DMatrix::from_iterator(2, 2, vec![1.0, 0.0, 0.0, gamma]);
132            FuncEvalMultivariate::new(f, g).with_hessian(hessian)
133        };
134
135        // Linesearch builder
136        let alpha = 1e-4;
137        let beta = 0.5;
138        let mut ls = BackTracking::new(alpha, beta);
139
140        // newton builder
141        let tol = 1e-8;
142        let x_0 = DVector::from(vec![1.0, 1.0]);
143        let mut nt = Newton::new(tol, x_0);
144
145        // Minimization
146        let max_iter_solver = 1000;
147        let max_iter_line_search = 100;
148
149        nt.minimize(&mut ls, oracle, max_iter_solver, max_iter_line_search, None)
150            .unwrap();
151
152        println!("Iterate: {:?}", nt.xk());
153
154        let eval = oracle(nt.xk());
155        println!("Function eval: {:?}", eval);
156        println!("Gradient norm: {:?}", eval.g().norm());
157        println!("tol: {:?}", tol);
158
159        let convergence = nt.has_converged(&eval);
160        println!("Convergence: {:?}", convergence);
161
162        assert!((eval.f() - 0.0).abs() < 1e-6);
163    }
164}
optimization_solvers/newton/mod.rs

optimization_solvers/newton/
mod.rs
