pointsource_algs: comparison src/pdps.rs

-:9738b51d90d7
+:4f468d35fa29
 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>
 */
+use crate::fb::{postprocess, prune_with_stats};
+use crate::forward_model::ForwardModel;
+use crate::measures::merging::SpikeMerging;
+use crate::measures::merging::SpikeMergingMethod;
+use crate::measures::{DiscreteMeasure, RNDM};
+use crate::plot::Plotter;
+pub use crate::prox_penalty::{InsertionConfig, ProxPenalty, StepLengthBoundPD};
+use crate::regularisation::RegTerm;
+use crate::types::*;
+use alg_tools::convex::{Conjugable, ConvexMapping, Prox};
+use alg_tools::error::DynResult;
+use alg_tools::iterate::AlgIteratorFactory;
+use alg_tools::linops::{Mapping, AXPY};
+use alg_tools::mapping::{DataTerm, Instance};
+use alg_tools::nalgebra_support::ToNalgebraRealField;
+use anyhow::ensure;
+use clap::ValueEnum;
 use numeric_literals::replace_float_literals;
-use serde::{Serialize, Deserialize};
+use serde::{Deserialize, Serialize};
-use nalgebra::DVector;
-use clap::ValueEnum;
-use alg_tools::iterate::AlgIteratorFactory;
-use alg_tools::euclidean::Euclidean;
-use alg_tools::linops::Mapping;
-use alg_tools::norms::{
-Linfinity,
-Projection,
-};
-use alg_tools::mapping::{RealMapping, Instance};
-use alg_tools::nalgebra_support::ToNalgebraRealField;
-use alg_tools::linops::AXPY;
-use crate::types::*;
-use crate::measures::{DiscreteMeasure, RNDM};
-use crate::measures::merging::SpikeMerging;
-use crate::forward_model::{
-ForwardModel,
-AdjointProductBoundedBy,
-};
-use crate::plot::{
-SeqPlotter,
-Plotting,
-PlotLookup
-};
-use crate::fb::{
-postprocess,
-prune_with_stats
-};
-pub use crate::prox_penalty::{
-FBGenericConfig,
-ProxPenalty
-};
-use crate::regularisation::RegTerm;
-use crate::dataterm::{
-DataTerm,
-L2Squared,
-L1
-};
-use crate::measures::merging::SpikeMergingMethod;
 /// Acceleration
 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, ValueEnum, Debug)]
 pub enum Acceleration {
 /// No acceleration
 None,
 /// Partial acceleration, $ω = 1/\sqrt{1+σ}$
 #[clap(name = "partial", help = "Partial acceleration, ω = 1/√(1+σ)")]
 Partial,
 /// Full acceleration, $ω = 1/\sqrt{1+2σ}$; no gap convergence guaranteed
-#[clap(name = "full", help = "Full acceleration, ω = 1/√(1+2σ); no gap convergence guaranteed")]
+#[clap(
-Full
+name = "full",
+help = "Full acceleration, ω = 1/√(1+2σ); no gap convergence guaranteed"
+)]
+Full,
 }
 #[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 {
+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> {
+pub struct PDPSConfig<F: Float> {
 /// Primal step length scaling. We must have `τ0 * σ0 < 1`.
-pub τ0 : F,
+pub τ0: F,
 /// Dual step length scaling. We must have `τ0 * σ0 < 1`.
-pub σ0 : F,
+pub σ0: F,
 /// Accelerate if available
-pub acceleration : Acceleration,
+pub acceleration: Acceleration,
 /// Generic parameters
-pub generic : FBGenericConfig<F>,
+pub generic: InsertionConfig<F>,
 }
 #[replace_float_literals(F::cast_from(literal))]
-impl<F : Float> Default for PDPSConfig<F> {
+impl<F: Float> Default for PDPSConfig<F> {
 fn default() -> Self {
 let τ0 = 5.0;
 PDPSConfig {
 τ0,
-σ0 : 0.99/τ0,
+σ0: 0.99 / τ0,
-acceleration : Acceleration::Partial,
+acceleration: Acceleration::Partial,
-generic : FBGenericConfig {
+generic: InsertionConfig {
-merging : SpikeMergingMethod { enabled : true, ..Default::default() },
+merging: SpikeMergingMethod { enabled: true, ..Default::default() },
-.. Default::default()
+..Default::default()
 },
 }
 }
 }
-/// Trait for data terms for the PDPS
-#[replace_float_literals(F::cast_from(literal))]
-pub trait PDPSDataTerm<F : Float, V, const N : usize> : DataTerm<F, V, N> {
-/// Calculate some subdifferential at `x` for the conjugate
-fn some_subdifferential(&self, x : V) -> V;
-/// Factor of strong convexity of the conjugate
-#[inline]
-fn factor_of_strong_convexity(&self) -> F {
-0.0
-}
-/// Perform dual update
-fn dual_update(&self, _y : &mut V, _y_prev : &V, _σ : F);
-}
-#[replace_float_literals(F::cast_from(literal))]
-impl<F, V, const N : usize> PDPSDataTerm<F, V, N>
-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 + σ));
-}
-}
-#[replace_float_literals(F::cast_from(literal))]
-impl<F : Float + nalgebra::RealField, const N : usize>
-PDPSDataTerm<F, DVector<F>, N>
-for L1 {
-fn some_subdifferential(&self, mut x : DVector<F>) -> DVector<F> {
-// nalgebra sucks for providing second copies of the same stuff that's elsewhere as well.
-x.iter_mut()
-.for_each(|v| if *v != F::ZERO { *v = *v/<F as NumTraitsFloat>::abs(*v) });
-x
-}
-#[inline]
-fn dual_update(&self, y : &mut DVector<F>, y_prev : &DVector<F>, σ : F) {
-y.axpy(1.0, y_prev, σ);
-y.proj_ball_mut(1.0, Linfinity);
-}
-}
 /// Iteratively solve the pointsource localisation problem using primal-dual proximal splitting.
 ///
