diff -r e2953ffd4e0b -r e9a460a0e638 src/subproblem/l1squared_unconstrained.rs
--- a/src/subproblem/l1squared_unconstrained.rs	Fri May 15 14:40:02 2026 -0500
+++ b/src/subproblem/l1squared_unconstrained.rs	Sun Jul 19 07:34:39 2026 +0200
@@ -3,13 +3,12 @@
 */
 
 use itertools::izip;
-use nalgebra::DVector;
+use nalgebra::{constraint::ShapeConstraint, DVector, Dyn, Storage, StorageMut, Vector, U1};
 use numeric_literals::replace_float_literals;
 use std::cmp::Ordering::*;
 
 use alg_tools::iterate::{AlgIteratorFactory, AlgIteratorState};
-use alg_tools::nalgebra_support::ToNalgebraRealField;
-use alg_tools::nanleast::NaNLeast;
+use alg_tools::nalgebra_support::{StridesOk, ToNalgebraRealField};
 use alg_tools::norms::{Dist, L1};
 use std::iter::zip;
 
@@ -44,31 +43,47 @@
 /// Clearly, if this condition fails for $x\_i$, it will fail for all the components
 /// already exluced. While, if it holds, it will hold for all components not excluded.
 #[replace_float_literals(F::cast_from(literal))]
