--- a/src/kernels/gaussian.rs Fri Apr 28 13:15:19 2023 +0300
+++ b/src/kernels/gaussian.rs Tue Dec 31 09:34:24 2024 -0500
@@ -17,9 +17,10 @@
Weighted,
Bounded,
};
-use alg_tools::mapping::Apply;
+use alg_tools::mapping::{Apply, Differentiable};
use alg_tools::maputil::array_init;
+use crate::types::Lipschitz;
use crate::fourier::Fourier;
use super::base::*;
use super::ball_indicator::CubeIndicator;
@@ -75,14 +76,45 @@
impl<S, const N : usize> Apply<Loc<S::Type, N>> for Gaussian<S, N>
where S : Constant {
type Output = S::Type;
- // This is not normalised to neither to have value 1 at zero or integral 1
- // (unless the cut-off ε=0).
#[inline]
fn apply(&self, x : Loc<S::Type, N>) -> Self::Output {
self.apply(&x)
}
}
+#[replace_float_literals(S::Type::cast_from(literal))]
+impl<'a, S, const N : usize> Differentiable<&'a Loc<S::Type, N>> for Gaussian<S, N>
+where S : Constant {
+ type Output = Loc<S::Type, N>;
+ #[inline]
+ fn differential(&self, x : &'a Loc<S::Type, N>) -> Self::Output {
+ x * (self.apply(x) / self.variance.value())
+ }
+}
+
+impl<S, const N : usize> Differentiable<Loc<S::Type, N>> for Gaussian<S, N>
+where S : Constant {
+ type Output = Loc<S::Type, N>;
+ // This is not normalised to neither to have value 1 at zero or integral 1
+ // (unless the cut-off ε=0).
+ #[inline]
+ fn differential(&self, x : Loc<S::Type, N>) -> Self::Output {
+ x * (self.apply(&x) / self.variance.value())
+ }
+}
+
+#[replace_float_literals(S::Type::cast_from(literal))]
+impl<S, const N : usize> Lipschitz<L2> for Gaussian<S, N>
+where S : Constant {
+ type FloatType = S::Type;
+ fn lipschitz_factor(&self, L2 : L2) -> Option<Self::FloatType> {
+ // f(x)=f_1(‖x‖_2/σ) * √(2π) / √(2πσ)^N, where f_1 is one-dimensional Gaussian with
+ // variance 1. The Lipschitz factor of f_1 is e^{-1/2}/√(2π), see, e.g.,
+ // https://math.stackexchange.com/questions/3630967/is-the-gaussian-density-lipschitz-continuous
+ // Thus the Lipschitz factor we want is e^{-1/2} / (√(2πσ)^N * σ).
+ Some((-0.5).exp() / (self.scale() * self.variance.value().sqrt()))
+ }
+}
#[replace_float_literals(S::Type::cast_from(literal))]
impl<'a, S, const N : usize> Gaussian<S, N>
@@ -169,8 +201,8 @@
Gaussian<S, N>>;
-/// This implements $χ\_{[-b, b]^n} \* (f χ\_{[-a, a]^n})$
-/// where $a,b>0$ and $f$ is a gaussian kernel on $ℝ^n$.
+/// This implements $g := χ\_{[-b, b]^n} \* (f χ\_{[-a, a]^n})$ where $a,b>0$ and $f$ is
+/// a gaussian kernel on $ℝ^n$. For an expression for $g$, see Lemma 3.9 in the manuscript.
#[replace_float_literals(F::cast_from(literal))]
impl<'a, F : Float, R, C, S, const N : usize> Apply<&'a Loc<F, N>>
for Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>>
@@ -224,6 +256,89 @@
}
}
+/// This implements the differential of $g := χ\_{[-b, b]^n} \* (f χ\_{[-a, a]^n})$ where $a,b>0$
+/// and $f$ is a gaussian kernel on $ℝ^n$. For an expression for the value of $g$, from which the
+/// derivative readily arises (at points of differentiability), see Lemma 3.9 in the manuscript.
+#[replace_float_literals(F::cast_from(literal))]
+impl<'a, F : Float, R, C, S, const N : usize> Differentiable<&'a Loc<F, N>>
+for Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>>
+where R : Constant<Type=F>,
+ C : Constant<Type=F>,
+ S : Constant<Type=F> {
+
+ type Output = Loc<F, N>;
+
+ #[inline]
+ fn differential(&self, y : &'a Loc<F, N>) -> Loc<F, N> {
+ let Convolution(ref ind,
+ SupportProductFirst(ref cut,
+ ref gaussian)) = self;
+ let a = cut.r.value();
+ let b = ind.r.value();
+ let σ = gaussian.variance.value().sqrt();
+ let π = F::PI;
+ let t = F::SQRT_2 * σ;
+ let c = σ * (8.0/π).sqrt();
+ let cd = (8.0).sqrt(); // σ * (8.0/π).sqrt() / t * (√2/π)
+
+ // Calculate the values for all component functions of the
+ // product. This is just the loop from apply above.
+ let unscaled_vs = y.map(|x| {
+ let c1 = -(a.min(b + x)); //(-a).max(-x-b);
+ let c2 = a.min(b - x);
+ if c1 >= c2 {
+ 0.0
+ } else {
+ let e1 = F::cast_from(erf((c1 / t).as_()));
+ let e2 = F::cast_from(erf((c2 / t).as_()));
+ debug_assert!(e2 >= e1);
+ c * (e2 - e1)
+ }
+ });
+ // This computes the gradient for each coordinate
+ product_differential(y, &unscaled_vs, |x| {
+ let c1 = -(a.min(b + x)); //(-a).max(-x-b);
+ let c2 = a.min(b - x);
+ if c1 >= c2 {
+ 0.0
+ } else {
+ // erf'(z) = (2/√π)*exp(-z^2), and we get extra factor -1/√(2*σ) = -1/t
+ // from the chain rule
+ let de1 = (-(c1/t).powi(2)).exp();
+ let de2 = (-(c2/t).powi(2)).exp();
+ cd * (de1 - de2)
+ }
+ }) / gaussian.scale()
+ }
+}
+
+impl<F : Float, R, C, S, const N : usize> Differentiable<Loc<F, N>>
+for Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>>
+where R : Constant<Type=F>,
+ C : Constant<Type=F>,
+ S : Constant<Type=F> {
+
+ type Output = Loc<F, N>;
+
+ #[inline]
+ fn differential(&self, y : Loc<F, N>) -> Loc<F, N> {
+ self.differential(&y)
+ }
+}
+
+#[replace_float_literals(F::cast_from(literal))]
+impl<'a, F : Float, R, C, S, const N : usize> Lipschitz<L2>
+for Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>>
+where R : Constant<Type=F>,
+ C : Constant<Type=F>,
+ S : Constant<Type=F> {
+ type FloatType = F;
+
+ fn lipschitz_factor(&self, L2 : L2) -> Option<F> {
+ todo!("This requirement some error function work.")
+ }
+}
+
impl<F : Float, R, C, S, const N : usize>
Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>>
where R : Constant<Type=F>,