Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Added derivative and Hessian for Goldstein-Price test function #440

Merged
merged 1 commit into from
Feb 9, 2024
Merged

Added derivative and Hessian for Goldstein-Price test function #440

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
278 changes: 276 additions & 2 deletions crates/argmin-testfunctions/src/goldsteinprice.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@
//!
//! Defined as
//!
//! `f(x_1, x_2) = [1 + (x_1 + x_2 + 1)^2 * (19 - 14*x_2 + + 3*x_1^2 - 14*x_2 6*x_1*x_2 + 3*x_2^2)]
//! `f(x_1, x_2) = [1 + (x_1 + x_2 + 1)^2 * (19 - 14*x_1 + 3*x_1^2 - 14*x_2 + 6*x_1*x_2 + 3*x_2^2)]
//! * [30 + (2*x_1 - 3*x_2)^2(18 - 32 * x_1 + 12* x_1^2 + 48 * x_2 -
//! 36 * x_1 * x_2 + 27 * x_2^2) ]`
//!
Expand All @@ -23,7 +23,7 @@ use num::{Float, FromPrimitive};
///
/// Defined as
///
/// `f(x_1, x_2) = [1 + (x_1 + x_2 + 1)^2 * (19 - 14*x_2 + 3*x_1^2 - 14*x_2 6*x_1*x_2 + 3*x_2^2)]
/// `f(x_1, x_2) = [1 + (x_1 + x_2 + 1)^2 * (19 - 14*x_1 + 3*x_1^2 - 14*x_2 + 6*x_1*x_2 + 3*x_2^2)]
/// * [30 + (2*x_1 - 3*x_2)^2(18 - 32 * x_1 + 12* x_1^2 + 48 * x_2 -
/// 36 * x_1 * x_2 + 27 * x_2^2) ]`
///
Expand Down Expand Up @@ -55,10 +55,203 @@ where
* (n18 - n32 * x1 + n12 * x1.powi(2) + n48 * x2 - n36 * x1 * x2 + n27 * x2.powi(2)))
}

/// Derivative of Goldstein-Price test function
pub fn goldsteinprice_derivative<T>(param: &[T; 2]) -> [T; 2]
where
T: Float + FromPrimitive,
{
let [x1, x2] = *param;

let n1 = T::from_f64(1.0).unwrap();
let n2 = T::from_f64(2.0).unwrap();
let n3 = T::from_f64(3.0).unwrap();
let n4 = T::from_f64(4.0).unwrap();
let n6 = T::from_f64(6.0).unwrap();
let n12 = T::from_f64(12.0).unwrap();
let n14 = T::from_f64(14.0).unwrap();
let n18 = T::from_f64(18.0).unwrap();
let n19 = T::from_f64(19.0).unwrap();
let n24 = T::from_f64(24.0).unwrap();
let n27 = T::from_f64(27.0).unwrap();
let n30 = T::from_f64(30.0).unwrap();
let n32 = T::from_f64(32.0).unwrap();
let n36 = T::from_f64(36.0).unwrap();
let n48 = T::from_f64(48.0).unwrap();
let n54 = T::from_f64(54.0).unwrap();

let x1s = x1.powi(2);
let x2s = x2.powi(2);

[
(n2 * (x1 + x2 + n1) * (n3 * x1s + n6 * x2 * x1 - n14 * x1 + n3 * x2s - n14 * x2 + n19)
+ (x1 + x2 + n1).powi(2) * (n6 * x1 + n6 * x2 - n14))
* ((n2 * x1 - n3 * x2).powi(2)
* (n12 * x1s - n36 * x2 * x1 - n32 * x1 + n27 * x2s + n48 * x2 + n18)
+ n30)
+ ((x1 + x2 + n1).powi(2)
* (n3 * x1s + n6 * x2 * x1 - n14 * x1 + n3 * x2s - n14 * x2 + n19)
+ n1)
* (n4
* (n2 * x1 - n3 * x2)
* (n12 * x1s - n36 * x2 * x1 - n32 * x1 + n27 * x2s + n48 * x2 + n18)
+ (n2 * x1 - n3 * x2).powi(2) * (n24 * x1 - n36 * x2 - n32)),
((x2 + x1 + n1).powi(2) * (n3 * x2s + n6 * x1 * x2 - n14 * x2 + n3 * x1s - n14 * x1 + n19)
+ n1)
* ((n2 * x1 - n3 * x2).powi(2) * (n54 * x2 - n36 * x1 + n48)
- n6 * (n2 * x1 - n3 * x2)
* (n27 * x2s - n36 * x1 * x2 + n48 * x2 + n12 * x1s - n32 * x1 + n18))
+ (n2
* (x2 + x1 + n1)
* (n3 * x2s + n6 * x1 * x2 - n14 * x2 + n3 * x1s - n14 * x1 + n19)
+ (x2 + x1 + n1).powi(2) * (n6 * x2 + n6 * x1 - n14))
* ((n2 * x1 - n3 * x2).powi(2)
* (n27 * x2s - n36 * x1 * x2 + n48 * x2 + n12 * x1s - n32 * x1 + n18)
+ n30),
]
}

