Mercurial > repos > pointsource_algs / file revision

/*!
Simple parametric forward model.
 */

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, 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;

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<E, F>,
    F: Float,
    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.write_observable(b, prefix)
    }

    /// Returns a zero observable
    fn zero_observable(&self) -> Self::Observable {
        self.0.zero_observable()
    }
}

#[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 𝒟, &'a IdOp<Z>>
where
    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,
{
    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))
    }
}

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

/*
/// This `impl` is bit of an abuse as the codomain of `Apre` is a [`Pair`] of a measure predual,
/// to which this `impl` applies, and another space.
impl<F, Apre, Z> LipschitzValues
for ColOp<Apre, IdOp<Z>>
where
    F : Float,
    Z : Clone + Space + ClosedAdd,
    Apre : LipschitzValues<FloatType = F>,
{
    type FloatType = F;
    /// Return (if one exists) a factor $L$ such that $A_*z$ is $L$-Lipschitz for all
    /// $z$ in the unit ball.
    fn value_unit_lipschitz_factor(&self) -> Option<Self::FloatType> {
        self.0.value_unit_lipschitz_factor()
    }

    /// Return (if one exists) a factor $L$ such that $∇A_*z$ is $L$-Lipschitz for all
    /// $z$ in the unit ball.
    fn value_diff_unit_lipschitz_factor(&self) -> Option<Self::FloatType> {
        self.0.value_diff_unit_lipschitz_factor()
    }
}
*/

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 + ClosedEuclidean<F>,
    A: Mapping<RNDM<N, F>, Codomain = Z>,
    A: BasicCurvatureBoundEstimates<F>,
{
    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)
            }
        }
    }
}