Mercurial > repos > pointsource_algs / file diff

--- a/src/subproblem/nonneg.rs	Tue Aug 01 10:32:12 2023 +0300
+++ b/src/subproblem/nonneg.rs	Thu Aug 29 00:00:00 2024 -0500
@@ -6,6 +6,7 @@
 use numeric_literals::replace_float_literals;
 use itertools::{izip, Itertools};
 use colored::Colorize;
+use std::cmp::Ordering::*;
 
 use alg_tools::iter::Mappable;
 use alg_tools::error::NumericalError;
@@ -22,38 +23,42 @@
     InnerMethod,
     InnerSettings
 };
+use super::l1squared_nonneg::max_interval_dist_to_zero;
 
 /// Compute the proximal operator of $x \mapsto x + \delta\_{[0, \infty)}$, i.e.,
 /// the non-negativity contrained soft-thresholding operator.
 #[inline]
 #[replace_float_literals(F::cast_from(literal))]
-fn nonneg_soft_thresholding<F : Float>(v : F, λ : F) -> F {
+pub(super) fn nonneg_soft_thresholding<F : Float>(v : F, λ : F) -> F {
     (v - λ).max(0.0)
 }
 
-/// Returns the ∞-norm minimal subdifferential of $x ↦ x^⊤Ax - g^⊤ x + λ\vec 1^⊤ x δ_{≥ 0}(x)$
+/// Returns the ∞-norm minimal subdifferential of $x ↦ x^⊤Ax/2 - g^⊤ x + λ\vec 1^⊤ x + δ_{≥ 0}(x)$
 /// at $x$.
 ///
 /// `v` will be modified and cannot be trusted to contain useful values afterwards.
-#[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())]
-fn min_subdifferential<F : Float + ToNalgebraRealField>(
-    v : &mut DVector<F::MixedType>,
-    mA : &DMatrix<F::MixedType>,
-    x : &DVector<F::MixedType>,
-    g : &DVector<F::MixedType>,
-    λ : F::MixedType
+#[replace_float_literals(F::cast_from(literal))]
+fn min_subdifferential<F : Float + nalgebra::RealField>(
+    v : &mut DVector<F>,
+    mA : &DMatrix<F>,
+    x : &DVector<F>,
+    g : &DVector<F>,
+    λ : F
 ) -> F {
     v.copy_from(g);
     mA.gemv(v, 1.0, x, -1.0);   // v =  Ax - g
     let mut val = 0.0;
     for (&v_i, &x_i) in izip!(v.iter(), x.iter()) {
-        // The subdifferential of the objective is $Ax - g + λ + ∂ δ_{≥ 0}(x)$.
-        let d = v_i + λ;
-        if x_i > 0.0 || d < 0.0 {
-            val = val.max(d.abs());
+        let (mut lb, mut ub) = (v_i, v_i);
+        match x_i.partial_cmp(&0.0) {
+            Some(Greater) => { lb += λ; ub += λ },
+            // Less should not happen
+            Some(Less|Equal) => { lb = F::NEG_INFINITY; ub += λ },
+            None => {},
         }
+        val = max_interval_dist_to_zero(val, lb, ub);
     }
-    F::from_nalgebra_mixed(val)
+    val
 }
 
 /// Forward-backward splitting implementation of [`quadratic_nonneg`].
@@ -98,7 +103,7 @@
         iters +=1;
 
         backup.map(|_| {
-            min_subdifferential(&mut v, mA, x, g, λ)
+            F::from_nalgebra_mixed(min_subdifferential(&mut v, mA, x, g, λ))
         })
     });
 
@@ -281,7 +286,7 @@
         // 4. Report solution quality
         state.if_verbose(|| {
             // Calculate subdifferential at the FB step `x` that hasn't yet had `s` yet added.
-            min_subdifferential(&mut v, mA, x, g, λ)
+            F::from_nalgebra_mixed(min_subdifferential(&mut v, mA, x, g, λ))
         })
     });
 
@@ -291,7 +296,7 @@
 /// This function applies an iterative method for the solution of the quadratic non-negativity
 /// constrained problem
 /// <div>$$
-///     \min_{x ∈ ℝ^n} \frac{1}{2} x^⊤Ax - g^⊤ x + λ{\vec 1}^⊤ x + c + δ_{≥ 0}(x).
+///     \min_{x ∈ ℝ^n} \frac{1}{2} x^⊤Ax - g^⊤ x + λ{\vec 1}^⊤ x + δ_{≥ 0}(x).
 /// $$</div>
 /// Semismooth Newton or forward-backward are supported based on the setting in `method`.
 /// The parameter `mA` is matrix $A$, and `g` and `λ` are as in the mathematical formulation.
@@ -307,27 +312,28 @@
 ///
 /// This function returns the number of iterations taken.
 pub fn quadratic_nonneg<F, I>(
-    method : InnerMethod,
     mA : &DMatrix<F::MixedType>,
     g : &DVector<F::MixedType>,
-    //c_ : F,
     λ : F,
     x : &mut DVector<F::MixedType>,
-    τ : F,
+    mA_normest : F,
+    inner : &InnerSettings<F>,
     iterator : I
 ) -> usize
 where F : Float + ToNalgebraRealField,
       I : AlgIteratorFactory<F>
 {
-    
-    match method {
-        InnerMethod::FB =>
-            quadratic_nonneg_fb(mA, g, λ, x, τ, iterator),
+    let inner_τ = inner.fb_τ0 / mA_normest;
+
+    match inner.method {
+        InnerMethod::FB | InnerMethod::Default =>
+            quadratic_nonneg_fb(mA, g, λ, x, inner_τ, iterator),
         InnerMethod::SSN =>
-            quadratic_nonneg_ssn(mA, g, λ, x, τ, iterator).unwrap_or_else(|e| {
+            quadratic_nonneg_ssn(mA, g, λ, x, inner_τ, iterator).unwrap_or_else(|e| {
                 println!("{}", format!("{e}. Using FB fallback.").red());
                 let ins = InnerSettings::<F>::default();
-                quadratic_nonneg_fb(mA, g, λ, x, τ, ins.iterator_options)
-            })
+                quadratic_nonneg_fb(mA, g, λ, x, inner_τ, ins.iterator_options)
+            }),
+        InnerMethod::PDPS => unimplemented!(),
     }
 }