/// Hessian of Goldstein-Price test function
pub fn goldsteinprice_hessian<T>(param: &[T; 2]) -> [[T; 2]; 2]
where
T: Float + FromPrimitive,
{
let [x1, x2] = *param;

let n840 = T::from_f64(840.0).unwrap();
let n1296 = T::from_f64(1296.0).unwrap();
let n2016 = T::from_f64(2016.0).unwrap();
let n2520 = T::from_f64(2520.0).unwrap();
let n2916 = T::from_f64(2916.0).unwrap();
let n3360 = T::from_f64(3360.0).unwrap();
let n4680 = T::from_f64(4680.0).unwrap();
let n5184 = T::from_f64(5184.0).unwrap();
let n5940 = T::from_f64(5940.0).unwrap();
let n6120 = T::from_f64(6120.0).unwrap();
let n6432 = T::from_f64(6432.0).unwrap();
let n6804 = T::from_f64(6804.0).unwrap();
let n7344 = T::from_f64(7344.0).unwrap();
let n7440 = T::from_f64(7440.0).unwrap();
let n7776 = T::from_f64(7776.0).unwrap();
let n8064 = T::from_f64(8064.0).unwrap();
let n10080 = T::from_f64(10080.0).unwrap();
let n10740 = T::from_f64(10740.0).unwrap();
let n11016 = T::from_f64(11016.0).unwrap();
let n11160 = T::from_f64(11160.0).unwrap();
let n11664 = T::from_f64(11664.0).unwrap();
let n12096 = T::from_f64(12096.0).unwrap();
let n14688 = T::from_f64(14688.0).unwrap();
let n15552 = T::from_f64(15552.0).unwrap();
let n15660 = T::from_f64(15660.0).unwrap();
let n17352 = T::from_f64(17352.0).unwrap();
let n17460 = T::from_f64(17460.0).unwrap();
let n17496 = T::from_f64(17496.0).unwrap();
let n18360 = T::from_f64(18360.0).unwrap();
let n19440 = T::from_f64(19440.0).unwrap();
let n19680 = T::from_f64(19680.0).unwrap();
let n20880 = T::from_f64(20880.0).unwrap();
let n23760 = T::from_f64(23760.0).unwrap();
let n24480 = T::from_f64(24480.0).unwrap();
let n25920 = T::from_f64(25920.0).unwrap();
let n26880 = T::from_f64(26880.0).unwrap();
let n27216 = T::from_f64(27216.0).unwrap();
let n27540 = T::from_f64(27540.0).unwrap();
let n28560 = T::from_f64(28560.0).unwrap();
let n29448 = T::from_f64(29448.0).unwrap();
let n30240 = T::from_f64(30240.0).unwrap();
let n30720 = T::from_f64(30720.0).unwrap();
let n31104 = T::from_f64(31104.0).unwrap();
let n32256 = T::from_f64(32256.0).unwrap();
let n34704 = T::from_f64(34704.0).unwrap();
let n36720 = T::from_f64(36720.0).unwrap();
let n38592 = T::from_f64(38592.0).unwrap();
let n38880 = T::from_f64(38880.0).unwrap();
let n40320 = T::from_f64(40320.0).unwrap();
let n40824 = T::from_f64(40824.0).unwrap();
let n41760 = T::from_f64(41760.0).unwrap();
let n42960 = T::from_f64(42960.0).unwrap();
let n43740 = T::from_f64(43740.0).unwrap();
let n47520 = T::from_f64(47520.0).unwrap();
let n48960 = T::from_f64(48960.0).unwrap();
let n51840 = T::from_f64(51840.0).unwrap();
let n58320 = T::from_f64(58320.0).unwrap();
let n59040 = T::from_f64(59040.0).unwrap();
let n64440 = T::from_f64(64440.0).unwrap();
let n69840 = T::from_f64(69840.0).unwrap();
let n70848 = T::from_f64(70848.0).unwrap();
let n73440 = T::from_f64(73440.0).unwrap();
let n73728 = T::from_f64(73728.0).unwrap();
let n92160 = T::from_f64(92160.0).unwrap();
let n104760 = T::from_f64(104760.0).unwrap();
let n132840 = T::from_f64(132840.0).unwrap();
let n142560 = T::from_f64(142560.0).unwrap();
let n141696 = T::from_f64(141696.0).unwrap();
let n172152 = T::from_f64(172152.0).unwrap();

let x1p2 = x1.powi(2);
let x1p3 = x1.powi(3);
let x1p4 = x1.powi(4);
let x1p5 = x1.powi(5);
let x1p6 = x1.powi(6);
let x2p2 = x2.powi(2);
let x2p3 = x2.powi(3);
let x2p4 = x2.powi(4);
let x2p5 = x2.powi(5);
let x2p6 = x2.powi(6);

let a = n8064 * x1p6
+ (-n12096 * x2 - n32256) * x1p5
+ (-n19440 * x2p2 + n40320 * x2 + n28560) * x1p4
+ (n24480 * x2p3 + n51840 * x2p2 - n3360 * x2 + n26880) * x1p3
+ (n15660 * x2p4 - n48960 * x2p3 - n64440 * x2p2 - n92160 * x2 - n29448) * x1p2
+ (-n11016 * x2p5 - n20880 * x2p4 + n7440 * x2p3 + n59040 * x2p2 + n34704 * x2 - n6432)
* x1
- n2916 * x2p6
+ n7344 * x2p5
+ n17460 * x2p4
+ n10080 * x2p3
+ n15552 * x2p2
+ n14688 * x2
+ n2520;

let b = n40824 * x2p6
+ (n40824 * x1 - n27216) * x2p5
+ (-n43740 * x1p2 + n58320 * x1 - n132840) * x2p4
+ (-n36720 * x1p3 + n73440 * x1p2 - n23760 * x1 + n38880) * x2p3
+ (n15660 * x1p4 - n41760 * x1p3 + n104760 * x1p2 - n142560 * x1 + n172152) * x2p2
+ (n7344 * x1p5 - n24480 * x1p4 + n7440 * x1p3 + n30240 * x1p2 - n141696 * x1 + n73728)
* x2
- n1296 * x1p6
+ n5184 * x1p5
- n10740 * x1p4
+ n19680 * x1p3
+ n15552 * x1p2
- n38592 * x1
+ n6120;

let offdiag = n6804 * x2p6
+ (n11664 - n17496 * x1) * x2p5
+ (-n27540 * x1p2 + n36720 * x1 - n5940) * x2p4
+ (n20880 * x1p3 - n41760 * x1p2 + n69840 * x1 - n47520) * x2p3
+ (n18360 * x1p4 - n48960 * x1p3 + n11160 * x1p2 + n30240 * x1 - n70848) * x2p2
+ (-n7776 * x1p5 + n25920 * x1p4 - n42960 * x1p3 + n59040 * x1p2 + n31104 * x1 - n38592)
* x2
- n2016 * x1p6
+ n8064 * x1p5
- n840 * x1p4
- n30720 * x1p3
+ n17352 * x1p2
+ n14688 * x1
- n4680;

[[a, offdiag], [offdiag, b]]
}

