src/mapping/quadratic.rs

branch
dev
changeset 99
9e5b9fc81c52
child 102
aead46a19767
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/mapping/quadratic.rs	Mon Apr 28 08:26:04 2025 -0500
@@ -0,0 +1,98 @@
+/*!
+Quadratic functions  of the form $\frac{1}{2}\|Ax-b\|_2^2$ for an operator $A$
+to a [`Euclidean`] space.
+*/
+
+#![allow(non_snake_case)]
+
+use super::{DifferentiableImpl, Differential, Lipschitz, Mapping};
+use crate::convex::ConvexMapping;
+use crate::euclidean::Euclidean;
+use crate::instance::{Instance, Space};
+use crate::linops::{BoundedLinear, Linear, Preadjointable};
+use crate::norms::{Norm, NormExponent, L2};
+use crate::types::Float;
+use std::marker::PhantomData;
+
+/// Functions of the form $\frac{1}{2}\|Ax-b\|_2^2$ for an operator $A$
+/// to a [`Euclidean`] space.
+#[derive(Clone, Copy)]
+pub struct Quadratic<'a, F: Float, Domain: Space, A: Mapping<Domain>> {
+    opA: &'a A,
+    b: &'a <A as Mapping<Domain>>::Codomain,
+    _phantoms: PhantomData<F>,
+}
+
+#[allow(non_snake_case)]
+impl<'a, F: Float, Domain: Space, A: Mapping<Domain>> Quadratic<'a, F, Domain, A> {
+    pub fn new(opA: &'a A, b: &'a A::Codomain) -> Self {
+        Quadratic {
+            opA,
+            b,
+            _phantoms: PhantomData,
+        }
+    }
+
+    pub fn operator(&self) -> &'a A {
+        self.opA
+    }
+
+    pub fn data(&self) -> &'a <A as Mapping<Domain>>::Codomain {
+        self.b
+    }
+}
+
+//+ AdjointProductBoundedBy<RNDM<F, N>, P, FloatType = F>,
+
+impl<'a, F: Float, X: Space, A: Mapping<X>> Mapping<X> for Quadratic<'a, F, X, A>
+where
+    A::Codomain: Euclidean<F>,
+{
+    type Codomain = F;
+
+    fn apply<I: Instance<X>>(&self, x: I) -> F {
+        // TODO: possibly (if at all more effcient) use GEMV once generalised
+        // to not require preallocation. However, Rust should be pretty efficient
+        // at not doing preallocations or anything here, as the result of self.opA.apply()
+        // can be consumed, so maybe GEMV is no use.
+        (self.opA.apply(x) - self.b).norm2_squared() / F::TWO
+    }
+}
+
+impl<'a, F: Float, X: Space, A: Linear<X>> ConvexMapping<X, F> for Quadratic<'a, F, X, A> where
+    A::Codomain: Euclidean<F>
+{
+}
+
+impl<'a, F, X, A> DifferentiableImpl<X> for Quadratic<'a, F, X, A>
+where
+    F: Float,
+    X: Space,
+    <A as Mapping<X>>::Codomain: Euclidean<F>,
+    A: Linear<X> + Preadjointable<X>,
+    <<A as Mapping<X>>::Codomain as Euclidean<F>>::Output: Instance<<A as Mapping<X>>::Codomain>,
+{
+    type Derivative = A::PreadjointCodomain;
+
+    fn differential_impl<I: Instance<X>>(&self, x: I) -> Self::Derivative {
+        // TODO: possibly (if at all more effcient) use GEMV once generalised
+        // to not require preallocation. However, Rust should be pretty efficient
+        // at not doing preallocations or anything here, as the result of self.opA.apply()
+        // can be consumed, so maybe GEMV is no use.
+        self.opA.preadjoint().apply(self.opA.apply(x) - self.b)
+    }
+}
+
+impl<'a, 'b, F, X, ExpX, A> Lipschitz<ExpX> for Differential<'b, X, Quadratic<'a, F, X, A>>
+where
+    F: Float,
+    X: Space + Clone + Norm<F, ExpX>,
+    ExpX: NormExponent,
+    A: Clone + BoundedLinear<X, ExpX, L2, F>,
+{
+    type FloatType = F;
+
+    fn lipschitz_factor(&self, seminorm: ExpX) -> Option<F> {
+        Some((*self.g).opA.opnorm_bound(seminorm, L2).powi(2))
+    }
+}

mercurial