1 Introduction↩︎

Deep neural networks have demonstrated remarkable success across a wide range of applications, yet their vulnerability to adversarial perturbations remains a significant obstacle to their deployment in safety-critical systems. These perturbations, often imperceptible to humans, can cause dramatic failures in model predictions—a phenomenon extensively documented in classification tasks [1]–[3]. Consequently, the field has pivoted toward defenses with formal guarantees, with randomized smoothing emerging as a leading scalable method for certifying robustness without making strong assumptions about the model architecture [4]. At a high level, randomized smoothing constructs a smoothed predictor by averaging the base model’s predictions under random perturbations of the input, and then uses the induced smoothness to certify that predictions remain stable within an input neighborhood. However, while certification for classification is well-established, the robustness of regression models—which predict continuous, real-valued outputs—remains less explored. This gap is critical; in domains like autonomous navigation and medical imaging, certifying that a continuous prediction will not deviate beyond a safe margin is a fundamental safety requirement.

Prominent recent approaches for regression, most notably \(\alpha\)-smoothing and RS-Reg [5], [6], operate by defining a probabilistic “safe” region and estimating the probability mass \(p_A\) that the base model’s outputs fall within this region under noise. This formulation is advantageous for its flexibility, as it allows for certification over arbitrary output distributions without assuming parametric forms. However, these approaches effectively reduces the complex behavior of the function to a simple binary outcome (inside or outside the region), discarding critical information regarding the local geometry of the smoothed landscape. Consequently, it often yields loose certificates that do not fully exploit the smoothness properties induced by the noise.

Furthermore, the certification criteria in \(\alpha\)-smoothing and RS-Reg are typically defined by ensuring that for any perturbed input, the smoothed model’s output remains within a neighborhood of the base model’s output [5], [6]. This formulation creates a misalignment between certification and deployment, as the actual model used for inference is the smoothed regressor, not the base model. While it is empirically possible to configure such methods to center on the smoothed prediction, doing so invalidates their theoretical guarantees.

In this work, we introduce a framework that addresses these limitations by leveraging higher-order information to derive certificates explicitly centered on the smoothed prediction. By incorporating the gradient norm of the smoothed regressor alongside the variance of the base output evaluated with smoothing noise, we characterize the “worst-case” base function analytically. This approach yields closed-form analytical bounds on the adversarial shift. Unlike methods that rely on certifying probabilistic acceptance regions, our theoretical bounds are deterministic given the population statistics, though practical certification utilizes high-probability estimates of these values. This distinction allows for significantly tighter radii by exploiting the function’s geometry, while simultaneously resolving the anchor inconsistency by guaranteeing the stability of the actual deployed prediction \(g(x)\). We empirically validate our framework on the MNIST rotation angle prediction task, demonstrating that incorporating gradient information yields significantly tighter robustness radii compared to baselines relying on variance or quantile-based inclusion bounds alone.

2 Related Work↩︎

2.0.0.1 Randomized smoothing and higher-order certification.

Randomized smoothing has established itself as a scalable standard for certifying \(\ell_2\) robustness in deep classifiers [4], [7], often relying on the Neyman-Pearson lemma to derive optimal certification regions from zeroth-order statistics. While subsequent work investigated incorporating first-order information [8], [9], adapting these insights to regression requires a fundamental shift in methodology. Unlike the classification setting, which focuses on decision boundary geometry or likelihood ratios, regression necessitates bounding the magnitude of continuous function variation. To achieve this, we depart from the Neyman-Pearson framework and instead employ a variational analysis to characterize the worst-case base function. This methodological shift is crucial: it allows us to make gradient information practically useful by explicitly incorporating the gradient as a constraint on the worst-case output shift.

2.0.0.2 Approaches to robust regression.

Prior to randomized smoothing, certification for regression was primarily explored through Lipschitz continuity analysis [10]. However, these techniques typically require white-box access to model parameters and scale poorly to deep architectures. In the black-box setting, recent works have attempted to bridge the gap by reducing regression to classification. For instance, [11] discretize the output space for object detection, while [12] reduce regression poisoning to a classification task.

Most relevant to our work, [6] and [5] developed the first direct smoothing frameworks for continuous regression. However, these methods rely on statistical inclusion probabilities rather than function geometry. First, as evidenced by the proof of Theorem 1 in [6], they treat certification as a binary problem: determining whether an output falls within a probabilistic acceptance region based on zeroth-order counts. Second, regarding \(\alpha\)-smoothing [5], while they successfully utilize \(\alpha\)-trimming to stabilize predictions—quantifying this benefit via the regularized beta function—their certificate characterizes local stability solely through the distribution of perturbed evaluations. It exploits the trimming to increase the probability of inclusion, but it does not utilize the gradient to analytically bound the worst-case magnitude of the output shift. In contrast, our approach is geometric: we utilize the gradient norm to directly constrain the variation of the smoothed function itself.

3 Preliminaries↩︎

3.0.0.1 Randomized Smoothing

Randomized smoothing is a post-processing method for improving the robustness of a given predictor. The method is motivated by the observation that large local variation can make a predictor vulnerable to small input perturbations: nearby inputs may receive substantially different predictions. Randomized smoothing reduces this local variation by replacing the predictor at each input with an average of its values over a distribution concentrated near that input. This averaging yields a smoother predictor whose change over a prescribed neighborhood can be precisely quantified, leading to certified robustness guarantees. Prior work has developed this framework primarily for classification; in this paper, we formulate a corresponding smoothing procedure for regression.

3.0.0.2 Notation and Setup

We consider a base regression model (e.g., a neural network) \(f: \mathbb{R}^d \to \mathbb{R}\) that maps an input \(\boldsymbol{x}\in \mathbb{R}^d\) to a \(1\)-dimensional output. Our analysis focuses on the certification of a smoothed regressor \(g: \mathbb{R}^d \to \mathbb{R}\), which is defined as the expectation of the base regressor \(f\) under isotropic Gaussian noise: \[\begin{align} \label{def:g95defn} g(\boldsymbol{x}) = \mathbb{E}[f(\boldsymbol{x}+\boldsymbol{e})], \quad \text{where } \boldsymbol{e}\sim \mathcal{N}(\boldsymbol{0}, \sigma^2 \boldsymbol{I}) \end{align}\tag{1}\]

We denote the \(\ell_2\) norm by \(\|\cdot\|_2\). We use \(\mathbb{E}[\cdot]\) and \(\mathrm{Var}(\cdot)\) to denote the expectation and variance operators respectively. We use \(y \in \mathbb{R}\) to represent the (potentially unknown) ground-truth value associated with an input \(\boldsymbol{x}\).

3.0.0.3 Threat Model and Certification Goal

We operate under the standard \(\ell_2\)-norm threat model. Consider a fixed nominal input \(\boldsymbol{z}\). An adversary may perturb \(\boldsymbol{z}\) by adding any perturbation \(\boldsymbol{\delta}\in \mathbb{R}^d\) as long as the perturbation is bounded, \(\|\boldsymbol{\delta}\|_2 \le R\). The resulting adversarial example is \(\boldsymbol{z}+ \boldsymbol{\delta}\).

Our objective is to certify the local stability of the smoothed regressor \(\boldsymbol{g}\) around \(\boldsymbol{z}\). We aim to find the largest possible radius \(R > 0\) such that for any perturbation \(\boldsymbol{\delta}\) satisfying \(\|\boldsymbol{\delta}\|_2 \le R\), the output of the smoothed function is guaranteed to remain within an acceptable distance, \(\epsilon\), of its prediction at clean input: \[\begin{align} \label{eq:cert95goal} |g(\boldsymbol{z}+ \boldsymbol{\delta}) - g(\boldsymbol{z})| \le \epsilon \end{align}\tag{2}\] This formulation provides a practical guarantee, as the certificate is centered on the model’s own prediction \(g(\boldsymbol{z})\), which is computable at test time.

Formally, certifying this radius requires bounding the worst-case shift of the smoothed function output under perturbation. Denote the worst-case shift as \(\Delta_{\max}\). When underlying base function is assumed to come from a function class \(\mathcal{F}\), \(\Delta_{\max}\) can be seen as the maximizing objective value achieved by this nested maximization problem: \[\begin{align} \label{obj:delta95max95nested95opt} \Delta_{\max}(R) = \max_{\|\boldsymbol{\delta}\|_2 \le R} \max_{f \in \mathcal{F}} \left|\mathbb{E}[f(\boldsymbol{z}+\boldsymbol{\delta}+\boldsymbol{e})] - \mathbb{E}[f(\boldsymbol{z}+\boldsymbol{e})] \right| \end{align}\tag{3}\]

To constrain \(\mathcal{F}\), we compute fundamental summary statistics describing \(f\), such as its mean and variance under gaussian noise. We then optimize over all base functions consistent with our observed statistics (e.g., variance \(C\) and gradient \(\mathbf{G}\)). We certify a radius by finding the largest \(R\) for which \(\Delta_{\max}(R) < \epsilon\).

4 Certification via Worst-Case Analysis↩︎

To derive the certified radius, we solve the nested maximization problem defined in Equation 3 to bound worst case shift \(\Delta_{\max}(R)\). To facilitate the analysis using variational calculus, we switch from the noise variable \(\boldsymbol{e}\) to the input variable \(\boldsymbol{x}\). Let \(\boldsymbol{x}= \boldsymbol{z}+ \boldsymbol{e}\). Then, under the nominal setting, \(\boldsymbol{x}\sim \mathcal{N}(\boldsymbol{z}, \sigma^2 \boldsymbol{I})\). Similarly, for the perturbed setting, \(\boldsymbol{x}\sim \mathcal{N}(\boldsymbol{z}+\boldsymbol{\delta}, \sigma^2 \boldsymbol{I})\). We define the Gaussian densities centered at the nominal and perturbed points as \(p_{\boldsymbol{z}} = \mathcal{N}(\boldsymbol{z}, \sigma^2 \boldsymbol{I})\) and \(p_{\boldsymbol{z}+\boldsymbol{\delta}} = \mathcal{N}(\boldsymbol{z}+\boldsymbol{\delta}, \sigma^2 \boldsymbol{I})\) respectively. With this notation, we can express the expectations in problem 3 as integrals over the input space \(\boldsymbol{x}\): \[\begin{align} \label{eq:objective95general} \Delta_{\max}(R) = \max _{\|\boldsymbol{\delta}\|_2 \leq R}\max _{f \in \mathcal{F}} \left|\mathbb{E}_{p_{\boldsymbol{z}+\boldsymbol{\delta}}}[f(\boldsymbol{x})] - \mathbb{E}_{p_{\boldsymbol{z}}}[f(\boldsymbol{x})]\right| \end{align}\tag{4}\] To obtain a certificate that holds for any base regressor consistent with our observations, we maximize this objective over the space of functions \(f\) satisfying specific empirical constraints.

4.1 Constraint Specifications↩︎

To make the optimization well-posed, we impose constraints based on properties of the base model \(f\) that can be empirically estimated.

We distinguish between two regimes regarding the mean of the function. In the unbounded setting (where the output \(f(\boldsymbol{x})\) can take any real value), the objective function 4 is invariant to constant shifts: replacing \(f(\boldsymbol{x})\) with \(f(\boldsymbol{x}) + c\) changes both expectations by \(c\), leaving the difference \(\Delta(\boldsymbol{\delta})\) unchanged (see Lemma 3 in Appendix 11). Consequently, we do not require knowledge of the specific mean value \(g(\boldsymbol{z})\) to certify the radius in unbounded cases. However, in the bounded setting, the specific value of the mean determines the relative position of the feasible bounds, breaking this invariance. We define optimization problems for both settings below.

4.1.0.1 Unbounded Setting: Variance Constraint.

Since the mean is uninformative in the unbounded setting due to shift invariance, we rely on the variance of the base function to constrain its scale. We assume the variance under the nominal distribution is bounded by an empirically measurable constant \(C\). This yields the following optimization problem: \[\begin{align} \label{opt:var95only} \underset{\|\boldsymbol{\delta}\|\le R, f}{\text{maximize}} \quad &\left|\mathbb{E}_{p_{\boldsymbol{z}+\boldsymbol{\delta}}}[f(\boldsymbol{x})] - \mathbb{E}_{p_{\boldsymbol{z}}}[f(\boldsymbol{x})]\right| \\ \text{subject to} \quad &\mathrm{Var}_{p_{\boldsymbol{z}}}(f(\boldsymbol{x})) \le C \nonumber \end{align}\tag{5}\]

4.1.0.2 Unbounded Setting: Variance and Gradient Constraints.

To capture local sensitivity, we additionally constrain the gradient of the smoothed function. We assume we have an estimate of the gradient vector \(\mathbf{G}\), which fixes both the magnitude and direction of the local trend. We refer to this setting as \((C, G)\) in our experiments. \[\begin{align} \label{opt:var95and95grad} \underset{\|\boldsymbol{\delta}\|\le R, f}{\text{maximize}} \quad &\left|\mathbb{E}_{p_{\boldsymbol{z}+\boldsymbol{\delta}}}[f(\boldsymbol{x})] - \mathbb{E}_{p_{\boldsymbol{z}}}[f(\boldsymbol{x})]\right| \\ \text{subject to} \quad &\mathrm{Var}_{p_{\boldsymbol{z}}}(f(\boldsymbol{x})) \le C \nonumber \\ &\nabla \mathbb{E}_{p_{\boldsymbol{z}}}[f(\boldsymbol{x})] = \mathbf{G} \nonumber \end{align}\tag{6}\] Our results later reduce the certified robustness bounds to depend only on \(\|\mathbf{G}\|\). This reduction is crucial, as empirically estimating a full high-dimensional gradient would require prohibitively many samples.

Bounded Setting: The Role of the Mean. In many regression tasks, the output is naturally bounded (e.g., steering angles). We impose a pointwise bound constraint \(f(\boldsymbol{x}) \in [-M, M]\). As noted above, the specific value of the mean \(E = \mathbb{E}_{p_{\boldsymbol{z}}}[f(\boldsymbol{x})]\) determines where the feasible “box” \([-M, M]\) is located relative to the distribution’s center (see Appendix 11). Therefore, we must explicitly constrain the mean \(E\).

We consider two variations in this setting. As a baseline, we define the \((E, C) + M\) setting (detailed in Appendix 14), which constrains the mean, variance, and bounds. However, our primary focus is the tightest setting which utilizes all available statistics, denoted as \((E, C, G) + M\):

Mean, Variance, Gradient, and Boundedness. We formulate the full optimization problem over functions \(f(\boldsymbol{x})\) as follows: \[\label{eq:opt95bounded95full} \begin{align} & \underset{\|\boldsymbol{\delta}\|\le R, f}{\text{maximize}} & & \left| \mathbb{E}_{p_{\boldsymbol{z}+\boldsymbol{\delta}}}[f(\boldsymbol{x})] - \mathbb{E}_{p_{\boldsymbol{z}}}[f(\boldsymbol{x})] \right| \\ & \text{subject to} & & \mathbb{E}_{p_{\boldsymbol{z}}}[f(\boldsymbol{x})] = E \\ & & & \mathrm{Var}_{p_{\boldsymbol{z}}}(f(\boldsymbol{x})) \le C \\ & & & \nabla \mathbb{E}_{p_{\boldsymbol{z}}}[f(\boldsymbol{x})] = \mathbf{G}\\ & & & -M \le f(\boldsymbol{x}) \le M \end{align}\tag{7}\]

