Generalizing Stochastic Smoothing for Differentiation and Gradient Estimation

## Abstract

Abstract We deal with the problem of gradient estimation for stochastic differentiable relaxations of algorithms, operators, simulators, and other non-differentiable functions.
Stochastic smoothing conventionally perturbs the input of a non-differentiable function with a differentiable density distribution with full support, smoothing it and enabling gradient estimation.
Our theory starts at first principles to derive stochastic smoothing with reduced assumptions, without requiring a differentiable density nor full support, and we present a general framework for relaxation and gradient estimation of non-differentiable black-box functions f : ℝ n → ℝ m f:\mathbb{R}^{n}\to\mathbb{R}^{m} .
We develop variance reduction for gradient estimation from 3 orthogonal perspectives.
Empirically, we benchmark 6 distributions and up to 24 variance reduction strategies for differentiable sorting and ranking, differentiable shortest-paths on graphs, differentiable rendering for pose estimation, as well as differentiable cryo-ET simulations.

## 1 Introduction

The differentiation of algorithms, operators, and other non-differentiable functions has been a topic of rapidly increasing interest in the machine learning community [, , , , , , ].
In particular, whenever we want to integrate a non-differentiable operation (such as ranking) into a machine learning pipeline, we need to relax it into a differentiable form in order to allow for backpropagation.
To give a concrete example, a body of recent work considered continuously relaxing the sorting and ranking operators for tasks like learning-to-rank [, , , , , , , , , ].
These works can be categorized into either casting sorting and ranking as a related problem (e.g., optimal transport []) and differentiably relaxing it (e.g., via entropy-regularized OT []) or by considering a sorting algorithm and continuously relaxing it on the level of individual operations or program statements [, , , ].
To give another example, in the space of differentiable graph algorithms and clustering, popular directions either relax algorithms on a statement level [] or cast the algorithm as a convex optimization problem and differentiate the solution of the optimization problem under perturbed parameterization [, , ].

Complementary to these directions of research, in this work, we consider algorithms, operators, simulators, and other non-differentiable functions directly as black-box functions and differentiably relax them via stochastic smoothing [], i.e., via stochastic perturbations of the inputs and via multiple function evaluation.
This is challenging as, so far, gradient estimators have come with large variance and supported only a restrictive set of smoothing distributions.
More concretely, for a black-box function $f:\mathbb{R}^{n}\to\mathbb{R}^{m}$, we consider the problem of estimating the derivative (or gradient) of the relaxation

$$ $f_{\epsilon}(x)=\mathbb{E}_{\epsilon\sim\mu}\big[f(x+\epsilon)\big]={\int}f(x+\epsilon)\,\mu(\epsilon)\,\mathrm{d}\epsilon$ (1) $$

where $\epsilon$ is a sample from a probability distribution with an (absolutely) continuous
density $\mu(\epsilon)$.
$f_{\epsilon}$ is a differentiable function (regardless of differentiability properties of $f$ itself, see Section 2 for details) and its gradient is well defined.
Under limiting restrictions on the probability distribution $\mu$ used for smoothing, gradient estimators exist in the literature [, , , ].

The contribution of this work lies in providing more generalized gradient estimators (reducing assumptions on $\mu$, Lemma 3) that exhibit reduced variances (Sec. 2.2) for the application of differentiably relaxing conventionally non-differentiable algorithms.
Moreover, we enable smoothing with and differentiation wrt. non-diagonal covariances (Thm. 7), characterize formal requirements for $f$ (Lem. 9+10), discriminate smoothing of algorithms and losses (Sec. 2.3), and provide a $k$-sample median extension (Apx. C).
The proposed approach is applicable for differentiation of (i) arbitrary(^1^1footnotemark: 1) functions, which are (ii) considered as a black-box, which (iii) can be called many times (primarily) at low cost, and (iv) should be smoothed with any distribution with absolutely continuous density on $\mathbb{R}$.
This contrasts prior work, which smoothed (i) convex optimizers [, ], (ii) used first-order gradients [, , , ], (iii) allowed calling an environment only once or few times in RL [], and/or (iv) smoothed with fully supported differentiable density distributions [, , , ].
In machine learning, many other subfields also utilize the ideas underlying stochastic smoothing; stochastic smoothing and similar methods can be found, e.g., in REINFORCE [], the score function estimator [], the CatLog-Derivative trick [], perturbed optimizers [, ], among others.

## 2 Differentiation via Stochastic Smoothing

We begin by recapitulating the core of the stochastic smoothing method.
The idea behind smoothing is that given a (potentially non-differentiable) function $f:\mathbb{R}^{n}\to\mathbb{R}^{m}$ (^1^11Traditionally and formally, only functions with compact range have been considered for $f$ (i.e., $f:\mathbb{R}^{n}\to[a,b]$) [, ]. More recently, i.a., Abernethy et al. [] have considered the general case of function with the real range. While this is very helpful, e.g., enabling linear functions for $f$, this is (even without our generalizations) not always finitely defined as we discuss with the help of degenerate examples in Appendix B. There, we characterize the set of valid $f$ leading to finitely defined $f_{\epsilon}$ in Lemma 9 as well as $\nabla f_{\epsilon}$ in Lemma 10 in dependence on $\mu$. We remark that, beyond discussions in Appendix B, we assume this to be satisfied by $f$.), we can relax the function to a differentiable function by perturbing its argument with a probability distribution:
if $\epsilon\in\mathbb{R}^{n}$ follows a distribution with a differentiable density $\mu(\epsilon)$, then $f_{\epsilon}(x)=\mathbb{E}_{\epsilon}[f(x+\epsilon)]$ is differentiable.

For the case of $m=1$, i.e., for a scalar function $f$, we can compute the gradient of $f_{\epsilon}$ by following and extending part of Lemma 1.5 in Abernethy et al. [] as follows:

###### Lemma 1 (Differentiable Density Smoothing) .

Given a function $f:\mathbb{R}^{n}\to\mathbb{R}$ (^1^1footnotemark: 1) and a differentiable probability density function $\mu(\epsilon)$ with full support on $\mathbb{R}^{n}$, then $f_{\epsilon}$ is differentiable and

$$ $\nabla_{\!x}f_{\epsilon}(x)=\nabla_{\!x}\mathbb{E}_{\epsilon\sim\mu}\big[f(x+\epsilon)\big]=\mathbb{E}_{\epsilon\sim\mu}\big[f(x+\epsilon)\cdot\nabla_{\!\epsilon}{-}\log\mu(\epsilon)\big]\,.$ (2) $$

###### Proof.

For didactic reasons, we include a full proof in the paper to support the reader’s understanding of the core of the method.
Via a change of variables, replacing $x+\epsilon$ by $u$, we obtain ($d\epsilon/du=1$)

$$ $\displaystyle f_{\epsilon}(x)=\int f(x+\epsilon)\mu(\epsilon)\,d\epsilon=\int f(u)\mu(u-x)\,du\,.$ (3) $$

Now,

$$ $\displaystyle\nabla_{\!x}f_{\epsilon}(x)=\nabla_{\!x}\int f(u)\mu(u-x)\,du=\int f(u)\,\nabla_{\!x}\mu(u-x)\,du\,.$ (4) $$

Using $\nabla_{\!x}\mu(u-x)=-\nabla_{\!\epsilon}\mu(\epsilon)$,

$$ $\kern 112.00014pt\nabla_{\!x}f_{\epsilon}(x)=-\int f(x+\epsilon)\,\nabla_{\!\epsilon}\mu(\epsilon)\,d\epsilon\,.$ (5) $$

Because $\frac{\partial\mu(\epsilon)}{\partial\epsilon}=\mu(\epsilon)\cdot\frac{\partial\log\mu(\epsilon)}{\partial\epsilon}\,$,
we can simplify the expression to

$$ $\displaystyle\nabla_{\!x}f_{\epsilon}(x)=-\int f(x+\epsilon)\,\mu(\epsilon)\,\nabla_{\!\epsilon}\log\mu(\epsilon)\,d\epsilon\,=\,\mathbb{E}_{\epsilon\sim\mu}\,\Big[f(x+\epsilon)\,\nabla_{\!\epsilon}{-}\log\mu(\epsilon)\Big]\,.$ (6) $$

∎

Empirically, for a number of samples $s$, this gradient estimator can be evaluated without bias via

$$ $\displaystyle\nabla_{\!x}f_{\epsilon}(x)\,\triangleq\,\frac{1}{s}\sum_{i=1}^{s}\,\Big[f(x+\epsilon_{i})\,\nabla_{\!\epsilon_{i}}{-}\log\mu(\epsilon_{i})\Big]\qquad\epsilon_{1},...,\epsilon_{s}\sim\mu\,.$ (7) $$

###### Corollary 2 (Differentiable Density Smoothing for Vector-valued Functions) .

We can extend Lemma 1 to vector-valued functions $f:\mathbb{R}^{n}\to\mathbb{R}^{m}$, allowing to compute the Jacobian matrix $\mathbf{J}_{f_{\epsilon}}\!\!\in\mathbb{R}^{m\times n}$ as

$$ $\mathbf{J}_{f_{\epsilon}}(x)=\mathbb{E}_{\epsilon\sim\mu}\Big[f(x+\epsilon)\cdot\big(\nabla_{\!\epsilon}{-}\log\mu(\epsilon)\big)^{\top}\Big]\,.$ (8) $$

We remark that prior work (e.g., []) limits $\mu$ to be a differentiable density with full support on $\mathbb{R}$, typically of exponential family, whereas we generalize it to any absolutely continuous density, and include additional generalizations.
This has important implications for distributions such as Cauchy, Laplace, and Triangular, which we show to have considerable practical relevance.

###### Lemma 3 (Requirement of Continuity of $\mu$ ) .

If $\mu(\epsilon)$ is absolutely continuous (and not necessarily differentiable), then $f_{\epsilon}$ is continuous and differentiable everywhere.

$$ $\nabla_{\!x}f_{\epsilon}(x)=\mathbb{E}_{\epsilon\sim\mu}\,\Big[f(x+\epsilon)\,\cdot\mathbf{1}_{\epsilon\notin\Omega}\cdot\,\nabla_{\!\epsilon}{-}\log\mu(\epsilon)\Big]\,.$ (9) $$

$\Omega$ is the zero-measure set of points with undefined gradient. We provide the proof in Appendix A.1.

Lemma 3 has important implications.
In particular, it enables smoothing with non-differentiable density distributions such as the Laplace distribution, the triangular distribution, and the Wigner Semicircle distribution [, ] while maintaining differentiability of $f_{\epsilon}$.

###### Remark 4 (Requirement of Continuity of $\mu$ ) .

However, it is crucial to mention that, for stochastic smoothing (Lemmas 1, 3, Corollary 2), $\mu$ has to be continuous.
For example, the uniform distribution is not a valid choice because it
does not have a continuous density on $\mathbb{R}$. ($\mathcal{U}(a,b)$ has discontinuities at $a,b$ where it jumps between 0 and $1/(b-a)$.)
With other formulations, e.g., [, ], it is possible to perform smoothing with a uniform distribution over a ball; however, if $f$ is discontinuous, uniform smoothing may not lead to a differentiable function.
Continuity is a requirement but not a sufficient condition, and absolutely continuous is a sufficient condition; however, the difference to continuity corresponds only to non-practical and adversarial examples, e.g., the Cantor or Weierstrass functions.

###### Remark 5 (Gaussian Smoothing) .

A popular special case of differentiable stochastic smoothing is smoothing with a Gaussian distribution $\mu_{\mathcal{N}}=N(\mathbf{0}_{n},\mathbf{I}_{n})$.
Here, due to the nature of the probability density function of a Gaussian, $\nabla_{\!\epsilon}{-}\log\mu_{\mathcal{N}}(\epsilon)=\epsilon$.
Further, when $\mu_{\mathcal{N}_{\sigma}}\!\!=N(\mathbf{0}_{n},\sigma^{2}\mathbf{I}_{n})$, then $\nabla_{\!\epsilon}{-}\log\mu_{\mathcal{N}_{\sigma}}(\epsilon)=\epsilon/\sigma$.
We emphasize that this equality only holds for the Gaussian distribution.

Equipped with the core idea behind stochastic smoothing, we can differentiate any function $f$ via perturbation with a probability distribution with (absolutely) continuous density.

Typically, probability distributions that we consider for smoothing are parameterized via a scale parameter, vi
7
., the standard deviation $\sigma$ in a Gaussian distribution or the scale $\gamma$ in a Cauchy distribution.
Extending the formalism above, we may be interested in differentiating with respect to the scale parameter $\gamma$ of our distribution $\mu$.
This becomes especially attractive when optimizing the scale and, thereby, degree of relaxation of our probability distribution.
While our formalism allows reparameterization to express $\gamma$ within $\mu$, we can also explicitly write it as

$$ $\nabla_{\!x}f_{\gamma\epsilon}(x)=\nabla_{\!x}\mathbb{E}_{\epsilon\sim\mu}\big[f(x+\gamma\cdot\epsilon)\big]=\mathbb{E}_{\epsilon\sim\mu}\big[f(x+\gamma\cdot\epsilon)\cdot\big(\nabla_{\!\epsilon}{-}\log\mu(\epsilon)\big)\,/\,\gamma\Big]\,.$ (10) $$

Now, we can differentiate wrt. $\gamma$, i.e., we can compute $\nabla_{\!\gamma}f_{\gamma\epsilon}(x)$.

###### Lemma 6 (Differentiation wrt. $\gamma$ ) .

Extending Lemma 1, Corollary 2, and Lemma 3, we
have

$$ $\nabla_{\!\gamma}f_{\gamma\epsilon}(x)=\nabla_{\!\gamma}\mathbb{E}_{\epsilon\sim\mu}\Big[f(x+\gamma\cdot\epsilon)\big]=\mathbb{E}_{\epsilon\sim\mu}\big[f(x+\gamma\cdot\epsilon)\cdot\big(-1+(\nabla_{\!\epsilon}{-}\log\mu(\epsilon))^{\top}\cdot\epsilon\big)\,/\,\gamma\Big]\,.$ (11) $$

The proof is deferred to Appendix A.2.

We can extend $\gamma$ for multivariate distributions to a scale matrix $\mathbf{\Sigma}/\mathbf{L}$ (e.g., a covariance matrix).

###### Theorem 7 (Multivariate Smoothing with Covariance Matrix) .

