diff --git a/CHANGELOG.md b/CHANGELOG.md index 08f21fb0..076bda35 100644 --- a/CHANGELOG.md +++ b/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) + + + + -## [v0.9.0] - Unreleased +## [v0.9.0] - 2024-03-20 ### Added +- Rust implementation of epigraph of squared Euclidean norm (constraint) - Implementation of `AffineSpace` ### Fixed diff --git a/Cargo.toml b/Cargo.toml index 4e4ea073..ae50f398 100644 --- a/Cargo.toml +++ b/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,10 +86,10 @@ 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 } @@ -97,11 +97,17 @@ rpmalloc = { version = "0.2.0", features = [ # 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. diff --git a/appveyor.yml b/appveyor.yml index 182dfc5f..e216a306 100644 --- a/appveyor.yml +++ b/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 diff --git a/docs/python-interface.md b/docs/python-interface.md index a332ea29..77983d8a 100644 --- a/docs/python-interface.md +++ b/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. | diff --git a/src/constraints/epigraph_squared_norm.rs b/src/constraints/epigraph_squared_norm.rs new file mode 100644 index 00000000..560b34bd --- /dev/null +++ b/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 + } +} diff --git a/src/constraints/mod.rs b/src/constraints/mod.rs index 2a02046d..dd0b7038 100644 --- a/src/constraints/mod.rs +++ b/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; diff --git a/src/constraints/sphere2.rs b/src/constraints/sphere2.rs index 2befc82b..8ea4d344 100644 --- a/src/constraints/sphere2.rs +++ b/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(¢er); + x.copy_from_slice(center); x[0] += self.radius; return; } diff --git a/src/constraints/tests.rs b/src/constraints/tests.rs index eaab661c..26343688 100644 --- a/src/constraints/tests.rs +++ b/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![