-pub(super) fn l1squared_prox<F: Float + nalgebra::RealField>(
+pub(super) fn l1squared_prox<F: Float + nalgebra::RealField, S1, S2>(
     sorted_abs: &mut DVector<F>,
-    x: &mut DVector<F>,
-    y: &DVector<F>,
+    x: &mut Vector<F, Dyn, S1>,
+    y: &Vector<F, Dyn, S2>,
     β: F,
-) {
+) where
+    S2: Storage<F, Dyn>,
+    S1: StorageMut<F, Dyn>,
+{
+    //let orig_x = x.clone();
     sorted_abs.copy_from(x);
     sorted_abs.axpy(-1.0, y, 1.0);
     sorted_abs.apply(|z_i| *z_i = num_traits::abs(*z_i));
-    sorted_abs
-        .as_mut_slice()
-        .sort_unstable_by(|a, b| NaNLeast(*a).cmp(&NaNLeast(*b)));
+    sorted_abs.as_mut_slice().sort_unstable_by(F::total_cmp);
 
     let mut n = sorted_abs.sum();
     for (m, az_i) in zip((1..=x.len() as u32).rev(), sorted_abs) {
-        // test first
-        let tmp = β * n / (1.0 + β * F::cast_from(m));
-        if *az_i <= tmp {
+        // test first. This is just *az_i <= tmp, for tmp defined below, without the division.
+        if *az_i <= β * (n - F::cast_from(m) * *az_i) {
             // Fail
             n -= *az_i;
         } else {
             // Success
+            let tmp = β * n / (1.0 + β * F::cast_from(m));
             x.zip_apply(y, |x_i, y_i| {
                 *x_i = y_i + soft_thresholding(*x_i - y_i, tmp)
             });
+            // //Check 0 ∈ w-x + β\norm{w-y}\_1\sign (w-y).
+            // let n: F = izip!(x.iter(), y)
+            //     .map(|(&w_i, &y_i)| NumTraitsFloat::abs(w_i - y_i))
+            //     .sum();
+            // for (&mut w_i, &x_i, &y_i) in izip!(x, &orig_x, y) {
+            //     if w_i > y_i {
+            //         assert_lt!(NumTraitsFloat::abs(w_i - x_i + n * β), 10.0 * F::EPSILON);
+            //     } else if w_i < y_i {
+            //         assert_lt!(NumTraitsFloat::abs(w_i - x_i - n * β), 10.0 * F::EPSILON);
+            //     } else {
+            //         assert_lt!(-n * β - 10.0 * F::EPSILON, w_i - x_i);
+            //         assert_lt!(w_i - x_i, n * β + 10.0 * F::EPSILON);
+            //     }
+            // }
             return;
         }
     }
@@ -80,15 +95,20 @@
 ///
 /// `v` will be modified and cannot be trusted to contain useful values afterwards.
 #[replace_float_literals(F::cast_from(literal))]
-fn min_subdifferential<F: Float + nalgebra::RealField>(
-    y: &DVector<F>,
-    x: &DVector<F>,
-    g: &DVector<F>,
+fn min_subdifferential<F: Float + nalgebra::RealField, S1, S2, S3>(
+    y: &Vector<F, Dyn, S1>,
+    x: &Vector<F, Dyn, S2>,
+    g: &Vector<F, Dyn, S3>,
     λ: F,
-    β: F,
-) -> F {
+) -> F
+where
+    S1: Storage<F, Dyn>,
+    S2: Storage<F, Dyn>,
+    S3: Storage<F, Dyn>,
+    ShapeConstraint: StridesOk<F, Dyn, U1, S2>,
+{
     let mut val = 0.0;
-    let tmp = β * y.dist(x, L1);
+    let tmp = y.dist(x, L1);
     for (&g_i, &x_i, y_i) in izip!(g.iter(), x.iter(), y.iter()) {
         let (mut lb, mut ub) = (-g_i, -g_i);
         match x_i.partial_cmp(y_i) {
@@ -132,12 +152,11 @@
 /// The `λ` component of the model is handled in the proximal step instead of the gradient step
 /// for potential performance improvements.
 #[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())]
-pub fn l1squared_unconstrained_pdps<F, I>(
-    y: &DVector<F::MixedType>,
-    g: &DVector<F::MixedType>,
+pub fn l1squared_unconstrained_pdps<F, I, S1, S2, S3>(
+    y: &Vector<F::MixedType, Dyn, S1>,
+    g: &Vector<F::MixedType, Dyn, S2>,
     λ_: F,
-    β_: F,
-    x: &mut DVector<F::MixedType>,
+    x: &mut Vector<F::MixedType, Dyn, S3>,
     τ_: F,
     σ_: F,
     iterator: I,
@@ -145,21 +164,24 @@
 where
     F: Float + ToNalgebraRealField,
     I: AlgIteratorFactory<F>,
+    S1: Storage<F::MixedType, Dyn>,
+    S2: Storage<F::MixedType, Dyn>,
+    S3: StorageMut<F::MixedType, Dyn>,
+    ShapeConstraint: StridesOk<F::MixedType, Dyn, U1, S3>,
 {
     let λ = λ_.to_nalgebra_mixed();
-    let β = β_.to_nalgebra_mixed();
     let τ = τ_.to_nalgebra_mixed();
     let σ = σ_.to_nalgebra_mixed();
     let mut w = DVector::zeros(x.len());
     let mut tmp = DVector::zeros(x.len());
-    let mut xprev = x.clone();
+    let mut xprev = x.clone_owned();
     let mut iters = 0;
 
     iterator.iterate(|state| {
         // Primal step: x^{k+1} = prox_{τ|.-y|_1^2}(x^k - τ (w^k - g))
         x.axpy(-τ, &w, 1.0);
         x.axpy(τ, g, 1.0);
-        l1squared_prox(&mut tmp, x, y, τ * β);
+        l1squared_prox(&mut tmp, x, y, τ);
 
         // Dual step: w^{k+1} = proj_{[-λ,λ]}(w^k + σ(2x^{k+1}-x^k))
         w.axpy(2.0 * σ, x, 1.0);
@@ -169,7 +191,7 @@
 
         iters += 1;
 
-        state.if_verbose(|| F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ, β)))
+        state.if_verbose(|| F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ)))
     });
 
     iters
@@ -193,12 +215,11 @@
 ///     \end{split}
 /// $$</div>
 #[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())]
-pub fn l1squared_unconstrained_pdps_alt<F, I>(
-    y: &DVector<F::MixedType>,
-    g: &DVector<F::MixedType>,
+pub fn l1squared_unconstrained_pdps_alt<F, I, S1, S2, S3>(
+    y: &Vector<F::MixedType, Dyn, S1>,
+    g: &Vector<F::MixedType, Dyn, S2>,
     λ_: F,
-    β_: F,
-    x: &mut DVector<F::MixedType>,
+    x: &mut Vector<F::MixedType, Dyn, S3>,
     τ_: F,
     σ_: F,
     θ_: F,
@@ -207,17 +228,20 @@
 where
     F: Float + ToNalgebraRealField,
     I: AlgIteratorFactory<F>,
+    S1: Storage<F::MixedType, Dyn>,
+    S2: Storage<F::MixedType, Dyn>,
+    S3: StorageMut<F::MixedType, Dyn>,
+    ShapeConstraint: StridesOk<F::MixedType, Dyn, U1, S3>,
 {
     let λ = λ_.to_nalgebra_mixed();
     let τ = τ_.to_nalgebra_mixed();
     let σ = σ_.to_nalgebra_mixed();
     let θ = θ_.to_nalgebra_mixed();
-    let β = β_.to_nalgebra_mixed();
     let σθ = σ * θ;
     let τθ = τ * θ;
     let mut w = DVector::zeros(x.len());
     let mut tmp = DVector::zeros(x.len());
-    let mut xprev = x.clone();
+    let mut xprev = x.clone_owned();
     let mut iters = 0;
 
     iterator.iterate(|state| {
@@ -236,14 +260,14 @@
         w.axpy(2.0, x, 1.0);
         w.axpy(-1.0, &xprev, 1.0);
         xprev.copy_from(&w); // use xprev as temporary variable
-        l1squared_prox(&mut tmp, &mut xprev, y, β / σθ);
+        l1squared_prox(&mut tmp, &mut xprev, y, 1.0 / σθ);
         w -= &xprev;
         w *= σ;
         xprev.copy_from(x);
 
         iters += 1;
 
-        state.if_verbose(|| F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ, β)))
+        state.if_verbose(|| F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ)))
     });
 
     iters
@@ -257,18 +281,21 @@
 ///
 /// This function returns the number of iterations taken.
 #[replace_float_literals(F::cast_from(literal))]
-pub fn l1squared_unconstrained<F, I>(
-    y: &DVector<F::MixedType>,
-    g: &DVector<F::MixedType>,
+pub fn l1squared_unconstrained<F, I, S1, S2, S3>(
+    y: &Vector<F::MixedType, Dyn, S1>,
+    g: &Vector<F::MixedType, Dyn, S2>,
     λ: F,
-    β: F,
-    x: &mut DVector<F::MixedType>,
+    x: &mut Vector<F::MixedType, Dyn, S3>,
     inner: &InnerSettings<F>,
     iterator: I,
 ) -> usize
 where
     F: Float + ToNalgebraRealField,
     I: AlgIteratorFactory<F>,
+    S1: Storage<F::MixedType, Dyn>,
+    S2: Storage<F::MixedType, Dyn>,
+    S3: StorageMut<F::MixedType, Dyn>,
+    ShapeConstraint: StridesOk<F::MixedType, Dyn, U1, S3>,
 {
     // Estimate of ‖K‖ for K=θ Id.
     let inner_θ = 1.0;
@@ -278,7 +305,7 @@
 
     match inner.method {
         InnerMethod::PDPS => {
-            l1squared_unconstrained_pdps_alt(y, g, λ, β, x, inner_τ, inner_σ, inner_θ, iterator)
+            l1squared_unconstrained_pdps_alt(y, g, λ, x, inner_τ, inner_σ, inner_θ, iterator)
         }
         other => unimplemented!("${other:?} is unimplemented"),
     }