diff -r 9738b51d90d7 -r 4f468d35fa29 src/forward_model/bias.rs
--- a/src/forward_model/bias.rs	Sun Apr 27 15:03:51 2025 -0500
+++ b/src/forward_model/bias.rs	Thu Feb 26 11:38:43 2026 -0500
@@ -2,12 +2,18 @@
 Simple parametric forward model.
  */
 
-use super::{AdjointProductBoundedBy, AdjointProductPairBoundedBy, BoundedCurvature, ForwardModel};
-use crate::measures::RNDM;
+use super::{BasicCurvatureBoundEstimates, BoundedCurvature, BoundedCurvatureGuess, ForwardModel};
+use crate::dataterm::QuadraticDataTerm;
+use crate::measures::{Radon, RNDM};
+use crate::prox_penalty::{RadonSquared, StepLengthBoundPair};
+use crate::seminorms::DiscreteMeasureOp;
 use alg_tools::direct_product::Pair;
-use alg_tools::error::DynError;
-use alg_tools::linops::{IdOp, Linear, RowOp, ZeroOp, AXPY};
-use alg_tools::mapping::Space;
+use alg_tools::error::{DynError, DynResult};
+use alg_tools::euclidean::ClosedEuclidean;
+use alg_tools::linops::{BoundedLinear, IdOp, RowOp, AXPY};
+use alg_tools::loc::Loc;
+use alg_tools::mapping::{Mapping, Space};
+use alg_tools::nalgebra_support::ToNalgebraRealField;
 use alg_tools::norms::{Norm, NormExponent, PairNorm, L2};
 use alg_tools::types::{ClosedAdd, Float};
 use numeric_literals::replace_float_literals;
@@ -16,9 +22,9 @@
     for RowOp<A, IdOp<A::Observable>>
 where
     E: NormExponent,
-    Domain: Space + Norm<F, E>,
+    Domain: Space + Norm<E, F>,
     F: Float,
-    A::Observable: ClosedAdd + Norm<F, L2> + 'static,
+    A::Observable: ClosedAdd + Norm<L2, F> + AXPY<Field = F> + 'static,
     A: ForwardModel<Domain, F, E> + 'static,
 {
     type Observable = A::Observable;
@@ -34,23 +40,46 @@
 }
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<Domain, F, A, D, Z> AdjointProductPairBoundedBy<Pair<Domain, Z>, D, IdOp<Z>>
-    for RowOp<A, IdOp<Z>>
+impl<'a, F, A, 𝒟, Z, const N: usize>
+    StepLengthBoundPair<F, QuadraticDataTerm<F, Pair<RNDM<N, F>, Z>, RowOp<A, IdOp<Z>>>>
+    for Pair<&'a 𝒟, &'a IdOp<Z>>
 where
-    Domain: Space,
-    F: Float,
-    Z: Clone + Space + ClosedAdd,
-    A: AdjointProductBoundedBy<Domain, D, FloatType = F, Codomain = Z>,
-    A::Codomain: ClosedAdd,
+    RNDM<N, F>: Space + for<'b> Norm<&'b 𝒟, F>,
+    F: Float + ToNalgebraRealField,
+    𝒟: DiscreteMeasureOp<Loc<N, F>, F>,
+    Z: Clone + ClosedEuclidean<F>,
+    A: for<'b> BoundedLinear<RNDM<N, F>, &'b 𝒟, L2, F, Codomain = Z>,
+    for<'b> &'b 𝒟: NormExponent,
 {
-    type FloatType = F;
+    fn step_length_bound_pair(
+        &self,
+        f: &QuadraticDataTerm<F, Pair<RNDM<N, F>, Z>, RowOp<A, IdOp<Z>>>,
+    ) -> DynResult<(F, F)> {
+        let l_0 = f.operator().0.opnorm_bound(self.0, L2)?.powi(2);
+        // [A_*; B_*][A, B] = [A_*A, A_* B; B_* A, B_* B] ≤ diag(2A_*A, 2B_*B)
+        // ≤ diag(2l_A𝒟_A, 2l_B𝒟_B), where now 𝒟_B=Id and l_B=1.
+        Ok((2.0 * l_0, 2.0))
+    }
+}
 
