pointsource_algs: comparison src/pdps.rs

-:efa60bc4f743
+:b087e3eab191
 This corresponds to the manuscript
 * Valkonen T. - _Proximal methods for point source localisation_,
 [arXiv:2212.02991](https://arxiv.org/abs/2212.02991).
-The main routine is [`pointsource_pdps`]. It is based on specilisatinn of
+The main routine is [`pointsource_pdps_reg`].
-[`generic_pointsource_fb_reg`] through relevant [`FBSpecialisation`] implementations.
 Both norm-2-squared and norm-1 data terms are supported. That is, implemented are solvers for
 <div>
 $$
 \min_{μ ∈ ℳ(Ω)}~ F_0(Aμ - b) + α \|μ\|_{ℳ(Ω)} + δ_{≥ 0}(μ),
 $$
 This is the task of <code>generic_pointsource_fb</code>, where we use <code>FBSpecialisation</code>
 to replace the specific residual $Aμ-b$ by $y$.
 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)$.
 </p>
-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;
 use serde::{Serialize, Deserialize};
 use nalgebra::DVector;
 use clap::ValueEnum;
-use alg_tools::iterate::{
+use alg_tools::iterate::AlgIteratorFactory;
-AlgIteratorFactory,
-AlgIteratorState,
-};
 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,
 };
 use alg_tools::bisection_tree::{
 BTNode,
 BTSearch,
 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,
 Plotting,
 PlotLookup
 };
 use crate::fb::{
 FBGenericConfig,
 insert_and_reweigh,
 postprocess,
-prune_and_maybe_simple_merge
+prune_with_stats
 };
 use crate::regularisation::RegTerm;
 use crate::dataterm::{
 DataTerm,
 L2Squared,
 /// Full acceleration, $ω = 1/\sqrt{1+2σ}$; no gap convergence guaranteed
 #[clap(name = "full", help = "Full acceleration, ω = 1/√(1+2σ); no gap convergence guaranteed")]
 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<F : Float>(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<F : Float> {
 /// Primal step length scaling. We must have `τ0 * σ0 < 1`.
 pub τ0 : F,
 fn dual_update(&self, _y : &mut V, _y_prev : &V, _σ : F);
 }
 #[replace_float_literals(F::cast_from(literal))]
-impl<F : Float, V :  Euclidean<F> + AXPY<F>, const N : usize>
+impl<F, V, const N : usize> PDPSDataTerm<F, V, N>
-PDPSDataTerm<F, V, N>
+for L2Squared
-for L2Squared {
+where
+F : Float,
+V :  Euclidean<F> + AXPY<F>,
+for<'b> &'b V : Instance<V>,
+{
 fn some_subdifferential(&self, x : V) -> V { x }
 fn factor_of_strong_convexity(&self) -> F {
 1.0
 }
 #[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 + σ));
 }
 }
 #[replace_float_literals(F::cast_from(literal))]
 impl<F : Float + nalgebra::RealField, const N : usize>
 op𝒟 : &'a 𝒟,
 pdpsconfig : &PDPSConfig<F>,
 iterator : I,
 mut plotter : SeqPlotter<F, N>,
 dataterm : D,
-) -> DiscreteMeasure<Loc<F, N>, F>
+) -> RNDM<F, N>
 where F : Float + ToNalgebraRealField,
 I : AlgIteratorFactory<IterInfo<F, N>>,
-for<'b> &'b A::Observable : std::ops::Neg<Output=A::Observable>
+for<'b> &'b A::Observable : std::ops::Neg<Output=A::Observable> + Instance<A::Observable>,
-+ std::ops::Add<A::Observable, Output=A::Observable>,
-//+ std::ops::Mul<F, Output=A::Observable>, // <-- FIXME: compiler overflow
-A::Observable : std::ops::MulAssign<F>,
 GA : SupportGenerator<F, N, SupportType = S, Id = usize> + Clone,
-A : ForwardModel<Loc<F, N>, F, PreadjointCodomain = BTFN<F, GA, BTA, N>>
+A : ForwardModel<RNDM<F, N>, F, PreadjointCodomain = BTFN<F, GA, BTA, N>>
-+ Lipschitz<&'a 𝒟, FloatType=F>,
++ AdjointProductBoundedBy<RNDM<F, N>, 𝒟, FloatType=F>,
 BTA : BTSearch<F, N, Data=usize, Agg=Bounds<F>>,
 G𝒟 : SupportGenerator<F, N, SupportType = K, Id = usize> + Clone,
 𝒟 : DiscreteMeasureOp<Loc<F, N>, F, PreCodomain = PreBTFN<F, G𝒟, N>>,
 𝒟::Codomain : RealMapping<F, N>,
 S: RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>,
 K: RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>,
 BTNodeLookup: BTNode<F, usize, Bounds<F>, N>,
 PlotLookup : Plotting<N>,