We have a function $f:\mathbb{R}^{n}\to\mathbb{R}^{m}$.
We assume $\epsilon$ is drawn from a multivariate distribution with absolutely continuous density in $\mathbb{R}^{n}$.
We have an invertible scale matrix $\mathbf{L}\in\mathbb{R}^{n\times n}$ (e.g., for a covariance matrix $\mathbf{\Sigma}$, $\mathbf{L}$ is based on its Cholesky decomposition $\mathbf{L}\mathbf{L}^{\top}=\mathbf{\Sigma}$).
We define $f_{\mathbf{L}\epsilon}(x)=\mathbb{E}_{\epsilon\sim\mu}\big[f(x+\mathbf{L}\cdot\epsilon)\big]\,$.
Then, our derivatives $\partial f_{\mathbf{L}\epsilon}(x)\,/\,\partial x~~(\in\mathbb{R}^{m\times n})$ and $\partial f_{\mathbf{L}\epsilon}(x)\,/\,\partial\mathbf{L}~~(\in\mathbb{R}^{m\times n\times n})$ can be computed as

$$ $\displaystyle\nabla_{\!x}\big(f_{\mathbf{L}\epsilon}(x)\big)_{i}$ $\displaystyle=\mathbb{E}_{\epsilon\sim\mu}\Big[f(x+\mathbf{L}\cdot\epsilon)_{i}\cdot\mathbf{L}^{-1}\cdot\big(\nabla_{\!\epsilon}{-}\log\mu(\epsilon)\big)\Big]\,,$ (12) $\displaystyle\nabla_{\!\mathbf{L}}\big(f_{\mathbf{L}\epsilon}(x)\big)_{i}$ $\displaystyle=\mathbb{E}_{\epsilon\sim\mu}\Big[f(x+\mathbf{L}\cdot\epsilon)_{i}\cdot\mathbf{L}^{-\top}\cdot\big(-1+(\nabla_{\!\epsilon}{-}\log\mu(\epsilon))\cdot\epsilon^{\top}\big)\Big]\,.$ (13) $$

Above, the indicator (from (9)) is omitted for a simplified exposition.
Proofs are deferred to Apx. A.3.

Before we continue with examples for distributions, we discuss two extensions, vi
7
. differentiating output covariances of smoothing, and differentiating the expected $k$-sample median.

For uncertainty quantification (UQ) applications, e.g., in the context of propagating distributions [], we may further be interested in computing the derivative of the output covariance matrix, as illustrated by the following theorem.

###### Theorem 8 (Output Covariance of Multivariate Smoothing for UQ) .

Given the assumptions of Theorem 7, we may also be interested in computing the output covariance $G_{\mathbf{L}\epsilon}$:

$$ $G_{\mathbf{L}\epsilon}(x)=\operatorname{Cov}_{\epsilon\sim\mu}\big[f(x+\mathbf{L}\cdot\epsilon)\big]\,.$ (14) $$

We can compute the derivative of the output covariance wrt. the input $\partial G_{\mathbf{L}\epsilon}(x)\,/\,\partial x(\in\mathbb{R}^{m\times m\times n})$ as

$$ $\displaystyle\nabla_{x}\,G_{\mathbf{L}\epsilon}(x)_{i,j}=$ $\displaystyle\mathbb{E}_{\epsilon\sim\mu}\Big[f(x+\mathbf{L}\epsilon)_{i}\cdot f(x+\mathbf{L}\epsilon)_{j}\,\cdot\mathbf{L}^{-1}\cdot\nabla_{\epsilon}{-}\log\mu(\epsilon)\Big]$ (15) $\displaystyle-f_{\mathbf{L}\epsilon}(x)_{i}\cdot\nabla_{x}\,f_{\mathbf{L}\epsilon}(x)_{j}~-f_{\mathbf{L}\epsilon}(x)_{j}\cdot\nabla_{x}\,f_{\mathbf{L}\epsilon}(x)_{i}$ $$

Further, we can compute the derivative wrt. the input scale matrix $\partial G_{\mathbf{L}\epsilon}(x)\,/\,\partial\mathbf{L}(\in\mathbb{R}^{m\times m\times n\times n})$ as

$$ $\displaystyle\nabla_{\mathbf{L}}\,G_{\mathbf{L}\epsilon}(x)_{i,j}=$ $\displaystyle-\mathbb{E}_{\epsilon\sim\mu}\Big[f(x+\mathbf{L}\epsilon)_{i}\cdot f(x+\mathbf{L}\epsilon)_{j}\cdot\big(\mathbf{L}^{-\top}\cdot\nabla_{\epsilon}\mu(\epsilon)\cdot\epsilon^{\top}\,/\mu(\epsilon)+\mathbf{L}^{-\top}\big)\Big]$ (16) $\displaystyle-f_{\mathbf{L}\epsilon}(x)_{i}\cdot\nabla_{\mathbf{L}}\,f_{\mathbf{L}\epsilon}(x)_{j}~-f_{\mathbf{L}\epsilon}(x)_{j}\cdot\nabla_{\mathbf{L}}\,f_{\mathbf{L}\epsilon}(x)_{i}$ $$

The proofs are deferred to Appendix A.4.

In Appendix C, we additionally extend stochastic smoothing to differentiating the expected $k$-sample median, show that it is differentiable, and provide an unbiased gradient estimator in Lemma 13.

###### Lemma 1 (Differentiable Density Smoothing) .

###### Proof.

###### Corollary 2 (Differentiable Density Smoothing for Vector-valued Functions) .

###### Lemma 3 (Requirement of Continuity of μ \mu ) .

###### Remark 4 (Requirement of Continuity of μ \mu ) .

###### Remark 5 (Gaussian Smoothing) .

###### Lemma 6 (Differentiation wrt. γ \gamma ) .

###### Theorem 7 (Multivariate Smoothing with Covariance Matrix) .

###### Theorem 8 (Output Covariance of Multivariate Smoothing for UQ) .

### 2.1 Distribution Examples

After covering the underlying theory of generalized stochastic smoothing, in this section, we provide examples of specific distributions that our theory applies to.
We illustrate the distributions in Table 1.

Before delving into individual choices for distributions, we provide a clarification for multivariate densities $\mu:\mathbb{R}^{n}\to\mathbb{R}_{\geq 0}$:
We consider the $n$-dimensional multivariate form of a distribution as the concatenation of $n$ independent univariate distributions.
Thus, for $\mu_{1}$ as the univariate formulation of the density, we have the proportionality $\mu(\epsilon)\simeq\prod_{i=1}^{n}\mu_{1}(\epsilon_{i})$.
We remark that the distribution by which we smooth ($\mathbf{L}\epsilon$) is not an isotropic (per-dimension independent) distribution.
Instead, through transformation by the scale matrix $\mathbf{L}$, e.g., in the case of the Gaussian distribution, $\mathbf{L}\epsilon$ covers the entire space of multivariate Gaussian distributions with arbitrary covariance matrices.

Beyond the Gaussian distribution, the logistic distribution offers heavier tails, and the Gumbel distribution provides max-stability, which can be important for specific tasks.
The Cauchy distribution [], with its undefined mean and infinite variance, also has important implications in smoothing: e.g., the Cauchy distribution is shown to provide monotonicity in differentiable sorting networks [].
While prior art [] heuristically utilized the Cauchy distribution for stochastic smoothing of argmax, this had been, thus far, without a general formal justification.

In this work, for the first time, we consider Laplace and triangular distributions.
First, the Laplace distribution, as the symmetric extension of the exponential distribution, does not lie in the space of exponential family distributions, and is not differentiable at 0.
Via Lemma 3, we show that stochastic smoothing can still be applied and is exactly correct despite non-differentiablity of the distribution.
A benefit of the Laplace distribution is that all samples contribute equally to the gradient computation ($|\nabla_{\epsilon}{-}\log\mu(\epsilon)|=1$), reducing variance.
Second, with the triangular distribution, we illustrate, for the first time, that stochastic smoothing can be performed even with a non-differentiable distribution with compact support ($[-1,1]$).
This is crucial if the domain of $f$ has to be limited to a compact set rather than the real domain, or in applications where smoothing beyond a limited distance to the original point is not meaningful.
A hypothetical application for this could be differentiating a physical motor controlled robot in reinforcement learning where we may not want to support an infinite range for safety considerations.

**Table 1: Probability distributions considered for generalized stochastic smoothing. Displayed is (from left to right) the density of the distribution $\mu(\epsilon)$ (plot + equation), the derivative of the NLL (equation), and the product between the density and the derivative of the NLL (plot). The latter plot corresponds to the kernel that $f$ is effectively convolved by to estimate the gradient. $(*)$: applies to $\epsilon\in(-1,1)\setminus\{0\}$, otherwise 0 or undefined.**
| Distribution | Density / PDF  $\mu(\epsilon)$ | $\nabla_{\!\epsilon}{-}\log\mu(\epsilon)$ | $\kern-0.50003pt\mu(\epsilon)\cdot\nabla_{\!\epsilon}{-}\log\mu(\epsilon)\kern-1.99997pt$ |  |
| --- | --- | --- | --- | --- |
| Gaussian |  | $\displaystyle\frac{1}{\sqrt{2\pi}}\exp\big(-\nicefrac{{1}}{{2}}\cdot\epsilon^{2}\big)$ | $\epsilon$ |  |
| Logistic |  | $\displaystyle\frac{\exp(-\epsilon)}{(1+\exp(-\epsilon))^{2}}$ | $\tanh(\epsilon/2)$ |  |
| Gumbel |  | $\exp(-\epsilon-\exp(-\epsilon))$ | $1-\exp(-\epsilon)$ |  |
| Cauchy |  | $\displaystyle\frac{1}{\pi\cdot(1+\epsilon^{2})}$ | $\displaystyle\frac{2\cdot\epsilon}{1+\epsilon^{2}}$ |  |
| Laplace |  | $\nicefrac{{1}}{{2}}\cdot\exp(-|\epsilon|)$ | $\operatorname{sign}(\epsilon)$ |  |
| Triangular |  | $\max(0,1-|\epsilon|)$ | $\displaystyle\frac{\operatorname{sign}(\epsilon)}{1-|\epsilon|}~^{(*)}$ |  |

### 2.2 Variance Reduction

Given an unbiased estimator of the gradient, e.g., in its simplest form (2), we desire reducing its variance, or, in other words, improve the quality of the gradient estimate for a given number of samples.
For this, we consider 3 orthogonal perspectives of variance reduction:
covariates, antithetic samples, and (randomized) quasi-Monte Carlo.

Figure: Figure 1: Comparison of covariates: a non-differentiable function (dark blue) is smoothed with a logistic distribution (light blue). The original gradient (dark red) is not everywhere defined, and does not meaningfully represent the gradient. The gradient of the smoothed function is shown in pink. Grey illustrates the variance of a gradient estimate with $5$ samples via the $[25\%,75\%]$ (dark grey) and $[10\%,90\%]$ (light grey) percentiles. Using $f(x)$ as a covariate, instead of using none reduces the gradient variance, in particular whenever $f(x)$ is large. Leave-one-out (LOO) further improves over $f(x)$ at discontinuities of the original function $f$ (i.e., at $x{=}1$), but has slightly higher variance than $f(x)$ where $f$ is continuous and has large values (i.e., at $x{=}{-}2$.)
Refer to caption: x13.png

To illustratively derive the first two variance reductions, let us consider the case of smoothing a constant function $f(x)=v$ for some large constant $v\gg 0$.
Naturally, $f_{\epsilon}(x)=v$ and $\nabla_{\!x}f_{\epsilon}(x)=0$.
However, for a finite number of samples $s$, our empirical estimate (e.g., (7)) will differ from 0 almost surely.
As the gradient of $f_{\epsilon}(x)-c$ wrt. $x$ does not depend on $c$, we have $\nabla_{\!x}f_{\epsilon}(x)=\nabla_{\!x}(f_{\epsilon}(x)-c)$.
If we choose $c=v$, the variance of the gradient estimator is reduced to 0.
For general and non-constant $f$, we can estimate the optimal choice of $c$ via $c=f(x)$ or via the leave-one-out estimator [, ] of $f_{\epsilon}(x)$.
In the fields of stochastic smoothing of optimizers and reinforcement learning this is known as the method of covariates.
$f(x)$ and LOO were previously considered for smoothing, e.g., in [] and [], respectively.
We illustrate the effects of both choices of covariates in Figure 1.

From an orthogonal perspective, we observe that $\mathbb{E}_{\epsilon\sim\mu}\!\big[\nabla_{\!\epsilon}{-}\!\log\mu(\epsilon))\big]{=}\,0$, which follows, e.g., from $\nabla_{\!x}f_{\epsilon}(x)=0\,{=}\,\mathbb{E}_{\epsilon\sim\mu}\!\big[v{\cdot}\nabla_{\!\epsilon}{-}\!\log\mu(\epsilon)\big]$.
For symmetric distributions, we can guarantee an empirical estimate to be 0 by always using pairs of antithetic samples [], i.e., complementary $\epsilon$s.
Using $\epsilon^{\prime}=-\epsilon$, we have $\nabla_{\!\epsilon}\log\mu(\epsilon)+\nabla_{\!\epsilon^{\prime}}\log\mu({\epsilon^{\prime}})=0$.
This is illustrated in Figure 2 (2).
In our experiments in the next section, we observe antithetic sampling to generally perform poorly in comparison to other variance reduction techniques.

Figure: Figure 2: Sampling strategies. Left to right: Monte-Carlo (MC), Antithetic Monte-Carlo, Cartesian Quasi-Monte-Carlo (QMC), Cartesian Randomized-Quasi-Monte-Carlo (RQMC), Latin-Hypercube Sampled QMC and RQMC. Samples can be transformed via the inverse CDF of a respective distribution.
Refer to caption: x16.png

A third perspective considers that points sampled with standard Monte Carlo (MC) methods (see Fig. 2 (1)), due to the random nature of the sampling, often form (accidental) clumps while other areas are void of samples.
To counteract this, quasi-Monte Carlo (QMC) methods [] spread out sampled points as evenly as possible by foregoing randomness and choosing points from a regular grid, e.g., a simple Cartesian grid, taking the grid cell centers as samples (see Fig. 2 (3)).
Via the inverse CDF of the respective distribution, the points can be mapped from the unit hypercube to samples from a respective distribution.
However, discarding randomness makes the sampling process deterministic, limits the dispersion introduced by the smoothing distribution to concrete points, and hence makes the estimator biased.
Randomized quasi-Monte Carlo (RQMC) [] methods overcome this difficulty by reintroducing some randomness.
Like QMC, RQMC uses a grid to subdivide $[0,1]^{n}$ into cells, but then samples a point from each cell (see Fig. 2 (4)) instead of taking the grid cell center.
While regular MC sampling leads to variances of $\mathcal{O}(1/\sqrt{s})$, RQMC reduces them
to
$\mathcal{O}(1/s^{1+2/n})$
for a number $s$ of samples and an input dimension of $n$ [].
For large $n$, we still have a rate of at least
$\mathcal{O}(1/s)\supset\mathcal{O}(1/s^{1+2/n})$
, which constitutes a substantial improvement over the regular reduction in $\mathcal{O}(1/\sqrt{s})$.
However, the default (i.e., Cartesian) QMC and RQMC methods require numbers of samples $s=k^{n}$ for $k\in\mathbb{N}_{+}$, which can become infeasible for large input dimensionalities.
Therefore, we also consider Latin-Hypercube Sampling (LHS) [], which uses a subset of grid cells such that each interval in every dimension is covered exactly once (see Fig. 2 (5+6)).
Finally, we remark that, to our knowledge, QMC and RQMC sampling strategies have not been considered in the field of gradient estimation.

