|
| 1 | +use super::Constraint; |
| 2 | + |
| 3 | +#[derive(Copy, Clone)] |
| 4 | +/// An $\\ell_p$ ball, that is, |
| 5 | +/// $$B_p^r = \\{ x \\in R^n : \Vert x \Vert_p \\leq r \\}$$ |
| 6 | +/// or a translated ball |
| 7 | +/// $$B_p^{x_c, r} = \\{ x \\in \mathbb{R}^n : \Vert x - x_c \Vert_p \\leq r \\},$$ |
| 8 | +/// with $1 < p < \infty$. |
| 9 | +/// |
| 10 | +/// # Projection problem |
| 11 | +/// |
| 12 | +/// Given a vector $x$, projection onto the ball means solving |
| 13 | +/// |
| 14 | +/// $$\Pi_{B_p^r}(x) = \mathrm{argmin}_{z \in B_p^r} \frac{1}{2}\Vert z - x \Vert_2^2.$$ |
| 15 | +/// |
| 16 | +/// If $x$ already belongs to the ball, the projection is $x$ itself. |
| 17 | +/// Otherwise, the projection lies on the boundary of the ball. |
| 18 | +/// |
| 19 | +/// # Numerical method |
| 20 | +/// |
| 21 | +/// For general $1 < p < \infty$, the projection does not admit a simple |
| 22 | +/// closed-form expression. This implementation computes it numerically using |
| 23 | +/// the optimality conditions of the constrained problem. |
| 24 | +/// |
| 25 | +/// For a ball centered at the origin, the projection is characterized by |
| 26 | +/// |
| 27 | +/// $$z_i = \operatorname{sign}(x_i)\,u_i(\lambda),$$ |
| 28 | +/// |
| 29 | +/// where each scalar $u_i(\lambda) \geq 0$ solves |
| 30 | +/// |
| 31 | +/// $$u_i + \lambda p\,u_i^{p-1} = |x_i|,$$ |
| 32 | +/// |
| 33 | +/// and the Lagrange multiplier $\lambda \geq 0$ is chosen so that the projected |
| 34 | +/// vector satisfies the feasibility condition $\|z\|_p = r.$ |
| 35 | +/// |
| 36 | +/// The implementation uses: |
| 37 | +/// |
| 38 | +/// - an *outer bisection loop* to determine the multiplier $\lambda$, |
| 39 | +/// - an *inner safeguarded Newton method* to solve the scalar nonlinear |
| 40 | +/// equation for each coordinate. |
| 41 | +/// |
| 42 | +/// The Newton iteration is combined with bracketing and a bisection fallback, |
| 43 | +/// which improves robustness. |
| 44 | +/// |
| 45 | +/// # Translated balls |
| 46 | +/// |
| 47 | +/// If the ball has a center $x_c$, projection is performed by translating the |
| 48 | +/// input vector to the origin, projecting onto the origin-centered ball, and |
| 49 | +/// translating the result back: |
| 50 | +/// |
| 51 | +/// $$\Pi_{B_p^{x_c, r}}(x) = x_c + \Pi_{B_p^r}(x - x_c).$$ |
| 52 | +/// |
| 53 | +/// # Convexity |
| 54 | +/// |
| 55 | +/// Since this module assumes `p > 1.0`, every [`BallP`] is convex, and therefore |
| 56 | +/// [`Constraint::project`] computes a unique Euclidean projection. |
| 57 | +/// |
| 58 | +/// # Examples |
| 59 | +/// |
| 60 | +/// Project onto the unit \(\ell_p\)-ball centered at the origin: |
| 61 | +/// |
| 62 | +/// ``` |
| 63 | +/// use optimization_engine::constraints::{BallP, Constraint}; |
| 64 | +/// |
| 65 | +/// let ball = BallP::new(None, 1.0, 1.5, 1e-10, 100); |
| 66 | +/// let mut x = vec![3.0, -1.0, 2.0]; |
| 67 | +/// ball.project(&mut x); |
| 68 | +/// ``` |
| 69 | +/// |
| 70 | +/// Project onto a translated \(\ell_p\)-ball: |
| 71 | +/// |
| 72 | +/// ``` |
| 73 | +/// use optimization_engine::constraints::{BallP, Constraint}; |
| 74 | +/// |
| 75 | +/// let center = vec![1.0, 1.0, 1.0]; |
| 76 | +/// let ball = BallP::new(Some(¢er), 2.0, 3.0, 1e-10, 100); |
| 77 | +/// let mut x = vec![4.0, -1.0, 2.0]; |
| 78 | +/// ball.project(&mut x); |
| 79 | +/// ``` |
| 80 | +/// |
| 81 | +/// # Notes |
| 82 | +/// |
| 83 | +/// - The projection is with respect to the *Euclidean norm* |
| 84 | +/// - The implementation is intended for general finite $p > 1.0$. If you need |
| 85 | +/// to project on a $\Vert{}\cdot{}\Vert_1$-ball or an $\Vert{}\cdot{}\Vert_\infty$-ball, |
| 86 | +/// use the implementations in [`Ball1`](crate::constraints::Ball1) |
| 87 | +/// and [`BallInf`](crate::constraints::BallInf). |
| 88 | +/// - Do not use this struct to project on a Euclidean ball; the implementation |
| 89 | +/// in [`Ball2`](crate::constraints::Ball2) is more efficient |
| 90 | +/// - The quality and speed of the computation depend on the chosen numerical |
| 91 | +/// tolerance and iteration limit. |
| 92 | +pub struct BallP<'a> { |
| 93 | + /// Optional center of the ball. |
| 94 | + /// |
| 95 | + /// If `None`, the ball is centered at the origin. |
| 96 | + /// If `Some(center)`, the ball is centered at `center`. |
| 97 | + center: Option<&'a [f64]>, |
| 98 | + |
| 99 | + /// Radius of the ball. |
| 100 | + /// |
| 101 | + /// Must be strictly positive. |
| 102 | + radius: f64, |
| 103 | + |
| 104 | + /// Exponent of the norm. |
| 105 | + /// |
| 106 | + /// Must satisfy `p > 1.0` and be finite. |
| 107 | + p: f64, |
| 108 | + |
| 109 | + /// Numerical tolerance used by the outer bisection on the Lagrange |
| 110 | + /// multiplier and by the inner Newton solver. |
| 111 | + tolerance: f64, |
| 112 | + |
| 113 | + /// Maximum number of iterations used by the outer bisection and |
| 114 | + /// the inner Newton solver. |
| 115 | + max_iter: usize, |
| 116 | +} |
| 117 | + |
| 118 | +impl<'a> BallP<'a> { |
| 119 | + /// Construct a new l_p ball with given center, radius, and exponent. |
| 120 | + /// |
| 121 | + /// - `center`: if `None`, the ball is centered at the origin |
| 122 | + /// - `radius`: radius of the ball |
| 123 | + /// - `p`: norm exponent, must satisfy `p > 1.0` and be finite |
| 124 | + /// - `tolerance`: tolerance for the numerical solvers |
| 125 | + /// - `max_iter`: maximum number of iterations for the numerical solvers |
| 126 | + pub fn new( |
| 127 | + center: Option<&'a [f64]>, |
| 128 | + radius: f64, |
| 129 | + p: f64, |
| 130 | + tolerance: f64, |
| 131 | + max_iter: usize, |
| 132 | + ) -> Self { |
| 133 | + assert!(radius > 0.0); |
| 134 | + assert!(p > 1.0 && p.is_finite()); |
| 135 | + assert!(tolerance > 0.0); |
| 136 | + assert!(max_iter > 0); |
| 137 | + |
| 138 | + BallP { |
| 139 | + center, |
| 140 | + radius, |
| 141 | + p, |
| 142 | + tolerance, |
| 143 | + max_iter, |
| 144 | + } |
| 145 | + } |
| 146 | + |
| 147 | + #[inline] |
| 148 | + /// Computes the $p$-norm of a given vector |
| 149 | + /// |
| 150 | + /// The $p$-norm of a vector $x\in \mathbb{R}^n$ is given by |
| 151 | + /// $$\Vert x \Vert_p = \left(\sum_{i=1}^{n} |x_i|^p\right)^{1/p},$$ |
| 152 | + /// for $p > 1$. |
| 153 | + fn lp_norm(&self, x: &[f64]) -> f64 { |
| 154 | + x.iter() |
| 155 | + .map(|xi| xi.abs().powf(self.p)) |
| 156 | + .sum::<f64>() |
| 157 | + .powf(1.0 / self.p) |
| 158 | + } |
| 159 | + |
| 160 | + fn project_lp_ball(&self, x: &mut [f64]) { |
| 161 | + let p = self.p; |
| 162 | + let r = self.radius; |
| 163 | + let tol = self.tolerance; |
| 164 | + let max_iter = self.max_iter; |
| 165 | + |
| 166 | + let current_norm = self.lp_norm(x); |
| 167 | + if current_norm <= r { |
| 168 | + return; |
| 169 | + } |
| 170 | + |
| 171 | + let abs_x: Vec<f64> = x.iter().map(|xi| xi.abs()).collect(); |
| 172 | + let target = r.powf(p); |
| 173 | + |
| 174 | + let radius_error = |lambda: f64| -> f64 { |
| 175 | + abs_x |
| 176 | + .iter() |
| 177 | + .map(|&a| { |
| 178 | + let u = Self::solve_coordinate_newton(a, lambda, p, tol, max_iter); |
| 179 | + u.powf(p) |
| 180 | + }) |
| 181 | + .sum::<f64>() |
| 182 | + - target |
| 183 | + }; |
| 184 | + |
| 185 | + let mut lambda_lo = 0.0_f64; |
| 186 | + let mut lambda_hi = 1.0_f64; |
| 187 | + |
| 188 | + while radius_error(lambda_hi) > 0.0 { |
| 189 | + lambda_hi *= 2.0; |
| 190 | + if lambda_hi > 1e20 { |
| 191 | + panic!("Failed to bracket the Lagrange multiplier"); |
| 192 | + } |
| 193 | + } |
| 194 | + |
| 195 | + for _ in 0..max_iter { |
| 196 | + let lambda_mid = 0.5 * (lambda_lo + lambda_hi); |
| 197 | + let err = radius_error(lambda_mid); |
| 198 | + |
| 199 | + if err.abs() <= tol { |
| 200 | + lambda_lo = lambda_mid; |
| 201 | + lambda_hi = lambda_mid; |
| 202 | + break; |
| 203 | + } |
| 204 | + |
| 205 | + if err > 0.0 { |
| 206 | + lambda_lo = lambda_mid; |
| 207 | + } else { |
| 208 | + lambda_hi = lambda_mid; |
| 209 | + } |
| 210 | + } |
| 211 | + |
| 212 | + let lambda_star = 0.5 * (lambda_lo + lambda_hi); |
| 213 | + |
| 214 | + x.iter_mut().zip(abs_x.iter()).for_each(|(xi, &a)| { |
| 215 | + let u = Self::solve_coordinate_newton(a, lambda_star, p, tol, max_iter); |
| 216 | + *xi = xi.signum() * u; |
| 217 | + }); |
| 218 | + } |
| 219 | + |
| 220 | + /// Solve for u >= 0 the equation u + lambda * p * u^(p-1) = a |
| 221 | + /// using safeguarded Newton iterations. |
| 222 | + /// |
| 223 | + /// The solution always belongs to [0, a], so Newton is combined with |
| 224 | + /// bracketing and a bisection fallback. |
| 225 | + fn solve_coordinate_newton(a: f64, lambda: f64, p: f64, tol: f64, max_iter: usize) -> f64 { |
| 226 | + if a == 0.0 { |
| 227 | + return 0.0; |
| 228 | + } |
| 229 | + |
| 230 | + if lambda == 0.0 { |
| 231 | + return a; |
| 232 | + } |
| 233 | + |
| 234 | + let mut lo = 0.0_f64; |
| 235 | + let mut hi = a; |
| 236 | + |
| 237 | + // Heuristic initial guess: |
| 238 | + // exact when p = 2, and usually in the right scale for general p. |
| 239 | + let mut u = (a / (1.0 + lambda * p)).clamp(lo, hi); |
| 240 | + |
| 241 | + for _ in 0..max_iter { |
| 242 | + let upm1 = u.powf(p - 1.0); |
| 243 | + let f = u + lambda * p * upm1 - a; |
| 244 | + |
| 245 | + if f.abs() <= tol { |
| 246 | + return u; |
| 247 | + } |
| 248 | + |
| 249 | + if f > 0.0 { |
| 250 | + hi = u; |
| 251 | + } else { |
| 252 | + lo = u; |
| 253 | + } |
| 254 | + |
| 255 | + let df = 1.0 + lambda * p * (p - 1.0) * u.powf(p - 2.0); |
| 256 | + let mut candidate = u - f / df; |
| 257 | + |
| 258 | + if !candidate.is_finite() || candidate <= lo || candidate >= hi { |
| 259 | + candidate = 0.5 * (lo + hi); |
| 260 | + } |
| 261 | + |
| 262 | + if (candidate - u).abs() <= tol * (1.0 + u.abs()) { |
| 263 | + return candidate; |
| 264 | + } |
| 265 | + |
| 266 | + u = candidate; |
| 267 | + } |
| 268 | + |
| 269 | + 0.5 * (lo + hi) |
| 270 | + } |
| 271 | +} |
| 272 | + |
| 273 | +impl<'a> Constraint for BallP<'a> { |
| 274 | + fn project(&self, x: &mut [f64]) { |
| 275 | + if let Some(center) = &self.center { |
| 276 | + assert_eq!(x.len(), center.len()); |
| 277 | + |
| 278 | + let mut shifted = vec![0.0; x.len()]; |
| 279 | + shifted |
| 280 | + .iter_mut() |
| 281 | + .zip(x.iter().zip(center.iter())) |
| 282 | + .for_each(|(s, (xi, ci))| *s = *xi - *ci); |
| 283 | + |
| 284 | + self.project_lp_ball(&mut shifted); |
| 285 | + |
| 286 | + x.iter_mut() |
| 287 | + .zip(shifted.iter().zip(center.iter())) |
| 288 | + .for_each(|(xi, (si, ci))| *xi = *ci + *si); |
| 289 | + } else { |
| 290 | + self.project_lp_ball(x); |
| 291 | + } |
| 292 | + } |
| 293 | + |
| 294 | + fn is_convex(&self) -> bool { |
| 295 | + true |
| 296 | + } |
| 297 | +} |
0 commit comments