numeris 0.5.14

Pure-Rust numerical algorithms library — high performance with SIMD support while also supporting no-std for embedded and WASM targets.
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

use crate::matrix::vector::Vector;
use crate::traits::FloatScalar;
use crate::Matrix;

/// Approximate the Jacobian of `f: R^N → R^M` using forward finite differences.
///
/// Uses step size `h_j = sqrt(ε) * max(|x_j|, 1)` for each component,
/// requiring `N + 1` function evaluations (one base evaluation + N perturbed).
///
/// # Example
///
/// ```
/// use numeris::optim::finite_difference_jacobian;
/// use numeris::Vector;
///
/// // f(x) = [x0^2, x0*x1], Jacobian = [[2*x0, 0], [x1, x0]]
/// let x = Vector::from_array([3.0_f64, 4.0]);
/// let j = finite_difference_jacobian(|x: &Vector<f64, 2>| {
///     Vector::from_array([x[0] * x[0], x[0] * x[1]])
/// }, &x);
/// assert!((j[(0, 0)] - 6.0).abs() < 1e-6);
/// assert!((j[(0, 1)] - 0.0).abs() < 1e-6);
/// assert!((j[(1, 0)] - 4.0).abs() < 1e-6);
/// assert!((j[(1, 1)] - 3.0).abs() < 1e-6);
/// ```
pub fn finite_difference_jacobian<T: FloatScalar, const M: usize, const N: usize>(
    mut f: impl FnMut(&Vector<T, N>) -> Vector<T, M>,
    x: &Vector<T, N>,
) -> Matrix<T, M, N> {
    let f0 = f(x);
    crate::fdiff::forward_diff_jacobian(x, &f0, f)
}

/// Approximate the gradient of `f: R^N → R` using forward finite differences.
///
/// Uses step size `h_j = sqrt(ε) * max(|x_j|, 1)` for each component,
/// requiring `N + 1` function evaluations.
///
/// # Example
///
/// ```
/// use numeris::optim::finite_difference_gradient;
/// use numeris::Vector;
///
/// // f(x) = x0^2 + 2*x1^2, grad = [2*x0, 4*x1]
/// let x = Vector::from_array([3.0_f64, 4.0]);
/// let g = finite_difference_gradient(|x: &Vector<f64, 2>| {
///     x[0] * x[0] + x[1] * x[1] * 2.0
/// }, &x);
/// assert!((g[0] - 6.0).abs() < 1e-5);
/// assert!((g[1] - 16.0).abs() < 1e-5);
/// ```
pub fn finite_difference_gradient<T: FloatScalar, const N: usize>(
    mut f: impl FnMut(&Vector<T, N>) -> T,
    x: &Vector<T, N>,
) -> Vector<T, N> {
    let f0 = f(x);
    let mut grad = Vector::<T, N>::zeros();

    for j in 0..N {
        let h = crate::fdiff::fd_step(x[j]);
        let mut x_pert = *x;
        x_pert[j] = x_pert[j] + h;
        grad[j] = (f(&x_pert) - f0) / h;
    }

    grad
}