--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/weight_optim.rs Tue Nov 29 15:36:12 2022 +0200
@@ -0,0 +1,131 @@
+//! Weight optimisation routines
+
+use numeric_literals::replace_float_literals;
+use alg_tools::nalgebra_support::ToNalgebraRealField;
+use alg_tools::loc::Loc;
+use alg_tools::iterate::AlgIteratorFactory;
+use crate::types::*;
+use crate::subproblem::*;
+use crate::measures::*;
+use crate::forward_model::ForwardModel;
+
+/// Helper struct for pre-initialising the finite-dimensional subproblems solver
+/// [`prepare_optimise_weights_l2`].
+///
+/// The pre-initialisation is done by [`prepare_optimise_weights_l2`].
+pub struct FindimData<F : Float> {
+ opAnorm_squared : F
+}
+
+/// Return a pre-initialisation struct for [`prepare_optimise_weights_l2`].
+///
+/// The parameter `opA` is the forward operator $A$.
+pub fn prepare_optimise_weights<F, A, const N : usize>(opA : &A) -> FindimData<F>
+where F : Float + ToNalgebraRealField,
+ A : ForwardModel<Loc<F, N>, F> {
+ FindimData{
+ opAnorm_squared : opA.opnorm_bound().powi(2)
+ }
+}
+
+/// 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\|_2^2 + α\|μ\|_ℳ + δ_{≥ 0}(μ)
+/// $$</div>
+/// only with respect to the weights of $μ$.
+///
+/// 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 [`prepare_optimise_weights`]:
+///
+/// Returns the number of iterations taken by the method configured in `inner`.
+pub fn optimise_weights_l2<'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>
+{
+ // Form and solve finite-dimensional subproblem.
+ let (Ã, 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, &Ã, &g̃, α, &mut x, inner_τ, iterator);
+ // Update masses of μ based on solution of finite-dimensional subproblem.
+ μ.set_masses_dvector(&x);
+
+ iters
+}
+
+/// Solve the finite-dimensional weight optimisation problem for the 2-norm-squared data fidelity
+/// point source localisation problem.
+///
+/// That is, we minimise
+/// <div>$$
+/// μ ↦ \|Aμ-b\|_1 + α\|μ\|_ℳ + δ_{≥ 0}(μ)
+/// $$</div>
+/// only with respect to the weights of $μ$.
+///
+/// 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.
+///
+/// Returns the number of iterations taken by the method configured in `inner`.
+#[replace_float_literals(F::cast_from(literal))]
+pub fn optimise_weights_l1<'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>
+{
+ // Form and solve finite-dimensional subproblem.
+ let à = opA.findim_model(&μ);
+ 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 l = (findim_data.opAnorm_squared * F::cast_from(μ.len())).sqrt();
+ let inner_σ = (0.99 / inner.τ0) / l;
+ let inner_τ = inner.τ0 / l;
+ let iters = l1_nonneg(L1InnerMethod::PDPS, &Ã, b, α, &mut x, inner_τ, inner_σ, iterator);
+ // Update masses of μ based on solution of finite-dimensional subproblem.
+ μ.set_masses_dvector(&x);
+
+ iters
+}