More thorough testing for lp-ball projection

About:
- test projection onto lp-ball using random points on the sphere
- introduce rand_distr as a test dependency
- unit test of test function that generated random points
- test projection when x is inside the ball

Author: alphaville

Cargo.toml

@@ -133,6 +133,7 @@ unit_test_utils = "0.1.4"
 icasadi_test = "0.0.3"
 # Random number generators for unit tests:
 rand = "0.10"
+rand_distr = "0.6"
 
 
 # --------------------------------------------------------------------------

src/constraints/ballp.rs

@@ -6,7 +6,7 @@ use super::Constraint;
 /// or a translated ball
 /// $$B_p^{x_c, r} = \\{ x \\in \mathbb{R}^n : \Vert x - x_c \Vert_p  \\leq r \\},$$
 /// with $1 < p < \infty$.
-/// 
+///
 /// # Projection problem
 ///
 /// Given a vector $x$, projection onto the ball means solving
@@ -81,10 +81,10 @@ use super::Constraint;
 /// # Notes
 ///
 /// - The projection is with respect to the *Euclidean norm*
-/// - The implementation is intended for general finite $p > 1.0$. If you need 
+/// - The implementation is intended for general finite $p > 1.0$. If you need
 /// to project on a $\Vert{}\cdot{}\Vert_1$-ball or an $\Vert{}\cdot{}
 /// \Vert_\infty$-ball, use the implementations in [`Ball1`](crate::constraints::Ball1) and [`BallInf`](crate::constraints::BallInf).
-/// - Do not use this struct to project on a Euclidean ball; the implementation 
+/// - Do not use this struct to project on a Euclidean ball; the implementation
 /// in [`Ball2`](crate::constraints::Ball2) is more efficient
 /// - The quality and speed of the computation depend on the chosen numerical
 ///   tolerance and iteration limit.
@@ -144,11 +144,17 @@ impl<'a> BallP<'a> {
     }
 
     #[inline]
