Function hessian 

pub fn hessian<F, G>(
    f: &G,
    x: &Array1<F>,
    epsilon: F,
) -> LinalgResult<Array2<F>>
where
    F: Float + NumAssign + Sum + ScalarOperand,
    G: Fn(&Array1<F>) -> F,

Calculate the Hessian matrix for a scalar-valued function

Computes a numerical approximation of the Hessian matrix (second derivatives) for a function that takes an n-dimensional input and produces a scalar output.

Arguments

• f - A function that maps from R^n to R
• x - The point at which to evaluate the Hessian
• epsilon - The step size for the finite difference approximation

Returns

• The Hessian matrix of shape (n, n)

Examples

use scirs2_core::ndarray::{array, Array1};
use scirs2_linalg::gradient::hessian;

// Define a simple quadratic function: f(x,y) = x^2 + xy + 2y^2
let f = |v: &Array1<f64>| -> f64 {
    let x = v[0];
    let y = v[1];
    x*x + x*y + 2.0*y*y
};

let x = array![1.0, 2.0];  // Point at which to evaluate the Hessian
let epsilon = 1e-5;

let hess = hessian(&f, &x, epsilon).unwrap();

// Analytical Hessian is:
// [∂²f/∂x², ∂²f/∂x∂y]   [2, 1]
// [∂²f/∂y∂x, ∂²f/∂y²] = [1, 4]

assert!((hess[[0, 0]] - 2.0).abs() < 1e-4);
assert!((hess[[0, 1]] - 1.0).abs() < 1e-4);
assert!((hess[[1, 0]] - 1.0).abs() < 1e-4);
assert!((hess[[1, 1]] - 4.0).abs() < 1e-4);