Mercurial > repos > pointsource_algs / file comparison

comparison: src/frank_wolfe.rs

src/frank_wolfe.rs

changeset 2
7a953a87b6c1
parent 0
eb3c7813b67a
equal deleted inserted replaced
0:eb3c7813b67a 2:7a953a87b6c1
50 SpikeMergingMethod, 50 SpikeMergingMethod,
51 SpikeMerging, 51 SpikeMerging,
52 }; 52 };
53 use crate::forward_model::ForwardModel; 53 use crate::forward_model::ForwardModel;
54 #[allow(unused_imports)] // Used in documentation 54 #[allow(unused_imports)] // Used in documentation
55 use crate::subproblem::{ 55 use crate::subproblem::InnerSettings;
56 quadratic_nonneg, 56 use crate::weight_optim::{
57 InnerSettings, 57 prepare_optimise_weights,
58 InnerMethod, 58 optimise_weights_l2
59 }; 59 };
60 use crate::tolerance::Tolerance; 60 use crate::tolerance::Tolerance;
61 use crate::plot::{ 61 use crate::plot::{
62 SeqPlotter, 62 SeqPlotter,
63 Plotting, 63 Plotting,
100 inner : Default::default(), 100 inner : Default::default(),
101 variant : FWVariant::FullyCorrective, 101 variant : FWVariant::FullyCorrective,
102 merging : Default::default(), 102 merging : Default::default(),
103 } 103 }
104 } 104 }
105 }
106
107 /// Helper struct for pre-initialising the finite-dimensional subproblems solver
108 /// [`prepare_optimise_weights`].
109 ///
110 /// The pre-initialisation is done by [`prepare_optimise_weights`].
111 pub struct FindimData<F : Float> {
112 opAnorm_squared : F
113 }
114
115 /// Return a pre-initialisation struct for [`prepare_optimise_weights`].
116 ///
117 /// The parameter `opA` is the forward operator $A$.
118 pub fn prepare_optimise_weights<F, A, const N : usize>(opA : &A) -> FindimData<F>
119 where F : Float + ToNalgebraRealField,
120 A : ForwardModel<Loc<F, N>, F> {
121 FindimData{
122 opAnorm_squared : opA.opnorm_bound().powi(2)
123 }
124 }
125
126 /// Solve the finite-dimensional weight optimisation problem for the 2-norm-squared data fidelity
127 /// point source localisation problem.
128 ///
129 /// That is, we minimise
130 /// <div>$$
131 /// μ ↦ \frac{1}{2}\|Aμ-b\|_w^2 + α\|μ\|_ℳ + δ_{≥ 0}(μ)
132 /// $$</div>
133 /// only with respect to the weights of $μ$.
134 ///
135 /// The parameter `μ` is the discrete measure whose weights are to be optimised.
136 /// The `opA` parameter is the forward operator $A$, while `b`$ and `α` are as in the
137 /// objective above. The method parameter are set in `inner` (see [`InnerSettings`]), while
138 /// `iterator` is used to iterate the steps of the method, and `plotter` may be used to
139 /// save intermediate iteration states as images. The parameter `findim_data` should be
140 /// prepared using [`prepare_optimise_weights`]:
141 ///
142 /// Returns the number of iterations taken by the method configured in `inner`.
143 pub fn optimise_weights<'a, F, A, I, const N : usize>(
144 μ : &mut DiscreteMeasure<Loc<F, N>, F>,
145 opA : &'a A,
146 b : &A::Observable,
147 α : F,
148 findim_data : &FindimData<F>,
149 inner : &InnerSettings<F>,
150 iterator : I
151 ) -> usize
152 where F : Float + ToNalgebraRealField,
153 I : AlgIteratorFactory<F>,
154 A : ForwardModel<Loc<F, N>, F>
155 {
156 // Form and solve finite-dimensional subproblem.
157 let (Ã, g̃) = opA.findim_quadratic_model(&μ, b);
158 let mut x = μ.masses_dvector();
159
160 // `inner_τ1` is based on an estimate of the operator norm of $A$ from ℳ(Ω) to
161 // ℝ^n. This estimate is a good one for the matrix norm from ℝ^m to ℝ^n when the
162 // former is equipped with the 1-norm. We need the 2-norm. To pass from 1-norm to
163 // 2-norm, we estimate
164 // ‖A‖_{2,2} := sup_{‖x‖_2 ≤ 1} ‖Ax‖_2 ≤ sup_{‖x‖_1 ≤ C} ‖Ax‖_2
165 // = C sup_{‖x‖_1 ≤ 1} ‖Ax‖_2 = C ‖A‖_{1,2},
166 // where C = √m satisfies ‖x‖_1 ≤ C ‖x‖_2. Since we are intested in ‖A_*A‖, no
167 // square root is needed when we scale:
168 let inner_τ = inner.τ0 / (findim_data.opAnorm_squared * F::cast_from(μ.len()));
169 let iters = quadratic_nonneg(inner.method, &Ã, &g̃, α, &mut x, inner_τ, iterator);
170 // Update masses of μ based on solution of finite-dimensional subproblem.
171 μ.set_masses_dvector(&x);
172
173 iters
174 } 105 }
175 106
176 /// Solve point source localisation problem using a conditional gradient method 107 /// Solve point source localisation problem using a conditional gradient method
177 /// for the 2-norm-squared data fidelity, i.e., the problem 108 /// for the 2-norm-squared data fidelity, i.e., the problem
178 /// <div>$$ 109 /// <div>$$
277 // The stop_target is only needed for the type system. 208 // The stop_target is only needed for the type system.
278 AlgIteratorOptions{ max_iter : 1, .. config.inner.iterator_options}.stop_target(0.0) 209 AlgIteratorOptions{ max_iter : 1, .. config.inner.iterator_options}.stop_target(0.0)
279 } 210 }
280 }; 211 };
281 212
282 inner_iters += optimise_weights(&mut μ, opA, b, α, &findim_data, &config.inner, inner_it); 213 inner_iters += optimise_weights_l2(&mut μ, opA, b, α, &findim_data, &config.inner, inner_it);
283 214
284 // Merge spikes and update residual for next step and `if_verbose` below. 215 // Merge spikes and update residual for next step and `if_verbose` below.
285 let n_before_merge = μ.len(); 216 let n_before_merge = μ.len();
286 residual = μ.merge_spikes_fitness(config.merging, 217 residual = μ.merge_spikes_fitness(config.merging,
287 |μ̃| opA.apply(μ̃) - b, 218 |μ̃| opA.apply(μ̃) - b,

Mercurial Repository: pointsource_algs

mercurial