diff -r 6105b5cd8d89 -r f0e8704d3f0e src/subproblem/unconstrained.rs
--- a/src/subproblem/unconstrained.rs	Tue Aug 01 10:25:09 2023 +0300
+++ b/src/subproblem/unconstrained.rs	Mon Feb 17 13:54:53 2025 -0500
@@ -23,11 +23,12 @@
     InnerMethod,
     InnerSettings
 };
+use super::l1squared_nonneg::max_interval_dist_to_zero;
 
 /// Compute the proximal operator of $x \mapsto |x|$, i.e., the soft-thresholding operator.
 #[inline]
 #[replace_float_literals(F::cast_from(literal))]
-fn soft_thresholding<F : Float>(v : F, λ : F) -> F {
+pub(crate) fn soft_thresholding<F : Float>(v : F, λ : F) -> F {
     if v > λ {
         v - λ
     } else if v < -λ {
@@ -40,27 +41,28 @@
 /// Returns the ∞-norm minimal subdifferential of $x ↦  x^⊤Ax - g^⊤ x + λ\|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 at x is $Ax - g + λ ∂‖·‖₁(x)$.
-        val = val.max(match x_i.partial_cmp(&0.0) {
-            Some(Greater) => v_i + λ,
-            Some(Less) => v_i - λ,
-            Some(Equal) => soft_thresholding(v_i, λ),
-            None => F::MixedType::nan(),
-        })
+        let (mut lb, mut ub) = (v_i, v_i);
+        match x_i.partial_cmp(&0.0) {
+            Some(Greater) => { lb += λ; ub += λ },
+            Some(Less) => { lb -= λ; ub -= λ },
+            Some(Equal) => { lb -= λ; ub += λ },
+            None => {},
+        }
+        val = max_interval_dist_to_zero(val, lb, ub);
     }
-    F::from_nalgebra_mixed(val)
+    val
 }
 
 
@@ -106,7 +108,7 @@
         iters +=1;
 
         backup.map(|_| {
-            min_subdifferential(&mut v, mA, x, g, λ)
+            F::from_nalgebra_mixed(min_subdifferential(&mut v, mA, x, g, λ))
         })
     });
 
@@ -242,7 +244,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, λ))
         })
     });
 
@@ -251,7 +253,7 @@
 
 /// This function applies an iterative method for the solution of the problem
 /// <div>$$
-///     \min_{x ∈ ℝ^n} \frac{1}{2} x^⊤Ax - g^⊤ x + λ\|x\|₁ + c.
+///     \min_{x ∈ ℝ^n} \frac{1}{2} x^⊤Ax - g^⊤ x + λ\|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.
@@ -262,27 +264,28 @@
 ///
 /// This function returns the number of iterations taken.
 pub fn quadratic_unconstrained<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>
 {
+    let inner_τ = inner.fb_τ0 / mA_normest;
     
-    match method {
+    match inner.method {
         InnerMethod::FB =>
-            quadratic_unconstrained_fb(mA, g, λ, x, τ, iterator),
+            quadratic_unconstrained_fb(mA, g, λ, x, inner_τ, iterator),
         InnerMethod::SSN =>
-            quadratic_unconstrained_ssn(mA, g, λ, x, τ, iterator).unwrap_or_else(|e| {
+            quadratic_unconstrained_ssn(mA, g, λ, x, inner_τ, iterator).unwrap_or_else(|e| {
                 println!("{}", format!("{e}. Using FB fallback.").red());
                 let ins = InnerSettings::<F>::default();
-                quadratic_unconstrained_fb(mA, g, λ, x, τ, ins.iterator_options)
-            })
+                quadratic_unconstrained_fb(mA, g, λ, x, inner_τ, ins.iterator_options)
+            }),
+        other => unimplemented!("${other:?} is unimplemented"),
     }
 }