#[cfg(test)]
mod tests {
use super::*;
use approx::assert_relative_eq;
use finitediff::FiniteDiff;
use proptest::prelude::*;
use std::{f32, f64};

#[test]
Expand All @@ -73,5 +266,86 @@ mod tests {
3_f64,
epsilon = f64::EPSILON
);

let deriv = goldsteinprice_derivative(&[0.0_f64, -1.0_f64]);
for i in 0..2 {
assert_relative_eq!(deriv[i], 0.0, epsilon = f64::EPSILON);
}
}

proptest! {
#[test]
fn test_goldsteinprice_derivative_finitediff(a in -2.0..2.0, b in -2.0..2.0) {
let param = [a, b];
let derivative = goldsteinprice_derivative(&param);
let derivative_fd = Vec::from(param).central_diff(&|x| goldsteinprice(&[x[0], x[1]]));
println!("1: {derivative:?} at {a}/{b}");
println!("2: {derivative_fd:?} at {a}/{b}");
for i in 0..derivative.len() {
assert_relative_eq!(derivative[i], derivative_fd[i], epsilon = 1e-1);
}
}
}

proptest! {
#[test]
fn test_goldsteinprice_derivative_finitediff_narrow(a in -0.5..0.5, b in -0.5..0.5) {
// This evaluates the function on a narrower domain, which allows us to have a lower
// epsilon, as the function is pretty steep at the boundary, which isn't great for
// accuracy when using finite differentiation.
let param = [a, b];
let derivative = goldsteinprice_derivative(&param);
let derivative_fd = Vec::from(param).central_diff(&|x| goldsteinprice(&[x[0], x[1]]));
println!("1: {derivative:?} at {a}/{b}");
println!("2: {derivative_fd:?} at {a}/{b}");
for i in 0..derivative.len() {
assert_relative_eq!(derivative[i], derivative_fd[i], epsilon = 1e-3);
}
}
}

