pointsource_algs: comparison src/frank

-:d29d1fcf5423
+:79943be70720
 use alg_tools::iterate::{
 AlgIteratorFactory,
 AlgIteratorState,
 AlgIteratorOptions,
+ValueIteratorFactory,
 };
 use alg_tools::euclidean::Euclidean;
 use alg_tools::norms::Norm;
 use alg_tools::linops::Apply;
 use alg_tools::sets::Cube;
 };
 use crate::forward_model::ForwardModel;
 #[allow(unused_imports)] // Used in documentation
 use crate::subproblem::{
 unconstrained::quadratic_unconstrained,
+nonneg::quadratic_nonneg,
 InnerSettings,
 InnerMethod,
 };
 use crate::tolerance::Tolerance;
 use crate::plot::{
 SeqPlotter,
 Plotting,
 PlotLookup
 };
+use crate::regularisation::{
+NonnegRadonRegTerm,
+RadonRegTerm,
+};
+use crate::fb::RegTerm;
 /// Settings for [`pointsource_fw`].
 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
 #[serde(default)]
 pub struct FWConfig<F : Float> {
 /// Helper struct for pre-initialising the finite-dimensional subproblems solver
 /// [`prepare_optimise_weights`].
 ///
 /// The pre-initialisation is done by [`prepare_optimise_weights`].
 pub struct FindimData<F : Float> {
-opAnorm_squared : F
+/// ‖A‖^2
-}
+opAnorm_squared : F,
+/// Bound $M_0$ from the Bredies–Pikkarainen article.
-/// Return a pre-initialisation struct for [`prepare_optimise_weights`].
+m0 : F
-///
+}
-/// The parameter `opA` is the forward operator $A$.
-pub fn prepare_optimise_weights<F, A, const N : usize>(opA : &A) -> FindimData<F>
+/// Trait for finite dimensional weight optimisation.
-where F : Float + ToNalgebraRealField,
+pub trait WeightOptim<
+F : Float + ToNalgebraRealField,
+A : ForwardModel<Loc<F, N>, F>,
+I : AlgIteratorFactory<F>,
+const N : usize
+> {
+/// Return a pre-initialisation struct for [`Self::optimise_weights`].
+///
+/// The parameter `opA` is the forward operator $A$.
+fn prepare_optimise_weights(&self, opA : &A, b : &A::Observable) -> FindimData<F>;
+/// Solve the finite-dimensional weight optimisation problem for the 2-norm-squared data fidelity
+/// point source localisation problem.
+///
+/// That is, we minimise
+/// <div>$$
+///     μ ↦ \frac{1}{2}\|Aμ-b\|_w^2 + G(μ)
+/// $$</div>
+/// only with respect to the weights of $μ$. Here $G$ is a regulariser modelled by `Self`.
+///
+/// The parameter `μ` is the discrete measure whose weights are to be optimised.
+/// The `opA` parameter is the forward operator $A$, while `b`$ and `α` are as in the
+/// objective above. The method parameter are set in `inner` (see [`InnerSettings`]), while
+/// `iterator` is used to iterate the steps of the method, and `plotter` may be used to
+/// save intermediate iteration states as images. The parameter `findim_data` should be
+/// prepared using [`Self::prepare_optimise_weights`]:
+///
+/// Returns the number of iterations taken by the method configured in `inner`.
+fn optimise_weights<'a>(
+&self,
+μ : &mut DiscreteMeasure<Loc<F, N>, F>,
+opA : &'a A,
+b : &A::Observable,
+findim_data : &FindimData<F>,
+inner : &InnerSettings<F>,
+iterator : I
+) -> usize;
+}
+/// Trait for regularisation terms supported by [`pointsource_fw_reg`].
+pub trait RegTermFW<
+F : Float + ToNalgebraRealField,
+A : ForwardModel<Loc<F, N>, F>,
+I : AlgIteratorFactory<F>,
+const N : usize
+> : RegTerm<F, N>
++ WeightOptim<F, A, I, N>
++ for<'a> Apply<&'a DiscreteMeasure<Loc<F, N>, F>, Output = F> {
+/// With $g = A\_\*(Aμ-b)$, returns $(x, g(x))$ for $x$ a new point to be inserted
+/// into $μ$, as determined by the regulariser.
+///
+/// The parameters `refinement_tolerance` and `max_steps` are passed to relevant
+/// [`BTFN`] minimisation and maximisation routines.
+fn find_insertion(
+&self,
+g : &mut A::PreadjointCodomain,
+refinement_tolerance : F,
+max_steps : usize
+) -> (Loc<F, N>, F);
+/// Insert point `ξ` into `μ` for the relaxed algorithm from Bredies–Pikkarainen.
+fn relaxed_insert<'a>(
+&self,
+μ : &mut DiscreteMeasure<Loc<F, N>, F>,
+g : &A::PreadjointCodomain,
+opA : &'a A,
+ξ : Loc<F, N>,
+v_ξ : F,
+findim_data : &FindimData<F>
+);
+}
+#[replace_float_literals(F::cast_from(literal))]
+impl<F : Float + ToNalgebraRealField, A, I, const N : usize> WeightOptim<F, A, I, N>
+for RadonRegTerm<F>
+where I : AlgIteratorFactory<F>,
 A : ForwardModel<Loc<F, N>, F> {
-FindimData{
-opAnorm_squared : opA.opnorm_bound().powi(2)
+fn prepare_optimise_weights(&self, opA : &A, b : &A::Observable) -> FindimData<F> {
-}
+FindimData{
-}
+opAnorm_squared : opA.opnorm_bound().powi(2),
+m0 : b.norm2_squared() / (2.0 * self.α()),
-/// Solve the finite-dimensional weight optimisation problem for the 2-norm-squared data fidelity
+}
-/// point source localisation problem.
+}
-///
-/// That is, we minimise
+fn optimise_weights<'a>(
-/// <div>$$
+&self,
-///     μ ↦ \frac{1}{2}\|Aμ-b\|_w^2 + α\|μ\|_ℳ + δ_{≥ 0}(μ)
+μ : &mut DiscreteMeasure<Loc<F, N>, F>,
-/// $$</div>
+opA : &'a A,
-/// only with respect to the weights of $μ$.
+b : &A::Observable,
-///
+findim_data : &FindimData<F>,
-/// The parameter `μ` is the discrete measure whose weights are to be optimised.
+inner : &InnerSettings<F>,
-/// The `opA` parameter is the forward operator $A$, while `b`$ and `α` are as in the
+iterator : I
-/// objective above. The method parameter are set in `inner` (see [`InnerSettings`]), while
+) -> usize {
-/// `iterator` is used to iterate the steps of the method, and `plotter` may be used to
-/// save intermediate iteration states as images. The parameter `findim_data` should be
+// Form and solve finite-dimensional subproblem.
-/// prepared using [`prepare_optimise_weights`]:
+let (Ã, g̃) = opA.findim_quadratic_model(&μ, b);
-///
+let mut x = μ.masses_dvector();
-/// Returns the number of iterations taken by the method configured in `inner`.
-pub fn optimise_weights<'a, F, A, I, const N : usize>(
+// `inner_τ1` is based on an estimate of the operator norm of $A$ from ℳ(Ω) to
-μ : &mut DiscreteMeasure<Loc<F, N>, F>,
+// ℝ^n. This estimate is a good one for the matrix norm from ℝ^m to ℝ^n when the
-opA : &'a A,
+// former is equipped with the 1-norm. We need the 2-norm. To pass from 1-norm to
-b : &A::Observable,
+// 2-norm, we estimate
-α : F,
+//      ‖A‖_{2,2} := sup_{‖x‖_2 ≤ 1} ‖Ax‖_2 ≤ sup_{‖x‖_1 ≤ C} ‖Ax‖_2
-findim_data : &FindimData<F>,
+//                 = C sup_{‖x‖_1 ≤ 1} ‖Ax‖_2 = C ‖A‖_{1,2},
-inner : &InnerSettings<F>,
+// where C = √m satisfies ‖x‖_1 ≤ C ‖x‖_2. Since we are intested in ‖A_*A‖, no
-iterator : I
+// square root is needed when we scale:
-) -> usize
+let inner_τ = inner.τ0 / (findim_data.opAnorm_squared * F::cast_from(μ.len()));
-where F : Float + ToNalgebraRealField,
+let iters = quadratic_unconstrained(inner.method, &Ã, &g̃, self.α(),
+&mut x, inner_τ, iterator);
+// Update masses of μ based on solution of finite-dimensional subproblem.
+μ.set_masses_dvector(&x);
+iters
+}
+}
+#[replace_float_literals(F::cast_from(literal))]
+impl<F : Float + ToNalgebraRealField, A, I, S, GA, BTA, const N : usize> RegTermFW<F, A, I, N>
+for RadonRegTerm<F>
+where Cube<F, N> : P2Minimise<Loc<F, N>, F>,
 I : AlgIteratorFactory<F>,
