macro_rules! get_jacobian {
($f:ident, $func_name:ident) => { ... };
}
Expand description
Get a function that returns the Jacobian of the provided multivariate, vector-valued function.
The Jacobian is computed using forward-mode automatic differentiation.
§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 Jacobian of $\mathbf{f}(\mathbf{x})$ at any point $\mathbf{x}\in\mathbb{R}^{n}$.
§Defining f
The multivariate, vector-valued function f
must have the following function signature:
fn f<S: Scalar, V: Vector<S>>(x: &V) -> 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::DVector
(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 Jacobian function generated by this macro will always return a dynamically-sized matrix,
even if the function f
uses statically-sized vectors. This is to avoid needing to pass a const
generic to this function to define the number of rows ($m$) of the Jacobian. Instead, the number
of rows is determined at runtime.
§Note
The function produced by this macro will perform $n$ evaluations of $\mathbf{f}(\mathbf{x})$ to evaluate its Jacobian.
§Example
Compute the Jacobian of
$$ \mathbf{f}(\mathbf{x})= \begin{bmatrix} x_{0} \\ 5x_{2} \\ 4x_{1}^{2}-2x_{2} \\ x_{2}\sin{x_{0}} \end{bmatrix} $$
at $\mathbf{x}=(5,6,7)^{T}$, and compare the result to the true result of
$$ \mathbf{J}\left((5,6,7)^{T}\right)= \begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & 5 \\ 0 & 48 & -2 \\ 7\cos{(5)} & 0 & \sin{(5)} \end{bmatrix} $$
§Using standard vectors
use linalg_traits::{Mat, Matrix, Scalar, Vector};
use numtest::*;
use numdiff::{get_jacobian, Dual, DualVector};
// Define the function, f(x).
fn f<S: Scalar, V: Vector<S>>(x: &V) -> V::DVectorT<S> {
V::DVectorT::from_slice(&[
x.vget(0),
x.vget(2) * 5.0,
x.vget(1).powi(2) * 4.0 - x.vget(2) * 2.0,
x.vget(2) * x.vget(0).sin(),
])
}
// Define the evaluation point.
let x0 = vec![5.0, 6.0, 7.0];
// Autogenerate the function "jac" that can be used to compute the Jacobian of f(x) at any point
// x.
get_jacobian!(f, jac);
// Evaluate the Jacobian using "jac".
let jac_eval: Mat<f64> = jac(&x0);
// True Jacobian of f(x) at the evaluation point.
let jac_eval_true: Mat<f64> = Mat::from_row_slice(
4,
3,
&[
1.0,
0.0,
0.0,
0.0,
0.0,
5.0,
0.0,
48.0,
-2.0,
7.0 * 5.0_f64.cos(),
0.0,
5.0_f64.sin(),
],
);
// Verify that the Jacobian function obtained using get_jacobian! computes the Jacobian
// correctly.
assert_eq!(jac_eval, jac_eval_true);
§Using other vector types
The function produced by get_jacobian!
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, SVector};
use ndarray::{array, Array1, Array2};
use numdiff::{get_jacobian, Dual, DualVector};
// Define the function, f(x).
fn f<S: Scalar, V: Vector<S>>(x: &V) -> V::DVectorT<S> {
V::DVectorT::from_slice(&[
x.vget(0),
x.vget(2) * 5.0,
x.vget(1).powi(2) * 4.0 - x.vget(2) * 2.0,
x.vget(2) * x.vget(0).sin(),
])
}
// Autogenerate the function "jac" that can be used to compute the Jacobian of f(x) at any point
// x.
get_jacobian!(f, jac);
// nalgebra::DVector
let x0: DVector<f64> = dvector![5.0, 6.0, 7.0];
let jac_eval: DMatrix<f64> = jac(&x0);
// nalgebra::SVector
let x0: SVector<f64, 3> = SVector::from_slice(&[5.0, 6.0, 7.0]);
let jac_eval: DMatrix<f64> = jac(&x0);
// ndarray::Array1
let x0: Array1<f64> = array![5.0, 6.0, 7.0];
let jac_eval: Array2<f64> = jac(&x0);
// faer::Mat
let x0: Mat<f64> = Mat::from_slice(&[5.0, 6.0, 7.0]);
let jac_eval: Mat<f64> = jac(&x0);