### 2.3 Smoothing of the Algorithm vs. the Objective

In many learning problems, we can write our training objective as $\ell(h(y))$ where $y$ is the output of a neural network model, $h$ is the algorithm, and the scalar function $\ell$ is the training objective (loss function) applied to the output of the algorithms.
In such cases, we can distinguish between smoothing the algorithms ($f=h$) and smoothing the loss ($f=\ell\circ h$).

When smoothing the algorithm, we compute the value and derivative of $\ell\big(\mathbb{E}_{\epsilon}[h(y+\epsilon)]\big)$.
This requires our loss function $\ell$ to be differentiable and capable of receiving relaxed inputs.
(For example, if the output of $h$ is binary, then $\ell$ has to be able to operate on real-valued inputs from $(0,1)$.)
In this case, the derivative of $\mathbb{E}_{\epsilon}[h(y+\epsilon)]$ is a Jacobian matrix (see Corollary 2).

When smoothing the objective / loss function, we compute the value and derivative of $\mathbb{E}_{\epsilon}[\ell(h(y+\epsilon))]$.
Here, the objective / loss $\ell$ does not need to be differentiable and can be limited to operate on discrete outputs of the algorithm $h$.
Here, the derivative of $\mathbb{E}_{\epsilon}[\ell(h(y+\epsilon))]$ is a gradient.

The optimal choice between smoothing the algorithm and smoothing the objective depends on different factors including the problem setting and algorithm, the availability of a real-variate and real-valued $\ell$, and the number of samples that can be afforded.
In practice, we observe that, whenever we can afford large numbers of samples, smoothing of the algorithm performs better.

## 3 Related Work

In the theoretical literature of gradient-free optimization, stochastic smoothing has been extensively studied [, , , ].
Our work extends existing results, generalizing the set of allowed distributions, considering vector-valued functions, anisotropic scale matrices, enabling $k$-sample median differentiation, and a characterization of finite definedness of expectations and their gradients based on the relationship between characteristics of the density and smoothed functions.

From a more applied perspective, stochastic smoothing has been applied for relaxing convex optimization problems [, , ]. In particular, convex optimization formulations of argmax [, ], the shortest-path problem [], and the clustering problem [] have been considered.
We remark that the perspective of smoothing any function or algorithm $f$, as in this work, differs from the perspective of perturbed optimizers.
In particular, optimizers are a special case of the functions we consider.

While we consider smoothing functions with real-valued inputs, there is also a rich literature of differentiating stochastic discrete programs [, , ].
These works typically use the inherent stochasticity from discrete random variables in programs and explicitly model the internals of the programs.
We consider real-variate black-box functions and smooth them with added input noise.

In the literature of reinforcement learning, a special case or analogous idea to stochastic smoothing can be found in the REINFORCE formulation where the (scalar) score function is smoothed via a policy [, , , ].
Compared to the literature, we enable new distributions and respective characterizations of requirements for the score functions.
We hope our results will pave their way into future RL research directions as they are also applicable to RL without major modification.

## 4 Experiments

For the experiments, we consider 4 experimental domains: sorting & ranking, graph algorithms, 3D mesh rendering, and cryo-electron tomography (cryoET) simulations.
The primary objective of the empirical evaluations is to compare different distributions as well as different variance reduction techniques.
We begin our evaluations by measuring the variance of the gradient estimators, and then continue with optimizations and using the differentiable relaxations in deep learning tasks.
We remark that, in each of the 4 experiments, $f$ does not have any non-zero gradients, and thus using first-order or path-wise gradients or gradient estimators is not possible.

Figure: Figure 3: Average $L_{2}$ norms between ground truth (oracle) and estimated gradient for different numbers of elements to sort and rank $n$, and different distributions. Each plot compares different variance reduction strategies as indicated in the legend to the right of the caption. Darker is better (smaller values). Colors are only comparable within each subplot. We use $1\,024$ samples, except for Cartesian and $n=3$ where we use $10^{3}=1\,000$ samples. An extension with $n\in\{7,10\}$ can be found in Figure 11 in the appendix. Absolute values are reported in Table 4.
[Uncaptioned image]: x22.png

### 4.1 Variance of Gradient Estimators

We evaluate the gradient variances for different variance reduction techniques in Figures 4 and 4.
For differentiable sorting and ranking, we smooth the (hard) permutation matrix that sorts an input vector ($f:\mathbb{R}^{n}\to\{0,1\}^{n\times n}$).
For diff. shortest-paths, we smooth the function that maps from a 2D cost-map to a binary encoding of the shortest-path under 8-neighborhood ($f:\mathbb{R}^{n\times n}\to\{0,1\}^{n\times n}$).
Both functions are not only non-differentiable, but also have no non-zero gradients anywhere.
For each distribution, we compare all combinations of the 3 complementary variance reduction techniques.

On the axis of sampling strategy, we can observe that, whenever available, Cartesian RQMC delivers the lowest variance.
The only exception is the triangular distribution, where latin QMC provides the lowest uncentered gradient variance (despite being a biased estimator) because of large contributions to the gradient for samples close to $-1$ and $1$.
Between latin QMC and RQMC, we can observe that their variance is equal except for the high-dimension cases of the Cauchy distribution and a few cases of the Gumbel distribution, where QMC is of lower variance.
However, due to the bias in QMC, RQMC would typically still be preferable over QMC.
We do not consider Cartesian QMC due to its substantially greater bias.
In heuristic conclusion, $\operatorname{RQMC\,(c.)}\succ\operatorname{RQMC\,(l.)}\succeq\operatorname{QMC\,(l.)}\succ\operatorname{MC}$.

On the axis of using antithetic sampling (left vs. right in each subplot), we observe that it consistently performs worse than the regular counterpart, except for vanilla MC without a covariate.
The reason for this is that antithetic sampling does not lead to a good sample-utilization trade-off once we consider quasi Monte-Carlo strategies.
For vanilla Monte-Carlo, antithetic sampling improves the results as long as we do not use the LOO covariate.
Thus, in the following, we consider antithetic only for MC.

On the axis of the covariate, we observe that LOO consistently provides the lowest gradient variances.
This aligns with intuition from Figure 1 where LOO provides the lowest variance at discontinuities (in this subsection, $f$ is discontinuous or constant everywhere).
Comparing no covariate and $f(x)$, the better choice has a strong dependence on the individual setting, which makes sense considering the binary outputs of the algorithms.
$f(x)$ would perform well for functions that attain large values while having fewer discontinuities.

Figure: Figure 5: Sorting benchmark ($n{=}5$). Exact match (EM) accuracy. Brighter is better (greater values). Values between subplots are compara- ble. IQM over 12 seeds and dis- played range of $[75\%,85.5\%]$.
Refer to caption: x24.png

In conclusion, the best setting is Cartesian RQMC with the LOO covariate and without antithetic sampling whenever available (only for $s=k^{n}$ samples for $k\in\mathbb{N}$).
The next best choice is typically RQMC with Latin hypercube sampling.

### 4.2 Differentiable Sorting & Ranking

After investigating the choices of variance reduction techniques wrt. the variance alone, in this section, we explore the utility of stochastic smoothing on the 4-digit MNIST sorting benchmark [].
Here, at each step, a set of $n{=}5$ 4-digit MNIST images (such as ) is presented to a CNN, which predicts the displayed scalar value for each of the $n$ images independently.
For training the model, no absolute information about the displayed value is provided, and only the ordering or ranking of the $n$ images according to their ground truth value is supervised.
The goal is to learn an order-preserving CNN, and the evaluation metric is the fraction of correctly inferred orders from the CNN (exact match accuracy).
Training the CNN requires a differentiable ranking operator (that maps from a vector to a differentiable permutation matrix) for the ranking loss.
Previous work has considered NeuralSort [], SoftSort [], casting sorting as a regularized OT problem [], and differentiable sorting networks (DSNs) [, ].
The state-of-the-art is monotonic DSN [], which utilizes a relaxation based on Cauchy distributions to provide monotonic differentiable sorting, which has strong theoretical and empirical advantages.

**Table 2: Sorting benchmark results ($n=5$), avg. over 12 seeds. ‘best (cv)’ refers to the best sampling strategy, as determined via cross-validation (thus, there is no bias from the selection of the strategy). Table 3 includes additional num. of samples and stds. Baselines are NeuralSort [], SoftSort [], Logistic DSNs [], Cauchy and Error-optimal DSNs [], and OT Sort [], avg. over at least 5 seeds each.**
| Baselines |  | Neu.S. | Soft.S. | L. DSN | C. DSN | E. DSN | OT. S. |
| --- | --- | --- | --- | --- | --- | --- | --- |
| — |  | 71.3 | 70.7 | 77.2 | 84.9 | 85.0 | 81.1 |
| Sampling | #s | Gauss. | Logis. | Gumbel | Cauchy | Laplace | Trian. |
| vanilla | 256 | 82.3 | 82.8 | 79.2 | 68.1 | 82.6 | 81.3 |
| best (cv) | 256 | 83.1 | 82.7 | 81.6 | 55.6 | 83.7 | 82.7 |
| vanilla | 1k | 81.3 | 83.7 | 82.0 | 68.5 | 80.6 | 82.8 |
| best (cv) | 1k | 83.9 | 84.0 | 84.2 | 73.0 | 84.3 | 82.4 |
| vanilla | 32k | 84.2 | 84.1 | 84.5 | 84.9 | 84.4 | 83.4 |
| best (cv) | 32k | 84.4 | 84.4 | 84.8 | 85.1 | 84.4 | 84.0 |

In Figure 5, we evaluate the performance of generalized stochastic smoothing with different distributions and different numbers of samples for each variance reduction technique.
We observe that, while the Cauchy distribution performs poorly for small numbers of samples, for large numbers of samples, the Cauchy distribution performs best.
This makes sense as the Cauchy distribution has infinite variance and, for DSNs, provides monotonicity.
We remark that large numbers of samples can easily be afforded in many applications (when comparing the high cost of neural networks to the vanishing cost of sorting/ranking within a loss function).
(Nevertheless, for $32\,768$ samples, the sorting operation starts to become the bottleneck.)
The
Laplace distribution is the best choice for smaller numbers of samples, which aligns with the characterization of it having the lowest variance because all samples contribute equally to the gradient.
Wrt. variance reduction, we continue to observe that vanilla MC performs worst.
RQMC performs best, except for Triangular, where QMC is best.
For the Gumbel distribution, we observe reduced performance for latin sampling.
Generally, we observe that $f(x)$ is the worst choice of covariate, but the effect lies within standard deviations.
In Table 2, we provide a numerical comparison to other differentiable sorting approaches.
We can observe that all choices of distributions improve over all baselines except for the monotonic DSNs, even at smaller numbers of samples (i.e., without measurable impact on training speed).
Finally, the Cauchy distribution leads to a minor improvement over the SOTA, without requiring a manually designed differentiable sorting algorithm; however, only at the computational cost of $32\,768$ samples.

Figure: Figure 6: Warcraft shortest-path experiment with 1000 samples. Brighter is better (larger values). Values between subplots are comparable. Exact match accuracy avg. over 5 seeds and displayed range $[70\%,96\%]$. Additional settings in Figures 13 and 14.
Refer to caption: x25.png

### 4.3 Differentiable Shortest-Paths

The Warcraft shortest-path benchmark [] is the established benchmark for differentiable shortest-path algorithms (e.g., [, , ]).
Here, a Warcraft pixel map is provided, a CNN predicts a $12\times 12$ cost matrix, a differentiable algorithm computes the shortest-path, and the supervision is only the ground truth shortest-path.
Berthet et al. [] considered stochastic smoothing with Fenchel-Young (FY) losses, which improves sample efficiency for small numbers of samples.
However, the FY loss does not improve for larger numbers of samples (e.g., Tab. 7.5 in []).

Figure: Figure 7: Warcraft shortest-path experiment using Gaussian smoothing of the algorithm (RQMC with latin hypercube-sampling and LOO covariate). Comparing the effects between the inverse temperature $\beta$ and the number of samples. We observe that with growing numbers of samples, the optimal inverse temperature increases, i.e., the optimal standard deviation for the Gaussian noise decreases. Averaged over 5 seeds.
Refer to caption: x26.png

As computing the shortest-path is computationally efficient and parallelizable (our implementation $\approx 5\,000\times$ faster than the Dijkstra implementation used in previous work [, ]), we can afford substantially larger numbers of samples, improving the quality of gradient estimation.
In Figure 6, we compare the performance of different smoothing strategies.
The logistic distribution performs best, and smoothing of the algorithm (top) performs better than smoothing of the loss (bottom).
Variance reduction via sampling strategies (antithetic, QMC, or RQMC) improves performance, and the best covariate is LOO.
For reference, the FY loss [] leads to an accuracy of $80.6\%$, regardless of the number of samples.
GSS consistently achieves $90\%+$ using 100 samples (see Fig. 13 right).
Using 10 000 samples, and variance reduction, we achieve $96.6\%$ in the best setting (Fig. 14) compared to the SOTA of $95.8\%$ []. In Fig. 7, we illustrate that smaller standard deviations (larger $\beta$) are better for more samples.

### 4.4 Differentiable Rendering

Figure: Figure 8: Utah teapot camera pose optimization. The metric is fraction of camera poses recovered; the initialization an- gle errors are uniformly distributed in $[15^{\circ},75^{\circ}]$. Brighter is better. Avg. over $768$ seeds. The dis- played range is $[0\%,90\%]$.
Refer to caption: x27.png

For differentiable rendering [, , , , , , ], we smooth a non-differentiable hard renderer via sampling.
This differs from DRPO [], which uses stochastic smoothing to relax the Heaviside and Argmax functions within an already differentiable renderer.
Instead, we consider the renderer as a black-box function.
This has the advantage of noise parameterized in the coordinate space rather than the image space.

