Mercurial > repos > pointsource_algs / changeset

--- a/src/subproblem/l1squared_nonneg.rs	Thu Feb 26 11:36:22 2026 -0500
+++ b/src/subproblem/l1squared_nonneg.rs	Mon Feb 23 18:18:02 2026 -0500
@@ -100,10 +100,43 @@
 /// which is also used internally.
 /// The third output indicates whether the output is “locked” to either $0$ (resp. $-y$ in
 /// the reformulation), or $y$ ($0$ in the reformulation).
-/// The fourth output indicates the next highest $λ$ where such locking would happen.
+/// The fourth output indicates $λ$ where such locking to y would happen. It may be negative,
+/// in which case the prox is always locked to y.
 #[replace_float_literals(F::cast_from(literal))]
 #[inline]
-fn shifted_nonneg_soft_thresholding<F: Float>(x: F, y: F, λ: F) -> (F, F, bool, Option<F>) {
+fn shifted_nonneg_soft_thresholding<F: Float>(x: F, y: F, λ: F) -> (F, F, bool, F) {
+    // The OC is
+    // x ∈ w + λsign(w-y) + N_{≥ 0}(w)
+    // x-y ∈ w' + λsign(w') + N_{≥ -y}(w')
+    // if w' > -y, then soft-thresholding. So if soft-thresholding stays above -y.
+    // If w' = -y, and w'≠0, then x ≤ λsign(-y)
+    // If w' = -y and w'=0, then x ≤ λ.
+    // w_i=max{sign(x-y)*max(|x-y| - λ, 0), -y} + y = max{sign(x-y)*max(|x-y| - λ, 0)+y, 0}
+    let (w, locked, lock) = if x >= y {
+        // w_i=max{max(x-y - λ, 0)+y, 0}=max{max(x - λ, y), 0} = max{x-λ, max(y, 0)}
+        let l = 0.0.max(y);
+        // x - l = x + min(0.0, y) = min(x, x + y)
+        if x - λ >= l {
+            (x - λ, false, x - l)
+        } else {
+            (l, true, F::NEG_INFINITY) // λ is increasing, so NEG_INFINITY is ok. So is x-l < λ
+        }
+    } else {
+        // w_i=max{-max(-(x-y) - λ, 0)+y, 0}=max{min((x-y) + λ, 0)+y, 0}=max{min(x + λ, y), 0}
+        if x + λ <= y {
+            // = max(x+λ, 0)
+            if x >= -λ {
+                (x + λ, false, y - x)
+            } else {
+                (0.0, true, 0.0.min(y) - x)
+            }
+        } else {
+            (0.0.max(y), true, F::NEG_INFINITY)
+        }
+    };
+    return (w, w - y, locked, lock);
+
+    /*
     let z = x - y; // The shift to f_0
     if -y >= 0.0 {
         if x > λ {
@@ -123,7 +156,7 @@
         (x - λ, z - λ, false, Some(z))
     } else {
         (y, 0.0, true, None)
-    }
+    }*/
 }
 
 /// Calculate $\prox\_f(x)$ for $f(x)=\frac{β}{2}\norm{x-y}\_1\^2 + δ\_{≥0}(x)$.
@@ -180,14 +213,15 @@
         let mut x_prime = 0.0;
         let mut max_shift = F::INFINITY;
         for (&x_i, &y_i) in izip!(x.iter(), y.iter()) {
-            let (_, a_shifted, locked, next_lock) = shifted_nonneg_soft_thresholding(x_i, y_i, λ);
+            let (_, w_i, locked, lock) = shifted_nonneg_soft_thresholding(x_i, y_i, λ);
 
-            if let Some(t) = next_lock {
-                max_shift = min(max_shift, t);
-                assert!(max_shift > λ);
+            // If not already locked (lock <= y), consider this coordinate in next level for λ.
+            if lock > λ {
+                max_shift = min(max_shift, lock);
             }
+
             if locked {
-                w_locked += abs(a_shifted);
+                w_locked += abs(w_i);
             } else {
                 n_unlocked += 1;
                 x_prime += abs(x_i - y_i);
@@ -203,9 +237,7 @@
             assert!(λ_new >= λ);
             // success
             x.zip_apply(y, |x_i, y_i| {
-                let (a, _, _, _) = shifted_nonneg_soft_thresholding(*x_i, y_i, λ_new);
-                //*x_i = y_i + lbd_soft_thresholding(*x_i, λ_new, -y_i)
-                *x_i = a;
+                (*x_i, _, _, _) = shifted_nonneg_soft_thresholding(*x_i, y_i, λ_new);
             });
             return;
         }