3 releases

✓ Uses Rust 2018 edition

0.1.2 Oct 14, 2019
0.1.1 Jan 11, 2019
0.1.0 Jan 11, 2019

#76 in Math

Download history 7/week @ 2019-06-29 13/week @ 2019-07-06 13/week @ 2019-07-13 9/week @ 2019-07-20 2/week @ 2019-07-27 26/week @ 2019-08-03 16/week @ 2019-08-10 10/week @ 2019-08-17 14/week @ 2019-08-24 11/week @ 2019-08-31 10/week @ 2019-09-07 7/week @ 2019-09-14 39/week @ 2019-09-21 14/week @ 2019-09-28 25/week @ 2019-10-05

65 downloads per month
Used in 3 crates (2 directly)

MIT/Apache

82KB
2K SLoC

Build Status

finitediff: Finite Differentiation

This crate contains a wide range of methods for the calculation of gradients, Jacobians and Hessians using forward and central differences. The methods have been implemented for input vectors of the type Vec<f64> and ndarray::Array1<f64>.

See the Documentation for details.

License

Licensed under either of

at your option.

Contribution

Unless you explicitly state otherwise, any contribution intentionally submitted for inclusion in the work by you, as defined in the Apache-2.0 license, shall be dual licensed as above, without any additional terms or conditions.


lib.rs:

This crate contains a wide range of methods for the calculation of gradients, Jacobians and Hessians using forward and central differences. The methods have been implemented for input vectors of the type Vec<f64> and ndarray::Array1<f64>. Central differences are more accurate but require more evaluations of the cost function and are therefore computationally more expensive.

References

Jorge Nocedal and Stephen J. Wright (2006). Numerical Optimization. Springer. ISBN 0-387-30303-0.

Usage

Add this to your dependencies section of Cargo.toml:

[dependencies]
finitediff = "0.1.1"

To use the FiniteDiff trait implementations on the ndarray types, please activate the ndarray feature:

[dependencies]
finitediff = { version = "0.1.1", features = ["ndarray"] }

Examples

Calculation of the gradient

This section illustrates the computation of a gradient of a fuction f at a point x of type Vec<f64>. Using forward differences requires n+1 evaluations of the function f while central differences requires 2*n evaluations.

For Vec<f64>

use finitediff::FiniteDiff;

// Define cost function `f(x)`
let f = |x: &Vec<f64>| -> f64 {
    // ...
#     x[0] + x[1].powi(2)
};

// Point at which gradient should be calculated
let x = vec![1.0f64, 1.0];

// Calculate gradient of `f` at `x` using forward differences
let grad_forward = x.forward_diff(&f);

// Calculate gradient of `f` at `x` using central differences
let grad_central = x.central_diff(&f);
#
#  // Desired solution
#  let res = vec![1.0f64, 2.0];
#
#  // Check result
#  for i in 0..2 {
#      assert!((res[i] - grad_forward[i]).abs() < 1e-6);
#      assert!((res[i] - grad_central[i]).abs() < 1e-6);
#  }

For ndarray::Array1<f64>

use ndarray::{array, Array1};
use finitediff::FiniteDiff;

// Define cost function `f(x)`
let f = |x: &Array1<f64>| -> f64 {
    // ...
#     x[0] + x[1].powi(2)
};

// Point at which gradient should be calculated
let x = array![1.0f64, 1.0];

// Calculate gradient of `f` at `x` using forward differences
let grad_forward = x.forward_diff(&f);

// Calculate gradient of `f` at `x` using central differences
let grad_central = x.central_diff(&f);
#
#  // Desired solution
#  let res = vec![1.0f64, 2.0];
#
#  // Check result
#  for i in 0..2 {
#      assert!((res[i] - grad_forward[i]).abs() < 1e-6);
#      assert!((res[i] - grad_central[i]).abs() < 1e-6);
#  }

Calculation of the Jacobian

Note that the same interface is also implemented for ndarray::Array1<f64> (not shown).

Full Jacobian

The calculation of the full Jacobian requires n+1 evaluations of the function f.

use finitediff::FiniteDiff;

let f = |x: &Vec<f64>| -> Vec<f64> {
    // ...
#      vec![
#          2.0 * (x[1].powi(3) - x[0].powi(2)),
#          3.0 * (x[1].powi(3) - x[0].powi(2)) + 2.0 * (x[2].powi(3) - x[1].powi(2)),
#          3.0 * (x[2].powi(3) - x[1].powi(2)) + 2.0 * (x[3].powi(3) - x[2].powi(2)),
#          3.0 * (x[3].powi(3) - x[2].powi(2)) + 2.0 * (x[4].powi(3) - x[3].powi(2)),
#          3.0 * (x[4].powi(3) - x[3].powi(2)) + 2.0 * (x[5].powi(3) - x[4].powi(2)),
#          3.0 * (x[5].powi(3) - x[4].powi(2)),
#      ]
};

