Mercurial > repos > pointsource_algs / file comparison

comparison: src/weight_optim.rs

src/weight_optim.rs

changeset 2
7a953a87b6c1
equal deleted inserted replaced
0:eb3c7813b67a 2:7a953a87b6c1
1 //! Weight optimisation routines
2
3 use numeric_literals::replace_float_literals;
4 use alg_tools::nalgebra_support::ToNalgebraRealField;
5 use alg_tools::loc::Loc;
6 use alg_tools::iterate::AlgIteratorFactory;
7 use crate::types::*;
8 use crate::subproblem::*;
9 use crate::measures::*;
10 use crate::forward_model::ForwardModel;
11
12 /// Helper struct for pre-initialising the finite-dimensional subproblems solver
13 /// [`prepare_optimise_weights_l2`].
14 ///
15 /// The pre-initialisation is done by [`prepare_optimise_weights_l2`].
16 pub struct FindimData<F : Float> {
17 opAnorm_squared : F
18 }
19
20 /// Return a pre-initialisation struct for [`prepare_optimise_weights_l2`].
21 ///
22 /// The parameter `opA` is the forward operator $A$.
23 pub fn prepare_optimise_weights<F, A, const N : usize>(opA : &A) -> FindimData<F>
24 where F : Float + ToNalgebraRealField,
25 A : ForwardModel<Loc<F, N>, F> {
26 FindimData{
27 opAnorm_squared : opA.opnorm_bound().powi(2)
28 }
29 }
30
31 /// Solve the finite-dimensional weight optimisation problem for the 2-norm-squared data fidelity
32 /// point source localisation problem.
33 ///
34 /// That is, we minimise
35 /// <div>$$
36 /// μ ↦ \frac{1}{2}\|Aμ-b\|_2^2 + α\|μ\|_ℳ + δ_{≥ 0}(μ)
37 /// $$</div>
38 /// only with respect to the weights of $μ$.
39 ///
40 /// The parameter `μ` is the discrete measure whose weights are to be optimised.
41 /// The `opA` parameter is the forward operator $A$, while `b`$ and `α` are as in the
42 /// objective above. The method parameter are set in `inner` (see [`InnerSettings`]), while
43 /// `iterator` is used to iterate the steps of the method, and `plotter` may be used to
44 /// save intermediate iteration states as images. The parameter `findim_data` should be
45 /// prepared using [`prepare_optimise_weights`]:
46 ///
47 /// Returns the number of iterations taken by the method configured in `inner`.
48 pub fn optimise_weights_l2<'a, F, A, I, const N : usize>(
49 μ : &mut DiscreteMeasure<Loc<F, N>, F>,
50 opA : &'a A,
51 b : &A::Observable,
52 α : F,
53 findim_data : &FindimData<F>,
54 inner : &InnerSettings<F>,
55 iterator : I
56 ) -> usize
57 where F : Float + ToNalgebraRealField,
58 I : AlgIteratorFactory<F>,
59 A : ForwardModel<Loc<F, N>, F>
60 {
61 // Form and solve finite-dimensional subproblem.
62 let (Ã, g̃) = opA.findim_quadratic_model(&μ, b);
63 let mut x = μ.masses_dvector();
64
65 // `inner_τ1` is based on an estimate of the operator norm of $A$ from ℳ(Ω) to
66 // ℝ^n. This estimate is a good one for the matrix norm from ℝ^m to ℝ^n when the
67 // former is equipped with the 1-norm. We need the 2-norm. To pass from 1-norm to
68 // 2-norm, we estimate
69 // ‖A‖_{2,2} := sup_{‖x‖_2 ≤ 1} ‖Ax‖_2 ≤ sup_{‖x‖_1 ≤ C} ‖Ax‖_2
70 // = C sup_{‖x‖_1 ≤ 1} ‖Ax‖_2 = C ‖A‖_{1,2},
71 // where C = √m satisfies ‖x‖_1 ≤ C ‖x‖_2. Since we are intested in ‖A_*A‖, no
72 // square root is needed when we scale:
73 let inner_τ = inner.τ0 / (findim_data.opAnorm_squared * F::cast_from(μ.len()));
74 let iters = quadratic_nonneg(inner.method, &Ã, &g̃, α, &mut x, inner_τ, iterator);
75 // Update masses of μ based on solution of finite-dimensional subproblem.
76 μ.set_masses_dvector(&x);
77
78 iters
79 }
80
81 /// Solve the finite-dimensional weight optimisation problem for the 2-norm-squared data fidelity
82 /// point source localisation problem.
83 ///
84 /// That is, we minimise
85 /// <div>$$
86 /// μ ↦ \|Aμ-b\|_1 + α\|μ\|_ℳ + δ_{≥ 0}(μ)
87 /// $$</div>
88 /// only with respect to the weights of $μ$.
89 ///
90 /// The parameter `μ` is the discrete measure whose weights are to be optimised.
91 /// The `opA` parameter is the forward operator $A$, while `b`$ and `α` are as in the
92 /// objective above. The method parameter are set in `inner` (see [`InnerSettings`]), while
93 /// `iterator` is used to iterate the steps of the method, and `plotter` may be used to
94 /// save intermediate iteration states as images.
95 ///
96 /// Returns the number of iterations taken by the method configured in `inner`.
97 #[replace_float_literals(F::cast_from(literal))]
98 pub fn optimise_weights_l1<'a, F, A, I, const N : usize>(
99 μ : &mut DiscreteMeasure<Loc<F, N>, F>,
100 opA : &'a A,
101 b : &A::Observable,
102 α : F,
103 findim_data : &FindimData<F>,
104 inner : &InnerSettings<F>,
105 iterator : I
106 ) -> usize
107 where F : Float + ToNalgebraRealField,
108 I : AlgIteratorFactory<F>,
109 A : ForwardModel<Loc<F, N>, F>
110 {
111 // Form and solve finite-dimensional subproblem.
112 let à = opA.findim_model(&μ);
113 let mut x = μ.masses_dvector();
114
115 // `inner_τ1` is based on an estimate of the operator norm of $A$ from ℳ(Ω) to
116 // ℝ^n. This estimate is a good one for the matrix norm from ℝ^m to ℝ^n when the
117 // former is equipped with the 1-norm. We need the 2-norm. To pass from 1-norm to
118 // 2-norm, we estimate
119 // ‖A‖_{2,2} := sup_{‖x‖_2 ≤ 1} ‖Ax‖_2 ≤ sup_{‖x‖_1 ≤ C} ‖Ax‖_2
120 // = C sup_{‖x‖_1 ≤ 1} ‖Ax‖_2 = C ‖A‖_{1,2},
121 // where C = √m satisfies ‖x‖_1 ≤ C ‖x‖_2. Since we are intested in ‖A_*A‖, no
122 // square root is needed when we scale:
123 let l = (findim_data.opAnorm_squared * F::cast_from(μ.len())).sqrt();
124 let inner_σ = (0.99 / inner.τ0) / l;
125 let inner_τ = inner.τ0 / l;
126 let iters = l1_nonneg(L1InnerMethod::PDPS, &Ã, b, α, &mut x, inner_τ, inner_σ, iterator);
127 // Update masses of μ based on solution of finite-dimensional subproblem.
128 μ.set_masses_dvector(&x);
129
130 iters
131 }

Mercurial Repository: pointsource_algs

mercurial