src/fb.rs

Sun, 11 Dec 2022 23:19:17 +0200

author
Tuomo Valkonen <tuomov@iki.fi>
date
Sun, 11 Dec 2022 23:19:17 +0200
changeset 23
9869fa1e0ccd
parent 13
bdc57366d4f5
child 24
d29d1fcf5423
permissions
-rw-r--r--

Print out experiment information when running it

0
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
1 /*!
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
2 Solver for the point source localisation problem using a forward-backward splitting method.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
3
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
4 This corresponds to the manuscript
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
5
13
bdc57366d4f5 arXiv links, README beautification
Tuomo Valkonen <tuomov@iki.fi>
parents: 8
diff changeset
6 * Valkonen T. - _Proximal methods for point source localisation_,
bdc57366d4f5 arXiv links, README beautification
Tuomo Valkonen <tuomov@iki.fi>
parents: 8
diff changeset
7 [arXiv:2212.02991](https://arxiv.org/abs/2212.02991).
0
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
8
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
9 The main routine is [`pointsource_fb`]. It is based on [`generic_pointsource_fb`], which is also
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
10 used by our [primal-dual proximal splitting][crate::pdps] implementation.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
11
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
12 FISTA-type inertia can also be enabled through [`FBConfig::meta`].
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
13
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
14 ## Problem
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
15
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
16 <p>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
17 Our objective is to solve
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
18 $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
19 \min_{μ ∈ ℳ(Ω)}~ F_0(Aμ-b) + α \|μ\|_{ℳ(Ω)} + δ_{≥ 0}(μ),
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
20 $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
21 where $F_0(y)=\frac{1}{2}\|y\|_2^2$ and the forward operator $A \in 𝕃(ℳ(Ω); ℝ^n)$.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
22 </p>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
23
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
24 ## Approach
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
25
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
26 <p>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
27 As documented in more detail in the paper, on each step we approximately solve
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
28 $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
29 \min_{μ ∈ ℳ(Ω)}~ F(x) + α \|μ\|_{ℳ(Ω)} + δ_{≥ 0}(x) + \frac{1}{2}\|μ-μ^k|_𝒟^2,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
30 $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
31 where $𝒟: 𝕃(ℳ(Ω); C_c(Ω))$ is typically a convolution operator.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
32 </p>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
33
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
34 ## Finite-dimensional subproblems.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
35
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
36 With $C$ a projection from [`DiscreteMeasure`] to the weights, and $x^k$ such that $x^k=Cμ^k$, we
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
37 form the discretised linearised inner problem
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
38 <p>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
39 $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
40 \min_{x ∈ ℝ^n}~ τ\bigl(F(Cx^k) + [C^*∇F(Cx^k)]^⊤(x-x^k) + α {\vec 1}^⊤ x\bigr)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
41 + δ_{≥ 0}(x) + \frac{1}{2}\|x-x^k\|_{C^*𝒟C}^2,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
42 $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
43 equivalently
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
44 $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
45 \begin{aligned}
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
46 \min_x~ & τF(Cx^k) - τ[C^*∇F(Cx^k)]^⊤x^k + \frac{1}{2} (x^k)^⊤ C^*𝒟C x^k
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
47 \\
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
48 &
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
49 - [C^*𝒟C x^k - τC^*∇F(Cx^k)]^⊤ x
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
50 \\
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
51 &
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
52 + \frac{1}{2} x^⊤ C^*𝒟C x
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
53 + τα {\vec 1}^⊤ x + δ_{≥ 0}(x),
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
54 \end{aligned}
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
55 $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
56 In other words, we obtain the quadratic non-negativity constrained problem
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
57 $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
58 \min_{x ∈ ℝ^n}~ \frac{1}{2} x^⊤ Ã x - b̃^⊤ x + c + τα {\vec 1}^⊤ x + δ_{≥ 0}(x).
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
59 $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
60 where
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
61 $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
62 \begin{aligned}
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
63 Ã & = C^*𝒟C,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
64 \\
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
65 g̃ & = C^*𝒟C x^k - τ C^*∇F(Cx^k)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
66 = C^* 𝒟 μ^k - τ C^*A^*(Aμ^k - b)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
67 \\
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
68 c & = τ F(Cx^k) - τ[C^*∇F(Cx^k)]^⊤x^k + \frac{1}{2} (x^k)^⊤ C^*𝒟C x^k
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
69 \\
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
70 &
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
71 = \frac{τ}{2} \|Aμ^k-b\|^2 - τ[Aμ^k-b]^⊤Aμ^k + \frac{1}{2} \|μ_k\|_{𝒟}^2
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
72 \\
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
73 &
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
74 = -\frac{τ}{2} \|Aμ^k-b\|^2 + τ[Aμ^k-b]^⊤ b + \frac{1}{2} \|μ_k\|_{𝒟}^2.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
75 \end{aligned}
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
76 $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
77 </p>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
78
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
79 We solve this with either SSN or FB via [`quadratic_nonneg`] as determined by
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
80 [`InnerSettings`] in [`FBGenericConfig::inner`].
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
81 */
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
82
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
83 use numeric_literals::replace_float_literals;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
84 use serde::{Serialize, Deserialize};
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
85 use colored::Colorize;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
86 use nalgebra::DVector;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
87
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
88 use alg_tools::iterate::{
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
89 AlgIteratorFactory,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
90 AlgIteratorState,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
91 };
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
92 use alg_tools::euclidean::Euclidean;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
93 use alg_tools::norms::Norm;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
94 use alg_tools::linops::Apply;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
95 use alg_tools::sets::Cube;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
96 use alg_tools::loc::Loc;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
97 use alg_tools::bisection_tree::{
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
98 BTFN,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
99 PreBTFN,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
100 Bounds,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
101 BTNodeLookup,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
102 BTNode,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
103 BTSearch,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
104 P2Minimise,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
105 SupportGenerator,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
106 LocalAnalysis,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
107 Bounded,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
108 };
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
109 use alg_tools::mapping::RealMapping;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
110 use alg_tools::nalgebra_support::ToNalgebraRealField;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
111
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
112 use crate::types::*;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
113 use crate::measures::{
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
114 DiscreteMeasure,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
115 DeltaMeasure,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
116 Radon
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
117 };
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
118 use crate::measures::merging::{
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
119 SpikeMergingMethod,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
120 SpikeMerging,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
121 };
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
122 use crate::forward_model::ForwardModel;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
123 use crate::seminorms::{
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
124 DiscreteMeasureOp, Lipschitz
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
125 };
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
126 use crate::subproblem::{
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
127 quadratic_nonneg,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
128 InnerSettings,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
129 InnerMethod,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
130 };
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
131 use crate::tolerance::Tolerance;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
132 use crate::plot::{
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
133 SeqPlotter,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
134 Plotting,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
135 PlotLookup
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
136 };
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
137
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
138 /// Method for constructing $μ$ on each iteration
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
139 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
140 #[allow(dead_code)]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
141 pub enum InsertionStyle {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
142 /// Resuse previous $μ$ from previous iteration, optimising weights
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
143 /// before inserting new spikes.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
144 Reuse,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
145 /// Start each iteration with $μ=0$.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
146 Zero,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
147 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
148
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
149 /// Meta-algorithm type
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
150 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
151 #[allow(dead_code)]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
152 pub enum FBMetaAlgorithm {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
153 /// No meta-algorithm
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
154 None,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
155 /// FISTA-style inertia
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
156 InertiaFISTA,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
157 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
158
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
159 /// Settings for [`pointsource_fb`].
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
160 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
161 #[serde(default)]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
162 pub struct FBConfig<F : Float> {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
163 /// Step length scaling
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
164 pub τ0 : F,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
165 /// Meta-algorithm to apply
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
166 pub meta : FBMetaAlgorithm,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
167 /// Generic parameters
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
168 pub insertion : FBGenericConfig<F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
169 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
170
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
171 /// Settings for the solution of the stepwise optimality condition in algorithms based on
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
172 /// [`generic_pointsource_fb`].
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
173 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
174 #[serde(default)]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
175 pub struct FBGenericConfig<F : Float> {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
176 /// Method for constructing $μ$ on each iteration; see [`InsertionStyle`].
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
177 pub insertion_style : InsertionStyle,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
178 /// Tolerance for point insertion.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
179 pub tolerance : Tolerance<F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
180 /// Stop looking for predual maximum (where to isert a new point) below
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
181 /// `tolerance` multiplied by this factor.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
182 pub insertion_cutoff_factor : F,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
183 /// Settings for branch and bound refinement when looking for predual maxima
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
184 pub refinement : RefinementSettings<F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
185 /// Maximum insertions within each outer iteration
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
186 pub max_insertions : usize,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
187 /// Pair `(n, m)` for maximum insertions `m` on first `n` iterations.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
188 pub bootstrap_insertions : Option<(usize, usize)>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
189 /// Inner method settings
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
190 pub inner : InnerSettings<F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
191 /// Spike merging method
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
192 pub merging : SpikeMergingMethod<F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
193 /// Tolerance multiplier for merges
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
194 pub merge_tolerance_mult : F,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
195 /// Spike merging method after the last step
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
196 pub final_merging : SpikeMergingMethod<F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
197 /// Iterations between merging heuristic tries
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
198 pub merge_every : usize,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
199 /// Save $μ$ for postprocessing optimisation
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
200 pub postprocessing : bool
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
201 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
202
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
203 #[replace_float_literals(F::cast_from(literal))]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
204 impl<F : Float> Default for FBConfig<F> {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
205 fn default() -> Self {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
206 FBConfig {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
207 τ0 : 0.99,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
208 meta : FBMetaAlgorithm::None,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
209 insertion : Default::default()
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
210 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
211 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
212 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
213
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
214 #[replace_float_literals(F::cast_from(literal))]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
215 impl<F : Float> Default for FBGenericConfig<F> {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
216 fn default() -> Self {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
217 FBGenericConfig {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
218 insertion_style : InsertionStyle::Reuse,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
219 tolerance : Default::default(),
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
220 insertion_cutoff_factor : 1.0,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
221 refinement : Default::default(),
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
222 max_insertions : 100,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
223 //bootstrap_insertions : None,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
224 bootstrap_insertions : Some((10, 1)),
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
225 inner : InnerSettings {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
226 method : InnerMethod::SSN,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
227 .. Default::default()
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
228 },
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
229 merging : SpikeMergingMethod::None,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
230 //merging : Default::default(),
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
231 final_merging : Default::default(),
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
232 merge_every : 10,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
233 merge_tolerance_mult : 2.0,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
234 postprocessing : false,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
235 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
236 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
237 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
238
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
239 /// Trait for specialisation of [`generic_pointsource_fb`] to basic FB, FISTA.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
240 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
241 /// The idea is that the residual $Aμ - b$ in the forward step can be replaced by an arbitrary
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
242 /// value. For example, to implement [primal-dual proximal splitting][crate::pdps] we replace it
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
243 /// with the dual variable $y$. We can then also implement alternative data terms, as the
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
244 /// (pre)differential of $F(μ)=F\_0(Aμ-b)$ is $F\'(μ) = A\_*F\_0\'(Aμ-b)$. In the case of the
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
245 /// quadratic fidelity $F_0(y)=\frac{1}{2}\\|y\\|_2^2$ in a Hilbert space, of course,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
246 /// $F\_0\'(Aμ-b)=Aμ-b$ is the residual.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
247 pub trait FBSpecialisation<F : Float, Observable : Euclidean<F>, const N : usize> : Sized {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
248 /// Updates the residual and does any necessary pruning of `μ`.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
249 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
250 /// Returns the new residual and possibly a new step length.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
251 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
252 /// The measure `μ` may also be modified to apply, e.g., inertia to it.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
253 /// The updated residual should correspond to the residual at `μ`.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
254 /// See the [trait documentation][FBSpecialisation] for the use and meaning of the residual.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
255 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
256 /// The parameter `μ_base` is the base point of the iteration, typically the previous iterate,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
257 /// but for, e.g., FISTA has inertia applied to it.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
258 fn update(
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
259 &mut self,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
260 μ : &mut DiscreteMeasure<Loc<F, N>, F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
261 μ_base : &DiscreteMeasure<Loc<F, N>, F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
262 ) -> (Observable, Option<F>);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
263
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
264 /// Calculates the data term value corresponding to iterate `μ` and available residual.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
265 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
266 /// Inertia and other modifications, as deemed, necessary, should be applied to `μ`.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
267 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
268 /// The blanket implementation correspondsn to the 2-norm-squared data fidelity
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
269 /// $\\|\text{residual}\\|\_2^2/2$.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
270 fn calculate_fit(
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
271 &self,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
272 _μ : &DiscreteMeasure<Loc<F, N>, F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
273 residual : &Observable
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
274 ) -> F {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
275 residual.norm2_squared_div2()
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
276 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
277
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
278 /// Calculates the data term value at $μ$.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
279 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
280 /// Unlike [`Self::calculate_fit`], no inertia, etc., should be applied to `μ`.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
281 fn calculate_fit_simple(
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
282 &self,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
283 μ : &DiscreteMeasure<Loc<F, N>, F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
284 ) -> F;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
285
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
286 /// Returns the final iterate after any necessary postprocess pruning, merging, etc.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
287 fn postprocess(self, mut μ : DiscreteMeasure<Loc<F, N>, F>, merging : SpikeMergingMethod<F>)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
288 -> DiscreteMeasure<Loc<F, N>, F>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
289 where DiscreteMeasure<Loc<F, N>, F> : SpikeMerging<F> {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
290 μ.merge_spikes_fitness(merging,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
291 |μ̃| self.calculate_fit_simple(μ̃),
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
292 |&v| v);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
293 μ.prune();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
294 μ
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
295 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
296
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
297 /// Returns measure to be used for value calculations, which may differ from μ.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
298 fn value_μ<'c, 'b : 'c>(&'b self, μ : &'c DiscreteMeasure<Loc<F, N>, F>)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
299 -> &'c DiscreteMeasure<Loc<F, N>, F> {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
300 μ
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
301 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
302 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
303
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
304 /// Specialisation of [`generic_pointsource_fb`] to basic μFB.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
305 struct BasicFB<
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
306 'a,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
307 F : Float + ToNalgebraRealField,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
308 A : ForwardModel<Loc<F, N>, F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
309 const N : usize
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
310 > {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
311 /// The data
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
312 b : &'a A::Observable,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
313 /// The forward operator
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
314 opA : &'a A,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
315 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
316
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
317 /// Implementation of [`FBSpecialisation`] for basic μFB forward-backward splitting.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
318 #[replace_float_literals(F::cast_from(literal))]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
319 impl<'a, F : Float + ToNalgebraRealField , A : ForwardModel<Loc<F, N>, F>, const N : usize>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
320 FBSpecialisation<F, A::Observable, N> for BasicFB<'a, F, A, N> {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
321 fn update(
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
322 &mut self,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
323 μ : &mut DiscreteMeasure<Loc<F, N>, F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
324 _μ_base : &DiscreteMeasure<Loc<F, N>, F>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
325 ) -> (A::Observable, Option<F>) {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
326 μ.prune();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
327 //*residual = self.opA.apply(μ) - self.b;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
328 let mut residual = self.b.clone();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
329 self.opA.gemv(&mut residual, 1.0, μ, -1.0);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
330 (residual, None)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
331 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
332
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
333 fn calculate_fit_simple(
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
334 &self,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
335 μ : &DiscreteMeasure<Loc<F, N>, F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
336 ) -> F {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
337 let mut residual = self.b.clone();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
338 self.opA.gemv(&mut residual, 1.0, μ, -1.0);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
339 residual.norm2_squared_div2()
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
340 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
341 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
342
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
343 /// Specialisation of [`generic_pointsource_fb`] to FISTA.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
344 struct FISTA<
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
345 'a,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
346 F : Float + ToNalgebraRealField,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
347 A : ForwardModel<Loc<F, N>, F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
348 const N : usize
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
349 > {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
350 /// The data
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
351 b : &'a A::Observable,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
352 /// The forward operator
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
353 opA : &'a A,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
354 /// Current inertial parameter
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
355 λ : F,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
356 /// Previous iterate without inertia applied.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
357 /// We need to store this here because `μ_base` passed to [`FBSpecialisation::update`] will
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
358 /// have inertia applied to it, so is not useful to use.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
359 μ_prev : DiscreteMeasure<Loc<F, N>, F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
360 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
361
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
362 /// Implementation of [`FBSpecialisation`] for μFISTA inertial forward-backward splitting.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
363 #[replace_float_literals(F::cast_from(literal))]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
364 impl<'a, F : Float + ToNalgebraRealField, A : ForwardModel<Loc<F, N>, F>, const N : usize>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
365 FBSpecialisation<F, A::Observable, N> for FISTA<'a, F, A, N> {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
366 fn update(
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
367 &mut self,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
368 μ : &mut DiscreteMeasure<Loc<F, N>, F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
369 _μ_base : &DiscreteMeasure<Loc<F, N>, F>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
370 ) -> (A::Observable, Option<F>) {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
371 // Update inertial parameters
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
372 let λ_prev = self.λ;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
373 self.λ = 2.0 * λ_prev / ( λ_prev + (4.0 + λ_prev * λ_prev).sqrt() );
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
374 let θ = self.λ / λ_prev - self.λ;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
375 // Perform inertial update on μ.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
376 // This computes μ ← (1 + θ) * μ - θ * μ_prev, pruning spikes where both μ
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
377 // and μ_prev have zero weight. Since both have weights from the finite-dimensional
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
378 // subproblem with a proximal projection step, this is likely to happen when the
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
379 // spike is not needed. A copy of the pruned μ without artithmetic performed is
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
380 // stored in μ_prev.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
381 μ.pruning_sub(1.0 + θ, θ, &mut self.μ_prev);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
382
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
383 //*residual = self.opA.apply(μ) - self.b;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
384 let mut residual = self.b.clone();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
385 self.opA.gemv(&mut residual, 1.0, μ, -1.0);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
386 (residual, None)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
387 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
388
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
389 fn calculate_fit_simple(
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
390 &self,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
391 μ : &DiscreteMeasure<Loc<F, N>, F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
392 ) -> F {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
393 let mut residual = self.b.clone();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
394 self.opA.gemv(&mut residual, 1.0, μ, -1.0);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
395 residual.norm2_squared_div2()
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
396 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
397
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
398 fn calculate_fit(
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
399 &self,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
400 _μ : &DiscreteMeasure<Loc<F, N>, F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
401 _residual : &A::Observable
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
402 ) -> F {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
403 self.calculate_fit_simple(&self.μ_prev)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
404 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
405
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
406 // For FISTA we need to do a final pruning as well, due to the limited
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
407 // pruning that can be done on each step.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
408 fn postprocess(mut self, μ_base : DiscreteMeasure<Loc<F, N>, F>, merging : SpikeMergingMethod<F>)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
409 -> DiscreteMeasure<Loc<F, N>, F>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
410 where DiscreteMeasure<Loc<F, N>, F> : SpikeMerging<F> {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
411 let mut μ = self.μ_prev;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
412 self.μ_prev = μ_base;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
413 μ.merge_spikes_fitness(merging,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
414 |μ̃| self.calculate_fit_simple(μ̃),
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
415 |&v| v);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
416 μ.prune();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
417 μ
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
418 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
419
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
420 fn value_μ<'c, 'b : 'c>(&'c self, _μ : &'c DiscreteMeasure<Loc<F, N>, F>)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
421 -> &'c DiscreteMeasure<Loc<F, N>, F> {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
422 &self.μ_prev
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
423 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
424 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
425
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
426 /// Iteratively solve the pointsource localisation problem using forward-backward splitting
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
427 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
428 /// The settings in `config` have their [respective documentation](FBConfig). `opA` is the
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
429 /// forward operator $A$, $b$ the observable, and $\lambda$ the regularisation weight.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
430 /// The operator `op𝒟` is used for forming the proximal term. Typically it is a convolution
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
431 /// operator. Finally, the `iterator` is an outer loop verbosity and iteration count control
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
432 /// as documented in [`alg_tools::iterate`].
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
433 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
434 /// For details on the mathematical formulation, see the [module level](self) documentation.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
435 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
436 /// Returns the final iterate.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
437 #[replace_float_literals(F::cast_from(literal))]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
438 pub fn pointsource_fb<'a, F, I, A, GA, 𝒟, BTA, G𝒟, S, K, const N : usize>(
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
439 opA : &'a A,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
440 b : &A::Observable,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
441 α : F,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
442 op𝒟 : &'a 𝒟,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
443 config : &FBConfig<F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
444 iterator : I,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
445 plotter : SeqPlotter<F, N>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
446 ) -> DiscreteMeasure<Loc<F, N>, F>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
447 where F : Float + ToNalgebraRealField,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
448 I : AlgIteratorFactory<IterInfo<F, N>>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
449 for<'b> &'b A::Observable : std::ops::Neg<Output=A::Observable>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
450 //+ std::ops::Mul<F, Output=A::Observable>, <-- FIXME: compiler overflow
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
451 A::Observable : std::ops::MulAssign<F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
452 GA : SupportGenerator<F, N, SupportType = S, Id = usize> + Clone,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
453 A : ForwardModel<Loc<F, N>, F, PreadjointCodomain = BTFN<F, GA, BTA, N>>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
454 + Lipschitz<𝒟, FloatType=F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
455 BTA : BTSearch<F, N, Data=usize, Agg=Bounds<F>>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
456 G𝒟 : SupportGenerator<F, N, SupportType = K, Id = usize> + Clone,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
457 𝒟 : DiscreteMeasureOp<Loc<F, N>, F, PreCodomain = PreBTFN<F, G𝒟, N>>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
458 𝒟::Codomain : RealMapping<F, N>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
459 S: RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
460 K: RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
461 BTNodeLookup: BTNode<F, usize, Bounds<F>, N>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
462 Cube<F, N>: P2Minimise<Loc<F, N>, F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
463 PlotLookup : Plotting<N>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
464 DiscreteMeasure<Loc<F, N>, F> : SpikeMerging<F> {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
465
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
466 let initial_residual = -b;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
467 let τ = config.τ0/opA.lipschitz_factor(&op𝒟).unwrap();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
468
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
469 match config.meta {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
470 FBMetaAlgorithm::None => generic_pointsource_fb(
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
471 opA, α, op𝒟, τ, &config.insertion, iterator, plotter, initial_residual,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
472 BasicFB{ b, opA }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
473 ),
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
474 FBMetaAlgorithm::InertiaFISTA => generic_pointsource_fb(
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
475 opA, α, op𝒟, τ, &config.insertion, iterator, plotter, initial_residual,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
476 FISTA{ b, opA, λ : 1.0, μ_prev : DiscreteMeasure::new() }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
477 ),
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
478 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
479 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
480
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
481 /// Generic implementation of [`pointsource_fb`].
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
482 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
483 /// The method can be specialised to even primal-dual proximal splitting through the
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
484 /// [`FBSpecialisation`] parameter `specialisation`.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
485 /// The settings in `config` have their [respective documentation](FBGenericConfig). `opA` is the
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
486 /// forward operator $A$, $b$ the observable, and $\lambda$ the regularisation weight.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
487 /// The operator `op𝒟` is used for forming the proximal term. Typically it is a convolution
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
488 /// operator. Finally, the `iterator` is an outer loop verbosity and iteration count control
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
489 /// as documented in [`alg_tools::iterate`].
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
490 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
491 /// The implementation relies on [`alg_tools::bisection_tree::BTFN`] presentations of
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
492 /// sums of simple functions usign bisection trees, and the related
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
493 /// [`alg_tools::bisection_tree::Aggregator`]s, to efficiently search for component functions
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
494 /// active at a specific points, and to maximise their sums. Through the implementation of the
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
495 /// [`alg_tools::bisection_tree::BT`] bisection trees, it also relies on the copy-on-write features
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
496 /// of [`std::sync::Arc`] to only update relevant parts of the bisection tree when adding functions.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
497 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
498 /// Returns the final iterate.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
499 #[replace_float_literals(F::cast_from(literal))]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
500 pub fn generic_pointsource_fb<'a, F, I, A, GA, 𝒟, BTA, G𝒟, S, K, Spec, const N : usize>(
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
501 opA : &'a A,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
502 α : F,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
503 op𝒟 : &'a 𝒟,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
504 mut τ : F,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
505 config : &FBGenericConfig<F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
506 iterator : I,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
507 mut plotter : SeqPlotter<F, N>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
508 mut residual : A::Observable,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
509 mut specialisation : Spec,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
510 ) -> DiscreteMeasure<Loc<F, N>, F>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
511 where F : Float + ToNalgebraRealField,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
512 I : AlgIteratorFactory<IterInfo<F, N>>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
513 Spec : FBSpecialisation<F, A::Observable, N>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
514 A::Observable : std::ops::MulAssign<F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
515 GA : SupportGenerator<F, N, SupportType = S, Id = usize> + Clone,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
516 A : ForwardModel<Loc<F, N>, F, PreadjointCodomain = BTFN<F, GA, BTA, N>>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
517 + Lipschitz<𝒟, FloatType=F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
518 BTA : BTSearch<F, N, Data=usize, Agg=Bounds<F>>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
519 G𝒟 : SupportGenerator<F, N, SupportType = K, Id = usize> + Clone,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
520 𝒟 : DiscreteMeasureOp<Loc<F, N>, F, PreCodomain = PreBTFN<F, G𝒟, N>>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
521 𝒟::Codomain : RealMapping<F, N>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
522 S: RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
523 K: RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
524 BTNodeLookup: BTNode<F, usize, Bounds<F>, N>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
525 Cube<F, N>: P2Minimise<Loc<F, N>, F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
526 PlotLookup : Plotting<N>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
527 DiscreteMeasure<Loc<F, N>, F> : SpikeMerging<F> {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
528
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
529 // Set up parameters
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
530 let quiet = iterator.is_quiet();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
531 let op𝒟norm = op𝒟.opnorm_bound();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
532 // We multiply tolerance by τ for FB since
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
533 // our subproblems depending on tolerances are scaled by τ compared to the conditional
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
534 // gradient approach.
8
ea3ca78873e8 Clean up / remove various unused FB algorithm family hacks.
Tuomo Valkonen <tuomov@iki.fi>
parents: 7
diff changeset
535 let tolerance = config.tolerance * τ * α;
0
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
536 let mut ε = tolerance.initial();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
537
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
538 // Initialise operators
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
539 let preadjA = opA.preadjoint();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
540
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
541 // Initialise iterates
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
542 let mut μ = DiscreteMeasure::new();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
543
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
544 let mut inner_iters = 0;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
545 let mut this_iters = 0;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
546 let mut pruned = 0;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
547 let mut merged = 0;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
548
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
549 let μ_diff = |μ_new : &DiscreteMeasure<Loc<F, N>, F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
550 μ_base : &DiscreteMeasure<Loc<F, N>, F>| {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
551 let mut ν : DiscreteMeasure<Loc<F, N>, F> = match config.insertion_style {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
552 InsertionStyle::Reuse => {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
553 μ_new.iter_spikes()
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
554 .zip(μ_base.iter_masses().chain(std::iter::repeat(0.0)))
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
555 .map(|(δ, α_base)| (δ.x, α_base - δ.α))
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
556 .collect()
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
557 },
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
558 InsertionStyle::Zero => {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
559 μ_new.iter_spikes()
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
560 .map(|δ| -δ)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
561 .chain(μ_base.iter_spikes().copied())
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
562 .collect()
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
563 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
564 };
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
565 ν.prune(); // Potential small performance improvement
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
566 ν
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
567 };
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
568
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
569 // Run the algorithm
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
570 iterator.iterate(|state| {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
571 // Calculate subproblem tolerances, and update main tolerance for next iteration
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
572 let τα = τ * α;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
573 let target_bounds = Bounds(τα - ε, τα + ε);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
574 let merge_tolerance = config.merge_tolerance_mult * ε;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
575 let merge_target_bounds = Bounds(τα - merge_tolerance, τα + merge_tolerance);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
576 let inner_tolerance = ε * config.inner.tolerance_mult;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
577 let refinement_tolerance = ε * config.refinement.tolerance_mult;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
578 let maximise_above = τα + ε * config.insertion_cutoff_factor;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
579 let ε_prev = ε;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
580 ε = tolerance.update(ε, state.iteration());
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
581
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
582 // Maximum insertion count and measure difference calculation depend on insertion style.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
583 let (m, warn_insertions) = match (state.iteration(), config.bootstrap_insertions) {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
584 (i, Some((l, k))) if i <= l => (k, false),
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
585 _ => (config.max_insertions, !quiet),
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
586 };
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
587 let max_insertions = match config.insertion_style {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
588 InsertionStyle::Zero => {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
589 todo!("InsertionStyle::Zero does not currently work with FISTA, so diabled.");
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
590 // let n = μ.len();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
591 // μ = DiscreteMeasure::new();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
592 // n + m
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
593 },
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
594 InsertionStyle::Reuse => m,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
595 };
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
596
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
597 // Calculate smooth part of surrogate model.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
598 // Using `std::mem::replace` here is not ideal, and expects that `empty_observable`
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
599 // has no significant overhead. For some reosn Rust doesn't allow us simply moving
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
600 // the residual and replacing it below before the end of this closure.
7
c32171f7cce5 Remove ergodic tolerance; it's not useful.
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
601 residual *= -τ;
0
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
602 let r = std::mem::replace(&mut residual, opA.empty_observable());
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
603 let minus_τv = preadjA.apply(r); // minus_τv = -τA^*(Aμ^k-b)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
604 // TODO: should avoid a second copy of μ here; μ_base already stores a copy.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
605 let ω0 = op𝒟.apply(μ.clone()); // 𝒟μ^k
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
606 //let g = &minus_τv + ω0; // Linear term of surrogate model
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
607
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
608 // Save current base point
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
609 let μ_base = μ.clone();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
610
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
611 // Add points to support until within error tolerance or maximum insertion count reached.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
612 let mut count = 0;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
613 let (within_tolerances, d) = 'insertion: loop {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
614 if μ.len() > 0 {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
615 // Form finite-dimensional subproblem. The subproblem references to the original μ^k
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
616 // from the beginning of the iteration are all contained in the immutable c and g.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
617 let à = op𝒟.findim_matrix(μ.iter_locations());
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
618 let g̃ = DVector::from_iterator(μ.len(),
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
619 μ.iter_locations()
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
620 .map(|ζ| minus_τv.apply(ζ) + ω0.apply(ζ))
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
621 .map(F::to_nalgebra_mixed));
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
622 let mut x = μ.masses_dvector();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
623
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
624 // The gradient of the forward component of the inner objective is C^*𝒟Cx - g̃.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
625 // We have |C^*𝒟Cx|_2 = sup_{|z|_2 ≤ 1} ⟨z, C^*𝒟Cx⟩ = sup_{|z|_2 ≤ 1} ⟨Cz|𝒟Cx⟩
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
626 // ≤ sup_{|z|_2 ≤ 1} |Cz|_ℳ |𝒟Cx|_∞ ≤ sup_{|z|_2 ≤ 1} |Cz|_ℳ |𝒟| |Cx|_ℳ
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
627 // ≤ sup_{|z|_2 ≤ 1} |z|_1 |𝒟| |x|_1 ≤ sup_{|z|_2 ≤ 1} n |z|_2 |𝒟| |x|_2
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
628 // = n |𝒟| |x|_2, where n is the number of points. Therefore
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
629 let inner_τ = config.inner.τ0 / (op𝒟norm * F::cast_from(μ.len()));
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
630
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
631 // Solve finite-dimensional subproblem.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
632 let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
633 inner_iters += quadratic_nonneg(config.inner.method, &Ã, &g̃, τ*α, &mut x,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
634 inner_τ, inner_it);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
635
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
636 // Update masses of μ based on solution of finite-dimensional subproblem.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
637 μ.set_masses_dvector(&x);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
638 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
639
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
640 // Form d = ω0 - τv - 𝒟μ = -𝒟(μ - μ^k) - τv for checking the proximate optimality
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
641 // conditions in the predual space, and finding new points for insertion, if necessary.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
642 let mut d = &minus_τv + op𝒟.preapply(μ_diff(&μ, &μ_base));
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
643
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
644 // If no merging heuristic is used, let's be more conservative about spike insertion,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
645 // and skip it after first round. If merging is done, being more greedy about spike
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
646 // insertion also seems to improve performance.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
647 let may_break = if let SpikeMergingMethod::None = config.merging {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
648 false
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
649 } else {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
650 count > 0
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
651 };
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
652
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
653 // If preliminary check indicates that we are in bonds, and if it otherwise matches
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
654 // the insertion strategy, skip insertion.
7
c32171f7cce5 Remove ergodic tolerance; it's not useful.
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
655 if may_break && target_bounds.superset(&d.bounds()) {
0
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
656 break 'insertion (true, d)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
657 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
658
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
659 // If the rough check didn't indicate stopping, find maximising point, maintaining for
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
660 // the calculations in the beginning of the loop that v_ξ = (ω0-τv-𝒟μ)(ξ) = d(ξ),
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
661 // where 𝒟μ is now distinct from μ0 after the insertions already performed.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
662 // We do not need to check lower bounds, as a solution of the finite-dimensional
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
663 // subproblem should always satisfy them.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
664
8
ea3ca78873e8 Clean up / remove various unused FB algorithm family hacks.
Tuomo Valkonen <tuomov@iki.fi>
parents: 7
diff changeset
665 // If μ has some spikes, only find a maximum of d if it is above a threshold
ea3ca78873e8 Clean up / remove various unused FB algorithm family hacks.
Tuomo Valkonen <tuomov@iki.fi>
parents: 7
diff changeset
666 // defined by the refinment tolerance.
ea3ca78873e8 Clean up / remove various unused FB algorithm family hacks.
Tuomo Valkonen <tuomov@iki.fi>
parents: 7
diff changeset
667 let (ξ, v_ξ) = match d.maximise_above(maximise_above, refinement_tolerance,
ea3ca78873e8 Clean up / remove various unused FB algorithm family hacks.
Tuomo Valkonen <tuomov@iki.fi>
parents: 7
diff changeset
668 config.refinement.max_steps) {
ea3ca78873e8 Clean up / remove various unused FB algorithm family hacks.
Tuomo Valkonen <tuomov@iki.fi>
parents: 7
diff changeset
669 None => break 'insertion (true, d),
ea3ca78873e8 Clean up / remove various unused FB algorithm family hacks.
Tuomo Valkonen <tuomov@iki.fi>
parents: 7
diff changeset
670 Some(res) => res,
0
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
671 };
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
672
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
673 // Break if maximum insertion count reached
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
674 if count >= max_insertions {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
675 let in_bounds2 = target_bounds.upper() >= v_ξ;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
676 break 'insertion (in_bounds2, d)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
677 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
678
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
679 // No point in optimising the weight here; the finite-dimensional algorithm is fast.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
680 μ += DeltaMeasure { x : ξ, α : 0.0 };
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
681 count += 1;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
682 };
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
683
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
684 if !within_tolerances && warn_insertions {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
685 // Complain (but continue) if we failed to get within tolerances
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
686 // by inserting more points.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
687 let err = format!("Maximum insertions reached without achieving \
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
688 subproblem solution tolerance");
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
689 println!("{}", err.red());
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
690 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
691
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
692 // Merge spikes
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
693 if state.iteration() % config.merge_every == 0 {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
694 let n_before_merge = μ.len();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
695 μ.merge_spikes(config.merging, |μ_candidate| {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
696 let mut d = &minus_τv + op𝒟.preapply(μ_diff(&μ_candidate, &μ_base));
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
697
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
698 if merge_target_bounds.superset(&d.bounds()) {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
699 return Some(())
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
700 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
701
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
702 let d_min_supp = μ_candidate.iter_spikes().filter_map(|&DeltaMeasure{ α, ref x }| {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
703 (α != 0.0).then(|| d.apply(x))
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
704 }).reduce(F::min);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
705
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
706 if d_min_supp.map_or(true, |b| b >= merge_target_bounds.lower()) &&
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
707 d.has_upper_bound(merge_target_bounds.upper(), refinement_tolerance,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
708 config.refinement.max_steps) {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
709 Some(())
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
710 } else {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
711 None
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
712 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
713 });
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
714 debug_assert!(μ.len() >= n_before_merge);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
715 merged += μ.len() - n_before_merge;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
716 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
717
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
718 let n_before_prune = μ.len();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
719 (residual, τ) = match specialisation.update(&mut μ, &μ_base) {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
720 (r, None) => (r, τ),
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
721 (r, Some(new_τ)) => (r, new_τ)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
722 };
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
723 debug_assert!(μ.len() <= n_before_prune);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
724 pruned += n_before_prune - μ.len();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
725
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
726 this_iters += 1;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
727
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
728 // Give function value if needed
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
729 state.if_verbose(|| {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
730 let value_μ = specialisation.value_μ(&μ);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
731 // Plot if so requested
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
732 plotter.plot_spikes(
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
733 format!("iter {} end; {}", state.iteration(), within_tolerances), &d,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
734 "start".to_string(), Some(&minus_τv),
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
735 Some(target_bounds), value_μ,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
736 );
8
ea3ca78873e8 Clean up / remove various unused FB algorithm family hacks.
Tuomo Valkonen <tuomov@iki.fi>
parents: 7
diff changeset
737 // Calculate mean inner iterations and reset relevant counters.
0
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
738 // Return the statistics
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
739 let res = IterInfo {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
740 value : specialisation.calculate_fit(&μ, &residual) + α * value_μ.norm(Radon),
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
741 n_spikes : value_μ.len(),
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
742 inner_iters,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
743 this_iters,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
744 merged,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
745 pruned,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
746 ε : ε_prev,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
747 postprocessing: config.postprocessing.then(|| value_μ.clone()),
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
748 };
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
749 inner_iters = 0;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
750 this_iters = 0;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
751 merged = 0;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
752 pruned = 0;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
753 res
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
754 })
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
755 });
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
756
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
757 specialisation.postprocess(μ, config.final_merging)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
758 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
759
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
760
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
761
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
762

mercurial