numdiff 0.7.3

/// Get a function that returns the Hessian of the provided multivariate, vector-valued function.
///
/// The Hessian is computed using forward-mode automatic differentiation with hyper-dual numbers.
///
/// # Arguments
///
/// * `f` - Multivariate, vector-valued function, $\mathbf{f}:\mathbb{R}^{n}\to\mathbb{R}^{m}$.
/// * `func_name` - Name of the function that will return the Hessian of $\mathbf{f}(\mathbf{x})$ at
///   any point $\mathbf{x}\in\mathbb{R}^{n}$.
/// * `param_type` (optional) - Type of the extra runtime parameter `p` that is passed to `f`.
///   Defaults to `[f64]` (implying that `f` accepts `p: &[f64]`).
///
/// # Defining `f`
///
/// The multivariate, vector-valued function `f` must have the following function signature:
///
/// ```ignore
/// fn f<S: Scalar, V: Vector<S>>(x: &V, p: &[f64]) -> V::DVectorT<S> {
///     // place function contents here
/// }
/// ```
///
/// For the automatic differentiation to work, `f` must be fully generic over the types of scalars
/// and vectors used. Additionally, the function must return an instance of `V::DVectorT` (a
/// dynamically-sized vector type that is compatible with `V`, see the
/// [`linalg-traits` docs](https://docs.rs/linalg-traits/latest/linalg_traits/trait.Vector.html#associatedtype.DVectorT)
/// for more information) since we did not want to burden the user with having to specify the size
/// of the output vector (i.e. $m$, where $\mathbf{f}:\mathbb{R}^{n}\to\mathbb{R}^{m}$) at compile
/// time, especially since users may be using this crate exclusively with dynamically-sized types.
///
/// # Warning
///
/// `f` cannot be defined as closure. It must be defined as a function.
///
/// # Warning
///
/// This Hessian function generated by this macro will always return a vector of dynamically-sized
/// matrices, even if the function `f` uses statically-sized vectors. This is to avoid needing to
/// pass const generics to this function to define the sizes of the Hessian matrices. Instead, the
/// dimensions are determined at runtime.
///
/// # Note
///
/// The function produced by this macro will perform $\frac{n(n+1)}{2}$ evaluations of
/// $\mathbf{f}(\mathbf{x})$ to evaluate its Hessian. This is more efficient than using finite
/// differences, which would require $\frac{n(n+1)}{2} + 1$ evaluations.
///
/// # Returns
///
/// The returned result is stored as a length-$m$ [`Vec`] of $n\times n$ matrices, where each matrix
/// represents the Hessian of the corresponding component of the vector function.
///
/// # Examples
///
/// ## Basic Example
///
/// Compute the Hessian of
///
/// $$
/// \mathbf{f}(\mathbf{x})=
/// \begin{bmatrix}
///     f_{0}(\mathbf{x}) \\\\
///     f_{1}(\mathbf{x})
/// \end{bmatrix}=
/// \begin{bmatrix}
///     x_{0}^{5}x_{1}+x_{0}\sin^{3}{x_{1}} \\\\
///     x_{0}^{3}+x_{1}^{4}-3x_{0}^{2}x_{1}^{2}
/// \end{bmatrix}
/// $$
///
/// at $\mathbf{x}=(5,8)^{T}$, and compare the result to the true result.
///
/// #### Using standard vectors
///
/// ```
/// use linalg_traits::{Mat, Matrix, Scalar, Vector};
/// use numtest::*;
///
/// use numdiff::{get_vhessian, HyperDual, HyperDualVector};
///
/// // Define the function, f(x).
/// fn f<S: Scalar, V: Vector<S>>(x: &V, _p: &[f64]) -> V::DVectorT<S> {
///     V::DVectorT::from_slice(&[
///         x.vget(0).powi(5) * x.vget(1) + x.vget(0) * x.vget(1).sin().powi(3),
///         x.vget(0).powi(3) + x.vget(1).powi(4) - S::new(3.0) * x.vget(0).powi(2) * x.vget(1).powi(2),
///     ])
/// }
///
/// // Define the evaluation point.
/// let x0 = vec![5.0, 8.0];
///
/// // Parameter vector (empty for this example).
/// let p = [];
///
/// // Autogenerate the function "hess" that can be used to compute the Hessian of f(x) at any point
/// // x.
/// get_vhessian!(f, hess);
///
/// // Evaluate the Hessian using "hess".
/// let hess_eval: Vec<Mat<f64>> = hess(&x0, &p);
///
/// // True Hessian of f(x) at the evaluation point.
/// let hess_f0_true: Mat<f64> = Mat::from_row_slice(
///     2,
///     2,
///     &[
///         20000.0,
///         3125.0 + 3.0 * 8.0_f64.sin().powi(2) * 8.0_f64.cos(),
///         3125.0 + 3.0 * 8.0_f64.sin().powi(2) * 8.0_f64.cos(),
///         30.0 * 8.0_f64.sin() * 8.0_f64.cos().powi(2) - 15.0 * 8.0_f64.sin().powi(3),
///     ],
/// );
/// let hess_f1_true: Mat<f64> = Mat::from_row_slice(2, 2, &[-354.0, -480.0, -480.0, 618.0]);
/// let hess_true: Vec<Mat<f64>> = vec![hess_f0_true, hess_f1_true];
///
/// // Verify that the Hessian function obtained using get_vhessian! computes the Hessian correctly.
/// assert_arrays_equal_to_decimal!(hess_eval[0], hess_true[0], 16);
/// assert_arrays_equal_to_decimal!(hess_eval[1], hess_true[1], 16);
/// ```
///
/// #### Using other vector types
///
/// The function produced by `get_vhessian!` can accept _any_ type for `x0`, as long as it
/// implements the `linalg_traits::Vector` trait.
///
/// ```
/// use faer::Mat;
/// use linalg_traits::{Scalar, Vector};
/// use nalgebra::{dvector, DMatrix, DVector, SMatrix, SVector};
/// use ndarray::{array, Array1, Array2};
///
/// use numdiff::{get_vhessian, HyperDual, HyperDualVector};
///
/// // Define the function, f(x).
/// fn f<S: Scalar, V: Vector<S>>(x: &V, _p: &[f64]) -> V::DVectorT<S> {
///     V::DVectorT::from_slice(&[
///         x.vget(0).powi(5) * x.vget(1) + x.vget(0) * x.vget(1).sin().powi(3),
///         x.vget(0).powi(3) + x.vget(1).powi(4) - S::new(3.0) * x.vget(0).powi(2) * x.vget(1).powi(2),
///     ])
/// }
///
/// // Parameter vector (empty for this example).
/// let p = [];
///
/// // Autogenerate the function "hess" that can be used to compute the Hessian of f(x) at any point
/// // x.
/// get_vhessian!(f, hess);
///
/// // nalgebra::DVector
/// let x0: DVector<f64> = dvector![5.0, 8.0];
/// let hess_eval: Vec<DMatrix<f64>> = hess(&x0, &p);
///
/// // nalgebra::SVector
/// let x0: SVector<f64, 2> = SVector::from_slice(&[5.0, 8.0]);
/// let hess_eval: Vec<DMatrix<f64>> = hess(&x0, &p);
///
/// // ndarray::Array1
/// let x0: Array1<f64> = array![5.0, 8.0];
/// let hess_eval: Vec<Array2<f64>> = hess(&x0, &p);
///
/// // faer::Mat
/// let x0: Mat<f64> = Mat::from_slice(&[5.0, 8.0]);
/// let hess_eval: Vec<Mat<f64>> = hess(&x0, &p);
/// ```
///
/// ## Example Passing Runtime Parameters
///
/// Compute the Hessian of a parameterized vector function
///
/// $$\mathbf{f}(\mathbf{x})=\begin{bmatrix}ax_{0}^{2}x_{1}+bx_{1}^{2}\\\\cx_{0}x_{1}+dx_{0}^{2}\end{bmatrix}$$
///
/// where $a$, $b$, $c$, and $d$ are runtime parameters.
///
/// ```
/// use linalg_traits::{Mat, Matrix, Scalar, Vector};
/// use numtest::*;
///
/// use numdiff::{get_vhessian, HyperDual, HyperDualVector};
///
/// // Define the function, f(x).
/// fn f<S: Scalar, V: Vector<S>>(x: &V, p: &[f64]) -> V::DVectorT<S> {
///     let a = S::new(p[0]);
///     let b = S::new(p[1]);
///     let c = S::new(p[2]);
///     let d = S::new(p[3]);
///     V::DVectorT::from_slice(&[
///         a * x.vget(0).powi(2) * x.vget(1) + b * x.vget(1).powi(2),
///         c * x.vget(0) * x.vget(1) + d * x.vget(0).powi(2),
///     ])
/// }
///
/// // Runtime parameters.
/// let a = 1.5;
/// let b = 2.0;
/// let c = 0.8;
/// let d = 3.0;
/// let p = [a, b, c, d];
///
/// // Evaluation point.
/// let x0 = vec![1.0, -0.5];
///
/// // Autogenerate the Hessian function.
/// get_vhessian!(f, hess);
///
/// // True Hessian functions for both components.
/// let hess_f0_true = Mat::from_row_slice(2, 2, &[
///     2.0 * a * x0[1], 2.0 * a * x0[0],
///     2.0 * a * x0[0], 2.0 * b
/// ]);
/// let hess_f1_true = Mat::from_row_slice(2, 2, &[
///     2.0 * d, c,
///     c, 0.0
/// ]);
/// let hess_true = vec![hess_f0_true, hess_f1_true];
///
/// // Compute the Hessian using the automatically generated function and compare.
/// let hess_eval: Vec<Mat<f64>> = hess(&x0, &p);
/// assert_arrays_equal_to_decimal!(hess_eval[0], hess_true[0], 16);
/// assert_arrays_equal_to_decimal!(hess_eval[1], hess_true[1], 16);
/// ```
///
/// ## Example Passing Custom Parameter Types
///
/// Use a custom parameter struct instead of `f64` values.
///
/// ```
/// use linalg_traits::{Mat, Matrix, Scalar, Vector};
/// use numtest::*;
///
/// use numdiff::{get_vhessian, HyperDual, HyperDualVector};
///
/// struct Data {
///     a: f64,
///     b: f64,
///     c: f64,
///     d: f64,
/// }
///
/// // Define the function, f(x).
/// fn f<S: Scalar, V: Vector<S>>(x: &V, p: &Data) -> V::DVectorT<S> {
///     let a = S::new(p.a);
///     let b = S::new(p.b);
///     let c = S::new(p.c);
///     let d = S::new(p.d);
///     V::DVectorT::from_slice(&[
///         a * x.vget(0).powi(2) * x.vget(1) + b * x.vget(1).powi(2),
///         c * x.vget(0) * x.vget(1) + d * x.vget(0).powi(2),
///     ])
/// }
///
/// // Runtime parameter struct.
/// let p = Data {
///     a: 1.5,
///     b: 2.0,
///     c: 0.8,
///     d: 3.0,
/// };
///
/// // Evaluation point.
/// let x0 = vec![1.0, -0.5];
///
/// // Autogenerate the Hessian function, telling the macro to expect a runtime parameter of type
/// // &Data.
/// get_vhessian!(f, hess, Data);
///
/// // True Hessian functions.
/// let hess_f0_true = Mat::from_row_slice(2, 2, &[
///     2.0 * p.a * x0[1], 2.0 * p.a * x0[0],
///     2.0 * p.a * x0[0], 2.0 * p.b
/// ]);
/// let hess_f1_true = Mat::from_row_slice(2, 2, &[
///     2.0 * p.d, p.c,
///     p.c, 0.0
/// ]);
/// let hess_true = vec![hess_f0_true, hess_f1_true];
///
/// // Compute the Hessian using the automatically generated function and compare.
/// let hess_eval: Vec<Mat<f64>> = hess(&x0, &p);
/// assert_arrays_equal_to_decimal!(hess_eval[0], hess_true[0], 16);
/// assert_arrays_equal_to_decimal!(hess_eval[1], hess_true[1], 16);
/// ```
#[macro_export]
macro_rules! get_vhessian {
    ($f:ident, $func_name:ident) => {
        get_vhessian!($f, $func_name, [f64]);
    };
    ($f:ident, $func_name:ident, $param_type:ty) => {
        /// Hessian of a multivariate, vector-valued function `f: ℝⁿ → ℝᵐ`.
        ///
        /// This function is generated for a specific function `f` using the
        /// `numdiff::get_vhessian!` macro.
        ///
        /// # Arguments
        ///
        /// * `x0` - Evaluation point, `x₀ ∈ ℝⁿ`.
        /// * `p` - Extra runtime parameter. This is a parameter (can be of any arbitrary type)
        ///   defined at runtime that the function may depend on but is not differentiated with
        ///   respect to.
        ///
        /// # Returns
        ///
        /// Hessian of `f` with respect to `x`, evaluated at `x = x₀`.
        ///
        /// Returns a vector of length `m`, where each element is the Hessian matrix
        /// `Hᵢ(x₀) = (∂²fᵢ/∂x²)|ₓ₌ₓ₀ ∈ ℝⁿˣⁿ` for the i-th component of `f`.
        fn $func_name<S, V>(x0: &V, p: &$param_type) -> Vec<V::DMatrixMxNf64>
        where
            S: Scalar,
            V: Vector<S>,
        {
            // Promote the evaluation point to a vector of hyper-dual numbers.
            let x0_hyper_dual = x0.clone().to_hyper_dual_vector();

            // Determine the dimensions.
            let n = x0.len();

            // Evaluate the function once to determine the output dimension.
            let f_result = $f(&x0_hyper_dual, p);
            let m = f_result.len();

            // Initialize a vector of matrices to store the Hessian, where each matrix is filled
            // with zeros.
            let mut hess = vec![x0.new_dmatrix_m_by_n_f64(n); m];

            // Variable to store the evaluation point perturbed with hyper-dual directions.
            let mut x_perturbed;

            // Evaluate the Hessian, iterating over the upper triangular elements.
            for i in 0..n {
                for j in i..n {
                    // Make a copy of the hyper-dual evaluation point.
                    x_perturbed = x0_hyper_dual.clone();

                    // Diagonal element (∂²f/∂xᵢ²).
                    if i == j {
                        // Set both ε₁ and ε₂ coefficients for variable i to 1.0.
                        let original = x_perturbed.vget(i);
                        x_perturbed.vset(i, HyperDual::new(original.get_a(), 1.0, 1.0, 0.0));
                    }
                    // Off-diagonal element (∂²f/∂xᵢ∂xⱼ).
                    else {
                        // Set ε₁ coefficient for variable i and ε₂ coefficient for variable j to
                        // 1.0.
                        let original_i = x_perturbed.vget(i);
                        let original_j = x_perturbed.vget(j);
                        x_perturbed.vset(i, HyperDual::new(original_i.get_a(), 1.0, 0.0, 0.0));
                        x_perturbed.vset(j, HyperDual::new(original_j.get_a(), 0.0, 1.0, 0.0));
                    }

                    // Evaluate the function at the perturbed point.
                    let f_result = $f(&x_perturbed, p);

                    // Extract the second-order derivatives for each component and store them.
                    for k in 0..m {
                        let f_k = f_result.vget(k);

                        // Extract the second-order derivative from the hyper-dual result.
                        let second_derivative = f_k.get_d();

                        // Store the result in both (i,j) and (j,i) positions (Hessian is
                        // symmetric).
                        hess[k][(i, j)] = second_derivative;
                        hess[k][(j, i)] = second_derivative;
                    }
                }
            }

            hess
        }
    };
}

