optimization_solvers/quasi_newton/
bfgs_b.rs

1use super::*;
2
3#[derive(derive_getters::Getters)]
4pub struct BFGSB {
5    approx_inv_hessian: DMatrix<Floating>,
6    x: DVector<Floating>,
7    k: usize,
8    tol: Floating,
9    s_norm: Option<Floating>,
10    y_norm: Option<Floating>,
11    identity: DMatrix<Floating>,
12    lower_bound: DVector<Floating>,
13    upper_bound: DVector<Floating>,
14}
15impl HasBounds for BFGSB {
16    fn lower_bound(&self) -> &DVector<Floating> {
17        &self.lower_bound
18    }
19    fn set_lower_bound(&mut self, lower_bound: DVector<Floating>) {
20        self.lower_bound = lower_bound;
21    }
22    fn set_upper_bound(&mut self, upper_bound: DVector<Floating>) {
23        self.upper_bound = upper_bound;
24    }
25    fn upper_bound(&self) -> &DVector<Floating> {
26        &self.upper_bound
27    }
28}
29
30impl BFGSB {
31    pub fn next_iterate_too_close(&self) -> bool {
32        match self.s_norm() {
33            Some(s) => s < &self.tol,
34            None => false,
35        }
36    }
37    pub fn gradient_next_iterate_too_close(&self) -> bool {
38        match self.y_norm() {
39            Some(y) => y < &self.tol,
40            None => false,
41        }
42    }
43    pub fn new(
44        tol: Floating,
45        x0: DVector<Floating>,
46        lower_bound: DVector<Floating>,
47        upper_bound: DVector<Floating>,
48    ) -> Self {
49        let n = x0.len();
50        let x0 = x0.box_projection(&lower_bound, &upper_bound);
51        let identity = DMatrix::identity(n, n);
52        BFGSB {
53            approx_inv_hessian: identity.clone(),
54            x: x0,
55            k: 0,
56            tol,
57            s_norm: None,
58            y_norm: None,
59            identity,
60            lower_bound,
61            upper_bound,
62        }
63    }
64}
65
66impl ComputeDirection for BFGSB {
67    fn compute_direction(
68        &mut self,
69        eval: &FuncEvalMultivariate,
70    ) -> Result<DVector<Floating>, SolverError> {
71        // Ok(-&self.approx_inv_hessian * eval.g())
72        let direction = &self.x - &self.approx_inv_hessian * eval.g();
73        let direction = direction.box_projection(&self.lower_bound, &self.upper_bound);
74        let direction = direction - &self.x;
75        Ok(direction)
76    }
77}
78
79impl LineSearchSolver for BFGSB {
80    fn k(&self) -> &usize {
81        &self.k
82    }
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_mut(&mut self) -> &mut usize {
90        &mut self.k
91    }
92    fn has_converged(&self, eval: &FuncEvalMultivariate) -> bool {
93        // either the gradient is small or the difference between the iterates is small
94        // eval.g().norm() < self.tol || self.next_iterate_too_close()
95        if self.next_iterate_too_close() {
96            warn!(target: "BFGSB","Minimization completed: next iterate too close");
97            true
98        } else if self.gradient_next_iterate_too_close() {
99            warn!(target: "BFGSB","Minimization completed: gradient next iterate too close");
100            true
101        } else {
102            eval.g().norm() < self.tol
103        }
104    }
105
106    fn update_next_iterate<LS: LineSearch>(
107        &mut self,
108        line_search: &mut LS,
109        eval_x_k: &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 next_iterate = self.xk() + step * direction;
123        debug!(target: "BFGSB", "ITERATE: {} + {} * {} = {}", self.xk(), step, direction, next_iterate);
124
125        let s = &next_iterate - &self.x;
126        self.s_norm = Some(s.norm());
127        let y = oracle(&next_iterate).g() - eval_x_k.g();
128        self.y_norm = Some(y.norm());
129
130        //updating iterate here, and then we will update the inverse hessian (if corrections are not too small)
131        *self.xk_mut() = next_iterate;
132
133        // We update the inverse hessian and the corrections in this hook which is triggered just after the calculation of the next iterate
134
135        if self.next_iterate_too_close() {
136            return Ok(());
137        }
138
139        if self.gradient_next_iterate_too_close() {
140            return Ok(());
141        }
142
143        let ys = &y.dot(&s);
144        let rho = 1.0 / ys;
145        let w_a = &s * &y.transpose();
146        let w_b = &y * &s.transpose();
147        let innovation = &s * &s.transpose();
148        let left_term = self.identity() - (w_a * rho);
149        let right_term = self.identity() - (w_b * rho);
150        self.approx_inv_hessian =
151            (left_term * &self.approx_inv_hessian * right_term) + innovation * rho;
152
153        Ok(())
154    }
155}
156
157#[cfg(test)]
158mod bfgsb_test {
159    use super::*;
160    #[test]
161    pub fn constrained_grad_desc_backtracking() {
162        std::env::set_var("RUST_LOG", "info");
163
164        let _ = Tracer::default()
165            .with_stdout_layer(Some(LogFormat::Normal))
166            .build();
167        let gamma = 999.0;
168        let f_and_g = |x: &DVector<Floating>| -> FuncEvalMultivariate {
169            let f = 0.5 * (x[0].powi(2) + gamma * x[1].powi(2));
170            let g = DVector::from(vec![x[0], gamma * x[1]]);
171            (f, g).into()
172        };
173
174        //bounds p
175        let lower_bounds = DVector::from_vec(vec![-f64::INFINITY, -f64::INFINITY]);
176        let upper_oounds = DVector::from_vec(vec![f64::INFINITY, f64::INFINITY]);
177        // Linesearch builder
178        let alpha = 1e-4;
179        let beta = 0.5;
180        let mut ls = BackTrackingB::new(alpha, beta, lower_bounds.clone(), upper_oounds.clone());
181
182        // Gradient descent builder
183        let tol = 1e-12;
184        let x_0 = DVector::from(vec![180.0, 152.0]);
185        let mut gd = BFGSB::new(tol, x_0, lower_bounds, upper_oounds);
186
187        // Minimization
188        let max_iter_solver = 10000;
189        let max_iter_line_search = 1000;
190
191        gd.minimize(
192            &mut ls,
193            f_and_g,
194            max_iter_solver,
195            max_iter_line_search,
196            None,
197        )
198        .unwrap();
199
200        println!("Iterate: {:?}", gd.xk());
201
202        let eval = f_and_g(gd.xk());
203        println!("Function eval: {:?}", eval);
204        println!(
205            "projected Gradient norm: {:?}",
206            gd.projected_gradient(&eval).infinity_norm()
207        );
208        println!("tol: {:?}", tol);
209
210        let convergence = gd.has_converged(&eval);
211        println!("Convergence: {:?}", convergence);
212    }
213}
optimization_solvers/quasi_newton/bfgs_b.rs

optimization_solvers/quasi_newton/
bfgs_b.rs
