diff -r e2953ffd4e0b -r e9a460a0e638 src/regularisation.rs
--- a/src/regularisation.rs	Fri May 15 14:40:02 2026 -0500
+++ b/src/regularisation.rs	Sun Jul 19 07:34:39 2026 +0200
@@ -234,8 +234,11 @@
     where
         M: MinMaxMapping<Domain, F>;
 
-    /// Convert bound on the regulariser to a bond on the Radon norm
+    /// Convert bound on the regulariser to a bound on the Radon norm
     fn radon_norm_bound(&self, b: F) -> F;
+
+    /// Returns true if $v$ is within the pointwise range of the subdifferential
+    fn subdiff_range(&self) -> Bounds<F>;
 }
 
 #[replace_float_literals(F::cast_from(literal))]
@@ -386,8 +389,7 @@
             // Solve finite-dimensional subproblem.
             let inner_tolerance = ε * config.inner.tolerance_mult;
             let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
-            stats.inner_iters +=
-                l1squared_nonneg(&y, &g_na, τα, 1.0, &mut x, &config.inner, inner_it);
+            stats.inner_iters += l1squared_nonneg(&y, &g_na, τα, &mut x, &config.inner, inner_it);
 
             // Update masses of μ based on solution of finite-dimensional subproblem.
             μ.set_masses_dvector(&x);
@@ -416,12 +418,12 @@
             μ.both_matching(radon_μ).all(|(α, rα, x)| {
                 let v = -d.apply(x); // TODO: observe ad hoc negation here, after minus_τv
                                      // switch to τv.
-                let (l1, u1) = match α.partial_cmp(&0.0).unwrap_or(Equal) {
+                let (l1, u1) = match α.total_cmp(&0.0) {
                     Greater => (τα, τα),
                     _ => (F::NEG_INFINITY, τα),
                     // Less should not happen; treated as Equal
                 };
-                let (l2, u2) = match rα.partial_cmp(&0.0).unwrap_or(Equal) {
+                let (l2, u2) = match rα.total_cmp(&0.0) {
                     Greater => (slack, slack),
                     Equal => (-slack, slack),
                     Less => (-slack, -slack),
@@ -465,6 +467,10 @@
     fn radon_norm_bound(&self, b: F) -> F {
         b / self.α()
     }
+
+    fn subdiff_range(&self) -> Bounds<F> {
+        Bounds(F::NEG_INFINITY, self.α())
+    }
 }
 
 #[replace_float_literals(F::cast_from(literal))]
@@ -644,7 +650,7 @@
             let inner_tolerance = ε * config.inner.tolerance_mult;
             let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
             stats.inner_iters +=
-                l1squared_unconstrained(&y, &g_na, τα, 1.0, &mut x, &config.inner, inner_it);
+                l1squared_unconstrained(&y, &g_na, τα, &mut x, &config.inner, inner_it);
 
             // Update masses of μ based on solution of finite-dimensional subproblem.
             μ.set_masses_dvector(&x);
@@ -732,4 +738,9 @@
     fn radon_norm_bound(&self, b: F) -> F {
         b / self.α()
     }
+
+    fn subdiff_range(&self) -> Bounds<F> {
+        let α = self.α();
+        Bounds(-α, α)
+    }
 }