-    fn lp_norm(x: &[f64], p: f64) -> f64 {
+    #[must_use]
+    /// Computes the $p$-norm of a given vector
+    ///
+    /// The $p$-norm of a vector $x\in \mathbb{R}^n$ is given by
+    /// $$\Vert x \Vert_p = \left(\sum_{i=1}^{n} |x_i|^p\right)^{1/p},$$
+    /// for $p > 1$.
+    pub fn lp_norm(&self, x: &[f64]) -> f64 {
         x.iter()
-            .map(|xi| xi.abs().powf(p))
+            .map(|xi| xi.abs().powf(self.p))
             .sum::<f64>()
-            .powf(1.0 / p)
+            .powf(1.0 / self.p)
     }
 
     #[inline]
@@ -162,7 +168,7 @@ impl<'a> BallP<'a> {
         let tol = self.tolerance;
         let max_iter = self.max_iter;
 
-        let current_norm = Self::lp_norm(x, p);
+        let current_norm = self.lp_norm(x);
         if current_norm <= r {
             return;
         }
@@ -171,7 +177,8 @@ impl<'a> BallP<'a> {
         let target = r.powf(p);
 
         let radius_error = |lambda: f64| -> f64 {
-            abs_x.iter()
+            abs_x
+                .iter()
                 .map(|&a| {
                     let u = Self::solve_coordinate_newton(a, lambda, p, tol, max_iter);
                     u.powf(p)
@@ -220,13 +227,7 @@ impl<'a> BallP<'a> {
     ///
     /// The solution always belongs to [0, a], so Newton is combined with
     /// bracketing and a bisection fallback.
-    fn solve_coordinate_newton(
-        a: f64,
-        lambda: f64,
-        p: f64,
-        tol: f64,
-        max_iter: usize,
-    ) -> f64 {
+    fn solve_coordinate_newton(a: f64, lambda: f64, p: f64, tol: f64, max_iter: usize) -> f64 {
         if a == 0.0 {
             return 0.0;
         }
@@ -298,4 +299,4 @@ impl<'a> Constraint for BallP<'a> {
     fn is_convex(&self) -> bool {
         true
     }
-}
+}

src/constraints/mod.rs

@@ -11,8 +11,8 @@
 mod affine_space;
 mod ball1;
 mod ball2;
-mod ballp;
 mod ballinf;
+mod ballp;
 mod cartesian_product;
 mod epigraph_squared_norm;
 mod finite;
@@ -28,8 +28,8 @@ mod zero;
 pub use affine_space::AffineSpace;
 pub use ball1::Ball1;
 pub use ball2::Ball2;
-pub use ballp::BallP;
 pub use ballinf::BallInf;
+pub use ballp::BallP;
 pub use cartesian_product::CartesianProduct;
 pub use epigraph_squared_norm::EpigraphSquaredNorm;
 pub use finite::FiniteSet;

src/constraints/tests.rs

@@ -2,6 +2,8 @@ use crate::matrix_operations;
 
 use super::*;
 use rand;
+use rand::RngExt;
+use rand_distr::{Distribution, Gamma};
 
 #[test]
 fn t_zero_set() {
@@ -990,16 +992,152 @@ fn t_affine_space_wrong_dimensions() {
     let _ = AffineSpace::new(a, b);
 }
 
+/// Sample `n_points` random points on the l_p sphere in R^dim:
+///
+///     { x : ||x||_p = radius }
+///
+/// Assumes:
+/// - `p > 0.0`
+/// - `radius > 0.0`
+/// - `dim > 0`
+///
+/// Returns a vector of points, each point being a `Vec<f64>` of length `dim`.
+pub fn sample_lp_sphere(n_points: usize, dim: usize, p: f64, radius: f64) -> Vec<Vec<f64>> {
+    assert!(n_points > 0);
+    assert!(dim > 0);
+    assert!(p > 0.0);
+    assert!(radius > 0.0);
+
+    let mut rng = rand::rng();
+
+    // Y_i ~ Gamma(shape=1/p, scale=1)
+    let gamma = Gamma::new(1.0 / p, 1.0).unwrap();
+
+    let mut points = Vec::with_capacity(n_points);
+
+    for _ in 0..n_points {
+        let mut g = Vec::with_capacity(dim);
+
+        for _ in 0..dim {
+            let y: f64 = gamma.sample(&mut rng);
+            let sign = if rng.random_bool(0.5) { 1.0 } else { -1.0 };
+            let value = sign * y.powf(1.0 / p);
+            g.push(value);
+        }
+
+        let norm_p = g
+            .iter()
+            .map(|x: &f64| x.abs().powf(p))
+            .sum::<f64>()
+            .powf(1.0 / p);
+
+        let point: Vec<f64> = g.into_iter().map(|x| radius * x / norm_p).collect();
+
+        points.push(point);
+    }
+
+    points
+}
 
 #[test]
-fn t_ballp_at_origin() {
-    let radius = 1.0;
-    let mut x = [1.0, 1.0];
+fn t_sample_lp_sphere_points_have_correct_norm() {
+    let n_points = 100_000;
+    let dim = 5;
+    let p = 3.2;
+    let radius = 0.9;
+
+    let samples = sample_lp_sphere(n_points, dim, p, radius);
+
+    assert_eq!(samples.len(), n_points);
+
+    for point in samples.iter() {
+        assert_eq!(point.len(), dim);
+
+        let norm_p = point
+            .iter()
+            .map(|xi| xi.abs().powf(p))
+            .sum::<f64>()
+            .powf(1.0 / p);
+
+        unit_test_utils::assert_nearly_equal(
+            radius,
+            norm_p,
+            1e-10,
+            1e-12,
+            "sample_lp_sphere produced a point with incorrect p-norm",
+        );
+    }
+}
+
+/// Check if a given vector `x_candidate_proj` is actually the projection
+/// of `x` onto the p-norm ball centered at the origin with a given radius
+///
+/// This is based on taking `sample_points` on the sphere of the lp-ball.
+fn is_norm_p_projection(
+    x: &[f64],
+    x_candidate_proj: &[f64],
+    p: f64,
+    radius: f64,
+    sample_points: usize,
+) -> bool {
+    let n = x.len();
+    assert_eq!(n, x_candidate_proj.len());
+    // e = x - x_candidate_proj
+    let e: Vec<f64> = x
+        .iter()
+        .zip(x_candidate_proj.iter())
+        .map(|(xi, yi)| xi - yi)
+        .collect();
+    let samples = sample_lp_sphere(sample_points, n, p, radius);
+    for xi in samples.iter() {
+        // w = x_candidate_proj - xi
+        let w: Vec<f64> = x_candidate_proj
+            .iter()
+            .zip(xi.iter())
+            .map(|(xproj_i, xi_i)| xproj_i - xi_i)
+            .collect();
+        let inner = matrix_operations::inner_product(&w, &e);
+        if inner < 0. {
+            return false;
+        }
+    }
+    true
+}
+
+#[test]
+fn t_ballp_at_origin_projection() {
+    let radius = 0.8;
+    let mut x = [1.0, -1.0, 6.0];
+    let x0 = x.clone();
     let p = 3.;
-    let tol = 1e-10;
+    let tol = 1e-16;
     let max_iters: usize = 200;
     let ball = BallP::new(None, radius, p, tol, max_iters);
     ball.project(&mut x);
-    println!("x = {:?}", x);
+    let p_norm_projected = ball.lp_norm(&x);
+    // Test that the p-norm of x is equal to the radius
+    // (this is because x was not in the lp-ball, so its
+    //  projection is on the sphere)
+    unit_test_utils::assert_nearly_equal(radius, p_norm_projected, 1e-10, 1e-15, "lol");
+    assert!(is_norm_p_projection(&x0, &x, p, radius, 10_000));
+}
 
-}
+#[test]
+fn t_ballp_at_origin_x_already_inside() {
+    let radius = 1.5;
+    let mut x = [0.5, -0.2, 0.1];
+    let x0 = x.clone();
+    let p = 3.;
+    let tol = 1e-16;
+    let max_iters: usize = 1200;
+    let ball = BallP::new(None, radius, p, tol, max_iters);
+    ball.project(&mut x);
+    assert!(ball.lp_norm(&x) <= radius);
+    unit_test_utils::assert_nearly_equal_array(
+        &x0,
+        &x,
+        1e-12,
+        1e-12,
+        "wrong projection on lp-ball",
+    );
+}