optimization_solvers/newton/
projected_newton.rs

1use super::*;
2
3#[derive(derive_getters::Getters)]
4pub struct ProjectedNewton {
5    grad_tol: Floating,
6    x: DVector<Floating>,
7    k: usize,
8    lower_bound: DVector<Floating>,
9    upper_bound: DVector<Floating>,
10    s_norm: Option<Floating>,
11    y_norm: Option<Floating>,
12}
13
14impl ProjectedNewton {
15    pub fn next_iterate_too_close(&self) -> bool {
16        match self.s_norm() {
17            Some(s) => s < &self.grad_tol,
18            None => false,
19        }
20    }
21    pub fn gradient_next_iterate_too_close(&self) -> bool {
22        match self.y_norm() {
23            Some(y) => y < &self.grad_tol,
24            None => false,
25        }
26    }
27    pub fn new(
28        grad_tol: Floating,
29        x0: DVector<Floating>,
30        lower_bound: DVector<Floating>,
31        upper_bound: DVector<Floating>,
32    ) -> Self {
33        let x0 = x0.box_projection(&lower_bound, &upper_bound);
34        // let
35        // let pg = DVector::zeros(x0.len());
36        Self {
37            grad_tol,
38            x: x0,
39            k: 0,
40            lower_bound,
41            upper_bound,
42            s_norm: None,
43            y_norm: None,
44            // pg,
45        }
46    }
47}
48
49impl HasBounds for ProjectedNewton {
50    fn lower_bound(&self) -> &DVector<Floating> {
51        &self.lower_bound
52    }
53    fn upper_bound(&self) -> &DVector<Floating> {
54        &self.upper_bound
55    }
56    fn set_lower_bound(&mut self, lower_bound: DVector<Floating>) {
57        self.lower_bound = lower_bound;
58    }
59    fn set_upper_bound(&mut self, upper_bound: DVector<Floating>) {
60        self.upper_bound = upper_bound;
61    }
62}
63
64impl ComputeDirection for ProjectedNewton {
65    fn compute_direction(
66        &mut self,
67        eval: &FuncEvalMultivariate,
68    ) -> Result<DVector<Floating>, SolverError> {
69        // Ok(-eval.g())
70        let hessian = eval
71            .hessian()
72            .clone()
73            .expect("Hessian not available in the oracle");
74        // let direction = &self.x - eval.g();
75        let direction = &self.x - &hessian.cholesky().unwrap().solve(eval.g());
76        let direction = direction.box_projection(&self.lower_bound, &self.upper_bound);
77        let direction = direction - &self.x;
78        Ok(direction)
79    }
80}
81
82impl LineSearchSolver for ProjectedNewton {
83    fn xk(&self) -> &DVector<Floating> {
84        &self.x
85    }
86    fn xk_mut(&mut self) -> &mut DVector<Floating> {
87        &mut self.x
88    }
89    fn k(&self) -> &usize {
90        &self.k
91    }
92    fn k_mut(&mut self) -> &mut usize {
93        &mut self.k
94    }
95    fn has_converged(&self, eval: &FuncEvalMultivariate) -> bool {
96        // we verify that the norm of the gradient is below the tolerance. If the projected gradient is available, then it means that we are in a constrained optimization setting and we verify if it is zero since this is equivalent to first order conditions of optimality in the setting of optimization with simple bounds (Theorem 12.3 from [Neculai Andrei, 2022])
97
98        if self.next_iterate_too_close() {
99            warn!(target: "bfgs","Minimization completed: next iterate too close");
100            true
101        } else if self.gradient_next_iterate_too_close() {
102            warn!(target: "bfgs","Minimization completed: gradient next iterate too close");
103            true
104        } else {
105            let proj_grad = self.projected_gradient(eval);
106            // warn!(target: "projected_newton", "Projected gradient: {:?}", proj_grad);
107            // we compute the infinity norm of the projected gradient
108            proj_grad.infinity_norm() < self.grad_tol
109        }
110    }
111
112    fn update_next_iterate<LS: LineSearch>(
113        &mut self,
114        line_search: &mut LS,
115        eval_x_k: &FuncEvalMultivariate, //eval: &FuncEvalMultivariate,
116        oracle: &mut impl FnMut(&DVector<Floating>) -> FuncEvalMultivariate,
117        direction: &DVector<Floating>,
118        max_iter_line_search: usize,
119    ) -> Result<(), SolverError> {
120        let step = line_search.compute_step_len(
121            self.xk(),
122            eval_x_k,
123            direction,
124            oracle,
125            max_iter_line_search,
126        );
127
128        debug!(target: "projected_newton", "ITERATE: {} + {} * {} = {}", self.xk(), step, direction, self.xk() + step * direction);
129
130        let next_iterate = self.xk() + step * direction;
131
132        let s = &next_iterate - &self.x;
133        self.s_norm = Some(s.norm());
134        let y = oracle(&next_iterate).g() - eval_x_k.g();
135        self.y_norm = Some(y.norm());
136
137        *self.xk_mut() = next_iterate;
138
139        Ok(())
140    }
141}
142
143#[cfg(test)]
144mod projected_newton_tests {
145    use super::*;
146    #[test]
147    pub fn constrained_grad_desc_backtracking() {
148        std::env::set_var("RUST_LOG", "info");
149
150        let _ = Tracer::default()
151            .with_stdout_layer(Some(LogFormat::Normal))
152            .build();
153        let gamma = 90.0;
154        let f_and_g = |x: &DVector<Floating>| -> FuncEvalMultivariate {
155            let f = 0.5 * (x[0].powi(2) + gamma * x[1].powi(2));
156            let g = DVector::from(vec![x[0], gamma * x[1]]);
157            let hessian = DMatrix::from_iterator(2, 2, vec![1.0, 0.0, 0.0, gamma]);
158            FuncEvalMultivariate::from((f, g)).with_hessian(hessian)
159        };
160
161        //bounds p
162        let lower_bounds = DVector::from_vec(vec![-f64::INFINITY, -f64::INFINITY]);
163        let upper_oounds = DVector::from_vec(vec![f64::INFINITY, f64::INFINITY]);
164        // Linesearch builder
165        let alpha = 1e-4;
166        let mut ls = GLLQuadratic::new(alpha, 15);
167
168        // Gradient descent builder
169        let tol = 1e-6;
170        let x_0 = DVector::from(vec![180.0, 152.0]);
171        let mut gd = ProjectedNewton::new(tol, x_0, lower_bounds, upper_oounds);
172
173        // Minimization
174        let max_iter_solver = 10000;
175        let max_iter_line_search = 1000;
176
177        gd.minimize(
178            &mut ls,
179            f_and_g,
180            max_iter_solver,
181            max_iter_line_search,
182            None,
183        )
184        .unwrap();
185
186        println!("Iterate: {:?}", gd.xk());
187
188        let eval = f_and_g(gd.xk());
189        println!("Function eval: {:?}", eval);
190        println!(
191            "projected Gradient norm: {:?}",
192            gd.projected_gradient(&eval).infinity_norm()
193        );
194        println!("tol: {:?}", tol);
195
196        let convergence = gd.has_converged(&eval);
197        println!("Convergence: {:?}", convergence);
198    }
199}
optimization_solvers/newton/projected_newton.rs

optimization_solvers/newton/
projected_newton.rs