We benchmark stochastic smoothing for rendering by optimizing the camera-pose (4-DoF) for a Utah teapot, an experiment inspired by [, ].
We illustrate the results in Figure 8.
Here, the logistic distribution performs best, and QMC/RQMC as well as LOO lead to the largest improvements.
While Fig. 8 shows smoothing the rendering algorithm, Fig. 12 performs smoothing of the training objective / loss.
Smoothing the algorithm is better because the loss (MSE), while well-defined on discrete renderings, is less meaningful on discrete renderings.

### 4.5 Differentiable Cryo-Electron Tomography

Transmission Electron Microscopy (TEM) transmits electron beams through thin specimens to form images [].
Due to the small electron beam wavelength, TEM leads to higher resolutions of

Figure: Figure 9: (a) Simulated Transmission Electron micrograph, (b) TMV structure with RNA (orange) and protein stacks (blue).
Refer to caption: x28.png

up to single columns of atoms. Obtaining high resolution images from TEM involves adjustments of various experimental parameters.
We apply smoothing to a realistic black-box TEM simulator [], optimizing sets of parameters to approximate reference Tobacco Mosaic Virus (TMV) [] micrographs. In
Figure 10, we perform two experiments: a 2-parameter study optimizing the microscope acceleration voltage and $x$-position of the specimen, and a 4-parameter study with additional parameters of the particle’s $y$-position and the primary lens focal length.
The micrograph image sizes are $400\times 400$ pixels, and accordingly we use smoothing of the loss.

Figure: Figure 10: RMSE to Ground Truth parameters for the 2-parameter (left) and 4-parameter experiment (right). We optimize the $L_{2}$ loss between generated and GT images using loss smoothing. No marker lines correspond to Gaussian, $\times$ to Laplace and $\triangle$ to Triangular distributions. Laplace and Triangular perform best; LOO leads to the largest improvements. Add. results are in Figure 15.
Refer to caption: x29.png

##### Summary of Experimental Results

Generally, we observe that QMC and RQMC perform best, whereas antithetic sampling performs rather poorly.
In low-dimensional problems, it is advisable to use RQMC (cartesian), and in higher dimensional problems (R)QMC (latin), still works well.
As for the covariate, LOO typically performs best; however, the choice of sampling strategy (QMC/RQMC) is more important than choosing the covariate.
In sorting and ranking, the Cauchy distribution performs best for large numbers of samples and for smaller numbers of samples, the Laplace distribution performs best.
In the shortest-path case, the logistic distribution performs best, and Gaussian closely follows. Here, we also observe that with larger numbers of samples, the optimal standard deviation decreases.
For differentiable rendering, the logistic distribution performs best.

##### Limitations

A limitation of our work is that zeroth-order gradient estimators are generally only competitive if the first-order gradients do not exist (see [] for discussions on exceptions).
In this vein, in order to be competitive with custom designed continuous relaxations like a differentiable renderer, we may need a very large number of samples, which could become prohibitive for expensive functions $f$.
The optimal choice of distribution depends on the function to be smoothed, which means there is no singular distribution that is optimal for all $f$; however, if one wants to limit the distribution to a single choice, we recommend the logistic or Laplace distribution, as, with their simple exponential convergence, they give a good middle ground between heavy-tailed and light-tailed distributions.
Finally, the variance reduction techniques like QMC/RQMC are not immediately applicable in single sample settings, and the variance reduction techniques in this paper build on evaluating $f$ many times.

## 5 Conclusion

In this work, we derived stochastic smoothing with reduced assumptions and outline a general framework for relaxation and gradient estimation of non-differentiable black-box functions.
This enables an increased set of distributions for stochastic smoothing, e.g., enabling smoothing with the triangular distribution while maintaining full differentiablility of $f_{\epsilon}$.
We investigated variance reduction for stochastic smoothing–based gradient estimation from 3 orthogonal perspectives, finding that RQMC and LOO are generally the best methods, whereas the popular antithetic sampling method performs rather poorly.
Moreover, enabled by supporting vector-valued functions, we disentangled the algorithm and objective, thus smoothing $f$ while analytically backpropagating through the loss $\ell$, improving gradient estimation.
We applied stochastic smoothing to differentiable sorting and ranking, diff. shortest-paths on graphs, diff. rendering for pose estimation and diff. cryo-ET simulations.
We hope that our work inspires the community to develop their own stochastic relaxations for differentiating non-differentiable algorithms, operators, and simulators.

#### Acknowledgments

We would like to acknowledge helpful discussions with Michael Kagan, Daniel Ratner, and Terry Suh.
This work was in part supported by the Land Salzburg within the WISS 2025 project IDA-Lab (20102-F1901166-KZP and 20204-WISS/225/197-2019), the U.S. DOE Contract No. DE-AC02-76SF00515, Zoox Inc, ARO (W911NF-21-1-0125), ONR (N00014-23-1-2159), the CZ Biohub, and SPRIN-D.

## Appendix A Proofs

### A.1 Proof of Lemma 3

###### Proof of Lemma 3 .

Let $\Omega\subset\mathbb{R}^{n}$ be the set of values where $\nabla_{\!\epsilon}\mu(\epsilon)$ is undefined.
$\mu$ is differentiable a.e. and $\Omega$ has Lebesgue measure 0.

Recapitulating (5) from the proof of Lemma 1, we have

$$ $\nabla_{\!x}f_{\epsilon}(x)=-\int f(x+\epsilon)\,\nabla_{\!\epsilon}\mu(\epsilon)\,d\epsilon\,.$ (17) $$

Replacing $\nabla_{\!\epsilon}\mu(\epsilon)$ by any of the weak derivatives $\nu$ of $\mu$, which exists and is integrable due to absolute continuity, we have

$$ $\displaystyle\nabla_{\!x}f_{\epsilon}(x)$ $\displaystyle=-\int f(x+\epsilon)\,\nu(\epsilon)\,d\epsilon\,$ (18) $\displaystyle=-\int_{\mathbb{R}^{n}\setminus\Omega}f(x+\epsilon)\,\nu(\epsilon)\,d\epsilon\,-\int_{\Omega}f(x+\epsilon)\,\nu(\epsilon)\,d\epsilon\,.$ (19) $$

Because $\mu$ is absolutely continuous
and as the Lebesgue measure of $\Omega$ is 0, per Hölder’s inequality

$$ $\int_{\Omega}|f(x+\epsilon)\,\nu(\epsilon)|\,d\epsilon\leq\int_{\Omega}|f(x+\epsilon)|\,d\epsilon\,\cdot\int_{\Omega}|\nu(\epsilon)|\,d\epsilon=\int_{\Omega}|f(x+\epsilon)|\,d\epsilon\,\cdot 0=0$ (20) $$

where $\int_{\Omega}|\nu(\epsilon)|\,d\epsilon=0$ follows from absolute continuity of $\mu$.
Thus,

$$ $\int_{\Omega}f(x+\epsilon)\,\nu(\epsilon)\,d\epsilon=0\,.$ (21) $$

As $\nu=\nabla_{\!\epsilon}\mu(\epsilon)$ for all $\epsilon\in\mathbb{R}^{n}\setminus\Omega$

$$ $\nabla_{\!x}f_{\epsilon}(x)=-\int_{\mathbb{R}^{n}\setminus\Omega}f(x+\epsilon)\,\nu(\epsilon)\,d\epsilon\,-\int_{\Omega}f(x+\epsilon)\,\nu(\epsilon)\,d\epsilon\,=-\int_{\mathbb{R}^{n}\setminus\Omega}f(x+\epsilon)\,\nabla_{\!\epsilon}\mu(\epsilon)\,d\epsilon\,,$ (22) $$

showing that for all possible choices of $\nu$, the gradient estimator coincides.
Thus, we complete our proof via

$$ $\nabla_{\!x}f_{\epsilon}(x)=-\int_{\mathbb{R}^{n}\setminus\Omega}f(x+\epsilon)\,\mu(\epsilon)\,\nabla_{\!\epsilon}\log\mu(\epsilon)\,d\epsilon=\mathbb{E}_{\epsilon\sim\mu}\,\Big[f(x+\epsilon)\,\cdot\mathbf{1}_{\epsilon\notin\Omega}\cdot\,\nabla_{\!\epsilon}{-}\log\mu(\epsilon)\Big]\,.$ (23) $$

After completing the proof, we remark that, if the density was not continuous, e.g., uniform $\mathcal{U}([0,1])$, then $\int_{\{0\}}\nabla_{\!\epsilon}\mu(\epsilon)\,d\epsilon=\Big[\mu(\epsilon)\Big]_{\epsilon\nearrow 0}^{\epsilon\searrow 0}=1$.
This means that the weak derivative is not defined (or loosely speaking “the derivative is infinity”), thereby violating the assumptions of Hölder’s inequality (Eq. 20).
This concludes that continuity is required for the proof to hold. ∎

###### Proof of Lemma 3 .

### A.2 Proof of Lemma 6

###### Proof of Lemma 6 .

$$ $\displaystyle\nabla_{\!\gamma}f_{\gamma\epsilon}(x)$ $\displaystyle=$ $\displaystyle\nabla_{\!\gamma}\mathbb{E}_{\epsilon\sim\mu}\left[f(x+\gamma\cdot\epsilon)\right]$ (24) $\displaystyle=$ $\displaystyle\nabla_{\!\gamma}\int f(x+\gamma\cdot\epsilon)\mu(\epsilon)d\epsilon$ (26) $\displaystyle{\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\textstyle\big(u=x+\epsilon\cdot\gamma\Rightarrow\epsilon=\frac{u-x}{\gamma}~;~\frac{du}{d\epsilon}=\gamma\Rightarrow d\epsilon=\frac{1}{\gamma}du\big)}$ $\displaystyle=$ $\displaystyle\nabla_{\!\gamma}\int f(u)\mu(\epsilon){\textstyle\frac{1}{\gamma}}du$ (27) $\displaystyle=$ $\displaystyle\int f(u)\nabla_{\!\gamma}(\mu(\epsilon){\textstyle\frac{1}{\gamma}})du$ (28) $\displaystyle=$ $\displaystyle\int f(u)({\textstyle\frac{1}{\gamma}}\nabla_{\!\gamma}\mu(\epsilon)+\mu(\epsilon)\nabla_{\!\gamma}{\textstyle\frac{1}{\gamma}})du$ (29) $\displaystyle=$ $\displaystyle\int f(u)({\textstyle\frac{1}{\gamma}}(\nabla_{\!\epsilon}\mu(\epsilon))^{\!\top}\textstyle{\frac{\partial\epsilon}{\partial\gamma}}-\mu(\epsilon){\textstyle\frac{1}{\gamma^{2}}})du$ (30) $\displaystyle=$ $\displaystyle\int f(u)({\textstyle\frac{1}{\gamma}}(\nabla_{\!\epsilon}\mu(\epsilon))^{\!\top}\textstyle{\frac{\partial}{\partial\gamma}\frac{u-x}{\gamma}}-\mu(\epsilon){\textstyle\frac{1}{\gamma^{2}}})du$ (31) $\displaystyle=$ $\displaystyle\int f(u)({\textstyle\frac{1}{\gamma}}(\nabla_{\!\epsilon}\mu(\epsilon))^{\!\top}(-\textstyle{\frac{\epsilon}{\gamma}})-{\textstyle\frac{1}{\gamma^{2}}}\mu(\epsilon))du$ (32) $\displaystyle=$ $\displaystyle\int f(u)(-(\nabla_{\!\epsilon}\mu(\epsilon))^{\top}\epsilon-\mu(\epsilon)){\textstyle\frac{1}{\gamma^{2}}}du$ (34) $\displaystyle{\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\textstyle\left(\nabla_{\!\epsilon}\log\mu(\epsilon)=\frac{1}{\mu(\epsilon)}\nabla_{\!\epsilon}\mu(\epsilon)\Rightarrow\nabla_{\!\epsilon}\mu(\epsilon)=\mu(\epsilon)\nabla_{\!\epsilon}\log\mu(\epsilon)\right)}$ $\displaystyle=$ $\displaystyle\int f(u)(-(\mu(\epsilon)\nabla_{\!\epsilon}\log\mu(\epsilon))^{\top}\epsilon-\mu(\epsilon)){\textstyle\frac{1}{\gamma^{2}}}\cdot\underbrace{\gamma\,d\epsilon}_{=du}$ (35) $\displaystyle=$ $\displaystyle\int f(u)\cdot(-(\nabla_{\!\epsilon}\log\mu(\epsilon))^{\top}\epsilon-1)\cdot{\textstyle\frac{1}{\gamma}}\cdot\mu(\epsilon)\,d\epsilon$ (36) $\displaystyle=$ $\displaystyle\int f(u)\cdot(-1+(\nabla_{\!\epsilon}{-}\log\mu(\epsilon))^{\top}\epsilon)\cdot{\textstyle\frac{1}{\gamma}}\cdot\mu(\epsilon)\,d\epsilon$ (37) $\displaystyle=$ $\displaystyle\mathbb{E}_{\epsilon\sim\mu}\big[f(x+\gamma\cdot\epsilon)\cdot\big({-}1+(\nabla_{\!\epsilon}{-}\log\mu(\epsilon))^{\top}\cdot\epsilon\big)\,/\,\gamma\big]\,.$ (38) $$

∎

###### Proof of Lemma 6 .

### A.3 Proof of Theorem 7

###### Proof of Theorem 7 .

Part 1: $\partial f_{\mathbf{L}\epsilon}(x)\,/\,\partial x$

We perform a change of variables, $u=x+\mathbf{L}\epsilon\implies\epsilon=\mathbf{L}^{-1}(u-x)$ and

$$ $\displaystyle d\epsilon=\frac{du}{du}d\epsilon=\frac{d\epsilon}{du}du=\frac{d\mathbf{L}^{-1}(u-x)}{du}du=\frac{d\mathbf{L}^{-1}u}{du}du=\det(\mathbf{L}^{-1})\,du$ (39) $$

Thus,

$$ $\displaystyle f_{\mathbf{L}\epsilon}(x)=\int f(x+\mathbf{L}\epsilon)\mu(\epsilon)\,d\epsilon=\int f(u)\cdot\mu(\mathbf{L}^{-1}(u-x))\cdot\det(\mathbf{L}^{-1})\,du\,.$ (40) $$

Now,

