Sun, 25 Sep 2022 21:45:56 +0300
Attempt at directly referring to measure masses as SliceStorageMut
0 | 1 | //! Iterative algorithms for solving finite-dimensional subproblems. |
2 | ||
3 | use serde::{Serialize, Deserialize}; | |
4
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
4 | use nalgebra::{ |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
5 | Matrix, |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
6 | Vector, |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
7 | Storage, |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
8 | StorageMut, |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
9 | Const, |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
10 | Dim, |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
11 | DVector, |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
12 | DMatrix, |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
13 | DefaultAllocator, |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
14 | Dynamic, |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
15 | }; |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
16 | use nalgebra::base::allocator::Allocator; |
0 | 17 | use numeric_literals::replace_float_literals; |
18 | use itertools::{izip, Itertools}; | |
19 | use colored::Colorize; | |
20 | ||
21 | use alg_tools::iter::Mappable; | |
22 | use alg_tools::error::NumericalError; | |
23 | use alg_tools::iterate::{ | |
24 | AlgIteratorFactory, | |
25 | AlgIteratorState, | |
26 | AlgIteratorOptions, | |
27 | Verbose, | |
28 | Step, | |
29 | }; | |
30 | use alg_tools::linops::GEMV; | |
31 | use alg_tools::nalgebra_support::ToNalgebraRealField; | |
32 | ||
33 | use crate::types::*; | |
34 | ||
35 | /// Method for solving finite-dimensional subproblems | |
36 | #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)] | |
37 | #[allow(dead_code)] | |
38 | pub enum InnerMethod { | |
39 | /// Forward-backward | |
40 | FB, | |
41 | /// Semismooth Newton | |
42 | SSN, | |
43 | } | |
44 | ||
45 | /// Settings for the solution of finite-dimensional subproblems | |
46 | #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)] | |
47 | pub struct InnerSettings<F : Float> { | |
48 | /// Method | |
49 | pub method : InnerMethod, | |
50 | /// Proportional step length (∈ [0, 1) for `InnerMethod::FB`). | |
51 | pub τ0 : F, | |
52 | /// Fraction of `tolerance` given to inner algorithm | |
53 | pub tolerance_mult : F, | |
54 | /// Iterator options | |
55 | #[serde(flatten)] | |
56 | pub iterator_options : AlgIteratorOptions, | |
57 | } | |
58 | ||
59 | #[replace_float_literals(F::cast_from(literal))] | |
60 | impl<F : Float> Default for InnerSettings<F> { | |
61 | fn default() -> Self { | |
62 | InnerSettings { | |
63 | τ0 : 0.99, | |
64 | iterator_options : AlgIteratorOptions { | |
65 | // max_iter cannot be very small, as initially FB needs many iterations, although | |
66 | // on later invocations even one or two tends to be enough | |
67 | max_iter : 2000, | |
68 | // verbose_iter affects testing of sufficient convergence, so we set it to | |
69 | // a small value… | |
70 | verbose_iter : Verbose::Every(1), | |
71 | // … but don't print out anything | |
72 | quiet : true, | |
73 | .. Default::default() | |
74 | }, | |
75 | method : InnerMethod::FB, | |
76 | tolerance_mult : 0.01, | |
77 | } | |
78 | } | |
79 | } | |
80 | ||
81 | /// Compute the proximal operator of $x \mapsto x + \delta\_{[0, \infty)}$, i.e., | |
82 | /// the non-negativity contrained soft-thresholding operator. | |
83 | #[inline] | |
84 | #[replace_float_literals(F::cast_from(literal))] | |
85 | fn nonneg_soft_thresholding<F : Float>(v : F, λ : F) -> F { | |
86 | (v - λ).max(0.0) | |
87 | } | |
88 | ||
89 | /// Forward-backward splitting implementation of [`quadratic_nonneg`]. | |
90 | /// For detailed documentation of the inputs and outputs, refer to there. | |
91 | /// | |
92 | /// The `λ` component of the model is handled in the proximal step instead of the gradient step | |
93 | /// for potential performance improvements. | |
94 | #[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())] | |
4
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
95 | pub fn quadratic_nonneg_fb<F, I, S, SA, SG, D>( |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
96 | mA : &Matrix<F::MixedType, D, D, SA>, |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
97 | g : &Vector<F::MixedType, D, SG>, |
0 | 98 | //c_ : F, |
99 | λ_ : F, | |
4
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
100 | x : &mut Vector<F::MixedType, D, S>, |
0 | 101 | τ_ : F, |
102 | iterator : I | |
103 | ) -> usize | |
104 | where F : Float + ToNalgebraRealField, | |
4
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
105 | I : AlgIteratorFactory<F>, |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
106 | D : Dim, |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
107 | S : StorageMut<F::MixedType, D, Const::<1>>, |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
108 | SA : Storage<F::MixedType, D, D>, |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
109 | SG : Storage<F::MixedType, D, Const::<1>>, |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
110 | DefaultAllocator : Allocator<F::MixedType, D> |
0 | 111 | { |
4
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
112 | let mut xprev = x.clone_owned(); |
0 | 113 | //let c = c_.to_nalgebra_mixed(); |
114 | let λ = λ_.to_nalgebra_mixed(); | |
115 | let τ = τ_.to_nalgebra_mixed(); | |
116 | let τλ = τ * λ; | |
117 | let mut iters = 0; | |
4
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
118 | let mut v = { |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
119 | let (r, c) = x.shape_generic(); |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
120 | Vector::zeros_generic(r, c) |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
121 | }; |
0 | 122 | |
123 | iterator.iterate(|state| { | |
124 | // Replace `x` with $x - τ[Ax-g]= [x + τg]- τAx$ | |
125 | v.copy_from(g); // v = g | |
126 | v.axpy(1.0, x, τ); // v = x + τ*g | |
127 | v.sygemv(-τ, mA, x, 1.0); // v = [x + τg]- τAx | |
128 | let backup = state.if_verbose(|| { | |
129 | xprev.copy_from(x) | |
130 | }); | |
131 | // Calculate the proximal map | |
132 | x.iter_mut().zip(v.iter()).for_each(|(x_i, &v_i)| { | |
133 | *x_i = nonneg_soft_thresholding(v_i, τλ); | |
134 | }); | |
135 | ||
136 | iters +=1; | |
137 | ||
138 | backup.map(|_| { | |
139 | // The subdifferential of the objective is $Ax - g + λ + ∂ δ_{≥ 0}(x)$. | |
140 | // We return the minimal ∞-norm over all subderivatives. | |
141 | v.copy_from(g); // d = g | |
142 | mA.gemv(&mut v, 1.0, x, -1.0); // d = Ax - g | |
143 | let mut val = 0.0; | |
144 | for (&v_i, &x_i) in izip!(v.iter(), x.iter()) { | |
145 | let d = v_i + λ; | |
146 | if x_i > 0.0 || d < 0.0 { | |
147 | val = val.max(d.abs()); | |
148 | } | |
149 | } | |
150 | F::from_nalgebra_mixed(val) | |
151 | }) | |
152 | }); | |
153 | ||
154 | iters | |
155 | } | |
156 | ||
157 | /// Semismooth Newton implementation of [`quadratic_nonneg`]. | |
158 | /// | |
159 | /// For detailed documentation of the inputs, refer to there. | |
160 | /// This function returns the number of iterations taken if there was no inversion failure, | |
161 | /// | |
162 | /// ## Method derivation | |
163 | /// | |
164 | /// **The below may look like garbage. Sorry, but rustdoc is obsolete rubbish | |
165 | /// that doesn't directly support by-now standard-in-markdown LaTeX math. Instead it | |
166 | /// forces one into unreliable KaTeX autorender postprocessing andescape hell and that | |
167 | /// it doesn't even process correctly.** | |
168 | /// | |
169 | /// <p> | |
170 | /// For the objective | |
171 | /// $$ | |
172 | /// J(x) = \frac{1}{2} x^⊤Ax - g^⊤ x + λ{\vec 1}^⊤ x + c + δ_{≥ 0}(x), | |
173 | /// $$ | |
174 | /// we have the optimality condition | |
175 | /// $$ | |
176 | /// x - \mathop{\mathrm{prox}}_{τλ{\vec 1}^⊤ + δ_{≥ 0}}(x - τ[Ax-g^⊤]) = 0, | |
177 | /// $$ | |
178 | /// which we write as | |
179 | /// $$ | |
180 | /// x - [G ∘ F](x)=0 | |
181 | /// $$ | |
182 | /// for | |
183 | /// $$ | |
184 | /// G(x) = \mathop{\mathrm{prox}}_{λ{\vec 1}^⊤ + δ_{≥ 0}} | |
185 | /// \quad\text{and}\quad | |
186 | /// F(x) = x - τ Ax + τ g^⊤ | |
187 | /// $$ | |
188 | /// We can use Newton derivative chain rule to compute | |
189 | /// $D_N[G ∘ F](x) = D_N G(F(x)) D_N F(x)$, where | |
190 | /// $D_N F(x) = \mathop{\mathrm{Id}} - τ A$, | |
191 | /// and $[D_N G(F(x))]_i = 1$ for inactive coordinates and $=0$ for active coordinates. | |
192 | /// </p> | |
193 | /// | |
194 | /// <p> | |
195 | /// The method itself involves solving $D_N[Id - G ∘ F](x^k) s^k = - [Id - G ∘ F](x^k)$ and | |
196 | /// updating $x^{k+1} = x^k + s^k$. Consequently | |
197 | /// $$ | |
198 | /// s^k - D_N G(F(x^k)) [s^k - τ As^k] = - x^k + [G ∘ F](x^k) | |
199 | /// $$ | |
200 | /// For $𝒜$ the set of active coordinates and $ℐ$ the set of inactive coordinates, this | |
201 | /// expands as | |
202 | /// $$ | |
203 | /// [τ A_{ℐ × ℐ}]s^k_ℐ = - x^k_ℐ + [G ∘ F](x^k)_ℐ - [τ A_{ℐ × 𝒜}]s^k_𝒜 | |
204 | /// $$ | |
205 | /// and | |
206 | /// $$ | |
207 | /// s^k_𝒜 = - x^k_𝒜 + [G ∘ F](x^k)_𝒜. | |
208 | /// $$ | |
209 | /// Thus on $𝒜$ the update $[x^k + s^k]_𝒜 = [G ∘ F](x^k)_𝒜$ is just the forward-backward update. | |
210 | /// </p> | |
211 | /// | |
212 | /// <p> | |
213 | /// We need to detect stopping by a subdifferential and return $x$ satisfying $x ≥ 0$, | |
214 | /// which is in general not true for the SSN. We therefore use that $[G ∘ F](x^k)$ is a valid | |
215 | /// forward-backward step. | |
216 | /// </p> | |
217 | #[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())] | |
4
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
218 | pub fn quadratic_nonneg_ssn<F, I, S, SA, SG>( |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
219 | mA : &Matrix<F::MixedType, Dynamic, Dynamic, SA>, |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
220 | g : &Vector<F::MixedType, Dynamic, SG>, |
0 | 221 | //c_ : F, |
222 | λ_ : F, | |
4
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
223 | x : &mut Vector<F::MixedType, Dynamic, S>, |
0 | 224 | τ_ : F, |
225 | iterator : I | |
226 | ) -> Result<usize, NumericalError> | |
4
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
227 | where F : Float + ToNalgebraRealField<MixedType=F>, |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
228 | I : AlgIteratorFactory<F>, |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
229 | S : StorageMut<F::MixedType, Dynamic, Const::<1>>, |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
230 | SA : Storage<F::MixedType, Dynamic, Dynamic>, |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
231 | SG : Storage<F::MixedType, Dynamic, Const::<1>>, |
0 | 232 | { |
233 | let n = x.len(); | |
4
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
234 | let mut xprev = x.clone_owned(); |
0 | 235 | let mut v = DVector::zeros(n); |
236 | //let c = c_.to_nalgebra_mixed(); | |
237 | let λ = λ_.to_nalgebra_mixed(); | |
238 | let τ = τ_.to_nalgebra_mixed(); | |
239 | let τλ = τ * λ; | |
240 | let mut inact : Vec<bool> = Vec::from_iter(std::iter::repeat(false).take(n)); | |
4
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
241 | let mut v = { |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
242 | let (r, c) = x.shape_generic(); |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
243 | Vector::zeros_generic(r, c) |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
244 | }; |
0 | 245 | let mut s = DVector::zeros(0); |
246 | let mut decomp = nalgebra::linalg::LU::new(DMatrix::zeros(0, 0)); | |
247 | let mut iters = 0; | |
248 | ||
249 | let res = iterator.iterate_fallible(|state| { | |
250 | // 1. Perform delayed SSN-update based on previously computed step on active | |
251 | // coordinates. The step is delayed to the beginning of the loop because | |
252 | // the SSN step may violate constraints, so we arrange `x` to contain at the | |
253 | // end of the loop the valid FB step that forms part of the SSN step | |
254 | let mut si = s.iter(); | |
255 | for (&ast, x_i, xprev_i) in izip!(inact.iter(), x.iter_mut(), xprev.iter_mut()) { | |
256 | if ast { | |
257 | *x_i = *xprev_i + *si.next().unwrap() | |
258 | } | |
259 | *xprev_i = *x_i; | |
260 | } | |
261 | ||
262 | //xprev.copy_from(x); | |
263 | ||
264 | // 2. Calculate FB step. | |
265 | // 2.1. Replace `x` with $x⁻ - τ[Ax⁻-g]= [x⁻ + τg]- τAx⁻$ | |
266 | x.axpy(τ, g, 1.0); // x = x⁻ + τ*g | |
267 | x.sygemv(-τ, mA, &xprev, 1.0); // x = [x⁻ + τg]- τAx⁻ | |
268 | // 2.2. Calculate prox and set of active coordinates at the same time | |
269 | let mut act_changed = false; | |
270 | let mut n_inact = 0; | |
271 | for (x_i, ast) in izip!(x.iter_mut(), inact.iter_mut()) { | |
272 | if *x_i > τλ { | |
273 | *x_i -= τλ; | |
274 | if !*ast { | |
275 | act_changed = true; | |
276 | *ast = true; | |
277 | } | |
278 | n_inact += 1; | |
279 | } else { | |
280 | *x_i = 0.0; | |
281 | if *ast { | |
282 | act_changed = true; | |
283 | *ast = false; | |
284 | } | |
285 | } | |
286 | } | |
287 | ||
288 | // *** x now contains forward-backward step *** | |
289 | ||
290 | // 3. Solve SSN step `s`. | |
291 | // 3.1 Construct [τ A_{ℐ × ℐ}] if the set of inactive coordinates has changed. | |
292 | if act_changed { | |
293 | let decomp_iter = inact.iter().cartesian_product(inact.iter()).zip(mA.iter()); | |
294 | let decomp_constr = decomp_iter.filter_map(|((&i_inact, &j_inact), &mAij)| { | |
295 | //(i_inact && j_inact).then_some(mAij * τ) | |
296 | (i_inact && j_inact).then_some(mAij) // 🔺 below matches removal of τ | |
297 | }); | |
298 | let mat = DMatrix::from_iterator(n_inact, n_inact, decomp_constr); | |
299 | decomp = nalgebra::linalg::LU::new(mat); | |
300 | } | |
301 | ||
302 | // 3.2 Solve `s` = $s_ℐ^k$ from | |
303 | // $[τ A_{ℐ × ℐ}]s^k_ℐ = - x^k_ℐ + [G ∘ F](x^k)_ℐ - [τ A_{ℐ × 𝒜}]s^k_𝒜$. | |
304 | // With current variable setup we have $[G ∘ F](x^k) = $`x` and $x^k = x⁻$ = `xprev`, | |
305 | // so the system to solve is $[τ A_{ℐ × ℐ}]s^k_ℐ = (x-x⁻)_ℐ - [τ A_{ℐ × 𝒜}](x-x⁻)_𝒜$ | |
306 | // The matrix $[τ A_{ℐ × ℐ}]$ we have already LU-decomposed above into `decomp`. | |
307 | s = if n_inact > 0 { | |
308 | // 3.2.1 Construct `rhs` = $(x-x⁻)_ℐ - [τ A_{ℐ × 𝒜}](x-x⁻)_𝒜$ | |
309 | let inactfilt = inact.iter().copied(); | |
310 | let rhs_iter = izip!(x.iter(), xprev.iter(), mA.row_iter()).filter_zip(inactfilt); | |
311 | let rhs_constr = rhs_iter.map(|(&x_i, &xprev_i, mAi)| { | |
312 | // Calculate row i of [τ A_{ℐ × 𝒜}]s^k_𝒜 = [τ A_{ℐ × 𝒜}](x-xprev)_𝒜 | |
313 | let actfilt = inact.iter().copied().map(std::ops::Not::not); | |
314 | let actit = izip!(x.iter(), xprev.iter(), mAi.iter()).filter_zip(actfilt); | |
315 | let actpart = actit.map(|(&x_j, &xprev_j, &mAij)| { | |
316 | mAij * (x_j - xprev_j) | |
317 | }).sum(); | |
318 | // Subtract it from [x-prev]_i | |
319 | //x_i - xprev_i - τ * actpart | |
320 | (x_i - xprev_i) / τ - actpart // 🔺 change matches removal of τ above | |
321 | }); | |
322 | let mut rhs = DVector::from_iterator(n_inact, rhs_constr); | |
323 | assert_eq!(rhs.len(), n_inact); | |
324 | // Solve the system | |
325 | if !decomp.solve_mut(&mut rhs) { | |
326 | return Step::Failure(NumericalError( | |
327 | "Failed to solve linear system for subproblem SSN." | |
328 | )) | |
329 | } | |
330 | rhs | |
331 | } else { | |
332 | DVector::zeros(0) | |
333 | }; | |
334 | ||
335 | iters += 1; | |
336 | ||
337 | // 4. Report solution quality | |
338 | state.if_verbose(|| { | |
339 | // Calculate subdifferential at the FB step `x` that hasn't yet had `s` yet added. | |
340 | // The subdifferential of the objective is $Ax - g + λ + ∂ δ_{≥ 0}(x)$. | |
341 | // We return the minimal ∞-norm over all subderivatives. | |
342 | v.copy_from(g); // d = g | |
343 | mA.gemv(&mut v, 1.0, x, -1.0); // d = Ax - g | |
344 | let mut val = 0.0; | |
345 | for (&v_i, &x_i) in izip!(v.iter(), x.iter()) { | |
346 | let d = v_i + λ; | |
347 | if x_i > 0.0 || d < 0.0 { | |
348 | val = val.max(d.abs()); | |
349 | } | |
350 | } | |
351 | F::from_nalgebra_mixed(val) | |
352 | }) | |
353 | }); | |
354 | ||
355 | res.map(|_| iters) | |
356 | } | |
357 | ||
358 | /// This function applies an iterative method for the solution of the quadratic non-negativity | |
359 | /// constrained problem | |
360 | /// <div>$$ | |
361 | /// \min_{x ∈ ℝ^n} \frac{1}{2} x^⊤Ax - g^⊤ x + λ{\vec 1}^⊤ x + c + δ_{≥ 0}(x). | |
362 | /// $$</div> | |
363 | /// Semismooth Newton or forward-backward are supported based on the setting in `method`. | |
364 | /// The parameter `mA` is matrix $A$, and `g` and `λ` are as in the mathematical formulation. | |
365 | /// The constant $c$ does not need to be provided. The step length parameter is `τ` while | |
366 | /// `x` contains the initial iterate and on return the final one. The `iterator` controls | |
367 | /// stopping. The “verbose” value output by all methods is the $ℓ\_∞$ distance of some | |
368 | /// subdifferential of the objective to zero. | |
369 | /// | |
370 | /// Interior point methods could offer a further alternative, for example, the one in: | |
371 | /// | |
372 | /// * Valkonen T. - _A method for weighted projections to the positive definite | |
373 | /// cone_, <https://doi.org/10.1080/02331934.2014.929680>. | |
374 | /// | |
375 | /// This function returns the number of iterations taken. | |
4
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
376 | pub fn quadratic_nonneg<F, I, S, SA, SG>( |
0 | 377 | method : InnerMethod, |
4
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
378 | mA : &Matrix<F::MixedType, Dynamic, Dynamic, SA>, |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
379 | g : &Vector<F::MixedType, Dynamic, SG>, |
0 | 380 | //c_ : F, |
381 | λ : F, | |
4
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
382 | x : &mut Vector<F::MixedType, Dynamic, S>, |
0 | 383 | τ : F, |
384 | iterator : I | |
385 | ) -> usize | |
4
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
386 | where F : Float + ToNalgebraRealField<MixedType=F>, |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
387 | I : AlgIteratorFactory<F>, |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
388 | S : StorageMut<F::MixedType, Dynamic, Const::<1>>, |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
389 | SA : Storage<F::MixedType, Dynamic, Dynamic>, |
5aa5c279e341
Attempt at directly referring to measure masses as SliceStorageMut
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
390 | SG : Storage<F::MixedType, Dynamic, Const::<1>>, |
0 | 391 | { |
392 | ||
393 | match method { | |
394 | InnerMethod::FB => | |
395 | quadratic_nonneg_fb(mA, g, λ, x, τ, iterator), | |
396 | InnerMethod::SSN => | |
397 | quadratic_nonneg_ssn(mA, g, λ, x, τ, iterator).unwrap_or_else(|e| { | |
398 | println!("{}", format!("{e}. Using FB fallback.").red()); | |
399 | let ins = InnerSettings::<F>::default(); | |
400 | quadratic_nonneg_fb(mA, g, λ, x, τ, ins.iterator_options) | |
401 | }) | |
402 | } | |
403 | } |