--- 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)
- }
-}