Skip to content

Commit

Permalink
use h=eps.cbrt() instead of sqrt for central differences
Browse files Browse the repository at this point in the history
  • Loading branch information
@stefan-k
stefan-k committed Mar 4, 2024
1 parent fa9f6b4 commit d03f507
Show file tree
Hide file tree
Showing 6 changed files with 48 additions and 48 deletions.
8 changes: 4 additions & 4 deletions crates/finitediff/src/ndarray_m/diff.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,14 +35,14 @@ pub fn central_diff_ndarray<F>(
where
F: Float + FromPrimitive,
{
let eps_sqrt = F::epsilon().sqrt();
let eps_cbrt = F::epsilon().cbrt();

let mut xt = x.clone();
(0..x.len())
.map(|i| {
let fx1 = mod_and_calc(&mut xt, f, i, eps_sqrt);
let fx2 = mod_and_calc(&mut xt, f, i, -eps_sqrt);
(fx1 - fx2) / (F::from_f64(2.0).unwrap() * eps_sqrt)
let fx1 = mod_and_calc(&mut xt, f, i, eps_cbrt);
let fx2 = mod_and_calc(&mut xt, f, i, -eps_cbrt);
(fx1 - fx2) / (F::from_f64(2.0).unwrap() * eps_cbrt)
})
.collect()
}
Expand Down
16 changes: 8 additions & 8 deletions crates/finitediff/src/ndarray_m/hessian.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@ pub fn central_hessian_ndarray<F>(
where
F: Float + FromPrimitive,
{
let eps_sqrt = F::epsilon().sqrt();
let eps_cbrt = F::epsilon().cbrt();

let mut xt = x.clone();
// TODO: get rid of this!
Expand All @@ -53,10 +53,10 @@ where
let n = x.len();
let mut out = ndarray::Array2::zeros((n, rn));
for i in 0..n {
let fx1 = mod_and_calc(&mut xt, grad, i, eps_sqrt);
let fx2 = mod_and_calc(&mut xt, grad, i, -eps_sqrt);
let fx1 = mod_and_calc(&mut xt, grad, i, eps_cbrt);
let fx2 = mod_and_calc(&mut xt, grad, i, -eps_cbrt);
for j in 0..rn {
out[(i, j)] = (fx1[j] - fx2[j]) / (F::from_f64(2.0).unwrap() * eps_sqrt);
out[(i, j)] = (fx1[j] - fx2[j]) / (F::from_f64(2.0).unwrap() * eps_cbrt);
}
}
// restore symmetry
Expand Down Expand Up @@ -87,13 +87,13 @@ pub fn central_hessian_vec_prod_ndarray<F>(
where
F: Float + FromPrimitive + ScalarOperand,
{
let eps_sqrt = F::epsilon().sqrt();
let eps_cbrt = F::epsilon().cbrt();

let x1 = x + &(p.mapv(|pi| pi * eps_sqrt));
let x2 = x - &(p.mapv(|pi| pi * eps_sqrt));
let x1 = x + &(p.mapv(|pi| pi * eps_cbrt));
let x2 = x - &(p.mapv(|pi| pi * eps_cbrt));
let fx1 = (grad)(&x1);
let fx2 = (grad)(&x2);
(fx1 - fx2) / (F::from_f64(2.0).unwrap() * eps_sqrt)
(fx1 - fx2) / (F::from_f64(2.0).unwrap() * eps_cbrt)
}

pub fn forward_hessian_nograd_ndarray<F>(
Expand Down
24 changes: 12 additions & 12 deletions crates/finitediff/src/ndarray_m/jacobian.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ pub fn central_jacobian_ndarray<F>(
where
F: Float + FromPrimitive,
{
let eps_sqrt = F::epsilon().sqrt();
let eps_cbrt = F::epsilon().cbrt();

let mut xt = x.clone();

Expand All @@ -60,10 +60,10 @@ where

let mut out = Array2::zeros((n, rn));
for i in 0..n {
let fx1 = mod_and_calc(&mut xt, fs, i, eps_sqrt);
let fx2 = mod_and_calc(&mut xt, fs, i, -eps_sqrt);
let fx1 = mod_and_calc(&mut xt, fs, i, eps_cbrt);
let fx2 = mod_and_calc(&mut xt, fs, i, -eps_cbrt);
for j in 0..rn {
out[(i, j)] = (fx1[j] - fx2[j]) / (F::from_f64(2.0).unwrap() * eps_sqrt);
out[(i, j)] = (fx1[j] - fx2[j]) / (F::from_f64(2.0).unwrap() * eps_cbrt);
}
}
out
Expand Down Expand Up @@ -93,12 +93,12 @@ pub fn central_jacobian_vec_prod_ndarray<F>(
where
F: Float + FromPrimitive + ScalarOperand,
{
let eps_sqrt = F::epsilon().sqrt();
let x1 = x + &p.mapv(|pi| eps_sqrt * pi);
let x2 = x + &p.mapv(|pi| -eps_sqrt * pi);
let eps_cbrt = F::epsilon().sqrt();
let x1 = x + &p.mapv(|pi| eps_cbrt * pi);
let x2 = x + &p.mapv(|pi| -eps_cbrt * pi);
let fx1 = (fs)(&x1);
let fx2 = (fs)(&x2);
(fx1 - fx2) / (F::from_f64(2.0).unwrap() * eps_sqrt)
(fx1 - fx2) / (F::from_f64(2.0).unwrap() * eps_cbrt)
}

pub fn forward_jacobian_pert_ndarray<F>(
Expand Down Expand Up @@ -142,19 +142,19 @@ pub fn central_jacobian_pert_ndarray<F>(
where
F: Float + FromPrimitive + AddAssign,
{
let eps_sqrt = F::epsilon().sqrt();
let eps_cbrt = F::epsilon().cbrt();

let mut out = ndarray::Array2::zeros((1, 1));
let mut xt = x.clone();
for (i, pert_item) in pert.iter().enumerate() {
for j in pert_item.x_idx.iter() {
xt[*j] += eps_sqrt;
xt[*j] += eps_cbrt;
}

let fx1 = (fs)(&xt);

for j in pert_item.x_idx.iter() {
xt[*j] = x[*j] - eps_sqrt;
xt[*j] = x[*j] - eps_cbrt;
}

let fx2 = (fs)(&xt);
Expand All @@ -169,7 +169,7 @@ where

for (k, x_idx) in pert_item.x_idx.iter().enumerate() {
for j in pert_item.r_idx[k].iter() {
out[(*x_idx, *j)] = (fx1[*j] - fx2[*j]) / (F::from_f64(2.0).unwrap() * eps_sqrt);
out[(*x_idx, *j)] = (fx1[*j] - fx2[*j]) / (F::from_f64(2.0).unwrap() * eps_cbrt);
}
}
}
Expand Down
8 changes: 4 additions & 4 deletions crates/finitediff/src/vec/diff.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,12 +30,12 @@ where
F: Float + FromPrimitive,
{
let mut xt = x.to_owned();
let eps_sqrt = F::epsilon().sqrt();
let eps_cbrt = F::epsilon().cbrt();
(0..x.len())
.map(|i| {
let fx1 = mod_and_calc(&mut xt, f, i, eps_sqrt);
let fx2 = mod_and_calc(&mut xt, f, i, -eps_sqrt);
(fx1 - fx2) / (F::from_f64(2.0).unwrap() * eps_sqrt)
let fx1 = mod_and_calc(&mut xt, f, i, eps_cbrt);
let fx2 = mod_and_calc(&mut xt, f, i, -eps_cbrt);
(fx1 - fx2) / (F::from_f64(2.0).unwrap() * eps_cbrt)
})
.collect()
}
Expand Down
16 changes: 8 additions & 8 deletions crates/finitediff/src/vec/hessian.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,15 +36,15 @@ pub fn central_hessian_vec<F>(x: &[F], grad: &dyn Fn(&Vec<F>) -> Vec<F>) -> Vec<
where
F: Float + FromPrimitive,
{
let eps_sqrt = F::epsilon().sqrt();
let eps_cbrt = F::epsilon().cbrt();
let mut xt = x.to_owned();
let out: Vec<Vec<F>> = (0..x.len())
.map(|i| {
let fx1 = mod_and_calc(&mut xt, grad, i, eps_sqrt);
let fx2 = mod_and_calc(&mut xt, grad, i, -eps_sqrt);
let fx1 = mod_and_calc(&mut xt, grad, i, eps_cbrt);
let fx2 = mod_and_calc(&mut xt, grad, i, -eps_cbrt);
fx1.iter()
.zip(fx2.iter())
.map(|(&a, &b)| (a - b) / (F::from_f64(2.0).unwrap() * eps_sqrt))
.map(|(&a, &b)| (a - b) / (F::from_f64(2.0).unwrap() * eps_cbrt))
.collect::<Vec<F>>()
})
.collect();
Expand Down Expand Up @@ -82,23 +82,23 @@ pub fn central_hessian_vec_prod_vec<F>(x: &[F], grad: &dyn Fn(&Vec<F>) -> Vec<F>
where
F: Float + FromPrimitive,
{
let eps_sqrt = F::epsilon().sqrt();
let eps_cbrt = F::epsilon().cbrt();
let out: Vec<F> = {
let x1 = x
.iter()
.zip(p.iter())
.map(|(&xi, &pi)| xi + pi * eps_sqrt)
.map(|(&xi, &pi)| xi + pi * eps_cbrt)
.collect();
let x2 = x
.iter()
.zip(p.iter())
.map(|(&xi, &pi)| xi - pi * eps_sqrt)
.map(|(&xi, &pi)| xi - pi * eps_cbrt)
.collect();
let fx1 = (grad)(&x1);
let fx2 = (grad)(&x2);
fx1.iter()
.zip(fx2.iter())
.map(|(&a, &b)| (a - b) / (F::from_f64(2.0).unwrap() * eps_sqrt))
.map(|(&a, &b)| (a - b) / (F::from_f64(2.0).unwrap() * eps_cbrt))
.collect::<Vec<F>>()
};
out
Expand Down
24 changes: 12 additions & 12 deletions crates/finitediff/src/vec/jacobian.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,14 +35,14 @@ where
F: Float + FromPrimitive,
{
let mut xt = x.to_owned();
let eps_sqrt = F::epsilon().sqrt();
let eps_cbrt = F::epsilon().cbrt();
(0..x.len())
.map(|i| {
let fx1 = mod_and_calc(&mut xt, fs, i, eps_sqrt);
let fx2 = mod_and_calc(&mut xt, fs, i, -eps_sqrt);
let fx1 = mod_and_calc(&mut xt, fs, i, eps_cbrt);
let fx2 = mod_and_calc(&mut xt, fs, i, -eps_cbrt);
fx1.iter()
.zip(fx2.iter())
.map(|(&a, &b)| (a - b) / (F::from_f64(2.0).unwrap() * eps_sqrt))
.map(|(&a, &b)| (a - b) / (F::from_f64(2.0).unwrap() * eps_cbrt))
.collect::<Vec<F>>()
})
.collect()
Expand Down Expand Up @@ -74,22 +74,22 @@ pub fn central_jacobian_vec_prod_vec<F>(x: &[F], fs: &dyn Fn(&Vec<F>) -> Vec<F>,
where
F: Float + FromPrimitive,
{
let eps_sqrt = F::epsilon().sqrt();
let eps_cbrt = F::epsilon().cbrt();
let x1 = x
.iter()
.zip(p.iter())
.map(|(&xi, &pi)| xi + eps_sqrt * pi)
.map(|(&xi, &pi)| xi + eps_cbrt * pi)
.collect();
let x2 = x
.iter()
.zip(p.iter())
.map(|(&xi, &pi)| xi - eps_sqrt * pi)
.map(|(&xi, &pi)| xi - eps_cbrt * pi)
.collect();
let fx1 = (fs)(&x1);
let fx2 = (fs)(&x2);
fx1.iter()
.zip(fx2.iter())
.map(|(&a, &b)| (a - b) / (F::from_f64(2.0).unwrap() * eps_sqrt))
.map(|(&a, &b)| (a - b) / (F::from_f64(2.0).unwrap() * eps_cbrt))
.collect::<Vec<F>>()
}

Expand Down Expand Up @@ -134,17 +134,17 @@ where
F: Float + FromPrimitive + AddAssign,
{
let mut out = vec![];
let eps_sqrt = F::epsilon().sqrt();
let eps_cbrt = F::epsilon().cbrt();
let mut xt = x.to_owned();
for (i, pert_item) in pert.iter().enumerate() {
for j in pert_item.x_idx.iter() {
xt[*j] += eps_sqrt;
xt[*j] += eps_cbrt;
}

let fx1 = (fs)(&xt);

for j in pert_item.x_idx.iter() {
xt[*j] = x[*j] - eps_sqrt;
xt[*j] = x[*j] - eps_cbrt;
}

let fx2 = (fs)(&xt);
Expand All @@ -159,7 +159,7 @@ where

for (k, x_idx) in pert_item.x_idx.iter().enumerate() {
for j in pert_item.r_idx[k].iter() {
out[*x_idx][*j] = (fx1[*j] - fx2[*j]) / (F::from_f64(2.0).unwrap() * eps_sqrt);
out[*x_idx][*j] = (fx1[*j] - fx2[*j]) / (F::from_f64(2.0).unwrap() * eps_cbrt);
}
}
}
Expand Down

0 comments on commit d03f507

Please sign in to comment.