4.2 Methodology: Calculus of Variations↩︎

The optimization defined above requires finding a function \(f^*\) from an infinite-dimensional space. This is the domain of the Calculus of Variations, which generalizes standard optimization to find functions that extremize a functional (an integral depending on the function).

This connects directly to the certification techniques used in classification. The seminal work by [4] relies on the Neyman-Pearson Lemma to find the worst-case base classifier. The Neyman-Pearson Lemma is fundamentally a variational result: it solves for the optimal indicator function (a 0/1 test) that maximizes power subject to a constraint on size [13].

In regression, we generalize this principle. We are not restricted to boolean 0/1 functions; our base function \(f\) is real-valued. We employ the method of Lagrange multipliers within the Calculus of Variations framework (see Appendix 10 for the detailed variational setup). We construct a Lagrangian functional \(\mathcal{L}[f]\) that incorporates our observed constraints and solve the Euler-Lagrange equations \(\frac{\delta \mathcal{L}}{\delta f} = 0\) to derive the analytical form of the worst-case function \(f^*\).

5 Worst Case Function↩︎

5.1 Warm-up: Worst Case Function with Variance Constraint↩︎

We begin by deriving the certified radius under the simplest assumption: that the base function \(f\) has bounded variance under the smoothing distribution. This serves as a baseline for our method and illustrates the variational approach in an isotropic setting.

We note that the form of the worst case function identified in Proposition [prop:variance95only] is the direct analogue of the Neyman-Pearson lemma: the optimal function is determined by the likelihood ratio of the two distributions. And the closed form formula for \(R\) shows how \(R\) varies with \(\epsilon, C\) and \(\sigma\).

While Proposition [prop:variance95only] provides a valid certificate, it generates an isotropic bound, relying solely on the scalar variance \(C\). Since variance is a rotationally invariant statistic, this certificate cannot distinguish between directions of low and high sensitivity. In the next section, we tighten this bound by incorporating the gradient.

5.2 Incorporating Gradient Estimate↩︎

We remain in the unbounded setting. We now tighten the variance-only bound by incorporating the gradient \(\mathbf{G}= \nabla g(\boldsymbol{z})\), solving the problem defined in 6 .

To derive the certified radius, we decouple the nested maximization. We first characterize the worst-case base function \(f\) for a fixed perturbation \(\boldsymbol{\delta}\).

The constant \(k(\boldsymbol{\delta})\) is guaranteed to be real-valued by the variance lower bound \(C \ge \sigma^2\|\mathbf{G}\|_2^2\) established in Proposition 2 in Appendix 9.

Having characterized the worst-case function \(f^*\) for any fixed \(\boldsymbol{\delta}\), we now solve the outer optimization problem over the perturbation. The following proposition establishes that the worst-case perturbation is necessarily aligned with the gradient direction \(\mathbf{G}\), reducing the multi-dimensional search to a single dimension.

Combining these results yields the closed-form worst-case shift.

As with the previous proposition, the variance lower bound \(C \ge \sigma^2\|\mathbf{G}\|_2^2\) (Proposition 2) ensures this value is real. Notably, this bound depends only on the norm of the gradient \(\|\mathbf{G}\|_2\), not its direction.

We observe that \(\Delta_{\max}(R)\) is strictly monotonically increasing with respect to \(R\) for \(R \ge 0\). The first term is monotonic because the function \(e^x-1-x\) is strictly increasing for \(x>0\) and the second term \(R\|G\|_2\) is linear in \(R\). This monotonicity guarantees that the equation \(\Delta_{\max}(R) =\epsilon\) has a unique solution. We calculate the certified radius for this setting using the bisection search procedure, detailed in Algorithm 1. Notably, the bound explicitly depends on the norm \(\left\|\mathbf{G}\right\|_2\), confirming that incorporating gradient information tightens the certificate.

The optimal function \(f^*\) reveals a fundamental structure: it is a superposition of exponential and linear ramps determined by the projection of the input onto the gradient direction. This generalizes the discrete “slab" constructions of [9] to continuous regression.

5.3 Bounded Function Value Setting↩︎

We now address the bounded setting, where \(f(\boldsymbol{x}) \in [-M, M]\). As formulated in Section 4.1, we consider two cases: utilizing the Mean and Variance as constraints \(((E, C) + M)\), and full set of constraints: Mean, Variance, and Gradient \(((E, C, \mathbf{G}) + M)\).

5.3.0.1 Case 1: Mean and Variance Constraints \(((E, C) + M)\).

As derived in Appendix 14, we simplify the search for the worst-case \(d\)-dimensional function \(f(\boldsymbol{x})\) by showing it is equivalent to finding a worst-case univariate function \(\phi(t)\). Here, the scalar variable \(t \sim \mathcal{N}(0, \sigma^2)\) represents the projection of the Gaussian noise along the perturbation direction. Consequently, the optimal univariate function \(\phi^*(t)\) takes the form of a clipped affine transformation of the likelihood ratio. While structurally related to Proposition [prop:variance95only], the bounded solution is distinct in two ways: the pointwise bounds \([-M, M]\) induce clipping in the optimal function form, and the explicit mean constraint introduces an additional Lagrange multiplier. We numerically solve for the mean and variance multipliers to satisfy the integral constraints, while the clipping explicitly enforces the pointwise bounds (Algorithm 4 in Appendix 14.3).

5.3.0.2 Case 2: Mean, Variance, and Gradient Constraints \(((E, C, \mathbf{G}) + M)\).

We now consider the fully constrained setting, denoted as \((E, C, \mathbf{G}) + M\), which yields the tightest certification bounds by utilizing all available statistics. Directly solving the variational problem 7 over the space of multivariate functions \(f: \mathbb{R}^d \to \mathbb{R}\) is computationally intractable. A primary theoretical contribution of this work is proving that this high-dimensional optimization admits an exact reduction to a tractable optimization over univariate functions.

This reduction renders the variational problem tractable: instead of optimizing over multivariate functions, we optimize over a univariate function \(\phi(t)\). The problem thus depends only on the scalar shift magnitude \(\left\|\boldsymbol{\delta}\right\|_2\), with the pointwise bounds inducing a clipped solution structure.

Algorithmic Solution. While clipping precludes an analytical solution for the multipliers, the dual problem remains convex. We efficiently solve for \((\lambda, \mu, \nu)\) using a dual optimization procedure (Algorithm 6 in Appendix 15) to evaluate the exact worst-case shift \(\Delta(R)\), enabling the computation of the certified radius via bisection search (Algorithm 5).

6 Estimation of Variance and Gradient Norm↩︎

To compute the certified radius, we estimate the variance \(C\) and squared gradient norm \(\|\nabla g(z)\|_2^2\) using U-statistics on \(n\) independent Gaussian samples [14], [15]. U-statistics provide minimum-variance unbiased estimators for these quantities.

To ensure the joint validity of our certificate with probability at least \(1-\alpha\), we apply a union bound, allocating failure probability \(\alpha/2\) to the variance estimation and \(\alpha/2\) to the gradient norm estimation. We construct valid two-sided confidence intervals for both parameters at significance level \(\alpha/2\) based on their asymptotic normality. For certification, we utilize the conservative upper bound of the variance interval and optimize the worst-case gradient norm over its respective valid interval. Detailed definitions of the kernels and the derivation of these confidence intervals are provided in Appendix 16.

A key distinction between our framework and prior work in classification certification (e.g., [4]) is the nature of the statistical guarantees. While classification often relies on exact Clopper-Pearson bounds for Binomial counts, our regression setting requires estimating continuous moments. Deriving tight, exact finite-sample bounds for these quantities without making strong distributional assumptions is often intractable. Consequently, we rely on the asymptotic normality of U-statistics. This represents a limitation of our current work, as strict finite-sample validity is not guaranteed. However, our convergence analysis (Appendix 17.2) indicates that for practical sample sizes \((n \ge 10,000)\), the asymptotic approximation is highly accurate and the certificates remain empirically sound.

7 Experiments↩︎

7.1 Validating Certificate Tightness on Synthetic Data↩︎

We evaluate the tightness and soundness of our unbounded \((C, G)\) certifier on synthetic functions where the true worst-case radius is computable, comparing performance against \(\alpha\)-smoothing. We test on three synthetic functions: a smooth quadratic bowl, a one-sided slice that is flat on one side and increases linearly beyond a threshold, and a “sandwich” function with lower and upper plateaus separated by a narrow transition band. For each function and each output tolerance \(\epsilon \in \{0.2, 0.5\}\), we cross-validate over noise levels \(\sigma \in \{0.1, 0.2, 0.5\}\) to find the optimal \(\sigma\) that maximizes mean certified radius. For \(\alpha\)-smoothing, we cross-validate over both \(\sigma \in \{0.1, 0.2, 0.5\}\) and trimming rates \(\alpha \in \{0.35, 0.49\}\) to find the best combination. The optimal \(\alpha\) value is \(0.49\) for all cases. We then compare both methods at their respective optimal parameter settings. All experiments use \(N = 5,000\) samples, \(P = 0.9\) success probability for both methods (denoting the target acceptance probability for \(\alpha\)-smoothing vs. statistical confidence for ours), and certify on the same 10 randomly sampled test points per function. Detailed function definitions and ground truth computation methodology are provided in Appendix 17.1.

Table 1 shows the comparison at optimal parameter settings. Our method consistently achieves larger mean certified radii across all functions and tolerance levels.

2.5pt

7.2 MNIST Rotation Certification and Comparison↩︎

We evaluate our certification methods on the MNIST rotation prediction task, a practical high-dimensional regression problem where the goal is to predict the rotation angle of a rotated handwritten digit. We use an E(2)-equivariant neural network (E2CNN) [16] trained to predict rotation angles in the range \([-\pi, \pi]\) radians (\([-180^\circ, 180^\circ]\) degrees). The task is well-suited for certification evaluation because: (1) the output is naturally bounded (angles are periodic), (2) the input space is high-dimensional (784 pixels), and (3) small pixel perturbations can induce large changes in continuous prediction [17], [18], making robustness certification critical for deployment. We evaluate on 100 stratified test samples (10 per digit class) from the rotated MNIST test set, using multiple noise levels \(\sigma \in \{0.06, 0.12, 0.25, 0.50, 0.75\}\) and output tolerance \(\epsilon = 10^{\circ}\) (approximately \(0.175\) radians). All methods are evaluated at a success probability of \(0.9\) (denoting the target acceptance probability for \(\alpha\)-smoothing vs. statistical confidence for ours) to ensure a fair comparison . Additional experimental details are provided in Appendix 17.3.

Figure 2 (a) shows the cumulative distribution functions (CDFs) of certified radii for our two methods—\((E, C) + M\) and \((E, C, G) + M\)—and the \(\alpha\)-smoothing baseline, using the optimal noise level \(\sigma\) for each method (determined by maximizing mean certified radius). At their respective best \(\sigma\) values, our \((E, C, G) + M\) method achieves a mean certified radius of \(0.210\) pixels with a median of \(0.208\) pixels at \(\sigma = 0.75\), certifying all 100 samples (0% failure rate at \(\epsilon = 10^\circ\)). In contrast, the \((E, C) + M\) method achieves a mean radius of \(0.090\) pixels at its optimal \(\sigma = 0.06\), demonstrating that gradient constraints provide substantial improvement (2.33\(\times\) mean improvement, 99/100 samples). The \(\alpha\)-smoothing baseline achieves a mean radius of \(0.120\) pixels at \(\sigma = 0.06\) with \(P = 0.9\), but fails to certify 6 out of 100 samples (6% abstain rate, certifying zero radius). Our \((E, C, G) + M\) method outperforms \(\alpha\)-smoothing by 1.76\(\times\) in mean certified radius and achieves a larger certified radius on 94 out of 100 samples (including all 6 samples where \(\alpha\)-smoothing fails to certify).

Our method demonstrates robust performance across a wide range of noise levels, while \(\alpha\)-smoothing exhibits catastrophic failure at large \(\sigma\) values (see Appendix 17.3 for detailed per-sigma results). At \(\sigma \geq 0.5\), \(\alpha\)-smoothing certifies zero radius on 100% of test points (failing to certify any robustness within \(\epsilon = 10^\circ\)), while our method maintains robust performance with mean radii of \(0.156\) and \(0.210\) pixels at \(\sigma = 0.5\) and \(\sigma = 0.75\), respectively. This robustness stems from our method’s ability to leverage variance and gradient information simultaneously, whereas \(\alpha\)-smoothing relies solely on probability mass estimation which becomes unreliable when the probability mass \(p_A\) becomes too small at large noise levels. The complete failure of \(\alpha\)-smoothing at moderate-to-large noise levels (\(\sigma \geq 0.5\)) limits its practical utility, as it restricts the method to a narrow range of noise levels and makes it fragile to hyperparameter selection, whereas our method can leverage larger noise levels to achieve better certified radii.

Although MNIST rotation has a periodic output space, our certificate does not exploit this periodicity: it treats the predictor as a bounded scalar-valued function and uses only moment and gradient information. To verify that the improvement is not specific to periodic targets, we next evaluate an aperiodic scalar regression task.

7.3 Aperiodic Age Estimation↩︎

To check that the improvement is not specific to periodic angle prediction, we also certify a face age-estimation task with scalar, aperiodic outputs. We use the UTKFace dataset and evaluate on 100 fixed test images. The base predictor is a pretrained MiVOLO-v2 face-age regressor. We use the model in inference mode without retraining or fine-tuning. We resize each image to \(64 \times 64\) RGB, scale pixel values to \([0,1]\), and add Gaussian smoothing noise in this normalized pixel space. Each noisy image is then bilinearly resized to the model’s \(384 \times 384\) input resolution and normalized channel-wise using the pretrained model’s image normalization constants before inference. We use the known valid age range \([0,116]\) as the output bound for our \((E,C)+M\) and \((E,C,G)+M\) certificates, with \(M=116\).

We use the same fixed-grid comparison protocol as for MNIST rotation. For our methods, we sweep \(\sigma \in \{0.06,0.12,0.25,0.5,0.75\}\); for \(\alpha\)-smoothing, we additionally sweep trimming levels \(\alpha \in \{0.35,0.49\}\). The output tolerance is \(\epsilon=6\) years. All methods are evaluated at success probability \(0.9\), using \(5{,}000\) Monte Carlo samples for our statistical estimates and \(5{,}000\) samples for the \(\alpha\)-smoothing probability estimate. Appendix 17.4 gives the full details.

Figure 2 (b) shows the CDFs using the best fixed configuration for each method, selected by mean certified radius. Both \((E,C)+M\) and \((E,C,G)+M\) select \(\sigma=0.75\), achieving mean certified radii \(1.651\) and \(2.152\), respectively, while \(\alpha\)-smoothing selects \((\sigma,\alpha)=(0.75,0.49)\) and achieves mean radius \(0.817\). As in the MNIST rotation task, \(\alpha\)-smoothing often fails to certify a nonzero radius: it returns radius zero on 59/100 UTKFace test images at its best fixed configuration. Our \((E,C,G)+M\) certificate is larger than \(\alpha\)-smoothing on 76/100 points. On those points, the average improvement is \(1.830\); on the 24 points where \(\alpha\)-smoothing is larger, its average advantage is only \(0.234\), giving a total gain/loss ratio (the sum of our positive radius gains divided by the sum of \(\alpha\)-smoothing’s positive gains) of \(24.74\). The appendix reports the corresponding radius–accuracy tradeoff for the smoothed age predictor.

