opensrdk_optimization/
lbfgs.rs

1use crate::vec::Vector as _;
2use crate::{line_search::*, Status};
3use opensrdk_linear_algebra::*;
4use std::collections::LinkedList;
5
6/// # Limited memory BFGS
7///
8/// https://en.wikipedia.org/wiki/Limited-memory_BFGS
9///
10/// - `memory`: Array length of memory
11/// - `delta`: Threshold of the difference between the value of f(x) and it of previous iteration for finishing
12/// - `epsilon`: Threshold of the norm of gradients for finishing
13/// - `max_iter`: Max limitation of iteration count
14/// - `line_search`: Configurations of line search
15pub struct Lbfgs {
16    memory: usize,
17    delta: f64,
18    epsilon: f64,
19    max_iter: usize,
20    line_search: LineSearch,
21}
22
23impl Default for Lbfgs {
24    fn default() -> Self {
25        Self {
26            memory: 8,
27            delta: 1e-6,
28            epsilon: 1e-6,
29            max_iter: 0,
30            line_search: LineSearch::default(),
31        }
32    }
33}
34
35impl Lbfgs {
36    pub fn new(
37        memory: usize,
38        delta: f64,
39        epsilon: f64,
40        max_iter: usize,
41        line_search: LineSearch,
42    ) -> Self {
43        Self {
44            memory,
45            delta,
46            epsilon,
47            max_iter,
48            line_search,
49        }
50    }
51
52    pub fn with_memory(mut self, memory: usize) -> Self {
53        self.memory = memory;
54
55        self
56    }
57
58    pub fn with_delta(mut self, delta: f64) -> Self {
59        assert!(delta.is_sign_positive(), "delta must be positive");
60
61        self.delta = delta;
62
63        self
64    }
65
66    pub fn with_epsilon(mut self, epsilon: f64) -> Self {
67        assert!(epsilon.is_sign_positive(), "epsilon must be positive");
68
69        self.epsilon = epsilon;
70
71        self
72    }
73
74    pub fn with_max_iter(mut self, max_iter: usize) -> Self {
75        self.max_iter = max_iter;
76
77        self
78    }
79
80    pub fn with_line_search(mut self, line_search: LineSearch) -> Self {
81        self.line_search = line_search;
82
83        self
84    }
85
86    pub fn minimize(&self, x: &mut [f64], fx_gfx: &dyn Fn(&[f64]) -> (f64, Vec<f64>)) -> Status {
87        let mut x_prev = Matrix::new(x.len(), 1);
88        let fx_prev = 0.0;
89        let mut gfx_prev = Matrix::new(x.len(), 1);
90
91        let mut s_y_rho_inv = LinkedList::<(Matrix, Matrix, f64)>::new();
92
93        for k in 1.. {
94            if self.max_iter != 0 && self.max_iter <= k {
95                return Status::MaxIter;
96            }
97
98            let (fx, gfx) = fx_gfx(x);
99            let dfx = fx - fx_prev;
100
101            if dfx.abs() / fx < self.delta {
102                return Status::Delta;
103            }
104
105            if gfx.l2_norm() < self.epsilon {
106                return Status::Epsilon;
107            }
108
109            let sk = x.to_vec().col_mat() - x_prev;
110            let yk = gfx.to_vec().col_mat() - gfx_prev;
111            let rhok_inv = (yk.t() * &sk)[0][0];
112
113            let mut h = gfx.clone().col_mat();
114
115            let mut alpha = vec![0.0; k as usize];
116
117            for (i, (si, yi, rhoi_inv)) in s_y_rho_inv.iter().enumerate().rev() {
118                alpha[i] = (si.t() * &h)[0][0] / rhoi_inv;
119
120                let alphai = alpha[i];
121
122                h = h - alphai * yi.clone();
123            }
124
125            let gamma = (sk.t() * &yk)[0][0] / (yk.t() * &yk)[0][0];
126
127            let mut z: Matrix = gamma * h.clone();
128
129            for (i, (si, yi, rhoi_inv)) in s_y_rho_inv.iter().enumerate() {
130                let alphai = alpha[i];
131                let betai = (yi.t() * &z)[0][0] / rhoi_inv;
132
133                z = z + (alphai - betai) * si.clone();
134            }
135
136            z = -1.0 * z;
137
138            s_y_rho_inv.push_back((sk, yk, rhok_inv));
139            if s_y_rho_inv.len() > self.memory {
140                s_y_rho_inv.pop_front();
141            }
142
143            let step_size = self.line_search.search(fx_gfx, x, z.slice());
144            let dx = step_size * z;
145
146            if !dx.slice().l2_norm().is_finite() {
147                return Status::NaN;
148            }
149
150            gfx_prev = h;
151            x_prev = x.to_vec().col_mat();
152            x.clone_from_slice((&x_prev + dx).slice());
153        }
154
155        Status::Success
156    }
157}
158
159#[cfg(test)]
160mod tests {
161    use crate::{lbfgs::*, numerical_diff};
162    use std::{error::Error, f64::consts::E, f64::consts::PI};
163    #[test]
164    fn it_works() {
165        result().unwrap();
166    }
167
168    #[test]
169    fn result() -> Result<(), Box<dyn Error>> {
170        let mut x = vec![1.2, 3.456789];
171
172        let fx_gfx = |x: &[f64]| {
173            let fx = 1.0 + 2.0 * x[0].powi(2) + 3.0 * x[1].powi(4);
174            let gx = vec![4.0 * x[0], 12.0 * x[1].powi(3)];
175
176            (fx, gx)
177        };
178        Lbfgs::default().minimize(&mut x, &fx_gfx);
179
180        println!("x: {:#?}", x);
181
182        Ok(())
183    }
184
185    #[test]
186    fn result_() -> Result<(), Box<dyn Error>> {
187        let mut x = vec![30.0, -20.0];
188
189        let fx = |x: &[f64]| {
190            let fx = 20.0
191                - 20.0
192                    * (-0.2
193                        * (1.0 / x.len() as f64 * x.into_iter().map(|xi| xi.powi(2)).sum::<f64>())
194                            .sqrt())
195                    .exp()
196                + E
197                - (1.0 / x.len() as f64
198                    * x.into_iter().map(|xi| (2.0 * PI * xi).cos()).sum::<f64>())
199                .exp();
200
201            fx
202        };
203        let fx_gfx = |x: &[f64]| {
204            let gfx = numerical_diff(&fx, x, 1e-3);
205            println!("{:#?}", gfx);
206
207            (fx(x), gfx)
208        };
209        Lbfgs::default().minimize(&mut x, &fx_gfx);
210
211        println!("x: {:#?}", x);
212
213        Ok(())
214    }
215}
opensrdk_optimization/lbfgs.rs

opensrdk_optimization/
lbfgs.rs