proptest! {
#[test]
fn test_goldsteinprice_hessian_finitediff(a in -2.0..2.0, b in -2.0..2.0) {
let param = [a, b];
let hessian = goldsteinprice_hessian(&param);
let hessian_fd =
Vec::from(param).central_hessian(&|x| goldsteinprice_derivative(&[x[0], x[1]]).to_vec());
let n = hessian.len();
println!("1: {hessian:?} at {a}/{b}");
println!("2: {hessian_fd:?} at {a}/{b}");
for i in 0..n {
assert_eq!(hessian[i].len(), n);
for j in 0..n {
if hessian_fd[i][j].is_finite() {
assert_relative_eq!(hessian[i][j], hessian_fd[i][j], epsilon = 1e-1);
}
}
}
}
}

proptest! {
#[test]
fn test_goldsteinprice_hessian_finitediff_narrow(a in -0.5..0.5, b in -0.5..0.5) {
// This evaluates the function on a narrower domain, which allows us to have a lower
// epsilon, as the function is pretty steep at the boundary, which isn't great for
// accuracy when using finite differentiation.
let param = [a, b];
let hessian = goldsteinprice_hessian(&param);
let hessian_fd =
Vec::from(param).central_hessian(&|x| goldsteinprice_derivative(&[x[0], x[1]]).to_vec());
let n = hessian.len();
println!("1: {hessian:?} at {a}/{b}");
println!("2: {hessian_fd:?} at {a}/{b}");
for i in 0..n {
assert_eq!(hessian[i].len(), n);
for j in 0..n {
if hessian_fd[i][j].is_finite() {
assert_relative_eq!(hessian[i][j], hessian_fd[i][j], epsilon = 1e-3);
}
}
}
}
}
}
Loading