pointsource_algs: comparison src/forward

-:9738b51d90d7
+:4f468d35fa29
 /*!
 Simple parametric forward model.
 */
-use super::{AdjointProductBoundedBy, AdjointProductPairBoundedBy, BoundedCurvature, ForwardModel};
+use super::{BasicCurvatureBoundEstimates, BoundedCurvature, BoundedCurvatureGuess, ForwardModel};
-use crate::measures::RNDM;
+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::error::{DynError, DynResult};
-use alg_tools::linops::{IdOp, Linear, RowOp, ZeroOp, AXPY};
+use alg_tools::euclidean::ClosedEuclidean;
-use alg_tools::mapping::Space;
+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;
 impl<Domain, F, A, E> ForwardModel<Pair<Domain, A::Observable>, F, PairNorm<E, L2, L2>>
 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;
 fn write_observable(&self, b: &Self::Observable, prefix: String) -> DynError {
 self.0.zero_observable()
 }
 }
 #[replace_float_literals(F::cast_from(literal))]
-impl<Domain, F, A, D, Z> AdjointProductPairBoundedBy<Pair<Domain, Z>, D, IdOp<Z>>
+impl<'a, F, A, 𝒟, Z, const N: usize>
-for RowOp<A, IdOp<Z>>
+StepLengthBoundPair<F, QuadraticDataTerm<F, Pair<RNDM<N, F>, Z>, RowOp<A, IdOp<Z>>>>
+for Pair<&'a 𝒟, &'a IdOp<Z>>
 where
-Domain: Space,
+RNDM<N, F>: Space + for<'b> Norm<&'b 𝒟, F>,
-F: Float,
+F: Float + ToNalgebraRealField,
-Z: Clone + Space + ClosedAdd,
+𝒟: DiscreteMeasureOp<Loc<N, F>, F>,
-A: AdjointProductBoundedBy<Domain, D, FloatType = F, Codomain = Z>,
+Z: Clone + ClosedEuclidean<F>,
-A::Codomain: ClosedAdd,
+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)> {
+#[replace_float_literals(F::cast_from(literal))]
-self.0.adjoint_product_bound(d).map(|l_0| {
+impl<'a, F, A, Z, const N: usize>
-// [A_*; B_*][A, B] = [A_*A, A_* B; B_* A, B_* B] ≤ diag(2A_*A, 2B_*B)
+StepLengthBoundPair<F, QuadraticDataTerm<F, Pair<RNDM<N, F>, Z>, RowOp<A, IdOp<Z>>>>
-// ≤ diag(2l_A𝒟_A, 2l_B𝒟_B), where now 𝒟_B=Id and l_B=1.
+for Pair<&'a RadonSquared, &'a IdOp<Z>>
-(2.0 * l_0, 2.0)
+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))
 }
 }
 /*
 /// This `impl` is bit of an abuse as the codomain of `Apre` is a [`Pair`] of a measure predual,
 self.0.value_diff_unit_lipschitz_factor()
 }
 }
 */
-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,
+Z: Clone + ClosedEuclidean<F>,
-A: BoundedCurvature<FloatType = F>,
+A: Mapping<RNDM<N, F>, Codomain = Z>,
+A: BasicCurvatureBoundEstimates<F>,
 {
-type FloatType = F;
+fn curvature_bound_components(
+&self,
-fn curvature_bound_components(&self) -> (Option<Self::FloatType>, Option<Self::FloatType>) {
+guess: BoundedCurvatureGuess,
-self.0.curvature_bound_components()
+) -> (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)
-}
-}

comparison: src/forward_model/bias.rs

src/forward_model/bias.rs

Mercurial > repos > pointsource_algs / file comparison

comparison: src/forward_model/bias.rs

src/forward_model/bias.rs