diff -r bd13c2ae3450 -r 56c8adc32b09 src/kernels/base.rs
--- a/src/kernels/base.rs	Fri Apr 28 13:15:19 2023 +0300
+++ b/src/kernels/base.rs	Tue Dec 31 09:34:24 2024 -0500
@@ -14,11 +14,12 @@
     GlobalAnalysis,
     Bounded,
 };
-use alg_tools::mapping::Apply;
-use alg_tools::maputil::{array_init, map2};
+use alg_tools::mapping::{Apply, Differentiable};
+use alg_tools::maputil::{array_init, map2, map1_indexed};
 use alg_tools::sets::SetOrd;
 
 use crate::fourier::Fourier;
+use crate::types::Lipschitz;
 
 /// Representation of the product of two kernels.
 ///
@@ -55,6 +56,33 @@
     }
 }
 
+impl<A, B, F : Float, const N : usize> Differentiable<Loc<F, N>>
+for SupportProductFirst<A, B>
+where A : for<'a> Apply<&'a Loc<F, N>, Output=F>
+          + for<'a> Differentiable<&'a Loc<F, N>, Output=Loc<F, N>>,
+      B : for<'a> Apply<&'a Loc<F, N>, Output=F>
+          + for<'a> Differentiable<&'a Loc<F, N>, Output=Loc<F, N>> {
+    type Output = Loc<F, N>;
+    #[inline]
+    fn differential(&self, x : Loc<F, N>) -> Self::Output {
+        self.0.differential(&x) * self.1.apply(&x) + self.1.differential(&x) * self.0.apply(&x)
+    }
+}
+
+impl<'a, A, B, F : Float, const N : usize> Differentiable<&'a Loc<F, N>>
+for SupportProductFirst<A, B>
+where A : Apply<&'a Loc<F, N>, Output=F>
+          + Differentiable<&'a Loc<F, N>, Output=Loc<F, N>>,
+      B : Apply<&'a Loc<F, N>, Output=F>
+          + Differentiable<&'a Loc<F, N>, Output=Loc<F, N>> {
+    type Output = Loc<F, N>;
+    #[inline]
+    fn differential(&self, x : &'a Loc<F, N>) -> Self::Output {
+        self.0.differential(&x) * self.1.apply(&x) + self.1.differential(&x) * self.0.apply(&x)
+    }
+}
+
+
 impl<'a, A, B, F : Float, const N : usize> Support<F, N>
 for SupportProductFirst<A, B>
 where A : Support<F, N>,
@@ -130,6 +158,28 @@
     }
 }
 
+impl<'a, A, B, F : Float, const N : usize> Differentiable<&'a Loc<F, N>>
+for SupportSum<A, B>
+where A : Differentiable<&'a Loc<F, N>, Output=Loc<F, N>>,
+      B : Differentiable<&'a Loc<F, N>, Output=Loc<F, N>> {
+    type Output = Loc<F, N>;
+    #[inline]
+    fn differential(&self, x : &'a Loc<F, N>) -> Self::Output {
+        self.0.differential(x) + self.1.differential(x)
+    }
+}
+
+impl<A, B, F : Float, const N : usize> Differentiable<Loc<F, N>>
+for SupportSum<A, B>
+where A : for<'a> Differentiable<&'a Loc<F, N>, Output=Loc<F, N>>,
+      B : for<'a> Differentiable<&'a Loc<F, N>, Output=Loc<F, N>> {
+    type Output = Loc<F, N>;
+    #[inline]
+    fn differential(&self, x : Loc<F, N>) -> Self::Output {
+        self.0.differential(&x) + self.1.differential(&x)
+    }
+}
+
 impl<'a, A, B, F : Float, const N : usize> Support<F, N>
 for SupportSum<A, B>
 where A : Support<F, N>,
@@ -174,6 +224,20 @@
     }
 }
 
