diff -r e98e1da2530d -r e7f1cb4bec78 src/convex.rs --- a/src/convex.rs Mon Apr 28 23:16:56 2025 -0500 +++ b/src/convex.rs Tue Apr 29 00:03:12 2025 -0500 @@ -2,27 +2,29 @@ Some convex analysis basics */ -use std::marker::PhantomData; -use crate::types::*; +use crate::euclidean::Euclidean; +use crate::instance::{DecompositionMut, Instance, InstanceMut}; +use crate::linops::{IdOp, Scaled}; use crate::mapping::{Mapping, Space}; -use crate::linops::IdOp; -use crate::instance::{Instance, InstanceMut, DecompositionMut}; +use crate::norms::*; use crate::operator_arithmetic::{Constant, Weighted}; -use crate::norms::*; +use crate::types::*; +use std::marker::PhantomData; /// Trait for convex mappings. Has no features, just serves as a constraint /// /// TODO: should constrain `Mapping::Codomain` to implement a partial order, /// but this makes everything complicated with little benefit. -pub trait ConvexMapping : Mapping -{} +pub trait ConvexMapping: Mapping {} /// Trait for mappings with a Fenchel conjugate /// /// The conjugate type has to implement [`ConvexMapping`], but a `Conjugable` mapping need /// not be convex. -pub trait Conjugable, F : Num = f64> : Mapping { - type Conjugate<'a> : ConvexMapping where Self : 'a; +pub trait Conjugable, F: Num = f64>: Mapping { + type Conjugate<'a>: ConvexMapping + where + Self: 'a; fn conjugate(&self) -> Self::Conjugate<'_>; } @@ -31,12 +33,14 @@ /// /// In contrast to [`Conjugable`], the preconjugate need not implement [`ConvexMapping`], /// but a `Preconjugable` mapping has to be convex. -pub trait Preconjugable : ConvexMapping +pub trait Preconjugable: ConvexMapping where - Domain : Space, - Predual : HasDual + Domain: Space, + Predual: HasDual, { - type Preconjugate<'a> : Mapping where Self : 'a; + type Preconjugate<'a>: Mapping + where + Self: 'a; fn preconjugate(&self) -> Self::Preconjugate<'_>; } @@ -45,53 +49,54 @@ /// /// The conjugate type has to implement [`ConvexMapping`], but a `Conjugable` mapping need /// not be convex. -pub trait Prox : Mapping { - type Prox<'a> : Mapping where Self : 'a; +pub trait Prox: Mapping { + type Prox<'a>: Mapping + where + Self: 'a; /// Returns a proximal mapping with weight τ - fn prox_mapping(&self, τ : Self::Codomain) -> Self::Prox<'_>; + fn prox_mapping(&self, τ: Self::Codomain) -> Self::Prox<'_>; /// Calculate the proximal mapping with weight τ - fn prox>(&self, τ : Self::Codomain, z : I) -> Domain { + fn prox>(&self, τ: Self::Codomain, z: I) -> Domain { self.prox_mapping(τ).apply(z) } /// Calculate the proximal mapping with weight τ in-place - fn prox_mut<'b>(&self, τ : Self::Codomain, y : &'b mut Domain) + fn prox_mut<'b>(&self, τ: Self::Codomain, y: &'b mut Domain) where - &'b mut Domain : InstanceMut, - Domain:: Decomp : DecompositionMut, - for<'a> &'a Domain : Instance, + &'b mut Domain: InstanceMut, + Domain::Decomp: DecompositionMut, + for<'a> &'a Domain: Instance, { *y = self.prox(τ, &*y); } } - /// Constraint to the unit ball of the norm described by `E`. -pub struct NormConstraint { - radius : F, - norm : NormMapping, +pub struct NormConstraint { + radius: F, + norm: NormMapping, } impl ConvexMapping for NormMapping where - Domain : Space, - E : NormExponent, - F : Float, - Self : Mapping -{} - + Domain: Space, + E: NormExponent, + F: Float, + Self: Mapping, +{ +} impl Mapping for NormConstraint where - Domain : Space + Norm, - F : Float, - E : NormExponent, + Domain: Space + Norm, + F: Float, + E: NormExponent, { type Codomain = F; - fn apply>(&self, d : I) -> F { + fn apply>(&self, d: I) -> F { if d.eval(|x| x.norm(self.norm.exponent)) <= self.radius { F::ZERO } else { @@ -102,68 +107,80 @@ impl ConvexMapping for NormConstraint where - Domain : Space, - E : NormExponent, - F : Float, - Self : Mapping -{} - + Domain: Space, + E: NormExponent, + F: Float, + Self: Mapping, +{ +} impl Conjugable for NormMapping where - E : HasDualExponent, - F : Float, - Domain : HasDual + Norm + Space, - >::DualSpace : Norm + E: HasDualExponent, + F: Float, + Domain: HasDual + Norm + Space, + >::DualSpace: Norm, { - type Conjugate<'a> = NormConstraint where Self : 'a; + type Conjugate<'a> + = NormConstraint + where + Self: 'a; fn conjugate(&self) -> Self::Conjugate<'_> { NormConstraint { - radius : F::ONE, - norm : self.exponent.dual_exponent().as_mapping() + radius: F::ONE, + norm: self.exponent.dual_exponent().as_mapping(), } } } impl Conjugable for Weighted, C> where - C : Constant, - E : HasDualExponent, - F : Float, - Domain : HasDual + Norm + Space, - >::DualSpace : Norm + C: Constant, + E: HasDualExponent, + F: Float, + Domain: HasDual + Norm + Space, + >::DualSpace: Norm, { - type Conjugate<'a> = NormConstraint where Self : 'a; + type Conjugate<'a> + = NormConstraint + where + Self: 'a; fn conjugate(&self) -> Self::Conjugate<'_> { NormConstraint { - radius : self.weight.value(), - norm : self.base_fn.exponent.dual_exponent().as_mapping() + radius: self.weight.value(), + norm: self.base_fn.exponent.dual_exponent().as_mapping(), } } } impl Prox for NormConstraint where - Domain : Space + Norm, - E : NormExponent, - F : Float, - NormProjection : Mapping, + Domain: Space + Norm, + E: NormExponent, + F: Float, + NormProjection: Mapping, { - type Prox<'a> = NormProjection where Self : 'a; + type Prox<'a> + = NormProjection + where + Self: 'a; #[inline] - fn prox_mapping(&self, _τ : Self::Codomain) -> Self::Prox<'_> { + fn prox_mapping(&self, _τ: Self::Codomain) -> Self::Prox<'_> { assert!(self.radius >= F::ZERO); - NormProjection{ radius : self.radius, exponent : self.norm.exponent } + NormProjection { + radius: self.radius, + exponent: self.norm.exponent, + } } } /// Projection to the unit ball of the norm described by `E`. -pub struct NormProjection { - radius : F, - exponent : E, +pub struct NormProjection { + radius: F, + exponent: E, } /* @@ -182,41 +199,42 @@ impl Mapping for NormProjection where - Domain : Space + Projection, - F : Float, - E : NormExponent, + Domain: Space + Projection, + F: Float, + E: NormExponent, { type Codomain = Domain; - fn apply>(&self, d : I) -> Domain { + fn apply>(&self, d: I) -> Domain { d.own().proj_ball(self.radius, self.exponent) } } +/// The zero mapping +pub struct Zero(PhantomData<(Domain, F)>); -/// The zero mapping -pub struct Zero(PhantomData<(Domain, F)>); - -impl Zero { +impl Zero { pub fn new() -> Self { Zero(PhantomData) } } -impl Mapping for Zero { +impl Mapping for Zero { type Codomain = F; /// Compute the value of `self` at `x`. - fn apply>(&self, _x : I) -> Self::Codomain { + fn apply>(&self, _x: I) -> Self::Codomain { F::ZERO } } -impl ConvexMapping for Zero { } - +impl ConvexMapping for Zero {} -impl, F : Float> Conjugable for Zero { - type Conjugate<'a> = ZeroIndicator where Self : 'a; +impl, F: Float> Conjugable for Zero { + type Conjugate<'a> + = ZeroIndicator + where + Self: 'a; #[inline] fn conjugate(&self) -> Self::Conjugate<'_> { @@ -224,12 +242,15 @@ } } -impl Preconjugable for Zero +impl Preconjugable for Zero where - Domain : Space, - Predual : HasDual + Domain: Space, + Predual: HasDual, { - type Preconjugate<'a> = ZeroIndicator where Self : 'a; + type Preconjugate<'a> + = ZeroIndicator + where + Self: 'a; #[inline] fn preconjugate(&self) -> Self::Preconjugate<'_> { @@ -237,38 +258,43 @@ } } -impl Prox for Zero { - type Prox<'a> = IdOp where Self : 'a; +impl Prox for Zero { + type Prox<'a> + = IdOp + where + Self: 'a; #[inline] - fn prox_mapping(&self, _τ : Self::Codomain) -> Self::Prox<'_> { + fn prox_mapping(&self, _τ: Self::Codomain) -> Self::Prox<'_> { IdOp::new() } } +/// The zero indicator +pub struct ZeroIndicator(PhantomData<(Domain, F)>); -/// The zero indicator -pub struct ZeroIndicator(PhantomData<(Domain, F)>); - -impl ZeroIndicator { +impl ZeroIndicator { pub fn new() -> Self { ZeroIndicator(PhantomData) } } -impl, F : Float> Mapping for ZeroIndicator { +impl, F: Float> Mapping for ZeroIndicator { type Codomain = F; /// Compute the value of `self` at `x`. - fn apply>(&self, x : I) -> Self::Codomain { + fn apply>(&self, x: I) -> Self::Codomain { x.eval(|x̃| if x̃.is_zero() { F::ZERO } else { F::INFINITY }) } } -impl, F : Float> ConvexMapping for ZeroIndicator { } +impl, F: Float> ConvexMapping for ZeroIndicator {} -impl, F : Float> Conjugable for ZeroIndicator { - type Conjugate<'a> = Zero where Self : 'a; +impl, F: Float> Conjugable for ZeroIndicator { + type Conjugate<'a> + = Zero + where + Self: 'a; #[inline] fn conjugate(&self) -> Self::Conjugate<'_> { @@ -276,15 +302,77 @@ } } -impl Preconjugable for ZeroIndicator +impl Preconjugable for ZeroIndicator where - Domain : Normed, - Predual : HasDual + Domain: Normed, + Predual: HasDual, { - type Preconjugate<'a> = Zero where Self : 'a; + type Preconjugate<'a> + = Zero + where + Self: 'a; #[inline] fn preconjugate(&self) -> Self::Preconjugate<'_> { Zero::new() } } + +/// The squared Euclidean norm divided by two +pub struct Norm222(PhantomData<(Domain, F)>); + +impl, F: Float> Norm222 { + pub fn new() -> Self { + Norm222(PhantomData) + } +} + +impl, F: Float> Mapping for Norm222 { + type Codomain = F; + + /// Compute the value of `self` at `x`. + fn apply>(&self, x: I) -> Self::Codomain { + x.eval(|z| z.norm2_squared() / F::TWO) + } +} + +impl, F: Float> ConvexMapping for Norm222 {} + +impl, F: Float> Conjugable for Norm222 { + type Conjugate<'a> + = Self + where + Self: 'a; + + #[inline] + fn conjugate(&self) -> Self::Conjugate<'_> { + Self::new() + } +} + +impl, F: Float> Preconjugable for Norm222 { + type Preconjugate<'a> + = Self + where + Self: 'a; + + #[inline] + fn preconjugate(&self) -> Self::Preconjugate<'_> { + Self::new() + } +} + +impl Prox for Norm222 +where + F: Float, + Domain: Euclidean, +{ + type Prox<'a> + = Scaled + where + Self: 'a; + + fn prox_mapping(&self, τ: F) -> Self::Prox<'_> { + Scaled(F::ONE / (F::ONE + τ)) + } +}