From 72d6c2b66d565e51a3c34fa24707942eb22fd0fb Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Wed, 11 Mar 2026 15:45:35 +0000 Subject: [PATCH 1/9] First implementation of ball_p projection About: - New file ballp.rs: impl of BallP projection - Very basic unit testing - Documentation Issues: - Addresses #333 --- src/constraints/ballp.rs | 301 +++++++++++++++++++++++++++++++++++++++ src/constraints/mod.rs | 2 + src/constraints/tests.rs | 14 ++ 3 files changed, 317 insertions(+) create mode 100644 src/constraints/ballp.rs diff --git a/src/constraints/ballp.rs b/src/constraints/ballp.rs new file mode 100644 index 00000000..c3bbaf95 --- /dev/null +++ b/src/constraints/ballp.rs @@ -0,0 +1,301 @@ +use super::Constraint; + +#[derive(Copy, Clone)] +/// An $\\ell_p$ ball, that is, +/// $$B_p^r = \\{ x \\in R^n : \Vert x \Vert_p \\leq r \\}$$ +/// 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 +/// +/// $$\Pi_{B_p^r}(x) = \mathrm{argmin}_{z \in B_p^r} \frac{1}{2}\Vert z - x \Vert_2^2.$$ +/// +/// If $x$ already belongs to the ball, the projection is $x$ itself. +/// Otherwise, the projection lies on the boundary of the ball. +/// +/// # Numerical method +/// +/// For general $1 < p < \infty$, the projection does not admit a simple +/// closed-form expression. This implementation computes it numerically using +/// the optimality conditions of the constrained problem. +/// +/// For a ball centered at the origin, the projection is characterized by +/// +/// $$z_i = \operatorname{sign}(x_i)\,u_i(\lambda),$$ +/// +/// where each scalar $u_i(\lambda) \geq 0$ solves +/// +/// $$u_i + \lambda p\,u_i^{p-1} = |x_i|,$$ +/// +/// and the Lagrange multiplier $\lambda \geq 0$ is chosen so that the projected +/// vector satisfies the feasibility condition $\|z\|_p = r.$ +/// +/// The implementation uses: +/// +/// - an *outer bisection loop* to determine the multiplier $\lambda$, +/// - an *inner safeguarded Newton method* to solve the scalar nonlinear +/// equation for each coordinate. +/// +/// The Newton iteration is combined with bracketing and a bisection fallback, +/// which improves robustness. +/// +/// # Translated balls +/// +/// If the ball has a center $x_c$, projection is performed by translating the +/// input vector to the origin, projecting onto the origin-centered ball, and +/// translating the result back: +/// +/// $$\Pi_{B_p^{x_c, r}}(x) = x_c + \Pi_{B_p^r}(x - x_c).$$ +/// +/// # Convexity +/// +/// Since this module assumes `p > 1.0`, every [`BallP`] is convex, and therefore +/// [`Constraint::project`] computes a unique Euclidean projection. +/// +/// # Examples +/// +/// Project onto the unit \(\ell_p\)-ball centered at the origin: +/// +/// ``` +/// use optimization_engine::constraints::{BallP, Constraint}; +/// +/// let ball = BallP::new(None, 1.0, 1.5, 1e-10, 100); +/// let mut x = vec![3.0, -1.0, 2.0]; +/// ball.project(&mut x); +/// ``` +/// +/// Project onto a translated \(\ell_p\)-ball: +/// +/// ``` +/// use optimization_engine::constraints::{BallP, Constraint}; +/// +/// let center = vec![1.0, 1.0, 1.0]; +/// let ball = BallP::new(Some(¢er), 2.0, 3.0, 1e-10, 100); +/// let mut x = vec![4.0, -1.0, 2.0]; +/// ball.project(&mut x); +/// ``` +/// +/// # Notes +/// +/// - The projection is with respect to the *Euclidean norm* +/// - 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 +/// in [`Ball2`](crate::constraints::Ball2) is more efficient +/// - The quality and speed of the computation depend on the chosen numerical +/// tolerance and iteration limit. +pub struct BallP<'a> { + /// Optional center of the ball. + /// + /// If `None`, the ball is centered at the origin. + /// If `Some(center)`, the ball is centered at `center`. + center: Option<&'a [f64]>, + + /// Radius of the ball. + /// + /// Must be strictly positive. + radius: f64, + + /// Exponent of the norm. + /// + /// Must satisfy `p > 1.0` and be finite. + p: f64, + + /// Numerical tolerance used by the outer bisection on the Lagrange + /// multiplier and by the inner Newton solver. + tolerance: f64, + + /// Maximum number of iterations used by the outer bisection and + /// the inner Newton solver. + max_iter: usize, +} + +impl<'a> BallP<'a> { + /// Construct a new l_p ball with given center, radius, and exponent. + /// + /// - `center`: if `None`, the ball is centered at the origin + /// - `radius`: radius of the ball + /// - `p`: norm exponent, must satisfy `p > 1.0` and be finite + /// - `tolerance`: tolerance for the numerical solvers + /// - `max_iter`: maximum number of iterations for the numerical solvers + pub fn new( + center: Option<&'a [f64]>, + radius: f64, + p: f64, + tolerance: f64, + max_iter: usize, + ) -> Self { + assert!(radius > 0.0); + assert!(p > 1.0 && p.is_finite()); + assert!(tolerance > 0.0); + assert!(max_iter > 0); + + BallP { + center, + radius, + p, + tolerance, + max_iter, + } + } + + #[inline] + fn lp_norm(x: &[f64], p: f64) -> f64 { + x.iter() + .map(|xi| xi.abs().powf(p)) + .sum::() + .powf(1.0 / p) + } + + #[inline] + fn project_origin(&self, x: &mut [f64]) { + self.project_lp_ball_general(x); + } + + fn project_lp_ball_general(&self, x: &mut [f64]) { + let p = self.p; + let r = self.radius; + let tol = self.tolerance; + let max_iter = self.max_iter; + + let current_norm = Self::lp_norm(x, p); + if current_norm <= r { + return; + } + + let abs_x: Vec = x.iter().map(|xi| xi.abs()).collect(); + let target = r.powf(p); + + let radius_error = |lambda: f64| -> f64 { + abs_x.iter() + .map(|&a| { + let u = Self::solve_coordinate_newton(a, lambda, p, tol, max_iter); + u.powf(p) + }) + .sum::() + - target + }; + + let mut lambda_lo = 0.0_f64; + let mut lambda_hi = 1.0_f64; + + while radius_error(lambda_hi) > 0.0 { + lambda_hi *= 2.0; + if lambda_hi > 1e20 { + panic!("Failed to bracket the Lagrange multiplier"); + } + } + + for _ in 0..max_iter { + let lambda_mid = 0.5 * (lambda_lo + lambda_hi); + let err = radius_error(lambda_mid); + + if err.abs() <= tol { + lambda_lo = lambda_mid; + lambda_hi = lambda_mid; + break; + } + + if err > 0.0 { + lambda_lo = lambda_mid; + } else { + lambda_hi = lambda_mid; + } + } + + let lambda_star = 0.5 * (lambda_lo + lambda_hi); + + x.iter_mut().zip(abs_x.iter()).for_each(|(xi, &a)| { + let u = Self::solve_coordinate_newton(a, lambda_star, p, tol, max_iter); + *xi = xi.signum() * u; + }); + } + + /// Solve for u >= 0 the equation u + lambda * p * u^(p-1) = a + /// using safeguarded Newton iterations. + /// + /// 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 { + if a == 0.0 { + return 0.0; + } + + if lambda == 0.0 { + return a; + } + + let mut lo = 0.0_f64; + let mut hi = a; + + // Heuristic initial guess: + // exact when p = 2, and usually in the right scale for general p. + let mut u = (a / (1.0 + lambda * p)).clamp(lo, hi); + + for _ in 0..max_iter { + let upm1 = u.powf(p - 1.0); + let f = u + lambda * p * upm1 - a; + + if f.abs() <= tol { + return u; + } + + if f > 0.0 { + hi = u; + } else { + lo = u; + } + + let df = 1.0 + lambda * p * (p - 1.0) * u.powf(p - 2.0); + let mut candidate = u - f / df; + + if !candidate.is_finite() || candidate <= lo || candidate >= hi { + candidate = 0.5 * (lo + hi); + } + + if (candidate - u).abs() <= tol * (1.0 + u.abs()) { + return candidate; + } + + u = candidate; + } + + 0.5 * (lo + hi) + } +} + +impl<'a> Constraint for BallP<'a> { + fn project(&self, x: &mut [f64]) { + if let Some(center) = &self.center { + assert_eq!(x.len(), center.len()); + + let mut shifted = vec![0.0; x.len()]; + shifted + .iter_mut() + .zip(x.iter().zip(center.iter())) + .for_each(|(s, (xi, ci))| *s = *xi - *ci); + + self.project_origin(&mut shifted); + + x.iter_mut() + .zip(shifted.iter().zip(center.iter())) + .for_each(|(xi, (si, ci))| *xi = *ci + *si); + } else { + self.project_origin(x); + } + } + + fn is_convex(&self) -> bool { + true + } +} \ No newline at end of file diff --git a/src/constraints/mod.rs b/src/constraints/mod.rs index dd0b7038..a0296aca 100644 --- a/src/constraints/mod.rs +++ b/src/constraints/mod.rs @@ -11,6 +11,7 @@ mod affine_space; mod ball1; mod ball2; +mod ballp; mod ballinf; mod cartesian_product; mod epigraph_squared_norm; @@ -27,6 +28,7 @@ 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 cartesian_product::CartesianProduct; pub use epigraph_squared_norm::EpigraphSquaredNorm; diff --git a/src/constraints/tests.rs b/src/constraints/tests.rs index 26343688..52064656 100644 --- a/src/constraints/tests.rs +++ b/src/constraints/tests.rs @@ -989,3 +989,17 @@ fn t_affine_space_wrong_dimensions() { let b = vec![1., 2., -0.5]; let _ = AffineSpace::new(a, b); } + + +#[test] +fn t_ballp_at_origin() { + let radius = 1.0; + let mut x = [1.0, 1.0]; + let p = 3.; + let tol = 1e-10; + let max_iters: usize = 200; + let ball = BallP::new(None, radius, p, tol, max_iters); + ball.project(&mut x); + println!("x = {:?}", x); + +} \ No newline at end of file From a11702bced4eeabd4d40eab6d9b014f08edc55fe Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Wed, 11 Mar 2026 17:13:42 +0000 Subject: [PATCH 2/9] 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 --- Cargo.toml | 1 + src/constraints/ballp.rs | 33 ++++----- src/constraints/mod.rs | 4 +- src/constraints/tests.rs | 150 +++++++++++++++++++++++++++++++++++++-- 4 files changed, 164 insertions(+), 24 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 70da22a4..1deeadb7 100644 --- a/Cargo.toml +++ b/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" # -------------------------------------------------------------------------- diff --git a/src/constraints/ballp.rs b/src/constraints/ballp.rs index c3bbaf95..89909007 100644 --- a/src/constraints/ballp.rs +++ b/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::() - .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 } -} \ No newline at end of file +} diff --git a/src/constraints/mod.rs b/src/constraints/mod.rs index a0296aca..d4ce2462 100644 --- a/src/constraints/mod.rs +++ b/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; diff --git a/src/constraints/tests.rs b/src/constraints/tests.rs index 52064656..5bf98629 100644 --- a/src/constraints/tests.rs +++ b/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` of length `dim`. +pub fn sample_lp_sphere(n_points: usize, dim: usize, p: f64, radius: f64) -> Vec> { + 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::() + .powf(1.0 / p); + + let point: Vec = 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::() + .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 = 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 = 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)); +} -} \ No newline at end of file +#[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", + ); +} From bc8316b18e36da52da3c95ce5346aa355008badc Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Wed, 11 Mar 2026 17:28:22 +0000 Subject: [PATCH 3/9] Better unit tests (lp-ball) About: - testing that p-norm of projection is <= radius - made lp_norm non-public function in BallP --- src/constraints/ballp.rs | 3 +-- src/constraints/tests.rs | 25 +++++++++++++++---------- 2 files changed, 16 insertions(+), 12 deletions(-) diff --git a/src/constraints/ballp.rs b/src/constraints/ballp.rs index 89909007..6a332ece 100644 --- a/src/constraints/ballp.rs +++ b/src/constraints/ballp.rs @@ -144,13 +144,12 @@ impl<'a> BallP<'a> { } #[inline] - #[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 { + fn lp_norm(&self, x: &[f64]) -> f64 { x.iter() .map(|xi| xi.abs().powf(self.p)) .sum::() diff --git a/src/constraints/tests.rs b/src/constraints/tests.rs index 5bf98629..12ba01e2 100644 --- a/src/constraints/tests.rs +++ b/src/constraints/tests.rs @@ -992,9 +992,7 @@ 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 } +/// Sample `n_points` random points on the `l_p` sphere in `R^dim` /// /// Assumes: /// - `p > 0.0` @@ -1082,6 +1080,18 @@ fn is_norm_p_projection( ) -> bool { let n = x.len(); assert_eq!(n, x_candidate_proj.len()); + + // Make sure ||x||_p ≤ radius + ε, where ε is a small tolerance + let feasibility_tol = 1e-10; + let norm_proj = x_candidate_proj + .iter() + .map(|xi| xi.abs().powf(p)) + .sum::() + .powf(1.0 / p); + if norm_proj > radius + feasibility_tol { + return false; + } + // e = x - x_candidate_proj let e: Vec = x .iter() @@ -1113,12 +1123,7 @@ fn t_ballp_at_origin_projection() { let tol = 1e-16; let max_iters: usize = 200; let ball = BallP::new(None, radius, p, tol, max_iters); - ball.project(&mut 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"); + ball.project(&mut x); assert!(is_norm_p_projection(&x0, &x, p, radius, 10_000)); } @@ -1132,7 +1137,6 @@ fn t_ballp_at_origin_x_already_inside() { 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, @@ -1141,3 +1145,4 @@ fn t_ballp_at_origin_x_already_inside() { "wrong projection on lp-ball", ); } + From ebc479300d073d8d0b72722822cbbfd6a7bc2ed9 Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Wed, 11 Mar 2026 17:28:42 +0000 Subject: [PATCH 4/9] Run cargo-fmt --- src/constraints/tests.rs | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/constraints/tests.rs b/src/constraints/tests.rs index 12ba01e2..e81d7ef2 100644 --- a/src/constraints/tests.rs +++ b/src/constraints/tests.rs @@ -1123,7 +1123,7 @@ fn t_ballp_at_origin_projection() { let tol = 1e-16; let max_iters: usize = 200; let ball = BallP::new(None, radius, p, tol, max_iters); - ball.project(&mut x); + ball.project(&mut x); assert!(is_norm_p_projection(&x0, &x, p, radius, 10_000)); } @@ -1145,4 +1145,3 @@ fn t_ballp_at_origin_x_already_inside() { "wrong projection on lp-ball", ); } - From 35ae945febd0e24fcc6dec016ebe8de3f9dfe641 Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Wed, 11 Mar 2026 17:35:07 +0000 Subject: [PATCH 5/9] Numerical stability of unit tests About: - Assert that inner product is < tolerance rather than <0 --- src/constraints/tests.rs | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/constraints/tests.rs b/src/constraints/tests.rs index e81d7ef2..7b926d82 100644 --- a/src/constraints/tests.rs +++ b/src/constraints/tests.rs @@ -1039,7 +1039,7 @@ pub fn sample_lp_sphere(n_points: usize, dim: usize, p: f64, radius: f64) -> Vec #[test] fn t_sample_lp_sphere_points_have_correct_norm() { - let n_points = 100_000; + let n_points = 10_000; let dim = 5; let p = 3.2; let radius = 0.9; @@ -1083,6 +1083,7 @@ fn is_norm_p_projection( // Make sure ||x||_p ≤ radius + ε, where ε is a small tolerance let feasibility_tol = 1e-10; + let inner_prod_tol = 1e-10; let norm_proj = x_candidate_proj .iter() .map(|xi| xi.abs().powf(p)) @@ -1107,7 +1108,7 @@ fn is_norm_p_projection( .map(|(xproj_i, xi_i)| xproj_i - xi_i) .collect(); let inner = matrix_operations::inner_product(&w, &e); - if inner < 0. { + if inner < inner_prod_tol { return false; } } From 39886d4648faaede1f15c1caf6126b6933028d9b Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Wed, 11 Mar 2026 23:03:29 +0000 Subject: [PATCH 6/9] New test t_ballp_at_xc_projection About: - Unit test with lp-ball centered not at origin - Clippy issues fixed --- src/constraints/ballp.rs | 18 +++++++----------- src/constraints/tests.rs | 39 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 46 insertions(+), 11 deletions(-) diff --git a/src/constraints/ballp.rs b/src/constraints/ballp.rs index 6a332ece..7510ba8a 100644 --- a/src/constraints/ballp.rs +++ b/src/constraints/ballp.rs @@ -82,10 +82,11 @@ use super::Constraint; /// /// - The projection is with respect to the *Euclidean norm* /// - 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). +/// 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 -/// in [`Ball2`](crate::constraints::Ball2) is more efficient +/// in [`Ball2`](crate::constraints::Ball2) is more efficient /// - The quality and speed of the computation depend on the chosen numerical /// tolerance and iteration limit. pub struct BallP<'a> { @@ -156,12 +157,7 @@ impl<'a> BallP<'a> { .powf(1.0 / self.p) } - #[inline] - fn project_origin(&self, x: &mut [f64]) { - self.project_lp_ball_general(x); - } - - fn project_lp_ball_general(&self, x: &mut [f64]) { + fn project_lp_ball(&self, x: &mut [f64]) { let p = self.p; let r = self.radius; let tol = self.tolerance; @@ -285,13 +281,13 @@ impl<'a> Constraint for BallP<'a> { .zip(x.iter().zip(center.iter())) .for_each(|(s, (xi, ci))| *s = *xi - *ci); - self.project_origin(&mut shifted); + self.project_lp_ball(&mut shifted); x.iter_mut() .zip(shifted.iter().zip(center.iter())) .for_each(|(xi, (si, ci))| *xi = *ci + *si); } else { - self.project_origin(x); + self.project_lp_ball(x); } } diff --git a/src/constraints/tests.rs b/src/constraints/tests.rs index 7b926d82..7e224a12 100644 --- a/src/constraints/tests.rs +++ b/src/constraints/tests.rs @@ -1071,6 +1071,9 @@ fn t_sample_lp_sphere_points_have_correct_norm() { /// 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. +/// +/// Note that this test is stochastic, so it only guarantees the correctness +/// of the projection in probability. fn is_norm_p_projection( x: &[f64], x_candidate_proj: &[f64], @@ -1146,3 +1149,39 @@ fn t_ballp_at_origin_x_already_inside() { "wrong projection on lp-ball", ); } + +#[test] +fn t_ballp_at_xc_projection() { + let radius = 0.8; + let mut x = [0.0, 0.1]; + let x_center = [1.0, 3.0]; + let p = 4.; + let tol = 1e-16; + let max_iters: usize = 200; + let ball = BallP::new( + Some(&x_center), + radius, + p, + tol, + max_iters); + ball.project(&mut x); + + let nrm = (x.iter() + .zip(x_center.iter()) + .fold(0.0, |s, (x, y)| { + (*x - *y).abs().powf(p) + s + })).powf(1./p); + unit_test_utils::assert_nearly_equal( + radius, nrm, + 1e-10, 1e-12, "wrong distance to lp-ball"); + + let proj_expected = [0.5178727276722618, 2.2277981662325224]; + unit_test_utils::assert_nearly_equal_array( + &proj_expected, + &x, + 1e-12, + 1e-12, + "wrong projection on lp-ball centered at xc != 0", + ); + +} From d8839a386ff60170fb6efbfe9c224cd020a58ba1 Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Wed, 11 Mar 2026 23:13:42 +0000 Subject: [PATCH 7/9] Cargo fmt (last time) Cargo fmt set up in VS Code --- src/constraints/ballp.rs | 4 ++-- src/constraints/tests.rs | 24 ++++++++---------------- 2 files changed, 10 insertions(+), 18 deletions(-) diff --git a/src/constraints/ballp.rs b/src/constraints/ballp.rs index 7510ba8a..46d52d80 100644 --- a/src/constraints/ballp.rs +++ b/src/constraints/ballp.rs @@ -82,8 +82,8 @@ use super::Constraint; /// /// - The projection is with respect to the *Euclidean norm* /// - 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) +/// 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 /// in [`Ball2`](crate::constraints::Ball2) is more efficient diff --git a/src/constraints/tests.rs b/src/constraints/tests.rs index 7e224a12..0bb32c43 100644 --- a/src/constraints/tests.rs +++ b/src/constraints/tests.rs @@ -1071,7 +1071,7 @@ fn t_sample_lp_sphere_points_have_correct_norm() { /// 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. -/// +/// /// Note that this test is stochastic, so it only guarantees the correctness /// of the projection in probability. fn is_norm_p_projection( @@ -1158,22 +1158,15 @@ fn t_ballp_at_xc_projection() { let p = 4.; let tol = 1e-16; let max_iters: usize = 200; - let ball = BallP::new( - Some(&x_center), - radius, - p, - tol, - max_iters); + let ball = BallP::new(Some(&x_center), radius, p, tol, max_iters); ball.project(&mut x); - - let nrm = (x.iter() + + let nrm = (x + .iter() .zip(x_center.iter()) - .fold(0.0, |s, (x, y)| { - (*x - *y).abs().powf(p) + s - })).powf(1./p); - unit_test_utils::assert_nearly_equal( - radius, nrm, - 1e-10, 1e-12, "wrong distance to lp-ball"); + .fold(0.0, |s, (x, y)| (*x - *y).abs().powf(p) + s)) + .powf(1. / p); + unit_test_utils::assert_nearly_equal(radius, nrm, 1e-10, 1e-12, "wrong distance to lp-ball"); let proj_expected = [0.5178727276722618, 2.2277981662325224]; unit_test_utils::assert_nearly_equal_array( @@ -1183,5 +1176,4 @@ fn t_ballp_at_xc_projection() { 1e-12, "wrong projection on lp-ball centered at xc != 0", ); - } From 6f6a547506d7490a3bfc5d1eded0deb8a0bf4aba Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Fri, 13 Mar 2026 17:19:24 +0000 Subject: [PATCH 8/9] [ci skip] documentation on website for BallP (rust) --- docs/openrust-basic.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/openrust-basic.md b/docs/openrust-basic.md index 412f7bd0..d737c80a 100644 --- a/docs/openrust-basic.md +++ b/docs/openrust-basic.md @@ -65,6 +65,7 @@ Constraints implement the namesake trait, [`Constraint`]. Implementations of [`C | [`AffineSpace`] | $U {}={} \\{u\in\mathbb{R}^n : Au = b\\}$ | | [`Ball1`] | $U {}={} \\{u\in\mathbb{R}^n : \Vert u-u^0\Vert_1 \leq r\\}$ | | [`Ball2`] | $U {}={} \\{u\in\mathbb{R}^n : \Vert u-u^0\Vert_2 \leq r\\}$ | +| [`BallP`] | $U {}={} \\{u\in\mathbb{R}^n : \Vert u-u^0\Vert_p \leq r\\}$ | | [`Sphere2`] | $U {}={} \\{u\in\mathbb{R}^n : \Vert u-u^0\Vert_2 = r\\}$ | | [`BallInf`] | $U {}={} \\{u\in\mathbb{R}^n : \Vert u-u^0\Vert_\infty \leq r\\}$ | | [`Halfspace`] | $U {}={} \\{u\in\mathbb{R}^n : \langle c, u\rangle \leq b\\}$ | @@ -362,6 +363,7 @@ the imposition of a maximum allowed duration, the exit status will be [`Sphere2`]: https://docs.rs/optimization_engine/*/optimization_engine/constraints/struct.Sphere2.html [`Ball1`]: https://docs.rs/optimization_engine/*/optimization_engine/constraints/struct.Ball1.html [`Ball2`]: https://docs.rs/optimization_engine/*/optimization_engine/constraints/struct.Ball2.html +[`BallP`]: https://docs.rs/optimization_engine/*/optimization_engine/constraints/struct.BallP.html [`BallInf`]: https://docs.rs/optimization_engine/*/optimization_engine/constraints/struct.BallInf.html [`Halfspace`]: https://docs.rs/optimization_engine/*/optimization_engine/constraints/struct.Halfspace.html [`Hyperplane`]: https://docs.rs/optimization_engine/*/optimization_engine/constraints/struct.Hyperplane.html From 351748d21fd97ed6ad239acba90f52c84204381e Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Fri, 13 Mar 2026 17:23:31 +0000 Subject: [PATCH 9/9] update Cargo.toml and CHANGELOG.md (v0.11.0) --- CHANGELOG.md | 10 ++++++++++ Cargo.toml | 2 +- 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 83c1cb8f..385e576e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,15 @@ 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.11.0] - Unreleased + +### Added + +- Implementation of `BallP` in Rust: projection on lp-ball + +[v0.11.0]: https://github.com/alphaville/optimization-engine/compare/v0.10.0...v0.11.0 [v0.10.0]: https://github.com/alphaville/optimization-engine/compare/v0.9.1...v0.10.0 [v0.9.1]: https://github.com/alphaville/optimization-engine/compare/v0.9.0...v0.9.1 [v0.9.0]: https://github.com/alphaville/optimization-engine/compare/v0.8.1...v0.9.0 diff --git a/Cargo.toml b/Cargo.toml index 1deeadb7..450ec90f 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -42,7 +42,7 @@ homepage = "https://alphaville.github.io/optimization-engine/" repository = "https://github.com/alphaville/optimization-engine" # Version of this crate (SemVer) -version = "0.10.0" +version = "0.11.0" edition = "2018"