Mercurial > repos > alg_tools / file comparison

comparison: src/convex.rs

src/convex.rs

changeset 90
b3c35d16affe
parent 88
086a59b3a2b4
equal deleted inserted replaced
25:d14c877e14b7 90:b3c35d16affe
1 /*!
2 Some convex analysis basics
3 */
4
5 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 /// Trait for convex mappings. Has no features, just serves as a constraint
14 ///
15 /// TODO: should constrain `Mapping::Codomain` to implement a partial order,
16 /// but this makes everything complicated with little benefit.
17 pub trait ConvexMapping<Domain : Space, F : Num = f64> : Mapping<Domain, Codomain = F>
18 {}
19
20 /// Trait for mappings with a Fenchel conjugate
21 ///
22 /// The conjugate type has to implement [`ConvexMapping`], but a `Conjugable` mapping need
23 /// not be convex.
24 pub trait Conjugable<Domain : HasDual<F>, F : Num = f64> : Mapping<Domain> {
25 type Conjugate<'a> : ConvexMapping<Domain::DualSpace, F> where Self : 'a;
26
27 fn conjugate(&self) -> Self::Conjugate<'_>;
28 }
29
30 /// Trait for mappings with a Fenchel preconjugate
31 ///
32 /// In contrast to [`Conjugable`], the preconjugate need not implement [`ConvexMapping`],
33 /// but a `Preconjugable` mapping has to be convex.
34 pub trait Preconjugable<Domain, Predual, F : Num = f64> : ConvexMapping<Domain, F>
35 where
36 Domain : Space,
37 Predual : HasDual<F>
38 {
39 type Preconjugate<'a> : Mapping<Predual> where Self : 'a;
40
41 fn preconjugate(&self) -> Self::Preconjugate<'_>;
42 }
43
44 /// Trait for mappings with a proximap map
45 ///
46 /// The conjugate type has to implement [`ConvexMapping`], but a `Conjugable` mapping need
47 /// not be convex.
48 pub trait Prox<Domain : Space> : Mapping<Domain> {
49 type Prox<'a> : Mapping<Domain, Codomain=Domain> where Self : 'a;
50
51 /// Returns a proximal mapping with weight τ
52 fn prox_mapping(&self, τ : Self::Codomain) -> Self::Prox<'_>;
53
54 /// Calculate the proximal mapping with weight τ
55 fn prox<I : Instance<Domain>>(&self, τ : Self::Codomain, z : I) -> Domain {
56 self.prox_mapping(τ).apply(z)
57 }
58
59 /// Calculate the proximal mapping with weight τ in-place
60 fn prox_mut<'b>(&self, τ : Self::Codomain, y : &'b mut Domain)
61 where
62 &'b mut Domain : InstanceMut<Domain>,
63 Domain:: Decomp : DecompositionMut<Domain>,
64 for<'a> &'a Domain : Instance<Domain>,
65 {
66 *y = self.prox(τ, &*y);
67 }
68 }
69
70
71 /// Constraint to the unit ball of the norm described by `E`.
72 pub struct NormConstraint<F : Float, E : NormExponent> {
73 radius : F,
74 norm : NormMapping<F, E>,
75 }
76
77 impl<Domain, E, F> ConvexMapping<Domain, F> for NormMapping<F, E>
78 where
79 Domain : Space,
80 E : NormExponent,
81 F : Float,
82 Self : Mapping<Domain, Codomain=F>
83 {}
84
85
86 impl<F, E, Domain> Mapping<Domain> for NormConstraint<F, E>
87 where
88 Domain : Space + Norm<F, E>,
89 F : Float,
90 E : NormExponent,
91 {
92 type Codomain = F;
93
94 fn apply<I : Instance<Domain>>(&self, d : I) -> F {
95 if d.eval(|x| x.norm(self.norm.exponent)) <= self.radius {
96 F::ZERO
97 } else {
98 F::INFINITY
99 }
100 }
101 }
102
103 impl<Domain, E, F> ConvexMapping<Domain, F> for NormConstraint<F, E>
104 where
105 Domain : Space,
106 E : NormExponent,
107 F : Float,
108 Self : Mapping<Domain, Codomain=F>
109 {}
110
111
112 impl<E, F, Domain> Conjugable<Domain, F> for NormMapping<F, E>
113 where
114 E : HasDualExponent,
115 F : Float,
116 Domain : HasDual<F> + Norm<F, E> + Space,
117 <Domain as HasDual<F>>::DualSpace : Norm<F, E::DualExp>
118 {
119 type Conjugate<'a> = NormConstraint<F, E::DualExp> where Self : 'a;
120
121 fn conjugate(&self) -> Self::Conjugate<'_> {
122 NormConstraint {
123 radius : F::ONE,
124 norm : self.exponent.dual_exponent().as_mapping()
125 }
126 }
127 }
128
129 impl<C, E, F, Domain> Conjugable<Domain, F> for Weighted<NormMapping<F, E>, C>
130 where
131 C : Constant<Type = F>,
132 E : HasDualExponent,
133 F : Float,
134 Domain : HasDual<F> + Norm<F, E> + Space,
135 <Domain as HasDual<F>>::DualSpace : Norm<F, E::DualExp>
136 {
137 type Conjugate<'a> = NormConstraint<F, E::DualExp> where Self : 'a;
138
139 fn conjugate(&self) -> Self::Conjugate<'_> {
140 NormConstraint {
141 radius : self.weight.value(),
142 norm : self.base_fn.exponent.dual_exponent().as_mapping()
143 }
144 }
145 }
146
147 impl<Domain, E, F> Prox<Domain> for NormConstraint<F, E>
148 where
149 Domain : Space + Norm<F, E>,
150 E : NormExponent,
151 F : Float,
152 NormProjection<F, E> : Mapping<Domain, Codomain=Domain>,
153 {
154 type Prox<'a> = NormProjection<F, E> where Self : 'a;
155
156 #[inline]
157 fn prox_mapping(&self, _τ : Self::Codomain) -> Self::Prox<'_> {
158 assert!(self.radius >= F::ZERO);
159 NormProjection{ radius : self.radius, exponent : self.norm.exponent }
160 }
161 }
162
163 /// Projection to the unit ball of the norm described by `E`.
164 pub struct NormProjection<F : Float, E : NormExponent> {
165 radius : F,
166 exponent : E,
167 }
168
169 /*
170 impl<F, Domain> Mapping<Domain> for NormProjection<F, L2>
171 where
172 Domain : Space + Euclidean<F> + std::ops::MulAssign<F>,
173 F : Float,
174 {
175 type Codomain = Domain;
176
177 fn apply<I : Instance<Domain>>(&self, d : I) -> Domain {
178 d.own().proj_ball2(self.radius)
179 }
180 }
181 */
182
183 impl<F, E, Domain> Mapping<Domain> for NormProjection<F, E>
184 where
185 Domain : Space + Projection<F, E>,
186 F : Float,
187 E : NormExponent,
188 {
189 type Codomain = Domain;
190
191 fn apply<I : Instance<Domain>>(&self, d : I) -> Domain {
192 d.own().proj_ball(self.radius, self.exponent)
193 }
194 }
195
196
197 /// The zero mapping
198 pub struct Zero<Domain : Space, F : Num>(PhantomData<(Domain, F)>);
199
200 impl<Domain : Space, F : Num> Zero<Domain, F> {
201 pub fn new() -> Self {
202 Zero(PhantomData)
203 }
204 }
205
206 impl<Domain : Space, F : Num> Mapping<Domain> for Zero<Domain, F> {
207 type Codomain = F;
208
209 /// Compute the value of `self` at `x`.
210 fn apply<I : Instance<Domain>>(&self, _x : I) -> Self::Codomain {
211 F::ZERO
212 }
213 }
214
215 impl<Domain : Space, F : Num> ConvexMapping<Domain, F> for Zero<Domain, F> { }
216
217
218 impl<Domain : HasDual<F>, F : Float> Conjugable<Domain, F> for Zero<Domain, F> {
219 type Conjugate<'a> = ZeroIndicator<Domain::DualSpace, F> where Self : 'a;
220
221 #[inline]
222 fn conjugate(&self) -> Self::Conjugate<'_> {
223 ZeroIndicator::new()
224 }
225 }
226
227 impl<Domain, Predual, F : Float> Preconjugable<Domain, Predual, F> for Zero<Domain, F>
228 where
229 Domain : Space,
230 Predual : HasDual<F>
231 {
232 type Preconjugate<'a> = ZeroIndicator<Predual, F> where Self : 'a;
233
234 #[inline]
235 fn preconjugate(&self) -> Self::Preconjugate<'_> {
236 ZeroIndicator::new()
237 }
238 }
239
240 impl<Domain : Space + Clone, F : Num> Prox<Domain> for Zero<Domain, F> {
241 type Prox<'a> = IdOp<Domain> where Self : 'a;
242
243 #[inline]
244 fn prox_mapping(&self, _τ : Self::Codomain) -> Self::Prox<'_> {
245 IdOp::new()
246 }
247 }
248
249
250 /// The zero indicator
251 pub struct ZeroIndicator<Domain : Space, F : Num>(PhantomData<(Domain, F)>);
252
253 impl<Domain : Space, F : Num> ZeroIndicator<Domain, F> {
254 pub fn new() -> Self {
255 ZeroIndicator(PhantomData)
256 }
257 }
258
259 impl<Domain : Normed<F>, F : Float> Mapping<Domain> for ZeroIndicator<Domain, F> {
260 type Codomain = F;
261
262 /// Compute the value of `self` at `x`.
263 fn apply<I : Instance<Domain>>(&self, x : I) -> Self::Codomain {
264 x.eval(|x̃| if x̃.is_zero() { F::ZERO } else { F::INFINITY })
265 }
266 }
267
268 impl<Domain : Normed<F>, F : Float> ConvexMapping<Domain, F> for ZeroIndicator<Domain, F> { }
269
270 impl<Domain : HasDual<F>, F : Float> Conjugable<Domain, F> for ZeroIndicator<Domain, F> {
271 type Conjugate<'a> = Zero<Domain::DualSpace, F> where Self : 'a;
272
273 #[inline]
274 fn conjugate(&self) -> Self::Conjugate<'_> {
275 Zero::new()
276 }
277 }
278
279 impl<Domain, Predual, F : Float> Preconjugable<Domain, Predual, F> for ZeroIndicator<Domain, F>
280 where
281 Domain : Normed<F>,
282 Predual : HasDual<F>
283 {
284 type Preconjugate<'a> = Zero<Predual, F> where Self : 'a;
285
286 #[inline]
287 fn preconjugate(&self) -> Self::Preconjugate<'_> {
288 Zero::new()
289 }
290 }

Mercurial Repository: alg_tools

mercurial