macro_rules! get_vhessian {
($f:ident, $func_name:ident) => { ... };
($f:ident, $func_name:ident, $param_type:ty) => { ... };
}Expand description
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 parameterpthat is passed tof. Defaults to[f64](implying thatfacceptsp: &[f64]).
§Defining f
The multivariate, vector-valued function f must have the following function signature:
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
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);