+impl<F : Float, M : Copy, A, B> Lipschitz<M> for SupportSum<A, B>
+where A : Lipschitz<M, FloatType = F>,
+      B : Lipschitz<M, FloatType = F> {
+    type FloatType = F;
+
+    fn lipschitz_factor(&self, m : M) -> Option<F> {
+        match (self.0.lipschitz_factor(m), self.1.lipschitz_factor(m)) {
+            (Some(l0), Some(l1)) => Some(l0 + l1),
+            _ => None
+        }
+    }
+}
+
+
 /// Representation of the convolution of two kernels.
 ///
 /// The kernels typically implement [`Support`]s and [`Mapping`][alg_tools::mapping::Mapping].
@@ -187,6 +251,16 @@
     pub B
 );
 
+impl<F : Float, M, A, B> Lipschitz<M> for Convolution<A, B>
+where A : Bounded<F> ,
+      B : Lipschitz<M, FloatType = F> {
+    type FloatType = F;
+
+    fn lipschitz_factor(&self, m : M) -> Option<F> {
+        self.1.lipschitz_factor(m).map(|l| l * self.0.bounds().uniform())
+    }
+}
+
 /// Representation of the autoconvolution of a kernel.
 ///
 /// The kernel typically implements [`Support`] and [`Mapping`][alg_tools::mapping::Mapping].
@@ -198,6 +272,16 @@
     pub A
 );
 
+impl<F : Float, M, C> Lipschitz<M> for AutoConvolution<C>
+where C : Lipschitz<M, FloatType = F> + Bounded<F> {
+    type FloatType = F;
+
+    fn lipschitz_factor(&self, m : M) -> Option<F> {
+        self.0.lipschitz_factor(m).map(|l| l * self.0.bounds().uniform())
+    }
+}
+
+
 /// Representation a multi-dimensional product of a one-dimensional kernel.
 ///
 /// For $G: ℝ → ℝ$, this is the function $F(x\_1, …, x\_n) := \prod_{i=1}^n G(x\_i)$.
@@ -229,6 +313,45 @@
     }
 }
 
+impl<'a, G, F : Float, const N : usize> Differentiable<&'a Loc<F, N>>
+for UniformProduct<G, N>
+where G : Apply<Loc<F, 1>, Output=F> + Differentiable<Loc<F, 1>, Output=F> {
+    type Output = Loc<F, N>;
+    #[inline]
+    fn differential(&self, x : &'a Loc<F, N>) -> Loc<F, N> {
+        let vs = x.map(|y| self.0.apply(Loc([y])));
+        product_differential(x, &vs, |y| self.0.differential(Loc([y])))
+    }
+}
+
+/// Helper function to calulate the differential of $f(x)=∏_{i=1}^N g(x_i)$.
+///
+/// The vector `x` is the location, `vs` consists of the values `g(x_i)`, and
+/// `gd` calculates the derivative `g'`.
+#[inline]
+pub(crate) fn product_differential<F : Float, G : Fn(F) -> F, const N : usize>(
+    x : &Loc<F, N>,
+    vs : &Loc<F, N>,
+    gd : G
+) -> Loc<F, N> {
+    map1_indexed(x, |i, &y| {
+        gd(y) * vs.iter()
+                  .zip(0..)
+                  .filter_map(|(v, j)| (j != i).then_some(*v))
+                  .product()
+    }).into()
+}
+
+impl<G, F : Float, const N : usize> Differentiable<Loc<F, N>>
+for UniformProduct<G, N>
+where G : Apply<Loc<F, 1>, Output=F> + Differentiable<Loc<F, 1>, Output=F> {
+    type Output = Loc<F, N>;
+    #[inline]
+    fn differential(&self, x : Loc<F, N>) -> Loc<F, N> {
+        self.differential(&x)
+    }
+}
+
 impl<G, F : Float, const N : usize> Support<F, N>
 for UniformProduct<G, N>
 where G : Support<F, 1> {