Mercurial > repos > pointsource_algs / file comparison

comparison: src/frank_wolfe.rs

src/frank_wolfe.rs

branch
dev
changeset 61
4f468d35fa29
parent 39
6316d68b58af
child 63
7a8a55fd41c0
equal deleted inserted replaced
60:9738b51d90d7 61:4f468d35fa29
11 11
12 * Bredies K., Pikkarainen H. - _Inverse problems in spaces of measures_, 12 * Bredies K., Pikkarainen H. - _Inverse problems in spaces of measures_,
13 DOI: [10.1051/cocv/2011205](https://doi.org/0.1051/cocv/2011205). 13 DOI: [10.1051/cocv/2011205](https://doi.org/0.1051/cocv/2011205).
14 */ 14 */
15 15
16 use nalgebra::{DMatrix, DVector};
16 use numeric_literals::replace_float_literals; 17 use numeric_literals::replace_float_literals;
17 use nalgebra::{DMatrix, DVector}; 18 use serde::{Deserialize, Serialize};
18 use serde::{Serialize, Deserialize};
19 //use colored::Colorize; 19 //use colored::Colorize;
20 20 use crate::dataterm::QuadraticDataTerm;
21 use alg_tools::iterate::{
22 AlgIteratorFactory,
23 AlgIteratorOptions,
24 ValueIteratorFactory,
25 };
26 use alg_tools::euclidean::Euclidean;
27 use alg_tools::norms::Norm;
28 use alg_tools::linops::Mapping;
29 use alg_tools::sets::Cube;
30 use alg_tools::loc::Loc;
31 use alg_tools::bisection_tree::{
32 BTFN,
33 Bounds,
34 BTNodeLookup,
35 BTNode,
36 BTSearch,
37 P2Minimise,
38 SupportGenerator,
39 LocalAnalysis,
40 };
41 use alg_tools::mapping::RealMapping;
42 use alg_tools::nalgebra_support::ToNalgebraRealField;
43 use alg_tools::norms::L2;
44
45 use crate::types::*;
46 use crate::measures::{
47 RNDM,
48 DiscreteMeasure,
49 DeltaMeasure,
50 Radon,
51 };
52 use crate::measures::merging::{
53 SpikeMergingMethod,
54 SpikeMerging,
55 };
56 use crate::forward_model::ForwardModel; 21 use crate::forward_model::ForwardModel;
22 use crate::measures::merging::{SpikeMerging, SpikeMergingMethod};
23 use crate::measures::{DeltaMeasure, DiscreteMeasure, Radon, RNDM};
24 use crate::plot::Plotter;
25 use crate::regularisation::{NonnegRadonRegTerm, RadonRegTerm, RegTerm};
57 #[allow(unused_imports)] // Used in documentation 26 #[allow(unused_imports)] // Used in documentation
58 use crate::subproblem::{ 27 use crate::subproblem::{
59 unconstrained::quadratic_unconstrained, 28 nonneg::quadratic_nonneg, unconstrained::quadratic_unconstrained, InnerMethod, InnerSettings,
60 nonneg::quadratic_nonneg,
61 InnerSettings,
62 InnerMethod,
63 }; 29 };
64 use crate::tolerance::Tolerance; 30 use crate::tolerance::Tolerance;
65 use crate::plot::{ 31 use crate::types::*;
66 SeqPlotter, 32 use alg_tools::bisection_tree::P2Minimise;
67 Plotting, 33 use alg_tools::bounds::MinMaxMapping;
68 PlotLookup 34 use alg_tools::error::DynResult;
69 }; 35 use alg_tools::euclidean::Euclidean;
70 use crate::regularisation::{ 36 use alg_tools::instance::Instance;
71 NonnegRadonRegTerm, 37 use alg_tools::iterate::{AlgIteratorFactory, AlgIteratorOptions, ValueIteratorFactory};
72 RadonRegTerm, 38 use alg_tools::linops::Mapping;
73 RegTerm 39 use alg_tools::loc::Loc;
74 }; 40 use alg_tools::nalgebra_support::ToNalgebraRealField;
41 use alg_tools::norms::Norm;
42 use alg_tools::norms::L2;
43 use alg_tools::sets::Cube;
75 44
76 /// Settings for [`pointsource_fw_reg`]. 45 /// Settings for [`pointsource_fw_reg`].
77 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)] 46 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
78 #[serde(default)] 47 #[serde(default)]
79 pub struct FWConfig<F : Float> { 48 pub struct FWConfig<F: Float> {
80 /// Tolerance for branch-and-bound new spike location discovery 49 /// Tolerance for branch-and-bound new spike location discovery
81 pub tolerance : Tolerance<F>, 50 pub tolerance: Tolerance<F>,
82 /// Inner problem solution configuration. Has to have `method` set to [`InnerMethod::FB`] 51 /// Inner problem solution configuration. Has to have `method` set to [`InnerMethod::FB`]
83 /// as the conditional gradient subproblems' optimality conditions do not in general have an 52 /// as the conditional gradient subproblems' optimality conditions do not in general have an
84 /// invertible Newton derivative for SSN. 53 /// invertible Newton derivative for SSN.
85 pub inner : InnerSettings<F>, 54 pub inner: InnerSettings<F>,
86 /// Variant of the conditional gradient method 55 /// Variant of the conditional gradient method
87 pub variant : FWVariant, 56 pub variant: FWVariant,
88 /// Settings for branch and bound refinement when looking for predual maxima 57 /// Settings for branch and bound refinement when looking for predual maxima
89 pub refinement : RefinementSettings<F>, 58 pub refinement: RefinementSettings<F>,
90 /// Spike merging heuristic 59 /// Spike merging heuristic
91 pub merging : SpikeMergingMethod<F>, 60 pub merging: SpikeMergingMethod<F>,
92 } 61 }
93 62
94 /// Conditional gradient method variant; see also [`FWConfig`]. 63 /// Conditional gradient method variant; see also [`FWConfig`].
95 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)] 64 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
96 #[allow(dead_code)] 65 #[allow(dead_code)]
99 FullyCorrective, 68 FullyCorrective,
100 /// Bredies–Pikkarainen. Forces `FWConfig.inner.max_iter = 1`. 69 /// Bredies–Pikkarainen. Forces `FWConfig.inner.max_iter = 1`.
101 Relaxed, 70 Relaxed,
102 } 71 }
103 72
104 impl<F : Float> Default for FWConfig<F> { 73 impl<F: Float> Default for FWConfig<F> {
105 fn default() -> Self { 74 fn default() -> Self {
106 FWConfig { 75 FWConfig {
107 tolerance : Default::default(), 76 tolerance: Default::default(),
108 refinement : Default::default(), 77 refinement: Default::default(),
109 inner : Default::default(), 78 inner: Default::default(),
110 variant : FWVariant::FullyCorrective, 79 variant: FWVariant::FullyCorrective,
111 merging : SpikeMergingMethod { enabled : true, ..Default::default() }, 80 merging: SpikeMergingMethod { enabled: true, ..Default::default() },
112 } 81 }
113 } 82 }
114 } 83 }
115 84
116 pub trait FindimQuadraticModel<Domain, F> : ForwardModel<DiscreteMeasure<Domain, F>, F> 85 pub trait FindimQuadraticModel<Domain, F>: ForwardModel<DiscreteMeasure<Domain, F>, F>
117 where 86 where
118 F : Float + ToNalgebraRealField, 87 F: Float + ToNalgebraRealField,
119 Domain : Clone + PartialEq, 88 Domain: Clone + PartialEq,
120 { 89 {
121 /// Return A_*A and A_* b 90 /// Return A_*A and A_* b
122 fn findim_quadratic_model( 91 fn findim_quadratic_model(
123 &self, 92 &self,
124 μ : &DiscreteMeasure<Domain, F>, 93 μ: &DiscreteMeasure<Domain, F>,
125 b : &Self::Observable 94 b: &Self::Observable,
126 ) -> (DMatrix<F::MixedType>, DVector<F::MixedType>); 95 ) -> (DMatrix<F::MixedType>, DVector<F::MixedType>);
127 } 96 }
128 97
129 /// Helper struct for pre-initialising the finite-dimensional subproblem solver. 98 /// Helper struct for pre-initialising the finite-dimensional subproblem solver.
130 pub struct FindimData<F : Float> { 99 pub struct FindimData<F: Float> {
131 /// ‖A‖^2 100 /// ‖A‖^2
132 opAnorm_squared : F, 101 opAnorm_squared: F,
133 /// Bound $M_0$ from the Bredies–Pikkarainen article. 102 /// Bound $M_0$ from the Bredies–Pikkarainen article.
134 m0 : F 103 m0: F,
135 } 104 }
136 105
137 /// Trait for finite dimensional weight optimisation. 106 /// Trait for finite dimensional weight optimisation.
138 pub trait WeightOptim< 107 pub trait WeightOptim<
139 F : Float + ToNalgebraRealField, 108 F: Float + ToNalgebraRealField,
140 A : ForwardModel<RNDM<F, N>, F>, 109 A: ForwardModel<RNDM<N, F>, F>,
141 I : AlgIteratorFactory<F>, 110 I: AlgIteratorFactory<F>,
142 const N : usize 111 const N: usize,
143 > { 112 >
144 113 {
145 /// Return a pre-initialisation struct for [`Self::optimise_weights`]. 114 /// Return a pre-initialisation struct for [`Self::optimise_weights`].
146 /// 115 ///
147 /// The parameter `opA` is the forward operator $A$. 116 /// The parameter `opA` is the forward operator $A$.
148 fn prepare_optimise_weights(&self, opA : &A, b : &A::Observable) -> FindimData<F>; 117 fn prepare_optimise_weights(&self, opA: &A, b: &A::Observable) -> DynResult<FindimData<F>>;
149 118
150 /// Solve the finite-dimensional weight optimisation problem for the 2-norm-squared data fidelity 119 /// Solve the finite-dimensional weight optimisation problem for the 2-norm-squared data fidelity
151 /// point source localisation problem. 120 /// point source localisation problem.
152 /// 121 ///
153 /// That is, we minimise 122 /// That is, we minimise
164 /// prepared using [`Self::prepare_optimise_weights`]: 133 /// prepared using [`Self::prepare_optimise_weights`]:
165 /// 134 ///
166 /// Returns the number of iterations taken by the method configured in `inner`. 135 /// Returns the number of iterations taken by the method configured in `inner`.
167 fn optimise_weights<'a>( 136 fn optimise_weights<'a>(
168 &self, 137 &self,
169 μ : &mut RNDM<F, N>, 138 μ: &mut RNDM<N, F>,
170 opA : &'a A, 139 opA: &'a A,
171 b : &A::Observable, 140 b: &A::Observable,
172 findim_data : &FindimData<F>, 141 findim_data: &FindimData<F>,
173 inner : &InnerSettings<F>, 142 inner: &InnerSettings<F>,
174 iterator : I 143 iterator: I,
175 ) -> usize; 144 ) -> usize;
176 } 145 }
177 146
178 /// Trait for regularisation terms supported by [`pointsource_fw_reg`]. 147 /// Trait for regularisation terms supported by [`pointsource_fw_reg`].
179 pub trait RegTermFW< 148 pub trait RegTermFW<
180 F : Float + ToNalgebraRealField, 149 F: Float + ToNalgebraRealField,
181 A : ForwardModel<RNDM<F, N>, F>, 150 A: ForwardModel<RNDM<N, F>, F>,
182 I : AlgIteratorFactory<F>, 151 I: AlgIteratorFactory<F>,
183 const N : usize 152 const N: usize,
184 > : RegTerm<F, N> 153 >: RegTerm<Loc<N, F>, F> + WeightOptim<F, A, I, N> + Mapping<RNDM<N, F>, Codomain = F>
185 + WeightOptim<F, A, I, N> 154 {
186 + Mapping<RNDM<F, N>, Codomain = F> {
187
188 /// With $g = A\_\*(Aμ-b)$, returns $(x, g(x))$ for $x$ a new point to be inserted 155 /// With $g = A\_\*(Aμ-b)$, returns $(x, g(x))$ for $x$ a new point to be inserted
189 /// into $μ$, as determined by the regulariser. 156 /// into $μ$, as determined by the regulariser.
190 /// 157 ///
191 /// The parameters `refinement_tolerance` and `max_steps` are passed to relevant 158 /// The parameters `refinement_tolerance` and `max_steps` are passed to relevant
192 /// [`BTFN`] minimisation and maximisation routines. 159 /// [`MinMaxMapping`] minimisation and maximisation routines.
193 fn find_insertion( 160 fn find_insertion(
194 &self, 161 &self,
195 g : &mut A::PreadjointCodomain, 162 g: &mut A::PreadjointCodomain,
196 refinement_tolerance : F, 163 refinement_tolerance: F,
197 max_steps : usize 164 max_steps: usize,
198 ) -> (Loc<F, N>, F); 165 ) -> (Loc<N, F>, F);
199 166
200 /// Insert point `ξ` into `μ` for the relaxed algorithm from Bredies–Pikkarainen. 167 /// Insert point `ξ` into `μ` for the relaxed algorithm from Bredies–Pikkarainen.
201 fn relaxed_insert<'a>( 168 fn relaxed_insert<'a>(
202 &self, 169 &self,
203 μ : &mut RNDM<F, N>, 170 μ: &mut RNDM<N, F>,
204 g : &A::PreadjointCodomain, 171 g: &A::PreadjointCodomain,
205 opA : &'a A, 172 opA: &'a A,
206 ξ : Loc<F, N>, 173 ξ: Loc<N, F>,
207 v_ξ : F, 174 v_ξ: F,
208 findim_data : &FindimData<F> 175 findim_data: &FindimData<F>,
209 ); 176 );
210 } 177 }
211 178
212 #[replace_float_literals(F::cast_from(literal))] 179 #[replace_float_literals(F::cast_from(literal))]
213 impl<F : Float + ToNalgebraRealField, A, I, const N : usize> WeightOptim<F, A, I, N> 180 impl<F: Float + ToNalgebraRealField, A, I, const N: usize> WeightOptim<F, A, I, N>
214 for RadonRegTerm<F> 181 for RadonRegTerm<F>
215 where I : AlgIteratorFactory<F>, 182 where
216 A : FindimQuadraticModel<Loc<F, N>, F> { 183 I: AlgIteratorFactory<F>,
217 184 A: FindimQuadraticModel<Loc<N, F>, F>,
218 fn prepare_optimise_weights(&self, opA : &A, b : &A::Observable) -> FindimData<F> { 185 {
219 FindimData{ 186 fn prepare_optimise_weights(&self, opA: &A, b: &A::Observable) -> DynResult<FindimData<F>> {
220 opAnorm_squared : opA.opnorm_bound(Radon, L2).powi(2), 187 Ok(FindimData {
221 m0 : b.norm2_squared() / (2.0 * self.α()), 188 opAnorm_squared: opA.opnorm_bound(Radon, L2)?.powi(2),
222 } 189 m0: b.norm2_squared() / (2.0 * self.α()),
190 })
223 } 191 }
224 192
225 fn optimise_weights<'a>( 193 fn optimise_weights<'a>(
226 &self, 194 &self,
227 μ : &mut RNDM<F, N>, 195 μ: &mut RNDM<N, F>,
228 opA : &'a A, 196 opA: &'a A,
229 b : &A::Observable, 197 b: &A::Observable,
230 findim_data : &FindimData<F>, 198 findim_data: &FindimData<F>,
231 inner : &InnerSettings<F>, 199 inner: &InnerSettings<F>,
232 iterator : I 200 iterator: I,
233 ) -> usize { 201 ) -> usize {
234
235 // Form and solve finite-dimensional subproblem. 202 // Form and solve finite-dimensional subproblem.
236 let (Ã, g̃) = opA.findim_quadratic_model(&μ, b); 203 let (Ã, g̃) = opA.findim_quadratic_model(&μ, b);
237 let mut x = μ.masses_dvector(); 204 let mut x = μ.masses_dvector();
238 205
239 // `inner_τ1` is based on an estimate of the operator norm of $A$ from ℳ(Ω) to 206 // `inner_τ1` is based on an estimate of the operator norm of $A$ from ℳ(Ω) to
243 // ‖A‖_{2,2} := sup_{‖x‖_2 ≤ 1} ‖Ax‖_2 ≤ sup_{‖x‖_1 ≤ C} ‖Ax‖_2 210 // ‖A‖_{2,2} := sup_{‖x‖_2 ≤ 1} ‖Ax‖_2 ≤ sup_{‖x‖_1 ≤ C} ‖Ax‖_2
244 // = C sup_{‖x‖_1 ≤ 1} ‖Ax‖_2 = C ‖A‖_{1,2}, 211 // = C sup_{‖x‖_1 ≤ 1} ‖Ax‖_2 = C ‖A‖_{1,2},
245 // where C = √m satisfies ‖x‖_1 ≤ C ‖x‖_2. Since we are intested in ‖A_*A‖, no 212 // where C = √m satisfies ‖x‖_1 ≤ C ‖x‖_2. Since we are intested in ‖A_*A‖, no
246 // square root is needed when we scale: 213 // square root is needed when we scale:
247 let normest = findim_data.opAnorm_squared * F::cast_from(μ.len()); 214 let normest = findim_data.opAnorm_squared * F::cast_from(μ.len());
248 let iters = quadratic_unconstrained(&Ã, &g̃, self.α(), &mut x, 215 let iters = quadratic_unconstrained(&Ã, &g̃, self.α(), &mut x, normest, inner, iterator);
249 normest, inner, iterator);
250 // Update masses of μ based on solution of finite-dimensional subproblem. 216 // Update masses of μ based on solution of finite-dimensional subproblem.
251 μ.set_masses_dvector(&x); 217 μ.set_masses_dvector(&x);
252 218
253 iters 219 iters
254 } 220 }
255 } 221 }
256 222
257 #[replace_float_literals(F::cast_from(literal))] 223 #[replace_float_literals(F::cast_from(literal))]
258 impl<F : Float + ToNalgebraRealField, A, I, S, GA, BTA, const N : usize> RegTermFW<F, A, I, N> 224 impl<F: Float + ToNalgebraRealField, A, I, const N: usize> RegTermFW<F, A, I, N> for RadonRegTerm<F>
259 for RadonRegTerm<F>
260 where 225 where
261 Cube<F, N> : P2Minimise<Loc<F, N>, F>, 226 Cube<N, F>: P2Minimise<Loc<N, F>, F>,
262 I : AlgIteratorFactory<F>, 227 I: AlgIteratorFactory<F>,
263 S: RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>, 228 A: FindimQuadraticModel<Loc<N, F>, F>,
264 GA : SupportGenerator<F, N, SupportType = S, Id = usize> + Clone, 229 A::PreadjointCodomain: MinMaxMapping<Loc<N, F>, F>,
265 A : FindimQuadraticModel<Loc<F, N>, F, PreadjointCodomain = BTFN<F, GA, BTA, N>>, 230 for<'a> &'a A::PreadjointCodomain: Instance<A::PreadjointCodomain>,
266 BTA : BTSearch<F, N, Data=usize, Agg=Bounds<F>>,
267 // FIXME: the following *should not* be needed, they are already implied 231 // FIXME: the following *should not* be needed, they are already implied
268 RNDM<F, N> : Mapping<A::PreadjointCodomain, Codomain = F>, 232 RNDM<N, F>: Mapping<A::PreadjointCodomain, Codomain = F>,
269 DeltaMeasure<Loc<F, N>, F> : Mapping<A::PreadjointCodomain, Codomain = F>, 233 DeltaMeasure<Loc<N, F>, F>: Mapping<A::PreadjointCodomain, Codomain = F>,
270 //A : Mapping<RNDM<F, N>, Codomain = A::Observable>, 234 {
271 //A : Mapping<DeltaMeasure<Loc<F, N>, F>, Codomain = A::Observable>,
272 {
273
274 fn find_insertion( 235 fn find_insertion(
275 &self, 236 &self,
276 g : &mut A::PreadjointCodomain, 237 g: &mut A::PreadjointCodomain,
277 refinement_tolerance : F, 238 refinement_tolerance: F,
278 max_steps : usize 239 max_steps: usize,
279 ) -> (Loc<F, N>, F) { 240 ) -> (Loc<N, F>, F) {
280 let (ξmax, v_ξmax) = g.maximise(refinement_tolerance, max_steps); 241 let (ξmax, v_ξmax) = g.maximise(refinement_tolerance, max_steps);
281 let (ξmin, v_ξmin) = g.minimise(refinement_tolerance, max_steps); 242 let (ξmin, v_ξmin) = g.minimise(refinement_tolerance, max_steps);
282 if v_ξmin < 0.0 && -v_ξmin > v_ξmax { 243 if v_ξmin < 0.0 && -v_ξmin > v_ξmax {
283 (ξmin, v_ξmin) 244 (ξmin, v_ξmin)
284 } else { 245 } else {
286 } 247 }
287 } 248 }
288 249
289 fn relaxed_insert<'a>( 250 fn relaxed_insert<'a>(
290 &self, 251 &self,
291 μ : &mut RNDM<F, N>, 252 μ: &mut RNDM<N, F>,
292 g : &A::PreadjointCodomain, 253 g: &A::PreadjointCodomain,
293 opA : &'a A, 254 opA: &'a A,
294 ξ : Loc<F, N>, 255 ξ: Loc<N, F>,
295 v_ξ : F, 256 v_ξ: F,
296 findim_data : &FindimData<F> 257 findim_data: &FindimData<F>,
297 ) { 258 ) {
298 let α = self.0; 259 let α = self.0;
299 let m0 = findim_data.m0; 260 let m0 = findim_data.m0;
300 let φ = |t| if t <= m0 { α * t } else { α / (2.0 * m0) * (t*t + m0 * m0) }; 261 let φ = |t| {
301 let v = if v_ξ.abs() <= α { 0.0 } else { m0 / α * v_ξ }; 262 if t <= m0 {
302 let δ = DeltaMeasure { x : ξ, α : v }; 263 α * t
264 } else {
265 α / (2.0 * m0) * (t * t + m0 * m0)
266 }
267 };
268 let v = if v_ξ.abs() <= α {
269 0.0
270 } else {
271 m0 / α * v_ξ
272 };
273 let δ = DeltaMeasure { x: ξ, α: v };
303 let dp = μ.apply(g) - δ.apply(g); 274 let dp = μ.apply(g) - δ.apply(g);
304 let d = opA.apply(&*μ) - opA.apply(δ); 275 let d = opA.apply(&*μ) - opA.apply(δ);
305 let r = d.norm2_squared(); 276 let r = d.norm2_squared();
306 let s = if r == 0.0 { 277 let s = if r == 0.0 {
307 1.0 278 1.0
308 } else { 279 } else {
309 1.0.min( (α * μ.norm(Radon) - φ(v.abs()) - dp) / r) 280 1.0.min((α * μ.norm(Radon) - φ(v.abs()) - dp) / r)
310 }; 281 };
311 *μ *= 1.0 - s; 282 *μ *= 1.0 - s;
312 *μ += δ * s; 283 *μ += δ * s;
313 } 284 }
314 } 285 }
315 286
316 #[replace_float_literals(F::cast_from(literal))] 287 #[replace_float_literals(F::cast_from(literal))]
317 impl<F : Float + ToNalgebraRealField, A, I, const N : usize> WeightOptim<F, A, I, N> 288 impl<F: Float + ToNalgebraRealField, A, I, const N: usize> WeightOptim<F, A, I, N>
318 for NonnegRadonRegTerm<F> 289 for NonnegRadonRegTerm<F>
319 where I : AlgIteratorFactory<F>, 290 where
320 A : FindimQuadraticModel<Loc<F, N>, F> { 291 I: AlgIteratorFactory<F>,
321 292 A: FindimQuadraticModel<Loc<N, F>, F>,
322 fn prepare_optimise_weights(&self, opA : &A, b : &A::Observable) -> FindimData<F> { 293 {
323 FindimData{ 294 fn prepare_optimise_weights(&self, opA: &A, b: &A::Observable) -> DynResult<FindimData<F>> {
324 opAnorm_squared : opA.opnorm_bound(Radon, L2).powi(2), 295 Ok(FindimData {
325 m0 : b.norm2_squared() / (2.0 * self.α()), 296 opAnorm_squared: opA.opnorm_bound(Radon, L2)?.powi(2),
326 } 297 m0: b.norm2_squared() / (2.0 * self.α()),
298 })
327 } 299 }
328 300
329 fn optimise_weights<'a>( 301 fn optimise_weights<'a>(
330 &self, 302 &self,
331 μ : &mut RNDM<F, N>, 303 μ: &mut RNDM<N, F>,
332 opA : &'a A, 304 opA: &'a A,
333 b : &A::Observable, 305 b: &A::Observable,
334 findim_data : &FindimData<F>, 306 findim_data: &FindimData<F>,
335 inner : &InnerSettings<F>, 307 inner: &InnerSettings<F>,
336 iterator : I 308 iterator: I,
337 ) -> usize { 309 ) -> usize {
338
339 // Form and solve finite-dimensional subproblem. 310 // Form and solve finite-dimensional subproblem.
340 let (Ã, g̃) = opA.findim_quadratic_model(&μ, b); 311 let (Ã, g̃) = opA.findim_quadratic_model(&μ, b);
341 let mut x = μ.masses_dvector(); 312 let mut x = μ.masses_dvector();
342 313
343 // `inner_τ1` is based on an estimate of the operator norm of $A$ from ℳ(Ω) to 314 // `inner_τ1` is based on an estimate of the operator norm of $A$ from ℳ(Ω) to
347 // ‖A‖_{2,2} := sup_{‖x‖_2 ≤ 1} ‖Ax‖_2 ≤ sup_{‖x‖_1 ≤ C} ‖Ax‖_2 318 // ‖A‖_{2,2} := sup_{‖x‖_2 ≤ 1} ‖Ax‖_2 ≤ sup_{‖x‖_1 ≤ C} ‖Ax‖_2
348 // = C sup_{‖x‖_1 ≤ 1} ‖Ax‖_2 = C ‖A‖_{1,2}, 319 // = C sup_{‖x‖_1 ≤ 1} ‖Ax‖_2 = C ‖A‖_{1,2},
349 // where C = √m satisfies ‖x‖_1 ≤ C ‖x‖_2. Since we are intested in ‖A_*A‖, no 320 // where C = √m satisfies ‖x‖_1 ≤ C ‖x‖_2. Since we are intested in ‖A_*A‖, no
350 // square root is needed when we scale: 321 // square root is needed when we scale:
351 let normest = findim_data.opAnorm_squared * F::cast_from(μ.len()); 322 let normest = findim_data.opAnorm_squared * F::cast_from(μ.len());
352 let iters = quadratic_nonneg(&Ã, &g̃, self.α(), &mut x, 323 let iters = quadratic_nonneg(&Ã, &g̃, self.α(), &mut x, normest, inner, iterator);
353 normest, inner, iterator);
354 // Update masses of μ based on solution of finite-dimensional subproblem. 324 // Update masses of μ based on solution of finite-dimensional subproblem.
355 μ.set_masses_dvector(&x); 325 μ.set_masses_dvector(&x);
356 326
357 iters 327 iters
358 } 328 }
359 } 329 }
360 330
361 #[replace_float_literals(F::cast_from(literal))] 331 #[replace_float_literals(F::cast_from(literal))]
362 impl<F : Float + ToNalgebraRealField, A, I, S, GA, BTA, const N : usize> RegTermFW<F, A, I, N> 332 impl<F: Float + ToNalgebraRealField, A, I, const N: usize> RegTermFW<F, A, I, N>
363 for NonnegRadonRegTerm<F> 333 for NonnegRadonRegTerm<F>
364 where 334 where
365 Cube<F, N> : P2Minimise<Loc<F, N>, F>, 335 Cube<N, F>: P2Minimise<Loc<N, F>, F>,
366 I : AlgIteratorFactory<F>, 336 I: AlgIteratorFactory<F>,
367 S: RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>, 337 A: FindimQuadraticModel<Loc<N, F>, F>,
368 GA : SupportGenerator<F, N, SupportType = S, Id = usize> + Clone, 338 A::PreadjointCodomain: MinMaxMapping<Loc<N, F>, F>,
369 A : FindimQuadraticModel<Loc<F, N>, F, PreadjointCodomain = BTFN<F, GA, BTA, N>>, 339 for<'a> &'a A::PreadjointCodomain: Instance<A::PreadjointCodomain>,
370 BTA : BTSearch<F, N, Data=usize, Agg=Bounds<F>>,
371 // FIXME: the following *should not* be needed, they are already implied 340 // FIXME: the following *should not* be needed, they are already implied
372 RNDM<F, N> : Mapping<A::PreadjointCodomain, Codomain = F>, 341 RNDM<N, F>: Mapping<A::PreadjointCodomain, Codomain = F>,
373 DeltaMeasure<Loc<F, N>, F> : Mapping<A::PreadjointCodomain, Codomain = F>, 342 DeltaMeasure<Loc<N, F>, F>: Mapping<A::PreadjointCodomain, Codomain = F>,
374 { 343 {
375
376 fn find_insertion( 344 fn find_insertion(
377 &self, 345 &self,
378 g : &mut A::PreadjointCodomain, 346 g: &mut A::PreadjointCodomain,
379 refinement_tolerance : F, 347 refinement_tolerance: F,
380 max_steps : usize 348 max_steps: usize,
381 ) -> (Loc<F, N>, F) { 349 ) -> (Loc<N, F>, F) {
382 g.maximise(refinement_tolerance, max_steps) 350 g.maximise(refinement_tolerance, max_steps)
383 } 351 }
384 352
385
386 fn relaxed_insert<'a>( 353 fn relaxed_insert<'a>(
387 &self, 354 &self,
388 μ : &mut RNDM<F, N>, 355 μ: &mut RNDM<N, F>,
389 g : &A::PreadjointCodomain, 356 g: &A::PreadjointCodomain,
390 opA : &'a A, 357 opA: &'a A,
391 ξ : Loc<F, N>, 358 ξ: Loc<N, F>,
392 v_ξ : F, 359 v_ξ: F,
393 findim_data : &FindimData<F> 360 findim_data: &FindimData<F>,
394 ) { 361 ) {
395 // This is just a verbatim copy of RadonRegTerm::relaxed_insert. 362 // This is just a verbatim copy of RadonRegTerm::relaxed_insert.
396 let α = self.0; 363 let α = self.0;
397 let m0 = findim_data.m0; 364 let m0 = findim_data.m0;
398 let φ = |t| if t <= m0 { α * t } else { α / (2.0 * m0) * (t*t + m0 * m0) }; 365 let φ = |t| {
399 let v = if v_ξ.abs() <= α { 0.0 } else { m0 / α * v_ξ }; 366 if t <= m0 {
400 let δ = DeltaMeasure { x : ξ, α : v }; 367 α * t
368 } else {
369 α / (2.0 * m0) * (t * t + m0 * m0)
370 }
371 };
372 let v = if v_ξ.abs() <= α {
373 0.0
374 } else {
375 m0 / α * v_ξ
376 };
377 let δ = DeltaMeasure { x: ξ, α: v };
401 let dp = μ.apply(g) - δ.apply(g); 378 let dp = μ.apply(g) - δ.apply(g);
402 let d = opA.apply(&*μ) - opA.apply(&δ); 379 let d = opA.apply(&*μ) - opA.apply(&δ);
403 let r = d.norm2_squared(); 380 let r = d.norm2_squared();
404 let s = if r == 0.0 { 381 let s = if r == 0.0 {
405 1.0 382 1.0
406 } else { 383 } else {
407 1.0.min( (α * μ.norm(Radon) - φ(v.abs()) - dp) / r) 384 1.0.min((α * μ.norm(Radon) - φ(v.abs()) - dp) / r)
408 }; 385 };
409 *μ *= 1.0 - s; 386 *μ *= 1.0 - s;
410 *μ += δ * s; 387 *μ += δ * s;
411 } 388 }
412 } 389 }
413
414 390
415 /// Solve point source localisation problem using a conditional gradient method 391 /// Solve point source localisation problem using a conditional gradient method
416 /// for the 2-norm-squared data fidelity, i.e., the problem 392 /// for the 2-norm-squared data fidelity, i.e., the problem
417 /// <div>$$ 393 /// <div>$$
418 /// \min_μ \frac{1}{2}\|Aμ-b\|_w^2 + G(μ), 394 /// \min_μ \frac{1}{2}\|Aμ-b\|_w^2 + G(μ),
423 /// The `opA` parameter is the forward operator $A$, while `b`$ is as in the 399 /// The `opA` parameter is the forward operator $A$, while `b`$ is as in the
424 /// objective above. The method parameter are set in `config` (see [`FWConfig`]), while 400 /// objective above. The method parameter are set in `config` (see [`FWConfig`]), while
425 /// `iterator` is used to iterate the steps of the method, and `plotter` may be used to 401 /// `iterator` is used to iterate the steps of the method, and `plotter` may be used to
426 /// save intermediate iteration states as images. 402 /// save intermediate iteration states as images.
427 #[replace_float_literals(F::cast_from(literal))] 403 #[replace_float_literals(F::cast_from(literal))]
428 pub fn pointsource_fw_reg<F, I, A, GA, BTA, S, Reg, const N : usize>( 404 pub fn pointsource_fw_reg<'a, F, I, A, Reg, Plot, const N: usize>(
429 opA : &A, 405 f: &'a QuadraticDataTerm<F, RNDM<N, F>, A>,
430 b : &A::Observable, 406 reg: &Reg,
431 reg : Reg, 407 //domain : Cube<N, F>,
432 //domain : Cube<F, N>, 408 config: &FWConfig<F>,
433 config : &FWConfig<F>, 409 iterator: I,
434 iterator : I, 410 mut plotter: Plot,
435 mut plotter : SeqPlotter<F, N>, 411 μ0 : Option<RNDM<N, F>>,
436 ) -> RNDM<F, N> 412 ) -> DynResult<RNDM<N, F>>
437 where F : Float + ToNalgebraRealField, 413 where
438 I : AlgIteratorFactory<IterInfo<F, N>>, 414 F: Float + ToNalgebraRealField,
439 for<'b> &'b A::Observable : std::ops::Neg<Output=A::Observable>, 415 I: AlgIteratorFactory<IterInfo<F>>,
440 GA : SupportGenerator<F, N, SupportType = S, Id = usize> + Clone, 416 A: ForwardModel<RNDM<N, F>, F>,
441 A : ForwardModel<RNDM<F, N>, F, PreadjointCodomain = BTFN<F, GA, BTA, N>>, 417 A::PreadjointCodomain: MinMaxMapping<Loc<N, F>, F>,
442 BTA : BTSearch<F, N, Data=usize, Agg=Bounds<F>>, 418 &'a A::PreadjointCodomain: Instance<A::PreadjointCodomain>,
443 S: RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>, 419 Cube<N, F>: P2Minimise<Loc<N, F>, F>,
444 BTNodeLookup: BTNode<F, usize, Bounds<F>, N>, 420 RNDM<N, F>: SpikeMerging<F>,
445 Cube<F, N>: P2Minimise<Loc<F, N>, F>, 421 Reg: RegTermFW<F, A, ValueIteratorFactory<F, AlgIteratorOptions>, N>,
446 PlotLookup : Plotting<N>, 422 Plot: Plotter<A::PreadjointCodomain, A::PreadjointCodomain, RNDM<N, F>>,
447 RNDM<F, N> : SpikeMerging<F>, 423 {
448 Reg : RegTermFW<F, A, ValueIteratorFactory<F, AlgIteratorOptions>, N> { 424 let opA = f.operator();
425 let b = f.data();
449 426
450 // Set up parameters 427 // Set up parameters
451 // We multiply tolerance by α for all algoritms. 428 // We multiply tolerance by α for all algoritms.
452 let tolerance = config.tolerance * reg.tolerance_scaling(); 429 let tolerance = config.tolerance * reg.tolerance_scaling();
453 let mut ε = tolerance.initial(); 430 let mut ε = tolerance.initial();
454 let findim_data = reg.prepare_optimise_weights(opA, b); 431 let findim_data = reg.prepare_optimise_weights(opA, b)?;
455 432
456 // Initialise operators 433 // Initialise operators
457 let preadjA = opA.preadjoint(); 434 let preadjA = opA.preadjoint();
458 435
459 // Initialise iterates 436 // Initialise iterates
460 let mut μ = DiscreteMeasure::new(); 437 let mut μ = μ0.unwrap_or_else(|| DiscreteMeasure::new());
461 let mut residual = -b; 438 let mut residual = f.residual(&μ);
462 439
463 // Statistics 440 // Statistics
464 let full_stats = |residual : &A::Observable, 441 let full_stats = |residual: &A::Observable, ν: &RNDM<N, F>, ε, stats| IterInfo {
465 ν : &RNDM<F, N>, 442 value: residual.norm2_squared_div2() + reg.apply(ν),
466 ε, stats| IterInfo { 443 n_spikes: ν.len(),
467 value : residual.norm2_squared_div2() + reg.apply(ν),
468 n_spikes : ν.len(),
469 ε, 444 ε,
470 .. stats 445 ..stats
471 }; 446 };
472 let mut stats = IterInfo::new(); 447 let mut stats = IterInfo::new();
473 448
474 // Run the algorithm 449 // Run the algorithm
475 for state in iterator.iter_init(|| full_stats(&residual, &μ, ε, stats.clone())) { 450 for state in iterator.iter_init(|| full_stats(&residual, &μ, ε, stats.clone())) {
478 453
479 // Calculate smooth part of surrogate model. 454 // Calculate smooth part of surrogate model.
480 let mut g = preadjA.apply(residual * (-1.0)); 455 let mut g = preadjA.apply(residual * (-1.0));
481 456
482 // Find absolute value maximising point 457 // Find absolute value maximising point
483 let (ξ, v_ξ) = reg.find_insertion(&mut g, refinement_tolerance, 458 let (ξ, v_ξ) =
484 config.refinement.max_steps); 459 reg.find_insertion(&mut g, refinement_tolerance, config.refinement.max_steps);
485 460
486 let inner_it = match config.variant { 461 let inner_it = match config.variant {
487 FWVariant::FullyCorrective => { 462 FWVariant::FullyCorrective => {
488 // No point in optimising the weight here: the finite-dimensional algorithm is fast. 463 // No point in optimising the weight here: the finite-dimensional algorithm is fast.
489 μ += DeltaMeasure { x : ξ, α : 0.0 }; 464 μ += DeltaMeasure { x: ξ, α: 0.0 };
490 stats.inserted += 1; 465 stats.inserted += 1;
491 config.inner.iterator_options.stop_target(inner_tolerance) 466 config.inner.iterator_options.stop_target(inner_tolerance)
492 }, 467 }
493 FWVariant::Relaxed => { 468 FWVariant::Relaxed => {
494 // Perform a relaxed initialisation of μ 469 // Perform a relaxed initialisation of μ
495 reg.relaxed_insert(&mut μ, &g, opA, ξ, v_ξ, &findim_data); 470 reg.relaxed_insert(&mut μ, &g, opA, ξ, v_ξ, &findim_data);
496 stats.inserted += 1; 471 stats.inserted += 1;
497 // The stop_target is only needed for the type system. 472 // The stop_target is only needed for the type system.
498 AlgIteratorOptions{ max_iter : 1, .. config.inner.iterator_options}.stop_target(0.0) 473 AlgIteratorOptions { max_iter: 1, ..config.inner.iterator_options }.stop_target(0.0)
499 } 474 }
500 }; 475 };
501 476
502 stats.inner_iters += reg.optimise_weights(&mut μ, opA, b, &findim_data, 477 stats.inner_iters +=
503 &config.inner, inner_it); 478 reg.optimise_weights(&mut μ, opA, b, &findim_data, &config.inner, inner_it);
504 479
505 // Merge spikes and update residual for next step and `if_verbose` below. 480 // Merge spikes and update residual for next step and `if_verbose` below.
506 let (r, count) = μ.merge_spikes_fitness(config.merging, 481 let (r, count) = μ.merge_spikes_fitness(
507 |μ̃| opA.apply(μ̃) - b, 482 config.merging,
508 A::Observable::norm2_squared); 483 |μ̃| f.residual(μ̃),
484 A::Observable::norm2_squared,
485 );
509 residual = r; 486 residual = r;
510 stats.merged += count; 487 stats.merged += count;
511 488
512 // Prune points with zero mass 489 // Prune points with zero mass
513 let n_before_prune = μ.len(); 490 let n_before_prune = μ.len();
518 stats.this_iters += 1; 495 stats.this_iters += 1;
519 let iter = state.iteration(); 496 let iter = state.iteration();
520 497
521 // Give statistics if needed 498 // Give statistics if needed
522 state.if_verbose(|| { 499 state.if_verbose(|| {
523 plotter.plot_spikes(iter, Some(&g), Option::<&S>::None, &μ); 500 plotter.plot_spikes(iter, Some(&g), Option::<&A::PreadjointCodomain>::None, &μ);
524 full_stats(&residual, &μ, ε, std::mem::replace(&mut stats, IterInfo::new())) 501 full_stats(
502 &residual,
503 &μ,
504 ε,
505 std::mem::replace(&mut stats, IterInfo::new()),
506 )
525 }); 507 });
526 508
527 // Update tolerance 509 // Update tolerance
528 ε = tolerance.update(ε, iter); 510 ε = tolerance.update(ε, iter);
529 } 511 }
530 512
531 // Return final iterate 513 // Return final iterate
532 μ 514 Ok(μ)
533 } 515 }

Mercurial Repository: pointsource_algs

mercurial