$$ $\displaystyle\nabla_{\!x}f_{\mathbf{L}\epsilon}(x)_{i}$ $\displaystyle=\nabla_{\!x}\int f(u)_{i}\cdot\mu(\mathbf{L}^{-1}(u-x))\cdot\det(\mathbf{L}^{-1})\,du$ (41) $\displaystyle=\int f(u)_{i}\,\cdot\nabla_{\!x}\big(\mu(\mathbf{L}^{-1}(u-x))\big)\cdot\det(\mathbf{L}^{-1})\,du$ (42) $\displaystyle=\int f(u)_{i}\,\cdot\mathbf{L}^{-1}\cdot\big(\nabla_{\!\epsilon}{-}\mu(\epsilon)\big)\cdot\det(\mathbf{L}^{-1})\,du$ (43) $\displaystyle=\int f(x+\mathbf{L}\epsilon)_{i}\,\cdot\mathbf{L}^{-1}\cdot\nabla_{\!\epsilon}{-}\mu(\epsilon)\,d\epsilon$ (44) $\displaystyle=\int f(x+\mathbf{L}\epsilon)_{i}\,\cdot\mathbf{L}^{-1}\cdot\mu(\epsilon)\cdot\nabla_{\!\epsilon}{-}\log\mu(\epsilon)\,d\epsilon$ (45) $\displaystyle=\mathbb{E}_{\epsilon\sim\mu}\Big[f(x+\mathbf{L}\epsilon)_{i}\,\cdot\mathbf{L}^{-1}\cdot\nabla_{\!\epsilon}{-}\log\mu(\epsilon)\Big]$ (46) $$

Part 2: $\partial f_{\mathbf{L}\epsilon}(x)\,/\,\partial\mathbf{L}$

We use the same change of variables as above.

$$ $\displaystyle\nabla_{\!\mathbf{L}}\,\mathbb{E}_{\epsilon\sim\mu}\big[f(x+\mathbf{L}\cdot\epsilon)_{i}\big]$ (47) $\displaystyle=\nabla_{\!\mathbf{L}}\int f(x+\mathbf{L}\epsilon)_{i}\cdot\mu(\epsilon)\,d\epsilon$ (48) $\displaystyle=\nabla_{\!\mathbf{L}}\int f(u)_{i}\cdot\mu(\mathbf{L}^{-1}(u-x))\cdot\det(\mathbf{L}^{-1})\,du$ (49) $\displaystyle=\int f(u)_{i}\cdot\nabla_{\!\mathbf{L}}\Big(\mu(\mathbf{L}^{-1}(u-x))\cdot\det(\mathbf{L^{-1}})\Big)\,du$ (50) $\displaystyle=\int f(x+\mathbf{L}\epsilon)_{i}\cdot\nabla_{\!\mathbf{L}}\Big(\mu(\mathbf{L}^{-1}(u-x))\cdot\det(\mathbf{L^{-1}})\Big)\,/\det(\mathbf{L^{-1}})\,d\epsilon$ (51) $\displaystyle=\mathbb{E}_{\epsilon\sim\mu}\Big[f(x+\mathbf{L}\epsilon)_{i}\cdot\nabla_{\!\mathbf{L}}\Big(\mu(\mathbf{L}^{-1}(u-x))\cdot\det(\mathbf{L^{-1}})\Big)\,\cdot\det(\mathbf{L})\,/\mu(\epsilon)\Big]$ (52) $$

Now, while $\nabla_{\!\mathbf{L}}\Big(\mu(\mathbf{L}^{-1}(u-x))\cdot\det(\mathbf{L^{-1}})\Big)$ may be computed via automatic differentiation, we can also solve it in closed-form.
Firstly, we can observe that

$$ $\displaystyle\nabla_{\!\mathbf{L}}\mu(\mathbf{L}^{-1}(u-x))$ $\displaystyle=\nabla_{\!\epsilon^{\top}}\mu(\epsilon)\cdot\nabla_{\!\mathbf{L}}(\mathbf{L}^{-1}\cdot(u-x))$ (53) $\displaystyle=\nabla_{\!\mathbf{L}}(\nabla_{\!\epsilon^{\top}}\mu(\epsilon)\cdot\mathbf{L}^{-1}\cdot(u-x))$ (54) $\displaystyle=-\,\mathbf{L}^{-\top}\cdot\nabla_{\!\epsilon}\mu(\epsilon)\cdot(u-x)^{\top}\cdot\mathbf{L}^{-\top}$ (55) $$

and

$$ $\displaystyle\nabla_{\!\mathbf{L}}\det(\mathbf{L}^{-1})$ $\displaystyle=-\det(\mathbf{L})^{-1}\cdot\mathbf{L}^{-\top}\,.$ (56) $$

We can combine this to resolve it in closed form to:

$$ $\displaystyle\nabla_{\!\mathbf{L}}\Big(\mu(\mathbf{L}^{-1}(u-x))\cdot\det(\mathbf{L^{-1}})\Big)$ $\displaystyle=-\,\mathbf{L}^{-\top}\cdot\nabla_{\!\epsilon}\mu(\epsilon)\cdot(u-x)^{\top}\cdot\mathbf{L}^{-\top}\cdot\det(\mathbf{L^{-1}})$ $\displaystyle~~~~\,-\mu(\mathbf{L}^{-1}(u-x))\cdot\det(\mathbf{L})^{-1}\,\cdot\,\mathbf{L}^{-\top}$ (57) $\displaystyle=-\,\mathbf{L}^{-\top}\cdot\nabla_{\!\epsilon}\mu(\epsilon)\cdot\big(\mathbf{L}^{-1}(u-x)\big)^{\top}\cdot\det(\mathbf{L^{-1}})$ $\displaystyle~~~~\,-\mu(\epsilon)\cdot\det(\mathbf{L})^{-1}\,\cdot\,\mathbf{L}^{-\top}$ (58) $\displaystyle=-\,\mathbf{L}^{-\top}\cdot\nabla_{\!\epsilon}\mu(\epsilon)\cdot\epsilon^{\top}\cdot\det(\mathbf{L^{-1}})$ $\displaystyle~~~~\,-\mu(\epsilon)\cdot\det(\mathbf{L})^{-1}\cdot\mathbf{L}^{-\top}$ (59) $\displaystyle=-\det(\mathbf{L^{-1}})\cdot\big(\mathbf{L}^{-\top}\cdot\nabla_{\!\epsilon}\mu(\epsilon)\cdot\epsilon^{\top}+\mu(\epsilon)\cdot\mathbf{L}^{-\top}\big)\,.$ (60) $$

Combing this with equation (52), we have

$$ $\displaystyle\nabla_{\!\mathbf{L}}\,\mathbb{E}_{\epsilon\sim\mu}\big[f(x+\mathbf{L}\cdot\epsilon)_{i}\big]$ $\displaystyle=\mathbb{E}_{\epsilon\sim\mu}\Bigg[f(x+\mathbf{L}\epsilon)_{i}\cdot\nabla_{\!\mathbf{L}}\Big(\mu(\mathbf{L}^{-1}(u-x))\cdot\det(\mathbf{L^{-1}})\Big)\,\cdot\det(\mathbf{L})\,/\mu(\epsilon)\Bigg]$ $\displaystyle=\mathbb{E}_{\epsilon\sim\mu}\Bigg[f(x+\mathbf{L}\epsilon)_{i}\,\cdot\,-\det(\mathbf{L^{-1}})\cdot\big(\mathbf{L}^{-\top}\cdot\nabla_{\!\epsilon}\mu(\epsilon)\cdot\epsilon^{\top}+\mu(\epsilon)\cdot\mathbf{L}^{-\top}\big)\,\cdot\,\det(\mathbf{L})\,/\,\mu(\epsilon)\Bigg]$ (61) $\displaystyle=\mathbb{E}_{\epsilon\sim\mu}\Bigg[f(x+\mathbf{L}\epsilon)_{i}\,\cdot\,-\big(\mathbf{L}^{-\top}\cdot\nabla_{\!\epsilon}\mu(\epsilon)\cdot\epsilon^{\top}+\mu(\epsilon)\cdot\mathbf{L}^{-\top}\big)\,/\,\mu(\epsilon)\Bigg]$ (62) $\displaystyle=\mathbb{E}_{\epsilon\sim\mu}\Bigg[f(x+\mathbf{L}\epsilon)_{i}\,\cdot\,-\big(\mathbf{L}^{-\top}\cdot\nabla_{\!\epsilon}\mu(\epsilon)\cdot\epsilon^{\top}\,/\,\mu(\epsilon)+\mathbf{L}^{-\top}\big)\Bigg]$ (63) $\displaystyle=\mathbb{E}_{\epsilon\sim\mu}\Big[f(x+\mathbf{L}\epsilon)_{i}\,\cdot\,\mathbf{L}^{-\top}\,\cdot\,\big(-1+\nabla_{\!\epsilon}{-}\log\mu(\epsilon)\cdot\epsilon^{\top}\big)\Big]\,.$ (64) $$

∎

###### Proof of Theorem 7 .

### A.4 Proof of Theorem 8

###### Proof of Theorem 8 .

Part 1: $\partial G_{\mathbf{L}\epsilon}(x)\,/\,\partial x$
We start by stating some helpful preliminaries:

$$ $\displaystyle G_{\mathbf{L}\epsilon}(x)$ $\displaystyle=\operatorname{Cov}_{\epsilon\sim\mu}(f(x+\mathbf{L}\epsilon))=\mathbb{E}_{\epsilon\sim\mu}\Big[\big(f(x+\mathbf{L}\epsilon)-f_{\mathbf{L}\epsilon}(x)\big)\,\big(f(x+\mathbf{L}\epsilon)-f_{\mathbf{L}\epsilon}(x)\big)^{\top}\Big]$ (65) $\displaystyle=\mathbb{E}_{\epsilon\sim\mu}\Big[f(x+\mathbf{L}\epsilon)\,f(x+\mathbf{L}\epsilon)^{\top}\Big]-f_{\mathbf{L}\epsilon}(x)\,f_{\mathbf{L}\epsilon}(x)^{\top}\,,$ (66) $\displaystyle G_{\mathbf{L}\epsilon}(x)_{i,j}$ $\displaystyle=\operatorname{Cov}_{\epsilon\sim\mu}(f(x+\mathbf{L}\epsilon)_{i},f(x+\mathbf{L}\epsilon)_{j})$ (67) $\displaystyle=\mathbb{E}_{\epsilon\sim\mu}\Big[\big(f(x+\mathbf{L}\epsilon)-f_{\mathbf{L}\epsilon}(x)\big)_{i}\,\cdot\,\big(f(x+\mathbf{L}\epsilon)-f_{\mathbf{L}\epsilon}(x)\big)_{j}\Big]$ (68) $\displaystyle=\mathbb{E}_{\epsilon\sim\mu}\Big[f(x+\mathbf{L}\epsilon)_{i}\cdot f(x+\mathbf{L}\epsilon)_{j}\Big]-f_{\mathbf{L}\epsilon}(x)_{i}\cdot f_{\mathbf{L}\epsilon}(x)_{j}\,.$ (69) $$

We proceed by computing the derivative of the left part of Equation 69.

$$ $\displaystyle\nabla_{x}\,\mathbb{E}_{\epsilon\sim\mu}\Big[f(x+\mathbf{L}\epsilon)_{i}\cdot f(x+\mathbf{L}\epsilon)_{j}\Big]$ (70) $\displaystyle=\mathbb{E}_{\epsilon\sim\mu}\Big[f(x+\mathbf{L}\epsilon)_{i}\cdot f(x+\mathbf{L}\epsilon)_{j}\,\cdot\mathbf{L}^{-1}\cdot\nabla_{\epsilon}\big(-\log\mu(\epsilon)\big)\Big]\,,$ (71) $$

which is analogous to Equations 41 until 46.
Now,

$$ $\displaystyle\nabla_{x}\,G_{\mathbf{L}\epsilon}(x)_{i,j}$ $\displaystyle=\nabla_{x}\,\mathbb{E}_{\epsilon\sim\mu}\Big[f(x+\mathbf{L}\epsilon)_{i}\cdot f(x+\mathbf{L}\epsilon)_{j}\Big]-\nabla_{x}\,\big(f_{\mathbf{L}\epsilon}(x)_{i}\cdot f_{\mathbf{L}\epsilon}(x)_{j}\big)$ (72) $\displaystyle=\mathbb{E}_{\epsilon\sim\mu}\Big[f(x+\mathbf{L}\epsilon)_{i}\cdot f(x+\mathbf{L}\epsilon)_{j}\,\cdot\mathbf{L}^{-1}\cdot\nabla_{\epsilon}\big(-\log\mu(\epsilon)\big)\Big]$ (73) $\displaystyle\quad-f_{\mathbf{L}\epsilon}(x)_{i}\cdot\nabla_{x}\,f_{\mathbf{L}\epsilon}(x)_{j}~-f_{\mathbf{L}\epsilon}(x)_{j}\cdot\nabla_{x}\,f_{\mathbf{L}\epsilon}(x)_{i}$ (74) $$

where $\nabla_{x}\,f_{\mathbf{L}\epsilon}(x)$ is defined as in part 1 of the proof of Theorem 7.

Part 2: $\partial G_{\mathbf{L}\epsilon}(x)\,/\,\partial\mathbf{L}$

$$ $\displaystyle\nabla_{\mathbf{L}}\,G_{\mathbf{L}\epsilon}(x)_{i,j}$ (75) $\displaystyle=\nabla_{\mathbf{L}}\,\mathbb{E}_{\epsilon\sim\mu}\Big[f(x+\mathbf{L}\epsilon)_{i}\cdot f(x+\mathbf{L}\epsilon)_{j}\Big]-\nabla_{\mathbf{L}}\,\big(f_{\mathbf{L}\epsilon}(x)_{i}\cdot f_{\mathbf{L}\epsilon}(x)_{j}\big)$ (76) $\displaystyle=\mathbb{E}_{\epsilon\sim\mu}\Big[f(x+\mathbf{L}\epsilon)_{i}\cdot f(x+\mathbf{L}\epsilon)_{j}\cdot\nabla_{\mathbf{L}}\Big(\mu(\mathbf{L}^{-1}(u-x))\cdot\det(\mathbf{L^{-1}})\Big)\,\cdot\det(\mathbf{L})\,/\mu(\epsilon)\Big]$ (77) $\displaystyle\quad-f_{\mathbf{L}\epsilon}(x)_{i}\cdot\nabla_{\mathbf{L}}\,f_{\mathbf{L}\epsilon}(x)_{j}~-f_{\mathbf{L}\epsilon}(x)_{j}\cdot\nabla_{\mathbf{L}}\,f_{\mathbf{L}\epsilon}(x)_{i}$ (78) $$

where $\nabla_{\mathbf{L}}\,f_{\mathbf{L}\epsilon}(x)$ is defined as in Theorem 7.

