optimization_solvers/steepest_descent/
spg.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 SpectralProjectedGradient {
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 SpectralProjectedGradient {
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 SpectralProjectedGradient {
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 SpectralProjectedGradient {
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 direction = direction.box_projection(&self.lower_bound, &self.upper_bound);
83        let direction = direction - &self.x;
84        Ok(direction)
85    }
86}
87
88impl LineSearchSolver for SpectralProjectedGradient {
89    fn has_converged(&self, eval: &FuncEvalMultivariate) -> bool {
90        let projected_gradient = self.projected_gradient(eval);
91        projected_gradient.infinity_norm() < self.grad_tol
92    }
93    fn xk(&self) -> &DVector<Floating> {
94        &self.x
95    }
96    fn xk_mut(&mut self) -> &mut DVector<Floating> {
97        &mut self.x
98    }
99    fn k(&self) -> &usize {
100        &self.k
101    }
102    fn k_mut(&mut self) -> &mut usize {
103        &mut self.k
104    }
105
106    fn update_next_iterate<LS: LineSearch>(
107        &mut self,
108        line_search: &mut LS,
109        eval_x_k: &FuncEvalMultivariate, //eval: &FuncEvalMultivariate,
110        oracle: &mut impl FnMut(&DVector<Floating>) -> FuncEvalMultivariate,
111        direction: &DVector<Floating>,
112        max_iter_line_search: usize,
113    ) -> Result<(), SolverError> {
114        let step = line_search.compute_step_len(
115            self.xk(),
116            eval_x_k,
117            direction,
118            oracle,
119            max_iter_line_search,
120        );
121
122        let xk = self.xk(); //immutable borrow
123
124        debug!(target: "spectral_projected_gradient", "ITERATE: {} + {} * {} = {}", xk, step, direction, xk + step * direction);
125
126        let next_iterate = xk + step * direction;
127
128        // we compute the correction terms:
129        let s_k = &next_iterate - xk;
130        let y_k = oracle(&next_iterate).g() - eval_x_k.g();
131
132        *self.xk_mut() = next_iterate;
133
134        // we update the Barzilai-Borwein scalar to be used in the next iteration for computing the descent direction
135        let skyk = s_k.dot(&y_k);
136        if skyk <= 0. {
137            debug!(target: "spectral_projected_gradient", "skyk = {} <= 0. Resetting lambda to lambda_max", skyk);
138            self.lambda = self.lambda_max;
139            return Ok(());
140        }
141        let sksk = s_k.dot(&s_k);
142        self.lambda = (sksk / skyk).min(self.lambda_max).max(self.lambda_min);
143        Ok(())
144    }
145}
146
147#[cfg(test)]
148mod spg_test {
149    use super::*;
150    #[test]
151    pub fn constrained_spg_backtracking() {
152        std::env::set_var("RUST_LOG", "info");
153
154        let _ = Tracer::default()
155            .with_stdout_layer(Some(LogFormat::Normal))
156            .build();
157        let gamma = 1e9;
158        let mut f_and_g = |x: &DVector<Floating>| -> FuncEvalMultivariate {
159            let f = 0.5 * (x[0].powi(2) + gamma * x[1].powi(2));
160            let g = DVector::from(vec![x[0], gamma * x[1]]);
161            (f, g).into()
162        };
163
164        //bounds
165        let lower_bounds = DVector::from_vec(vec![-1., 47.]);
166        let upper_oounds = DVector::from_vec(vec![f64::INFINITY, f64::INFINITY]);
167        // Linesearch builder
168        let c1 = 1e-4;
169        let m = 10;
170        // let ls = BackTrackingB::new(alpha, beta, lower_bounds.clone(), upper_oounds.clone());
171        let mut ls = GLLQuadratic::new(c1, m);
172
173        // Gradient descent builder
174        let tol = 1e-12;
175        let x_0 = DVector::from(vec![180.0, 152.0]);
176        let mut gd =
177            SpectralProjectedGradient::new(tol, x_0, &mut f_and_g, lower_bounds, upper_oounds);
178
179        // Minimization
180        let max_iter_solver = 10000;
181        let max_iter_line_search = 1000;
182
183        gd.minimize(
184            &mut ls,
185            f_and_g,
186            max_iter_solver,
187            max_iter_line_search,
188            None,
189        )
190        .unwrap();
191
192        println!("Iterate: {:?}", gd.xk());
193
194        let eval = f_and_g(gd.xk());
195        println!("Function eval: {:?}", eval);
196        println!(
197            "Projected Gradient norm: {:?}",
198            gd.projected_gradient(&eval).norm()
199        );
200        println!("tol: {:?}", tol);
201
202        let convergence = gd.has_converged(&eval);
203        println!("Convergence: {:?}", convergence);
204    }
205}
optimization_solvers/steepest_descent/spg.rs

optimization_solvers/steepest_descent/
spg.rs