-A : ForwardModel<Loc<F, N>, F>
+S: RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>,
-{
+GA : SupportGenerator<F, N, SupportType = S, Id = usize> + Clone,
-// Form and solve finite-dimensional subproblem.
+A : ForwardModel<Loc<F, N>, F, PreadjointCodomain = BTFN<F, GA, BTA, N>>,
-let (Ã, g̃) = opA.findim_quadratic_model(&μ, b);
+BTA : BTSearch<F, N, Data=usize, Agg=Bounds<F>> {
-let mut x = μ.masses_dvector();
+fn find_insertion(
-// `inner_τ1` is based on an estimate of the operator norm of $A$ from ℳ(Ω) to
+&self,
-// ℝ^n. This estimate is a good one for the matrix norm from ℝ^m to ℝ^n when the
+g : &mut A::PreadjointCodomain,
-// former is equipped with the 1-norm. We need the 2-norm. To pass from 1-norm to
+refinement_tolerance : F,
-// 2-norm, we estimate
+max_steps : usize
-//      ‖A‖_{2,2} := sup_{‖x‖_2 ≤ 1} ‖Ax‖_2 ≤ sup_{‖x‖_1 ≤ C} ‖Ax‖_2
+) -> (Loc<F, N>, F) {
-//                 = C sup_{‖x‖_1 ≤ 1} ‖Ax‖_2 = C ‖A‖_{1,2},
+let (ξmax, v_ξmax) = g.maximise(refinement_tolerance, max_steps);
-// where C = √m satisfies ‖x‖_1 ≤ C ‖x‖_2. Since we are intested in ‖A_*A‖, no
+let (ξmin, v_ξmin) = g.minimise(refinement_tolerance, max_steps);
-// square root is needed when we scale:
+if v_ξmin < 0.0 && -v_ξmin > v_ξmax {
-let inner_τ = inner.τ0 / (findim_data.opAnorm_squared * F::cast_from(μ.len()));
+(ξmin, v_ξmin)
-let iters = quadratic_unconstrained(inner.method, &Ã, &g̃, α, &mut x, inner_τ, iterator);
+} else {
-// Update masses of μ based on solution of finite-dimensional subproblem.
+(ξmax, v_ξmax)
-μ.set_masses_dvector(&x);
+}
+}
-iters
-}
+fn relaxed_insert<'a>(
+&self,
+μ : &mut DiscreteMeasure<Loc<F, N>, F>,
+g : &A::PreadjointCodomain,
+opA : &'a A,
+ξ : Loc<F, N>,
+v_ξ : F,
+findim_data : &FindimData<F>
+) {
+let α = self.0;
+let m0 = findim_data.m0;
+let φ = |t| if t <= m0 { α * t } else { α / (2.0 * m0) * (t*t + m0 * m0) };
+let v = if v_ξ.abs() <= α { 0.0 } else { m0 / α * v_ξ };
+let δ = DeltaMeasure { x : ξ, α : v };
+let dp = μ.apply(g) - δ.apply(g);
+let d = opA.apply(&*μ) - opA.apply(&δ);
+let r = d.norm2_squared();
+let s = if r == 0.0 {
+1.0
+} else {
+1.0.min( (α * μ.norm(Radon) - φ(v.abs()) - dp) / r)
+};
+*μ *= 1.0 - s;
+*μ += δ * s;
+}
+}
+#[replace_float_literals(F::cast_from(literal))]
+impl<F : Float + ToNalgebraRealField, A, I, const N : usize> WeightOptim<F, A, I, N>
+for NonnegRadonRegTerm<F>
+where I : AlgIteratorFactory<F>,
+A : ForwardModel<Loc<F, N>, F> {
+fn prepare_optimise_weights(&self, opA : &A, b : &A::Observable) -> FindimData<F> {
+FindimData{
+opAnorm_squared : opA.opnorm_bound().powi(2),
+m0 : b.norm2_squared() / (2.0 * self.α()),
+}
+}
+fn optimise_weights<'a>(
+&self,
+μ : &mut DiscreteMeasure<Loc<F, N>, F>,
+opA : &'a A,
+b : &A::Observable,
+findim_data : &FindimData<F>,
+inner : &InnerSettings<F>,
+iterator : I
+) -> usize {
+// Form and solve finite-dimensional subproblem.
+let (Ã, g̃) = opA.findim_quadratic_model(&μ, b);
+let mut x = μ.masses_dvector();
+// `inner_τ1` is based on an estimate of the operator norm of $A$ from ℳ(Ω) to
+// ℝ^n. This estimate is a good one for the matrix norm from ℝ^m to ℝ^n when the
+// former is equipped with the 1-norm. We need the 2-norm. To pass from 1-norm to
+// 2-norm, we estimate
+//      ‖A‖_{2,2} := sup_{‖x‖_2 ≤ 1} ‖Ax‖_2 ≤ sup_{‖x‖_1 ≤ C} ‖Ax‖_2
+//                 = C sup_{‖x‖_1 ≤ 1} ‖Ax‖_2 = C ‖A‖_{1,2},
+// where C = √m satisfies ‖x‖_1 ≤ C ‖x‖_2. Since we are intested in ‖A_*A‖, no
+// square root is needed when we scale:
+let inner_τ = inner.τ0 / (findim_data.opAnorm_squared * F::cast_from(μ.len()));
+let iters = quadratic_nonneg(inner.method, &Ã, &g̃, self.α(),
+&mut x, inner_τ, iterator);
+// Update masses of μ based on solution of finite-dimensional subproblem.
+μ.set_masses_dvector(&x);
+iters
+}
+}
+#[replace_float_literals(F::cast_from(literal))]
+impl<F : Float + ToNalgebraRealField, A, I, S, GA, BTA, const N : usize> RegTermFW<F, A, I, N>
+for NonnegRadonRegTerm<F>
+where Cube<F, N> : P2Minimise<Loc<F, N>, F>,
+I : AlgIteratorFactory<F>,
+S: RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>,
+GA : SupportGenerator<F, N, SupportType = S, Id = usize> + Clone,
+A : ForwardModel<Loc<F, N>, F, PreadjointCodomain = BTFN<F, GA, BTA, N>>,
+BTA : BTSearch<F, N, Data=usize, Agg=Bounds<F>> {
+fn find_insertion(
+&self,
+g : &mut A::PreadjointCodomain,
+refinement_tolerance : F,
+max_steps : usize
+) -> (Loc<F, N>, F) {
+g.maximise(refinement_tolerance, max_steps)
+}
+fn relaxed_insert<'a>(
+&self,
+μ : &mut DiscreteMeasure<Loc<F, N>, F>,
+g : &A::PreadjointCodomain,
+opA : &'a A,
+ξ : Loc<F, N>,
+v_ξ : F,
+findim_data : &FindimData<F>
+) {
+// This is just a verbatim copy of RadonRegTerm::relaxed_insert.
+let α = self.0;
+let m0 = findim_data.m0;
+let φ = |t| if t <= m0 { α * t } else { α / (2.0 * m0) * (t*t + m0 * m0) };
+let v = if v_ξ.abs() <= α { 0.0 } else { m0 / α * v_ξ };
+let δ = DeltaMeasure { x : ξ, α : v };
+let dp = μ.apply(g) - δ.apply(g);
+let d = opA.apply(&*μ) - opA.apply(&δ);
+let r = d.norm2_squared();
+let s = if r == 0.0 {
+1.0
+} else {
+1.0.min( (α * μ.norm(Radon) - φ(v.abs()) - dp) / r)
+};
+*μ *= 1.0 - s;
+*μ += δ * s;
+}
+}
 /// Solve point source localisation problem using a conditional gradient method
 /// for the 2-norm-squared data fidelity, i.e., the problem
 /// <div>$$
