--- a/src/fb.rs Sun Dec 11 23:19:17 2022 +0200
+++ b/src/fb.rs Sun Dec 11 23:25:53 2022 +0200
@@ -6,8 +6,8 @@
* Valkonen T. - _Proximal methods for point source localisation_,
[arXiv:2212.02991](https://arxiv.org/abs/2212.02991).
-The main routine is [`pointsource_fb`]. It is based on [`generic_pointsource_fb`], which is also
-used by our [primal-dual proximal splitting][crate::pdps] implementation.
+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`].
@@ -83,17 +83,17 @@
use numeric_literals::replace_float_literals;
use serde::{Serialize, Deserialize};
use colored::Colorize;
-use nalgebra::DVector;
+use nalgebra::{DVector, DMatrix};
use alg_tools::iterate::{
AlgIteratorFactory,
AlgIteratorState,
};
use alg_tools::euclidean::Euclidean;
-use alg_tools::norms::Norm;
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,
@@ -113,7 +113,6 @@
use crate::measures::{
DiscreteMeasure,
DeltaMeasure,
- Radon
};
use crate::measures::merging::{
SpikeMergingMethod,
@@ -124,7 +123,8 @@
DiscreteMeasureOp, Lipschitz
};
use crate::subproblem::{
- quadratic_nonneg,
+ nonneg::quadratic_nonneg,
+ unconstrained::quadratic_unconstrained,
InnerSettings,
InnerMethod,
};
@@ -134,6 +134,10 @@
Plotting,
PlotLookup
};
+use crate::regularisation::{
+ NonnegRadonRegTerm,
+ RadonRegTerm,
+};
/// Method for constructing $μ$ on each iteration
#[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
@@ -156,7 +160,7 @@
InertiaFISTA,
}
-/// Settings for [`pointsource_fb`].
+/// Settings for [`pointsource_fb_reg`].
#[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
#[serde(default)]
pub struct FBConfig<F : Float> {
@@ -169,7 +173,7 @@
}
/// Settings for the solution of the stepwise optimality condition in algorithms based on
-/// [`generic_pointsource_fb`].
+/// [`generic_pointsource_fb_reg`].
#[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
#[serde(default)]
pub struct FBGenericConfig<F : Float> {
@@ -236,7 +240,7 @@
}
}
-/// Trait for specialisation of [`generic_pointsource_fb`] to basic FB, FISTA.
+/// 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
@@ -301,7 +305,7 @@
}
}
-/// Specialisation of [`generic_pointsource_fb`] to basic μFB.
+/// Specialisation of [`generic_pointsource_fb_reg`] to basic μFB.
struct BasicFB<
'a,
F : Float + ToNalgebraRealField,
@@ -340,7 +344,7 @@
}
}
-/// Specialisation of [`generic_pointsource_fb`] to FISTA.
+/// Specialisation of [`generic_pointsource_fb_reg`] to FISTA.
struct FISTA<
'a,
F : Float + ToNalgebraRealField,
@@ -423,62 +427,290 @@
}
}
-/// Iteratively solve the pointsource localisation problem using 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.
-/// 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.
-///
-/// Returns the final iterate.
+
+/// 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;
+}
+
#[replace_float_literals(F::cast_from(literal))]
-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>,
- //+ 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> {
+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;
- let initial_residual = -b;
- let τ = config.τ0/opA.lipschitz_factor(&op𝒟).unwrap();
+ // 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))
+ }
+ }
- match config.meta {
- FBMetaAlgorithm::None => generic_pointsource_fb(
- opA, α, op𝒟, τ, &config.insertion, iterator, plotter, initial_residual,
- BasicFB{ b, opA }
- ),
- FBMetaAlgorithm::InertiaFISTA => generic_pointsource_fb(
- opA, α, op𝒟, τ, &config.insertion, iterator, plotter, initial_residual,
- FISTA{ b, opA, λ : 1.0, μ_prev : DiscreteMeasure::new() }
- ),
+ 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.α()
}
}
-/// Generic implementation of [`pointsource_fb`].
+#[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`].
///
/// The method can be specialised to even primal-dual proximal splitting through the
/// [`FBSpecialisation`] parameter `specialisation`.
@@ -497,16 +729,18 @@
///
/// Returns the final iterate.
#[replace_float_literals(F::cast_from(literal))]
-pub fn generic_pointsource_fb<'a, F, I, A, GA, 𝒟, BTA, G𝒟, S, K, Spec, const N : usize>(
+pub fn generic_pointsource_fb_reg<
+ 'a, F, I, A, GA, 𝒟, BTA, G𝒟, S, K, Spec, Reg, const N : usize
+>(
opA : &'a A,
- α : F,
+ reg : Reg,
op𝒟 : &'a 𝒟,
mut τ : F,
config : &FBGenericConfig<F>,
iterator : I,
mut plotter : SeqPlotter<F, N>,
mut residual : A::Observable,
- mut specialisation : Spec,
+ mut specialisation : Spec
) -> DiscreteMeasure<Loc<F, N>, F>
where F : Float + ToNalgebraRealField,
I : AlgIteratorFactory<IterInfo<F, N>>,
@@ -522,17 +756,16 @@
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> {
+ DiscreteMeasure<Loc<F, N>, F> : SpikeMerging<F>,
+ Reg : RegTerm<F, N> {
// Set up parameters
let quiet = iterator.is_quiet();
let op𝒟norm = op𝒟.opnorm_bound();
- // 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 * τ * α;
+ // 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
@@ -568,17 +801,6 @@
// Run the algorithm
iterator.iterate(|state| {
- // Calculate subproblem tolerances, and update main tolerance for next iteration
- let τα = τ * α;
- let target_bounds = Bounds(τα - ε, τα + ε);
- let merge_tolerance = config.merge_tolerance_mult * ε;
- let merge_target_bounds = Bounds(τα - merge_tolerance, τα + merge_tolerance);
- let inner_tolerance = ε * config.inner.tolerance_mult;
- let refinement_tolerance = ε * config.refinement.tolerance_mult;
- let maximise_above = τα + ε * config.insertion_cutoff_factor;
- let ε_prev = ε;
- ε = tolerance.update(ε, state.iteration());
-
// 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),
@@ -626,12 +848,10 @@
// ≤ 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 inner_τ = config.inner.τ0 / (op𝒟norm * F::cast_from(μ.len()));
+ let Ã_normest = op𝒟norm * F::cast_from(μ.len());
// Solve finite-dimensional subproblem.
- let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
- inner_iters += quadratic_nonneg(config.inner.method, &Ã, &g̃, τ*α, &mut x,
- inner_τ, inner_it);
+ inner_iters += reg.solve_findim(&Ã, &g̃, τ, &mut x, Ã_normest, ε, config);
// Update masses of μ based on solution of finite-dimensional subproblem.
μ.set_masses_dvector(&x);
@@ -644,36 +864,23 @@
// 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 may_break = if let SpikeMergingMethod::None = config.merging {
+ let skip_by_rough_check = if let SpikeMergingMethod::None = config.merging {
false
} else {
count > 0
};
- // If preliminary check indicates that we are in bonds, and if it otherwise matches
- // the insertion strategy, skip insertion.
- if may_break && target_bounds.superset(&d.bounds()) {
- break 'insertion (true, d)
- }
-
- // If the rough check didn't indicate stopping, find maximising point, maintaining for
- // the calculations in the beginning of the loop that v_ξ = (ω0-τv-𝒟μ)(ξ) = d(ξ),
- // where 𝒟μ is now distinct from μ0 after the insertions already performed.
- // We do not need to check lower bounds, as a solution of the finite-dimensional
- // subproblem should always satisfy them.
-
- // If μ has some spikes, only find a maximum of d if it is above a threshold
- // defined by the refinment tolerance.
- let (ξ, v_ξ) = match d.maximise_above(maximise_above, refinement_tolerance,
- config.refinement.max_steps) {
+ // 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 {
- let in_bounds2 = target_bounds.upper() >= v_ξ;
- break 'insertion (in_bounds2, d)
+ break 'insertion (in_bounds, d)
}
// No point in optimising the weight here; the finite-dimensional algorithm is fast.
@@ -695,21 +902,8 @@
μ.merge_spikes(config.merging, |μ_candidate| {
let mut d = &minus_τv + op𝒟.preapply(μ_diff(&μ_candidate, &μ_base));
- if merge_target_bounds.superset(&d.bounds()) {
- return Some(())
- }
-
- let d_min_supp = μ_candidate.iter_spikes().filter_map(|&DeltaMeasure{ α, ref x }| {
- (α != 0.0).then(|| d.apply(x))
- }).reduce(F::min);
-
- if d_min_supp.map_or(true, |b| b >= merge_target_bounds.lower()) &&
- d.has_upper_bound(merge_target_bounds.upper(), refinement_tolerance,
- config.refinement.max_steps) {
- Some(())
- } else {
- None
- }
+ reg.verify_merge_candidate(&mut d, μ_candidate, τ, ε, &config)
+ .then_some(())
});
debug_assert!(μ.len() >= n_before_merge);
merged += μ.len() - n_before_merge;
@@ -725,6 +919,10 @@
this_iters += 1;
+ // Update main tolerance for next iteration
+ let ε_prev = ε;
+ ε = tolerance.update(ε, state.iteration());
+
// Give function value if needed
state.if_verbose(|| {
let value_μ = specialisation.value_μ(&μ);
@@ -732,12 +930,12 @@
plotter.plot_spikes(
format!("iter {} end; {}", state.iteration(), within_tolerances), &d,
"start".to_string(), Some(&minus_τv),
- Some(target_bounds), value_μ,
+ reg.target_bounds(τ, ε_prev), value_μ,
);
// Calculate mean inner iterations and reset relevant counters.
// Return the statistics
let res = IterInfo {
- value : specialisation.calculate_fit(&μ, &residual) + α * value_μ.norm(Radon),
+ value : specialisation.calculate_fit(&μ, &residual) + reg.apply(value_μ),
n_spikes : value_μ.len(),
inner_iters,
this_iters,
@@ -757,6 +955,128 @@
specialisation.postprocess(μ, config.final_merging)
}
+/// Iteratively solve the pointsource localisation problem using 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.
+/// 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.
+///
+/// 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> {
+
+ let initial_residual = -b;
+ let τ = config.τ0/opA.lipschitz_factor(&op𝒟).unwrap();
+
+ 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() },
+ ),
+ }
+}
+
+//
+// Deprecated interfaces
+//
+
+#[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> {
+
+ pointsource_fb_reg(opA, b, NonnegRadonRegTerm(α), op𝒟, config, iterator, plotter)
+}
+#[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> {
+ generic_pointsource_fb_reg(opA, NonnegRadonRegTerm(α), op𝒟, τ, config, iterator, plotter,
+ residual, specialisation)
+}