basin 3.0.0

Numerical optimization in pure Rust, with pluggable linear-algebra backends and WASM support.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125

use crate::core::math::{MatrixIdentity, MatrixIndex, Scalar, VectorLen};

/// Constructs BFGS's inverse Hessian at the start of each fresh solve.
///
/// [`Bfgs`](super::Bfgs) calls [`initialize`](Self::initialize) once, after
/// evaluating the starting point's cost and gradient, including on point warm
/// starts. Exact continuation skips this call and retains the current model
/// and the first-update rescaling flag.
///
/// Implementations must return a finite, symmetric positive-definite `n × n`
/// matrix, where `n` is the parameter dimension. They must construct a fresh
/// model on each call rather than reuse an evolving approximation. An empty
/// `0 × 0` matrix is valid for an empty parameter vector.
///
/// Closures and functions `Fn(&V, &V) -> M` implement this trait when
/// `V: VectorLen` and `M: MatrixIndex<F>`. They receive the point and its
/// gradient, in that order. Their returned matrices are checked for shape,
/// finite entries, exact symmetry, and positive definiteness using a pure-Rust
/// Cholesky factorization. Invalid matrices panic during initialization.
/// A named implementation can support custom storage without entry access,
/// and is responsible for enforcing the matrix contract itself.
///
/// With `serde`, a named implementation may derive serialization. Closures
/// support owned checkpoints without serialization or cloning.
pub trait InitialInverseHessian<V, M, F: Scalar = f64> {
    /// Construct a fresh model from the evaluated starting point and gradient.
    fn initialize(&self, param: &V, gradient: &V) -> M;

    /// Whether to replace the model with `(sᵀy / yᵀy) I` immediately before
    /// the first accepted curvature update. Defaults to `false` to preserve
    /// custom scaling. The default identity initializer returns `true`.
    fn rescale_on_first_update(&self) -> bool {
        false
    }
}

/// BFGS's default identity seed with rescaling on the first curvature update.
///
/// Requires only the matrix's existing [`MatrixIdentity`] capability.
#[derive(Clone, Copy, Debug, Default)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct IdentityInverseHessian;

impl<V: VectorLen, M: MatrixIdentity, F: Scalar> InitialInverseHessian<V, M, F>
    for IdentityInverseHessian
{
    fn initialize(&self, param: &V, _: &V) -> M {
        M::identity(param.vec_len())
    }

    fn rescale_on_first_update(&self) -> bool {
        true
    }
}

/// A fixed seed configured by [`Bfgs::with_initial_inverse_hessian`](super::Bfgs::with_initial_inverse_hessian).
///
/// Each fresh solve checks the seed and clones it into the working model.
/// BFGS updates never modify the stored seed. Requires `M: Clone +
/// MatrixIndex<F>` only when this initializer is selected.
#[derive(Clone, Debug)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct FixedInverseHessian<M>(pub(super) M);

impl<V: VectorLen, M: Clone + MatrixIndex<F>, F: Scalar>
    InitialInverseHessian<V, M, F> for FixedInverseHessian<M>
{
    fn initialize(&self, param: &V, _: &V) -> M {
        validate::<M, F>(&self.0, param.vec_len());
        self.0.clone()
    }
}

impl<V: VectorLen, M: MatrixIndex<F>, F: Scalar, Init: Fn(&V, &V) -> M>
    InitialInverseHessian<V, M, F> for Init
{
    fn initialize(&self, param: &V, gradient: &V) -> M {
        let matrix = self(param, gradient);
        validate::<M, F>(&matrix, param.vec_len());
        matrix
    }
}

fn validate<M: MatrixIndex<F>, F: Scalar>(matrix: &M, n: usize) {
    assert!(
        matrix.matrix_rows() == n && matrix.matrix_cols() == n,
        "BFGS initial inverse Hessian must match the parameter dimension"
    );
    // Check both triangles because matvec uses the whole matrix, unlike a
    // factorization that can treat just one triangle as authoritative.
    for i in 0..n {
        for j in 0..=i {
            let a = matrix.matrix_entry(i, j);
            let b = matrix.matrix_entry(j, i);
            assert!(
                a.is_finite() && b.is_finite() && a == b,
                "BFGS initial inverse Hessian must be finite and symmetric"
            );
        }
    }

    // Validate custom seeds without requiring backend factorization support.
    let mut l = vec![F::zero(); n * n];
    for i in 0..n {
        for j in 0..=i {
            let mut value = matrix.matrix_entry(i, j);
            for k in 0..j {
                value = value - l[i * n + k] * l[j * n + k];
            }
            assert!(
                value.is_finite() && (i != j || value > F::zero()),
                "BFGS initial inverse Hessian must be positive definite"
            );
            l[i * n + j] = if i == j {
                value.sqrt()
            } else {
                value / l[j * n + j]
            };
            assert!(
                l[i * n + j].is_finite(),
                "BFGS initial inverse Hessian factorization must be finite"
            );
        }
    }
}