-///     \min_μ \frac{1}{2}\|Aμ-b\|_w^2 + α\|μ\|_ℳ + δ_{≥ 0}(μ).
+///     \min_μ \frac{1}{2}\|Aμ-b\|_w^2 + G(μ),
-/// $$</div>
+/// $$
+/// where $G$ is the regularisation term modelled by `reg`.
+/// </div>
 ///
-/// The `opA` parameter is the forward operator $A$, while `b`$ and `α` are as in the
+/// The `opA` parameter is the forward operator $A$, while `b`$ is as in the
 /// objective above. The method parameter are set in `config` (see [`FWConfig`]), while
 /// `iterator` is used to iterate the steps of the method, and `plotter` may be used to
 /// save intermediate iteration states as images.
 #[replace_float_literals(F::cast_from(literal))]
-pub fn pointsource_fw<'a, F, I, A, GA, BTA, S, const N : usize>(
+pub fn pointsource_fw_reg<'a, F, I, A, GA, BTA, S, Reg, const N : usize>(
 opA : &'a A,
 b : &A::Observable,
-α : F,
+reg : Reg,
 //domain : Cube<F, N>,
 config : &FWConfig<F>,
 iterator : I,
 mut plotter : SeqPlotter<F, N>,
 ) -> DiscreteMeasure<Loc<F, N>, F>
 BTA : BTSearch<F, N, Data=usize, Agg=Bounds<F>>,
 S: 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 : RegTermFW<F, A, ValueIteratorFactory<F, AlgIteratorOptions>, N> {
 // Set up parameters
 // We multiply tolerance by α for all algoritms.
-let tolerance = config.tolerance * α;
+let tolerance = config.tolerance * reg.tolerance_scaling();
 let mut ε = tolerance.initial();
-let findim_data = prepare_optimise_weights(opA);
+let findim_data = reg.prepare_optimise_weights(opA, b);
-let m0 = b.norm2_squared() / (2.0 * α);
-let φ = |t| if t <= m0 { α * t } else { α / (2.0 * m0) * (t*t + m0 * m0) };
 // Initialise operators
 let preadjA = opA.preadjoint();
 // Initialise iterates
 // the residual and replacing it below before the end of this closure.
 let r = std::mem::replace(&mut residual, opA.empty_observable());
 let mut g = -preadjA.apply(r);
 // Find absolute value maximising point
-let (ξmax, v_ξmax) = g.maximise(refinement_tolerance,
+let (ξ, v_ξ) = reg.find_insertion(&mut g, refinement_tolerance,
 config.refinement.max_steps);
-let (ξmin, v_ξmin) = g.minimise(refinement_tolerance,
-config.refinement.max_steps);
-let (ξ, v_ξ) = if v_ξmin < 0.0 && -v_ξmin > v_ξmax {
-(ξmin, v_ξmin)
-} else {
-(ξmax, v_ξmax)
-};
 let inner_it = match config.variant {
 FWVariant::FullyCorrective => {
 // No point in optimising the weight here: the finite-dimensional algorithm is fast.
 μ += DeltaMeasure { x : ξ, α : 0.0 };
 config.inner.iterator_options.stop_target(inner_tolerance)
 },
 FWVariant::Relaxed => {
 // Perform a relaxed initialisation of μ
-let v = if v_ξ.abs() <= α { 0.0 } else { m0 / α * v_ξ };
+reg.relaxed_insert(&mut μ, &g, opA, ξ, v_ξ, &findim_data);
-let δ = DeltaMeasure { x : ξ, α : v };
-let dp = μ.apply(&g) - δ.apply(&g);
-let d = opA.apply(&μ) - opA.apply(&δ);
-let r = d.norm2_squared();
-let s = if r == 0.0 {
-1.0
-} else {
-1.0.min( (α * μ.norm(Radon) - φ(v.abs()) - dp) / r)
-};
-μ *= 1.0 - s;
-μ += δ * s;
 // The stop_target is only needed for the type system.
 AlgIteratorOptions{ max_iter : 1, .. config.inner.iterator_options}.stop_target(0.0)
 }
 };
-inner_iters += optimise_weights(&mut μ, opA, b, α, &findim_data, &config.inner, inner_it);
+inner_iters += reg.optimise_weights(&mut μ, opA, b, &findim_data, &config.inner, inner_it);
 // Merge spikes and update residual for next step and `if_verbose` below.
 let n_before_merge = μ.len();
 residual = μ.merge_spikes_fitness(config.merging,
 |μ̃| opA.apply(μ̃) - b,
 format!("iter {} start", state.iteration()), &g,
 "".to_string(), None::<&A::PreadjointCodomain>,
 None, &μ
 );
 let res = IterInfo {
-value : residual.norm2_squared_div2() + α * μ.norm(Radon),
+value : residual.norm2_squared_div2() + reg.apply(&μ),
 n_spikes : μ.len(),
 inner_iters,
 this_iters,
 merged,
 pruned,
 // Return final iterate
 μ
 }
+//
+// Deprecated interface
+//
+#[deprecated(note = "Use `pointsource_fw_reg`")]
+pub fn pointsource_fw<'a, F, I, A, GA, BTA, S, const N : usize>(
+opA : &'a A,
+b : &A::Observable,
+α : F,
+//domain : Cube<F, N>,
+config : &FWConfig<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>>,
+BTA : BTSearch<F, N, Data=usize, Agg=Bounds<F>>,
+S: 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_fw_reg(opA, b, NonnegRadonRegTerm(α), config, iterator, plotter)
+}
+#[deprecated(note = "Use `WeightOptim::optimise_weights`")]
+pub fn optimise_weights<'a, F, A, I, const N : usize>(
+μ : &mut DiscreteMeasure<Loc<F, N>, F>,
+opA : &'a A,
+b : &A::Observable,
+α : F,
+findim_data : &FindimData<F>,
+inner : &InnerSettings<F>,
+iterator : I
+) -> usize
+where F : Float + ToNalgebraRealField,
+I : AlgIteratorFactory<F>,
+A : ForwardModel<Loc<F, N>, F>
+{
+NonnegRadonRegTerm(α).optimise_weights(μ, opA, b, findim_data, inner, iterator)
+}

comparison: src/frank_wolfe.rs

src/frank_wolfe.rs

Mercurial > repos > pointsource_algs / file comparison

comparison: src/frank_wolfe.rs

src/frank_wolfe.rs