Merge pull request #334 from alphaville/feature/projection-epigraph-norm

Projection on the epigraph of the squared Euclidean norm

CHANGELOG.md

@@ -7,13 +7,20 @@ and this project adheres to [Semantic Versioning](http://semver.org/).
 
 Note: This is the main Changelog file for the Rust solver. The Changelog file for the Python interface (`opengen`) can be found in [/open-codegen/CHANGELOG.md](open-codegen/CHANGELOG.md)
 
+
+<!-- ---------------------
+      Unreleased
+     --------------------- -->
+
+
 <!-- ---------------------
-      Not released
+      v0.9.0
      --------------------- -->
-## [v0.9.0] - Unreleased
+## [v0.9.0] - 2024-03-20
 
 ### Added
 
+- Rust implementation of epigraph of squared Euclidean norm (constraint) 
 - Implementation of `AffineSpace`
 
 ### Fixed

Cargo.toml

@@ -77,7 +77,7 @@ rustdoc-args = ["--html-in-header", "katex-header.html"]
 # D.E.P.E.N.D.E.N.C.I.E.S
 # --------------------------------------------------------------------------
 [dependencies]
-num = "0.4.0"
+num = "0.4"
 
 # Our own stuff - L-BFGS: limited-memory BFGS directions
 lbfgs = "0.2"
@@ -86,22 +86,28 @@ lbfgs = "0.2"
 instant = { version = "0.1" }
 
 # Wasm-bindgen is only activated if OpEn is compiled with `--features wasm`
-wasm-bindgen = { version = "0.2.74", optional = true }
+wasm-bindgen = { version = "0.2", optional = true }
 
 # sc-allocator provides an implementation of a bump allocator
-rpmalloc = { version = "0.2.0", features = [
+rpmalloc = { version = "0.2", features = [
     "guards",
     "statistics",
 ], optional = true }
 
 # jemallocator is an optional feature; it will only be loaded if the feature 
 # `jem` is used (i.e., if we compile with `cargo build --features jem`)
 [target.'cfg(not(target_env = "msvc"))'.dependencies]
-jemallocator = { version = "0.5.0", optional = true }
+jemallocator = { version = "0.5", optional = true }
+
+
+# computation of roots of cubic equation needed for the projection on the 
+# epigraph of the squared Euclidean norm
+roots = "0.0.8"
 
 # Least squares solver
 ndarray = { version = "0.15", features = ["approx"] }
-modcholesky = "0.1.3"
+modcholesky = "0.1"
+
 
 # --------------------------------------------------------------------------
 # F.E.A.T.U.R.E.S.

appveyor.yml

@@ -68,6 +68,7 @@ build: false
 #directly or perform other testing commands. Rust will automatically be placed in the PATH
 # environment variable.
 test_script:
+  - cargo add roots
   - cargo add ndarray --features approx
   - cargo add modcholesky
   - cargo build

docs/python-interface.md

@@ -87,6 +87,7 @@ following types of constraints:
 | `NoConstraints`    | No constraints - the whole $\mathbb{R}^{n}$|
 | `Rectangle`        | Rectangle, $$R = \\{u \in \mathbb{R}^{n_u} {}:{} f_{\min} \leq u \leq f_{\max}\\},$$ for example, `Rectangle(fmin, fmax)` |
 | `SecondOrderCone`  | Second-order aka "ice cream" aka "Lorenz" cone |
+| `EpigraphSquaredNorm`| The epigraph of the squared Eucliden norm is a set of the form $X = \\{(z, t) \in \mathbb{R}^{n+1}:  \Vert z \Vert  \leq t\\}$. | 
 | `CartesianProduct` | Cartesian product of any of the above. See more information below. |
 
 

src/constraints/epigraph_squared_norm.rs

@@ -0,0 +1,96 @@
+use crate::matrix_operations;
+
+use super::Constraint;
+
+#[derive(Copy, Clone, Default)]
+/// The epigraph of the squared Eucliden norm is a set of the form
+/// $X = \\{x = (z, t) \in \mathbb{R}^{n}\times \mathbb{R} {}:{} \\|z\\|^2 \leq t \\}.$
+pub struct EpigraphSquaredNorm {}
+
+impl EpigraphSquaredNorm {
+    /// Create a new instance of the epigraph of the squared norm.
+    ///
+    /// Note that you do not need to specify the dimension.
+    pub fn new() -> Self {
+        EpigraphSquaredNorm {}
+    }
+}
+
+impl Constraint for EpigraphSquaredNorm {
+    ///Project on the epigraph of the squared Euclidean norm.
+    ///
+    /// The projection is computed as detailed
+    /// [here](https://mathematix.wordpress.com/2017/05/02/projection-on-the-epigraph-of-the-squared-euclidean-norm/).
+    ///
+    /// ## Arguments
+    /// - `x`: The given vector $x$ is updated with the projection on the set
+    ///
+    /// ## Example
+    ///
+    /// ```rust
+    /// use optimization_engine::constraints::*;
+    ///
+    /// let epi = EpigraphSquaredNorm::new();
+    /// let mut x = [1., 2., 3., 4.];    
+    /// epi.project(&mut x);
+    /// ```
+    fn project(&self, x: &mut [f64]) {
+        let nx = x.len() - 1;
+        assert!(nx > 0, "x must have a length of at least 2");
+        let z: &[f64] = &x[..nx];
+        let t: f64 = x[nx];
+        let norm_z_sq = matrix_operations::norm2_squared(z);
+        if norm_z_sq <= t {
+            return;
+        }
+
+        let theta = 1. - 2. * t;
+        let a3 = 4.;
+        let a2 = 4. * theta;
+        let a1 = theta * theta;
+        let a0 = -norm_z_sq;
+
+        let cubic_poly_roots = roots::find_roots_cubic(a3, a2, a1, a0);
+        let mut right_root = f64::NAN;
+        let mut scaling = f64::NAN;
+
+        // Find right root
+        cubic_poly_roots.as_ref().iter().for_each(|ri| {
+            if *ri > 0. {
+                let denom = 1. + 2. * (*ri - t);
+                if ((norm_z_sq / (denom * denom)) - *ri).abs() < 1e-6 {
+                    right_root = *ri;
+                    scaling = denom;
+                }
+            }
+        });
+
+        // Refinement of root with Newton-Raphson
+        let mut refinement_error = 1.;
+        let newton_max_iters: usize = 5;
+        let newton_eps = 1e-14;
+        let mut zsol = right_root;
+        let mut iter = 0;
+        while refinement_error > newton_eps && iter < newton_max_iters {
+            let zsol_sq = zsol * zsol;
+            let zsol_cb = zsol_sq * zsol;
+            let p_z = a3 * zsol_cb + a2 * zsol_sq + a1 * zsol + a0;
+            let dp_z = 3. * a3 * zsol_sq + 2. * a2 * zsol + a1;
+            zsol -= p_z / dp_z;
+            refinement_error = p_z.abs();
+            iter += 1;
+        }
+        right_root = zsol;
+
+        // Projection
+        for xi in x.iter_mut().take(nx) {
+            *xi /= scaling;
+        }
+        x[nx] = right_root;
+    }
+
+    /// This is a convex set, so this function returns `True`
+    fn is_convex(&self) -> bool {
+        true
+    }
+}

src/constraints/mod.rs

@@ -13,6 +13,7 @@ mod ball1;
 mod ball2;
 mod ballinf;
 mod cartesian_product;
+mod epigraph_squared_norm;
 mod finite;
 mod halfspace;
 mod hyperplane;
@@ -28,6 +29,7 @@ pub use ball1::Ball1;
 pub use ball2::Ball2;
 pub use ballinf::BallInf;
 pub use cartesian_product::CartesianProduct;
+pub use epigraph_squared_norm::EpigraphSquaredNorm;
 pub use finite::FiniteSet;
 pub use halfspace::Halfspace;
 pub use hyperplane::Hyperplane;

src/constraints/sphere2.rs

@@ -38,7 +38,7 @@ impl<'a> Constraint for Sphere2<'a> {
         if let Some(center) = &self.center {
             let norm_difference = crate::matrix_operations::norm2_squared_diff(x, center).sqrt();
             if norm_difference <= epsilon {
-                x.copy_from_slice(&center);
+                x.copy_from_slice(center);
                 x[0] += self.radius;
                 return;
             }

src/constraints/tests.rs

@@ -1,5 +1,6 @@
+use crate::matrix_operations;
+
 use super::*;
-use modcholesky::ModCholeskySE99;
 use rand;
 
 #[test]
@@ -874,6 +875,53 @@ fn t_ball1_alpha_negative() {
     let _ = Ball1::new(None, -1.);
 }
 
+#[test]
+fn t_epigraph_squared_norm_inside() {
+    let epi = EpigraphSquaredNorm::new();
+    let mut x = [1., 2., 10.];
+    let x_correct = x.clone();
+    epi.project(&mut x);
+    unit_test_utils::assert_nearly_equal_array(
+        &x_correct,
+        &x,
+        1e-12,
+        1e-14,
+        "wrong projection on epigraph of squared norm",
+    );
+}
+
+#[test]
+fn t_epigraph_squared_norm() {
+    let epi = EpigraphSquaredNorm::new();
+    for i in 0..100 {
+        let t = 0.01 * i as f64;
+        let mut x = [1., 2., 3., t];
+        epi.project(&mut x);
+        let err = (matrix_operations::norm2_squared(&x[..3]) - x[3]).abs();
+        assert!(err < 1e-10, "wrong projection on epigraph of squared norm");
+    }
+}
+
+#[test]
+fn t_epigraph_squared_norm_correctness() {
+    let epi = EpigraphSquaredNorm::new();
+    let mut x = [1., 2., 3., 4.];
+    let x_correct = [
+        0.560142228903570,
+        1.120284457807140,
+        1.680426686710711,
+        4.392630432414829,
+    ];
+    epi.project(&mut x);
+    unit_test_utils::assert_nearly_equal_array(
+        &x_correct,
+        &x,
+        1e-12,
+        1e-14,
+        "wrong projection on epigraph of squared norm",
+    );
+}
+
 #[test]
 fn t_affine_space() {
     let a = vec![