7.4 Certified Accuracy Analysis↩︎

To evaluate the practical utility of our certificates, we employ three complementary metrics that disentangle the quantity (coverage) from the quality (precision) of the certified predictions.

First, following standard protocols [4], we define absolute certified accuracy (ACR) at radius \(R\). This measures the overall utility of the system: what fraction of the test set is both robust and correct? \[\text{cert\_acc}_{\text{abs}}(R) = \frac{1}{n} \sum_{i=1}^{n} \mathbf{1}\{d(\hat{y}_i, y_i) \leq \epsilon_{\text{corr}} \wedge r_i \geq R\}, \label{eq:certified95accuracy95absolute}\tag{9}\] where \(n\) is the total number of samples, \(d(\cdot, \cdot)\) is the circular distance, and \(\epsilon_{\text{corr}} = 10^\circ\) is the correctness tolerance.

Second, to assess the reliability of valid certificates, we measure conditional certified accuracy. This asks: among the subset of samples \(\mathcal{S}_R = \{i : r_i \geq R\}\) that represent valid certificates, what fraction are actually correct? \[\text{cert\_acc}_{\text{cond}}(R) = \frac{1}{|\mathcal{S}_R|} \sum_{i \in \mathcal{S}_R} \mathbf{1}\{d(\hat{y}_i, y_i) \leq \epsilon_{\text{corr}}\}. \label{eq:certified95accuracy95conditional}\tag{10}\]

Third, to capture the continuous nature of regression beyond binary correctness, we measure the certified mean circular absolute error (MAE). This is the average circular distance among samples that meet the robustness requirement: \[\mathrm{cert\;circular\;MAE}(R) = \frac{1}{|S_R|}\sum_{i\in S_R} d(\hat{y}_i,y_i). \label{eq:certified95mean95distance}\tag{11}\]

Table 2 presents these metrics on the MNIST rotation task. To ensure a fair comparison across radii, for each method we select a single fixed noise level \(\sigma\) that maximizes the area under the certified accuracy curve (or minimizes mean distance error).

Results and Analysis. Our \((E, C, G) + M\) method demonstrates superior performance across all metrics. In terms of utility (Absolute Accuracy), it maintains \(95\%\) accuracy up to \(R=0.15\) pixels, whereas the baseline \((E, C) + M\) drops to \(2\%\) and \(\alpha\)-smoothing drops to \(39\%\). In terms of reliability (Conditional Accuracy), all methods achieve high scores (near \(100\%\)), indicating that when a certificate is issued, the prediction is reliably within the \(10^\circ\) tolerance. Finally, regarding precision (mean circular MAE), our method achieves a mean error of \(3.90^\circ\) at \(R=0.15\), significantly lower than the \(4.65^\circ\) of \(\alpha\)-smoothing and \(5.80^\circ\) of \((E, C) + M\). Notably, at the largest radius threshold \(R=0.25\), our method not only certifies \(5\%\) of the data (where others certify roughly \(0\%\)) but does so with exceptional precision, achieving a mean distance of just \(1.94^\circ\).

3pt

7.5 Certificate Validation: Tightness and Soundness↩︎

We validate our certificates through two complementary analyses: tightness (comparing certified radii to optimization-based radii) and soundness (verifying that certified radii do not exceed empirically-found worst-case radii).

7.5.0.1 Tightness Analysis.

We compare certified radii with optimization-based radii computed using projected gradient descent (PGD) [19] to find worst-case adversarial perturbations. We evaluate on 100 MNIST rotation samples at \(\sigma = 0.5\) using \((E, C, G) + M\). Of these, 18 samples hit the PGD search bound \(R_{\max} = 5.0\) pixels and are excluded from ratio analysis (these samples have optimization-based radius \(\geq 5.0\) but unknown exact value). For the remaining 82 samples, Table 3 and Figure 3 show a mean ratio of \(3.41\times\) (median \(3.48\times\)), indicating optimization-based radii are \(3.4\times\) larger than certified radii on average. Most samples have ratios between \(2\times\) and \(5\times\), demonstrating consistent conservatism. The \(3.4\times\) average gap is comparable to the discrepancy between certified and empirical radii observed in robust classification, where true robustness often exceeds certified bounds by factors of two or more [4]. Detailed methodology is in Appendix 17.6.

7.5.0.2 Soundness Validation.

Samples with ratio \(\geq 1.0\) are empirically sound (certified radius \(\leq\) optimization-based radius). Of 100 samples, 98 (98%) are sound: 80 of the 82 analyzed samples have ratio \(\geq 1.0\), and all 18 samples hitting the search bound are definitively sound (optimization-based \(\geq 5.0 >\) any certified radius). The 2 samples (2%) with ratio \(< 1.0\) (smallest: \(0.39\times\)) reflect finite-sample estimation error, expected with \(95\%\) confidence intervals (theoretical \(5\%\) failure rate vs. observed \(2\%\)). This confirms our method behaves as theoretically predicted. See Appendix 17.6 for detailed boundary analysis.

Impact Statement↩︎

This paper presents work whose goal is to advance the field of Machine Learning, specifically by enhancing the reliability and safety of regression models. Our framework for certifying robustness has potential positive societal impacts in safety-critical domains such as autonomous navigation and medical imaging, where guaranteeing that continuous predictions remain within safe margins is essential. We do not feel there are specific negative ethical or societal consequences that must be highlighted here.

The Appendix is organized as follows. Appendix 8 provides an extended discussion of related work, placing our contributions in the context of existing adversarial defense and certification literature. Appendix 9 establishes fundamental lemmas and propositions, including the application of Stein’s Lemma and variance lower bounds, which serve as building blocks for our proofs. Appendix 10 introduces the Calculus of Variations framework and details how our certification objectives are formulated as variational problems. Appendix 11 clarifies the role of the mean value constraint in bounded versus unbounded settings. We then provide the complete derivations for our certified radii: Appendix 12 covers the unbounded variance-only setting (Section 5.1), Appendix 13 details the unbounded gradient-constrained setting (Section 5.2), Appendix 14 derives the solution for the bounded setting with mean and variance constraints, and Appendix 15 presents the proofs for our fully constrained setting (Section 5.3), including the dimensionality reduction theorem. Finally, Appendix 16 details the statistical estimation procedures using U-statistics, and Appendix 17 provides comprehensive experimental details, including model architectures, hyperparameters, and additional validation results.

8 Extended Related Work↩︎

8.0.0.1 Historical context: From attacks to certification.

Early work on adversarial robustness focused on empirical defenses and attack strategies [20], [21]. As adaptive attacks systematically defeated many empirical defenses, the field pivoted toward provable robustness. This produced a spectrum of certification paradigms, including mixed-integer formulations and convex relaxations. However, randomized smoothing [4] eventually emerged as the standard for scalable \(\ell_2\) certification in modern deep networks, trading off the tightness of exact verification for the ability to scale to high-dimensional datasets like ImageNet.

8.0.0.2 Regularization and local sensitivity.

Complementary to certification, several works use input-gradient regularization to control local sensitivity during training. [22] derives per-example robustness bounds from the loss gradient norm and a modulus-of-continuity term. While not a smoothing method per se, its objective—shrinking local change by limiting first-order behavior—aligns with our use of gradient constraints. Our work differs by providing post-hoc certification for smoothed models rather than a training-time penalty, though the geometric intuition is similar.

8.0.0.3 Detailed comparison with existing regression methods.

Recent works by [6] and [5] formulate robustness as a probability of inclusion. Specifically, [6] defines a neighborhood around the ground-truth value and estimates a lower confidence bound on the probability that the noisy observation falls within this box. [5] improves upon this for unbounded outputs using \(\alpha\)-trimming. A theoretical limitation of these approaches is their reliance on a fixed anchor (such as the base prediction \(f(x)\) or ground truth \(y\)) to define the “safe” region. While it is possible to experimentally center these regions on the smoothed prediction, doing so violates the independence assumptions required for their guarantees, as the center becomes a stochastic variable. By contrast, our method analytically certifies the stability of the value \(g(x)\) itself, ensuring the certificate is theoretically valid for the actual deployed prediction.

8.0.0.4 Theoretical limitations of smoothing.

A parallel body of work maps out the inherent limitations of smoothing-based approaches. [23] and [24] show that for \(\ell_2\) certification, the certified radius must inevitably decrease as \(\sigma/\sqrt{d}\) in high dimensions to maintain accuracy. Similarly, [25] argue that for \(\ell_\infty\) perturbations, randomized smoothing may be fundamentally incapable of nontrivial certification for natural image distributions. Our work focuses on tightening the constant factor of these bounds via gradient information, but is subject to the same fundamental dimensional scaling laws inherent to the smoothing framework.

9 Useful Lemmas and Propositions↩︎

10 Our Problem as Calculus of Variations Problem↩︎