Using the closed-form solution for $\nabla_{\mathbf{L}}\Big(\mu(\mathbf{L}^{-1}(u-x))\cdot\det(\mathbf{L^{-1}})\Big)$ from the proof of Theorem 7, we can simplify it to

$$ $\displaystyle\nabla_{\mathbf{L}}\,$ $\displaystyle G_{\mathbf{L}\epsilon}(x)_{i,j}$ (79) $\displaystyle=~$ $\displaystyle\mathbb{E}_{\epsilon\sim\mu}\Big[f(x+\mathbf{L}\epsilon)_{i}\cdot f(x+\mathbf{L}\epsilon)_{j}\cdot\nabla_{\mathbf{L}}\Big(\mu(\mathbf{L}^{-1}(u-x))\cdot\det(\mathbf{L^{-1}})\Big)\,\cdot\det(\mathbf{L})\,/\mu(\epsilon)\Big]$ $\displaystyle-f_{\mathbf{L}\epsilon}(x)_{i}\cdot\nabla_{\mathbf{L}}\,f_{\mathbf{L}\epsilon}(x)_{j}~-f_{\mathbf{L}\epsilon}(x)_{j}\cdot\nabla_{\mathbf{L}}\,f_{\mathbf{L}\epsilon}(x)_{i}$ (80) $\displaystyle=~$ $\displaystyle-\mathbb{E}_{\epsilon\sim\mu}\Big[f(x+\mathbf{L}\epsilon)_{i}\cdot f(x+\mathbf{L}\epsilon)_{j}\cdot\big(\mathbf{L}^{-\top}\cdot\nabla_{\epsilon}\mu(\epsilon)\cdot\epsilon^{\top}+\mu(\epsilon)\cdot\mathbf{L}^{-\top}\big)\,/\mu(\epsilon)\Big]$ $\displaystyle-f_{\mathbf{L}\epsilon}(x)_{i}\cdot\nabla_{\mathbf{L}}\,f_{\mathbf{L}\epsilon}(x)_{j}~-f_{\mathbf{L}\epsilon}(x)_{j}\cdot\nabla_{\mathbf{L}}\,f_{\mathbf{L}\epsilon}(x)_{i}$ (81) $\displaystyle=~$ $\displaystyle-\mathbb{E}_{\epsilon\sim\mu}\Big[f(x+\mathbf{L}\epsilon)_{i}\cdot f(x+\mathbf{L}\epsilon)_{j}\cdot\big(\mathbf{L}^{-\top}\cdot\nabla_{\epsilon}\mu(\epsilon)\cdot\epsilon^{\top}\,/\mu(\epsilon)+\mathbf{L}^{-\top}\big)\Big]$ $\displaystyle-f_{\mathbf{L}\epsilon}(x)_{i}\cdot\nabla_{\mathbf{L}}\,f_{\mathbf{L}\epsilon}(x)_{j}~-f_{\mathbf{L}\epsilon}(x)_{j}\cdot\nabla_{\mathbf{L}}\,f_{\mathbf{L}\epsilon}(x)_{i}$ (82) $$

∎

###### Proof of Theorem 8 .

## Appendix B Discussion of Properties of f f for Finitely Defined f ϵ f_{\epsilon} and ∇ f ϵ \nabla f_{\epsilon}

When we have a function $f$ that is not defined with a compact range with $f:\mathbb{R}^{n}\to\mathbb{R}$, and have a density $\mu$ with unbounded support (e.g., Gaussian or Cauchy), we may experience $f_{\epsilon}$ or even $\nabla f_{\epsilon}$ to not be finitely defined.
For example, virtually any distribution with full support on $\mathbb{R}$ leads to the smoothing $f_{\epsilon}$ of the degenerate function $f:x\mapsto\exp(\exp(\exp(\exp(x^{2}))))$ to not be finitely defined.

We say a function, as described via an expectation, is finitely defined iff it is defined (i.e., the expectation has a value) and its value is finite (i.e., not infinity).
For example, the first moment of the Cauchy distribution is undefined, and the second moment is infinite; thus, both moments are not finitely defined.

We remark that the considerations in this appendix also apply to prior works that enable the real plane as the output space of $f$.
We further remark that writing an expression for smoothing and the gradient of a arbitrary function with non-compact range is not necessarily false; however, e.g., any claim that smoothness is guaranteed if the gradient jumps from $-\infty$ to $\infty$ (e.g., the power tower in the first paragraph) is not formally correct.
We remark that characterizing valid $f$s via a Lipschitz or other continuity requirement is not applicable because this would defeat the goal of differentiating non-differentiable and discontinuous $f$.

In the following, we discuss when $f_{\epsilon}$ or $\nabla f_{\epsilon}$ are finitely defined.
For this, let us cover a few preliminaries:

Let a function $f(x)$ be called $\mathcal{O}(b(x))$ bounded if there exist $c,v\in\mathcal{O}(b(x))$ and $\bar{c},\bar{v}\in\mathbb{R}$ such that

$$ $\bar{c}+c(x)\leq f(x)\leq\bar{v}+v(x)\qquad\forall x\,.$ (83) $$

For example, a function may be called polynomially bounded (wrt. a polynomial $b(x)$) if (but not only if) $-b(x)\leq f(x)\leq b(x)$.

Moreover, let a density $\mu$ with support $\mathbb{R}$ be called decaying faster than $b(x)$ if $\mu(x)\in{o}(b(x))$.
For example, the standard Gaussian density decays faster than $\exp(-|x|)$, i.e., $\mu(x)\in o(\exp(-|x|))$.
Additionally, we can say that Gaussian density decays at rate $\exp(-x^{2})$, i.e., $\mu(x)\in\theta(\exp(-x^{2}))$.

Now, we can formally characterize finite definedness of $f_{\epsilon}$ and $\nabla f_{\epsilon}$:

###### Lemma 9 (Finite Definedness of $f_{\epsilon}$ ) .

$f_{\epsilon}$ is finitely defined if there exists an increasing function $b(\cdot)$ such that

$$ $f(x)\text{ is bounded by }\mathcal{O}(b(x))\qquad\text{and}\qquad\mu(\epsilon)\in\mathcal{O}(1/b(\epsilon+\alpha\epsilon)/\epsilon^{(1+\alpha)})$ (84) $$

for some $\alpha>0$.

###### Proof.

To show that $f_{\epsilon}$ exists, we need to show that

$$ $\int_{\mathbb{R}}\big|f(x+\epsilon)\cdot\mu(\epsilon)\big|\,d\epsilon$ (85) $$

is finite for all $x$.
Let $\tilde{f}$ be an absolutely upper bound of $f$, and w.l.o.g. let us choose $\tilde{f}(y)=b(y)+\bar{b}$ with $b(y)>1$ for $y\in\mathbb{R}$.
Further, as per the assumptions $\mu(\epsilon)<\frac{1}{\epsilon^{(1+\alpha)}\cdot b(\epsilon+\alpha\epsilon)}\cdot w$ for all $\epsilon<\omega_{1}$ as well as all $\epsilon>\omega_{2}$ for some $w,\omega_{1},\omega_{2}$.
Let us restrict $\omega_{1},\omega_{2}$ to $\omega_{1}<-|x|/\alpha$ and $\omega_{2}>|x|/\alpha$.
It is trivial to see that

$$ $\int_{\omega_{1}}^{\omega_{2}}\big|f(x+\epsilon)\cdot\mu(\epsilon)\big|\,d\epsilon<\infty\,.$ (86) $$

W.l.o.g., let us consider the upper remainder:

$$ $\displaystyle\int_{\omega_{2}}^{\infty}\big|f(x+\epsilon)\cdot\mu(\epsilon)\big|\,d\epsilon$ $\displaystyle\leq\int_{\omega_{2}}^{\infty}\big|\tilde{f}(x+\epsilon)\cdot\mu(\epsilon)\big|\,d\epsilon$ (87) $\displaystyle\leq\int_{\omega_{2}}^{\infty}\left|(b(x+\epsilon)+\bar{b})\cdot\frac{1}{\epsilon^{(1+\alpha)}\cdot b(\epsilon+\alpha\epsilon)}\cdot w\right|\,d\epsilon$ (88) $\displaystyle=\int_{\omega_{2}}^{\infty}\left|\left(\frac{{b(x+\epsilon)}}{\epsilon^{(1+\alpha)}\cdot{b(\epsilon+\alpha\epsilon)}}+\frac{\bar{b}}{\epsilon^{(1+\alpha)}\cdot b(\epsilon+\alpha\epsilon)}\right)\cdot w\right|\,d\epsilon$ (89) $\displaystyle\leq\int_{\omega_{2}}^{\infty}\left|\left(\frac{b(x+\epsilon)}{\epsilon^{(1+\alpha)}\cdot b(\epsilon+|x|)}+\frac{\bar{b}}{\epsilon^{(1+\alpha)}\cdot b(\epsilon+\alpha\epsilon)}\right)\cdot w\right|\,d\epsilon$ (90) $\displaystyle\leq\int_{\omega_{2}}^{\infty}\left|\left(\frac{1}{\epsilon^{(1+\alpha)}}+\frac{\bar{b}}{\epsilon^{(1+\alpha)}\cdot b(\epsilon+\alpha\epsilon)}\right)\cdot w\right|\,d\epsilon$ (91) $\displaystyle<\int_{\omega_{2}}^{\infty}\left|\frac{1}{\epsilon^{(1+\alpha)}}+\frac{\bar{b}}{\epsilon^{(1+\alpha)}}\right|\,d\epsilon\cdot w$ (92) $\displaystyle=\int_{\omega_{2}}^{\infty}\left|\frac{1}{\epsilon^{(1+\alpha)}}\right|\,d\epsilon\cdot w\cdot(1+\bar{b})<\infty\,.$ (93) $$

That $\int_{\omega_{2}}^{\infty}\frac{1}{\epsilon^{(1+\alpha)}}\,d\epsilon$ is finite for the step in (93) can be shown via

$$ $\int_{\omega_{2}}^{\infty}\frac{1}{\epsilon^{(1+\alpha)}}\,d\epsilon=\int_{\omega_{2}}^{\infty}\epsilon^{-1-\alpha}\,d\epsilon=\left[-\frac{1}{\alpha}\epsilon^{-\alpha}\right]_{\omega_{2}}^{\infty}=\left[-\frac{1}{\alpha}\lim_{\epsilon\to\infty}\epsilon^{-\alpha}+\frac{1}{\alpha}\omega_{2}^{-\alpha}\right]=\frac{1}{\alpha}\omega_{2}^{-\alpha}\,.$ $$

The same can be shown analogously for the integral $\int_{-\infty}^{\omega_{1}}$.
This completes the proof.
∎

###### Lemma 10 (Finite Definedness of $\nabla f_{\epsilon}$ ) .

$\nabla f_{\epsilon}$ is finitely defined if there exists an increasing function $b(\cdot)$ such that

$$ $f(x)\text{ is bounded by }\mathcal{O}(b(x))\qquad\text{and}\qquad\big|\mu(\epsilon)\cdot\nabla_{\!\epsilon}{-}\log\mu(\epsilon)\big|\in\mathcal{O}(1/b(\epsilon+\alpha\epsilon)/\epsilon^{(1+\alpha)})$ (94) $$

for some $\alpha>0$.

###### Proof.

The proof of Lemma 9 also applies here, but with $\big|\mu(\epsilon)\cdot\nabla_{\!\epsilon}{-}\log\mu(\epsilon)\big|<\frac{1}{\epsilon^{(1+\alpha)}\cdot b(\epsilon+\alpha\epsilon)}\cdot w$ for all $\epsilon<\omega_{1}$ as well as all $\epsilon>\omega_{2}$ for some $w,\omega_{1},\omega_{2}$.
∎

###### Example 11 (Cauchy and the Identity) .

Let $\mu$ be the density of a Cauchy distribution and let $f(x)=x$.
The tightest $b$ for $f(x)\in\mathcal{O}(b(x))$ is $b(x)=x$.

We have $\mu(\epsilon)\in\theta(1/\epsilon^{2})$ and thus $\mu(\epsilon)\notin o(1/\epsilon^{2})$.
$f_{\epsilon}$, i.e., the mean of the Cauchy distribution is not defined.

However, its gradient $\nabla f_{\epsilon}=1$ is indeed finitely defined.
In particular, we can see that

$$ $\mu(\epsilon)\cdot\nabla_{\!\epsilon}{-}\log\mu(\epsilon)=\frac{2\epsilon}{\pi\cdot(1+\epsilon^{2})\cdot(1+\epsilon^{2})}\in\mathcal{O}(1/\epsilon^{3})\,.$ (95) $$

This is an intriguing property of the Cauchy distribution (or other edge cases) where $f_{\epsilon}$ is undefined whereas $\nabla f_{\epsilon}$ is finitely and well-defined.
In practice, we often only require the gradient for stochastic gradient descent, which means that we often only require $\nabla f_{\epsilon}$ to be well defined and do not necessarily need to evaluate $f_{\epsilon}$ depending on the application.

Additional discussions for the Cauchy distribution and an extension of stochastic smoothing to the $k$-sample median can be found in the next appendix.

###### Lemma 9 (Finite Definedness of f ϵ f_{\epsilon} ) .

###### Proof.

###### Lemma 10 (Finite Definedness of ∇ f ϵ \nabla f_{\epsilon} ) .

###### Proof.

###### Example 11 (Cauchy and the Identity) .

## Appendix C Stochastic Smoothing, Medians, and the Cauchy Distribution

In this section, we provide a discussion of a special case of stochastic smoothing with the Cauchy distribution, and provide an extension of stochastic smoothing to the $k$-sample median.
This becomes important if the range of $f$ is not subset of a compact set, and thus $\mathbb{E}_{\epsilon\sim\mu}\big[f(x+\epsilon)\big]$ becomes undefined for some choice of distribution $\mu$.
For example, for $f(x+\epsilon)=\epsilon$ and $\mu$ being the density of a Cauchy distribution, $\mathbb{E}_{\epsilon\sim\mu}\big[f(x+\epsilon)\big]=\mathbb{E}_{\epsilon\sim\mu}\big[\epsilon\big]$ is undefined.
Nevertheless, even in this case, the gradient estimators discussed in this paper for $\nabla_{\!x}\mathbb{E}_{\epsilon\sim\mu}\big[f(x+\epsilon)\big]$ remain well defined.
This is practically relevant because $\mathbb{E}_{\epsilon\sim\mu}\big[f(x+\epsilon)\big]$ does not need to be finitely defined as long as $\nabla_{\!x}\mathbb{E}_{\epsilon\sim\mu}\big[f(x+\epsilon)\big]$ is well defined.
Further, we remark that the undefinedness of $\mathbb{E}_{\epsilon\sim\mu}\big[f(x+\epsilon)\big]$ requires the range of $f$ to be unbounded, i.e., if there exists a maximum / minimum possible output, then it is well defined.
Moreover, there exist $f$ with unbounded range for which $\mathbb{E}_{\epsilon\sim\mu}\big[f(x+\epsilon)\big]$ also remains well defined.