-DiscreteMeasure<Loc<F, N>, F> : SpikeMerging<F>,
+RNDM<F, N> : SpikeMerging<F>,
 D : PDPSDataTerm<F, A::Observable, N>,
 Reg : RegTerm<F, N> {
+// 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 op𝒟norm = op𝒟.opnorm_bound(Radon, Linfinity);
-let l = opA.lipschitz_factor(&op𝒟).unwrap().sqrt();
+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();
 // We multiply tolerance by τ for FB since our subproblems depending on tolerances are scaled
 // Initialise iterates
 let mut μ = DiscreteMeasure::new();
 let mut y = dataterm.some_subdifferential(-b);
 let mut y_prev = y.clone();
+let full_stats = |μ : &RNDM<F, N>, ε, 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`
+let τv = opA.preadjoint().apply(y * τ);
-// 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);
 // Save current base point
 let μ_base = μ.clone();
 // Insert and reweigh
-let (d, within_tolerances) = insert_and_reweigh(
+let (d, _within_tolerances) = insert_and_reweigh(
-&mut μ, &minus_τv, &μ_base, None,
+&mut μ, &τv, &μ_base, None,
 op𝒟, op𝒟norm,
 τ, ε,
-config, &reg, state, &mut stats
+config, &reg, &state, &mut stats
 );
 // Prune and possibly merge spikes
-prune_and_maybe_simple_merge(
+if config.merge_now(&state) {
-&mut μ, &minus_τv, &μ_base,
+stats.merged += μ.merge_spikes(config.merging, |μ_candidate| {
-op𝒟,
+let mut d = &τv + op𝒟.preapply(μ_candidate.sub_matching(&μ_base));
-τ, ε,
+reg.verify_merge_candidate(&mut d, μ_candidate, τ, ε, &config)
-config, &reg, state, &mut stats
+});
-);
+}
+stats.pruned += prune_with_stats(&mut μ);
 // Update step length parameters
-let ω = match pdpsconfig.acceleration {
+let ω = pdpsconfig.acceleration.accelerate(&mut τ, &mut σ, γ);
-Acceleration::None => 1.0,
-Acceleration::Partial => {
-let ω = 1.0 / (1.0 + γ * σ).sqrt();
-σ = σ * ω;
-τ = τ / ω;
-ω
-},
-Acceleration::Full => {
-let ω = 1.0 / (1.0 + 2.0 * γ * σ).sqrt();
-σ = σ * ω;
-τ = τ / ω;
-ω
-},
-};
 // Do dual update
 y = b.clone();                          // y = b
 opA.gemv(&mut y, 1.0 + ω, &μ, -1.0);    // y = A[(1+ω)μ^{k+1}]-b
 opA.gemv(&mut y, -ω, &μ_base, 1.0);     // y = A[(1+ω)μ^{k+1} - ω μ^k]-b
 dataterm.dual_update(&mut y, &y_prev, σ);
 y_prev.copy_from(&y);
-// Update main tolerance for next iteration
+// Give statistics if requested
-let ε_prev = ε;
+let iter = state.iteration();
-ε = tolerance.update(ε, state.iteration());
 stats.this_iters += 1;
-// Give function value if needed
 state.if_verbose(|| {
-// Plot if so requested
+plotter.plot_spikes(iter, Some(&d), Some(&τv), &μ);
-plotter.plot_spikes(
+full_stats(&μ, ε, std::mem::replace(&mut stats, IterInfo::new()))
-format!("iter {} end; {}", state.iteration(), within_tolerances), &d,
+});
-"start".to_string(), Some(&minus_τv),
-reg.target_bounds(τ, ε_prev), &μ,
+ε = tolerance.update(ε, iter);
-);
+}
-// 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
-})
-});
 postprocess(μ, config, dataterm, opA, b)
 }

Mercurial > repos > pointsource_algs / file comparison

comparison: src/pdps.rs

src/pdps.rs