The most basic problem in the calculus of variations is to find a function, let’s call it \(y(\boldsymbol{x})\), that extremizes (maximizes or minimizes) a functional \(I[y]\) of the form: \[I[y] = \int L(\boldsymbol{x}, y, y') \,dx\] where \(L\) is some function of the independent variable \(x\), the function \(y(\boldsymbol{x})\), and its derivative \(y'(\boldsymbol{x})\).

When we have constraints, we use the method of Lagrange multipliers. For a problem of the form:

• Extremize: \(I[y] = \int L(\boldsymbol{x}, y, y') \,d\boldsymbol{x}\)

• Subject to: \(J[y] = \int K(\boldsymbol{x}, y, y') \,d\boldsymbol{x}= \text{Constant}\)

We form a new Lagrangian functional, \(\mathcal{L}[y] = I[y] - \lambda J[y]\), and find the function \(y(\boldsymbol{x})\) that makes it stationary by solving its Euler-Lagrange equation.

\[\frac{\partial L}{\partial f} - \frac{d}{d\boldsymbol{x}} \frac{\partial L}{\partial f'}=0\]

10.1 Mapping Our Problem to the Standard Form↩︎

In our specific problem, the unknown function we are trying to find is \(f(\boldsymbol{x})\). We now map our objective and constraints to the standard form described above.

10.2 The Functional to Maximize, \(I[f]\)↩︎

Our objective is to maximize the change in the smoothed function \(g(\boldsymbol{z})\), which is given by: \[I[f] = \int f(\boldsymbol{x}) p_{z+\delta}(\boldsymbol{x}) \,d\boldsymbol{x}- \int f(\boldsymbol{x}) p_z(\boldsymbol{x}) \,d\boldsymbol{x}\] We can combine this into a single integral to identify the integrand \(L\): \[I[f] = \int \underbrace{f(\boldsymbol{x}) \left( p_{\boldsymbol{z}+\delta}(\boldsymbol{x}) - p_{\boldsymbol{z}}(\boldsymbol{x}) \right)}_{L(\boldsymbol{x}, f(\boldsymbol{x}))} \,d\boldsymbol{x}\] So, the integrand of our main functional is \(L(\boldsymbol{x}, f(\boldsymbol{x})) = f(\boldsymbol{x}) (p_{\boldsymbol{z}+\boldsymbol{\delta}}(\boldsymbol{x}) - p_{\boldsymbol{z}}(\boldsymbol{x}))\). Notice that this function does not depend on the derivative, \(f'(\boldsymbol{x})\).

10.3 The Constraint Functionals, \(J_i[f]\)↩︎

We have two types of constraints that we can express as functionals.

• Variance Constraint: The constraint is \(\text{Var}_{p_{\boldsymbol{z}}}(f(\boldsymbol{x})) \le C\). We can write this functional as: \[J_1[f] = \int \underbrace{(f(\boldsymbol{x}) - g(\boldsymbol{z}))^2 p_{\boldsymbol{z}}(\boldsymbol{x})}_{K_1(\boldsymbol{x}, f(\boldsymbol{x}))} \,d\boldsymbol{x}\le C,\] where \(g(\boldsymbol{z}) = \int f(\boldsymbol{x}) p_{\boldsymbol{z}}(\boldsymbol{x}) d\boldsymbol{x}\) and depends on \(f\). The integrand for this constraint is \(K_1(\boldsymbol{x}, f(\boldsymbol{x})) = (f(\boldsymbol{x}) - g(\boldsymbol{z}))^2 p_{\boldsymbol{z}}(\boldsymbol{x})\). Again, this does not depend on \(f'(\boldsymbol{x})\).

• Gradient Constraint: This is a vector of \(d\) constraints, where the \(i\)-th component is: \[J_{2,i}[f] = \int \underbrace{f(\boldsymbol{x}) \frac{\boldsymbol{x}_i - \boldsymbol{z}_i}{\sigma^2} p_{\boldsymbol{z}}(\boldsymbol{x})}_{K_{2,i}(\boldsymbol{x}, f(\boldsymbol{x}))} \,d\boldsymbol{x}= G_i\] The integrand is \(K_{2,i}(\boldsymbol{x}, f(\boldsymbol{x})) = f(\boldsymbol{x}) \frac{\boldsymbol{x}_i - \boldsymbol{z}_i}{\sigma^2} p_{\boldsymbol{z}}(\boldsymbol{x})\). This also does not depend on \(f'(\boldsymbol{x})\).

10.4 The Full Lagrangian Functional↩︎

We combine the objective and constraints using our Lagrange multipliers, the scalar \(\lambda\) and the vector \(\boldsymbol{\mu}\): \[\mathcal{L}[f] = I[f] - \lambda J_1[f] - \boldsymbol{\mu}^T \mathbf{J_2}[f]\] Writing this with a single integrand , \(L_{\text{total}}\), we get: \[\mathcal{L}[f] = \int L_{\text{total}}(\boldsymbol{x}, f(\boldsymbol{x}) \, d\boldsymbol{x}= \int \left( L(\boldsymbol{x},f) - \lambda K_1(\boldsymbol{x},f) - \boldsymbol{\mu}^T \mathbf{K_2}(\boldsymbol{x},f) \right)\,d\boldsymbol{x}\]

10.5 The Simplified Euler-Lagrange Equation↩︎

The full Euler-Lagrange equation for the functional \(\mathcal{L}[f]\) is: \[\frac{\partial L_{\text{total}}}{\partial f} - \frac{d}{d \boldsymbol{x}}\frac{\partial L_{\text{total}}}{\partial f'} = 0\]

Since none of our integrands (\(L\), \(K_1\), or \(\mathbf{K_2}\)) depend on the derivative \(f'(\boldsymbol{x})\), the term \(\frac{\partial L_{\text{total}}}{\partial f'}\) is identically zero. This makes the entire second part of the Euler-Lagrange equation vanish.

Therefore, the necessary condition for finding the optimal function \(f^*(\boldsymbol{x})\) simplifies from a potentially complex differential equation to a simple algebraic equation: \[\frac{\partial L_{\text{total}}}{\partial f} = 0\] This final equation is precisely the “functional derivative" step to find the analytical form of the worst-case function \(f^*(\boldsymbol{x})\).

10.6 Details of Derivation of Stationary Condition from Lagrangian Functional↩︎

In this section, we show how to derive the stationary condition from the Lagrangian functional, given that \(g(\boldsymbol{z})\) is dependent on \(f\). We use the optimization problem with only variance constraint 5 to exemplify this. Note other optimization problems which have their stationary condition’s corresponding term follow similarly.

Recall the Lagrangian for optimization problem 5 is \[\mathcal{L}(f, \lambda) = \underbrace{\int f(\boldsymbol{x})(p_{\boldsymbol{z}+\boldsymbol{\delta}}(\boldsymbol{x}) - p_{\boldsymbol{z}}(\boldsymbol{x}))d\boldsymbol{x}}_{\text{Term A}} - \lambda \underbrace{\left( \int (f(\boldsymbol{x}) - g(\boldsymbol{z}))^2 p_{\boldsymbol{z}}(\boldsymbol{x})d\boldsymbol{x}- C \right)}_{\text{Term B}}\] where \(g(\boldsymbol{z}) = \int f(\boldsymbol{x})p_{\boldsymbol{z}}(\boldsymbol{x})d\boldsymbol{x}\), and \(p_{\boldsymbol{z}}\) is defined as in Section 4 to be \(\mathcal{N}(\boldsymbol{z}, \sigma^2 I)\). We focus on the variance term. The functional derivative of Term A is straightforward: \[\frac{\delta}{\delta f(\boldsymbol{x})} \left( \text{Term A} \right) = p_{\boldsymbol{z}+\boldsymbol{\delta}}(\boldsymbol{x}) - p_{\boldsymbol{z}}(\boldsymbol{x})\]

For Term B, we utilize the fact that for a distribution centered at \(\boldsymbol{z}\): \[\mathrm{Var}_{\boldsymbol{e}}(f)=\frac{1}{2}\int \int(f(\boldsymbol{z}+\boldsymbol{e}')-f(\boldsymbol{z}+\boldsymbol{e}))^2 dP(\boldsymbol{e})dP(\boldsymbol{e}')\] where \(\boldsymbol{e}\) and \(\boldsymbol{e}'\) are independent copies of the same random variable, that is \(\boldsymbol{e}, \boldsymbol{e}' \sim \mathcal{N}(0, \sigma^2 I)\).

We calculate the functional derivative by perturbing \(f\) to \(f + \epsilon h\):

\[\begin{align} \frac{d}{d \epsilon}\mathrm{Var}[f + \epsilon h] \bigg|_{\epsilon=0} &= \int \int ( f (\boldsymbol{z}+\boldsymbol{e}') - f(\boldsymbol{z}+\boldsymbol{e})) (h(\boldsymbol{z}+\boldsymbol{e}') - h(\boldsymbol{z}+\boldsymbol{e})) dP(\boldsymbol{e})dP(\boldsymbol{e}') \\ &= 2 \int f (\boldsymbol{z}+\boldsymbol{e}') h(\boldsymbol{z}+\boldsymbol{e}')dP(\boldsymbol{e}') - 2 \int\int f(\boldsymbol{z}+\boldsymbol{e}) h(\boldsymbol{z}+\boldsymbol{e}') dP(\boldsymbol{e})dP(\boldsymbol{e}')\\ &= 2 \int \left[ f (\boldsymbol{z}+\boldsymbol{e}') - \underbrace{\int f(\boldsymbol{z}+\boldsymbol{e}) dP(\boldsymbol{e})}_{g(\boldsymbol{z})} \right] h(\boldsymbol{z}+\boldsymbol{e}')dP(\boldsymbol{e}') \end{align}\]

To recover the density form used in the Lagrangian, we perform a change of variables for the outer integral. Let \(\boldsymbol{x}= \boldsymbol{z}+ \boldsymbol{e}'\). Then \(dP(\boldsymbol{e}') = p_{\boldsymbol{z}}(\boldsymbol{x})d\boldsymbol{x}\). This yields: \[\int 2 (f(\boldsymbol{x}) - g(\boldsymbol{z})) h(\boldsymbol{x}) p_{\boldsymbol{z}}(\boldsymbol{x}) d\boldsymbol{x}\]

Thus, variance term’s functional derivative with respect to \(f(\boldsymbol{x})\) is \(2(f(\boldsymbol{x}) - g(\boldsymbol{z}))p_{\boldsymbol{z}}(\boldsymbol{x})\).

Putting functional derivative from Term A and Term B together we get the stationary condition shown in Equation 12 : \[(p_{\boldsymbol{z}+\boldsymbol{\delta}}(\boldsymbol{x}) - p_{\boldsymbol{x}}(\boldsymbol{x})) - 2\lambda (f(\boldsymbol{x}) - g(\boldsymbol{z})) p_{\boldsymbol{x}}(\boldsymbol{x}).\]

11 Mean Value Constraint↩︎

We compare the role of the mean function \(g(\boldsymbol{z})=\mathbb{E}_{p_{\boldsymbol{z}}}[f]\) in the unbounded and bounded settings. In both cases, the objective function (the harm) \(\Delta(\boldsymbol{\delta})=\mathbb{E}_{p_{\boldsymbol{z}+\boldsymbol{\delta}}}[f]-\mathbb{E}_{p_{\boldsymbol{z}}}[f]\) is invariant under adding constants to \(f\). However, the impact of the mean on the feasible set differs significantly.

In the unbounded setting, the constraints are shift-invariant; thus, the specific value of \(g(\boldsymbol{z})\) is immaterial. In the bounded setting, the pointwise box constraint couples the admissible centered function to the mean \(g(\boldsymbol{z})\) via a shift of the feasible interval. Consequently, an unknown \(g(\boldsymbol{z})\) forces us to consider the worst-case shift, leading to a looser (more conservative) certificate.

12 Proofs in Section 5.1↩︎

In this section, we restate and give proof for Proposition [prop:variance95only], which gives form for worst case function and closed form for certified radius in unbounded function value setting with variance constraint.

13 Details and Proofs of Section 5.2↩︎

In Section 5.2, we derive the worst-case certificate for the setting constrained by the variance \(C\) and the gradient vector \(\mathbf{G}\)  6 .

13.1 Lagrangian Formulation and Optimal Function Form↩︎

Since the problem involves maximizing a linear objective over the intersection of a strictly convex set (the variance ball) and a hyperplane (the gradient constraint), the resulting feasible set remains strictly convex. Consequently, the optimal solution shape is unique (up to a constant shift), and there is only one set of Lagrange multipliers that satisfies the stationarity and constraint conditions.

To solve for the worst-case function \(f^*\), we employ the Calculus of Variations. We introduce a scalar Lagrange multiplier \(\lambda\) for the variance constraint and a vector multiplier \(\boldsymbol{\mu}\in \mathbb{R}^d\) for the gradient constraint. Using the identity of gradient of \(g(\boldsymbol{z})\) (equivalently \(\mathbb{E}_{p_{\boldsymbol{z}}}[f]\)) derived in Lemma 2, \(\nabla_{\boldsymbol{z}} \mathbb{E}_{p_{\boldsymbol{z}}}[f] = \mathbb{E}_{p_{\boldsymbol{z}}}[f(\boldsymbol{x}) \frac{\boldsymbol{x}-\boldsymbol{z}}{\sigma^2}]\), we construct the full Lagrangian functional:

\[\begin{align} \mathcal{L}[f, \lambda, \boldsymbol{\mu}] = &\int f(\boldsymbol{x}) (p_{\boldsymbol{z}+\boldsymbol{\delta}}(\boldsymbol{x}) - p_{\boldsymbol{z}}(\boldsymbol{x})) d\boldsymbol{x}\\ &- \lambda \left( \int (f(\boldsymbol{x}) - g(\boldsymbol{z}))^2 p_{\boldsymbol{z}}(\boldsymbol{x}) d\boldsymbol{x}- C \right) \\ &- \boldsymbol{\mu}^\top \left( \int f(\boldsymbol{x}) \frac{\boldsymbol{x}-\boldsymbol{z}}{\sigma^2} p_{\boldsymbol{z}}(\boldsymbol{x}) d\boldsymbol{x}- \mathbf{G} \right) \end{align}.\]

We seek the function \(f^*\) that makes this functional stationary. Taking the functional derivative with respect to \(f(\boldsymbol{x})\) (as detailed in Appendix Section 10.6 regarding the variance term) and setting it to zero yields the Euler-Lagrange condition:

\[\frac{\delta \mathcal{L}}{\delta f(\boldsymbol{x})} = (p_{\boldsymbol{z}+\boldsymbol{\delta}}(\boldsymbol{x}) - p_{\boldsymbol{z}}(\boldsymbol{x})) - 2\lambda (f(\boldsymbol{x}) - g(\boldsymbol{z})) p_{\boldsymbol{z}}(\boldsymbol{x}) - \boldsymbol{\mu}^\top \frac{\boldsymbol{x}-\boldsymbol{z}}{\sigma^2} p_{\boldsymbol{z}}(\boldsymbol{x}) = 0\]

Since the Gaussian density \(p_{\boldsymbol{z}}(\boldsymbol{x})\) is strictly positive everywhere, we can divide through by \(p_{\boldsymbol{z}}(\boldsymbol{x})\) to isolate the optimal function form. This gives:

\[\left( \frac{p_{\boldsymbol{z}+\boldsymbol{\delta}}(\boldsymbol{x})}{p_{\boldsymbol{z}}(\boldsymbol{x})} - 1 \right) - 2\lambda (f(\boldsymbol{x}) - g(\boldsymbol{z})) - \boldsymbol{\mu}^\top \frac{\boldsymbol{x}-\boldsymbol{z}}{\sigma^2} = 0\]

To simplify the notation, we define two auxiliary terms: the likelihood ratio shift \(A(\boldsymbol{x}, \boldsymbol{\delta})\) and the function \(B(\boldsymbol{x})\):

\[\begin{align} \label{eq:likelihood95A} A(\boldsymbol{x}, \boldsymbol{\delta}) \triangleq \frac{p_{\boldsymbol{z}+\boldsymbol{\delta}}(\boldsymbol{x})}{p_{\boldsymbol{z}}(\boldsymbol{x})} - 1 = \exp\left( \frac{\boldsymbol{\delta}^\top(\boldsymbol{x}-\boldsymbol{z})}{\sigma^2} - \frac{\|\boldsymbol{\delta}\|^2}{2\sigma^2} \right) - 1 \end{align}\tag{13}\]

\[\begin{align} \label{eq:score95func95B} B(\boldsymbol{x}) \triangleq \frac{\boldsymbol{x}-\boldsymbol{z}}{\sigma^2} \end{align}\tag{14}\]

Substituting these definitions into the stationarity condition, we obtain:

\[A(\boldsymbol{x}, \boldsymbol{\delta}) - 2\lambda (f(\boldsymbol{x}) - g(\boldsymbol{z})) - \boldsymbol{\mu}^\top B(\boldsymbol{x}) = 0\]

We observe that the variance constraint must be active, i.e., \(\lambda > 0\). If we assume \(\lambda = 0\), the stationarity condition implies \(A(\boldsymbol{x}, \boldsymbol{\delta}) = \boldsymbol{\mu}^\top B(\boldsymbol{x})\). However, \(A(\boldsymbol{x}, \boldsymbol{\delta})\) involves the exponential term \(e^{\boldsymbol{\delta}^{\top}(\boldsymbol{x}-\boldsymbol{z})/\sigma^2}\), whereas \(\boldsymbol{\mu}^\top B(\boldsymbol{x})\) is linear in \(\boldsymbol{x}\). A non-constant exponential function cannot equal a linear function almost everywhere (they intersect at most on a set of measure zero). Thus, we must have \(\lambda \neq 0\). Since \(\lambda\) corresponds to an inequality constraint, dual feasibility requires \(\lambda \ge 0\), and therefore \(\lambda > 0\).

Solving for the centered function \(f^*(\boldsymbol{x}) - g(\boldsymbol{z})\):

\[f^*(\boldsymbol{x}) - g(\boldsymbol{z}) = \frac{1}{2\lambda} \left( A(\boldsymbol{x}, \boldsymbol{\delta}) - \boldsymbol{\mu}^\top B(\boldsymbol{x}) \right)\]

Again defining the scalar constant \(k \triangleq \frac{1}{2\lambda}\), we arrive at the general form of the worst-case function:

\[\begin{align} \label{eq:general95form95unbounded95c95g95constrained} f^*(\boldsymbol{x}) - g(\boldsymbol{z}) = k \left( A(\boldsymbol{x}, \boldsymbol{\delta}) - \boldsymbol{\mu}^\top B(\boldsymbol{x}) \right) \end{align}\tag{15}\]

This derivation reveals that the optimal function is a linear combination of the likelihood ratio (exponential in \(\boldsymbol{x}\)) and a function linear in \(\boldsymbol{x}\). The unknown constants \(k\) and \(\boldsymbol{\mu}\) must now be determined by satisfying the primal constraints.

13.2 Derivation of the Worst-Case Function↩︎

Having established the general structure of the optimal function in Section 13.1, we must now enforce the primal constraints to determine the Lagrange multipliers. By solving for the scalar \(k\) and vector \(\boldsymbol{\mu}\) in terms of the problem parameters (\(\boldsymbol{\delta}, C, \mathbf{G}\)) and substituting them back into the general form, we obtain the explicit worst-case function characterized in the following proposition.

To prove this, we first establish three key Gaussian expectations involving the likelihood ratio shift \(A(\boldsymbol{x}, \boldsymbol{\delta})\) and the function \(B(\boldsymbol{x})\) defined previously.

13.3 Worst-Case Harm and Direction↩︎

With the optimal function \(f^*\) determined for a fixed perturbation \(\boldsymbol{\delta}\), we can now compute the worst-case harm and optimize the perturbation direction.

14 Bounded Function Value with Mean and Variance Constraint↩︎

Let \(p_{\boldsymbol{z}} = \mathcal{N}(\boldsymbol{z}, \sigma^2 I_d)\). We consider the following optimization problem where the function \(f\) is constrained by a fixed mean \(E\), a variance budget \(C\), and pointwise bounds \([-M, M]\): \[\begin{align} \label{eq:primal95ec95m} \max_{f, \boldsymbol{\delta}} \quad & \Delta(f, \boldsymbol{\delta}) = \mathbb{E}_{p_{\boldsymbol{z}+\boldsymbol{\delta}}}[f(\boldsymbol{x})] - \mathbb{E}_{p_{\boldsymbol{z}}}[f(\boldsymbol{x})] \\ \text{s.t.} \quad & \left\|\boldsymbol{\delta}\right\|_2 \le R \notag \\ (E) \quad & \mathbb{E}_{p_{\boldsymbol{z}}}[f(\boldsymbol{x})] = E \notag \\ (C) \quad & \mathbb{E}_{p_{\boldsymbol{z}}}[(f(\boldsymbol{x})-E)^2] \le C \notag \\ (M) \quad & -M \le f(\boldsymbol{x}) \le M \quad \forall \boldsymbol{x}\in \mathbb{R}^d \notag \end{align}\tag{21}\]

14.0.0.1 High-Variance vs. Low-Variance Regimes

Before deriving the worst-case function, we distinguish between two fundamental regimes determined by the relationship between the variance budget \(C\) and the geometric bounds of the function space.

A function bounded by \([-M, M]\) with mean \(E\) achieves its maximum possible variance when it takes only the extreme values \(\{-M, M\}\) (a “bang-bang” solution). The variance of such a function is exactly \(M^2 - E^2\). This bound delineates two cases:

1. The High-Variance Regime (\(C \ge M^2 - E^2\)): The variance constraint is redundant. The adversary can construct a step function (Neyman-Pearson solution) using only the values \(\{-M, M\}\) that satisfies the variance constraint trivially. This reduces to the standard bounded certification problem (e.g., [4], [9]).

2. The Low-Variance Regime (\(C < M^2 - E^2\)): The variance constraint is strictly tighter than the natural bound imposed by \(M\). The adversary is precluded from using a simple step function and must construct a solution that transitions smoothly or is clipped to satisfy the stricter variance budget.

In the context of robust regression, we are specifically interested in the Low-Variance Regime, as it characterizes functions that possess limited volatility. Consequently, we proceed under the assumption that \(C < M^2 - E^2\).

14.1 Derivation of the 1-Dimensional Reduction↩︎

We claim that the \(d\)-dimensional optimization problem can be reduced to a 1-dimensional optimization problem over a scalar function \(\phi(t)\) by exploiting the symmetry of the Gaussian distribution and the isotropic nature of the constraints.

14.1.0.1 Symmetry and Projection

Since the prior distribution \(p_{\boldsymbol{z}}\) is isotropic and the constraints \((E)\), \((C)\), and \((M)\) are rotationally invariant, the problem is symmetric with respect to the direction of perturbation. Without loss of generality, we fix the perturbation direction. Let \(\boldsymbol{\delta}\) be an arbitrary perturbation with magnitude \(\alpha = \left\|\boldsymbol{\delta}\right\|_2\). We define the unit vector \(\boldsymbol{u}= \boldsymbol{\delta}/ \alpha\).

We decompose the input space relative to \(\boldsymbol{u}\). Let \(t = \boldsymbol{u}^\top (\boldsymbol{x}- \boldsymbol{z})\) be the projection of the input \(\boldsymbol{x}\) onto the direction of perturbation. We define a 1-dimensional function \(\phi: \mathbb{R}\to \mathbb{R}\) as the conditional expectation of \(f(\boldsymbol{x})\) given the projection \(t\): \[\phi(t) = \mathbb{E}_{p_{\boldsymbol{z}}}[f(\boldsymbol{x}) \mid \boldsymbol{u}^\top (\boldsymbol{x}- \boldsymbol{z}) = t]\]

14.1.0.2 Equivalence of the Objective Function

The likelihood ratio between the shifted distribution \(p_{\boldsymbol{z}+\boldsymbol{\delta}}\) and the base distribution \(p_{\boldsymbol{z}}\) depends only on the projection \(t\): \[w(t) \triangleq \frac{p_{\boldsymbol{z}+\boldsymbol{\delta}}(\boldsymbol{x})}{p_{\boldsymbol{z}}(\boldsymbol{x})} - 1 = \exp\left( \frac{\alpha t}{\sigma^2} - \frac{\alpha^2}{2\sigma^2} \right) - 1\] Using the tower property of conditional expectation, the objective function becomes: \[\begin{align} \mathbb{E}_{p_{\boldsymbol{z}+\boldsymbol{\delta}}}[f(\boldsymbol{x})] - \mathbb{E}_{p_{\boldsymbol{z}}}[f(\boldsymbol{x})] &= \mathbb{E}_{p_{\boldsymbol{z}}}\left[ f(\boldsymbol{x}) \left( \frac{p_{\boldsymbol{z}+\boldsymbol{\delta}}(\boldsymbol{x})}{p_{\boldsymbol{z}}(\boldsymbol{x})} - 1 \right) \right] \notag \\ &= \mathbb{E}_{p_{\boldsymbol{z}}}\left[ f(\boldsymbol{x}) w(t) \right] \notag \\ &= \mathbb{E}_{T \sim \mathcal{N}(0, \sigma^2)} \left[ \mathbb{E}[f(\boldsymbol{x}) \mid T=t] w(t) \right] \notag \\ &= \mathbb{E}_{T \sim \mathcal{N}(0, \sigma^2)} [\phi(T) w(T)] \end{align}\]

14.1.0.3 Preservation of Constraints

We verify that the projected function \(\phi(t)\) satisfies the constraints if \(f(\boldsymbol{x})\) does.

• Mean (E): \(\mathbb{E}_{p_{\boldsymbol{z}}}[f(\boldsymbol{x})] = \mathbb{E}_T[\mathbb{E}[f(\boldsymbol{x})|T]] = \mathbb{E}_T[\phi(T)] = E\).

• Boundedness (M): Since \(-M \le f(\boldsymbol{x}) \le M\) holds almost surely, the conditional expectation satisfies \(-M \le \phi(t) \le M\) almost surely.

• Variance (C): Applying Jensen’s Inequality to the convex function \(h(y) = (y-E)^2\): \[(\phi(t) - E)^2 = (\mathbb{E}[f(\boldsymbol{x}) - E \mid t])^2 \le \mathbb{E}[(f(\boldsymbol{x}) - E)^2 \mid t]\] Taking the expectation over \(T\) yields \(\mathbb{E}_T[(\phi(T) - E)^2] \le C\).

Thus, maximizing the original objective is equivalent to finding the worst-case 1-D function \(\phi(t)\) and shift \(\alpha \in [0, R]\).

14.2 The Lagrangian Formulation and Derivation of Optimal Form↩︎

We solve the inner optimization problem for the function \(\phi\) given a fixed shift \(\alpha\). Once the optimal \(\phi_\alpha^*\) is found, we maximize over \(\alpha \in [0, R]\). The primal optimization problem is: \[\begin{align} \max_{\phi} \quad & \int \phi(t) w(t) p_0(t) \, dt \\ \text{s.t.} \quad & \int (\phi(t) - E)^2 p_0(t) \, dt \le C \quad (\text{Variance})\\ & \int \phi(t) p_0(t) \, dt = E \quad (\text{Mean})\\ & -M \le \phi(t) \le M \quad \forall t \in \mathbb{R} \quad (\text{Boundedness}) \end{align}\]

We construct the Lagrangian \(\mathcal{L}\) with Lagrange multipliers \(\lambda\) (variance), \(\mu\) (mean), \(\eta_1(t)\) (upper bound), and \(\eta_2(t)\) (lower bound). \[\begin{align} \mathcal{L}(\phi, \lambda, \mu, \eta_1, \eta_2) &= \int \phi(t) w(t) p_0(t) \, dt \\ &\quad - \lambda \left( \int (\phi(t) - E)^2 p_0(t) \, dt - C \right) \\ &\quad - \mu \left( \int \phi(t) p_0(t) \, dt - E \right) \\ &\quad - \int \eta_1(t) (\phi(t) - M) p_0(t) \, dt \\ &\quad - \int \eta_2(t) (-\phi(t) - M) p_0(t) \, dt \end{align}\] Grouping terms, the Lagrangian becomes: \[\mathcal{L} = \int \left[ \phi(t) w(t) - \lambda (\phi(t)-E)^2 - \mu \phi(t) - \eta_1(t)(\phi(t) - M) + \eta_2(t)(\phi(t) + M) \right] p_0(t) \, dt + \text{const}\]

Optimality Conditions↩︎

1. Stationarity Condition
Taking the functional derivative of \(\mathcal{L}\) with respect to \(\phi(t)\) and setting it to zero gives the necessary condition for the optimal solution \(\phi^*\): \[\frac{\delta \mathcal{L}}{\delta \phi(t)} \bigg|_{\phi^*} = \left( w(t) - 2\lambda (\phi^*(t) - E) - \mu - \eta_1(t) + \eta_2(t) \right) p_0(t) = 0\] Rearranging for \(\phi^*(t)\): \[\begin{align} \label{eq:stationarity95mean95ec95m} 2\lambda (\phi^*(t) - E) = w(t) - \mu - \eta_1(t) + \eta_2(t) \end{align}\tag{22}\]

2. KKT Conditions
For \(\phi^*(t)\) to be optimal, it must satisfy:

• Primal Feasibility: \(-M \le \phi^*(t) \le M\).

• Dual Feasibility: \(\lambda \ge 0\), \(\eta_1(t) \ge 0\), and \(\eta_2(t) \ge 0\).

• Complementary Slackness: \[\lambda(\mathbb{E}[(\phi^*-E)^2]-C) = 0, \quad \eta_1(t)(\phi^*(t) - M) = 0, \quad \eta_2(t)(-\phi^*(t) - M) = 0\]

Proof that the Variance Constraint is Active (\(\lambda > 0\))↩︎

We prove that \(\lambda > 0\) by contradiction, given the assumption that \(C < M^2 - E^2\). Assume the variance constraint is inactive (\(\lambda = 0\)). The stationarity condition (Eq. 22 ) becomes \(\eta_1(t) - \eta_2(t) = w(t) - \mu\). Based on the sign of \(w(t) - \mu\), one of the multipliers must be positive, forcing \(\phi^*(t)\) to the boundary:

• If \(w(t) > \mu \implies \eta_1(t) > 0 \implies \phi^*(t) = M\).

• If \(w(t) < \mu \implies \eta_2(t) > 0 \implies \phi^*(t) = -M\).

Thus, if \(\lambda=0\), the optimal function \(\phi^*(t)\) takes only the extreme values \(\{-M, M\}\) almost everywhere. The variance of such a solution is maximal: \(\mathrm{Var}(\phi^*) = M^2 - E^2\). This contradicts the Low-Variance Regime assumption that \(\mathrm{Var}(\phi) \le C < M^2 - E^2\). Therefore, the assumption \(\lambda=0\) is false, and we must have \(\lambda > 0\).

Deriving the Optimal Form↩︎

With \(\lambda > 0\), we define the unconstrained candidate function \(h(t)\): \[h(t; \lambda, \mu) \triangleq E + \frac{w(t) - \mu}{2\lambda}\] Analyzing the KKT conditions for the active/inactive bound cases yields the clipped solution: \[\begin{align} \label{eq:optimal95phi95final95ec95m} \phi^*(t) = \text{clip}_{[-M, M]} \left( E + \frac{w(t) - \mu}{2\lambda} \right) \end{align}\tag{23}\]

14.3 Numerical Procedure↩︎

14.3.0.1 Solving for the Lagrange Multipliers

The multipliers \(\lambda \text{ and } \mu\) are determined by enforcing the variance and mean constraints. Let \(\mathbf{v} = [\lambda, \mu]^\top\) be the vector of unknowns. We define the residual vector \(\mathbf{F}(\lambda, \mu)\): \[\begin{align} F_1(\lambda, \mu) &\triangleq \mathbb{E}_T \left[ \left( \phi^*(T; \lambda, \mu) - E \right)^2 \right] - C = 0 \tag{24}\\ F_2(\lambda, \mu) &\triangleq \mathbb{E}_T \left[ \phi^*(T; \lambda, \mu) \right] - E = 0 \tag{25} \end{align}\]

14.3.0.2 Monotonicity and Bisection

Let \(W(R) \triangleq \max_{\|\boldsymbol{\delta}\|_2 \le R} \Delta(\boldsymbol{\delta})\). Due to rotational invariance, this is equivalent to maximizing over the shift magnitude \(\alpha \in [0, R]\). Since the set of allowable shifts grows with \(R\), \(W(R)\) is non-decreasing. This validates the use of bisection search for the certified radius.

14.3.0.3 Numerical Root Finding

The system \(\mathbf{F}(\lambda, \mu) = \mathbf{0}\) corresponds to the stationarity conditions of the dual problem. Since the primal problem is convex with linear constraints, the dual function is concave. This ensures that the solution \([\lambda^*, \mu^*]\) is unique and corresponds to the global optimum. We solve this system using a standard numerical root-finding algorithm (e.g., hybrid Powell method or Levenberg-Marquardt).

15 Bounded Function with Mean, Variance and Gradient Constraint↩︎

In this section, we show in the bounded function setting, with mean, variance and gradient constraint, the nested optimization of finding the worst case perturbation in \(d\)-dimension \(\|\boldsymbol{\delta}\|\le R\) while finding the worst case function \(f\) can be reduced to finding a \(1\)-dimensional perturbation \(\alpha \in \mathbb{R}\), and a function \(\phi\) which takes one dimensional input. We first show that the worst case perturbation \(\boldsymbol{\delta}\) has to align with gradient \(\mathbf{G}\)’s direction. This result allows us to reduce the problem on \(d\)-dimensional variable to be on \(1\)-dimension variable.

15.1 Proofs for Section 5.3: Bounded Function Value↩︎

Theorem 3 shows in the bounded function value setting with variance, gradient, and expected value constraints, the worst perturbation \(\boldsymbol{\delta}\) aligns with gradient \(\mathbf{G}\), we can reduce the optimization problem in \(d\)-dimensional \(\boldsymbol{\delta}\) to an optimization problem on \(1\)-dimensional problem input.

We build on this result. In below proof given to Proposition 4, we write out the optimization problem of finding the worst-case perturbation in \(1\)-dimensional variable \(\alpha \in \mathbb{R}\) and a function \(\phi\) of a 1-dimensional input.

15.2 General Properties of the Optimal Solution↩︎

Before deriving the explicit functional form of the solution, we establish a critical property regarding the variance constraint. Specifically, we show that for the class of functions relevant to neural networks (continuous or smooth functions), the variance constraint is strictly active (\(\lambda > 0\)).

15.3 The Lagrangian Formulation and Optimal Form↩︎

15.4 Numerical Procedure↩︎

Solving for the Lagrange Multipliers

With the analytic form of \(\phi^*(t)\) established in Corollary [cor:bounded95ecgm95optimal95form], we must determine the values of the Lagrange multipliers \(\lambda, \mu, \text{ and } \nu\) that satisfy the original constraints. We treat the multipliers as variables in a dual optimization problem.

Monotonicity of the Bounded Problem Before detailing the solver, we verify that the worst-case harm is monotonic in the radius \(R\), which validates the use of bisection search.

Let \(W(R)\) denote the worst-case function value change given a perturbation budget \(R\): \[W(R) \triangleq \max_{\|\boldsymbol{\delta}\|_2 \le R} \Delta(\boldsymbol{\delta}).\] Using our 1-D reduction, this is equivalent to: \[W(R)=\max _{\alpha \in \{-R, R\}} \Delta(\alpha).\] Observe that the set of feasible perturbations grows with \(R\). If \(R_1 < R_2\), then any perturbation valid for \(R_1\) is also valid for \(R_2\) (conceptually, though our reduction fixes \(\|\boldsymbol{\delta}\|=R\), the relaxation to inequality \(\|\boldsymbol{\delta}\| \le R\) yields the same optimum due to convexity/concavity arguments in the radius). Maximizing over a larger set yields a value greater than or equal to maximizing over a subset. Thus, \(W(R_2) \ge W(R_1)\). This guarantee ensures that a unique solution (or a single transition point) exists for \(W(R)=\epsilon\), making bisection search valid.

15.4.0.1 Lagrange Dual Method

The dual function \(g(\lambda, \mu, \nu)\) is derived by maximizing the Lagrangian with respect to \(\phi\) for fixed multipliers: \[g(\lambda, \mu, \nu) = \max_{\phi: -M-E \le \phi \le M-E} \mathcal{L}(\phi, \lambda, \mu, \nu)\] Recall the Lagrangian from Appendix 15.3 (excluding pointwise multipliers which are handled implicitly by the maximization): \[\begin{align} \mathcal{L} &= \int \phi(t) w(t) p_0(t) \, dt - \lambda \left( \int \phi(t)^2 p_0(t) \, dt - C \right) \\ &\quad - \mu \left( \int \phi(t) t p_0(t) \, dt - \sigma^2 \|\mathbf{G}\|_2 \right) - \nu \left( \int \phi(t) p_0(t) \, dt - 0 \right) \end{align}\] Collecting constant terms and grouping the integrals, we express the dual function as: \[g(\lambda, \mu, \nu) = \underbrace{\int \left( \phi^*(t)[w(t) - \mu t - \nu] - \lambda \phi^*(t)^2 \right) p_0(t) \, dt}_{\text{Optimized Integral Part}} + \underbrace{\lambda C + \mu \sigma^2 \|\mathbf{G}\|_2}_{\text{Constant Parts}}\] where \(\phi^*(t)\) is the optimal clipped form dependent on \((\lambda, \mu, \nu)\). We minimize this dual function to find the optimal multipliers: \[\min_{\lambda > 0, \mu \in \mathbb{R}, \nu \in \mathbb{R}} g(\lambda, \mu, \nu)\]

15.4.0.2 Gradients of the Dual

By the Envelope Theorem, the gradient of the dual function with respect to the multipliers is exactly the negative of the constraint residuals. Since the dual function is convex, we can use these gradients for optimization (e.g., via Gradient Descent): \[\begin{align} \frac{\partial g}{\partial \lambda} &= C - \int \phi^*(t)^2 p_0(t) \, dt \\ \frac{\partial g}{\partial \mu} &= \sigma^2 \|\mathbf{G}\|_2 - \int \phi^*(t) t p_0(t) \, dt \\ \frac{\partial g}{\partial \nu} &= 0 - \int \phi^*(t) p_0(t) \, dt \end{align}\]

16 Details of Estimation↩︎

In this section, we provide the theoretical details for the estimators used in Section 6. We first establish the validity of our joint certification procedure using the union bound, then state the general asymptotic properties of U-statistics, and finally apply these results to derive the specific confidence intervals for the variance and the gradient norm.

16.1 Joint Confidence Validity via Union Bound↩︎

To guarantee that the certified radius is valid with probability at least \(1-\alpha\), we must bound the estimation error for both the variance \(C\) and the gradient norm \(\|\mathbf{G}\|_2\). We use a union bound to split the failure budget, allocating \(\alpha/2\) to each term. We construct two-sided confidence intervals for both parameters at the significance level \(\alpha/2\) (i.e., using the critical value \(z_{1-\alpha/4}\)). This approach ensures validity as follows:

1. Variance (\(C\)): We calculate a two-sided interval \([C_{\text{low}}, C_{\text{high}}]\) and utilize the conservative upper bound \(C_{\text{high}}\). Due to the two-sided construction, the probability that the true variance exceeds this bound is at most \(\alpha/4\) (half of the \(\alpha/2\) budget).

2. Gradient Norm (\(\|\mathbf{G}\|_2\)): We calculate a two-sided interval \([G_{\text{low}}, G_{\text{high}}]\). Since the certified radius is non-monotonic with respect to the gradient norm (Corollary [cor:certified95radius]), we compute the radius by finding the worst-case gradient norm \(G^* \in [G_{\text{low}}, G_{\text{high}}]\) that minimizes the radius. The true gradient norm falls within this interval with probability \(1 - \alpha/2\).

The union bound guarantees that the joint certificate holds with probability at least \(1 - (\alpha/4 + \alpha/2) = 1 - 0.75\alpha \ge 1-\alpha\).

16.2 General Theory: Asymptotic Normality of U-Statistics↩︎

Both our variance and gradient norm estimators rely on U-statistics of order \(r=2\). We utilize the standard asymptotic theory for non-degenerate U-statistics [14], [15].

To construct confidence intervals, we estimate the asymptotic variance \(4\zeta_1\) using consistent sample estimates of the underlying moments.

16.3 Variance Estimation↩︎

We estimate the local variance \(C = \mathrm{Var}(f(\boldsymbol{x}))\) where \(\boldsymbol{x}\sim \mathcal{N}(\boldsymbol{z}, \sigma^2 I)\). Let \(Y = f(\boldsymbol{x})\) denote the random function output. We are given \(n\) i.i.d. samples \(Y_1, \dots, Y_n\). Using the identity \(\mathrm{Var}(Y) = \frac{1}{2}\mathbb{E}[(Y_1 - Y_2)^2]\), we employ the kernel \(h(y_i, y_j) = \frac{1}{2}(y_i - y_j)^2\).

Estimator equivalence. As noted in the main text, the U-statistic with this kernel simplifies exactly to the standard unbiased sample variance: \[\hat{\sigma}^2 = \binom{n}{2}^{-1} \sum_{i<j} \frac{1}{2}(Y_i - Y_j)^2 = \frac{1}{n-1}\sum_{i=1}^n (Y_i - \bar{Y})^2 = S^2.\]

Asymptotic Variance. We calculate \(\zeta_1\) by projecting the kernel onto one variable. To distinguish from the smoothing parameter \(\sigma\) and Lagrange multipliers used elsewhere, we explicitly denote the true moments of the output distribution \(Y\) as \(\mu_Y = \mathbb{E}[Y]\) and \(\sigma_Y^2 = \mathrm{Var}(Y)\). The projection is: \[h_1(t) = \mathbb{E}[h(t, Y_j)] = \frac{1}{2}\mathbb{E}[(t - Y_j)^2] = \frac{1}{2}((t - \mu_Y)^2 + \sigma_Y^2).\] The variance of this projection is: \[\zeta_1 = \mathrm{Var}(h_1(Y_i)) = \mathrm{Var}\left(\frac{1}{2}(Y_i - \mu_Y)^2\right) = \frac{1}{4}(m_4 - \sigma_Y^4),\] where \(m_4 = \mathbb{E}[(Y-\mu_Y)^4]\) is the fourth central moment of \(Y\).

Confidence Interval. We require a high-probability upper bound on the variance. We estimate the asymptotic variance using the sample fourth central moment \(\hat{m}_4\) and sample variance \(S^2\). Consistent with our union bound strategy, we use the critical value \(z_{1-\alpha/4}\) (corresponding to a two-sided interval with significance level \(\alpha/2\)). The upper bound is: \[C_{upper} = S^2 + z_{1-\alpha/4} \sqrt{\frac{\hat{m}_4 - (S^2)^4}{n}}.\]

16.4 Gradient Norm Estimation↩︎

We estimate the squared gradient norm \(\theta = \|\nabla g(\boldsymbol{z})\|^2\). Recall that \(\nabla g(\boldsymbol{z}) = \mathbb{E}[\mathbf{W}]\), where \(\mathbf{W} = \frac{1}{\sigma^2}(\boldsymbol{x}- \boldsymbol{z})f(\boldsymbol{x})\). We are given \(n\) i.i.d. samples \(\mathbf{W}_1, \dots, \mathbf{W}_n\). We use the kernel \(h(\mathbf{w}_i, \mathbf{w}_j) = \mathbf{w}_i^\top \mathbf{w}_j\).

Asymptotic Variance. The parameter to be estimated is \(\theta = \|\mathbb{E}[\mathbf{W}]\|^2\). We denote the true mean vector (the gradient) as \(\mathbf{m} = \mathbb{E}[\mathbf{W}]\) and the covariance matrix as \(\Sigma_{\mathbf{W}} = \mathrm{Cov}(\mathbf{W})\). The projection of the kernel is: \[h_1(\mathbf{w}) = \mathbb{E}[\mathbf{w}^\top \mathbf{W}_j] = \mathbf{w}^\top \mathbb{E}[\mathbf{W}_j] = \mathbf{w}^\top \mathbf{m}.\] The variance of this projection is \(\zeta_1 = \mathrm{Var}(\mathbf{W}_i^\top \mathbf{m}) = \mathbf{m}^\top \Sigma_{\mathbf{W}} \mathbf{m}\). Therefore, the U-statistic \(\hat{\theta}\) follows: \[\sqrt{n}(\hat{\theta} - \theta) \xrightarrow{d} \mathcal{N}(0, 4\mathbf{m}^\top \Sigma_{\mathbf{W}} \mathbf{m}).\]

Computational Stability in High Dimensions. Directly estimating the \(d \times d\) covariance matrix \(\Sigma_{\mathbf{W}}\) to compute \(\mathbf{m}^\top \Sigma_{\mathbf{W}} \mathbf{m}\) can be computationally expensive and unstable when \(d\) is large. However, we observe that \(\mathbf{m}^\top \Sigma_{\mathbf{W}} \mathbf{m}\) is simply the variance of the scalar projection \(\mathbf{W}^\top \mathbf{m}\). We estimate this scalar variance directly using the sample variance of the projected points \(\{ \mathbf{W}_i^\top \hat{\mathbf{m}} \}_{i=1}^n\), where \(\hat{\mathbf{m}} = \frac{1}{n}\sum \mathbf{W}_i\). This avoids explicit covariance matrix computation. Let \(\hat{v}^2\) denote this estimated scalar variance. We estimate the asymptotic variance scalar \(\mathbf{m}^\top \Sigma_{\mathbf{W}} \mathbf{m}\) using the sample variance of the projected points \(\mathbf{W}_i^\top \hat{\mathbf{m}}\). While this estimator may have finite-sample bias due to the dependency between \(\hat{\mathbf{m}}\) and \(\mathbf{W}_i\), it is a consistent estimator of the true asymptotic variance. By Slutsky’s theorem, consistency is sufficient to ensure the validity of the derived confidence intervals. While the estimation of the projection direction \(\mathbf{m}\) itself is subject to high-dimensional noise (scaling with \(d/n\)), the use of a conservative upper confidence bound on the scalar projection variance ensures that the resulting gradient norm estimate remains reliable for certification purposes.

Confidence Interval. We construct a two-sided confidence interval for \(\theta\) with significance level \(\alpha/2\) using the critical value \(z_{1-\alpha/4}\). Let \(\hat{\sigma}_{\theta} = \sqrt{4\hat{v}^2/n}\). The bounds for the squared norm \(\theta\) are: \[\theta_{low} = \hat{\theta} - z_{1-\alpha/4} \hat{\sigma}_{\theta}, \quad \theta_{high} = \hat{\theta} + z_{1-\alpha/4} \hat{\sigma}_{\theta}.\] To obtain the interval for the gradient norm \(G = \sqrt{\theta}\), we take the square root of the bounds (truncating the lower bound at zero): \[G_{low} = \sqrt{\max(0, \theta_{low})}, \quad G_{high} = \sqrt{\max(0, \theta_{high})}.\] These bounds define the search interval \([G_{low}, G_{high}]\) used in our certification procedure.

17 More Experiment Details↩︎

17.1 Synthetic Functions and Ground Truth Computation↩︎

We validate our certification methods on synthetic functions where the true worst-case radius can be computed analytically or via high-precision numerical methods. This allows us to verify soundness (certified radius \(\leq\) true radius) and evaluate tightness (how close certified radius is to true radius). All functions operate on 2-dimensional input space (\(d = 2\)), with input coordinates \(x = (x_1, x_2) \in \mathbb{R}^2\).

We evaluate three functions, each designed to test different aspects of the certification method.

The quadratic function is defined as \(f(x) = s \cdot \|x - c\|_2^2\) with center \(c = (0.0, 0.0)\) and scale \(s = 1.0\). This function is smooth and grows quadratically, making it ideal for testing tightness on well-behaved functions. This function has no output clipping, allowing us to test the method’s behavior when outputs can grow arbitrarily large.

The slice function is defined as \(f(x) = \max(0, x_1 - t)\) with threshold \(t = 0.0\). This function has a discontinuous jump at \(x_1 = t\) and grows linearly for \(x_1 > t\), testing the method’s robustness to non-smooth functions without explicit output bounds. The function is unbounded above but bounded below at zero.

The sandwich function is defined as \(f(x) = \max(0, \min(1, x_1 - t))\) with threshold \(t = 0.0\). This function is piecewise constant with a transition region, creating a challenging case for certification. While the function itself is bounded (outputs in \([0, 1]\)), we still use the unbounded \((C, G)\) certifier (without providing an explicit bound \(M\)) to test the method’s behavior when no output bound is supplied.

We compute the true worst-case radius \(R_{\text{true}}\) using optimization-based methods. Since these functions lack analytical solutions for the worst-case radius, we use a combination of Monte Carlo expectation estimation and differential evolution optimization [26] (implemented via scipy.optimize.differential_evolution). Differential evolution is a population-based global optimization algorithm that is well-suited for finding worst-case perturbations in non-convex optimization problems. For each test point \(z\), we first compute the baseline expectation \(g(\boldsymbol{z}) = E[f(\boldsymbol{z}+ \boldsymbol{e})]\), \(\boldsymbol{e}\sim \mathcal{N}(0,\sigma^2 I)\) using Monte Carlo estimation with \(50,000\) samples. We then use differential evolution to find the worst-case perturbation direction \(\boldsymbol{\delta}^*\) at a given radius \(r\) that maximizes the change in expectation \(|g(\boldsymbol{z}+ \boldsymbol{\delta}) - g(\boldsymbol{z})|\) subject to \(\|\boldsymbol{\delta}\|_2 \leq r\). The true radius is found via binary search: we search for the largest radius \(r\) such that the worst-case expectation change at radius \(r\) (computed using \(100,000\) samples for final evaluation) does not exceed the output tolerance \(\epsilon\). The binary search uses 30 iterations with tolerance \(10^{-4}\), and the differential evolution uses 50 iterations with population size 10. This approach provides high-precision estimates of the true worst-case radius while remaining computationally tractable.

All synthetic function experiments use consistent parameters: test points drawn uniformly from \([-1.0, 1.0] \times [-1.0, 1.0]\) using random seed 42, 10 randomly sampled test points per function, cross-validation over \(\sigma \in \{0.1, 0.2, 0.5\}\) for our certifier and over both \(\sigma \in \{0.1, 0.2, 0.5\}\) and \(\alpha \in \{0.35, 0.49\}\) for \(\alpha\)-smoothing to find optimal parameters, output tolerance levels \(\epsilon \in \{0.2, 0.5\}\), \(N = 5,000\) samples, and \(P = 0.9\) success probability for both our certifier and \(\alpha\)-smoothing.

17.2 Convergence Analysis of Radius Estimators↩︎

We validate that our statistical estimators converge correctly and that certified radii converge to their theoretical values as the sample size \(N\) increases. We perform this analysis on a single test sample from the MNIST rotation task, using sample sizes \(N \in \{100, 500, 1000, 5000, 10000\}\) with multiple independent trials per \(N\) value. Theoretical values for \(C\) and \(\|\mathbf{G}\|\) are estimated using large-sample Monte Carlo with \(N = 50,000\) samples.

17.2.1 Part 1: Estimator Convergence↩︎

For each \(N\), we compute point estimates \(\hat{C}\) and \(\hat{\theta}\) using U-statistic estimators, along with upper confidence bounds (UCBs) \(C_{\text{UCB}}\) and \(\theta_{\text{UCB}}\) using analytical confidence intervals based on asymptotic normality (z-critical or t-distribution). Figure 7 shows that both estimators converge to their true values with \(O(1/\sqrt{N})\) rate, with mean bias less than \(5\%\) at \(N = 10,000\) and confidence intervals shrinking proportionally.

17.2.2 Part 2: Radius Convergence↩︎

We compute the empirical certified radius \(R_{\text{empirical}}\) using estimated \(C_{\text{UCB}}\) and \(\theta_{\text{UCB}}\), and the theoretical radius \(R_{\text{theoretical}}\) using the large-sample estimates. Figure 8 shows that the empirical radius converges to the theoretical radius as \(N\) increases: at \(N = 5,000\), within \(10\%\) of the theoretical value, and at \(N = 10,000\), within \(5\%\). This validates that estimation error decreases with \(N\) and that the certification method is well-calibrated.

17.2.3 Experimental Setup↩︎

We use a single MNIST test image (index 0) with parameters \(\sigma = 0.5\), \(\epsilon = 10^\circ\) (0.175 radians), \(M = \pi\) radians, and confidence level \(95\%\) (failure probability \(\alpha = 0.05\)).

17.3 MNIST Rotation Experimental Details↩︎

17.3.1 Model Architecture↩︎

We use an E(2)-equivariant convolutional neural network (E2CNN) architecture for rotation angle prediction. The model uses rotation-equivariant convolutions that preserve the E(2) symmetry group structure with \(N=8\) rotational symmetries. The architecture consists of three convolutional layers with increasing channel dimensions (16, 32, 64), followed by a final layer mapping to invariant features (128 channels). The output layer is a regression head that predicts rotation angles in the range \([-\pi, \pi]\) radians via a \((\cos \theta, \sin \theta)\) representation. The model contains approximately 50,000 parameters. Full architecture details and implementation are available in our code repository.

17.3.2 Training Procedure↩︎

The model was trained on the rotated MNIST dataset using the Adam optimizer with an initial learning rate of \(0.001\), which was reduced by a factor of 0.5 every 20 epochs using a StepLR scheduler. Training was performed for 30 epochs with a batch size of 64, using mean squared error (MSE) loss on the \((\cos \theta, \sin \theta)\) representation of angles. The validation metric was circular mean absolute error (MAE) in degrees, and no additional data augmentation was applied since the dataset already contains full rotation range \([0^\circ, 360^\circ]\). During training, angles are converted to the range \([-180^\circ, 180^\circ]\) for consistency with the model’s output range. The data was split into 80% training and 20% validation sets.

Training was performed on GPU and converged within 30 epochs. The final model achieves a test set mean absolute error (MAE) of \(4.53^{\circ}\) (\(\pm 9.48^{\circ}\) standard deviation). Training and validation losses decreased steadily, with validation MAE stabilizing around epoch 20-25 (coinciding with the learning rate reduction at epoch 20). The model demonstrates good generalization, with test set performance matching validation performance.

17.3.3 Test Set Selection↩︎

We evaluate on 100 stratified test samples from the rotated MNIST test set, with 10 samples per digit class (0–9). Samples were selected using stratified random sampling with seed 42 for reproducibility. The test indices are stored in the test manifest file and match the indices used in all comparison experiments to ensure fair evaluation across methods.

17.3.4 Certification Parameters↩︎

For all certification experiments, we use a sample size of \(N = 10,000\) for statistical estimation, which was validated through convergence analysis, showing stable estimates with small confidence intervals. We employ a total confidence level of \(1-\alpha\) (total failure probability \(\alpha\)), where different experiments may use different values of \(\alpha\) depending on the desired success probability. The mean estimate \(E = \mathbb{E}_{\boldsymbol{e}\sim \mathcal{N}(0, \sigma^2 I)}[f(\boldsymbol{z}+\boldsymbol{e})]\) is computed as a point estimate (sample mean) without a confidence interval. Numerical integration is performed using 40 Gauss-Hermite quadrature points, and the dual optimization problem is solved with a tolerance of \(10^{-6}\) and a maximum of 1000 iterations. The output tolerance is set to \(\epsilon = 10^{\circ}\) (0.1745 radians).

Confidence level details: We use a total confidence level of \(1-\alpha\) (total failure probability \(\alpha\)). The allocation of \(\alpha\) depends on the method: For methods using variance and gradient constraints \((E, C, G) + M\), we allocate \(\alpha/2\) failure probability to each constraint (variance \(C\) and gradient norm \(\|\mathbf{G}\|\)) via the union bound, ensuring the total failure probability does not exceed \(\alpha\). For methods using only the variance constraint \((E, C) + M\), the entire failure probability \(\alpha\) is allocated to the variance constraint \(C\) (no split). For each constraint, we construct a two-sided confidence interval, which further reduces the per-tail failure probability. This ensures that with probability at least \(1-\alpha\), the relevant constraints are simultaneously satisfied within their respective confidence intervals. Note the per-sigma comparison table (Table 4) uses \(\alpha = 0.05\) (\(95\%\) confidence) for \((E, C, G) + M\) and \(\alpha = 0.025\) (\(97.5\%\) confidence) for \((E, C) + M\), while the main text comparison (Section 7.2) uses \(\alpha = 0.10\) (\(90\%\) confidence) for both methods. The mean estimate \(E\) is computed as the sample mean point estimate without a confidence interval. This is justified because with \(N=10,000\) samples, the standard error of the mean is \(\sqrt{C/N}\), which is small relative to the variance and gradient bounds that dominate the certified radius.

Quadrature points: We use 40 Gauss-Hermite quadrature points to numerically approximate Gaussian expectations of the form \(\mathbb{E}[f(T)]\) where \(T \sim \mathcal{N}(0, \sigma^2)\). Gauss-Hermite quadrature is a standard numerical integration technique that provides exact integration for polynomials up to degree \(2n-1\) using \(n\) quadrature points. We use this method because the dual optimization problem requires computing expectations of functions \(\phi^*(t)\) that may have discontinuities (due to clipping at \(\pm M\)), making analytical integration intractable. The choice of 40 points provides a good balance between accuracy and computational efficiency, as verified through sensitivity analysis.

Dual optimization tolerance: The dual optimization problem is solved using iterative methods (gradient-based optimization or root finding) to find the Lagrange multipliers \((\lambda, \mu, \nu)\) that satisfy the Karush-Kuhn-Tucker (KKT) conditions. The tolerance \(10^{-6}\) controls when the optimization algorithm terminates, ensuring that constraint violations (e.g., \(|\mathbb{E}[\phi^2] - C|\)) are below this threshold. This tolerance is chosen to be sufficiently tight to guarantee numerical accuracy while remaining computationally tractable.

17.3.5 Alpha-Smoothing Parameters↩︎

For the \(\alpha\)-smoothing baseline comparison, we evaluate both trimming parameters \(\alpha \in \{0.35, 0.49\}\) across all noise levels \(\sigma \in \{0.06, 0.12, 0.25, 0.50, 0.75\}\) and report results using the \(\alpha\) value that yields the best performance at each \(\sigma\). In practice, we found that \(\alpha = 0.49\) performs best for \(\sigma = 0.06\) (the optimal \(\sigma\) for \(\alpha\)-smoothing), while \(\alpha = 0.35\) may be preferred at other \(\sigma\) values. We use a probability threshold of \(P = 0.9\) and allocate \(n_{\text{tr}} = 10,000\) samples for \(p_A\) estimation (matching the \(N=10,000\) used in our methods). For the Binomial mapping (described below), we use \(n_{\text{sample}} = 500\) when \(\alpha=0.49\) or \(n_{\text{sample}} = 200\) when \(\alpha=0.35\). We enable circular distance for angle predictions, and all experiments use random seed 42 for reproducibility.

Note on sample sizes: The parameter \(n_{\text{tr}} = 10,000\) is the main Monte Carlo sample size used for estimating \(p_A\) (the probability that the output remains within tolerance), which directly corresponds to the \(N = 10,000\) samples used in our \((E, C) + M\) and \((E, C, G) + M\) methods. The parameter \(n_{\text{sample}}\) is used for a different purpose: computing the quantile \(q\) via Binomial distribution mapping, and does not affect the main statistical estimation. Therefore, the comparison is fair: all methods use 10,000 samples for the primary statistical estimation.

Binomial mapping: The \(\alpha\)-smoothing method uses a Binomial distribution mapping to convert the trimming parameter \(\alpha\) and sample size \(n_{\text{sample}}\) into a quantile threshold \(q\) that appears in the certified radius formula. Specifically, given \(\alpha\), \(n_{\text{sample}}\), and target probability \(P\), the method solves for \(q\) such that \(\text{BinomCDF}(\lfloor \alpha \cdot n_{\text{sample}} \rfloor, n_{\text{sample}}, 1-q) = P\), where \(\text{BinomCDF}\) is the cumulative distribution function of a Binomial random variable. This mapping accounts for the fact that \(\alpha\)-trimming removes the \(\alpha\) fraction of extreme values, and the quantile \(q\) represents the probability threshold that ensures the trimmed mean remains within the acceptable region with probability at least \(P\). For details, see [5].

Center method and comparison with our methods: The implementation of \(\alpha\)-smoothing method measures the output shift relative to the base function’s clean prediction \(f(z)\) rather than the smoothed regressor’s output \(g_\alpha(z)\). This choice is required by the theoretical guarantee of [5]. The acceptable region is defined as \([f(z) - \epsilon, f(z) + \epsilon]\) (with circular wrapping for angles). In contrast, our \((E, C) + M\) and \((E, C, G) + M\) methods certify output shift relative to the smoothed regressor’s expectation \(g(z) = \mathbb{E}_{\boldsymbol{e}}[f(\boldsymbol{z}+\boldsymbol{e})]\).

17.3.6 Noise Level Selection↩︎

We evaluate all methods across multiple noise levels \(\sigma \in \{0.06, 0.12, 0.25, 0.50, 0.75\}\) to understand performance sensitivity. Based on mean certified radius, the optimal \(\sigma\) values are: \(\sigma = 0.06\) for \((E, C) + M\) (mean radius = \(0.090\) pixels), \(\sigma = 0.75\) for \((E, C, G) + M\) (mean radius = \(0.207\) pixels), and \(\sigma = 0.06\) with \(\alpha = 0.49\) for \(\alpha\)-smoothing (mean radius = \(0.120\) pixels). The different optimal \(\sigma\) values reflect the different mechanisms: methods without gradient information benefit from smaller noise (less restrictive variance bounds), while methods with gradient information can leverage larger noise levels effectively.

17.3.7 Per-Sigma Performance↩︎

Table 4 shows the mean certified radius for each method across all noise levels. Our \((E, C, G) + M\) method demonstrates robust performance across all noise levels, maintaining non-zero certified radii at all \(\sigma\) values. In contrast, \(\alpha\)-smoothing exhibits catastrophic failure at large \(\sigma\) values (\(\sigma \geq 0.5\)), certifying zero radius for all samples. This highlights a critical limitation of methods that rely solely on probability mass estimation, which becomes unreliable when the probability mass becomes too small at large noise levels.

17.3.8 Computational Details↩︎

All certification experiments were performed on NVIDIA A100 GPUs using batched model inference for efficient Monte Carlo sampling. For each test sample, the certification process involves: (1) statistical estimation using \(N = 10,000\) Monte Carlo samples to compute variance \(C\) and gradient norm \(\|\mathbf{G}\|\) estimates, and (2) dual optimization to solve for the certified radius. The average wall-clock time per certificate computation is roughly 5–10 seconds on an A100 GPU, with the majority of time spent on statistical estimation (Monte Carlo sampling and U-statistic computation). The dual optimization step is relatively fast (typically under 1 second) due to the efficient numerical integration using Gauss-Hermite quadrature. For the full experimental evaluation across 100 test samples and 5 noise levels (\(\sigma \in \{0.06, 0.12, 0.25, 0.50, 0.75\}\)), the total compute time corresponds to roughly 1–2 GPU-hours. Parallelization across multiple test samples on the cluster further reduced the experimental turnaround time.

17.4 Age-Estimation Experimental Details↩︎

17.4.1 Task and Dataset↩︎

The age-estimation experiment evaluates robustness for an aperiodic scalar regression task. Each input is a face image and the target is the subject’s age in years. Unlike MNIST rotation, the output is not circular, so prediction error is measured by ordinary absolute error in years. We use the UTKFace dataset,¹ whose filenames encode age labels, and evaluate on 100 fixed test images. The same selected images are used for every method and every parameter setting.

17.4.2 Model and Preprocessing↩︎

The base regressor is a pretrained MiVOLO-v2 face-age model [27] used only for inference; ². We do not retrain or fine-tune it. The certification input is a \(64 \times 64\) RGB image with pixel values scaled to \([0,1]\). Gaussian smoothing noise is added at this resolution, so certified radii are \(\ell_2\) perturbation norms in the same normalized pixel space.

Before passing a noisy image to MiVOLO-v2, we resize it to the model’s \(384 \times 384\) input size using bilinear interpolation: each new pixel is computed as a weighted average of the four nearest pixels in the \(64 \times 64\) image. We then apply the pretrained model’s channel normalization constants and run inference. This resizing and normalization step is deterministic preprocessing applied after the randomized smoothing perturbation. We use the known valid age range \([0,116]\) as the output bound for our bounded certificates, with \(M=116\).

17.4.3 Test Set Selection↩︎

We evaluate on 100 fixed UTKFace test images selected from the held-out split and store the selected test indices for reuse across all sweeps. This is important because the UTKFace runs were split across separate jobs for our methods and for \(\alpha\)-smoothing; using the same selected indices ensures that all paired comparisons are pointwise comparisons on identical images.

17.4.4 Certification Parameters↩︎

For our \((E,C,G)+M\) certificate, we use \(N=5{,}000\) Monte Carlo samples per point to estimate the mean, variance, and gradient-norm quantities. The output tolerance is \(\epsilon=6\) years. We evaluate success probability \(0.9\) for all methods. For \((E,C,G)+M\), the failure probability is split across the estimated constraints used by the certificate. For \((E,C)+M\), we reuse the same saved Monte Carlo estimates but recompute the variance and mean confidence bounds with the two-way union-bound split for \((E,C)\) at success probability \(0.9\), avoiding another expensive model-inference sweep. The dual certificate optimization uses Gauss-Hermite quadrature as in the MNIST rotation experiment.

For each \(\sigma\), we also estimate the smoothed predictor \(g(x)=\mathbb{E}_{\eta}[f(x+\eta)]\) by Monte Carlo averaging over the same noisy model evaluations used for certification. This gives a real-valued non-adversarial accuracy measure for the smoothed model. The base regressor has clean mean absolute error \(6.73\) years on the 100 evaluated points; the appendix tradeoff plot reports how the smoothed MAE changes with \(\sigma\).

17.4.5 Single-Point Convergence Analysis↩︎

Mirroring the MNIST convergence analysis, we validate the statistical estimators on one representative UTKFace test image. The selected image has true age \(26\) years and clean model prediction \(26.55\) years. We use \(\sigma=0.06\) and output tolerance \(\epsilon=5\) years, and compare finite-sample estimates at \(N \in \{100,500,1000,2000,5000\}\) against large-sample pseudo-true values computed with \(N=10{,}000\). Each finite-sample setting is repeated over 10 independent trials.

The large-sample pseudo-true estimates are \(g(x)=28.24\), \(C=0.971\), \(\theta=\|G\|_2^2=768.73\), and \(\|G\|_2=27.73\). Figure 9 shows convergence of the variance estimator and squared-gradient estimator toward these reference values: point estimates are averaged over trials and shaded bands show the corresponding \(95\%\) confidence bands. The variance estimate is already within about \(2\%\) of the reference value across the tested sample sizes. The gradient estimator is noisier at small \(N\), as expected in the high-dimensional image input space, but moves toward the pseudo-true value as \(N\) increases.

Figure 9 reports confidence-interval coverage and certified-radius convergence. The nominal confidence level is \(95\%\) in this diagnostic run. Empirical coverage is at least \(90\%\) for \(C\) and \(100\%\) for \(\theta\) and \(\|G\|_2\) across the tested sample sizes. The certified radius converges quickly to the large-sample reference radius \(0.1088\): the mean radius is within \(1.8\%\) at \(N=500\), within \(0.8\%\) at \(N=2000\), and within \(0.5\%\) at \(N=5000\).

17.4.6 Alpha-Smoothing Parameters↩︎

For the \(\alpha\)-smoothing baseline, we evaluate trimming parameters \(\alpha \in \{0.35,0.49\}\) across the same noise levels \(\sigma \in \{0.06,0.12,0.25,0.5,0.75\}\). We use probability threshold \(P=0.9\) and \(n_{\mathrm{tr}}=5{,}000\) samples to estimate \(p_A\), matching the main Monte Carlo sample size used by our methods. As in the MNIST rotation experiment, the additional Binomial mapping sample size controls the quantile threshold used by \(\alpha\)-smoothing and is separate from the main \(p_A\) estimation sample size; we use \(n_{\text{sample}}=200\) for \(\alpha=0.35\) and \(n_{\text{sample}}=500\) for \(\alpha=0.49\).

17.4.7 Noise Level Selection and Per-Sigma Performance↩︎

Each method selects the fixed grid configuration with the largest mean certified radius over the 100 test images. Table 5 reports the per-sigma mean certified radius. For \(\alpha\)-smoothing, the table reports the better trimming parameter at each \(\sigma\); \(\alpha=0.49\) is best at all five noise levels in this sweep. All three methods select \(\sigma=0.75\) by mean certified radius.

4pt

4pt

17.4.8 Radius–Accuracy Tradeoff and Computation↩︎

The best-radius configuration for our methods uses the largest smoothing noise in the grid. As expected, larger \(\sigma\) improves certified radius but can degrade the smoothed predictor’s accuracy. Figure 10 reports this tradeoff using the Monte Carlo estimate of the smoothed age predictor’s MAE.

The experiment was run on NVIDIA A100 GPUs using batched model inference. We used mixed precision and tensor-native preprocessing for efficient Monte Carlo evaluation. The grid was sharded across test points and parameter settings on the cluster; the \((E,C)+M\) curve was then postprocessed from saved estimates without rerunning model inference.

17.5 Certified Accuracy Experimental Details↩︎

The metric definitions for absolute certified accuracy, conditional certified accuracy, and certified mean distance are provided in Section 7.4. Here we describe the experimental implementation details.

For rotation angle predictions, which are in the range \([-\pi, \pi]\) radians (\([-180^\circ, 180^\circ]\) degrees), we use circular distance to account for the periodic nature of angles. The circular distance between two angles \(\hat{y}\) and \(y\) (in degrees) is defined as \(d(\hat{y}, y) = \min(|\hat{y} - y|, 360^\circ - |\hat{y} - y|)\), which ensures that angles near the boundaries of \([-180^\circ, 180^\circ]\) are correctly compared. For example, if the true angle is \(-175^\circ\) and the prediction is \(175^\circ\), the circular distance is \(10^\circ\) (not \(350^\circ\)), correctly reflecting that these angles are close on the circle. Thus, the mean circular MAE reported in Table 2 averages this absolute circular error over certified points.

17.5.1 Sigma Selection and Evaluation Protocol↩︎

For each method and each radius threshold \(R\), we use different sigma selection strategies depending on the metric. For certified accuracy, we select the \(\sigma\) value that maximizes certified accuracy at that threshold. For mean distance, we select the \(\sigma\) value that minimizes mean distance (i.e., gives the best prediction quality) for certified samples at that threshold.

We evaluate certified accuracy on 100 stratified test samples (10 per digit class) from the rotated MNIST test set. For each sample, we compute the certified radius \(r_i\) using the certification method, compute the clean model prediction \(\hat{y}_i\) (smoothed model output), extract the true label \(y_i\) from the rotated MNIST dataset, and compute the prediction error \(d(\hat{y}_i, y_i)\) using circular distance. We then check if the sample is both correct (\(d(\hat{y}_i, y_i) \leq 10^\circ\)) and certified (\(r_i \geq R\)), and count the fraction of samples satisfying both conditions. This evaluation protocol uses the true label \(y_i\) (from the rotated dataset) rather than the model output, ensuring that we measure accuracy relative to ground truth rather than self-consistency.

17.6 Optimization-Based Radius Computation Methodology↩︎

The optimization-based radius provides an empirical estimate of the true worst-case certified radius by using white-box optimization to find adversarial perturbations. While this is not a provable guarantee, it provides a useful lower bound for evaluating certificate tightness.

17.6.1 Problem Formulation↩︎

For a smoothed function \(g_\sigma(x) = \mathbb{E}_{\eta \sim \mathcal{N}(0, \sigma^2 I)}[f(x + \eta)]\), the optimization-based radius is the largest radius \(R\) such that: \[\max_{\|\delta\|_2 \leq R} |g_\sigma(x + \delta) - g_\sigma(x)| \leq \epsilon,\] where \(\epsilon\) is the output tolerance. This is computed using projected gradient descent (PGD) to find worst-case perturbations, combined with Monte Carlo estimation to evaluate the smoothed function.

17.6.2 Algorithm Overview↩︎

The computation proceeds in three steps:

Step 1: Baseline Estimation. We compute the clean baseline \(g_\sigma(x)\) using Monte Carlo estimation with \(N_{\text{mc}} = 50,000\) samples: \[g_\sigma(x) \approx \frac{1}{N_{\text{mc}}} \sum_{i=1}^{N_{\text{mc}}} f(x + \eta_i),\] where \(\eta_i \sim \mathcal{N}(0, \sigma^2 I)\) are noise samples. For angle predictions, which are in the range \([-\pi, \pi]\) radians (\([-180^\circ, 180^\circ]\) degrees), we use circular averaging to handle the periodic nature of angles. Specifically, we compute the mean of cosine and sine values separately: \(\bar{\cos} = \frac{1}{N_{\text{mc}}}\sum_{i=1}^{N_{\text{mc}}} \cos(\theta_i)\) and \(\bar{\sin} = \frac{1}{N_{\text{mc}}}\sum_{i=1}^{N_{\text{mc}}} \sin(\theta_i)\), where \(\theta_i\) are the predicted angles in radians. We then use the two-argument arctangent function \(\text{arctan2}(\bar{\sin}, \bar{\cos})\) to recover the mean angle. The \(\text{arctan2}\) function computes the angle in the correct quadrant by taking both sine and cosine as inputs, ensuring that angles near the boundaries of \([-\pi, \pi]\) are correctly averaged (e.g., averaging \(-175^\circ\) and \(175^\circ\) gives \(180^\circ\) rather than \(0^\circ\)).

Step 2: Worst-Case Perturbation Search. For a candidate radius \(R\), we use PGD to find the worst-case perturbation \(\delta^*\): \[\begin{align} \delta^{(t+1)} &= \Pi_{\|\delta\|_2 \leq R}\left(\delta^{(t)} + \alpha \nabla_\delta |g_\sigma(x + \delta^{(t)}) - g_\sigma(x)|\right), \end{align}\] where \(\Pi_{\|\delta\|_2 \leq R}\) projects onto the \(\ell_2\) ball of radius \(R\), and \(\alpha\) is the step size. At each iteration, we estimate \(g_\sigma(x + \delta^{(t)})\) using Monte Carlo with the same noise samples (Common Random Numbers, CRN) to reduce variance.

Step 3: Binary Search. We perform binary search on \(R\) to find the largest radius where the worst-case output change \(\leq \epsilon\). We initialize \(R_{\text{low}} = 0\) and \(R_{\text{high}} = R_{\max}\) (typically \(5.0\) pixels for MNIST experiments). While \(R_{\text{high}} - R_{\text{low}} > \text{tol}\) (where \(\text{tol} = 10^{-3}\) pixels), we compute \(R_{\text{mid}} = (R_{\text{low}} + R_{\text{high}})/2\), evaluate the worst-case change at \(R_{\text{mid}}\) using PGD, and update the bounds: if the change \(\leq \epsilon\), we set \(R_{\text{low}} = R_{\text{mid}}\); otherwise, we set \(R_{\text{high}} = R_{\text{mid}}\). The algorithm returns \(R_{\text{low}}\) as the optimization-based radius. The binary search typically converges in 10-15 iterations.

17.6.3 Implementation Details↩︎

We use \(N_{\text{mc}} = 50,000\) Monte Carlo samples, validated to provide stable estimates of the smoothed function. To reduce variance, we employ Common Random Numbers (CRN): the same noise bank is reused across all evaluations for a given sample.

For the PGD optimization, we use \(n_{\text{restarts}} = 5\) restarts to avoid local minima, with \(n_{\text{steps}} = 50\) steps per restart and a step size of \(\alpha = 0.01\) (adaptive based on radius). The perturbation is projected onto the \(\ell_2\) ball and clipped to the \([0, 1]\) pixel range to ensure valid image values.

The binary search is initialized with range \([0, R_{\max}]\) where \(R_{\max} = 5.0\) pixels for MNIST experiments, using a tolerance of \(\text{tol} = 10^{-3}\) pixels. The search typically converges in 10-15 iterations, well below the maximum of 20 iterations.

17.6.4 Limitations and Assumptions↩︎

The optimization-based radius is an approximation with several limitations. First, since it relies on optimization, PGD may not find the global optimum, so the optimization-based radius may underestimate the true worst-case radius. Second, the smoothed function \(g_\sigma(x)\) is estimated with finite Monte Carlo samples, introducing estimation error. Third, the method requires white-box access, meaning we can query the model many times (typically \(O(N_{\text{mc}} \times n_{\text{restarts}} \times n_{\text{steps}})\) evaluations per sample) to compute gradients and evaluate the smoothed function at different perturbations. Note that this does not require access to the model’s internal parameters or weights; we only need the ability to evaluate the model function \(f(x)\) at arbitrary inputs, which is available during evaluation but not in deployment scenarios. Finally, the computational cost of \(O(N_{\text{mc}} \times n_{\text{restarts}} \times n_{\text{steps}})\) model evaluations per sample makes it expensive for large-scale evaluation. Despite these limitations, the optimization-based radius provides a useful empirical benchmark for evaluating certificate tightness, as it represents the best achievable estimate using optimization-based methods.

17.6.5 Validation↩︎

To validate the optimization-based radius computation, we verify three properties. First, we check soundness: certified radii should be \(\leq\) optimization-based radii (conservative), which is confirmed empirically. Second, we verify convergence: results are stable across different random seeds. Third, we assess sensitivity: results are robust to PGD hyperparameters (step size, restarts). Empirical validation shows that certified radii are consistently \(\leq\) optimization-based radii, confirming that our certificates are sound (conservative) as expected.

17.6.6 Experimental Setup for Tightness Analysis↩︎

We evaluate tightness on the MNIST rotation task using 100 stratified test samples from rotated MNIST (10 per digit class). These same 100 samples are used across all certification methods: variance/gradient estimation, certified radius computation (both \((E, C)+M\) and \((E, C, G)+M\) methods), and alpha-smoothing certification.

We evaluate at a single noise level \(\sigma = 0.5\) due to the high computational cost of optimization-based radius computation, with output tolerance \(\epsilon = 10^{\circ}\) (0.1745 radians). The binary search uses a maximum search radius of \(R_{\max} = 5.0\) pixels in raw pixel space, with a tolerance of \(10^{-3}\) pixels. Of the 100 samples, 18 (18.0%) hit the \(R_{\max}\) bound and are excluded from ratio analysis.

For the 82 samples that do not hit the bound, we compute the tightness ratio \(r_i = R_{\text{opt},i} / R_{\text{cert},i}\) for each sample \(i\). The mean ratio reported in Table 3 is the arithmetic mean of individual ratios: \(\bar{r} = \frac{1}{82}\sum_{i=1}^{82} r_i\). Note that this differs from the ratio of means \((\bar{R}_{\text{opt}} / \bar{R}_{\text{cert}})\) due to the nonlinearity of division; we report the mean of ratios as it better represents the typical per-sample tightness. All summary statistics in Table 3 (mean certified, mean optimization-based, mean ratio, median ratio) are computed over the same 82 samples for consistency.

For the optimization-based radius computation, we use \(N_{\text{mc}} = 50,000\) Monte Carlo samples for the baseline estimation and \(N_{\text{attack}} = 2,000\) samples per PGD step for evaluating the smoothed function during optimization. We use 5 PGD restarts with 100 steps per restart to avoid local minima, resulting in approximately \(50,000 + (5 \times 100 \times 2,000) = 1,050,000\) model evaluations per sample. All computations are performed on GPU.

17.6.7 Soundness Analysis: Samples Near the Boundary↩︎

A ratio \(r \geq 1.0\) indicates a sound certificate (certified radius does not exceed the empirically-found worst-case radius), while \(r < 1.0\) suggests potential finite-sample estimation error.

We find 2 samples (2.4% of uncapped samples) with ratio \(< 1.0\): Sample 77 with \(r = 0.39\) and Sample 13 with \(r = 0.65\). For these samples, the certified radius exceeds the optimization-based radius by factors of \(2.6\times\) and \(1.5\times\) respectively, indicating that finite-sample variance and gradient estimates led to overconfident certificates. This is consistent with the probabilistic nature of our confidence intervals: with \(95\%\) confidence level, we expect approximately \(5\%\) of samples to fall outside the bounds, and the observed \(2.4\%\) rate is within this tolerance.

Additionally, 4 samples (4.9% of uncapped) have ratios in the range \([1.1, 1.3]\): Samples 1, 51, 55, and 61. These samples are technically sound (ratio \(\geq 1.0\)) but represent cases where the certificates are relatively tight, with certified radii within \(14\%-30\%\) of the optimization-based radii. This is visible in the tightness histogram (Figure 3) as samples clustered near the \(r = 1.0\) threshold.

In total, 6 out of 82 uncapped samples (\(7.3\%\)) have tightness ratios below \(1.3\), representing the “edge cases” where statistical estimation is most challenging. The remaining 76 samples (\(92.7\%\)) have ratios \(\geq 1.3\), with a mean ratio of \(3.40\times\) and median of \(3.46\times\), indicating substantial conservatism in the certified radii. Including the 18 capped samples (which are definitively sound with optimization-based radius \(\geq 5.0\) pixels exceeding all certified radii), the overall empirical soundness rate is 98% (98 out of 100 samples).

References↩︎

---

1. UTKFace dataset page: https://github.com/aicip/UTKFace.↩︎

2. MiVOLO-v2 model checkpoint: https://huggingface.co/iitolstykh/mivolo_v2↩︎