To account for cases where $\mathbb{E}_{\epsilon\sim\mu}\big[f(x+\epsilon)\big]$ may not be well defined or not a robust statistic, we introduce an extension of smoothing to the median.
We begin by defining the $k$-sample median.

###### Definition 12 ( $k$ -Sample Median) .

For a number of samples $k>1$, and a distribution $\zeta$, we say that

$$ $\mathbb{E}_{z_{1},z_{2},...,z_{k}\sim\zeta}\Big[\operatorname{median}\left\{z_{1},z_{2},...,z_{k}\right\}\Big]$ (96) $$

is the $k$-sample median. For multivariate distributions, let $\operatorname{median}$ be the per-dimension median.

Indeed, for $k\geq 5$, the $k$-sample median estimator is shown to have finite variance for the Cauchy distribution (Theorem 3 and Example 2 in []), which implies a well defined $k$-sample median.
Moreover, for any distribution with a density of the median bounded away from 0, the first and second moments are guaranteed to be finitely defined for sufficiently large $k$.
This is important for non-trivial $f$ with $f(\epsilon)\neq\epsilon$ for at least one $\epsilon$ with $\epsilon\sim\mu$, which implies $\zeta\neq\mu$.
Thus, rather than computing and differentiating the expected value, we can differentiate the $k$-sample median.

###### Lemma 13 (Differentiation of the $k$ -Sample Median) .

With the $k$-sample median smoothing as

$$ $f^{(k)}_{\epsilon}(x)=\mathbb{E}_{\epsilon_{1},...,\epsilon_{k}\sim\mu}\Big[\operatorname{median}\left\{f(x+\epsilon_{1}),...,f(x+\epsilon_{k})\right\}\Big]\,,$ (97) $$

we can differentiate $f^{(k)}_{\epsilon}(x)$ as

$$ $\nabla_{\!x}f^{(k)}_{\epsilon}(x)=\mathbb{E}_{\epsilon_{1},...,\epsilon_{k}\sim\mu}\Big[f(x+\epsilon_{r(\epsilon)})\cdot\nabla_{\!\epsilon_{r(\epsilon)}}{-}\log\mu(\epsilon_{r(\epsilon)})\Big]$ (98) $$

where $r(\epsilon)$ is the arg-median of the set $\{f(x+\epsilon_{1}),...,f(x+\epsilon_{k})\}$, which is equivalent to the implicit definition via $f(x+\epsilon_{r(\epsilon)})=\operatorname{median}\left\{f(x+\epsilon_{1}),...,f(x+\epsilon_{k})\right\}$.

###### Proof.

We denote $\epsilon_{1:k}\sim\mu^{(1:k)}$ such that $\epsilon_{1:k}=\big[\epsilon^{\top}_{1},...,\epsilon^{\top}_{k}\big]^{\top}$ and $\epsilon_{i}\sim\mu\,\,\forall i\in\{1,...,k\}$.

$$ $\displaystyle\nabla_{\!x}f^{(k)}_{\epsilon}(x)$ $\displaystyle=\nabla_{\!x}\mathbb{E}_{\epsilon_{1},...,\epsilon_{k}\sim\mu}\Big[\operatorname{median}\left\{f(x+\epsilon_{1}),...,f(x+\epsilon_{k})\right\}\Big]$ (99) $\displaystyle=\nabla_{\!x}\mathbb{E}_{\epsilon_{1:k}\sim\mu^{(1:k)}}\Big[\operatorname{median}\left\{f(x+\epsilon_{1}),...,f(x+\epsilon_{k})\right\}\Big]$ (100) $\displaystyle=\nabla_{\!x}\int_{\mathbb{R}^{n\cdot k}}\operatorname{median}\left\{f(x+\epsilon_{1}),...,f(x+\epsilon_{k})\right\}\cdot\mu^{(1:k)}(\epsilon_{1:k})\,d\epsilon_{1:k}$ (101) $\big(x_{1},...,x_{k}=x\big)$ $\displaystyle=\sum_{j=1}^{k}\nabla_{\!x_{j}}\int_{\mathbb{R}^{n\cdot k}}\operatorname{median}\left\{f(x_{1}+\epsilon_{1}),...,f(x_{k}+\epsilon_{k})\right\}\cdot\mu^{(1:k)}(\epsilon_{1:k})\,d\epsilon_{1:k}$ (102) $$

As a shorthand, we abbreviate the indicator $\mathbb{1}_{f(x_{j}+\epsilon_{j})=\operatorname{median}\left\{f(x_{1}+\epsilon_{1}),...,f(x_{k}+\epsilon_{k})\right\}}$ as $\mathbb{1}_{j,\epsilon_{1:k}}$ and abbreviate $\mathbb{1}_{f(u_{j})=\operatorname{median}\left\{f(u_{1}),...,f(u_{k})\right\}}$ as $\mathbb{1}_{j,u_{1:k}}$:

$$ $\displaystyle\nabla_{\!x}f^{(k)}_{\epsilon}(x)$ $\displaystyle=\sum_{j=1}^{k}\nabla_{\!x_{j}}\int_{\mathbb{R}^{n\cdot k}}f(x_{j}+\epsilon_{j})\cdot\mathbb{1}_{j,\epsilon_{1:k}}\cdot\mu^{(1:k)}(\epsilon_{1:k})\,d\epsilon_{1:k}$ (103) $\displaystyle=\sum_{j=1}^{k}\nabla_{\!x_{j}}\int_{\mathbb{R}^{n\cdot k}}f(u)\cdot\mathbb{1}_{j,u{1:k}}\cdot\mu^{(1:k)}(u_{1:k}-x)\,du_{1:k}$ (104) $\displaystyle=\sum_{j=1}^{k}\int_{\mathbb{R}^{n\cdot k}}f(u)\cdot\mathbb{1}_{j,u{1:k}}\cdot\nabla_{\!x_{j}}\mu^{(1:k)}(u_{1:k}-x)\,du_{1:k}$ (105) $\displaystyle=\sum_{j=1}^{k}\int_{\mathbb{R}^{n\cdot k}}f(x+\epsilon_{j})\cdot\mathbb{1}_{j,\epsilon_{1:k}}\cdot-\nabla_{\!\epsilon_{j}}\mu^{(1:k)}(\epsilon_{1:k})\,d\epsilon_{1:k}$ (106) $$

We have

$$ $\displaystyle\nabla_{\!\epsilon_{j}}\mu^{(1:k)}(\epsilon_{1:k})=\mu^{(1:k)}(\epsilon_{1:k})\cdot\nabla_{\!\epsilon_{j}}\log\mu^{(1:k)}(\epsilon_{1:k})=\mu^{(1:k)}(\epsilon_{1:k})\cdot\nabla_{\!\epsilon_{j}}\log\mu(\epsilon_{j})\,.$ (107) $$

Thus,

$$ $\displaystyle\nabla_{\!x}f^{(k)}_{\epsilon}(x)$ $\displaystyle=\sum_{j=1}^{k}\int_{\mathbb{R}^{n\cdot k}}f(x+\epsilon_{j})\cdot\mathbb{1}_{j,\epsilon_{1:k}}\cdot-\mu^{(1:k)}(\epsilon_{1:k})\cdot\nabla_{\!\epsilon_{j}}\log\mu(\epsilon_{j})\,d\epsilon_{1:k}$ (108) $\displaystyle=\int_{\mathbb{R}^{n\cdot k}}\sum_{j=1}^{k}\left[\mathbb{1}_{j,\epsilon_{1:k}}\cdot f(x+\epsilon_{j})\cdot\nabla_{\!\epsilon_{j}}{-}\log\mu(\epsilon_{j})\right]\cdot\mu^{(1:k)}(\epsilon_{1:k})\,d\epsilon_{1:k}$ (109) $$

Indicating the choice of median in dependence of $\epsilon_{1:k}$, we define $r(\epsilon_{1:k})$ s.t. $\mathbb{1}_{r(\epsilon_{1:k}),\epsilon_{1:k}}=1$. Thus,

$$ $\displaystyle\nabla_{\!x}f^{(k)}_{\epsilon}(x)$ $\displaystyle=\int_{\mathbb{R}^{n\cdot k}}f(x+\epsilon_{r(\epsilon_{1:k})})\cdot\nabla_{\!\epsilon_{r(\epsilon_{1:k})}}{-}\log\mu(\epsilon_{r(\epsilon_{1:k})})\cdot\mu^{(1:k)}(\epsilon_{1:k})\,d\epsilon_{1:k}$ (110) $\displaystyle=\mathbb{E}_{\epsilon_{1:k}\sim\mu^{(1:k)}}\Big[f(x+\epsilon_{r(\epsilon_{1:k})})\cdot\nabla_{\!\epsilon_{r(\epsilon_{1:k})}}{-}\log\mu(\epsilon_{r(\epsilon_{1:k})})\Big]$ (111) $$

This concludes the proof.
∎

Empirically, we can estimate $\nabla_{\!x}f^{(k)}_{\epsilon}(x)$ for $s$ propagated samples ($s>k$) without bias as

$$ $\displaystyle\nabla_{\!x}f^{(k)}_{\epsilon}(x)\,\triangleq\,\sum_{i=1}^{s}\,\Big[q_{i}\cdot f(x+\epsilon_{i})\cdot\nabla_{\!\epsilon_{i}}{-}\log\mu(\epsilon_{i})\Big]\qquad\epsilon_{1},...,\epsilon_{s}\sim\mu$ (112) $$

where $q_{i}$ is the probability of $f(x+\epsilon_{i})$ being the median in a subset of $k$ samples, i.e., under uniqueness of $g_{i}$s, we have

$$ $\displaystyle q_{i}=\frac{\displaystyle\sum_{\{h_{1},...,h_{k}\}\subset\{g_{1},...,g_{s}\}}\kern-20.00003pt\mathbb{1}\big(g_{i}=\operatorname{median}\{h_{1},...,h_{k}\}\big)}{\displaystyle\binom{s}{k}}\qquad\qquad g_{i}:=f(x+\epsilon_{i})\,.$ (113) $$

We remark that, in case of non-uniqueness, it is adequate to split the probability among the candidates; however, under non-discreteness assumptions on $f$ (density of $\zeta<\infty$, the converse typically implies the range of $f$ being a subset of a compact set), this almost surely (with probability 1) does not occur.

We have shown that the $k$-sample median $f^{(k)}_{\epsilon}(x)$ is differentiable and demonstrated an unbiased gradient estimator for it.
A straightforward extension for the case of $f$ being differentiable is differentiating through the median via a $k\to\infty$-sample median, e.g., via setting $s=k^{2}$.
The $k\to\infty$ extension for differentiating through the median itself requires $f$ being differentiable because, for discontinuous $f$, $f^{(k)}_{\epsilon}(x)$ is differentiable only for $k<\infty$.
(As an illustration, the median of the Heaviside function under a symmetric perturbation $\mu$ with density at 0 bounded away from 0 is the exactly the Heaviside function.)

###### Definition 12 ( k k -Sample Median) .

###### Lemma 13 (Differentiation of the k k -Sample Median) .

###### Proof.

## Appendix D Experimental Details

##### MNIST Sorting Benchmark Experiments

We train for 100 000 steps at a learning rate of 0.001 with the Adam optimizer using a batch size of 100.
Following the requirements of the benchmark, we use the same model as previous works [, , ].
That is, two convolutional layers with a kernel size of $5\times 5$, 32 and 64 channels respectively, each followed by a ReLU and MaxPool layer; after flattening, this is followed by a fully connected layer with a size of 64, a ReLU layer, and a fully connected output layer mapping to a scalar.
For each distribution and number of samples, we choose the optimal $\gamma\in\{1,1/3,0.1\}$.

##### Warcraft Shortest-Path Benchmark Experiments

Following the established protocol [], we train for 50 epochs with the Adam optimizer at a batch size of 70 and an initial learning rate of 0.001.
The learning rate decays by a factor of 10 after 30 and 40 epochs each. The model is the first block of ResNet18.
The hyperparameter $\gamma=1/\beta$ as specified in Figures 13 and 14.

##### Utah Teapot Camera Pose Optimization Experiments

We initialize the pose to be perturbed by angles uniformly sampled from $[15^{\circ},75^{\circ}]$.
The ground truth orientation is randomly sampled from the sphere of possible orientations.
The ground truth camera angle is $20^{\circ}$, and the ground truth camera distance is uniformly sampled from $[2.5,4]$.
The initial camera distance is sampled as being uniformly offset by $[-0.5,6]$, thus the feasible set of initial camera distance guesses lies in $[2,10]$.
The initial camera angle is uniformly sampled from $[10^{\circ},30^{\circ}]$.
We optimize for 1 000 steps with the Adam optimizer [$(\beta_{1},\beta_{2})=(0.5,0.99)$] and the CosineAnnealingLR scheduler with an initial learning rate of $0.3$.
We schedule the diagonal of $\mathbf{L}$ to decay exponentially from $[0.1,5^{\circ},5^{\circ},0.25^{\circ}]\cdot 10^{0.75}$ to $[0.1,5^{\circ},5^{\circ},0.25^{\circ}]\cdot 10^{-1.75}$ (the dimensions are camera distance, 2 pose angles, and the camera angle).
As discussed, the success criterion is finding the angle within $5^{\circ}$ of the ground truth angle. There is typically no local minimum within $5^{\circ}$ and it is a reliable indicator for successful alignment.

##### Differentiable Cryo-Electron Tomography Experiments

The ground truth values of the parameters are set to $300$ kV for acceleration voltage, $3$ mm for the focal length, and the ground truth sample specimen is centered as $(x,y)=(0,0)$ nm units. For reporting the RMSE metric, the acceleration voltages are normalized by a factor of $100$ to ensure that all parameters vary over commensurate ranges. For the 2-parameter optimization, the feasible set of acceleration voltage varied over a range of $[0,1000]$ kV and the feasible set of the specimen’s $x$-position varied over the range $[-5,5]$. For the 4-parameter optimization, the feasible set of acceleration voltage varied over a range of $[0,600]$ kV, the focal length ranges over $[0,6]$ mm, the $x$- and $y$-positions range over $[-3,3]$. We use the Adam optimizer for both experiments, with [$(\beta_{1},\beta_{2})=(0.5,0.9)$]. For the MC Search baseline, we generate sets of $n$ uniform random points in the feasible region of the parameters, generate micrographs for these random parameter tuples using the TEM simulator [], and identify the parameter tuple in the set having the lowest mean squared error with respect to the ground truth image. The RMSE between this parameter tuple and the ground truth parameters is the metric for the specific set of $n$ randomly generated values. This is repeated $20$ times to obtain the mean and standard deviation of the RMSE metric at that $n$.