-    fn adjoint_product_pair_bound(&self, d: &D, _: &IdOp<Z>) -> Option<(F, F)> {
-        self.0.adjoint_product_bound(d).map(|l_0| {
-            // [A_*; B_*][A, B] = [A_*A, A_* B; B_* A, B_* B] ≤ diag(2A_*A, 2B_*B)
-            // ≤ diag(2l_A𝒟_A, 2l_B𝒟_B), where now 𝒟_B=Id and l_B=1.
-            (2.0 * l_0, 2.0)
-        })
+#[replace_float_literals(F::cast_from(literal))]
+impl<'a, F, A, Z, const N: usize>
+    StepLengthBoundPair<F, QuadraticDataTerm<F, Pair<RNDM<N, F>, Z>, RowOp<A, IdOp<Z>>>>
+    for Pair<&'a RadonSquared, &'a IdOp<Z>>
+where
+    RNDM<N, F>: Space + Norm<Radon, F>,
+    F: Float + ToNalgebraRealField,
+    Z: Clone + ClosedEuclidean<F>,
+    A: BoundedLinear<RNDM<N, F>, Radon, L2, F, Codomain = Z>,
+{
+    fn step_length_bound_pair(
+        &self,
+        f: &QuadraticDataTerm<F, Pair<RNDM<N, F>, Z>, RowOp<A, IdOp<Z>>>,
+    ) -> DynResult<(F, F)> {
+        let l_0 = f.operator().0.opnorm_bound(Radon, L2)?.powi(2);
+        // [A_*; B_*][A, B] = [A_*A, A_* B; B_* A, B_* B] ≤ diag(2A_*A, 2B_*B)
+        // ≤ diag(2l_A𝒟_A, 2l_B𝒟_B), where now 𝒟_B=Id and l_B=1.
+        Ok((2.0 * l_0, 2.0))
     }
 }
 
@@ -79,30 +108,44 @@
 }
 */
 
-impl<F, A, Z> BoundedCurvature for RowOp<A, IdOp<Z>>
+use BoundedCurvatureGuess::*;
+
+/// Curvature error control: helper bounds for (4.2d), (5.2a), (5.2b), (5.15a), and (5.16a).
+///
+/// Based on Lemma 5.11 and Example 5.12, the helper bound for (5.15a) and (5.16a) is (3.8).
+/// Due to Example 6.1, defining $v^k$ as the projection $F'$ to the predual space of the
+/// measures, returns, if possible, and subject to the guess being correct, factors $ℓ_F$ and
+/// $Θ²$ such that $B_{P_ℳ^* F'(μ, z)} dγ ≤ ℓ_F c_2$ and
+/// $⟨P_ℳ^*[F'(μ+Δ, z)-F'(μ, z)]|Δ⟩ ≤ Θ²|γ|(c_2)‖γ‖$, where $Δ=(π_♯^1-π_♯^0)γ$.
+/// For our $F(μ, z)=\frac{1}{2}\|Aμ+z-b\|^2$, we have $F'(μ, z)=A\_*(Aμ+z-b)$, so
+/// $F'(μ+Δ, z)-F'(μ, z)=A\_*AΔ$ is independent of $z$, and the bounding can be calculated
+/// as in the case without $z$, based on Lemma 3.8.
+///
+/// We use Remark 5.15 and Example 5.16 for (4.2d) and (5.2a) with the additional effect of $z$.
+/// This is based on a Lipschitz estimate for $∇v^k$, where we still, similarly to the Example,
+/// have $∇v^k(x)=∇A\_*(x)[Aμ^k+z^k-b]$. We estimate the final term similarly to the example,
+/// assuming for the guess [`BetterThanZero`] that every iterate is better than $(μ, z)=0$.
+/// This the final estimate is exactly as in the example, without $z$.
+/// Thus we can directly use [`BasicCurvatureBoundEstimates`] on the operator $A$.
+impl<F, A, Z, const N: usize> BoundedCurvature<F>
+    for QuadraticDataTerm<F, Pair<RNDM<N, F>, Z>, RowOp<A, IdOp<Z>>>
 where
     F: Float,
-    Z: Clone + Space + ClosedAdd,
-    A: BoundedCurvature<FloatType = F>,
+    Z: Clone + ClosedEuclidean<F>,
+    A: Mapping<RNDM<N, F>, Codomain = Z>,
+    A: BasicCurvatureBoundEstimates<F>,
 {
-    type FloatType = F;
-
-    fn curvature_bound_components(&self) -> (Option<Self::FloatType>, Option<Self::FloatType>) {
-        self.0.curvature_bound_components()
+    fn curvature_bound_components(
+        &self,
+        guess: BoundedCurvatureGuess,
+    ) -> (DynResult<F>, DynResult<F>) {
+        match guess {
+            BetterThanZero => {
+                let opA = &self.operator().0;
+                let b = self.data();
+                let (ℓ_F0, θ2) = opA.basic_curvature_bound_components();
+                (ℓ_F0.map(|l| l * b.norm2()), θ2)
+            }
+        }
     }
 }
-
-#[replace_float_literals(F::cast_from(literal))]
-impl<'a, F, D, XD, Y, const N: usize> AdjointProductBoundedBy<RNDM<F, N>, D>
-    for ZeroOp<'a, RNDM<F, N>, XD, Y, F>
-where
-    F: Float,
-    Y: AXPY<F> + Clone,
-    D: Linear<RNDM<F, N>>,
-{
-    type FloatType = F;
-    /// Return $L$ such that $A_*A ≤ L𝒟$ is bounded by some `other` operator $𝒟$.
-    fn adjoint_product_bound(&self, _: &D) -> Option<F> {
-        Some(0.0)
-    }
-}