#[cfg(test)]
mod tests {
    use crate::{HyperDual, HyperDualVector};
    use linalg_traits::{Mat, Matrix, Scalar, Vector};
    use nalgebra::{DMatrix, DVector, SMatrix, SVector, dvector};
    use ndarray::{Array1, Array2, array};
    use numtest::*;

    #[test]
    fn test_vhessian_1() {
        // Function to take the Hessian of.
        fn f<S: Scalar, V: Vector<S>>(x: &V, _p: &[f64]) -> V::DVectorT<S> {
            V::DVectorT::from_slice(&[x.vget(0).powi(3)])
        }

        // Evaluation point.
        let x0 = vec![2.0];

        // Parameter vector (empty for this test).
        let p = [];

        // True Hessian function.
        let hess = |x: &Vec<f64>| vec![Mat::from_row_slice(1, 1, &[6.0 * x[0]])];

        // Hessian function obtained via forward-mode automatic differentation.
        get_vhessian!(f, hess_autodiff);

        // Evaluate the Hessian using both functions.
        let hess_eval_autodiff: Vec<Mat<f64>> = hess_autodiff(&x0, &p);
        let hess_eval: Vec<Mat<f64>> = hess(&x0);

        // Test autodiff Hessian against true Hessian.
        assert_arrays_equal_to_decimal!(hess_eval_autodiff[0], hess_eval[0], 16);
    }