### D.1 Assets

List of assets:

- •
The sixth platonic solid (aka. Teapotahedron or Utah tea pot) []   [License N/A]
- •
Multi-digit MNIST [], which builds on MNIST []   [MIT License / CC License]
- •
Warcraft shortest-path data set []   [MIT License]
- •
PyTorch []   [BSD 3-Clause License]
- •
TEM-simulator []   [GNU General Public License]

### D.2 Runtimes

The runtimes for sorting and shortest-path experiments are for one full training on 1 GPU.
The pose optimization experiment runtimes are the total time for all 768 seeds on 1 GPU.
For the TEM-simulator, we report the CPU time per simulation sample, which is the dominant and only the measureable component of the total optimization routine time.
The choice of distribution, covariate, and choice of variance reduction does not have a measurable effect on training times.

- •
MNIST Sorting Benchmark Experiments [1 Nvidia V100 GPU]
–
Training w/ 256 samples:
65 min
–
Training w/ 1 024 samples:
67 min
–
Training w/ 2 048 samples:
68 min
–
Training w/ 8 192 samples:
77 min
–
Training w/ 32 768 samples:
118 min
- •
Warcraft Shortest-Path Benchmark Experiments [1 Nvidia V100 GPU]
–
Training w/ 10 samples:
9 min
–
Training w/ 100 samples:
19 min
–
Training w/ 1 000 samples:
26 min
–
Training w/ 10 000 samples:
101 min
- •
Utah Teapot Camera Pose Optimization Experiments [1 Nvidia A6000 GPU]
–
Optimization on 768 seeds w/ 16 samples: 25 min
–
Optimization on 768 seeds w/ 64 samples: 81 min
–
Optimization on 768 seeds w/ 256 samples: 362 min
- •
Differentiable Cryo-Electron Tomography Experiments [CPU: 44 Intel Xeon Gold 5118]
–
Simulator time per sample on 1 CPU core: 67 sec

## Appendix E Additional Experimental Results

**Table 3: Extension of Table 2 with additional numbers of samples and standard deviations.**
| Baselines |  | Neu.S. | Soft.S. | L. DSN | C. DSN | E. DSN | OT. S. |
| --- | --- | --- | --- | --- | --- | --- | --- |
| — |  | 71.3 | 70.7 | 77.2 | 84.9 | 85.0 | 81.1 |
| Sampling | #s | Gauss. | Logis. | Gumbel | Cauchy | Laplace | Trian. |
| vanilla | 256 | 82.3$\pm$2.0 | 82.8$\pm$0.9 | 79.2$\pm$9.7 | 68.1$\pm$19.3 | 82.6$\pm$0.8 | 81.3$\pm$1.2 |
| best (cv) | 256 | 83.1$\pm$1.6 | 82.7$\pm$1.8 | 81.6$\pm$3.6 | 55.6$\pm$13.3 | 83.7$\pm$0.8 | 82.7$\pm$1.1 |
| vanilla | 1024 | 81.3$\pm$9.1 | 83.7$\pm$0.7 | 82.0$\pm$1.6 | 68.5$\pm$24.8 | 80.6$\pm$9.0 | 82.8$\pm$1.0 |
| best (cv) | 1024 | 83.9$\pm$0.6 | 84.0$\pm$0.5 | 84.2$\pm$0.6 | 73.0$\pm$12.6 | 84.3$\pm$0.6 | 82.4$\pm$1.6 |
| vanilla | 2048 | 84.1$\pm$0.6 | 83.6$\pm$0.8 | 84.0$\pm$0.5 | 75.7$\pm$11.6 | 83.8$\pm$0.7 | 83.2$\pm$0.6 |
| best (cv) | 2048 | 84.2$\pm$0.5 | 84.2$\pm$0.6 | 84.6$\pm$0.4 | 82.0$\pm$2.2 | 84.8$\pm$0.5 | 83.4$\pm$0.5 |
| vanilla | 8192 | 84.0$\pm$0.6 | 84.2$\pm$0.8 | 84.0$\pm$0.6 | 83.6$\pm$1.0 | 83.9$\pm$1.0 | 83.6$\pm$0.7 |
| best (cv) | 8192 | 84.4$\pm$0.6 | 84.5$\pm$0.5 | 84.1$\pm$0.7 | 84.3$\pm$0.5 | 84.3$\pm$0.4 | 83.7$\pm$0.4 |
| vanilla | 32768 | 84.2$\pm$0.5 | 84.1$\pm$0.4 | 84.5$\pm$0.7 | 84.9$\pm$0.5 | 84.4$\pm$0.5 | 83.4$\pm$0.8 |
| best (cv) | 32768 | 84.4$\pm$0.4 | 84.4$\pm$0.4 | 84.8$\pm$0.5 | 85.1$\pm$0.4 | 84.4$\pm$0.4 | 84.0$\pm$0.3 |

Figure: Figure 11: Average $L_{2}$ norms between ground truth (oracle) and estimated gradient for different numbers of elements to sort and rank $n$, and different distributions. Each plot compares different variance reduction strategies as indicated in the legend to the right of the caption. Darker is better (smaller values). Colors are only comparable within each subplot. We use $1\,024$ samples, except for Cartesian and $n=3$ where we use $10^{3}=1\,000$ samples.
Refer to caption: x31.png

Figure: Figure 12: Utah teapot camera pose optimization with smoothing of the loss, compared to Figure 8, which performs smoothing of the algorithm. Smoothing the algorithm is consistently better, with the largest effect for larger numbers of samples. Results averaged over 768 seeds.
Refer to caption: x32.png

Figure: Figure 13: Warcraft shortest-path experiment. Left: 10 samples. Right: 100 samples. Averaged over 5 seeds. Brighter is better. Values between subplots are comparable. The displayed range is $[70\%,96.5\%]$.
Refer to caption: x33.png

Figure: Figure 14: Warcraft shortest-path experiment. Left: 1 000 samples. Right: 10 000 samples. Averaged over 5 seeds. Brighter is better. Values between subplots are comparable. The displayed range is $[70\%,96.5\%]$.
Refer to caption: x36.png

Figure: Figure 15: Cryo-Electron Tomography Experiments: RMSE with respect to Ground Truth parameters for different number of parameters optimized and for different number of samples per optimization step: (Top Left) 2-parameters & number of samples=9, (Top Right) 2-parameters & number of samples=25, (Bottom Left) 2-parameters & number of samples=36, (Bottom Right) 4-parameters. No marker lines correspond to Gaussian, $\times$ corresponds to Laplace, and $\triangle$ corresponds to Triangular distributions. Ascertaining optimal parameters with minimal evaluations is important not just for high resolution imaging, but also to minimize radiation damage to the specimen. In this light, of the covariate choices, LOO generally leads to best improvement and none consistently leads to deterioration in performance. The Laplace and Triangular distributions lead to best performance. For the Gaussian distribution, Cartesian RQMC is generally exhibiting best results.
Refer to caption: x39.png

**Table 4: Individual absolute values from the variance simulations for differentiable sorting in Figure 4. The minimum and values within 1% of the minimum are indicated as bold.**
|  | none | $f(x)$ | LOO | none | $f(x)$ | LOO |
| --- | --- | --- | --- | --- | --- | --- |
|  | regular | antithetic |  |  |  |  |
| MC | 0.0084 | 0.0079 | 0.0046 | 0.0055 | 0.0054 | 0.0053 |
| QMC (lat.) | 0.0029 | 0.0030 | 0.0030 | 0.0036 | 0.0036 | 0.0036 |
| RQMC (l.) | 0.0030 | 0.0030 | 0.0030 | 0.0036 | 0.0035 | 0.0036 |
| RQMC (c.) | 0.0012 | 0.0013 | 0.0012 | 0.0014 | 0.0014 | 0.0014 |

**Table 5: Individual absolute values from the variance simulations for differentiable shortest-paths in Figure 4. The minimum and values within 1% of the minimum are indicated as bold.**
|  | none | $f(x)$ | LOO | none | $f(x)$ | LOO |
| --- | --- | --- | --- | --- | --- | --- |
|  | regular | antithetic |  |  |  |  |
| MC | 1330.01 | 4.17 | 4.17 | 8.32 | 8.32 | 8.34 |
| QMC (lat.) | 4.04 | 4.04 | 4.04 | 8.04 | 8.04 | 8.07 |
| RQMC (l.) | 4.25 | 4.05 | 4.05 | 8.10 | 8.09 | 8.12 |

## Peer Reviews

confidence: 5
contribution: 1
initial_rating: 3
presentation: 4
questions: .
soundness: 4
strengths: .
summary: This work presents methods for stochastic smoothing of non-differentiable functions, motivated by applications where there is the need to optimize non-differentiable functions. Some experiments are provided.
weaknesses: The novelty is insufficient. The idea of smoothing using Gaussian perturbations is a well-established technique in 0-th order optimization or non-smooth optimization. The authors present further general smoothing conditions, requiring that the density function be absolutely continuous, but these derivations are mostly standard exercises in mathematical analysis. The variance reduction techniques are also not novel. Overall, my view is that there is very little theoretical novelty.

The authors provide applications and experiments, but I do not think they are strong enough to warrant publication on the basis of the experiments.

---

confidence: 3
contribution: 3
initial_rating: 8
presentation: 4
questions: Can this method be applied to problems where the $ L^1 $ loss is minimized?
soundness: 3
strengths: The paper is very well written, with ideas clearly conveyed. Mathematical descriptions are precise.

I like the fact that multiple density functions are considered and discussed.

The discussion on variance reduction is insightful.

Many different applications are presented. And at the end, some recommendations are made about the "optimal" densities, as well as on the sampling methods.
summary: This paper presents improved results on gradient estimation by stochastic smoothing. In particular, the density function is assumed to be absolutely continuous as opposed to being differentiable. The authors discuss three ways to reduce variance of the estimated gradients. Finally, the paper presents application of this smoothing on various problem where gradients of non-differentiable functions are required.

This paper builds on existing ideas, but it does present some strong and comprehensive results that are both insightful.
weaknesses: None.

---

confidence: 4
contribution: 1
initial_rating: 3
presentation: 2
questions: 1. Can the authors comment on how their theoretical or experimental results might be applicable for discrete-to-continuous relaxations? Particularly, how would quantities such as lipschitz smoothness, convexity, or general ease of optimization be affected by different choices of distributions?

2. Can the authors comment on the run-time of their algorithms compared to the baselines in Section 4?
soundness: 2
strengths: The experiments in Section 4 seem quite interesting, particularly the use of smoothing for relaxing discrete functions to continuous functions.
summary: The authors study stochastic gradients, with applications in stochastic smoothing and variance reduction of gradient estimation. 

I recommend rejection of this paper because I feel there is a lack of novelty and impact for machine learning.

I present a more detailed summary of the paper along with my criticisms under the weakness section.
weaknesses: The goals of the paper are somewhat diverse, so I will go over each result and state my criticisms.

1. In Lemma 3, the authors present a result for smoothing an objective $f$ using a distribution $\mu$ that is possibly non-differentiable. However, the non-differentiable set for $\mu$ is assumed to be $0$-measure. In any learning application I can think of, there seems to be no meaningful difference between "differentiable" and "almost everywhere differentiable". For instance, one can always convolve $\mu$ with a tiny Gaussian, of variance less than machine precision, and immediately get a differentiable density, while not affecting the outcome in any meaningful way. There are also no quantitative analysis which indicate that choosing a non-differentiable $\mu$ would offer any advantages. 
2. Theorem 7 discusses smoothing with an additional rescaling matrix involved. The result appears unsurprising, and the proof appears to use just calculus. I think the authors want to apply this to finding a optimal rescaling matrix $\mathbf{L}$, but the consequences of Theorem 7 are never discussed.
3. The authors discuss numerous variance reduction techniques in Section 2.1, but there is no attempt at any quantitative and theoretical comparison of these approaches. The authors did include some simple synthetic experiments, but it is questionable how well the observations will generalize to realistic scenarios.
4. In Section 2.3, the authors discuss smoothing of the algorithm vs the objective. I am not convinced by the setup of $\ell(h(y))$ where $\ell$ is the loss, $h$ is an algorithm, and $y$ is the model output. 
5. In Section 4, the authors present a series of experiments that compare the effectiveness of different distributions and different variance reduction techniques. I am not sure how well the conclusions here generalize.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 4
questions: Is it possible to demonstrate the advantage of your smoothing approach over the previous works in theory, not only in practice, in some way? For instance, by incorporating the proposed gradient estimator into the SGD or Adam optimizer and providing a better convergence result? 

There exist many general variance reduction techniques in the field of optimization (in such algorithms as SVRG, L-SVRG, Page, SAGA, Variance Reduced Adam, SPIDER etc.). Do the proposed variance reduction techniques generalize the techniques from the optimization field or orthogonal to them? Can they be combined together?

Smoothing should introduce some bias to the estimator. There are some works on the analysis of SGD with biased estimators [1, 2]. They both consider smoothing techniques. Does the proposed smoothing technique fit into the frameworks of assumptions in these papers?

[1] Ahmad Ajalloeian, Sebastian U. Stich, On the Convergence of SGD with Biased Gradients
[2] Yury Demidovich, Grigory Malinovsky, Igor Sokolov, Peter Richtárik, A Guide Through the Zoo of Biased SGD
soundness: 3
strengths: The work contains rigorous theorem statements and proofs. The work contains extensive empirical results.
summary: The work addresses gradient estimation for stochastic differentiable relaxations of non-differentiable functions. The developed theory does not require a differentiable density or full support as it was in prior works.
weaknesses: The work considers the problem of estimating the gradient of the relaxation of the algorithm or of the loss function. In Section 2.3 the authors formalize the algorithm smoothing and the objective smoothing. It seems that the paper does not provide the theoretical convergence analysis of the backbone stochastic optimization algorithms, but only introduces the techniques of smoothing and variance reduction and then tests the theoretical findings in practice.

Meta-review: While the paper aims to address interesting problems, several reviewers have found its theoretical contribution, in light of the existing literature, to be rather minimal, with its impact on machine learning deemed questionable. It has been noted that the theoretical derivations are at times trivial, and the results occasionally lack meaningful applications. Furthermore, the numerical examples fall short of being convincing and fail to complement the shortcomings of the theoretical results.

Decision: Reject