Mercurial > repos > pointsource_algs / file comparison

comparison: src/subproblem/l1squared_unconstrained.rs

src/subproblem/l1squared_unconstrained.rs

changeset 72
e9a460a0e638
parent 54
b3312eee105c
equal deleted inserted replaced
71:e2953ffd4e0b 72:e9a460a0e638
1 /*! 1 /*!
2 Iterative algorithms for solving the finite-dimensional subproblem without constraints. 2 Iterative algorithms for solving the finite-dimensional subproblem without constraints.
3 */ 3 */
4 4
5 use itertools::izip; 5 use itertools::izip;
6 use nalgebra::DVector; 6 use nalgebra::{constraint::ShapeConstraint, DVector, Dyn, Storage, StorageMut, Vector, U1};
7 use numeric_literals::replace_float_literals; 7 use numeric_literals::replace_float_literals;
8 use std::cmp::Ordering::*; 8 use std::cmp::Ordering::*;
9 9
10 use alg_tools::iterate::{AlgIteratorFactory, AlgIteratorState}; 10 use alg_tools::iterate::{AlgIteratorFactory, AlgIteratorState};
11 use alg_tools::nalgebra_support::ToNalgebraRealField; 11 use alg_tools::nalgebra_support::{StridesOk, ToNalgebraRealField};
12 use alg_tools::nanleast::NaNLeast;
13 use alg_tools::norms::{Dist, L1}; 12 use alg_tools::norms::{Dist, L1};
14 use std::iter::zip; 13 use std::iter::zip;
15 14
16 use super::l1squared_nonneg::max_interval_dist_to_zero; 15 use super::l1squared_nonneg::max_interval_dist_to_zero;
17 use super::unconstrained::soft_thresholding; 16 use super::unconstrained::soft_thresholding;
42 /// sequentially: just check that the smallest assumed-nonzero component $i$ satisfies the 41 /// sequentially: just check that the smallest assumed-nonzero component $i$ satisfies the
43 /// condition of soft-thresholding to remain non-zero: $|x\_i|>τ\norm{x'}/(1+τm)$. 42 /// condition of soft-thresholding to remain non-zero: $|x\_i|>τ\norm{x'}/(1+τm)$.
44 /// Clearly, if this condition fails for $x\_i$, it will fail for all the components 43 /// Clearly, if this condition fails for $x\_i$, it will fail for all the components
45 /// already exluced. While, if it holds, it will hold for all components not excluded. 44 /// already exluced. While, if it holds, it will hold for all components not excluded.
46 #[replace_float_literals(F::cast_from(literal))] 45 #[replace_float_literals(F::cast_from(literal))]
47 pub(super) fn l1squared_prox<F: Float + nalgebra::RealField>( 46 pub(super) fn l1squared_prox<F: Float + nalgebra::RealField, S1, S2>(
48 sorted_abs: &mut DVector<F>, 47 sorted_abs: &mut DVector<F>,
49 x: &mut DVector<F>, 48 x: &mut Vector<F, Dyn, S1>,
50 y: &DVector<F>, 49 y: &Vector<F, Dyn, S2>,
51 β: F, 50 β: F,
52 ) { 51 ) where
52 S2: Storage<F, Dyn>,
53 S1: StorageMut<F, Dyn>,
54 {
55 //let orig_x = x.clone();
53 sorted_abs.copy_from(x); 56 sorted_abs.copy_from(x);
54 sorted_abs.axpy(-1.0, y, 1.0); 57 sorted_abs.axpy(-1.0, y, 1.0);
55 sorted_abs.apply(|z_i| *z_i = num_traits::abs(*z_i)); 58 sorted_abs.apply(|z_i| *z_i = num_traits::abs(*z_i));
56 sorted_abs 59 sorted_abs.as_mut_slice().sort_unstable_by(F::total_cmp);
57 .as_mut_slice()
58 .sort_unstable_by(|a, b| NaNLeast(*a).cmp(&NaNLeast(*b)));
59 60
60 let mut n = sorted_abs.sum(); 61 let mut n = sorted_abs.sum();
61 for (m, az_i) in zip((1..=x.len() as u32).rev(), sorted_abs) { 62 for (m, az_i) in zip((1..=x.len() as u32).rev(), sorted_abs) {
62 // test first 63 // test first. This is just *az_i <= tmp, for tmp defined below, without the division.
63 let tmp = β * n / (1.0 + β * F::cast_from(m)); 64 if *az_i <= β * (n - F::cast_from(m) * *az_i) {
64 if *az_i <= tmp {
65 // Fail 65 // Fail
66 n -= *az_i; 66 n -= *az_i;
67 } else { 67 } else {
68 // Success 68 // Success
69 let tmp = β * n / (1.0 + β * F::cast_from(m));
69 x.zip_apply(y, |x_i, y_i| { 70 x.zip_apply(y, |x_i, y_i| {
70 *x_i = y_i + soft_thresholding(*x_i - y_i, tmp) 71 *x_i = y_i + soft_thresholding(*x_i - y_i, tmp)
71 }); 72 });
73 // //Check 0 ∈ w-x + β\norm{w-y}\_1\sign (w-y).
74 // let n: F = izip!(x.iter(), y)
75 // .map(|(&w_i, &y_i)| NumTraitsFloat::abs(w_i - y_i))
76 // .sum();
77 // for (&mut w_i, &x_i, &y_i) in izip!(x, &orig_x, y) {
78 // if w_i > y_i {
79 // assert_lt!(NumTraitsFloat::abs(w_i - x_i + n * β), 10.0 * F::EPSILON);
80 // } else if w_i < y_i {
81 // assert_lt!(NumTraitsFloat::abs(w_i - x_i - n * β), 10.0 * F::EPSILON);
82 // } else {
83 // assert_lt!(-n * β - 10.0 * F::EPSILON, w_i - x_i);
84 // assert_lt!(w_i - x_i, n * β + 10.0 * F::EPSILON);
85 // }
86 // }
72 return; 87 return;
73 } 88 }
74 } 89 }
75 // m = 0 should always work, but x is zero. 90 // m = 0 should always work, but x is zero.
76 x.fill(0.0); 91 x.fill(0.0);
78 93
79 /// Returns the ∞-norm minimal subdifferential of $x ↦ (β/2)|x-y|_1^2 - g^⊤ x + λ\|x\|₁$ at $x$. 94 /// Returns the ∞-norm minimal subdifferential of $x ↦ (β/2)|x-y|_1^2 - g^⊤ x + λ\|x\|₁$ at $x$.
80 /// 95 ///
81 /// `v` will be modified and cannot be trusted to contain useful values afterwards. 96 /// `v` will be modified and cannot be trusted to contain useful values afterwards.
82 #[replace_float_literals(F::cast_from(literal))] 97 #[replace_float_literals(F::cast_from(literal))]
83 fn min_subdifferential<F: Float + nalgebra::RealField>( 98 fn min_subdifferential<F: Float + nalgebra::RealField, S1, S2, S3>(
84 y: &DVector<F>, 99 y: &Vector<F, Dyn, S1>,
85 x: &DVector<F>, 100 x: &Vector<F, Dyn, S2>,
86 g: &DVector<F>, 101 g: &Vector<F, Dyn, S3>,
87 λ: F, 102 λ: F,
88 β: F, 103 ) -> F
89 ) -> F { 104 where
105 S1: Storage<F, Dyn>,
106 S2: Storage<F, Dyn>,
107 S3: Storage<F, Dyn>,
108 ShapeConstraint: StridesOk<F, Dyn, U1, S2>,
109 {
90 let mut val = 0.0; 110 let mut val = 0.0;
91 let tmp = β * y.dist(x, L1); 111 let tmp = y.dist(x, L1);
92 for (&g_i, &x_i, y_i) in izip!(g.iter(), x.iter(), y.iter()) { 112 for (&g_i, &x_i, y_i) in izip!(g.iter(), x.iter(), y.iter()) {
93 let (mut lb, mut ub) = (-g_i, -g_i); 113 let (mut lb, mut ub) = (-g_i, -g_i);
94 match x_i.partial_cmp(y_i) { 114 match x_i.partial_cmp(y_i) {
95 Some(Greater) => { 115 Some(Greater) => {
96 lb += tmp; 116 lb += tmp;
130 /// For detailed documentation of the inputs and outputs, refer to there. 150 /// For detailed documentation of the inputs and outputs, refer to there.
131 /// 151 ///
132 /// The `λ` component of the model is handled in the proximal step instead of the gradient step 152 /// The `λ` component of the model is handled in the proximal step instead of the gradient step
133 /// for potential performance improvements. 153 /// for potential performance improvements.
134 #[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())] 154 #[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())]
135 pub fn l1squared_unconstrained_pdps<F, I>( 155 pub fn l1squared_unconstrained_pdps<F, I, S1, S2, S3>(
136 y: &DVector<F::MixedType>, 156 y: &Vector<F::MixedType, Dyn, S1>,
137 g: &DVector<F::MixedType>, 157 g: &Vector<F::MixedType, Dyn, S2>,
138 λ_: F, 158 λ_: F,
139 β_: F, 159 x: &mut Vector<F::MixedType, Dyn, S3>,
140 x: &mut DVector<F::MixedType>,
141 τ_: F, 160 τ_: F,
142 σ_: F, 161 σ_: F,
143 iterator: I, 162 iterator: I,
144 ) -> usize 163 ) -> usize
145 where 164 where
146 F: Float + ToNalgebraRealField, 165 F: Float + ToNalgebraRealField,
147 I: AlgIteratorFactory<F>, 166 I: AlgIteratorFactory<F>,
167 S1: Storage<F::MixedType, Dyn>,
168 S2: Storage<F::MixedType, Dyn>,
169 S3: StorageMut<F::MixedType, Dyn>,
170 ShapeConstraint: StridesOk<F::MixedType, Dyn, U1, S3>,
148 { 171 {
149 let λ = λ_.to_nalgebra_mixed(); 172 let λ = λ_.to_nalgebra_mixed();
150 let β = β_.to_nalgebra_mixed();
151 let τ = τ_.to_nalgebra_mixed(); 173 let τ = τ_.to_nalgebra_mixed();
152 let σ = σ_.to_nalgebra_mixed(); 174 let σ = σ_.to_nalgebra_mixed();
153 let mut w = DVector::zeros(x.len()); 175 let mut w = DVector::zeros(x.len());
154 let mut tmp = DVector::zeros(x.len()); 176 let mut tmp = DVector::zeros(x.len());
155 let mut xprev = x.clone(); 177 let mut xprev = x.clone_owned();
156 let mut iters = 0; 178 let mut iters = 0;
157 179
158 iterator.iterate(|state| { 180 iterator.iterate(|state| {
159 // Primal step: x^{k+1} = prox_{τ|.-y|_1^2}(x^k - τ (w^k - g)) 181 // Primal step: x^{k+1} = prox_{τ|.-y|_1^2}(x^k - τ (w^k - g))
160 x.axpy(-τ, &w, 1.0); 182 x.axpy(-τ, &w, 1.0);
161 x.axpy(τ, g, 1.0); 183 x.axpy(τ, g, 1.0);
162 l1squared_prox(&mut tmp, x, y, τ * β); 184 l1squared_prox(&mut tmp, x, y, τ);
163 185
164 // Dual step: w^{k+1} = proj_{[-λ,λ]}(w^k + σ(2x^{k+1}-x^k)) 186 // Dual step: w^{k+1} = proj_{[-λ,λ]}(w^k + σ(2x^{k+1}-x^k))
165 w.axpy(2.0 * σ, x, 1.0); 187 w.axpy(2.0 * σ, x, 1.0);
166 w.axpy(-σ, &xprev, 1.0); 188 w.axpy(-σ, &xprev, 1.0);
167 w.apply(|w_i| *w_i = num_traits::clamp(*w_i, -λ, λ)); 189 w.apply(|w_i| *w_i = num_traits::clamp(*w_i, -λ, λ));
168 xprev.copy_from(x); 190 xprev.copy_from(x);
169 191
170 iters += 1; 192 iters += 1;
171 193
172 state.if_verbose(|| F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ, β))) 194 state.if_verbose(|| F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ)))
173 }); 195 });
174 196
175 iters 197 iters
176 } 198 }
177 199
191 /// & = \min_{x ∈ ℝ^n} \max_{w} ⟨θ w, x⟩ - g^⊤ x + λ\|x\|₁ 213 /// & = \min_{x ∈ ℝ^n} \max_{w} ⟨θ w, x⟩ - g^⊤ x + λ\|x\|₁
192 /// - \left(x ↦ \frac{β}{2θ} |x-y|_1^2 \right)^*(w). 214 /// - \left(x ↦ \frac{β}{2θ} |x-y|_1^2 \right)^*(w).
193 /// \end{split} 215 /// \end{split}
194 /// $$</div> 216 /// $$</div>
195 #[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())] 217 #[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())]
196 pub fn l1squared_unconstrained_pdps_alt<F, I>( 218 pub fn l1squared_unconstrained_pdps_alt<F, I, S1, S2, S3>(
197 y: &DVector<F::MixedType>, 219 y: &Vector<F::MixedType, Dyn, S1>,
198 g: &DVector<F::MixedType>, 220 g: &Vector<F::MixedType, Dyn, S2>,
199 λ_: F, 221 λ_: F,
200 β_: F, 222 x: &mut Vector<F::MixedType, Dyn, S3>,
201 x: &mut DVector<F::MixedType>,
202 τ_: F, 223 τ_: F,
203 σ_: F, 224 σ_: F,
204 θ_: F, 225 θ_: F,
205 iterator: I, 226 iterator: I,
206 ) -> usize 227 ) -> usize
207 where 228 where
208 F: Float + ToNalgebraRealField, 229 F: Float + ToNalgebraRealField,
209 I: AlgIteratorFactory<F>, 230 I: AlgIteratorFactory<F>,
231 S1: Storage<F::MixedType, Dyn>,
232 S2: Storage<F::MixedType, Dyn>,
233 S3: StorageMut<F::MixedType, Dyn>,
234 ShapeConstraint: StridesOk<F::MixedType, Dyn, U1, S3>,
210 { 235 {
211 let λ = λ_.to_nalgebra_mixed(); 236 let λ = λ_.to_nalgebra_mixed();
212 let τ = τ_.to_nalgebra_mixed(); 237 let τ = τ_.to_nalgebra_mixed();
213 let σ = σ_.to_nalgebra_mixed(); 238 let σ = σ_.to_nalgebra_mixed();
214 let θ = θ_.to_nalgebra_mixed(); 239 let θ = θ_.to_nalgebra_mixed();
215 let β = β_.to_nalgebra_mixed();
216 let σθ = σ * θ; 240 let σθ = σ * θ;
217 let τθ = τ * θ; 241 let τθ = τ * θ;
218 let mut w = DVector::zeros(x.len()); 242 let mut w = DVector::zeros(x.len());
219 let mut tmp = DVector::zeros(x.len()); 243 let mut tmp = DVector::zeros(x.len());
220 let mut xprev = x.clone(); 244 let mut xprev = x.clone_owned();
221 let mut iters = 0; 245 let mut iters = 0;
222 246
223 iterator.iterate(|state| { 247 iterator.iterate(|state| {
224 // Primal step: x^{k+1} = soft_τλ(x^k - τ(θ w^k -g)) 248 // Primal step: x^{k+1} = soft_τλ(x^k - τ(θ w^k -g))
225 x.axpy(-τθ, &w, 1.0); 249 x.axpy(-τθ, &w, 1.0);
234 // where q/σ = w^k/σ + (2x^{k+1}-x^k), 258 // where q/σ = w^k/σ + (2x^{k+1}-x^k),
235 w /= σ; 259 w /= σ;
236 w.axpy(2.0, x, 1.0); 260 w.axpy(2.0, x, 1.0);
237 w.axpy(-1.0, &xprev, 1.0); 261 w.axpy(-1.0, &xprev, 1.0);
238 xprev.copy_from(&w); // use xprev as temporary variable 262 xprev.copy_from(&w); // use xprev as temporary variable
239 l1squared_prox(&mut tmp, &mut xprev, y, β / σθ); 263 l1squared_prox(&mut tmp, &mut xprev, y, 1.0 / σθ);
240 w -= &xprev; 264 w -= &xprev;
241 w *= σ; 265 w *= σ;
242 xprev.copy_from(x); 266 xprev.copy_from(x);
243 267
244 iters += 1; 268 iters += 1;
245 269
246 state.if_verbose(|| F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ, β))) 270 state.if_verbose(|| F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ)))
247 }); 271 });
248 272
249 iters 273 iters
250 } 274 }
251 275
255 /// $$</div> 279 /// $$</div>
256 /// Only PDPS is supported. 280 /// Only PDPS is supported.
257 /// 281 ///
258 /// This function returns the number of iterations taken. 282 /// This function returns the number of iterations taken.
259 #[replace_float_literals(F::cast_from(literal))] 283 #[replace_float_literals(F::cast_from(literal))]
260 pub fn l1squared_unconstrained<F, I>( 284 pub fn l1squared_unconstrained<F, I, S1, S2, S3>(
261 y: &DVector<F::MixedType>, 285 y: &Vector<F::MixedType, Dyn, S1>,
262 g: &DVector<F::MixedType>, 286 g: &Vector<F::MixedType, Dyn, S2>,
263 λ: F, 287 λ: F,
264 β: F, 288 x: &mut Vector<F::MixedType, Dyn, S3>,
265 x: &mut DVector<F::MixedType>,
266 inner: &InnerSettings<F>, 289 inner: &InnerSettings<F>,
267 iterator: I, 290 iterator: I,
268 ) -> usize 291 ) -> usize
269 where 292 where
270 F: Float + ToNalgebraRealField, 293 F: Float + ToNalgebraRealField,
271 I: AlgIteratorFactory<F>, 294 I: AlgIteratorFactory<F>,
295 S1: Storage<F::MixedType, Dyn>,
296 S2: Storage<F::MixedType, Dyn>,
297 S3: StorageMut<F::MixedType, Dyn>,
298 ShapeConstraint: StridesOk<F::MixedType, Dyn, U1, S3>,
272 { 299 {
273 // Estimate of ‖K‖ for K=θ Id. 300 // Estimate of ‖K‖ for K=θ Id.
274 let inner_θ = 1.0; 301 let inner_θ = 1.0;
275 let normest = inner_θ; 302 let normest = inner_θ;
276 303
277 let (inner_τ, inner_σ) = (inner.pdps_τσ0.0 / normest, inner.pdps_τσ0.1 / normest); 304 let (inner_τ, inner_σ) = (inner.pdps_τσ0.0 / normest, inner.pdps_τσ0.1 / normest);
278 305
279 match inner.method { 306 match inner.method {
280 InnerMethod::PDPS => { 307 InnerMethod::PDPS => {
281 l1squared_unconstrained_pdps_alt(y, g, λ, β, x, inner_τ, inner_σ, inner_θ, iterator) 308 l1squared_unconstrained_pdps_alt(y, g, λ, x, inner_τ, inner_σ, inner_θ, iterator)
282 } 309 }
283 other => unimplemented!("${other:?} is unimplemented"), 310 other => unimplemented!("${other:?} is unimplemented"),
284 } 311 }
285 } 312 }

Mercurial Repository: pointsource_algs

mercurial