    #[test]
    fn test_vhessian_2() {
        // Function to take the Hessian of.
        fn f<S: Scalar, V: Vector<S>>(x: &V, _p: &[f64]) -> V::DVectorT<S> {
            V::DVectorT::from_slice(&[x.vget(0).powi(2) + x.vget(1).powi(3)])
        }

        // Evaluation point.
        let x0 = SVector::from_row_slice(&[1.0, 2.0]);

        // Parameter vector (empty for this test).
        let p = [];

        // True Hessian function.
        let hess = |x: &SVector<f64, 2>| {
            vec![SMatrix::<f64, 2, 2>::from_row_slice(&[
                2.0,
                0.0,
                0.0,
                6.0 * x[1],
            ])]
        };

        // Hessian function obtained via forward-mode automatic differentation.
        get_vhessian!(f, hess_autodiff);

        // Evaluate the Hessian using both functions.
        let hess_eval_autodiff: Vec<DMatrix<f64>> = hess_autodiff(&x0, &p);
        let hess_eval: Vec<SMatrix<f64, 2, 2>> = hess(&x0);

        // Test autodiff Hessian against true Hessian.
        assert_arrays_equal_to_decimal!(hess_eval_autodiff[0], hess_eval[0], 16);
    }

    #[test]
    fn test_vhessian_3() {
        // Function to take the Hessian of.
        fn f<S: Scalar, V: Vector<S>>(x: &V, _p: &[f64]) -> V::DVectorT<S> {
            V::DVectorT::from_slice(&[
                x.vget(0).powi(5) * x.vget(1) + x.vget(0) * x.vget(1).sin().powi(3)
            ])
        }

        // Evaluation point.
        let x0 = array![1.0, 2.0];

        // Parameter vector (empty for this test).
        let p = [];

        // True Hessian function.
        let hess = |x: &Array1<f64>| {
            vec![Array2::<f64>::from_row_slice(
                2,
                2,
                &[
                    20.0 * x[0].powi(3) * x[1],
                    5.0 * x[0].powi(4) + 3.0 * x[1].sin().powi(2) * x[1].cos(),
                    5.0 * x[0].powi(4) + 3.0 * x[1].sin().powi(2) * x[1].cos(),
                    6.0 * x[0] * x[1].sin() * x[1].cos().powi(2) - 3.0 * x[0] * x[1].sin().powi(3),
                ],
            )]
        };

        // Hessian function obtained via forward-mode automatic differentation.
        get_vhessian!(f, hess_autodiff);

        // Evaluate the Hessian using both functions.
        let hess_eval_autodiff: Vec<Array2<f64>> = hess_autodiff(&x0, &p);
        let hess_eval: Vec<Array2<f64>> = hess(&x0);

        // Test autodiff Hessian against true Hessian.
        assert_arrays_equal_to_decimal!(hess_eval_autodiff[0], hess_eval[0], 16);
    }

