$$ @@ -37,10 +36,6 @@ For $F_0(y)=\frac{1}{2}\|y\|_2^2$ the second part reads $y = Aμ -b$. For $F_0(y)=\|y\|_1$ the second part reads $y ∈ ∂\|·\|_1(Aμ - b)$.

- -Based on zero initialisation for $μ$, we use the [`Subdifferentiable`] trait to make an -initialisation corresponding to the second part of the optimality conditions. -In the algorithm itself, standard proximal steps are taking with respect to $F\_0^* + ⟨b, ·⟩$. */ use numeric_literals::replace_float_literals; @@ -48,13 +43,10 @@ use nalgebra::DVector; use clap::ValueEnum; -use alg_tools::iterate::{ - AlgIteratorFactory, - AlgIteratorState, -}; +use alg_tools::iterate::AlgIteratorFactory; use alg_tools::loc::Loc; use alg_tools::euclidean::Euclidean; -use alg_tools::linops::Apply; +use alg_tools::linops::Mapping; use alg_tools::norms::{ Linfinity, Projection, @@ -69,14 +61,17 @@ SupportGenerator, LocalAnalysis, }; -use alg_tools::mapping::RealMapping; +use alg_tools::mapping::{RealMapping, Instance}; use alg_tools::nalgebra_support::ToNalgebraRealField; use alg_tools::linops::AXPY; use crate::types::*; -use crate::measures::DiscreteMeasure; +use crate::measures::{DiscreteMeasure, RNDM, Radon}; use crate::measures::merging::SpikeMerging; -use crate::forward_model::ForwardModel; +use crate::forward_model::{ + AdjointProductBoundedBy, + ForwardModel +}; use crate::seminorms::DiscreteMeasureOp; use crate::plot::{ SeqPlotter, @@ -87,7 +82,7 @@ FBGenericConfig, insert_and_reweigh, postprocess, - prune_and_maybe_simple_merge + prune_with_stats }; use crate::regularisation::RegTerm; use crate::dataterm::{ @@ -110,7 +105,30 @@ Full } -/// Settings for [`pointsource_pdps`]. +#[replace_float_literals(F::cast_from(literal))] +impl Acceleration { + /// PDPS parameter acceleration. Updates τ and σ and returns ω. + /// This uses dual strong convexity, not primal. + fn accelerate(self, τ : &mut F, σ : &mut F, γ : F) -> F { + match self { + Acceleration::None => 1.0, + Acceleration::Partial => { + let ω = 1.0 / (1.0 + γ * (*σ)).sqrt(); + *σ *= ω; + *τ /= ω; + ω + }, + Acceleration::Full => { + let ω = 1.0 / (1.0 + 2.0 * γ * (*σ)).sqrt(); + *σ *= ω; + *τ /= ω; + ω + }, + } + } +} + +/// Settings for [`pointsource_pdps_reg`]. #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)] #[serde(default)] pub struct PDPSConfig { @@ -155,9 +173,13 @@ #[replace_float_literals(F::cast_from(literal))] -impl + AXPY, const N : usize> -PDPSDataTerm -for L2Squared { +impl PDPSDataTerm +for L2Squared +where + F : Float, + V : Euclidean + AXPY, + for<'b> &'b V : Instance, +{ fn some_subdifferential(&self, x : V) -> V { x } fn factor_of_strong_convexity(&self) -> F { @@ -166,7 +188,7 @@ #[inline] fn dual_update(&self, y : &mut V, y_prev : &V, σ : F) { - y.axpy(1.0 / (1.0 + σ), &y_prev, σ / (1.0 + σ)); + y.axpy(1.0 / (1.0 + σ), y_prev, σ / (1.0 + σ)); } } @@ -210,16 +232,13 @@ iterator : I, mut plotter : SeqPlotter, dataterm : D, -) -> DiscreteMeasure, F> +) -> RNDM where F : Float + ToNalgebraRealField, I : AlgIteratorFactory>, - for<'b> &'b A::Observable : std::ops::Neg - + std::ops::Add, - //+ std::ops::Mul, // <-- FIXME: compiler overflow - A::Observable : std::ops::MulAssign, + for<'b> &'b A::Observable : std::ops::Neg + Instance, GA : SupportGenerator + Clone, - A : ForwardModel, F, PreadjointCodomain = BTFN> - + Lipschitz<&'a 𝒟, FloatType=F>, + A : ForwardModel, F, PreadjointCodomain = BTFN> + + AdjointProductBoundedBy, 𝒟, FloatType=F>, BTA : BTSearch>, G𝒟 : SupportGenerator + Clone, 𝒟 : DiscreteMeasureOp, F, PreCodomain = PreBTFN>, @@ -228,14 +247,20 @@ K: RealMapping + LocalAnalysis, N>, BTNodeLookup: BTNode, N>, PlotLookup : Plotting, - DiscreteMeasure, F> : SpikeMerging, + RNDM : SpikeMerging, D : PDPSDataTerm, Reg : RegTerm { + // Check parameters + assert!(pdpsconfig.τ0 > 0.0 && + pdpsconfig.σ0 > 0.0 && + pdpsconfig.τ0 * pdpsconfig.σ0 <= 1.0, + "Invalid step length parameters"); + // Set up parameters let config = &pdpsconfig.generic; - let op𝒟norm = op𝒟.opnorm_bound(); - let l = opA.lipschitz_factor(&op𝒟).unwrap().sqrt(); + let op𝒟norm = op𝒟.opnorm_bound(Radon, Linfinity); + let l = opA.adjoint_product_bound(&op𝒟).unwrap().sqrt(); let mut τ = pdpsconfig.τ0 / l; let mut σ = pdpsconfig.σ0 / l; let γ = dataterm.factor_of_strong_convexity(); @@ -249,53 +274,42 @@ let mut μ = DiscreteMeasure::new(); let mut y = dataterm.some_subdifferential(-b); let mut y_prev = y.clone(); + let full_stats = |μ : &RNDM, ε, stats| IterInfo { + value : dataterm.calculate_fit_op(μ, opA, b) + reg.apply(μ), + n_spikes : μ.len(), + ε, + // postprocessing: config.postprocessing.then(|| μ.clone()), + .. stats + }; let mut stats = IterInfo::new(); // Run the algorithm - iterator.iterate(|state| { + for state in iterator.iter_init(|| full_stats(&μ, ε, stats.clone())) { // Calculate smooth part of surrogate model. - // Using `std::mem::replace` here is not ideal, and expects that `empty_observable` - // has no significant overhead. For some reosn Rust doesn't allow us simply moving - // the residual and replacing it below before the end of this closure. - y *= -τ; - let r = std::mem::replace(&mut y, opA.empty_observable()); - let minus_τv = opA.preadjoint().apply(r); + let τv = opA.preadjoint().apply(y * τ); // Save current base point let μ_base = μ.clone(); // Insert and reweigh - let (d, within_tolerances) = insert_and_reweigh( - &mut μ, &minus_τv, &μ_base, None, + let (d, _within_tolerances) = insert_and_reweigh( + &mut μ, &τv, &μ_base, None, op𝒟, op𝒟norm, τ, ε, - config, ®, state, &mut stats + config, ®, &state, &mut stats ); // Prune and possibly merge spikes - prune_and_maybe_simple_merge( - &mut μ, &minus_τv, &μ_base, - op𝒟, - τ, ε, - config, ®, state, &mut stats - ); + if config.merge_now(&state) { + stats.merged += μ.merge_spikes(config.merging, |μ_candidate| { + let mut d = &τv + op𝒟.preapply(μ_candidate.sub_matching(&μ_base)); + reg.verify_merge_candidate(&mut d, μ_candidate, τ, ε, &config) + }); + } + stats.pruned += prune_with_stats(&mut μ); // Update step length parameters - let ω = match pdpsconfig.acceleration { - Acceleration::None => 1.0, - Acceleration::Partial => { - let ω = 1.0 / (1.0 + γ * σ).sqrt(); - σ = σ * ω; - τ = τ / ω; - ω - }, - Acceleration::Full => { - let ω = 1.0 / (1.0 + 2.0 * γ * σ).sqrt(); - σ = σ * ω; - τ = τ / ω; - ω - }, - }; + let ω = pdpsconfig.acceleration.accelerate(&mut τ, &mut σ, γ); // Do dual update y = b.clone(); // y = b @@ -304,32 +318,17 @@ dataterm.dual_update(&mut y, &y_prev, σ); y_prev.copy_from(&y); - // Update main tolerance for next iteration - let ε_prev = ε; - ε = tolerance.update(ε, state.iteration()); + // Give statistics if requested + let iter = state.iteration(); stats.this_iters += 1; - // Give function value if needed state.if_verbose(|| { - // Plot if so requested - plotter.plot_spikes( - format!("iter {} end; {}", state.iteration(), within_tolerances), &d, - "start".to_string(), Some(&minus_τv), - reg.target_bounds(τ, ε_prev), &μ, - ); - // Calculate mean inner iterations and reset relevant counters. - // Return the statistics - let res = IterInfo { - value : dataterm.calculate_fit_op(&μ, opA, b) + reg.apply(&μ), - n_spikes : μ.len(), - ε : ε_prev, - postprocessing: config.postprocessing.then(|| μ.clone()), - .. stats - }; - stats = IterInfo::new(); - res - }) - }); + plotter.plot_spikes(iter, Some(&d), Some(&τv), &μ); + full_stats(&μ, ε, std::mem::replace(&mut stats, IterInfo::new())) + }); + + ε = tolerance.update(ε, iter); + } postprocess(μ, config, dataterm, opA, b) }