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 70da22a4..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" @@ -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/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 diff --git a/src/constraints/ballp.rs b/src/constraints/ballp.rs new file mode 100644 index 00000000..46d52d80 --- /dev/null +++ b/src/constraints/ballp.rs @@ -0,0 +1,297 @@ +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] + /// 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$. + fn lp_norm(&self, x: &[f64]) -> f64 { + x.iter() + .map(|xi| xi.abs().powf(self.p)) + .sum::() + .powf(1.0 / self.p) + } + + fn project_lp_ball(&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); + 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_lp_ball(&mut shifted); + + x.iter_mut() + .zip(shifted.iter().zip(center.iter())) + .for_each(|(xi, (si, ci))| *xi = *ci + *si); + } else { + self.project_lp_ball(x); + } + } + + fn is_convex(&self) -> bool { + true + } +} diff --git a/src/constraints/mod.rs b/src/constraints/mod.rs index dd0b7038..d4ce2462 100644 --- a/src/constraints/mod.rs +++ b/src/constraints/mod.rs @@ -12,6 +12,7 @@ mod affine_space; mod ball1; mod ball2; mod ballinf; +mod ballp; mod cartesian_product; mod epigraph_squared_norm; mod finite; @@ -28,6 +29,7 @@ pub use affine_space::AffineSpace; pub use ball1::Ball1; pub use ball2::Ball2; 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 26343688..0bb32c43 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() { @@ -989,3 +991,189 @@ fn t_affine_space_wrong_dimensions() { let b = vec![1., 2., -0.5]; let _ = AffineSpace::new(a, b); } + +/// Sample `n_points` random points on the `l_p` sphere in `R^dim` +/// +/// 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_sample_lp_sphere_points_have_correct_norm() { + let n_points = 10_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. +/// +/// 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], + p: f64, + radius: f64, + sample_points: usize, +) -> 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 inner_prod_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() + .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 < inner_prod_tol { + 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-16; + let max_iters: usize = 200; + let ball = BallP::new(None, radius, p, tol, max_iters); + ball.project(&mut x); + 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); + unit_test_utils::assert_nearly_equal_array( + &x0, + &x, + 1e-12, + 1e-12, + "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", + ); +}