    #[test]
    fn test_vhessian_4() {
        // Function to take the Hessian of.
        fn f<S: Scalar, V: Vector<S>>(x: &V, _p: &[f64]) -> V::DVectorT<S> {
            V::DVectorT::from_slice(&[
                x.vget(0).powi(5) * x.vget(1) + x.vget(0) * x.vget(1).sin().powi(3),
                x.vget(0).powi(3) + x.vget(1).powi(4)
                    - S::new(3.0) * x.vget(0).powi(2) * x.vget(1).powi(2),
            ])
        }

        // Evaluation point.
        let x0 = DVector::<f64>::from_slice(&[1.0, 2.0]);

        // Parameter vector (empty for this test).
        let p = [];

        // True Hessian function.
        let hess = |x: &DVector<f64>| {
            vec![
                DMatrix::<f64>::from_row_slice(
                    2,
                    2,
                    &[
                        20.0 * x[0].powi(3) * x[1],
                        5.0 * x[0].powi(4) + 3.0 * x[1].sin().powi(2) * x[1].cos(),
                        5.0 * x[0].powi(4) + 3.0 * x[1].sin().powi(2) * x[1].cos(),
                        6.0 * x[0] * x[1].sin() * x[1].cos().powi(2)
                            - 3.0 * x[0] * x[1].sin().powi(3),
                    ],
                ),
                DMatrix::<f64>::from_row_slice(
                    2,
                    2,
                    &[
                        6.0 * x[0] - 6.0 * x[1].powi(2),
                        -12.0 * x[0] * x[1],
                        -12.0 * x[0] * x[1],
                        12.0 * x[1].powi(2) - 6.0 * x[0].powi(2),
                    ],
                ),
            ]
        };

        // Hessian function obtained via forward-mode automatic differentation.
        get_vhessian!(f, hess_autodiff);

        // Evaluate the Hessian using both functions.
        let hess_eval_autodiff: Vec<DMatrix<f64>> = hess_autodiff(&x0, &p);
        let hess_eval: Vec<DMatrix<f64>> = hess(&x0);

        // Test autodiff Hessian against true Hessian.
        assert_arrays_equal_to_decimal!(hess_eval_autodiff[0], hess_eval[0], 16);
        assert_arrays_equal_to_decimal!(hess_eval_autodiff[1], hess_eval[1], 16);
    }

