diff -r 88d994f2d63d -r d9318317ccac src/forward_pdps.rs
--- a/src/forward_pdps.rs	Tue Aug 04 23:47:29 2026 -0500
+++ b/src/forward_pdps.rs	Wed Aug 05 17:50:38 2026 -0500
@@ -141,12 +141,8 @@
 {
     // Check parameters
     ensure!(
-        config.τ0 > 0.0
-            && config.τ0 < 1.0
-            && config.σp0 > 0.0
-            && config.σp0 < 1.0
-            && config.σd0 >= 0.0
-            && config.σp0 * config.σd0 <= 1.0,
+        config.τ0 > 0.0 && config.σp0 > 0.0,
+        //&& config.τ0 < 1.0 && config.σp0 < 1.0
         "Invalid step length parameters"
     );
 
@@ -160,23 +156,23 @@
     let idOpZ = IdOp::new();
     let opKz_adj = opKz.adjoint();
     let (l, l_z) = Pair(prox_penalty, &idOpZ).step_length_bound_pair(&f)?;
-    // We need to satisfy
-    //
-    //     τσ_dM(1-σ_p L_z)/(1 - τ L) + [σ_p L_z + σ_pσ_d‖K_z‖^2] < 1
-    //                                  ^^^^^^^^^^^^^^^^^^^^^^^^^
-    // with 1 > σ_p L_z and 1 > τ L.
-    //
-    // To do so, we first solve σ_p and σ_d from standard PDPS step length condition
-    // ^^^^^ < 1. then we solve τ from  the rest.
-    // If opKZ is the zero operator, then we set σ_d = 0 for τ to be calculated correctly below.
-    let σ_d = if is_fb { 0.0 } else { config.σd0 / nKz };
-    let σ_p = config.σp0 / (l_z + config.σd0 * nKz);
-    // Observe that = 1 - ^^^^^^^^^^^^^^^^^^^^^ = 1 - σ_{p,0}
-    // We get the condition τσ_d M (1-σ_p L_z) < (1-σ_{p,0})*(1-τ L)
-    // ⟺ τ [ σ_d M (1-σ_p L_z) + (1-σ_{p,0}) L ] < (1-σ_{p,0})
-    let φ = 1.0 - config.σp0;
-    let a = 1.0 - σ_p * l_z;
-    let τ = config.τ0 * φ / (σ_d * bigM * a + φ * l);
+
+    let (τ, σ_p, σ_d) = if is_fb {
+        (config.τ0 / l, config.σp0 / l_z, 0.0)
+    } else {
+        ensure!(
+            //config.τ0 < 1.0 &&
+            config.σp0 < 1.0 && config.σd0 > 0.0 && config.σp0 * config.σd0 <= 1.0,
+            "Invalid step length parameters"
+        );
+        let σ_d = config.σd0 / nKz;
+        let σ_p = config.σp0 / (l_z + config.σd0 * nKz);
+        let φ = 1.0 - config.σp0;
+        let a = 1.0 - σ_p * l_z;
+        let τ = config.τ0 * φ / (σ_d * bigM * a + φ * l);
+        (τ, σ_p, σ_d)
+    };
+
     // Acceleration is not currently supported
     // let γ = dataterm.factor_of_strong_convexity();
     let ω = 1.0;