-/// The `dataterm` should be either [`L1`] for norm-1 data term or [`L2Squared`] for norm-2-squared.
 /// The settings in `config` have their [respective documentation](PDPSConfig). `opA` is the
 /// forward operator $A$, $b$ the observable, and $\lambda$ the regularisation weight.
 /// The operator `op𝒟` is used for forming the proximal term. Typically it is a convolution
 /// operator. Finally, the `iterator` is an outer loop verbosity and iteration count control
 /// as documented in [`alg_tools::iterate`].
 ///
 /// For the mathematical formulation, see the [module level](self) documentation and the manuscript.
 ///
 /// Returns the final iterate.
 #[replace_float_literals(F::cast_from(literal))]
-pub fn pointsource_pdps_reg<F, I, A, D, Reg, P, const N : usize>(
+pub fn pointsource_pdps_reg<'a, F, I, A, Phi, Reg, Plot, P, const N: usize>(
-opA : &A,
+f: &'a DataTerm<F, RNDM<N, F>, A, Phi>,
-b : &A::Observable,
+reg: &Reg,
-reg : Reg,
+prox_penalty: &P,
-prox_penalty : &P,
+pdpsconfig: &PDPSConfig<F>,
-pdpsconfig : &PDPSConfig<F>,
+iterator: I,
-iterator : I,
+mut plotter: Plot,
-mut plotter : SeqPlotter<F, N>,
+μ0 : Option<RNDM<N, F>>,
-dataterm : D,
+) -> DynResult<RNDM<N, F>>
-) -> RNDM<F, N>
 where