let x = vec![1.0f64, 1.0, 1.0, 1.0, 1.0, 1.0];

// Using forward differences
let jacobian_forward = x.forward_jacobian(&f);

// Using central differences
let jacobian_central = x.central_jacobian(&f);

#  let res = vec![
#      vec![-4.0, -6.0, 0.0, 0.0, 0.0, 0.0],
#      vec![6.0, 5.0, -6.0, 0.0, 0.0, 0.0],
#      vec![0.0, 6.0, 5.0, -6.0, 0.0, 0.0],
#      vec![0.0, 0.0, 6.0, 5.0, -6.0, 0.0],
#      vec![0.0, 0.0, 0.0, 6.0, 5.0, -6.0],
#      vec![0.0, 0.0, 0.0, 0.0, 6.0, 9.0],
#  ];
#
#  // Check result
#  for i in 0..6 {
#      for j in 0..6 {
#          assert!((res[i][j] - jacobian_forward[i][j]).abs() < 1e-6);
#          assert!((res[i][j] - jacobian_central[i][j]).abs() < 1e-6);
#      }
#  }

Product of the Jacobian J(x) with a vector p

Directly computing J(x)*p can be much more efficient than computing J(x) first and then multiplying it with p. While computing the full Jacobian J(x) requires n+1 evaluations of f, J(x)*p only requires 2.

use finitediff::FiniteDiff;

let f = |x: &Vec<f64>| -> Vec<f64> {
    // ...
#      vec![
#          2.0 * (x[1].powi(3) - x[0].powi(2)),
#          3.0 * (x[1].powi(3) - x[0].powi(2)) + 2.0 * (x[2].powi(3) - x[1].powi(2)),
#          3.0 * (x[2].powi(3) - x[1].powi(2)) + 2.0 * (x[3].powi(3) - x[2].powi(2)),
#          3.0 * (x[3].powi(3) - x[2].powi(2)) + 2.0 * (x[4].powi(3) - x[3].powi(2)),
#          3.0 * (x[4].powi(3) - x[3].powi(2)) + 2.0 * (x[5].powi(3) - x[4].powi(2)),
#          3.0 * (x[5].powi(3) - x[4].powi(2)),
#      ]
};

let x = vec![1.0f64, 1.0, 1.0, 1.0, 1.0, 1.0];
let p = vec![1.0f64, 2.0, 3.0, 4.0, 5.0, 6.0];

// using forward differences
let jacobian_forward = x.forward_jacobian_vec_prod(&f, &p);

// using central differences
let jacobian_central = x.central_jacobian_vec_prod(&f, &p);
#
#  let res = vec![8.0, 22.0, 27.0, 32.0, 37.0, 24.0];
#
#  // Check result
#  for i in 0..6 {
#      assert!((res[i] - jacobian_forward[i]).abs() < 11.0*1e-6);
#      assert!((res[i] - jacobian_central[i]).abs() < 11.0*1e-6);
#  }

Sparse Jacobian

If the Jacobian is sparse its structure can be exploited using perturbation vectors. See Nocedal & Wright for details.

use finitediff::{FiniteDiff, PerturbationVector};

let f = |x: &Vec<f64>| -> Vec<f64> {
    // ...
#      vec![
#          2.0 * (x[1].powi(3) - x[0].powi(2)),
#          3.0 * (x[1].powi(3) - x[0].powi(2)) + 2.0 * (x[2].powi(3) - x[1].powi(2)),
#          3.0 * (x[2].powi(3) - x[1].powi(2)) + 2.0 * (x[3].powi(3) - x[2].powi(2)),
#          3.0 * (x[3].powi(3) - x[2].powi(2)) + 2.0 * (x[4].powi(3) - x[3].powi(2)),
#          3.0 * (x[4].powi(3) - x[3].powi(2)) + 2.0 * (x[5].powi(3) - x[4].powi(2)),
#          3.0 * (x[5].powi(3) - x[4].powi(2)),
#      ]
};

let x = vec![1.0f64, 1.0, 1.0, 1.0, 1.0, 1.0];

let pert = vec![
    PerturbationVector::new()
        .add(0, vec![0, 1])
        .add(3, vec![2, 3, 4]),
    PerturbationVector::new()
        .add(1, vec![0, 1, 2])
        .add(4, vec![3, 4, 5]),
    PerturbationVector::new()
        .add(2, vec![1, 2, 3])
        .add(5, vec![4, 5]),
];

// using forward differences
let jacobian_forward = x.forward_jacobian_pert(&f, &pert);