    #[test]
    fn test_vhessian_5() {
        // Function to take the Hessian of.
        #[allow(clippy::many_single_char_names)]
        fn f<S: Scalar, V: Vector<S>>(x: &V, p: &[f64]) -> V::DVectorT<S> {
            let a = S::new(p[0]);
            let b = S::new(p[1]);
            let c = S::new(p[2]);
            let d = S::new(p[3]);
            V::DVectorT::from_slice(&[
                a * x.vget(0).powi(2) * x.vget(1) + b * x.vget(1).powi(2),
                c * x.vget(0) * x.vget(1) + d * x.vget(0).powi(2),
            ])
        }

        // Parameter vector.
        let p = [1.5, 2.0, 0.8, 3.0];

        // Evaluation point.
        let x0 = dvector![1.0, -0.5];

        // True Hessian function.
        let hess = |x: &DVector<f64>, p: &[f64]| {
            vec![
                DMatrix::<f64>::from_row_slice(
                    2,
                    2,
                    &[
                        2.0 * p[0] * x[1],
                        2.0 * p[0] * x[0],
                        2.0 * p[0] * x[0],
                        2.0 * p[1],
                    ],
                ),
                DMatrix::<f64>::from_row_slice(2, 2, &[2.0 * p[3], p[2], p[2], 0.0]),
            ]
        };

        // Hessian function obtained via forward-mode automatic differentiation.
        get_vhessian!(f, hess_autodiff);

        // Evaluate the Hessian using both methods.
        let hess_eval_autodiff: Vec<DMatrix<f64>> = hess_autodiff(&x0, &p);
        let hess_eval: Vec<DMatrix<f64>> = hess(&x0, &p);

        // Test autodiff Hessian against true Hessian.
        assert_arrays_equal_to_decimal!(hess_eval_autodiff[0], hess_eval[0], 16);
        assert_arrays_equal_to_decimal!(hess_eval_autodiff[1], hess_eval[1], 16);
    }

