From 7e9b6261e045e5e841cfaa9db1dc56bd7b554aea Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Fri, 13 Mar 2026 13:00:51 +0000 Subject: [PATCH 01/12] New constraint (Rust): sq. norm epigraph About: - Implementaiton of epigraph of sq. norm as constraint - Basic unit test - Documentation --- src/constraints/mod.rs | 2 + src/constraints/squared_norm_epigraph.rs | 254 +++++++++++++++++++++++ src/constraints/tests.rs | 13 ++ 3 files changed, 269 insertions(+) create mode 100644 src/constraints/squared_norm_epigraph.rs diff --git a/src/constraints/mod.rs b/src/constraints/mod.rs index dd0b7038..8b29ef54 100644 --- a/src/constraints/mod.rs +++ b/src/constraints/mod.rs @@ -22,6 +22,7 @@ mod rectangle; mod simplex; mod soc; mod sphere2; +mod squared_norm_epigraph; mod zero; pub use affine_space::AffineSpace; @@ -38,6 +39,7 @@ pub use rectangle::Rectangle; pub use simplex::Simplex; pub use soc::SecondOrderCone; pub use sphere2::Sphere2; +pub use squared_norm_epigraph::EpiSqNorm; pub use zero::Zero; /// A set which can be used as a constraint diff --git a/src/constraints/squared_norm_epigraph.rs b/src/constraints/squared_norm_epigraph.rs new file mode 100644 index 00000000..7572137d --- /dev/null +++ b/src/constraints/squared_norm_epigraph.rs @@ -0,0 +1,254 @@ +use super::Constraint; + +#[derive(Copy, Clone)] +/// Epigraph of the squared Euclidean norm, that is, +/// $$\mathrm{epi}\Vert \cdot \Vert^2 = \\{(x, z)\in\mathbb{R}^{n+1}: \Vert x \Vert^2 \leq z \\}$$ +/// +/// A point $u = (x, z)$ is represented by a slice `u` of length at least 2: +/// +/// - `u[..u.len()-1]` is the vector part `x`, +/// - `u[u.len()-1]` is the scalar part `z`. +pub struct EpiSqNorm { + /// Tolerance used to accept a candidate root. + tol_root: f64, + + /// Tolerance used in Newton refinement. + tol_newton: f64, + + /// Maximum number of Newton iterations. + max_iter: usize, +} + +impl EpiSqNorm { + /// Constructor + /// + /// Creates a new instance of `EpiSqNorm` + pub fn new() -> Self { + Self { + tol_root: 1e-6, + tol_newton: 1e-10, + max_iter: 20, + } + } + + /// Set `tol_root` parameter + pub fn with_tol_root(mut self, tol_root: f64) -> Self { + self.tol_root = tol_root; + self + } + + /// Set `tol_newton` parameter + pub fn with_tol_newton(mut self, tol_newton: f64) -> Self { + self.tol_newton = tol_newton; + self + } + + /// Set `tol_newton` parameter + pub fn with_max_iter(mut self, max_iter: usize) -> Self { + self.max_iter = max_iter; + self + } + + #[inline] + fn norm2_squared(x: &[f64]) -> f64 { + x.iter().map(|xi| xi * xi).sum() + } + + /// Returns the real roots of $a r^3 + b r^2 + c r + d = 0$. + fn solve_cubic_real_roots(a: f64, b: f64, c: f64, d: f64) -> Vec { + assert!(a != 0.0); + + let p = (3.0 * a * c - b * b) / (3.0 * a * a); + let q = (2.0 * b.powi(3) - 9.0 * a * b * c + 27.0 * a * a * d) / (27.0 * a * a * a); + let shift = -b / (3.0 * a); + let delta = q * q / 4.0 + p * p * p / 27.0; + + if delta > 0.0 { + let sqrt_delta = delta.sqrt(); + let u = (-q / 2.0 + sqrt_delta).cbrt(); + let v = (-q / 2.0 - sqrt_delta).cbrt(); + vec![u + v + shift] + } else if delta.abs() <= 1e-14 { + let u = (-q / 2.0).cbrt(); + vec![2.0 * u + shift, -u + shift] + } else { + let rho = (-p * p * p / 27.0).sqrt(); + let phi = (-q / (2.0 * rho)).acos(); + let m = 2.0 * (-p / 3.0).sqrt(); + + let r1 = m * (phi / 3.0).cos() + shift; + let r2 = m * ((phi + 2.0 * std::f64::consts::PI) / 3.0).cos() + shift; + let r3 = m * ((phi + 4.0 * std::f64::consts::PI) / 3.0).cos() + shift; + + vec![r1, r2, r3] + } + } + + /// Computes candidate roots of + /// $$4 r^3 + 4 θ r^2 + θ^2 r - \Vert x \Vert^2 = 0.$$ + fn cubic_roots(theta: f64, x_norm_sq: f64) -> Vec { + let b = 4.0 * theta; + let c = theta * theta; + let d = -x_norm_sq; + + let discr = 72.0 * b * c * d - 4.0 * b.powi(3) * d + b * b * c * c + - 16.0 * c.powi(3) + - 432.0 * d * d; + let discr0 = b * b - 12.0 * c; + + if discr.abs() <= 1e-14 { + if discr0.abs() <= 1e-14 { + return vec![-b / 12.0]; + } else { + let single = (16.0 * b * c - 144.0 * d - b.powi(3)) / (4.0 * discr0); + let double = (36.0 * d - b * c) / (2.0 * discr0); + return vec![single, double]; + } + } + Self::solve_cubic_real_roots(4.0, b, c, d) + } + + /// Newton refinement for + /// 4 z^3 + 4 θ z^2 + θ^2 z - ||x||² = 0. + fn newton_solver(x_norm_sq: f64, theta: f64, z0: f64, max_iter: usize, tolx: f64) -> f64 { + let mut zsol = z0; + let mut zprev = z0 - 1.0; + + for _ in 0..max_iter { + let numerator = + 4.0 * zsol.powi(3) + 4.0 * theta * zsol.powi(2) + theta * theta * zsol - x_norm_sq; + let denominator = 12.0 * zsol.powi(2) + 8.0 * theta * zsol + theta * theta; + + let next = zsol - numerator / denominator; + let err = (next - zprev).abs(); + + zprev = zsol; + zsol = next; + + if err < tolx { + return zsol; + } + } + + zsol + } +} + +impl Constraint for EpiSqNorm { + /// Project a point onto the epigraph of the squared Euclidean norm. + /// + /// A point $(y, t)$ is represented by a slice `x` of length at least `2`: + /// + /// - `x[..x.len()-1]` is the vector part $y$, + /// - `x[x.len()-1]` is the scalar part $t$. + /// + /// The set is + /// $$C = \\{ (u, s) \in \mathbb{R}^n \times \mathbb{R} : \|u\|_2^2 \le s \\}.$$ + /// This method computes the Euclidean projection of $(y, t)$ onto $C$. + /// If the input is already feasible, that is, if $\Vert y \Vert_2^2 \leq t$, + /// then the input is left unchanged. + /// + /// # Method + /// + /// For an infeasible point $(y, t)$, the projection $(u, s)$ satisfies the + /// optimality conditions $y - u = 2 (s - t) u$, and $s = \Vert u\Vert_2^2$. + /// Solving for $u$, we have + /// $$u = \frac{y}{1 + 2(s - t)}.$$ + /// Substituting this into $s = \Vert u \Vert_2^2$ yields the scalar equation + /// $$s \left(1 + 2(s - t)\right)^2 = \Vert y \Vert_2^2.$$ + /// Defining $\theta = 1 - 2t$, this becomes the cubic polynomial + /// $$4s^3 + 4\theta s^2 + \theta^2 s - \|y\|_2^2 = 0.$$ + /// + /// The implementation follows this procedure: + /// + /// 1. Check whether the point is already in the epigraph. + /// 2. Form the cubic equation in the projected scalar variable $s$. + /// 3. Compute its real candidate roots. + /// 4. Select a root that is consistent with + /// + /// $$ + /// u = \frac{y}{1 + 2(s - t)} + /// $$ + /// + /// and + /// + /// $$ + /// s \approx \|u\|_2^2. + /// $$ + /// + /// 5. Refine the selected root with a few Newton iterations applied to + /// + /// $$ + /// f(s) = 4s^3 + 4\theta s^2 + \theta^2 s - \Vert y \Vert_2^2. + /// $$ + /// + /// 6. Recover the projected vector component from + /// + /// $$ + /// u = \frac{y}{1 + 2(s - t)}. + /// $$ + /// + /// # Notes + /// + /// - The projection is taken with respect to the standard Euclidean norm on + /// the product space. + /// - The last entry of `x` is interpreted as the epigraph variable. + /// - The set is convex, so the projection is uniquely defined. + /// + /// # Panics + /// + /// Panics if: + /// + /// - `x.len() < 2`, + /// - no admissible real root is found for the cubic equation. + fn project(&self, x: &mut [f64]) { + assert!( + x.len() >= 2, + "EpiSqNorm projection requires a slice of length at least 2" + ); + + let n = x.len(); + let z = x[n - 1]; + let x_vec = &x[..n - 1]; + let x_norm_sq = Self::norm2_squared(x_vec); + + // Already feasible + if x_norm_sq <= z { + return; + } + + let theta = 1.0 - 2.0 * z; + let roots = Self::cubic_roots(theta, x_norm_sq); + + let mut z_proj = None; + + for &r in roots.iter() { + let denom = 1.0 + 2.0 * (r - z); + if denom.abs() <= 1e-15 { + continue; + } + + let candidate_norm_sq = x_vec.iter().map(|xi| (xi / denom).powi(2)).sum::(); + + if (candidate_norm_sq - r).abs() <= self.tol_root { + z_proj = Some(r); + break; + } + } + + let mut z_sol = z_proj.expect("No admissible real root found in EpiSqNorm::project"); + + z_sol = Self::newton_solver(x_norm_sq, theta, z_sol, self.max_iter, self.tol_newton); + + let denom = 1.0 + 2.0 * (z_sol - z); + + x[..n - 1].iter_mut().for_each(|xi| *xi /= denom); + x[n - 1] = z_sol; + } + + /// The epigraph of the squared Euclidean norm is convex, + /// so this returns `true`. + fn is_convex(&self) -> bool { + true + } +} diff --git a/src/constraints/tests.rs b/src/constraints/tests.rs index 26343688..26b24fad 100644 --- a/src/constraints/tests.rs +++ b/src/constraints/tests.rs @@ -989,3 +989,16 @@ fn t_affine_space_wrong_dimensions() { let b = vec![1., 2., -0.5]; let _ = AffineSpace::new(a, b); } + +#[test] +fn t_epigraph_squared_norm_basic() { + let epi = EpiSqNorm::new() + .with_max_iter(20) + .with_tol_root(1e-6) + .with_tol_newton(1e-10); + let mut u = [2.0, -1.0, 0.5]; + epi.project(&mut u); + let x_norm_sq = u[..2].iter().map(|xi| xi * xi).sum::(); + let z = u[2]; + assert!(x_norm_sq <= z + 1e-8); +} From 7b09cb172ebc45fe364b7bb9fac7a7bbb4cfd651 Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Fri, 13 Mar 2026 13:24:12 +0000 Subject: [PATCH 02/12] Epigraph of sq. Eucl. norm: improvements About: - algorithmic improvements for better num. stability - better documentation --- src/constraints/epigraph_squared_norm.rs | 141 ++++++++++--- src/constraints/mod.rs | 2 - src/constraints/squared_norm_epigraph.rs | 254 ----------------------- src/constraints/tests.rs | 13 -- 4 files changed, 107 insertions(+), 303 deletions(-) delete mode 100644 src/constraints/squared_norm_epigraph.rs diff --git a/src/constraints/epigraph_squared_norm.rs b/src/constraints/epigraph_squared_norm.rs index 560b34bd..77d3b0ff 100644 --- a/src/constraints/epigraph_squared_norm.rs +++ b/src/constraints/epigraph_squared_norm.rs @@ -3,27 +3,69 @@ 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 \\}.$ +/// The epigraph of the squared Euclidean norm, that is, +/// $$ +/// X = \{x = (z, t) \in \mathbb{R}^{n}\times \mathbb{R} : \|z\|_2^2 \le t \}. +/// $$ +/// +/// A point is represented by a slice `x` whose last entry is the scalar +/// component `t`, while the preceding entries form the vector component `z`. 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. + #[must_use] pub fn new() -> Self { EpigraphSquaredNorm {} } } impl Constraint for EpigraphSquaredNorm { - ///Project on the epigraph of the squared Euclidean norm. + /// 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/). + /// Let the input be represented as $(z,t)$, where `z` is the vector formed + /// by the first `x.len() - 1` entries of `x`, and `t` is the last entry. + /// This method computes the Euclidean projection of $(z,t)$ onto + /// + /// $$ + /// \operatorname{epi}\|\cdot\|_2^2 + /// = + /// \{(u,s) \in \mathbb{R}^n \times \mathbb{R} : \|u\|_2^2 \le s\}. + /// $$ + /// + /// If the point is already feasible, that is, if + /// + /// $$ + /// \|z\|_2^2 \le t, + /// $$ + /// + /// then the input is left unchanged. + /// + /// Otherwise, the projection is computed using the methodology described + /// [here](https://mathematix.wordpress.com/2017/05/02/projection-on-the-epigraph-of-the-squared-euclidean-norm/): + /// + /// 1. Eliminate the projected vector variable to obtain a cubic equation + /// in the projected scalar variable. + /// 2. Compute the real roots of that cubic. + /// 3. Select the admissible root that is consistent with the projection + /// formula. + /// 4. Refine the selected root with a few Newton iterations. + /// 5. Recover the projected vector component using the refined root. /// /// ## Arguments - /// - `x`: The given vector $x$ is updated with the projection on the set + /// + /// - `x`: The given vector `x` is updated with the projection on the set + /// + /// ## Panics + /// + /// Panics if: + /// + /// - `x.len() < 2`, + /// - no admissible real root is found, + /// - the Newton derivative becomes too small, + /// - the final scaling factor is numerically singular. /// /// ## Example /// @@ -31,56 +73,87 @@ impl Constraint for EpigraphSquaredNorm { /// use optimization_engine::constraints::*; /// /// let epi = EpigraphSquaredNorm::new(); - /// let mut x = [1., 2., 3., 4.]; + /// + /// // Here, z = [1., 2., 3.] and t = 4. + /// let mut x = [1., 2., 3., 4.]; + /// /// epi.project(&mut x); /// ``` fn project(&self, x: &mut [f64]) { + assert!( + x.len() >= 2, + "EpigraphSquaredNorm::project requires x.len() >= 2" + ); + 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 z = &x[..nx]; + let t = x[nx]; let norm_z_sq = matrix_operations::norm2_squared(z); + + // Already feasible if norm_z_sq <= t { return; } - let theta = 1. - 2. * t; - let a3 = 4.; - let a2 = 4. * theta; + // Cubic: + // 4 r^3 + 4 theta r^2 + theta^2 r - ||z||^2 = 0 + let theta = 1.0 - 2.0 * t; + let a3 = 4.0; + let a2 = 4.0 * 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; + + let root_tol = 1e-6; + let mut right_root: Option = None; + + // Pick the first admissible real root + for &ri in cubic_poly_roots.as_ref().iter() { + let denom = 1.0 + 2.0 * (ri - t); + + // We need a valid scaling and consistency with ||z_proj||^2 = ri + if denom > 0.0 { + let candidate_norm_sq = norm_z_sq / (denom * denom); + if (candidate_norm_sq - ri).abs() <= root_tol { + right_root = Some(ri); + break; } } - }); + } + + let mut zsol = + right_root.expect("EpigraphSquaredNorm::project: no admissible real root found"); - // Refinement of root with Newton-Raphson - let mut refinement_error = 1.; + // Newton refinement 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 { + + for _ in 0..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; + if p_z.abs() <= newton_eps { + break; + } + + let dp_z = 3.0 * a3 * zsol_sq + 2.0 * a2 * zsol + a1; + assert!( + dp_z.abs() > 1e-15, + "EpigraphSquaredNorm::project: Newton derivative too small" + ); + zsol -= p_z / dp_z; - refinement_error = p_z.abs(); - iter += 1; } - right_root = zsol; + + let right_root = zsol; + let scaling = 1.0 + 2.0 * (right_root - t); + + assert!( + scaling.abs() > 1e-15, + "EpigraphSquaredNorm::project: scaling factor too small" + ); // Projection for xi in x.iter_mut().take(nx) { @@ -89,7 +162,7 @@ impl Constraint for EpigraphSquaredNorm { x[nx] = right_root; } - /// This is a convex set, so this function returns `True` + /// 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 8b29ef54..dd0b7038 100644 --- a/src/constraints/mod.rs +++ b/src/constraints/mod.rs @@ -22,7 +22,6 @@ mod rectangle; mod simplex; mod soc; mod sphere2; -mod squared_norm_epigraph; mod zero; pub use affine_space::AffineSpace; @@ -39,7 +38,6 @@ pub use rectangle::Rectangle; pub use simplex::Simplex; pub use soc::SecondOrderCone; pub use sphere2::Sphere2; -pub use squared_norm_epigraph::EpiSqNorm; pub use zero::Zero; /// A set which can be used as a constraint diff --git a/src/constraints/squared_norm_epigraph.rs b/src/constraints/squared_norm_epigraph.rs deleted file mode 100644 index 7572137d..00000000 --- a/src/constraints/squared_norm_epigraph.rs +++ /dev/null @@ -1,254 +0,0 @@ -use super::Constraint; - -#[derive(Copy, Clone)] -/// Epigraph of the squared Euclidean norm, that is, -/// $$\mathrm{epi}\Vert \cdot \Vert^2 = \\{(x, z)\in\mathbb{R}^{n+1}: \Vert x \Vert^2 \leq z \\}$$ -/// -/// A point $u = (x, z)$ is represented by a slice `u` of length at least 2: -/// -/// - `u[..u.len()-1]` is the vector part `x`, -/// - `u[u.len()-1]` is the scalar part `z`. -pub struct EpiSqNorm { - /// Tolerance used to accept a candidate root. - tol_root: f64, - - /// Tolerance used in Newton refinement. - tol_newton: f64, - - /// Maximum number of Newton iterations. - max_iter: usize, -} - -impl EpiSqNorm { - /// Constructor - /// - /// Creates a new instance of `EpiSqNorm` - pub fn new() -> Self { - Self { - tol_root: 1e-6, - tol_newton: 1e-10, - max_iter: 20, - } - } - - /// Set `tol_root` parameter - pub fn with_tol_root(mut self, tol_root: f64) -> Self { - self.tol_root = tol_root; - self - } - - /// Set `tol_newton` parameter - pub fn with_tol_newton(mut self, tol_newton: f64) -> Self { - self.tol_newton = tol_newton; - self - } - - /// Set `tol_newton` parameter - pub fn with_max_iter(mut self, max_iter: usize) -> Self { - self.max_iter = max_iter; - self - } - - #[inline] - fn norm2_squared(x: &[f64]) -> f64 { - x.iter().map(|xi| xi * xi).sum() - } - - /// Returns the real roots of $a r^3 + b r^2 + c r + d = 0$. - fn solve_cubic_real_roots(a: f64, b: f64, c: f64, d: f64) -> Vec { - assert!(a != 0.0); - - let p = (3.0 * a * c - b * b) / (3.0 * a * a); - let q = (2.0 * b.powi(3) - 9.0 * a * b * c + 27.0 * a * a * d) / (27.0 * a * a * a); - let shift = -b / (3.0 * a); - let delta = q * q / 4.0 + p * p * p / 27.0; - - if delta > 0.0 { - let sqrt_delta = delta.sqrt(); - let u = (-q / 2.0 + sqrt_delta).cbrt(); - let v = (-q / 2.0 - sqrt_delta).cbrt(); - vec![u + v + shift] - } else if delta.abs() <= 1e-14 { - let u = (-q / 2.0).cbrt(); - vec![2.0 * u + shift, -u + shift] - } else { - let rho = (-p * p * p / 27.0).sqrt(); - let phi = (-q / (2.0 * rho)).acos(); - let m = 2.0 * (-p / 3.0).sqrt(); - - let r1 = m * (phi / 3.0).cos() + shift; - let r2 = m * ((phi + 2.0 * std::f64::consts::PI) / 3.0).cos() + shift; - let r3 = m * ((phi + 4.0 * std::f64::consts::PI) / 3.0).cos() + shift; - - vec![r1, r2, r3] - } - } - - /// Computes candidate roots of - /// $$4 r^3 + 4 θ r^2 + θ^2 r - \Vert x \Vert^2 = 0.$$ - fn cubic_roots(theta: f64, x_norm_sq: f64) -> Vec { - let b = 4.0 * theta; - let c = theta * theta; - let d = -x_norm_sq; - - let discr = 72.0 * b * c * d - 4.0 * b.powi(3) * d + b * b * c * c - - 16.0 * c.powi(3) - - 432.0 * d * d; - let discr0 = b * b - 12.0 * c; - - if discr.abs() <= 1e-14 { - if discr0.abs() <= 1e-14 { - return vec![-b / 12.0]; - } else { - let single = (16.0 * b * c - 144.0 * d - b.powi(3)) / (4.0 * discr0); - let double = (36.0 * d - b * c) / (2.0 * discr0); - return vec![single, double]; - } - } - Self::solve_cubic_real_roots(4.0, b, c, d) - } - - /// Newton refinement for - /// 4 z^3 + 4 θ z^2 + θ^2 z - ||x||² = 0. - fn newton_solver(x_norm_sq: f64, theta: f64, z0: f64, max_iter: usize, tolx: f64) -> f64 { - let mut zsol = z0; - let mut zprev = z0 - 1.0; - - for _ in 0..max_iter { - let numerator = - 4.0 * zsol.powi(3) + 4.0 * theta * zsol.powi(2) + theta * theta * zsol - x_norm_sq; - let denominator = 12.0 * zsol.powi(2) + 8.0 * theta * zsol + theta * theta; - - let next = zsol - numerator / denominator; - let err = (next - zprev).abs(); - - zprev = zsol; - zsol = next; - - if err < tolx { - return zsol; - } - } - - zsol - } -} - -impl Constraint for EpiSqNorm { - /// Project a point onto the epigraph of the squared Euclidean norm. - /// - /// A point $(y, t)$ is represented by a slice `x` of length at least `2`: - /// - /// - `x[..x.len()-1]` is the vector part $y$, - /// - `x[x.len()-1]` is the scalar part $t$. - /// - /// The set is - /// $$C = \\{ (u, s) \in \mathbb{R}^n \times \mathbb{R} : \|u\|_2^2 \le s \\}.$$ - /// This method computes the Euclidean projection of $(y, t)$ onto $C$. - /// If the input is already feasible, that is, if $\Vert y \Vert_2^2 \leq t$, - /// then the input is left unchanged. - /// - /// # Method - /// - /// For an infeasible point $(y, t)$, the projection $(u, s)$ satisfies the - /// optimality conditions $y - u = 2 (s - t) u$, and $s = \Vert u\Vert_2^2$. - /// Solving for $u$, we have - /// $$u = \frac{y}{1 + 2(s - t)}.$$ - /// Substituting this into $s = \Vert u \Vert_2^2$ yields the scalar equation - /// $$s \left(1 + 2(s - t)\right)^2 = \Vert y \Vert_2^2.$$ - /// Defining $\theta = 1 - 2t$, this becomes the cubic polynomial - /// $$4s^3 + 4\theta s^2 + \theta^2 s - \|y\|_2^2 = 0.$$ - /// - /// The implementation follows this procedure: - /// - /// 1. Check whether the point is already in the epigraph. - /// 2. Form the cubic equation in the projected scalar variable $s$. - /// 3. Compute its real candidate roots. - /// 4. Select a root that is consistent with - /// - /// $$ - /// u = \frac{y}{1 + 2(s - t)} - /// $$ - /// - /// and - /// - /// $$ - /// s \approx \|u\|_2^2. - /// $$ - /// - /// 5. Refine the selected root with a few Newton iterations applied to - /// - /// $$ - /// f(s) = 4s^3 + 4\theta s^2 + \theta^2 s - \Vert y \Vert_2^2. - /// $$ - /// - /// 6. Recover the projected vector component from - /// - /// $$ - /// u = \frac{y}{1 + 2(s - t)}. - /// $$ - /// - /// # Notes - /// - /// - The projection is taken with respect to the standard Euclidean norm on - /// the product space. - /// - The last entry of `x` is interpreted as the epigraph variable. - /// - The set is convex, so the projection is uniquely defined. - /// - /// # Panics - /// - /// Panics if: - /// - /// - `x.len() < 2`, - /// - no admissible real root is found for the cubic equation. - fn project(&self, x: &mut [f64]) { - assert!( - x.len() >= 2, - "EpiSqNorm projection requires a slice of length at least 2" - ); - - let n = x.len(); - let z = x[n - 1]; - let x_vec = &x[..n - 1]; - let x_norm_sq = Self::norm2_squared(x_vec); - - // Already feasible - if x_norm_sq <= z { - return; - } - - let theta = 1.0 - 2.0 * z; - let roots = Self::cubic_roots(theta, x_norm_sq); - - let mut z_proj = None; - - for &r in roots.iter() { - let denom = 1.0 + 2.0 * (r - z); - if denom.abs() <= 1e-15 { - continue; - } - - let candidate_norm_sq = x_vec.iter().map(|xi| (xi / denom).powi(2)).sum::(); - - if (candidate_norm_sq - r).abs() <= self.tol_root { - z_proj = Some(r); - break; - } - } - - let mut z_sol = z_proj.expect("No admissible real root found in EpiSqNorm::project"); - - z_sol = Self::newton_solver(x_norm_sq, theta, z_sol, self.max_iter, self.tol_newton); - - let denom = 1.0 + 2.0 * (z_sol - z); - - x[..n - 1].iter_mut().for_each(|xi| *xi /= denom); - x[n - 1] = z_sol; - } - - /// The epigraph of the squared Euclidean norm is convex, - /// so this returns `true`. - fn is_convex(&self) -> bool { - true - } -} diff --git a/src/constraints/tests.rs b/src/constraints/tests.rs index 26b24fad..26343688 100644 --- a/src/constraints/tests.rs +++ b/src/constraints/tests.rs @@ -989,16 +989,3 @@ fn t_affine_space_wrong_dimensions() { let b = vec![1., 2., -0.5]; let _ = AffineSpace::new(a, b); } - -#[test] -fn t_epigraph_squared_norm_basic() { - let epi = EpiSqNorm::new() - .with_max_iter(20) - .with_tol_root(1e-6) - .with_tol_newton(1e-10); - let mut u = [2.0, -1.0, 0.5]; - epi.project(&mut u); - let x_norm_sq = u[..2].iter().map(|xi| xi * xi).sum::(); - let z = u[2]; - assert!(x_norm_sq <= z + 1e-8); -} From bec77e9eb74d1c06f9b967786d0fc36d3aa626e1 Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Fri, 13 Mar 2026 13:40:12 +0000 Subject: [PATCH 03/12] [ci skip] EpigraphSquaredNorm: improved docs --- src/constraints/epigraph_squared_norm.rs | 28 +++++------------------- 1 file changed, 6 insertions(+), 22 deletions(-) diff --git a/src/constraints/epigraph_squared_norm.rs b/src/constraints/epigraph_squared_norm.rs index 77d3b0ff..f91617c8 100644 --- a/src/constraints/epigraph_squared_norm.rs +++ b/src/constraints/epigraph_squared_norm.rs @@ -5,7 +5,7 @@ use super::Constraint; #[derive(Copy, Clone, Default)] /// The epigraph of the squared Euclidean norm, that is, /// $$ -/// X = \{x = (z, t) \in \mathbb{R}^{n}\times \mathbb{R} : \|z\|_2^2 \le t \}. +/// X = \\{x = (z, t) \in \mathbb{R}^{n}\times \mathbb{R} : \Vert z\Vert_2^2 \leq t \\}. /// $$ /// /// A point is represented by a slice `x` whose last entry is the scalar @@ -30,29 +30,13 @@ impl Constraint for EpigraphSquaredNorm { /// This method computes the Euclidean projection of $(z,t)$ onto /// /// $$ - /// \operatorname{epi}\|\cdot\|_2^2 - /// = - /// \{(u,s) \in \mathbb{R}^n \times \mathbb{R} : \|u\|_2^2 \le s\}. + /// \mathrm{epi}\Vert \cdot\Vert_2^2 = \\{(u,s) \in \mathbb{R}^n \times \mathbb{R} : \Vert u\Vert_2^2 \leq s \\}. /// $$ /// - /// If the point is already feasible, that is, if - /// - /// $$ - /// \|z\|_2^2 \le t, - /// $$ - /// - /// then the input is left unchanged. - /// - /// Otherwise, the projection is computed using the methodology described - /// [here](https://mathematix.wordpress.com/2017/05/02/projection-on-the-epigraph-of-the-squared-euclidean-norm/): - /// - /// 1. Eliminate the projected vector variable to obtain a cubic equation - /// in the projected scalar variable. - /// 2. Compute the real roots of that cubic. - /// 3. Select the admissible root that is consistent with the projection - /// formula. - /// 4. Refine the selected root with a few Newton iterations. - /// 5. Recover the projected vector component using the refined root. + /// If the point is already feasible, that is, if $\Vert z\Vert_2^2 \leq t,$ + /// then the input is left unchanged, otherwise, the projection is computed using + /// the methodology described + /// [here](https://mathematix.wordpress.com/2017/05/02/projection-on-the-epigraph-of-the-squared-euclidean-norm/). /// /// ## Arguments /// From 40087b2f88943d9fa646cf914b64c10c865cc22e Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Sat, 14 Mar 2026 13:10:09 +0000 Subject: [PATCH 04/12] [ci skip] update CHANGELOG --- CHANGELOG.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 385e576e..3d615814 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -16,6 +16,10 @@ Note: This is the main Changelog file for the Rust solver. The Changelog file fo - Implementation of `BallP` in Rust: projection on lp-ball +### Changed + +- Algorithmic improvements in `EpigraphSquaredNorm` (numerically stable Newton refinement) and more detailed docs +