optimization_solvers/newton/
spn.rs

1//implementation of Spectral Projected Gradient method from Birgin, Ernesto & Martínez, José Mario & Raydan, Marcos. (2014). Spectral Projected Gradient Methods: Review and Perspectives. Journal of Statistical Software. 60. 1-21. 10.18637/jss.v060.i03.
2
3// The spectral projected gradient follows the same approach of the projected gradient method, but:
4// - It rescales the descent direction via a safeguarded Barzila-Borwein scalar, that is between the min and max eigenvalues of the average hessian between x_k and x_k + ts_k (hence the ``spectral'' denomination)
5
6// - The algorithm is typically paired with a non-monotone line search (quadratic or cubic) as that one descibed in [Grippo, Lampariello, Lucidi, 1986] because sometimes enforcing the sufficient decrease condition, which is typical in armijo line search, can be too restrictive. Notice that this kind of line-search embeds the monotone line-search by simply setting to 1 the look-back parameter when evaluating the armijo condition.
7
8use super::*;
9
10#[derive(derive_getters::Getters)]
11pub struct SpectralProjectedNewton {
12    grad_tol: Floating,
13    x: DVector<Floating>,
14    k: usize,
15    lower_bound: DVector<Floating>,
16    upper_bound: DVector<Floating>,
17    lambda: Floating,
18    lambda_min: Floating,
19    lambda_max: Floating,
20}
21
22impl SpectralProjectedNewton {
23    pub fn with_lambdas(mut self, lambda_min: Floating, lambda_max: Floating) -> Self {
24        self.lambda_min = lambda_min;
25        self.lambda_max = lambda_max;
26        self
27    }
28    pub fn new(
29        grad_tol: Floating,
30        x0: DVector<Floating>,
31        oracle: &mut impl FnMut(&DVector<Floating>) -> FuncEvalMultivariate,
32        lower_bound: DVector<Floating>,
33        upper_bound: DVector<Floating>,
34    ) -> Self {
35        let x0 = x0.box_projection(&lower_bound, &upper_bound);
36        let lambda_min = 1e-3;
37        let lambda_max = 1e3;
38
39        // we initialize lambda0 as equation 8 from [Birgin, Martínez, Raydan, 2014]
40        let eval0 = oracle(&x0);
41        let direction0 = &x0 - eval0.g();
42        let direction0 = direction0.box_projection(&lower_bound, &upper_bound);
43        let direction0 = direction0 - &x0;
44        let lambda = (1. / direction0.infinity_norm())
45            .min(lambda_max)
46            .max(lambda_min);
47
48        Self {
49            grad_tol,
50            x: x0,
51            k: 0,
52            lower_bound,
53            upper_bound,
54            lambda,
55            lambda_min,
56            lambda_max,
57        }
58    }
59}
60
61impl HasBounds for SpectralProjectedNewton {
62    fn lower_bound(&self) -> &DVector<Floating> {
63        &self.lower_bound
64    }
65    fn set_lower_bound(&mut self, lower_bound: DVector<Floating>) {
66        self.lower_bound = lower_bound;
67    }
68    fn set_upper_bound(&mut self, upper_bound: DVector<Floating>) {
69        self.upper_bound = upper_bound;
70    }
71    fn upper_bound(&self) -> &DVector<Floating> {
72        &self.upper_bound
73    }
74}
75
76impl ComputeDirection for SpectralProjectedNewton {
77    fn compute_direction(
78        &mut self,
79        eval: &FuncEvalMultivariate,
80    ) -> Result<DVector<Floating>, SolverError> {
81        // let direction = &self.x - self.lambda * eval.g();
82        let hessian = eval
83            .hessian()
84            .clone()
85            .expect("Hessian not available in the oracle");
86        let direction = &self.x - self.lambda * &hessian.cholesky().unwrap().solve(eval.g());
87        let direction = direction.box_projection(&self.lower_bound, &self.upper_bound);
88        let direction = direction - &self.x;
89        Ok(direction)
90    }
91}
92
93impl LineSearchSolver for SpectralProjectedNewton {
94    fn has_converged(&self, eval: &FuncEvalMultivariate) -> bool {
95        let projected_gradient = self.projected_gradient(eval);
96        projected_gradient.infinity_norm() < self.grad_tol
97    }
98    fn xk(&self) -> &DVector<Floating> {
99        &self.x
100    }
101    fn xk_mut(&mut self) -> &mut DVector<Floating> {
102        &mut self.x
103    }
104    fn k(&self) -> &usize {
105        &self.k
106    }
107    fn k_mut(&mut self) -> &mut usize {
108        &mut self.k
109    }
110
111    fn update_next_iterate<LS: LineSearch>(
112        &mut self,
113        line_search: &mut LS,
114        eval_x_k: &FuncEvalMultivariate, //eval: &FuncEvalMultivariate,
115        oracle: &mut impl FnMut(&DVector<Floating>) -> FuncEvalMultivariate,
116        direction: &DVector<Floating>,
117        max_iter_line_search: usize,
118    ) -> Result<(), SolverError> {
119        let step = line_search.compute_step_len(
120            self.xk(),
121            eval_x_k,
122            direction,
123            oracle,
124            max_iter_line_search,
125        );
126
127        let xk = self.xk(); //immutable borrow
128
129        debug!(target: "spectral_projected_newton", "ITERATE: {} + {} * {} = {}", xk, step, direction, xk + step * direction);
130
131        let next_iterate = xk + step * direction;
132
133        // we compute the correction terms:
134        let s_k = &next_iterate - xk;
135        let y_k = oracle(&next_iterate).g() - eval_x_k.g();
136
137        *self.xk_mut() = next_iterate;
138
139        // we update the Barzilai-Borwein scalar to be used in the next iteration for computing the descent direction
140        let skyk = s_k.dot(&y_k);
141        if skyk <= 0. {
142            debug!(target: "spectral_projected_newton", "skyk = {} <= 0. Resetting lambda to lambda_max", skyk);
143            self.lambda = self.lambda_max;
144            return Ok(());
145        }
146        let sksk = s_k.dot(&s_k);
147        self.lambda = (sksk / skyk).min(self.lambda_max).max(self.lambda_min);
148        Ok(())
149    }
150}
151
152#[cfg(test)]
153mod spg_test {
154    use super::*;
155    #[test]
156    pub fn constrained_spg_backtracking() {
157        std::env::set_var("RUST_LOG", "info");
158
159        let _ = Tracer::default()
160            .with_stdout_layer(Some(LogFormat::Normal))
161            .build();
162        let gamma = 1e9;
163        let mut f_and_g = |x: &DVector<Floating>| -> FuncEvalMultivariate {
164            let f = 0.5 * (x[0].powi(2) + gamma * x[1].powi(2));
165            let g = DVector::from(vec![x[0], gamma * x[1]]);
166            let hessian = DMatrix::from_iterator(2, 2, vec![1.0, 0.0, 0.0, gamma]);
167            FuncEvalMultivariate::from((f, g)).with_hessian(hessian)
168        };
169
170        //bounds
171        let lower_bounds = DVector::from_vec(vec![-1., 47.]);
172        let upper_oounds = DVector::from_vec(vec![f64::INFINITY, f64::INFINITY]);
173        // Linesearch builder
174        let c1 = 1e-4;
175        let m = 10;
176        // let ls = BackTrackingB::new(alpha, beta, lower_bounds.clone(), upper_oounds.clone());
177        let mut ls = GLLQuadratic::new(c1, m);
178
179        // Gradient descent builder
180        let tol = 1e-12;
181        let x_0 = DVector::from(vec![180.0, 152.0]);
182        let mut gd =
183            SpectralProjectedNewton::new(tol, x_0, &mut f_and_g, lower_bounds, upper_oounds);
184
185        // Minimization
186        let max_iter_solver = 10000;
187        let max_iter_line_search = 1000;
188
189        gd.minimize(
190            &mut ls,
191            f_and_g,
192            max_iter_solver,
193            max_iter_line_search,
194            None,
195        )
196        .unwrap();
197
198        println!("Iterate: {:?}", gd.xk());
199
200        let eval = f_and_g(gd.xk());
201        println!("Function eval: {:?}", eval);
202        println!(
203            "Projected Gradient norm: {:?}",
204            gd.projected_gradient(&eval).norm()
205        );
206        println!("tol: {:?}", tol);
207
208        let convergence = gd.has_converged(&eval);
209        println!("Convergence: {:?}", convergence);
210    }
211}
optimization_solvers/newton/spn.rs

optimization_solvers/newton/
spn.rs