    #[test]
    fn test_vhessian_custom_params() {
        struct Data {
            a: f64,
            b: f64,
            c: f64,
            d: f64,
        }

        // Function to take the Hessian of.
        #[allow(clippy::many_single_char_names)]
        fn f<S: Scalar, V: Vector<S>>(x: &V, p: &Data) -> V::DVectorT<S> {
            let a = S::new(p.a);
            let b = S::new(p.b);
            let c = S::new(p.c);
            let d = S::new(p.d);
            V::DVectorT::from_slice(&[
                a * x.vget(0).powi(2) * x.vget(1) + b * x.vget(1).powi(2),
                c * x.vget(0) * x.vget(1) + d * x.vget(0).powi(2),
            ])
        }

        // Runtime parameter struct.
        let p = Data {
            a: 1.5,
            b: 2.0,
            c: 0.8,
            d: 3.0,
        };

        // Evaluation point.
        let x0 = vec![1.0, -0.5];

        // Hessian function obtained via forward-mode automatic differentiation.
        get_vhessian!(f, hess, Data);

        // True Hessian function.
        let hess_f0_true = Mat::from_row_slice(
            2,
            2,
            &[
                2.0 * p.a * x0[1],
                2.0 * p.a * x0[0],
                2.0 * p.a * x0[0],
                2.0 * p.b,
            ],
        );
        let hess_f1_true = Mat::from_row_slice(2, 2, &[2.0 * p.d, p.c, p.c, 0.0]);
        let hess_true = [hess_f0_true, hess_f1_true];

        // Evaluate the Hessian using both functions.
        let hess_eval: Vec<Mat<f64>> = hess(&x0, &p);

        // Test autodiff Hessian against true Hessian.
        assert_arrays_equal_to_decimal!(hess_eval[0], hess_true[0], 16);
        assert_arrays_equal_to_decimal!(hess_eval[1], hess_true[1], 16);
    }
}