Bayesian Optimization

Cargo.toml

@@ -20,6 +20,8 @@ num-integer = "0.1"
 opensrdk-linear-algebra = "0.8.4"
 opensrdk-kernel-method = "0.2.0"
 special = "0.8.1"
+optimize = "0.1.0"
+ndarray = "0.15.6"
 
 [dev-dependencies]
 blas-src = { version = "0.8", features = ["intel-mkl"] } 

src/nonparametric/elliptical_process/bayesian_optimization/expected_improvement.rs

@@ -0,0 +1,43 @@
+use super::AcquisitionFunctions;
+// use crate::{Normal, NormalParams};
+
+pub struct ExpectedImprovement {
+    f_vec: Vec<f64>,
+}
+
+impl AcquisitionFunctions for ExpectedImprovement {
+    fn value(&self, theta: &crate::NormalParams) -> f64 {
+        let mu = theta.mu();
+        let sigma = theta.sigma();
+        let tau = self.f_vec.iter().max().unwrap();
+        let t = (mu - tau) / sigma;
+        // let n = Normal;
+        // let phi_large = n.p_kernel(n, t, &NormalParams::new(0.0, 1.0).unwrap());
+        let phi_large = pdf(t);
+        let phi_small = cdf(t);
+
+        (mu - tau) * phi_large + sigma * phi_small
+    }
+}
+
+// Abramowitz and Stegun (1964) formula 26.2.17
+// precision: abs(err) < 7.5e-8
+
+fn pdf(x: f64) -> f64 {
+    ((-x * x) / 2.0).exp() / (2.0 * std::f64::consts::PI)
+}
+fn cdf(x: f64) -> f64 {
+    // constants
+    const p: f64 = 0.2316419;
+    const b1: f64 = 0.31938153;
+    const b2: f64 = -0.356563782;
+    const b3: f64 = 1.781477937;
+    const b4: f64 = -1.821255978;
+    const b5: f64 = 1.330274429;
+
+    let t = 1.0 / (1.0 + p * x.abs());
+    let z = (-x * x / 2.0).exp() / (2.0 * std::f64::consts::PI).sqrt();
+    let y = 1.0 - z * ((((b5 * t + b4) * t + b3) * t + b2) * t + b1) * t;
+
+    return if x > 0.0 { y } else { 1.0 - y };
+}

src/nonparametric/elliptical_process/bayesian_optimization/mod.rs

@@ -0,0 +1,120 @@
+pub mod expected_improvement;
+pub mod upper_confidence_bound;
+
+pub use expected_improvement::*;
+use ndarray::{Array, ArrayView1};
+use opensrdk_kernel_method::{Periodic, RBF};
+use rand::rngs::StdRng;
+pub use upper_confidence_bound::*;
+
+use crate::{nonparametric::GaussianProcessRegressor, NormalParams};
+use optimize::NelderMeadBuilder;
+
+use super::BaseEllipticalProcessParams;
+
+pub trait AcquisitionFunctions {
+    fn value(&self, theta: &NormalParams) -> f64;
+}
+
+struct Data {
+    x_data: Vec<f64>,
+    y_data: Vec<f64>,
+}
+
+#[test]
+fn test_main() {
+    let mut n: usize = 0;
+    let mut data = Data {
+        x_data: vec![],
+        y_data: vec![],
+    };
+
+    loop {
+        let mut rng = StdRng::from_seed([1; 32]);
+        let mut x: f64 = rng.gen();
+
+        sampling(&data, &x);
+
+        n += 1;
+
+        if n == 20 {
+            break;
+        }
+    }
+
+    loop {
+        let xs = maximize_ucb(&data, n);
+        // let xs = maximize_ei(&data);
+
+        sampling(&data, &xs);
+
+        n += 1;
+    }
+}
+
+fn objective(x: &f64) -> f64 {
+    x + x ^ 2.0
+}
+
+fn sampling(mut data: &Data, x: &f64) {
+    let y = objective(x);
+    data.x_data.push(x);
+    data.y_data.push(y);
+}
+
+fn gp_regression(x: &Vec<f64>, y: &Vec<f64>, xs: f64) -> NormalParams {
+    let kernel = RBF + Periodic;
+    let theta = vec![1.0; kernel.params_len()];
+    let sigma = 1.0;
+
+    let base_params = BaseEllipticalProcessParams::new(kernel, x, theta, sigma).unwrap();
+    let params_y = base_params.exact(&y).unwrap();
+    let mu = params_y.gp_predict(xs).unwrap().mu();
+    let sigma = params_y.gp_predict(xs).unwrap().sigma();
+
+    [mu, sigma]
+}
+
+fn maximize_ucb(data: &Data, n: usize) -> f64 {
+    let func_to_minimize = |xs: ArrayView1<f64>| {
+        let theta: NormalParams = gp_regression(&data.x_data, &data.y_data, xs);
+        let ucb = UpperConfidenceBound { trial: n };
+        -ucb.value(&theta)
+    };
+
+    let minimizer = NelderMeadBuilder::default()
+        .xtol(1e-6f64)
+        .ftol(1e-6f64)
+        .maxiter(50000)
+        .build()
+        .unwrap();
+
+    // Set the starting guess
+    let args = Array::from_vec(vec![3.0, -8.3]);
+    let xs = minimizer.minimize(&func_to_minimize, args.view());
+
+    xs
+}
+
+fn maximize_ei(data: &Data) -> f64 {
+    let func_to_minimize = |xs: ArrayView1<f64>| {
+        let theta: NormalParams = gp_regression(&data.x_data, &data.y_data, xs);
+        let ei = ExpectedImprovement {
+            f_vec: &data.y_data,
+        };
+        -ei.value(&theta)
+    };
+
+    let minimizer = NelderMeadBuilder::default()
+        .xtol(1e-6f64)
+        .ftol(1e-6f64)
+        .maxiter(50000)
+        .build()
+        .unwrap();
+
+    // Set the starting guess
+    let args = Array::from_vec(vec![3.0, -8.3]);
+    let xs = minimizer.minimize(&func_to_minimize, args.view());
+
+    xs
+}

src/nonparametric/elliptical_process/bayesian_optimization/upper_confidence_bound.rs

@@ -0,0 +1,18 @@
+use num_integer::sqrt;
+
+use super::AcquisitionFunctions;
+
+pub struct UpperConfidenceBound {
+    trial: f64,
+}
+
+impl AcquisitionFunctions for UpperConfidenceBound {
+    fn value(&self, theta: &crate::NormalParams) -> f64 {
+        let mu = theta.mu();
+        let sigma = theta.sigma();
+        let n = self.trial;
+        let k = sqrt(n.ln() / n);
+
+        mu + k * sigma
+    }
+}

src/nonparametric/elliptical_process/mod.rs

@@ -1,10 +1,12 @@
+pub mod bayesian_optimization;
 pub mod exact;
 pub mod ey;
 pub mod kernel_matrix;
 pub mod kiss_love;
 pub mod regressor;
 pub mod sparse;
 
+pub use bayesian_optimization::*;
 pub use exact::*;
 pub use ey::*;
 pub use kernel_matrix::*;