Projection on the epigraph of the squared Euclidean norm

CHANGELOG.md

@@ -14,6 +14,7 @@ Note: This is the main Changelog file for the Rust solver. The Changelog file fo
 
 ### Added
 
+- Rust implementation of epigraph of squared Euclidean norm (constraint) 
 - Implementation of `AffineSpace`
 
 ### Fixed

Cargo.toml

@@ -99,10 +99,16 @@ rpmalloc = { version = "0.2.0", features = [
 [target.'cfg(not(target_env = "msvc"))'.dependencies]
 jemallocator = { version = "0.5.0", optional = true }
 
+
+# computation of roots of cubic equation needed for the projection on the 
+# epigraph of the squared Eucliean norm
+roots = "0.0.8"
+
 # Least squares solver
 ndarray = { version = "0.15", features = ["approx"] }
 modcholesky = "0.1.3"
 
+
 # --------------------------------------------------------------------------
 # 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

src/constraints/epigraph_squared_norm.rs

@@ -0,0 +1,97 @@
+use crate::matrix_operations;
+
+use super::Constraint;
+
+#[derive(Copy, Clone)]
+/// 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;
+                }
+                return;
+            }
+        });
+
+        // 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 = zsol - p_z / dp_z;
+            refinement_error = p_z.abs();
+            iter += 1;
+        }
+        right_root = zsol;
+
+        // Projection
+        for i in 0..nx {
+            x[i] /= 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,7 @@
+use crate::matrix_operations;
+
 use super::*;
 use modcholesky::ModCholeskySE99;
 use rand;
 
 #[test]
@@ -875,7 +877,52 @@
 }
 
 #[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",
+    );
+}
 fn t_affine_space() {
     let a = vec![
         0.5, 0.1, 0.2, -0.3, -0.6, 0.3, 0., 0.5, 1.0, 0.1, -1.0, -0.4,
     ];