Skip to main content

lbfgs/
lib.rs

1//! # lbfgs
2//! ```
3//! use lbfgs::*;
4//!
5//! fn main() {
6//!     // Problem size and the number of stored vectors in L-BFGS cannot be zero
7//!     let problem_size = 3;
8//!     let lbfgs_memory_size = 5;
9//!
10//!     // Create the L-BFGS instance with curvature and C-BFGS checks enabled
11//!     let mut lbfgs = Lbfgs::<f64>::new(problem_size, lbfgs_memory_size)
12//!         .with_sy_epsilon(1e-8)     // L-BFGS acceptance condition on s'*y > sy_espsilon
13//!         .with_cbfgs_alpha(1.0)     // C-BFGS condition:
14//!         .with_cbfgs_epsilon(1e-4); // y'*s/||s||^2 > epsilon * ||grad(x)||^alpha
15//!
16//!     // Starting value is always accepted (no s or y vectors yet)
17//!     assert_eq!(
18//!         lbfgs.update_hessian(&[0.0, 0.0, 0.0], &[0.0, 0.0, 0.0]),
19//!         UpdateStatus::UpdateOk
20//!     );
21//!
22//!     // Rejected because of CBFGS condition
23//!     assert_eq!(
24//!         lbfgs.update_hessian(&[-0.838, 0.260, 0.479], &[-0.5, 0.6, -1.2]),
25//!         UpdateStatus::Rejection
26//!     );
27//!
28//!     // This will fail because y'*s == 0 (curvature condition)
29//!     assert_eq!(
30//!         lbfgs.update_hessian(
31//!             &[-0.5, 0.6, -1.2],
32//!             &[0.419058177461747, 0.869843029576958, 0.260313940846084]
33//!         ),
34//!         UpdateStatus::Rejection
35//!     );
36//!
37//!     // A proper update that will be accepted
38//!     assert_eq!(
39//!         lbfgs.update_hessian(&[-0.5, 0.6, -1.2], &[0.1, 0.2, -0.3]),
40//!         UpdateStatus::UpdateOk
41//!     );
42//!
43//!     // Apply Hessian approximation on a gradient
44//!     let mut g = [-3.1, 1.5, 2.1];
45//!     let correct_dir = [-1.100601247872944, -0.086568349404424, 0.948633011911515];
46//!
47//!     lbfgs.apply_hessian(&mut g);
48//!
49//!     assert!((g[0] - correct_dir[0]).abs() < 1e-12);
50//!     assert!((g[1] - correct_dir[1]).abs() < 1e-12);
51//!     assert!((g[2] - correct_dir[2]).abs() < 1e-12);
52//! }
53//! ```
54//!
55//! # Errors
56//!
57//! `update_hessian` will give errors if the C-BFGS or L-BFGS curvature conditions are not met.
58//!
59//! # Panics
60//!
61//! `with_sy_epsilon`, `with_cbfgs_alpha`, and `with_cbfgs_epsilon` will panic if given negative
62//! values.
63//!
64//! `update_hessian` and `apply_hessian` will panic if given slices which are not the same length
65//! as the `problem_size`.
66//!
67
68use num_traits::Float;
69
70pub mod vec_ops;
71
72#[cfg(test)]
73mod tests;
74
75/// Precision is a trait extending `num_traits::Float` to provide type-specific constants
76/// (the default sy tolerance).
77pub trait LbfgsPrecision: Float {
78    /// Default s-y tolerance for L-BFGS updates
79    const DEFAULT_SY_TOLERANCE: Self;
80    /// Absolute tolerance (e.g., for unit tests)
81    const ABS_TOL: Self;
82    /// Relative tolerance (e.g., for unit tests)
83    const REL_TOL: Self;
84}
85
86impl LbfgsPrecision for f64 {
87    const DEFAULT_SY_TOLERANCE: f64 = 1e-10;
88    const ABS_TOL: f64 = 1e-8;
89    const REL_TOL: f64 = 1e-10;
90}
91
92impl LbfgsPrecision for f32 {
93    const DEFAULT_SY_TOLERANCE: f32 = 1e-8;
94    const ABS_TOL: f32 = 1e-5;
95    const REL_TOL: f32 = 1e-5;
96}
97
98/// LBFGS Buffer
99///
100/// The Limited-memory BFGS algorithm is used to estimate curvature information for the
101/// gradient of a function as well as other operators and is often used in numerical
102/// optimization and numerical methods in general.
103///
104/// `Lbfgs` maintains a buffer of pairs `(s,y)` and values `rho` (inverse of inner products
105/// of `s` and `y`)
106///
107///
108#[derive(Debug)]
109pub struct Lbfgs<T = f64>
110where
111    T: LbfgsPrecision + std::iter::Sum<T>,
112{
113    /// The number of vectors in s and y that are currently in use
114    active_size: usize,
115    /// Used to warm-start the Hessian estimation with H_0 = gamma * I
116    gamma: T,
117    /// s holds the vectors of state difference s_k = x_{k+1} - x_k, s_0 holds the most recent s
118    s: Vec<Vec<T>>,
119    /// y holds the vectors of the function g (usually cost function gradient) difference:
120    /// y_k = g_{k+1} - g_k, y_0 holds the most recent y
121    y: Vec<Vec<T>>,
122    /// Intermediary storage for the forward L-BFGS pass
123    alpha: Vec<T>,
124    /// Intermediary storage for the forward L-BFGS pass
125    rho: Vec<T>,
126    /// The alpha parameter of the C-BFGS criterion
127    cbfgs_alpha: T,
128    /// The epsilon parameter of the C-BFGS criterion
129    cbfgs_epsilon: T,
130    /// Limit on the inner product s'*y for acceptance in the buffer
131    sy_epsilon: T,
132    /// Holds the state of the last `update_hessian`, used to calculate the `s_k` vectors
133    old_state: Vec<T>,
134    /// Holds the g of the last `update_hessian`, used to calculate the `y_k` vectors
135    old_g: Vec<T>,
136    /// Check to see if the `old_*` variables have valid data
137    first_old: bool,
138}
139
140#[derive(Debug, Copy, Clone, PartialEq)]
141pub enum UpdateStatus {
142    /// The g and state was accepted to update the Hessian estimate
143    UpdateOk,
144    /// The g and state was rejected by the C-BFGS criteria
145    Rejection,
146}
147
148impl<T> Lbfgs<T>
149where
150    T: LbfgsPrecision + std::iter::Sum<T>,
151{
152    /// Create a new L-BFGS instance with a specific problem and L-BFGS buffer size
153    pub fn new(problem_size: usize, buffer_size: usize) -> Lbfgs<T> {
154        debug_assert!(problem_size > 0);
155        debug_assert!(buffer_size > 0);
156
157        Lbfgs {
158            active_size: 0,
159            gamma: T::one(),
160            s: vec![vec![T::zero(); problem_size]; buffer_size + 1], // +1 for the temporary checking area
161            y: vec![vec![T::zero(); problem_size]; buffer_size + 1], // +1 for the temporary checking area
162            alpha: vec![T::zero(); buffer_size],
163            rho: vec![T::zero(); buffer_size + 1],
164            cbfgs_alpha: T::zero(),
165            cbfgs_epsilon: T::zero(),
166            sy_epsilon: T::DEFAULT_SY_TOLERANCE,
167            old_state: vec![T::zero(); problem_size],
168            old_g: vec![T::zero(); problem_size],
169            first_old: true,
170        }
171    }
172
173    /// Update the default C-BFGS alpha
174    pub fn with_cbfgs_alpha(mut self, alpha: T) -> Self {
175        debug_assert!(alpha >= T::zero(), "Negative alpha");
176
177        self.cbfgs_alpha = alpha;
178        self
179    }
180
181    /// Update the default C-BFGS epsilon
182    pub fn with_cbfgs_epsilon(mut self, epsilon: T) -> Self {
183        debug_assert!(epsilon >= T::zero(), "sy_epsilon must be non-negative");
184
185        self.cbfgs_epsilon = epsilon;
186        self
187    }
188
189    /// Update the default sy_epsilon
190    pub fn with_sy_epsilon(mut self, sy_epsilon: T) -> Self {
191        debug_assert!(sy_epsilon >= T::zero(), "sy_epsilon must be non-negative");
192
193        self.sy_epsilon = sy_epsilon;
194        self
195    }
196
197    /// "Empties" the buffer
198    ///
199    /// This is a cheap operation as it amounts to setting certain internal flags
200    pub fn reset(&mut self) {
201        self.active_size = 0;
202        self.first_old = true;
203    }
204
205    /// Apply the current Hessian estimate to an input vector
206    pub fn apply_hessian(&mut self, g: &mut [T]) {
207        debug_assert!(g.len() == self.old_g.len());
208
209        if self.active_size == 0 {
210            // No Hessian available, the g is the best we can do for now
211            return;
212        }
213
214        let active_s = &self.s[0..self.active_size];
215        let active_y = &self.y[0..self.active_size];
216        let rho = &self.rho[0..self.active_size];
217        let alpha = &mut self.alpha;
218
219        let q = g;
220
221        // Perform the forward L-BFGS algorithm
222        for (s_k, (y_k, (rho_k, alpha_k))) in active_s
223            .iter()
224            .zip(active_y.iter().zip(rho.iter().zip(alpha.iter_mut())))
225        {
226            let a = *rho_k * vec_ops::inner_product(s_k, q);
227
228            *alpha_k = a;
229
230            vec_ops::inplace_vec_add(q, y_k, -a);
231        }
232
233        // Apply the initial Hessian estimate and form r = H_0 * q, where H_0 = gamma * I
234        vec_ops::scalar_mult(q, self.gamma);
235        let r = q;
236
237        // Perform the backward L-BFGS algorithm
238        for (s_k, (y_k, (rho_k, alpha_k))) in active_s
239            .iter()
240            .zip(active_y.iter().zip(rho.iter().zip(alpha.iter())))
241            .rev()
242        {
243            let beta = *rho_k * vec_ops::inner_product(y_k, r);
244
245            vec_ops::inplace_vec_add(r, s_k, *alpha_k - beta);
246        }
247
248        // The g with the Hessian applied is available in the input g
249        // r = H_k * grad f
250    }
251
252    /// Check the validity of the newly added s and y vectors. Based on the condition in:
253    /// D.-H. Li and M. Fukushima, "On the global convergence of the BFGS method for nonconvex
254    /// unconstrained optimization problems," vol. 11, no. 4, pp. 1054–1064, jan 2001.
255    fn new_s_and_y_valid(&mut self, g: &[T]) -> bool {
256        let s = self.s.last().unwrap();
257        let y = self.y.last().unwrap();
258        let rho = self.rho.last_mut().unwrap();
259        let ys = vec_ops::inner_product(s, y);
260        let norm_s_squared = vec_ops::inner_product(s, s);
261
262        *rho = T::one() / ys;
263
264        if norm_s_squared <= T::min_positive_value()
265            || (self.sy_epsilon > T::zero() && ys <= self.sy_epsilon)
266        {
267            // In classic L-BFGS, the buffer should be updated only if
268            // y'*s is strictly positive and |s| is nonzero
269            false
270        } else if self.cbfgs_epsilon > T::zero() && self.cbfgs_alpha > T::zero() {
271            // Check the CBFGS condition of Li and Fukushima
272            // Condition: (y^T * s) / ||s||^2 > epsilon * ||grad(x)||^alpha
273            let lhs_cbfgs = ys / norm_s_squared;
274            let rhs_cbfgs = self.cbfgs_epsilon * vec_ops::norm2(g).powf(self.cbfgs_alpha);
275
276            lhs_cbfgs > rhs_cbfgs
277        } else {
278            // The standard L-BFGS conditions are satisfied and C-BFGS is
279            // not active (either cbfgs_epsilon <= 0.0 or cbfgs_alpha <= 0.0)
280            true
281        }
282    }
283
284    /// Saves vectors to update the Hessian estimate
285    pub fn update_hessian(&mut self, g: &[T], state: &[T]) -> UpdateStatus {
286        debug_assert!(g.len() == self.old_state.len());
287        debug_assert!(state.len() == self.old_state.len());
288
289        // First iteration, only save
290        if self.first_old {
291            self.first_old = false;
292
293            self.old_state.copy_from_slice(state);
294            self.old_g.copy_from_slice(g);
295
296            return UpdateStatus::UpdateOk;
297        }
298
299        // Form the new s_k in the temporary area
300        vec_ops::difference_and_save(self.s.last_mut().unwrap(), &state, &self.old_state);
301
302        // Form the new y_k in the temporary area
303        vec_ops::difference_and_save(self.y.last_mut().unwrap(), &g, &self.old_g);
304
305        // Check that the s and y are valid to use
306        if !self.new_s_and_y_valid(g) {
307            return UpdateStatus::Rejection;
308        }
309
310        self.old_state.copy_from_slice(state);
311        self.old_g.copy_from_slice(g);
312
313        // Move the new s_0,  y_0 and rho_0 to the front
314        self.s.rotate_right(1);
315        self.y.rotate_right(1);
316        self.rho.rotate_right(1);
317
318        // Update the Hessian estimate
319        self.gamma = (T::one() / self.rho[0]) / vec_ops::inner_product(&self.y[0], &self.y[0]);
320
321        // Update the indexes and number of active, -1 comes from the temporary area used in
322        // the end of s and y to check if they are valid
323        self.active_size = (self.s.len() - 1).min(self.active_size + 1);
324
325        UpdateStatus::UpdateOk
326    }
327}