// using central differences
let jacobian_central = x.central_jacobian_pert(&f, &pert);
#
#  let res = vec![
#      vec![-4.0, -6.0, 0.0, 0.0, 0.0, 0.0],
#      vec![6.0, 5.0, -6.0, 0.0, 0.0, 0.0],
#      vec![0.0, 6.0, 5.0, -6.0, 0.0, 0.0],
#      vec![0.0, 0.0, 6.0, 5.0, -6.0, 0.0],
#      vec![0.0, 0.0, 0.0, 6.0, 5.0, -6.0],
#      vec![0.0, 0.0, 0.0, 0.0, 6.0, 9.0],
#  ];
#
#  // Check result
#  for i in 0..6 {
#      for j in 0..6 {
#          assert!((res[i][j] - jacobian_forward[i][j]).abs() < 1e-6);
#          assert!((res[i][j] - jacobian_central[i][j]).abs() < 1e-6);
#      }
#  }

Calculation of the Hessian

Note that the same interface is also implemented for ndarray::Array1<f64> (not shown).

Full Hessian

use finitediff::FiniteDiff;

let g = |x: &Vec<f64>| -> Vec<f64> {
    // ...
#     vec![1.0, 2.0 * x[1], x[3].powi(2), 2.0 * x[3] * x[2]]
};

let x = vec![1.0f64, 1.0, 1.0, 1.0];

// using forward differences
let hessian_forward = x.forward_hessian(&g);

// using central differences
let hessian_central = x.central_hessian(&g);
#
#  let res = vec![
#      vec![0.0, 0.0, 0.0, 0.0],
#      vec![0.0, 2.0, 0.0, 0.0],
#      vec![0.0, 0.0, 0.0, 2.0],
#      vec![0.0, 0.0, 2.0, 2.0],
#  ];
#
#  // Check result
#  for i in 0..4 {
#      for j in 0..4 {
#          assert!((res[i][j] - hessian_forward[i][j]).abs() < 1e-6);
#          assert!((res[i][j] - hessian_central[i][j]).abs() < 1e-6);
#      }
#  }

Product of the Hessian H(x) with a vector p

use finitediff::FiniteDiff;

let g = |x: &Vec<f64>| -> Vec<f64> {
    // ...
#     vec![1.0, 2.0 * x[1], x[3].powi(2), 2.0 * x[3] * x[2]]
};

let x = vec![1.0f64, 1.0, 1.0, 1.0];
let p = vec![2.0, 3.0, 4.0, 5.0];

// using forward differences
let hessian_forward = x.forward_hessian_vec_prod(&g, &p);

// using forward differences
let hessian_central = x.central_hessian_vec_prod(&g, &p);
#
#  let res = vec![0.0, 6.0, 10.0, 18.0];
#
#  for i in 0..4 {
#      assert!((res[i] - hessian_forward[i]).abs() < 1e-6);
#      assert!((res[i] - hessian_central[i]).abs() < 1e-6);
#  }

Calculation of the Hessian without knowledge of the gradient

use finitediff::FiniteDiff;

let f = |x: &Vec<f64>| -> f64 {
    // ...
#     x[0] + x[1].powi(2) + x[2] * x[3].powi(2)
};

let x = vec![1.0f64, 1.0, 1.0, 1.0];

let hessian = x.forward_hessian_nograd(&f);
#
#  let res = vec![
#      vec![0.0, 0.0, 0.0, 0.0],
#      vec![0.0, 2.0, 0.0, 0.0],
#      vec![0.0, 0.0, 0.0, 2.0],
#      vec![0.0, 0.0, 2.0, 2.0],
#  ];
#
#  // Check result
#  for i in 0..4 {
#      for j in 0..4 {
#          assert!((res[i][j] - hessian[i][j]).abs() < 1e-6)
#      }
#  }

Calculation of the sparse Hessian without knowledge of the gradient

use finitediff::FiniteDiff;

let f = |x: &Vec<f64>| -> f64 {
    // ...
#     x[0] + x[1].powi(2) + x[2] * x[3].powi(2)
};

let x = vec![1.0f64, 1.0, 1.0, 1.0];

// Indices at which the Hessian should be evaluated. All other
// elements of the Hessian will be zero
let indices = vec![[1, 1], [2, 3], [3, 3]];

let hessian = x.forward_hessian_nograd_sparse(&f, indices);
#
#  let res = vec![
#      vec![0.0, 0.0, 0.0, 0.0],
#      vec![0.0, 2.0, 0.0, 0.0],
#      vec![0.0, 0.0, 0.0, 2.0],
#      vec![0.0, 0.0, 2.0, 2.0],
#  ];
#
#  // Check result
#  for i in 0..4 {
#      for j in 0..4 {
#          assert!((res[i][j] - hessian[i][j]).abs() < 1e-6)
#      }
#  }

Dependencies