Tue, 29 Nov 2022 15:36:12 +0200
fubar
//! 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 }