-F : Float + ToNalgebraRealField,
+F: Float + ToNalgebraRealField,
-I : AlgIteratorFactory<IterInfo<F, N>>,
+I: AlgIteratorFactory<IterInfo<F>>,
-A : ForwardModel<RNDM<F, N>, F>
+A: ForwardModel<RNDM<N, F>, F>,
-+ AdjointProductBoundedBy<RNDM<F, N>, P, FloatType=F>,
+for<'b> &'b A::Observable: Instance<A::Observable>,
-A::PreadjointCodomain : RealMapping<F, N>,
+A::Observable: AXPY<Field = F>,
-for<'b> &'b A::Observable : std::ops::Neg<Output=A::Observable> + Instance<A::Observable>,
+RNDM<N, F>: SpikeMerging<F>,
-PlotLookup : Plotting<N>,
+Reg: RegTerm<Loc<N, F>, F>,
-RNDM<F, N> : SpikeMerging<F>,
+Phi: Conjugable<A::Observable, F>,
-D : PDPSDataTerm<F, A::Observable, N>,
+for<'b> Phi::Conjugate<'b>: Prox<A::Observable>,
-Reg : RegTerm<F, N>,
+P: ProxPenalty<Loc<N, F>, A::PreadjointCodomain, Reg, F> + StepLengthBoundPD<F, A, RNDM<N, F>>,
-P : ProxPenalty<F, A::PreadjointCodomain, Reg, N>,
+Plot: Plotter<P::ReturnMapping, A::PreadjointCodomain, RNDM<N, F>>,
 {
 // Check parameters
-assert!(pdpsconfig.τ0 > 0.0 &&
+ensure!(
-pdpsconfig.σ0 > 0.0 &&
+pdpsconfig.τ0 > 0.0 && pdpsconfig.σ0 > 0.0 && pdpsconfig.τ0 * pdpsconfig.σ0 <= 1.0,
-pdpsconfig.τ0 * pdpsconfig.σ0 <= 1.0,
+"Invalid step length parameters"
-"Invalid step length parameters");
+);
+let opA = f.operator();
+let b = f.data();
+let phistar = f.fidelity().conjugate();
 // Set up parameters
 let config = &pdpsconfig.generic;
-let l = opA.adjoint_product_bound(prox_penalty).unwrap().sqrt();
+let l = prox_penalty.step_length_bound_pd(opA)?;
 let mut τ = pdpsconfig.τ0 / l;
 let mut σ = pdpsconfig.σ0 / l;
-let γ = dataterm.factor_of_strong_convexity();
+let γ = phistar.factor_of_strong_convexity();
 // We multiply tolerance by τ for FB since our subproblems depending on tolerances are scaled
 // by τ compared to the conditional gradient approach.
 let tolerance = config.tolerance * τ * reg.tolerance_scaling();
 let mut ε = tolerance.initial();
 // Initialise iterates
-let mut μ = DiscreteMeasure::new();
+let mut μ = μ0.unwrap_or_else(|| DiscreteMeasure::new());
-let mut y = dataterm.some_subdifferential(-b);
+let mut y = f.residual(&μ);
-let mut y_prev = y.clone();
+let full_stats = |μ: &RNDM<N, F>, ε, stats| IterInfo {
-let full_stats = |μ : &RNDM<F, N>, ε, stats| IterInfo {
+value: f.apply(μ) + reg.apply(μ),
-value : dataterm.calculate_fit_op(μ, opA, b) + reg.apply(μ),
+n_spikes: μ.len(),
-n_spikes : μ.len(),
 ε,
 // postprocessing: config.postprocessing.then(|| μ.clone()),
-.. stats
+..stats
 };
 let mut stats = IterInfo::new();
 // Run the algorithm
 for state in iterator.iter_init(|| full_stats(&μ, ε, stats.clone())) {
 // Calculate smooth part of surrogate model.
-let mut τv = opA.preadjoint().apply(y * τ);
+// FIXME: the clone is required to avoid compiler overflows with reference-Mul requirement above.
+let mut τv = opA.preadjoint().apply(y.clone() * τ);
 // Save current base point
 let μ_base = μ.clone();
 // Insert and reweigh
 let (maybe_d, _within_tolerances) = prox_penalty.insert_and_reweigh(
-&mut μ, &mut τv, &μ_base, None,
+&mut μ, &mut τv, &μ_base, None, τ, ε, config, &reg, &state, &mut stats,
-τ, ε,
+)?;
-config, &reg, &state, &mut stats
-);
 // Prune and possibly merge spikes
 if config.merge_now(&state) {
-stats.merged += prox_penalty.merge_spikes_no_fitness(
+stats.merged += prox_penalty
-&mut μ, &mut τv, &μ_base, None, τ, ε, config, &reg,
+.merge_spikes_no_fitness(&mut μ, &mut τv, &μ_base, None, τ, ε, config, &reg);
-);
 }
 stats.pruned += prune_with_stats(&mut μ);
 // Update step length parameters
 let ω = pdpsconfig.acceleration.accelerate(&mut τ, &mut σ, γ);
 // Do dual update
-y = b.clone();                          // y = b
+// y = y_prev + τb
-opA.gemv(&mut y, 1.0 + ω, &μ, -1.0);    // y = A[(1+ω)μ^{k+1}]-b
+y.axpy(τ, b, 1.0);
-opA.gemv(&mut y, -ω, &μ_base, 1.0);     // y = A[(1+ω)μ^{k+1} - ω μ^k]-b
+// y = y_prev - τ(A[(1+ω)μ^{k+1}]-b)
-dataterm.dual_update(&mut y, &y_prev, σ);
+opA.gemv(&mut y, -τ * (1.0 + ω), &μ, 1.0);
-y_prev.copy_from(&y);
+// y = y_prev - τ(A[(1+ω)μ^{k+1} - ω μ^k]-b)
+opA.gemv(&mut y, τ * ω, &μ_base, 1.0);
+y = phistar.prox(τ, y);
 // Give statistics if requested
 let iter = state.iteration();
 stats.this_iters += 1;
 });
 ε = tolerance.update(ε, iter);
 }
-postprocess(μ, config, dataterm, opA, b)
+postprocess(μ, config, |μ| f.apply(μ))
 }

Mercurial > repos > pointsource_algs / file comparison

comparison: src/pdps.rs

src/pdps.rs