--- a/src/regularisation.rs Tue Aug 01 10:32:12 2023 +0300
+++ b/src/regularisation.rs Thu Aug 29 00:00:00 2024 -0500
@@ -34,9 +34,13 @@
use crate::subproblem::{
nonneg::quadratic_nonneg,
unconstrained::quadratic_unconstrained,
+ l1squared_unconstrained::l1squared_unconstrained,
+ l1squared_nonneg::l1squared_nonneg
};
use alg_tools::iterate::AlgIteratorFactory;
+use std::cmp::Ordering::{Greater, Less, Equal};
+
/// The regularisation term $α\\|μ\\|\_{ℳ(Ω)} + δ_{≥ 0}(μ)$ for [`pointsource_fb_reg`] and other
/// algorithms.
///
@@ -131,6 +135,23 @@
config : &FBGenericConfig<F>
) -> usize;
+ /// Approximately solve the problem
+ /// <div>$$
+ /// \min_{x ∈ ℝ^n} \frac{1}{2} |x-y|_1^2 - g^⊤ x + τ G(x)
+ /// $$</div>
+ /// for $G$ depending on the trait implementation.
+ ///
+ /// Returns the number of iterations taken.
+ fn solve_findim_l1squared(
+ &self,
+ y : &DVector<F::MixedType>,
+ g : &DVector<F::MixedType>,
+ τ : F,
+ x : &mut DVector<F::MixedType>,
+ ε : 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
@@ -151,8 +172,37 @@
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> {
+ self.find_tolerance_violation_slack(d, τ, ε, skip_by_rough_check, config, F::ZERO)
+ }
+
+ /// Find a point where `d` may violate the tolerance `ε`.
+ ///
+ /// This version includes a `slack` parameter to expand the tolerances.
+ /// It is used for Radon-norm squared proximal term in [`crate::radon_fb`].
+ ///
+ /// 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_slack<G, BT>(
+ &self,
+ d : &mut BTFN<F, G, BT, N>,
+ τ : F,
+ ε : F,
+ skip_by_rough_check : bool,
+ config : &FBGenericConfig<F>,
+ slack : 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.
@@ -169,6 +219,28 @@
G::SupportType : Mapping<Loc<F, N>,Codomain=F>
+ LocalAnalysis<F, Bounds<F>, N>;
+ /// Verify that `d` is in bounds `ε` for a merge candidate `μ`
+ ///
+ /// This version is s used for Radon-norm squared proximal term in [`crate::radon_fb`].
+ /// The [`DiscreteMeasure`]s `μ` and `radon_μ` are supposed to have same coordinates at
+ /// same agreeing indices.
+ ///
+ /// `ε` is the current main tolerance and `τ` a scaling factor for the regulariser.
+ fn verify_merge_candidate_radonsq<G, BT>(
+ &self,
+ d : &mut BTFN<F, G, BT, N>,
+ μ : &DiscreteMeasure<Loc<F, N>, F>,
+ τ : F,
+ ε : F,
+ config : &FBGenericConfig<F>,
+ radon_μ :&DiscreteMeasure<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>;
+
+
/// TODO: document this
fn target_bounds(&self, τ : F, ε : F) -> Option<Bounds<F>>;
@@ -197,6 +269,9 @@
G : SupportGenerator<F, N, Id=BT::Data>,
G::SupportType : Mapping<Loc<F, N>,Codomain=F>
+ LocalAnalysis<F, Bounds<F>, N>;
+
+ /// Convert bound on the regulariser to a bond on the Radon norm
+ fn radon_norm_bound(&self, b : F) -> F;
}
#[replace_float_literals(F::cast_from(literal))]
@@ -215,30 +290,47 @@
) -> 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)
+ quadratic_nonneg(mA, g, τ * self.α(), x,
+ mA_normest, &config.inner, inner_it)
+
}
+ fn solve_findim_l1squared(
+ &self,
+ y : &DVector<F::MixedType>,
+ g : &DVector<F::MixedType>,
+ τ : F,
+ x : &mut DVector<F::MixedType>,
+ ε : F,
+ config : &FBGenericConfig<F>
+ ) -> usize {
+ let inner_tolerance = ε * config.inner.tolerance_mult;
+ let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
+ l1squared_nonneg(y, g, τ * self.α(), 1.0, x,
+ &config.inner, inner_it)
+ }
+
+
#[inline]
- fn find_tolerance_violation<G, BT>(
+ fn find_tolerance_violation_slack<G, BT>(
&self,
d : &mut BTFN<F, G, BT, N>,
τ : F,
ε : F,
skip_by_rough_check : bool,
config : &FBGenericConfig<F>,
+ slack : 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 keep_below = τα + slack + ε;
+ let maximise_above = τα + slack + ε * 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
+ // If preliminary check indicates that we are in bounds, and if it otherwise matches
// the insertion strategy, skip insertion.
if skip_by_rough_check && d.bounds().upper() <= keep_below {
None
@@ -271,9 +363,9 @@
return (
bnd.lower() >= keep_supp_above
||
- μ.iter_spikes().map(|&DeltaMeasure{ α : β, ref x }| {
- (β == 0.0) || d.apply(x) >= keep_supp_above
- }).all(std::convert::identity)
+ μ.iter_spikes().all(|&DeltaMeasure{ α, ref x }| {
+ (α == 0.0) || d.apply(x) >= keep_supp_above
+ })
) && (
bnd.upper() <= keep_below
||
@@ -281,6 +373,50 @@
)
}
+ fn verify_merge_candidate_radonsq<G, BT>(
+ &self,
+ d : &mut BTFN<F, G, BT, N>,
+ μ : &DiscreteMeasure<Loc<F, N>, F>,
+ τ : F,
+ ε : F,
+ config : &FBGenericConfig<F>,
+ radon_μ :&DiscreteMeasure<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 refinement_tolerance = ε * config.refinement.tolerance_mult;
+ let merge_tolerance = config.merge_tolerance_mult * ε;
+ let slack = radon_μ.norm(Radon);
+ let bnd = d.bounds();
+
+ return {
+ μ.both_matching(radon_μ)
+ .all(|(α, rα, x)| {
+ let v = d.apply(x);
+ let (l1, u1) = match α.partial_cmp(&0.0).unwrap_or(Equal) {
+ Greater => (τα, τα),
+ _ => (F::NEG_INFINITY, τα),
+ // Less should not happen; treated as Equal
+ };
+ let (l2, u2) = match rα.partial_cmp(&0.0).unwrap_or(Equal) {
+ Greater => (slack, slack),
+ Equal => (-slack, slack),
+ Less => (-slack, -slack),
+ };
+ // TODO: both fail.
+ (l1 + l2 - merge_tolerance <= v) && (v <= u1 + u2 + merge_tolerance)
+ })
+ } && {
+ let keep_below = τα + slack + merge_tolerance;
+ 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(τα - ε, τα + ε))
@@ -310,9 +446,12 @@
G : SupportGenerator<F, N, Id=BT::Data>,
G::SupportType : Mapping<Loc<F, N>,Codomain=F>
+ LocalAnalysis<F, Bounds<F>, N> {
- //let w = |x| 1.0.min((ε + d.apply(x))/(τ * self.α()));
- let τw = |x| τ.min((ε + d.apply(x))/self.α());
- τw(z) - τw(y)
+ let w = |x| 1.0.min((ε + d.apply(x))/(τ * self.α()));
+ w(z) - w(y)
+ }
+
+ fn radon_norm_bound(&self, b : F) -> F {
+ b / self.α()
}
}
@@ -331,28 +470,43 @@
) -> 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)
+ quadratic_unconstrained(mA, g, τ * self.α(), x,
+ mA_normest, &config.inner, inner_it)
}
- fn find_tolerance_violation<G, BT>(
+ fn solve_findim_l1squared(
+ &self,
+ y : &DVector<F::MixedType>,
+ g : &DVector<F::MixedType>,
+ τ : F,
+ x : &mut DVector<F::MixedType>,
+ ε : F,
+ config : &FBGenericConfig<F>
+ ) -> usize {
+ let inner_tolerance = ε * config.inner.tolerance_mult;
+ let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
+ l1squared_unconstrained(y, g, τ * self.α(), 1.0, x,
+ &config.inner, inner_it)
+ }
+
+ fn find_tolerance_violation_slack<G, BT>(
&self,
d : &mut BTFN<F, G, BT, N>,
τ : F,
ε : F,
skip_by_rough_check : bool,
config : &FBGenericConfig<F>,
+ slack : 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 keep_below = τα + slack + ε;
+ let keep_above = -(τα + slack) - ε;
+ let maximise_above = τα + slack + ε * config.insertion_cutoff_factor;
+ let minimise_below = -(τα + slack) - ε * 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
@@ -405,14 +559,13 @@
return (
(bnd.lower() >= keep_supp_pos_above && bnd.upper() <= keep_supp_neg_below)
||
- μ.iter_spikes().map(|&DeltaMeasure{ α : β, ref x }| {
- use std::cmp::Ordering::*;
+ μ.iter_spikes().all(|&DeltaMeasure{ α : β, ref x }| {
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
||
@@ -426,6 +579,56 @@
)
}
+ fn verify_merge_candidate_radonsq<G, BT>(
+ &self,
+ d : &mut BTFN<F, G, BT, N>,
+ μ : &DiscreteMeasure<Loc<F, N>, F>,
+ τ : F,
+ ε : F,
+ config : &FBGenericConfig<F>,
+ radon_μ : &DiscreteMeasure<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 refinement_tolerance = ε * config.refinement.tolerance_mult;
+ let merge_tolerance = config.merge_tolerance_mult * ε;
+ let slack = radon_μ.norm(Radon);
+ let bnd = d.bounds();
+
+ return {
+ μ.both_matching(radon_μ)
+ .all(|(α, rα, x)| {
+ let v = d.apply(x);
+ let (l1, u1) = match α.partial_cmp(&0.0).unwrap_or(Equal) {
+ Greater => (τα, τα),
+ Equal => (-τα, τα),
+ Less => (-τα, -τα),
+ };
+ let (l2, u2) = match rα.partial_cmp(&0.0).unwrap_or(Equal) {
+ Greater => (slack, slack),
+ Equal => (-slack, slack),
+ Less => (-slack, -slack),
+ };
+ (l1 + l2 - merge_tolerance <= v) && (v <= u1 + u2 + merge_tolerance)
+ })
+ } && {
+ let keep_below = τα + slack + merge_tolerance;
+ bnd.upper() <= keep_below
+ ||
+ d.has_upper_bound(keep_below, refinement_tolerance,
+ config.refinement.max_steps)
+ } && {
+ let keep_above = -τα - slack - merge_tolerance;
+ 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(-τα - ε, τα + ε))
@@ -457,14 +660,14 @@
+ LocalAnalysis<F, Bounds<F>, N> {
let α = self.α();
- // let w = |x| {
- // let dx = d.apply(x);
- // ((-ε + dx)/(τ * α)).max(1.0.min(ε + dx)/(τ * α))
- // };
- let τw = |x| {
+ let w = |x| {
let dx = d.apply(x);
- ((-ε + dx)/α).max(τ.min(ε + dx)/α)
+ ((-ε + dx)/(τ * α)).max(1.0.min(ε + dx)/(τ * α))
};
- τw(z) - τw(y)
+ w(z) - w(y)
+ }
+
+ fn radon_norm_bound(&self, b : F) -> F {
+ b / self.α()
}
}
\ No newline at end of file