--- a/src/fb.rs Tue Aug 01 10:25:09 2023 +0300
+++ b/src/fb.rs Mon Feb 17 13:54:53 2025 -0500
@@ -6,10 +6,7 @@
* Valkonen T. - _Proximal methods for point source localisation_,
[arXiv:2212.02991](https://arxiv.org/abs/2212.02991).
-The main routine is [`pointsource_fb_reg`]. It is based on [`generic_pointsource_fb_reg`], which is
-also used by our [primal-dual proximal splitting][crate::pdps] implementation.
-
-FISTA-type inertia can also be enabled through [`FBConfig::meta`].
+The main routine is [`pointsource_fb_reg`].
## Problem
@@ -76,650 +73,94 @@
$$
</p>
-We solve this with either SSN or FB via [`quadratic_nonneg`] as determined by
-[`InnerSettings`] in [`FBGenericConfig::inner`].
+We solve this with either SSN or FB as determined by
+[`crate::subproblem::InnerSettings`] in [`FBGenericConfig::inner`].
*/
+use colored::Colorize;
use numeric_literals::replace_float_literals;
-use serde::{Serialize, Deserialize};
-use colored::Colorize;
-use nalgebra::{DVector, DMatrix};
+use serde::{Deserialize, Serialize};
-use alg_tools::iterate::{
- AlgIteratorFactory,
- AlgIteratorState,
-};
use alg_tools::euclidean::Euclidean;
-use alg_tools::linops::Apply;
-use alg_tools::sets::Cube;
-use alg_tools::loc::Loc;
-use alg_tools::mapping::Mapping;
-use alg_tools::bisection_tree::{
- BTFN,
- PreBTFN,
- Bounds,
- BTNodeLookup,
- BTNode,
- BTSearch,
- P2Minimise,
- SupportGenerator,
- LocalAnalysis,
- Bounded,
-};
+use alg_tools::instance::Instance;
+use alg_tools::iterate::AlgIteratorFactory;
+use alg_tools::linops::{Mapping, GEMV};
use alg_tools::mapping::RealMapping;
use alg_tools::nalgebra_support::ToNalgebraRealField;
+use crate::dataterm::{calculate_residual, DataTerm, L2Squared};
+use crate::forward_model::{AdjointProductBoundedBy, ForwardModel};
+use crate::measures::merging::SpikeMerging;
+use crate::measures::{DiscreteMeasure, RNDM};
+use crate::plot::{PlotLookup, Plotting, SeqPlotter};
+pub use crate::prox_penalty::{FBGenericConfig, ProxPenalty};
+use crate::regularisation::RegTerm;
use crate::types::*;
-use crate::measures::{
- DiscreteMeasure,
- DeltaMeasure,
-};
-use crate::measures::merging::{
- SpikeMergingMethod,
- SpikeMerging,
-};
-use crate::forward_model::ForwardModel;
-use crate::seminorms::{
- DiscreteMeasureOp, Lipschitz
-};
-use crate::subproblem::{
- nonneg::quadratic_nonneg,
- unconstrained::quadratic_unconstrained,
- InnerSettings,
- InnerMethod,
-};
-use crate::tolerance::Tolerance;
-use crate::plot::{
- SeqPlotter,
- Plotting,
- PlotLookup
-};
-use crate::regularisation::{
- NonnegRadonRegTerm,
- RadonRegTerm,
-};
-
-/// Method for constructing $μ$ on each iteration
-#[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
-#[allow(dead_code)]
-pub enum InsertionStyle {
- /// Resuse previous $μ$ from previous iteration, optimising weights
- /// before inserting new spikes.
- Reuse,
- /// Start each iteration with $μ=0$.
- Zero,
-}
-
-/// Meta-algorithm type
-#[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
-#[allow(dead_code)]
-pub enum FBMetaAlgorithm {
- /// No meta-algorithm
- None,
- /// FISTA-style inertia
- InertiaFISTA,
-}
/// Settings for [`pointsource_fb_reg`].
#[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
#[serde(default)]
-pub struct FBConfig<F : Float> {
+pub struct FBConfig<F: Float> {
/// Step length scaling
- pub τ0 : F,
- /// Meta-algorithm to apply
- pub meta : FBMetaAlgorithm,
+ pub τ0: F,
/// Generic parameters
- pub insertion : FBGenericConfig<F>,
-}
-
-/// Settings for the solution of the stepwise optimality condition in algorithms based on
-/// [`generic_pointsource_fb_reg`].
-#[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
-#[serde(default)]
-pub struct FBGenericConfig<F : Float> {
- /// Method for constructing $μ$ on each iteration; see [`InsertionStyle`].
- pub insertion_style : InsertionStyle,
- /// Tolerance for point insertion.
- pub tolerance : Tolerance<F>,
- /// Stop looking for predual maximum (where to isert a new point) below
- /// `tolerance` multiplied by this factor.
- pub insertion_cutoff_factor : F,
- /// Settings for branch and bound refinement when looking for predual maxima
- pub refinement : RefinementSettings<F>,
- /// Maximum insertions within each outer iteration
- pub max_insertions : usize,
- /// Pair `(n, m)` for maximum insertions `m` on first `n` iterations.
- pub bootstrap_insertions : Option<(usize, usize)>,
- /// Inner method settings
- pub inner : InnerSettings<F>,
- /// Spike merging method
- pub merging : SpikeMergingMethod<F>,
- /// Tolerance multiplier for merges
- pub merge_tolerance_mult : F,
- /// Spike merging method after the last step
- pub final_merging : SpikeMergingMethod<F>,
- /// Iterations between merging heuristic tries
- pub merge_every : usize,
- /// Save $μ$ for postprocessing optimisation
- pub postprocessing : bool
+ pub generic: FBGenericConfig<F>,
}
#[replace_float_literals(F::cast_from(literal))]
-impl<F : Float> Default for FBConfig<F> {
+impl<F: Float> Default for FBConfig<F> {
fn default() -> Self {
FBConfig {
- τ0 : 0.99,
- meta : FBMetaAlgorithm::None,
- insertion : Default::default()
- }
- }
-}
-
-#[replace_float_literals(F::cast_from(literal))]
-impl<F : Float> Default for FBGenericConfig<F> {
- fn default() -> Self {
- FBGenericConfig {
- insertion_style : InsertionStyle::Reuse,
- tolerance : Default::default(),
- insertion_cutoff_factor : 1.0,
- refinement : Default::default(),
- max_insertions : 100,
- //bootstrap_insertions : None,
- bootstrap_insertions : Some((10, 1)),
- inner : InnerSettings {
- method : InnerMethod::SSN,
- .. Default::default()
- },
- merging : SpikeMergingMethod::None,
- //merging : Default::default(),
- final_merging : Default::default(),
- merge_every : 10,
- merge_tolerance_mult : 2.0,
- postprocessing : false,
+ τ0: 0.99,
+ generic: Default::default(),
}
}
}
-/// Trait for specialisation of [`generic_pointsource_fb_reg`] to basic FB, FISTA.
-///
-/// The idea is that the residual $Aμ - b$ in the forward step can be replaced by an arbitrary
-/// value. For example, to implement [primal-dual proximal splitting][crate::pdps] we replace it
-/// with the dual variable $y$. We can then also implement alternative data terms, as the
-/// (pre)differential of $F(μ)=F\_0(Aμ-b)$ is $F\'(μ) = A\_*F\_0\'(Aμ-b)$. In the case of the
-/// quadratic fidelity $F_0(y)=\frac{1}{2}\\|y\\|_2^2$ in a Hilbert space, of course,
-/// $F\_0\'(Aμ-b)=Aμ-b$ is the residual.
-pub trait FBSpecialisation<F : Float, Observable : Euclidean<F>, const N : usize> : Sized {
- /// Updates the residual and does any necessary pruning of `μ`.
- ///
- /// Returns the new residual and possibly a new step length.
- ///
- /// The measure `μ` may also be modified to apply, e.g., inertia to it.
- /// The updated residual should correspond to the residual at `μ`.
- /// See the [trait documentation][FBSpecialisation] for the use and meaning of the residual.
- ///
- /// The parameter `μ_base` is the base point of the iteration, typically the previous iterate,
- /// but for, e.g., FISTA has inertia applied to it.
- fn update(
- &mut self,
- μ : &mut DiscreteMeasure<Loc<F, N>, F>,
- μ_base : &DiscreteMeasure<Loc<F, N>, F>,
- ) -> (Observable, Option<F>);
-
- /// Calculates the data term value corresponding to iterate `μ` and available residual.
- ///
- /// Inertia and other modifications, as deemed, necessary, should be applied to `μ`.
- ///
- /// The blanket implementation correspondsn to the 2-norm-squared data fidelity
- /// $\\|\text{residual}\\|\_2^2/2$.
- fn calculate_fit(
- &self,
- _μ : &DiscreteMeasure<Loc<F, N>, F>,
- residual : &Observable
- ) -> F {
- residual.norm2_squared_div2()
- }
-
- /// Calculates the data term value at $μ$.
- ///
- /// Unlike [`Self::calculate_fit`], no inertia, etc., should be applied to `μ`.
- fn calculate_fit_simple(
- &self,
- μ : &DiscreteMeasure<Loc<F, N>, F>,
- ) -> F;
-
- /// Returns the final iterate after any necessary postprocess pruning, merging, etc.
- fn postprocess(self, mut μ : DiscreteMeasure<Loc<F, N>, F>, merging : SpikeMergingMethod<F>)
- -> DiscreteMeasure<Loc<F, N>, F>
- where DiscreteMeasure<Loc<F, N>, F> : SpikeMerging<F> {
- μ.merge_spikes_fitness(merging,
- |μ̃| self.calculate_fit_simple(μ̃),
- |&v| v);
- μ.prune();
- μ
- }
-
- /// Returns measure to be used for value calculations, which may differ from μ.
- fn value_μ<'c, 'b : 'c>(&'b self, μ : &'c DiscreteMeasure<Loc<F, N>, F>)
- -> &'c DiscreteMeasure<Loc<F, N>, F> {
- μ
- }
-}
-
-/// Specialisation of [`generic_pointsource_fb_reg`] to basic μFB.
-struct BasicFB<
- 'a,
- F : Float + ToNalgebraRealField,
- A : ForwardModel<Loc<F, N>, F>,
- const N : usize
-> {
- /// The data
- b : &'a A::Observable,
- /// The forward operator
- opA : &'a A,
-}
-
-/// Implementation of [`FBSpecialisation`] for basic μFB forward-backward splitting.
-#[replace_float_literals(F::cast_from(literal))]
-impl<'a, F : Float + ToNalgebraRealField , A : ForwardModel<Loc<F, N>, F>, const N : usize>
-FBSpecialisation<F, A::Observable, N> for BasicFB<'a, F, A, N> {
- fn update(
- &mut self,
- μ : &mut DiscreteMeasure<Loc<F, N>, F>,
- _μ_base : &DiscreteMeasure<Loc<F, N>, F>
- ) -> (A::Observable, Option<F>) {
- μ.prune();
- //*residual = self.opA.apply(μ) - self.b;
- let mut residual = self.b.clone();
- self.opA.gemv(&mut residual, 1.0, μ, -1.0);
- (residual, None)
- }
-
- fn calculate_fit_simple(
- &self,
- μ : &DiscreteMeasure<Loc<F, N>, F>,
- ) -> F {
- let mut residual = self.b.clone();
- self.opA.gemv(&mut residual, 1.0, μ, -1.0);
- residual.norm2_squared_div2()
- }
-}
-
-/// Specialisation of [`generic_pointsource_fb_reg`] to FISTA.
-struct FISTA<
- 'a,
- F : Float + ToNalgebraRealField,
- A : ForwardModel<Loc<F, N>, F>,
- const N : usize
-> {
- /// The data
- b : &'a A::Observable,
- /// The forward operator
- opA : &'a A,
- /// Current inertial parameter
- λ : F,
- /// Previous iterate without inertia applied.
- /// We need to store this here because `μ_base` passed to [`FBSpecialisation::update`] will
- /// have inertia applied to it, so is not useful to use.
- μ_prev : DiscreteMeasure<Loc<F, N>, F>,
-}
-
-/// Implementation of [`FBSpecialisation`] for μFISTA inertial forward-backward splitting.
-#[replace_float_literals(F::cast_from(literal))]
-impl<'a, F : Float + ToNalgebraRealField, A : ForwardModel<Loc<F, N>, F>, const N : usize>
-FBSpecialisation<F, A::Observable, N> for FISTA<'a, F, A, N> {
- fn update(
- &mut self,
- μ : &mut DiscreteMeasure<Loc<F, N>, F>,
- _μ_base : &DiscreteMeasure<Loc<F, N>, F>
- ) -> (A::Observable, Option<F>) {
- // Update inertial parameters
- let λ_prev = self.λ;
- self.λ = 2.0 * λ_prev / ( λ_prev + (4.0 + λ_prev * λ_prev).sqrt() );
- let θ = self.λ / λ_prev - self.λ;
- // Perform inertial update on μ.
- // This computes μ ← (1 + θ) * μ - θ * μ_prev, pruning spikes where both μ
- // and μ_prev have zero weight. Since both have weights from the finite-dimensional
- // subproblem with a proximal projection step, this is likely to happen when the
- // spike is not needed. A copy of the pruned μ without artithmetic performed is
- // stored in μ_prev.
- μ.pruning_sub(1.0 + θ, θ, &mut self.μ_prev);
-
- //*residual = self.opA.apply(μ) - self.b;
- let mut residual = self.b.clone();
- self.opA.gemv(&mut residual, 1.0, μ, -1.0);
- (residual, None)
- }
-
- fn calculate_fit_simple(
- &self,
- μ : &DiscreteMeasure<Loc<F, N>, F>,
- ) -> F {
- let mut residual = self.b.clone();
- self.opA.gemv(&mut residual, 1.0, μ, -1.0);
- residual.norm2_squared_div2()
- }
-
- fn calculate_fit(
- &self,
- _μ : &DiscreteMeasure<Loc<F, N>, F>,
- _residual : &A::Observable
- ) -> F {
- self.calculate_fit_simple(&self.μ_prev)
- }
-
- // For FISTA we need to do a final pruning as well, due to the limited
- // pruning that can be done on each step.
- fn postprocess(mut self, μ_base : DiscreteMeasure<Loc<F, N>, F>, merging : SpikeMergingMethod<F>)
- -> DiscreteMeasure<Loc<F, N>, F>
- where DiscreteMeasure<Loc<F, N>, F> : SpikeMerging<F> {
- let mut μ = self.μ_prev;
- self.μ_prev = μ_base;
- μ.merge_spikes_fitness(merging,
- |μ̃| self.calculate_fit_simple(μ̃),
- |&v| v);
- μ.prune();
- μ
- }
-
- fn value_μ<'c, 'b : 'c>(&'c self, _μ : &'c DiscreteMeasure<Loc<F, N>, F>)
- -> &'c DiscreteMeasure<Loc<F, N>, F> {
- &self.μ_prev
- }
-}
-
-
-/// Abstraction of regularisation terms for [`generic_pointsource_fb_reg`].
-pub trait RegTerm<F : Float + ToNalgebraRealField, const N : usize>
-: for<'a> Apply<&'a DiscreteMeasure<Loc<F, N>, F>, Output = F> {
- /// Approximately solve the problem
- /// <div>$$
- /// \min_{x ∈ ℝ^n} \frac{1}{2} x^⊤Ax - g^⊤ x + τ G(x)
- /// $$</div>
- /// for $G$ depending on the trait implementation.
- ///
- /// The parameter `mA` is $A$. An estimate for its opeator norm should be provided in
- /// `mA_normest`. The initial iterate and output is `x`. The current main tolerance is `ε`.
- ///
- /// Returns the number of iterations taken.
- fn solve_findim(
- &self,
- mA : &DMatrix<F::MixedType>,
- g : &DVector<F::MixedType>,
- τ : F,
- x : &mut DVector<F::MixedType>,
- mA_normest : F,
- ε : F,
- config : &FBGenericConfig<F>
- ) -> usize;
-
- /// Find a point where `d` may violate the tolerance `ε`.
- ///
- /// If `skip_by_rough_check` is set, do not find the point if a rough check indicates that we
- /// are in bounds. `ε` is the current main tolerance and `τ` a scaling factor for the
- /// regulariser.
- ///
- /// Returns `None` if `d` is in bounds either based on the rough check, or a more precise check
- /// terminating early. Otherwise returns a possibly violating point, the value of `d` there,
- /// and a boolean indicating whether the found point is in bounds.
- fn find_tolerance_violation<G, BT>(
- &self,
- d : &mut BTFN<F, G, BT, N>,
- τ : F,
- ε : F,
- skip_by_rough_check : bool,
- config : &FBGenericConfig<F>,
- ) -> Option<(Loc<F, N>, F, bool)>
- where BT : BTSearch<F, N, Agg=Bounds<F>>,
- G : SupportGenerator<F, N, Id=BT::Data>,
- G::SupportType : Mapping<Loc<F, N>,Codomain=F>
- + LocalAnalysis<F, Bounds<F>, N>;
-
- /// Verify that `d` is in bounds `ε` for a merge candidate `μ`
- ///
- /// `ε` is the current main tolerance and `τ` a scaling factor for the regulariser.
- fn verify_merge_candidate<G, BT>(
- &self,
- d : &mut BTFN<F, G, BT, N>,
- μ : &DiscreteMeasure<Loc<F, N>, F>,
- τ : F,
- ε : F,
- config : &FBGenericConfig<F>,
- ) -> bool
- where BT : BTSearch<F, N, Agg=Bounds<F>>,
- G : SupportGenerator<F, N, Id=BT::Data>,
- G::SupportType : Mapping<Loc<F, N>,Codomain=F>
- + LocalAnalysis<F, Bounds<F>, N>;
-
- fn target_bounds(&self, τ : F, ε : F) -> Option<Bounds<F>>;
-
- /// Returns a scaling factor for the tolerance sequence.
- ///
- /// Typically this is the regularisation parameter.
- fn tolerance_scaling(&self) -> F;
+pub(crate) fn prune_with_stats<F: Float, const N: usize>(μ: &mut RNDM<F, N>) -> usize {
+ let n_before_prune = μ.len();
+ μ.prune();
+ debug_assert!(μ.len() <= n_before_prune);
+ n_before_prune - μ.len()
}
#[replace_float_literals(F::cast_from(literal))]
-impl<F : Float + ToNalgebraRealField, const N : usize> RegTerm<F, N> for NonnegRadonRegTerm<F>
-where Cube<F, N> : P2Minimise<Loc<F, N>, F> {
- fn solve_findim(
- &self,
- mA : &DMatrix<F::MixedType>,
- g : &DVector<F::MixedType>,
- τ : F,
- x : &mut DVector<F::MixedType>,
- mA_normest : F,
- ε : F,
- config : &FBGenericConfig<F>
- ) -> usize {
- let inner_tolerance = ε * config.inner.tolerance_mult;
- let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
- let inner_τ = config.inner.τ0 / mA_normest;
- quadratic_nonneg(config.inner.method, mA, g, τ * self.α(), x,
- inner_τ, inner_it)
- }
-
- #[inline]
- fn find_tolerance_violation<G, BT>(
- &self,
- d : &mut BTFN<F, G, BT, N>,
- τ : F,
- ε : F,
- skip_by_rough_check : bool,
- config : &FBGenericConfig<F>,
- ) -> Option<(Loc<F, N>, F, bool)>
- where BT : BTSearch<F, N, Agg=Bounds<F>>,
- G : SupportGenerator<F, N, Id=BT::Data>,
- G::SupportType : Mapping<Loc<F, N>,Codomain=F>
- + LocalAnalysis<F, Bounds<F>, N> {
- let τα = τ * self.α();
- let keep_below = τα + ε;
- let maximise_above = τα + ε * config.insertion_cutoff_factor;
- let refinement_tolerance = ε * config.refinement.tolerance_mult;
-
- // If preliminary check indicates that we are in bonds, and if it otherwise matches
- // the insertion strategy, skip insertion.
- if skip_by_rough_check && d.bounds().upper() <= keep_below {
- None
- } else {
- // If the rough check didn't indicate no insertion needed, find maximising point.
- d.maximise_above(maximise_above, refinement_tolerance, config.refinement.max_steps)
- .map(|(ξ, v_ξ)| (ξ, v_ξ, v_ξ <= keep_below))
- }
- }
-
- fn verify_merge_candidate<G, BT>(
- &self,
- d : &mut BTFN<F, G, BT, N>,
- μ : &DiscreteMeasure<Loc<F, N>, F>,
- τ : F,
- ε : F,
- config : &FBGenericConfig<F>,
- ) -> bool
- where BT : BTSearch<F, N, Agg=Bounds<F>>,
- G : SupportGenerator<F, N, Id=BT::Data>,
- G::SupportType : Mapping<Loc<F, N>,Codomain=F>
- + LocalAnalysis<F, Bounds<F>, N> {
- let τα = τ * self.α();
- let refinement_tolerance = ε * config.refinement.tolerance_mult;
- let merge_tolerance = config.merge_tolerance_mult * ε;
- let keep_below = τα + merge_tolerance;
- let keep_supp_above = τα - merge_tolerance;
- let bnd = d.bounds();
-
- return (
- bnd.lower() >= keep_supp_above
- ||
- μ.iter_spikes().map(|&DeltaMeasure{ α : β, ref x }| {
- (β == 0.0) || d.apply(x) >= keep_supp_above
- }).all(std::convert::identity)
- ) && (
- bnd.upper() <= keep_below
- ||
- d.has_upper_bound(keep_below, refinement_tolerance, config.refinement.max_steps)
- )
- }
-
- fn target_bounds(&self, τ : F, ε : F) -> Option<Bounds<F>> {
- let τα = τ * self.α();
- Some(Bounds(τα - ε, τα + ε))
- }
-
- fn tolerance_scaling(&self) -> F {
- self.α()
- }
+pub(crate) fn postprocess<
+ F: Float,
+ V: Euclidean<F> + Clone,
+ A: GEMV<F, RNDM<F, N>, Codomain = V>,
+ D: DataTerm<F, V, N>,
+ const N: usize,
+>(
+ mut μ: RNDM<F, N>,
+ config: &FBGenericConfig<F>,
+ dataterm: D,
+ opA: &A,
+ b: &V,
+) -> RNDM<F, N>
+where
+ RNDM<F, N>: SpikeMerging<F>,
+ for<'a> &'a RNDM<F, N>: Instance<RNDM<F, N>>,
+{
+ μ.merge_spikes_fitness(
+ config.final_merging_method(),
+ |μ̃| dataterm.calculate_fit_op(μ̃, opA, b),
+ |&v| v,
+ );
+ μ.prune();
+ μ
}
-#[replace_float_literals(F::cast_from(literal))]
-impl<F : Float + ToNalgebraRealField, const N : usize> RegTerm<F, N> for RadonRegTerm<F>
-where Cube<F, N> : P2Minimise<Loc<F, N>, F> {
- fn solve_findim(
- &self,
- mA : &DMatrix<F::MixedType>,
- g : &DVector<F::MixedType>,
- τ : F,
- x : &mut DVector<F::MixedType>,
- mA_normest: F,
- ε : F,
- config : &FBGenericConfig<F>
- ) -> usize {
- let inner_tolerance = ε * config.inner.tolerance_mult;
- let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
- let inner_τ = config.inner.τ0 / mA_normest;
- quadratic_unconstrained(config.inner.method, mA, g, τ * self.α(), x,
- inner_τ, inner_it)
- }
-
- fn find_tolerance_violation<G, BT>(
- &self,
- d : &mut BTFN<F, G, BT, N>,
- τ : F,
- ε : F,
- skip_by_rough_check : bool,
- config : &FBGenericConfig<F>,
- ) -> Option<(Loc<F, N>, F, bool)>
- where BT : BTSearch<F, N, Agg=Bounds<F>>,
- G : SupportGenerator<F, N, Id=BT::Data>,
- G::SupportType : Mapping<Loc<F, N>,Codomain=F>
- + LocalAnalysis<F, Bounds<F>, N> {
- let τα = τ * self.α();
- let keep_below = τα + ε;
- let keep_above = -τα - ε;
- let maximise_above = τα + ε * config.insertion_cutoff_factor;
- let minimise_below = -τα - ε * config.insertion_cutoff_factor;
- let refinement_tolerance = ε * config.refinement.tolerance_mult;
-
- // If preliminary check indicates that we are in bonds, and if it otherwise matches
- // the insertion strategy, skip insertion.
- if skip_by_rough_check && Bounds(keep_above, keep_below).superset(&d.bounds()) {
- None
- } else {
- // If the rough check didn't indicate no insertion needed, find maximising point.
- let mx = d.maximise_above(maximise_above, refinement_tolerance,
- config.refinement.max_steps);
- let mi = d.minimise_below(minimise_below, refinement_tolerance,
- config.refinement.max_steps);
-
- match (mx, mi) {
- (None, None) => None,
- (Some((ξ, v_ξ)), None) => Some((ξ, v_ξ, keep_below >= v_ξ)),
- (None, Some((ζ, v_ζ))) => Some((ζ, v_ζ, keep_above <= v_ζ)),
- (Some((ξ, v_ξ)), Some((ζ, v_ζ))) => {
- if v_ξ - τα > τα - v_ζ {
- Some((ξ, v_ξ, keep_below >= v_ξ))
- } else {
- Some((ζ, v_ζ, keep_above <= v_ζ))
- }
- }
- }
- }
- }
-
- fn verify_merge_candidate<G, BT>(
- &self,
- d : &mut BTFN<F, G, BT, N>,
- μ : &DiscreteMeasure<Loc<F, N>, F>,
- τ : F,
- ε : F,
- config : &FBGenericConfig<F>,
- ) -> bool
- where BT : BTSearch<F, N, Agg=Bounds<F>>,
- G : SupportGenerator<F, N, Id=BT::Data>,
- G::SupportType : Mapping<Loc<F, N>,Codomain=F>
- + LocalAnalysis<F, Bounds<F>, N> {
- let τα = τ * self.α();
- let refinement_tolerance = ε * config.refinement.tolerance_mult;
- let merge_tolerance = config.merge_tolerance_mult * ε;
- let keep_below = τα + merge_tolerance;
- let keep_above = -τα - merge_tolerance;
- let keep_supp_pos_above = τα - merge_tolerance;
- let keep_supp_neg_below = -τα + merge_tolerance;
- let bnd = d.bounds();
-
- return (
- (bnd.lower() >= keep_supp_pos_above && bnd.upper() <= keep_supp_neg_below)
- ||
- μ.iter_spikes().map(|&DeltaMeasure{ α : β, ref x }| {
- use std::cmp::Ordering::*;
- match β.partial_cmp(&0.0) {
- Some(Greater) => d.apply(x) >= keep_supp_pos_above,
- Some(Less) => d.apply(x) <= keep_supp_neg_below,
- _ => true,
- }
- }).all(std::convert::identity)
- ) && (
- bnd.upper() <= keep_below
- ||
- d.has_upper_bound(keep_below, refinement_tolerance,
- config.refinement.max_steps)
- ) && (
- bnd.lower() >= keep_above
- ||
- d.has_lower_bound(keep_above, refinement_tolerance,
- config.refinement.max_steps)
- )
- }
-
- fn target_bounds(&self, τ : F, ε : F) -> Option<Bounds<F>> {
- let τα = τ * self.α();
- Some(Bounds(-τα - ε, τα + ε))
- }
-
- fn tolerance_scaling(&self) -> F {
- self.α()
- }
-}
-
-
-/// Generic implementation of [`pointsource_fb_reg`].
+/// Iteratively solve the pointsource localisation problem using forward-backward splitting.
///
-/// The method can be specialised to even primal-dual proximal splitting through the
-/// [`FBSpecialisation`] parameter `specialisation`.
-/// The settings in `config` have their [respective documentation](FBGenericConfig). `opA` is the
+/// The settings in `config` have their [respective documentation](FBConfig). `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 details on the mathematical formulation, see the [module level](self) documentation.
+///
/// The implementation relies on [`alg_tools::bisection_tree::BTFN`] presentations of
/// sums of simple functions usign bisection trees, and the related
/// [`alg_tools::bisection_tree::Aggregator`]s, to efficiently search for component functions
@@ -729,233 +170,103 @@
///
/// Returns the final iterate.
#[replace_float_literals(F::cast_from(literal))]
-pub fn generic_pointsource_fb_reg<
- 'a, F, I, A, GA, 𝒟, BTA, G𝒟, S, K, Spec, Reg, const N : usize
->(
- opA : &'a A,
- reg : Reg,
- op𝒟 : &'a 𝒟,
- mut τ : F,
- config : &FBGenericConfig<F>,
- iterator : I,
- mut plotter : SeqPlotter<F, N>,
- mut residual : A::Observable,
- mut specialisation : Spec
-) -> DiscreteMeasure<Loc<F, N>, F>
-where F : Float + ToNalgebraRealField,
- I : AlgIteratorFactory<IterInfo<F, N>>,
- Spec : FBSpecialisation<F, A::Observable, N>,
- 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>>
- + Lipschitz<𝒟, 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>,
- Reg : RegTerm<F, N> {
-
+pub fn pointsource_fb_reg<F, I, A, Reg, P, const N: usize>(
+ opA: &A,
+ b: &A::Observable,
+ reg: Reg,
+ prox_penalty: &P,
+ fbconfig: &FBConfig<F>,
+ iterator: I,
+ mut plotter: SeqPlotter<F, N>,
+) -> RNDM<F, N>
+where
+ F: Float + ToNalgebraRealField,
+ I: AlgIteratorFactory<IterInfo<F, N>>,
+ for<'b> &'b A::Observable: std::ops::Neg<Output = A::Observable>,
+ A: ForwardModel<RNDM<F, N>, F> + AdjointProductBoundedBy<RNDM<F, N>, P, FloatType = F>,
+ A::PreadjointCodomain: RealMapping<F, N>,
+ PlotLookup: Plotting<N>,
+ RNDM<F, N>: SpikeMerging<F>,
+ Reg: RegTerm<F, N>,
+ P: ProxPenalty<F, A::PreadjointCodomain, Reg, N>,
+{
// Set up parameters
- let quiet = iterator.is_quiet();
- let op𝒟norm = op𝒟.opnorm_bound();
+ let config = &fbconfig.generic;
+ let τ = fbconfig.τ0 / opA.adjoint_product_bound(prox_penalty).unwrap();
// 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 operators
- let preadjA = opA.preadjoint();
-
// Initialise iterates
let mut μ = DiscreteMeasure::new();
-
- let mut inner_iters = 0;
- let mut this_iters = 0;
- let mut pruned = 0;
- let mut merged = 0;
+ let mut residual = -b;
- let μ_diff = |μ_new : &DiscreteMeasure<Loc<F, N>, F>,
- μ_base : &DiscreteMeasure<Loc<F, N>, F>| {
- let mut ν : DiscreteMeasure<Loc<F, N>, F> = match config.insertion_style {
- InsertionStyle::Reuse => {
- μ_new.iter_spikes()
- .zip(μ_base.iter_masses().chain(std::iter::repeat(0.0)))
- .map(|(δ, α_base)| (δ.x, α_base - δ.α))
- .collect()
- },
- InsertionStyle::Zero => {
- μ_new.iter_spikes()
- .map(|δ| -δ)
- .chain(μ_base.iter_spikes().copied())
- .collect()
- }
- };
- ν.prune(); // Potential small performance improvement
- ν
+ // Statistics
+ let full_stats = |residual: &A::Observable, μ: &RNDM<F, N>, ε, stats| IterInfo {
+ value: residual.norm2_squared_div2() + reg.apply(μ),
+ n_spikes: μ.len(),
+ ε,
+ //postprocessing: config.postprocessing.then(|| μ.clone()),
+ ..stats
};
+ let mut stats = IterInfo::new();
// Run the algorithm
- iterator.iterate(|state| {
- // Maximum insertion count and measure difference calculation depend on insertion style.
- let (m, warn_insertions) = match (state.iteration(), config.bootstrap_insertions) {
- (i, Some((l, k))) if i <= l => (k, false),
- _ => (config.max_insertions, !quiet),
- };
- let max_insertions = match config.insertion_style {
- InsertionStyle::Zero => {
- todo!("InsertionStyle::Zero does not currently work with FISTA, so diabled.");
- // let n = μ.len();
- // μ = DiscreteMeasure::new();
- // n + m
- },
- InsertionStyle::Reuse => m,
- };
-
+ for state in iterator.iter_init(|| full_stats(&residual, &μ, ε, 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.
- residual *= -τ;
- let r = std::mem::replace(&mut residual, opA.empty_observable());
- let minus_τv = preadjA.apply(r); // minus_τv = -τA^*(Aμ^k-b)
- // TODO: should avoid a second copy of μ here; μ_base already stores a copy.
- let ω0 = op𝒟.apply(μ.clone()); // 𝒟μ^k
- //let g = &minus_τv + ω0; // Linear term of surrogate model
+ let mut τv = opA.preadjoint().apply(residual * τ);
// Save current base point
let μ_base = μ.clone();
-
- // Add points to support until within error tolerance or maximum insertion count reached.
- let mut count = 0;
- let (within_tolerances, d) = 'insertion: loop {
- if μ.len() > 0 {
- // Form finite-dimensional subproblem. The subproblem references to the original μ^k
- // from the beginning of the iteration are all contained in the immutable c and g.
- let à = op𝒟.findim_matrix(μ.iter_locations());
- let g̃ = DVector::from_iterator(μ.len(),
- μ.iter_locations()
- .map(|ζ| minus_τv.apply(ζ) + ω0.apply(ζ))
- .map(F::to_nalgebra_mixed));
- let mut x = μ.masses_dvector();
- // The gradient of the forward component of the inner objective is C^*𝒟Cx - g̃.
- // We have |C^*𝒟Cx|_2 = sup_{|z|_2 ≤ 1} ⟨z, C^*𝒟Cx⟩ = sup_{|z|_2 ≤ 1} ⟨Cz|𝒟Cx⟩
- // ≤ sup_{|z|_2 ≤ 1} |Cz|_ℳ |𝒟Cx|_∞ ≤ sup_{|z|_2 ≤ 1} |Cz|_ℳ |𝒟| |Cx|_ℳ
- // ≤ sup_{|z|_2 ≤ 1} |z|_1 |𝒟| |x|_1 ≤ sup_{|z|_2 ≤ 1} n |z|_2 |𝒟| |x|_2
- // = n |𝒟| |x|_2, where n is the number of points. Therefore
- let Ã_normest = op𝒟norm * F::cast_from(μ.len());
-
- // Solve finite-dimensional subproblem.
- inner_iters += reg.solve_findim(&Ã, &g̃, τ, &mut x, Ã_normest, ε, config);
-
- // Update masses of μ based on solution of finite-dimensional subproblem.
- μ.set_masses_dvector(&x);
- }
-
- // Form d = ω0 - τv - 𝒟μ = -𝒟(μ - μ^k) - τv for checking the proximate optimality
- // conditions in the predual space, and finding new points for insertion, if necessary.
- let mut d = &minus_τv + op𝒟.preapply(μ_diff(&μ, &μ_base));
+ // Insert and reweigh
+ let (maybe_d, _within_tolerances) = prox_penalty.insert_and_reweigh(
+ &mut μ, &mut τv, &μ_base, None, τ, ε, config, ®, &state, &mut stats,
+ );
- // If no merging heuristic is used, let's be more conservative about spike insertion,
- // and skip it after first round. If merging is done, being more greedy about spike
- // insertion also seems to improve performance.
- let skip_by_rough_check = if let SpikeMergingMethod::None = config.merging {
- false
- } else {
- count > 0
- };
-
- // Find a spike to insert, if needed
- let (ξ, _v_ξ, in_bounds) = match reg.find_tolerance_violation(
- &mut d, τ, ε, skip_by_rough_check, config
- ) {
- None => break 'insertion (true, d),
- Some(res) => res,
- };
-
- // Break if maximum insertion count reached
- if count >= max_insertions {
- break 'insertion (in_bounds, d)
- }
-
- // No point in optimising the weight here; the finite-dimensional algorithm is fast.
- μ += DeltaMeasure { x : ξ, α : 0.0 };
- count += 1;
- };
-
- if !within_tolerances && warn_insertions {
- // Complain (but continue) if we failed to get within tolerances
- // by inserting more points.
- let err = format!("Maximum insertions reached without achieving \
- subproblem solution tolerance");
- println!("{}", err.red());
+ // Prune and possibly merge spikes
+ if config.merge_now(&state) {
+ stats.merged += prox_penalty.merge_spikes(
+ &mut μ,
+ &mut τv,
+ &μ_base,
+ None,
+ τ,
+ ε,
+ config,
+ ®,
+ Some(|μ̃: &RNDM<F, N>| L2Squared.calculate_fit_op(μ̃, opA, b)),
+ );
}
- // Merge spikes
- if state.iteration() % config.merge_every == 0 {
- let n_before_merge = μ.len();
- μ.merge_spikes(config.merging, |μ_candidate| {
- let mut d = &minus_τv + op𝒟.preapply(μ_diff(&μ_candidate, &μ_base));
+ stats.pruned += prune_with_stats(&mut μ);
+
+ // Update residual
+ residual = calculate_residual(&μ, opA, b);
+
+ let iter = state.iteration();
+ stats.this_iters += 1;
- reg.verify_merge_candidate(&mut d, μ_candidate, τ, ε, &config)
- .then_some(())
- });
- debug_assert!(μ.len() >= n_before_merge);
- merged += μ.len() - n_before_merge;
- }
-
- let n_before_prune = μ.len();
- (residual, τ) = match specialisation.update(&mut μ, &μ_base) {
- (r, None) => (r, τ),
- (r, Some(new_τ)) => (r, new_τ)
- };
- debug_assert!(μ.len() <= n_before_prune);
- pruned += n_before_prune - μ.len();
-
- this_iters += 1;
+ // Give statistics if needed
+ state.if_verbose(|| {
+ plotter.plot_spikes(iter, maybe_d.as_ref(), Some(&τv), &μ);
+ full_stats(
+ &residual,
+ &μ,
+ ε,
+ std::mem::replace(&mut stats, IterInfo::new()),
+ )
+ });
// Update main tolerance for next iteration
- let ε_prev = ε;
- ε = tolerance.update(ε, state.iteration());
+ ε = tolerance.update(ε, iter);
+ }
- // Give function value if needed
- state.if_verbose(|| {
- let value_μ = specialisation.value_μ(&μ);
- // 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), value_μ,
- );
- // Calculate mean inner iterations and reset relevant counters.
- // Return the statistics
- let res = IterInfo {
- value : specialisation.calculate_fit(&μ, &residual) + reg.apply(value_μ),
- n_spikes : value_μ.len(),
- inner_iters,
- this_iters,
- merged,
- pruned,
- ε : ε_prev,
- postprocessing: config.postprocessing.then(|| value_μ.clone()),
- };
- inner_iters = 0;
- this_iters = 0;
- merged = 0;
- pruned = 0;
- res
- })
- });
-
- specialisation.postprocess(μ, config.final_merging)
+ postprocess(μ, config, L2Squared, opA, b)
}
-/// Iteratively solve the pointsource localisation problem using forward-backward splitting
+/// Iteratively solve the pointsource localisation problem using inertial forward-backward splitting.
///
/// The settings in `config` have their [respective documentation](FBConfig). `opA` is the
/// forward operator $A$, $b$ the observable, and $\lambda$ the regularisation weight.
@@ -965,118 +276,114 @@
///
/// For details on the mathematical formulation, see the [module level](self) documentation.
///
+/// The implementation relies on [`alg_tools::bisection_tree::BTFN`] presentations of
+/// sums of simple functions usign bisection trees, and the related
+/// [`alg_tools::bisection_tree::Aggregator`]s, to efficiently search for component functions
+/// active at a specific points, and to maximise their sums. Through the implementation of the
+/// [`alg_tools::bisection_tree::BT`] bisection trees, it also relies on the copy-on-write features
+/// of [`std::sync::Arc`] to only update relevant parts of the bisection tree when adding functions.
+///
/// Returns the final iterate.
#[replace_float_literals(F::cast_from(literal))]
-pub fn pointsource_fb_reg<'a, F, I, A, GA, 𝒟, BTA, G𝒟, S, K, Reg, const N : usize>(
- opA : &'a A,
- b : &A::Observable,
- reg : Reg,
- op𝒟 : &'a 𝒟,
- config : &FBConfig<F>,
- iterator : I,
- plotter : SeqPlotter<F, N>,
-) -> DiscreteMeasure<Loc<F, N>, F>
-where F : Float + ToNalgebraRealField,
- I : AlgIteratorFactory<IterInfo<F, N>>,
- for<'b> &'b A::Observable : std::ops::Neg<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>>
- + Lipschitz<𝒟, 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>,
- Cube<F, N>: P2Minimise<Loc<F, N>, F>,
- PlotLookup : Plotting<N>,
- DiscreteMeasure<Loc<F, N>, F> : SpikeMerging<F>,
- Reg : RegTerm<F, N> {
+pub fn pointsource_fista_reg<F, I, A, Reg, P, const N: usize>(
+ opA: &A,
+ b: &A::Observable,
+ reg: Reg,
+ prox_penalty: &P,
+ fbconfig: &FBConfig<F>,
+ iterator: I,
+ mut plotter: SeqPlotter<F, N>,
+) -> RNDM<F, N>
+where
+ F: Float + ToNalgebraRealField,
+ I: AlgIteratorFactory<IterInfo<F, N>>,
+ for<'b> &'b A::Observable: std::ops::Neg<Output = A::Observable>,
+ A: ForwardModel<RNDM<F, N>, F> + AdjointProductBoundedBy<RNDM<F, N>, P, FloatType = F>,
+ A::PreadjointCodomain: RealMapping<F, N>,
+ PlotLookup: Plotting<N>,
+ RNDM<F, N>: SpikeMerging<F>,
+ Reg: RegTerm<F, N>,
+ P: ProxPenalty<F, A::PreadjointCodomain, Reg, N>,
+{
+ // Set up parameters
+ let config = &fbconfig.generic;
+ let τ = fbconfig.τ0 / opA.adjoint_product_bound(prox_penalty).unwrap();
+ let mut λ = 1.0;
+ // 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();
- let initial_residual = -b;
- let τ = config.τ0/opA.lipschitz_factor(&op𝒟).unwrap();
+ // Initialise iterates
+ let mut μ = DiscreteMeasure::new();
+ let mut μ_prev = DiscreteMeasure::new();
+ let mut residual = -b;
+ let mut warned_merging = false;
- match config.meta {
- FBMetaAlgorithm::None => generic_pointsource_fb_reg(
- opA, reg, op𝒟, τ, &config.insertion, iterator, plotter, initial_residual,
- BasicFB{ b, opA },
- ),
- FBMetaAlgorithm::InertiaFISTA => generic_pointsource_fb_reg(
- opA, reg, op𝒟, τ, &config.insertion, iterator, plotter, initial_residual,
- FISTA{ b, opA, λ : 1.0, μ_prev : DiscreteMeasure::new() },
- ),
- }
-}
+ // Statistics
+ let full_stats = |ν: &RNDM<F, N>, ε, stats| IterInfo {
+ value: L2Squared.calculate_fit_op(ν, opA, b) + reg.apply(ν),
+ n_spikes: ν.len(),
+ ε,
+ // postprocessing: config.postprocessing.then(|| ν.clone()),
+ ..stats
+ };
+ let mut stats = IterInfo::new();
-//
-// Deprecated interfaces
-//
+ // 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(residual * τ);
-#[deprecated(note = "Use `pointsource_fb_reg`")]
-pub fn pointsource_fb<'a, F, I, A, GA, 𝒟, BTA, G𝒟, S, K, const N : usize>(
- opA : &'a A,
- b : &A::Observable,
- α : F,
- op𝒟 : &'a 𝒟,
- config : &FBConfig<F>,
- iterator : I,
- plotter : SeqPlotter<F, N>
-) -> DiscreteMeasure<Loc<F, N>, F>
-where F : Float + ToNalgebraRealField,
- I : AlgIteratorFactory<IterInfo<F, N>>,
- for<'b> &'b A::Observable : std::ops::Neg<Output=A::Observable>,
- 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>>
- + Lipschitz<𝒟, 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>,
- Cube<F, N>: P2Minimise<Loc<F, N>, F>,
- PlotLookup : Plotting<N>,
- DiscreteMeasure<Loc<F, N>, F> : SpikeMerging<F> {
+ // Save current base point
+ let μ_base = μ.clone();
+
+ // Insert new spikes and reweigh
+ let (maybe_d, _within_tolerances) = prox_penalty.insert_and_reweigh(
+ &mut μ, &mut τv, &μ_base, None, τ, ε, config, ®, &state, &mut stats,
+ );
- pointsource_fb_reg(opA, b, NonnegRadonRegTerm(α), op𝒟, config, iterator, plotter)
-}
+ // (Do not) merge spikes.
+ if config.merge_now(&state) && !warned_merging {
+ let err = format!("Merging not supported for μFISTA");
+ println!("{}", err.red());
+ warned_merging = true;
+ }
+ // Update inertial prameters
+ let λ_prev = λ;
+ λ = 2.0 * λ_prev / (λ_prev + (4.0 + λ_prev * λ_prev).sqrt());
+ let θ = λ / λ_prev - λ;
-#[deprecated(note = "Use `generic_pointsource_fb_reg`")]
-pub fn generic_pointsource_fb<'a, F, I, A, GA, 𝒟, BTA, G𝒟, S, K, Spec, const N : usize>(
- opA : &'a A,
- α : F,
- op𝒟 : &'a 𝒟,
- τ : F,
- config : &FBGenericConfig<F>,
- iterator : I,
- plotter : SeqPlotter<F, N>,
- residual : A::Observable,
- specialisation : Spec,
-) -> DiscreteMeasure<Loc<F, N>, F>
-where F : Float + ToNalgebraRealField,
- I : AlgIteratorFactory<IterInfo<F, N>>,
- Spec : FBSpecialisation<F, A::Observable, N>,
- 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>>
- + Lipschitz<𝒟, 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>,
- Cube<F, N>: P2Minimise<Loc<F, N>, F>,
- PlotLookup : Plotting<N>,
- DiscreteMeasure<Loc<F, N>, F> : SpikeMerging<F> {
+ // Perform inertial update on μ.
+ // This computes μ ← (1 + θ) * μ - θ * μ_prev, pruning spikes where both μ
+ // and μ_prev have zero weight. Since both have weights from the finite-dimensional
+ // subproblem with a proximal projection step, this is likely to happen when the
+ // spike is not needed. A copy of the pruned μ without artithmetic performed is
+ // stored in μ_prev.
+ let n_before_prune = μ.len();
+ μ.pruning_sub(1.0 + θ, θ, &mut μ_prev);
+ //let μ_new = (&μ * (1.0 + θ)).sub_matching(&(&μ_prev * θ));
+ // μ_prev = μ;
+ // μ = μ_new;
+ debug_assert!(μ.len() <= n_before_prune);
+ stats.pruned += n_before_prune - μ.len();
- generic_pointsource_fb_reg(opA, NonnegRadonRegTerm(α), op𝒟, τ, config, iterator, plotter,
- residual, specialisation)
+ // Update residual
+ residual = calculate_residual(&μ, opA, b);
+
+ let iter = state.iteration();
+ stats.this_iters += 1;
+
+ // Give statistics if needed
+ state.if_verbose(|| {
+ plotter.plot_spikes(iter, maybe_d.as_ref(), Some(&τv), &μ_prev);
+ full_stats(&μ_prev, ε, std::mem::replace(&mut stats, IterInfo::new()))
+ });
+
+ // Update main tolerance for next iteration
+ ε = tolerance.update(ε, iter);
+ }
+
+ postprocess(μ_prev, config, L2Squared, opA, b)
}