Mercurial > repos > alg_tools / file comparison

comparison: src/convex.rs

src/convex.rs

branch
dev
changeset 104
e7f1cb4bec78
parent 88
086a59b3a2b4
child 105
103aa137fcb2
equal deleted inserted replaced
103:e98e1da2530d 104:e7f1cb4bec78
1 /*! 1 /*!
2 Some convex analysis basics 2 Some convex analysis basics
3 */ 3 */
4 4
5 use crate::euclidean::Euclidean;
6 use crate::instance::{DecompositionMut, Instance, InstanceMut};
7 use crate::linops::{IdOp, Scaled};
8 use crate::mapping::{Mapping, Space};
9 use crate::norms::*;
10 use crate::operator_arithmetic::{Constant, Weighted};
11 use crate::types::*;
5 use std::marker::PhantomData; 12 use std::marker::PhantomData;
6 use crate::types::*;
7 use crate::mapping::{Mapping, Space};
8 use crate::linops::IdOp;
9 use crate::instance::{Instance, InstanceMut, DecompositionMut};
10 use crate::operator_arithmetic::{Constant, Weighted};
11 use crate::norms::*;
12 13
13 /// Trait for convex mappings. Has no features, just serves as a constraint 14 /// Trait for convex mappings. Has no features, just serves as a constraint
14 /// 15 ///
15 /// TODO: should constrain `Mapping::Codomain` to implement a partial order, 16 /// TODO: should constrain `Mapping::Codomain` to implement a partial order,
16 /// but this makes everything complicated with little benefit. 17 /// but this makes everything complicated with little benefit.
17 pub trait ConvexMapping<Domain : Space, F : Num = f64> : Mapping<Domain, Codomain = F> 18 pub trait ConvexMapping<Domain: Space, F: Num = f64>: Mapping<Domain, Codomain = F> {}
18 {}
19 19
20 /// Trait for mappings with a Fenchel conjugate 20 /// Trait for mappings with a Fenchel conjugate
21 /// 21 ///
22 /// The conjugate type has to implement [`ConvexMapping`], but a `Conjugable` mapping need 22 /// The conjugate type has to implement [`ConvexMapping`], but a `Conjugable` mapping need
23 /// not be convex. 23 /// not be convex.
24 pub trait Conjugable<Domain : HasDual<F>, F : Num = f64> : Mapping<Domain> { 24 pub trait Conjugable<Domain: HasDual<F>, F: Num = f64>: Mapping<Domain> {
25 type Conjugate<'a> : ConvexMapping<Domain::DualSpace, F> where Self : 'a; 25 type Conjugate<'a>: ConvexMapping<Domain::DualSpace, F>
26 where
27 Self: 'a;
26 28
27 fn conjugate(&self) -> Self::Conjugate<'_>; 29 fn conjugate(&self) -> Self::Conjugate<'_>;
28 } 30 }
29 31
30 /// Trait for mappings with a Fenchel preconjugate 32 /// Trait for mappings with a Fenchel preconjugate
31 /// 33 ///
32 /// In contrast to [`Conjugable`], the preconjugate need not implement [`ConvexMapping`], 34 /// In contrast to [`Conjugable`], the preconjugate need not implement [`ConvexMapping`],
33 /// but a `Preconjugable` mapping has to be convex. 35 /// but a `Preconjugable` mapping has to be convex.
34 pub trait Preconjugable<Domain, Predual, F : Num = f64> : ConvexMapping<Domain, F> 36 pub trait Preconjugable<Domain, Predual, F: Num = f64>: ConvexMapping<Domain, F>
35 where 37 where
36 Domain : Space, 38 Domain: Space,
37 Predual : HasDual<F> 39 Predual: HasDual<F>,
38 { 40 {
39 type Preconjugate<'a> : Mapping<Predual> where Self : 'a; 41 type Preconjugate<'a>: Mapping<Predual>
42 where
43 Self: 'a;
40 44
41 fn preconjugate(&self) -> Self::Preconjugate<'_>; 45 fn preconjugate(&self) -> Self::Preconjugate<'_>;
42 } 46 }
43 47
44 /// Trait for mappings with a proximap map 48 /// Trait for mappings with a proximap map
45 /// 49 ///
46 /// The conjugate type has to implement [`ConvexMapping`], but a `Conjugable` mapping need 50 /// The conjugate type has to implement [`ConvexMapping`], but a `Conjugable` mapping need
47 /// not be convex. 51 /// not be convex.
48 pub trait Prox<Domain : Space> : Mapping<Domain> { 52 pub trait Prox<Domain: Space>: Mapping<Domain> {
49 type Prox<'a> : Mapping<Domain, Codomain=Domain> where Self : 'a; 53 type Prox<'a>: Mapping<Domain, Codomain = Domain>
54 where
55 Self: 'a;
50 56
51 /// Returns a proximal mapping with weight τ 57 /// Returns a proximal mapping with weight τ
52 fn prox_mapping(&self, τ : Self::Codomain) -> Self::Prox<'_>; 58 fn prox_mapping(&self, τ: Self::Codomain) -> Self::Prox<'_>;
53 59
54 /// Calculate the proximal mapping with weight τ 60 /// Calculate the proximal mapping with weight τ
55 fn prox<I : Instance<Domain>>(&self, τ : Self::Codomain, z : I) -> Domain { 61 fn prox<I: Instance<Domain>>(&self, τ: Self::Codomain, z: I) -> Domain {
56 self.prox_mapping(τ).apply(z) 62 self.prox_mapping(τ).apply(z)
57 } 63 }
58 64
59 /// Calculate the proximal mapping with weight τ in-place 65 /// Calculate the proximal mapping with weight τ in-place
60 fn prox_mut<'b>(&self, τ : Self::Codomain, y : &'b mut Domain) 66 fn prox_mut<'b>(&self, τ: Self::Codomain, y: &'b mut Domain)
61 where 67 where
62 &'b mut Domain : InstanceMut<Domain>, 68 &'b mut Domain: InstanceMut<Domain>,
63 Domain:: Decomp : DecompositionMut<Domain>, 69 Domain::Decomp: DecompositionMut<Domain>,
64 for<'a> &'a Domain : Instance<Domain>, 70 for<'a> &'a Domain: Instance<Domain>,
65 { 71 {
66 *y = self.prox(τ, &*y); 72 *y = self.prox(τ, &*y);
67 } 73 }
68 } 74 }
69 75
70
71 /// Constraint to the unit ball of the norm described by `E`. 76 /// Constraint to the unit ball of the norm described by `E`.
72 pub struct NormConstraint<F : Float, E : NormExponent> { 77 pub struct NormConstraint<F: Float, E: NormExponent> {
73 radius : F, 78 radius: F,
74 norm : NormMapping<F, E>, 79 norm: NormMapping<F, E>,
75 } 80 }
76 81
77 impl<Domain, E, F> ConvexMapping<Domain, F> for NormMapping<F, E> 82 impl<Domain, E, F> ConvexMapping<Domain, F> for NormMapping<F, E>
78 where 83 where
79 Domain : Space, 84 Domain: Space,
80 E : NormExponent, 85 E: NormExponent,
81 F : Float, 86 F: Float,
82 Self : Mapping<Domain, Codomain=F> 87 Self: Mapping<Domain, Codomain = F>,
83 {} 88 {
84 89 }
85 90
86 impl<F, E, Domain> Mapping<Domain> for NormConstraint<F, E> 91 impl<F, E, Domain> Mapping<Domain> for NormConstraint<F, E>
87 where 92 where
88 Domain : Space + Norm<F, E>, 93 Domain: Space + Norm<F, E>,
89 F : Float, 94 F: Float,
90 E : NormExponent, 95 E: NormExponent,
91 { 96 {
92 type Codomain = F; 97 type Codomain = F;
93 98
94 fn apply<I : Instance<Domain>>(&self, d : I) -> F { 99 fn apply<I: Instance<Domain>>(&self, d: I) -> F {
95 if d.eval(|x| x.norm(self.norm.exponent)) <= self.radius { 100 if d.eval(|x| x.norm(self.norm.exponent)) <= self.radius {
96 F::ZERO 101 F::ZERO
97 } else { 102 } else {
98 F::INFINITY 103 F::INFINITY
99 } 104 }
100 } 105 }
101 } 106 }
102 107
103 impl<Domain, E, F> ConvexMapping<Domain, F> for NormConstraint<F, E> 108 impl<Domain, E, F> ConvexMapping<Domain, F> for NormConstraint<F, E>
104 where 109 where
105 Domain : Space, 110 Domain: Space,
106 E : NormExponent, 111 E: NormExponent,
107 F : Float, 112 F: Float,
108 Self : Mapping<Domain, Codomain=F> 113 Self: Mapping<Domain, Codomain = F>,
109 {} 114 {
110 115 }
111 116
112 impl<E, F, Domain> Conjugable<Domain, F> for NormMapping<F, E> 117 impl<E, F, Domain> Conjugable<Domain, F> for NormMapping<F, E>
113 where 118 where
114 E : HasDualExponent, 119 E: HasDualExponent,
115 F : Float, 120 F: Float,
116 Domain : HasDual<F> + Norm<F, E> + Space, 121 Domain: HasDual<F> + Norm<F, E> + Space,
117 <Domain as HasDual<F>>::DualSpace : Norm<F, E::DualExp> 122 <Domain as HasDual<F>>::DualSpace: Norm<F, E::DualExp>,
118 { 123 {
119 type Conjugate<'a> = NormConstraint<F, E::DualExp> where Self : 'a; 124 type Conjugate<'a>
125 = NormConstraint<F, E::DualExp>
126 where
127 Self: 'a;
120 128
121 fn conjugate(&self) -> Self::Conjugate<'_> { 129 fn conjugate(&self) -> Self::Conjugate<'_> {
122 NormConstraint { 130 NormConstraint {
123 radius : F::ONE, 131 radius: F::ONE,
124 norm : self.exponent.dual_exponent().as_mapping() 132 norm: self.exponent.dual_exponent().as_mapping(),
125 } 133 }
126 } 134 }
127 } 135 }
128 136
129 impl<C, E, F, Domain> Conjugable<Domain, F> for Weighted<NormMapping<F, E>, C> 137 impl<C, E, F, Domain> Conjugable<Domain, F> for Weighted<NormMapping<F, E>, C>
130 where 138 where
131 C : Constant<Type = F>, 139 C: Constant<Type = F>,
132 E : HasDualExponent, 140 E: HasDualExponent,
133 F : Float, 141 F: Float,
134 Domain : HasDual<F> + Norm<F, E> + Space, 142 Domain: HasDual<F> + Norm<F, E> + Space,
135 <Domain as HasDual<F>>::DualSpace : Norm<F, E::DualExp> 143 <Domain as HasDual<F>>::DualSpace: Norm<F, E::DualExp>,
136 { 144 {
137 type Conjugate<'a> = NormConstraint<F, E::DualExp> where Self : 'a; 145 type Conjugate<'a>
146 = NormConstraint<F, E::DualExp>
147 where
148 Self: 'a;
138 149
139 fn conjugate(&self) -> Self::Conjugate<'_> { 150 fn conjugate(&self) -> Self::Conjugate<'_> {
140 NormConstraint { 151 NormConstraint {
141 radius : self.weight.value(), 152 radius: self.weight.value(),
142 norm : self.base_fn.exponent.dual_exponent().as_mapping() 153 norm: self.base_fn.exponent.dual_exponent().as_mapping(),
143 } 154 }
144 } 155 }
145 } 156 }
146 157
147 impl<Domain, E, F> Prox<Domain> for NormConstraint<F, E> 158 impl<Domain, E, F> Prox<Domain> for NormConstraint<F, E>
148 where 159 where
149 Domain : Space + Norm<F, E>, 160 Domain: Space + Norm<F, E>,
150 E : NormExponent, 161 E: NormExponent,
151 F : Float, 162 F: Float,
152 NormProjection<F, E> : Mapping<Domain, Codomain=Domain>, 163 NormProjection<F, E>: Mapping<Domain, Codomain = Domain>,
153 { 164 {
154 type Prox<'a> = NormProjection<F, E> where Self : 'a; 165 type Prox<'a>
155 166 = NormProjection<F, E>
156 #[inline] 167 where
157 fn prox_mapping(&self, _τ : Self::Codomain) -> Self::Prox<'_> { 168 Self: 'a;
169
170 #[inline]
171 fn prox_mapping(&self, _τ: Self::Codomain) -> Self::Prox<'_> {
158 assert!(self.radius >= F::ZERO); 172 assert!(self.radius >= F::ZERO);
159 NormProjection{ radius : self.radius, exponent : self.norm.exponent } 173 NormProjection {
174 radius: self.radius,
175 exponent: self.norm.exponent,
176 }
160 } 177 }
161 } 178 }
162 179
163 /// Projection to the unit ball of the norm described by `E`. 180 /// Projection to the unit ball of the norm described by `E`.
164 pub struct NormProjection<F : Float, E : NormExponent> { 181 pub struct NormProjection<F: Float, E: NormExponent> {
165 radius : F, 182 radius: F,
166 exponent : E, 183 exponent: E,
167 } 184 }
168 185
169 /* 186 /*
170 impl<F, Domain> Mapping<Domain> for NormProjection<F, L2> 187 impl<F, Domain> Mapping<Domain> for NormProjection<F, L2>
171 where 188 where
180 } 197 }
181 */ 198 */
182 199
183 impl<F, E, Domain> Mapping<Domain> for NormProjection<F, E> 200 impl<F, E, Domain> Mapping<Domain> for NormProjection<F, E>
184 where 201 where
185 Domain : Space + Projection<F, E>, 202 Domain: Space + Projection<F, E>,
186 F : Float, 203 F: Float,
187 E : NormExponent, 204 E: NormExponent,
188 { 205 {
189 type Codomain = Domain; 206 type Codomain = Domain;
190 207
191 fn apply<I : Instance<Domain>>(&self, d : I) -> Domain { 208 fn apply<I: Instance<Domain>>(&self, d: I) -> Domain {
192 d.own().proj_ball(self.radius, self.exponent) 209 d.own().proj_ball(self.radius, self.exponent)
193 } 210 }
194 } 211 }
195 212
196
197 /// The zero mapping 213 /// The zero mapping
198 pub struct Zero<Domain : Space, F : Num>(PhantomData<(Domain, F)>); 214 pub struct Zero<Domain: Space, F: Num>(PhantomData<(Domain, F)>);
199 215
200 impl<Domain : Space, F : Num> Zero<Domain, F> { 216 impl<Domain: Space, F: Num> Zero<Domain, F> {
201 pub fn new() -> Self { 217 pub fn new() -> Self {
202 Zero(PhantomData) 218 Zero(PhantomData)
203 } 219 }
204 } 220 }
205 221
206 impl<Domain : Space, F : Num> Mapping<Domain> for Zero<Domain, F> { 222 impl<Domain: Space, F: Num> Mapping<Domain> for Zero<Domain, F> {
207 type Codomain = F; 223 type Codomain = F;
208 224
209 /// Compute the value of `self` at `x`. 225 /// Compute the value of `self` at `x`.
210 fn apply<I : Instance<Domain>>(&self, _x : I) -> Self::Codomain { 226 fn apply<I: Instance<Domain>>(&self, _x: I) -> Self::Codomain {
211 F::ZERO 227 F::ZERO
212 } 228 }
213 } 229 }
214 230
215 impl<Domain : Space, F : Num> ConvexMapping<Domain, F> for Zero<Domain, F> { } 231 impl<Domain: Space, F: Num> ConvexMapping<Domain, F> for Zero<Domain, F> {}
216 232
217 233 impl<Domain: HasDual<F>, F: Float> Conjugable<Domain, F> for Zero<Domain, F> {
218 impl<Domain : HasDual<F>, F : Float> Conjugable<Domain, F> for Zero<Domain, F> { 234 type Conjugate<'a>
219 type Conjugate<'a> = ZeroIndicator<Domain::DualSpace, F> where Self : 'a; 235 = ZeroIndicator<Domain::DualSpace, F>
236 where
237 Self: 'a;
220 238
221 #[inline] 239 #[inline]
222 fn conjugate(&self) -> Self::Conjugate<'_> { 240 fn conjugate(&self) -> Self::Conjugate<'_> {
223 ZeroIndicator::new() 241 ZeroIndicator::new()
224 } 242 }
225 } 243 }
226 244
227 impl<Domain, Predual, F : Float> Preconjugable<Domain, Predual, F> for Zero<Domain, F> 245 impl<Domain, Predual, F: Float> Preconjugable<Domain, Predual, F> for Zero<Domain, F>
228 where 246 where
229 Domain : Space, 247 Domain: Space,
230 Predual : HasDual<F> 248 Predual: HasDual<F>,
231 { 249 {
232 type Preconjugate<'a> = ZeroIndicator<Predual, F> where Self : 'a; 250 type Preconjugate<'a>
251 = ZeroIndicator<Predual, F>
252 where
253 Self: 'a;
233 254
234 #[inline] 255 #[inline]
235 fn preconjugate(&self) -> Self::Preconjugate<'_> { 256 fn preconjugate(&self) -> Self::Preconjugate<'_> {
236 ZeroIndicator::new() 257 ZeroIndicator::new()
237 } 258 }
238 } 259 }
239 260
240 impl<Domain : Space + Clone, F : Num> Prox<Domain> for Zero<Domain, F> { 261 impl<Domain: Space + Clone, F: Num> Prox<Domain> for Zero<Domain, F> {
241 type Prox<'a> = IdOp<Domain> where Self : 'a; 262 type Prox<'a>
242 263 = IdOp<Domain>
243 #[inline] 264 where
244 fn prox_mapping(&self, _τ : Self::Codomain) -> Self::Prox<'_> { 265 Self: 'a;
266
267 #[inline]
268 fn prox_mapping(&self, _τ: Self::Codomain) -> Self::Prox<'_> {
245 IdOp::new() 269 IdOp::new()
246 } 270 }
247 } 271 }
248 272
249
250 /// The zero indicator 273 /// The zero indicator
251 pub struct ZeroIndicator<Domain : Space, F : Num>(PhantomData<(Domain, F)>); 274 pub struct ZeroIndicator<Domain: Space, F: Num>(PhantomData<(Domain, F)>);
252 275
253 impl<Domain : Space, F : Num> ZeroIndicator<Domain, F> { 276 impl<Domain: Space, F: Num> ZeroIndicator<Domain, F> {
254 pub fn new() -> Self { 277 pub fn new() -> Self {
255 ZeroIndicator(PhantomData) 278 ZeroIndicator(PhantomData)
256 } 279 }
257 } 280 }
258 281
259 impl<Domain : Normed<F>, F : Float> Mapping<Domain> for ZeroIndicator<Domain, F> { 282 impl<Domain: Normed<F>, F: Float> Mapping<Domain> for ZeroIndicator<Domain, F> {
260 type Codomain = F; 283 type Codomain = F;
261 284
262 /// Compute the value of `self` at `x`. 285 /// Compute the value of `self` at `x`.
263 fn apply<I : Instance<Domain>>(&self, x : I) -> Self::Codomain { 286 fn apply<I: Instance<Domain>>(&self, x: I) -> Self::Codomain {
264 x.eval(|x̃| if x̃.is_zero() { F::ZERO } else { F::INFINITY }) 287 x.eval(|x̃| if x̃.is_zero() { F::ZERO } else { F::INFINITY })
265 } 288 }
266 } 289 }
267 290
268 impl<Domain : Normed<F>, F : Float> ConvexMapping<Domain, F> for ZeroIndicator<Domain, F> { } 291 impl<Domain: Normed<F>, F: Float> ConvexMapping<Domain, F> for ZeroIndicator<Domain, F> {}
269 292
270 impl<Domain : HasDual<F>, F : Float> Conjugable<Domain, F> for ZeroIndicator<Domain, F> { 293 impl<Domain: HasDual<F>, F: Float> Conjugable<Domain, F> for ZeroIndicator<Domain, F> {
271 type Conjugate<'a> = Zero<Domain::DualSpace, F> where Self : 'a; 294 type Conjugate<'a>
295 = Zero<Domain::DualSpace, F>
296 where
297 Self: 'a;
272 298
273 #[inline] 299 #[inline]
274 fn conjugate(&self) -> Self::Conjugate<'_> { 300 fn conjugate(&self) -> Self::Conjugate<'_> {
275 Zero::new() 301 Zero::new()
276 } 302 }
277 } 303 }
278 304
279 impl<Domain, Predual, F : Float> Preconjugable<Domain, Predual, F> for ZeroIndicator<Domain, F> 305 impl<Domain, Predual, F: Float> Preconjugable<Domain, Predual, F> for ZeroIndicator<Domain, F>
280 where 306 where
281 Domain : Normed<F>, 307 Domain: Normed<F>,
282 Predual : HasDual<F> 308 Predual: HasDual<F>,
283 { 309 {
284 type Preconjugate<'a> = Zero<Predual, F> where Self : 'a; 310 type Preconjugate<'a>
311 = Zero<Predual, F>
312 where
313 Self: 'a;
285 314
286 #[inline] 315 #[inline]
287 fn preconjugate(&self) -> Self::Preconjugate<'_> { 316 fn preconjugate(&self) -> Self::Preconjugate<'_> {
288 Zero::new() 317 Zero::new()
289 } 318 }
290 } 319 }
320
321 /// The squared Euclidean norm divided by two
322 pub struct Norm222<Domain: Space, F: Float>(PhantomData<(Domain, F)>);
323
324 impl<Domain: Euclidean<F>, F: Float> Norm222<Domain, F> {
325 pub fn new() -> Self {
326 Norm222(PhantomData)
327 }
328 }
329
330 impl<Domain: Euclidean<F>, F: Float> Mapping<Domain> for Norm222<Domain, F> {
331 type Codomain = F;
332
333 /// Compute the value of `self` at `x`.
334 fn apply<I: Instance<Domain>>(&self, x: I) -> Self::Codomain {
335 x.eval(|z| z.norm2_squared() / F::TWO)
336 }
337 }
338
339 impl<Domain: Euclidean<F>, F: Float> ConvexMapping<Domain, F> for Norm222<Domain, F> {}
340
341 impl<Domain: Euclidean<F>, F: Float> Conjugable<Domain, F> for Norm222<Domain, F> {
342 type Conjugate<'a>
343 = Self
344 where
345 Self: 'a;
346
347 #[inline]
348 fn conjugate(&self) -> Self::Conjugate<'_> {
349 Self::new()
350 }
351 }
352
353 impl<Domain: Euclidean<F>, F: Float> Preconjugable<Domain, Domain, F> for Norm222<Domain, F> {
354 type Preconjugate<'a>
355 = Self
356 where
357 Self: 'a;
358
359 #[inline]
360 fn preconjugate(&self) -> Self::Preconjugate<'_> {
361 Self::new()
362 }
363 }
364
365 impl<Domain, F> Prox<Domain> for Norm222<Domain, F>
366 where
367 F: Float,
368 Domain: Euclidean<F, Output = Domain>,
369 {
370 type Prox<'a>
371 = Scaled<F>
372 where
373 Self: 'a;
374
375 fn prox_mapping(&self, τ: F) -> Self::Prox<'_> {
376 Scaled(F::ONE / (F::ONE + τ))
377 }
378 }

Mercurial Repository: alg_tools

mercurial