stochastic method, smoothing method, minEmax problem, Wasserstein distributionally robust optimization, adversarially robust training.

1 Introduction↩︎

This paper studies the following minimization–expectation–maximization (minEmax) problem: \[\label{eq:minmax} \min_{{\mathbf{y}}\in\mathbb{R}^{m_1}}\left\{g({\mathbf{y}}):=\varphi({\mathbf{y}})+\mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}}\left[\max_{{\mathbf{z}}\in\mathcal{Z}} \Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}})\right]\right\}.\tag{1}\] Here \(\mathcal{Z}\subseteq\mathbb{R}^{m_2}\) is a nonempty compact set, \(\mathbb{P}\) is a probability distribution on \(\mathcal{X}\subseteq\mathbb{R}^{m_3}\), and independent samples can be drawn from \(\mathbb{P}\). For each \({\mathbf{x}}\in\mathcal{X}\), define \[\label{eq:def-Phi-i} \Phi({\mathbf{y}};{\mathbf{x}}):=\max_{{\mathbf{z}}\in\mathcal{Z}}\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}}).\tag{2}\] Throughout the paper, we make the following assumptions on problem 1 .

Assumption 1 allows nonsmooth convex regularizers, such as indicator functions of compact convex sets. The feasible set \(\mathcal{Z}\) may be a connected compact set or a finite discrete set; in either case it is compact. By Berge’s maximum theorem [1], Assumption 1 implies that \(\Phi(\cdot;{\mathbf{x}})\) is well-defined and continuous for every \({\mathbf{x}}\in\mathcal{X}\). In addition, Assumption 1(i) imposes differentiability only with respect to the outer variable \({\mathbf{y}}\); no differentiability with respect to \({\mathbf{z}}\) is required, while Assumption 2 gives boundedness of the \({\mathbf{y}}\)-gradient and Lipschitz continuity in the inner variable.

Problem 1 encompasses a wide variety of modern machine learning applications, particularly those involving challenging inner maximization problems that may be nonsmooth in \({\mathbf{z}}\) or may involve a finite feasible set. Prominent examples include adversarially robust training (e.g., [2]–[5]) and distributionally robust optimization (DRO) (e.g., [6], [7]); see Section 1.1 for more details.

When \(\mathbb{P}\) is the uniform distribution on a finite dataset \(\{{\mathbf{x}}_1,\ldots,{\mathbf{x}}_n\}\), problem 1 can be written as \[\min_{{\mathbf{y}}}\;\max_{{\mathbf{z}}_1,\ldots,{\mathbf{z}}_n\in\mathcal{Z}} \left\{\varphi({\mathbf{y}})+\frac{1}{n}\sum_{i=1}^n\Psi({\mathbf{y}},{\mathbf{z}}_i;{\mathbf{x}}_i)\right\}.\] Minimax methods (see, e.g., [8]–[14]) can then be applied in principle, although the reformulation may be expensive for a large \(n\); see, e.g., [15]. For a general distribution \(\mathbb{P}\), an analogous reformulation would require an infinite-dimensional measurable decision rule \({\mathbf{x}}\mapsto {\mathbf{z}}({\mathbf{x}})\), rather than finitely many variables \({\mathbf{z}}_1,\ldots,{\mathbf{z}}_n\). This fundamental distinction underscores the necessity for developing a new optimization framework that can efficiently handle the inherent complexity of expectation-over-maximization problems with general probability distributions.

We therefore view problem 1 as nonsmooth minimization of an expected-value objective. Though a subgradient method is standard for nonsmooth stochastic minimization, its convergence analysis typically requires access to high-quality (stochastic) subgradient evaluations. For problem 1 , accurate inner objective values are not sufficient for this purpose: Example 1 in Appendix 6 shows that a nearly optimal inner maximizer can yield a vector far from the Clarke subdifferential of the outer objective. We avoid this issue by applying the log-mean-exp (LME) smoothing [16] to \(\Phi(\cdot;{\mathbf{x}})\). Specifically, for every \({\mathbf{x}}\in\mathcal{X}\) and \(\mu>0\), define the LME smoothing function \(\widetilde{\Phi}\) as follows: \[\label{eq:smoothphi} \widetilde{\Phi}({\mathbf{y}},\mu;{\mathbf{x}}): =\mu \log \mathbb{E}_{{\mathbf{z}}\sim \zeta} [e^{\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}})/\mu}] =\mu\log\int_{\mathcal{Z}}\exp(\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}})/\mu)\,\zeta(d{\mathbf{z}}),\tag{3}\] where \(\zeta\) denotes a fixed probability measure on \(\mathcal{Z}\). In particular, if \(\mathcal{Z}\) is finite, \(\zeta\) corresponds to the uniform measure; if \(\mathcal{Z}\) is a connected compact set (with positive Lebesgue measure), \(\zeta\) is the normalized Lebesgue measure on \(\mathcal{Z}\) (i.e., \(\zeta({\mathcal{Z}})=1\)).

By leveraging the properties of the LME smoothing function, we effectively mitigate the computational challenges discussed above and develop a stochastic smoothing proximal gradient method (SSPG). Our analysis in Section 3 establishes that every almost sure accumulation point of the sequence generated by SSPG is a Clarke stationary point of problem 1 . Given that the primal objective function is Clarke regular under our assumptions, such accumulation points are additionally directional stationary [17]–[19]. We further establish the iteration complexity results of SSPG to find an approximate stationary point under two different notions.

The remainder of this section proceeds as follows: we first describe several relevant applications structured as 1 , then summarize our main technical contributions, and finally introduce necessary notation and definitions.

1.1 Applications↩︎

Two representative applications with the structure of 1 are Wasserstein distributionally robust optimization (WDRO) and adversarially robust training.

1.1.1 WDRO↩︎

The WDRO problem [6], [20] can be formulated as \[\label{eq:model} \min_{{\boldsymbol{\theta}}\in \Theta} \max _{\widehat{\mathbb{P}} \in \mathcal{B}_\delta(\mathbb{P})} \mathbb{E}_{{\mathbf{x}}\sim \widehat{\mathbb{P}}}\left[\ell({\boldsymbol{\theta}},{\mathbf{x}})\right],\tag{4}\] where \(\mathbb{P}\) is the nominal distribution on the uncertainty set \(\mathcal{Z}\), \(\widehat{\mathbb{P}}\) ranges over alternative distributions in the Wasserstein ambiguity set, \(\ell:\mathbb{R}^{m_1}\times \mathbb{R}^{m_2}\mapsto \mathbb{R}\) is a given (possibly nonconvex) loss function, \(\Theta\) is a closed convex set, and the ambiguity set \(\mathcal{B}_\delta(\mathbb{P})\) is defined as the \(\delta\)-ball in the \(p\)-th Wasserstein distance centered at \(\mathbb{P}\), i.e., \(\mathcal{B}_\delta(\mathbb{P})=\{\widehat{\mathbb{P}} \in \mathcal{P}(\mathcal{Z})\mid d_{\mathcal{W}_p}(\widehat{\mathbb{P}}, \mathbb{P}) \leq \delta\}.\) Here, \(\mathcal{P}(\mathcal{Z})\) is the space of probability distributions \(\widehat{\mathbb{P}}\) supported on \(\mathcal{Z}\) with \(\mathbb{E}_{{\mathbf{x}}\sim \widehat{\mathbb{P}}} [\|{\mathbf{x}}\|] <\infty\), and the \(p\)-th Wasserstein distance [21] between distributions \(\mathbb{Q}_1, \mathbb{Q}_2 \in \mathcal{P}(\mathcal{Z})\) is defined by \[\begin{align} &d_{\mathcal{W}_p}\left(\mathbb{Q}_1, \mathbb{Q}_2\right):= \\ &\inf \left\{\left(\int_{\mathcal{Z}\times \mathcal{Z}}\left\|{\mathbf{z}}_1-{\mathbf{z}}_2\right\|_p^p \mathrm{\Pi}\left(\mathrm{d} {\mathbf{z}}_1, \mathrm{d} {\mathbf{z}}_2\right)\right)^{1/p} \bigg| \begin{array}{l} \mathrm{\Pi} \text{ is a joint distribution of } {\mathbf{z}}_1 \text{ and } {\mathbf{z}}_2 \text{ with}\\ \text{marginal distributions } \mathbb{Q}_1 \text{ and } \mathbb{Q}_2 \text{, respectively} \end{array}\right\}. \end{align}\] For a fixed decision \({\boldsymbol{\theta}}\), the inner problem in 4 is the worst-case risk: \[\label{eq:worstrisk} \max_{\widehat{\mathbb{P}} \in \mathcal{B}_\delta(\mathbb{P})} \mathbb{E}_{{\mathbf{x}}\sim \widehat{\mathbb{P}}}\left[\ell({\boldsymbol{\theta}},{\mathbf{x}})\right].\tag{5}\] Under some mild conditions [22], [23], the worst-case risk in 5 is finite and attainable, and furthermore, strong duality holds. The dual problem of 5 is given by [7], [22], [23] \[\label{eq:worstriskdual} \min_{\lambda\geq 0} \lambda \delta^p+\mathbb{E}_{{\mathbf{x}}\sim \mathbb{P}}\left[\max _{{\mathbf{z}}\in \mathcal{Z}}\{\ell({\boldsymbol{\theta}},{\mathbf{z}})-\lambda d({\mathbf{x}},{\mathbf{z}})\}\right].\tag{6}\] Here, \(d: \mathbb{R}^{m_2}\times \mathbb{R}^{m_2}\mapsto \mathbb{R}\) denotes the transport cost of the Wasserstein metric of order \(p\in \mathbb{N}_+\) defined as \(d({\mathbf{z}}_1,{\mathbf{z}}_2):=\|{\mathbf{z}}_1-{\mathbf{z}}_2\|_p^p\) and \(\lambda\in\mathbb{R}_+\) is the Lagrangian multiplier with respect to the inequality constraint \(d_{\mathcal{W}_p}^p(\widehat{\mathbb{P}}, \mathbb{P}) \leq \delta^p\). By strong duality, the optimal values of the aforementioned two models are the same. The following lemma gives the compact interval for the dual multiplier; its proof is deferred to Appendix 7.

Hence, problem 6 can be restricted to \(\lambda\in[0,B_\lambda]\) without loss of optimality, where \(B_{\lambda}= L\,C_{p,m_2}\,\delta^{-(p-1)}\). Substituting this dual representation into 4 , and explicitly imposing the constraint \(\lambda\in[0,B_{\lambda}]\), we arrive at the following reformulation: \[\label{eq:model2} \min_{{\boldsymbol{\theta}}\in \Theta, \lambda\in[0,B_{\lambda}]} g({\boldsymbol{\theta}}, \lambda):= \lambda \delta^p+\mathbb{E}_{{\mathbf{x}}\sim \mathbb{P}}\left[\max _{{\mathbf{z}}\in \mathcal{Z}}\{\ell({\boldsymbol{\theta}},{\mathbf{z}})-\lambda d({\mathbf{x}},{\mathbf{z}})\}\right].\tag{7}\] This problem has the form 1 , with \({\mathbf{y}}=({\boldsymbol{\theta}},\lambda)\) and \(\Psi(({\boldsymbol{\theta}},\lambda),{\mathbf{z}};{\mathbf{x}})=\ell({\boldsymbol{\theta}},{\mathbf{z}})-\lambda d({\mathbf{x}},{\mathbf{z}}).\) It satisfies Assumption 1 when \(\Theta\), \([0,B_\lambda]\), and \(\mathcal{Z}\) are compact.

1.1.2 Adversarially robust training↩︎

In a classification setting, the adversarially robust training [2]–[5] problem takes the form of \[\label{eq:adt}\min _{{{\boldsymbol{\theta}}\in \Theta}}\frac{1}{n}\sum_{i=1}^n \max _{{\mathbf{z}}\in \Delta({{\mathbf{x}}}_i)} \ell\left({\bar{{\mathbf{x}}}_i}, f({\boldsymbol{\theta}}, {\mathbf{z}})\right).\tag{8}\] Here, \(\{({\mathbf{x}}_i,{\bar{{\mathbf{x}}}}_i)\}_{i=1}^n\subset \mathbb{R}^{m_2}\times \mathbb{R}^{m_3}\) are sampled ordered pairs of data points and their corresponding labels (respectively), \(\Delta({{\mathbf{x}}})=\left\{{\mathbf{z}}\in\mathbb{R}^{m_2}\big| d\left({{\mathbf{x}}}, {\mathbf{z}}\right) \leq \epsilon\right\}\) represents a set of feasible perturbations applied to a given data point \({\mathbf{x}}\) based on some distance function \(d\), \(f:\mathbb{R}^{m_1}\times \mathbb{R}^{m_2}\mapsto \mathbb{R}^{m_3}\) is a prediction function, and \(\ell:~\mathbb{R}^{m_3}\times \mathbb{R}^{m_3}\mapsto \mathbb{R}\) is a (possibly nonconvex) loss function. In the context of computer vision, this model seeks to identify the worst-case perturbations for images within a prescribed radius [3], [4]. Let \(\mathbb{P}_n\) be the empirical distribution on the datasets. Then \[\frac{1}{n}\sum_{i=1}^n \max _{{\mathbf{z}}\in \Delta({{\mathbf{x}}}_i)} \ell\left({\bar{{\mathbf{x}}}_i}, f({\boldsymbol{\theta}}, {\mathbf{z}})\right) = \mathbb{E}_{({\mathbf{x}},\bar{\mathbf{x}})\sim \mathbb{P}}\max _{{\mathbf{z}}\in \Delta({{\mathbf{x}}})} \ell\left({\bar{{\mathbf{x}}}}, f({\boldsymbol{\theta}}, {\mathbf{z}})\right),\] and thus 8 can be formulated into 1 . This formulation satisfies Assumption 1 whenever \(\Theta\) is compact, \(\ell(\cdot,\cdot)\) is continuous, and \(\ell(\bar{\mathbf{x}},\cdot)\) and \(f(\cdot,{\mathbf{z}})\) are smooth for all \(\bar{\mathbf{x}}\) and \({\mathbf{z}}\). Indeed, adversarially robust training can be viewed as a dual problem of the \(\infty\)-WDRO problem whose ambiguity set is characterized by the Wasserstein distance with \(p=+\infty\); see [24]. Notice that the dual formulation in problem 7 applies for \(p\in [1,\infty)\) but not for \(p=\infty\).

1.2 Contributions↩︎

This paper develops a stochastic smoothing framework for the minEmax problem 1 , a nonsmooth expected-value problem in which the outer variable is to minimize the expectation of a pointwise maximum. This formulation covers, among other examples, Wasserstein distributionally robust optimization (WDRO) after dualization of the worst-case risk and adversarially robust training. Our main contributions are summarized below.

1. We introduce the log-mean-exp (LME) smoothing in 3 for each \({\mathbf{x}}\in{\mathcal{X}}\). We prove that \(\widetilde{\Phi}(\cdot,\mu;{\mathbf{x}})\) is a valid smoothing function of \(\Phi(\cdot;{\mathbf{x}})\) defined in 2 , preserves the basic Lipschitz bound in \({\mathbf{y}}\), and has a \({\mathbf{y}}\)-gradient Lipschitz constant of order \(L_\Psi+l_\Psi^2/\mu\). We also derive approximation-gap bounds for both finite \({\mathcal{Z}}\) and connected compact \({\mathcal{Z}}\). These bounds make it possible to translate stationarity of the smoothed objectives into Goldstein stationarity of problem 1 .

2. Using the LME smoothing, we propose SSPG, a stochastic smoothing proximal gradient method. At iteration \(k\), SSPG uses a stochastic estimator of \(\nabla_{{\mathbf{y}}}\mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}}[\widetilde{\Phi}({\mathbf{y}},\mu_k;{\mathbf{x}})]\) and performs one proximal gradient step with respect to the nonsmooth convex term \(\varphi\), while the smoothing parameter \(\mu_k\) is kept fixed or decreased according to a prescribed rule. We establish a nonasymptotic stationarity bound for a randomly selected iterate. In particular, SSPG obtains an \(\epsilon\)-scaled stationary point in expectation with iteration complexity \(O(\epsilon^{-3})\) and sample complexity \(O(\epsilon^{-5})\). After converting scaled stationarity to Goldstein stationarity, SSPG obtains an \((O(\epsilon),O(\epsilon))\)-Goldstein stationary point with iteration complexity \(\widetilde{O}(\epsilon^{-4})\) and sample complexity \(\widetilde{O}(\epsilon^{-6})\). With a summable refinement schedule, every almost-sure cluster point of the resulting sequence of approximate stationary points is Clarke stationary, and hence directional stationary, for 1 .

3. For WDRO, we prove an explicit global upper bound on the optimal dual multiplier in the Wasserstein worst-case-risk reformulation. This allows the dual variable to be restricted to a compact interval without loss of optimality, making the compact-domain assumptions compatible with WDRO analysis and implementation. We then test SSPG on newsvendor, robust regression, \(\infty\)-Wasserstein/adversarially learning, and adversarially robust deep-learning tasks. The experiments show that SSPG is competitive with the tested baselines and, in several settings, improves the accuracy–runtime tradeoff relative to nonsmoothed or fixed-entropic alternatives.

1.3 Definitions and Notations↩︎

Let \(\|\cdot\|\) be the 2-norm of a vector. We use \([M]\) to represent \(\{1,2,\ldots,M\}\) for an integer \(M\). Given a compact set \(\mathcal{Z}\), denote \(|\mathcal{Z}|\) as the cardinality of \(\cal{Z}\) if it is a finite set and volume of \(\cal{Z}\) if it is a continuous set. The indicator function is denoted as \(\mathbb{I}.\) We refer to \(\mathbf{1}_m\) as the \(m\)-dimensional vector of all ones. We write \(\operatorname{Proj}\), \(\operatorname{Prox}\), and \(\operatorname{dist}\) for the projection, proximal, and distance operators, respectively. The normal cone of a convex set \(\mathcal{Y}\) at \({\mathbf{y}}^*\) is defined as \(\mathcal{N}_{\mathcal{Y}}({\mathbf{y}}^*)=\{{\boldsymbol{\gamma}}:\langle {\boldsymbol{\gamma}}, {\mathbf{y}}-{\mathbf{y}}^*\rangle \leq 0, \forall {\mathbf{y}}\in \mathcal{Y}\}.\)

Let \(h:\mathbb{R}^d\to\mathbb{R}\) be locally Lipschitz continuous. We denote by \(h'(\bar{{\mathbf{y}}};{\mathbf{d}})\) the directional derivative of \(h\) at \(\bar{{\mathbf{y}}}\) in direction \({\mathbf{d}}\): \[h^{\prime}\left(\bar{{\mathbf{y}}} ; {\mathbf{d}}\right)=\lim _{ t \downarrow 0} \frac{1}{t}(h(\bar{{\mathbf{y}}}+t {\mathbf{d}})-h(\bar{{\mathbf{y}}})).\] If \(h^{\prime}\left(\bar{{\mathbf{y}}} ; {\mathbf{d}}\right)\) is well-defined for any unit vector \({\mathbf{d}}\), we say \(h\) is directionally differentiable at \(\bar{{\mathbf{y}}}\). It is known that if \(h\) is piecewise differentiable and Lipschitz continuous, then it is directionally differentiable [19], [25]. Denote \(h^{\circ}(\bar{ {\mathbf{y}}} ; {\mathbf{d}})\) as the generalized directional derivative at \(\bar{{\mathbf{y}}}\) along the direction \({\mathbf{d}}\) [26], by \[h^{\circ}(\bar{ {\mathbf{y}}} ; {\mathbf{d}}):=\limsup _{{{\mathbf{y}}\rightarrow \bar{ {\mathbf{y}}}}, {t \downarrow 0}} \frac{1}{t}(h({\mathbf{y}}+t {\mathbf{d}})-h({\mathbf{y}})).\] In general, it holds \(h^{\circ}(\bar{ {\mathbf{y}}} ; {\mathbf{d}})\geq h^{\prime}\left(\bar{{\mathbf{y}}} ; {\mathbf{d}}\right)\) [26], and \(h^{\prime}\left(\bar{{\mathbf{y}}} ; {\mathbf{d}}\right)=h^{\circ}(\bar{ {\mathbf{y}}} ; {\mathbf{d}})={\mathbf{d}}^{\top}\nabla h(\bar{ {\mathbf{y}}})\) if \(h\) is differentiable. We use \(\partial h\) to denote the (Clarke) subdifferential [26] of a continuous function \(h\), i.e., \[\partial h(\bar{{\mathbf{y}}}) := \left\{ {\boldsymbol{\xi}}\in\mathbb{R}^d: \langle {\boldsymbol{\xi}},{\mathbf{d}}\rangle\le h^\circ(\bar{{\mathbf{y}}};{\mathbf{d}}) \;\text{for all }{\mathbf{d}}\in\mathbb{R}^d \right\}.\] For \(\mu\ge 0\), we define the Goldstein \(\mu\)-subdifferential [27] of \(h\) at \({\mathbf{y}}\) by \(\partial^{\mu} h({\mathbf{y}}) := \operatorname{conv} \left( \bigcup_{{\mathbf{u}}\in \mathbb{B}({\mathbf{y}},\mu)} \partial h({\mathbf{u}}) \right),\) where \(\mathbb{B}({\mathbf{y}},\mu)=\{{\mathbf{u}}:\|{\mathbf{u}}-{\mathbf{y}}\|\le \mu\}\).

From Definition 1, if \(h\) is convex or differentiable, then it is Clarke regular [26] and directionally differentiable.

For a locally Lipschitz function, directional stationarity implies Clarke stationarity. If \(h\) is Clarke regular, the converse also holds. Thus the two notions coincide for Clarke regular objectives. The \((\mu,\epsilon)\)-Goldstein condition is a finite-accuracy relaxation: it aggregates Clarke subgradients in a \(\mu\)-neighborhood of the candidate point and requires the resulting convex hull to contain a vector of norm at most \(\epsilon\). It is particularly suitable for nonasymptotic analysis of nonsmooth nonconvex stochastic optimization, when exact Clarke stationarity is difficult to certify. As \(\mu\downarrow 0\) and \(\epsilon\downarrow 0\), Goldstein stationarity recovers Clarke stationarity under outer-semicontinuity conditions.

Approximate stationarity conditions based on Goldstein subdifferentials have recently emerged as tools in the nonasymptotic analysis of nonsmooth nonconvex optimization. For instance, [29] established connections between “uniform smoothing" and Goldstein stationarity, and developed both deterministic and stochastic gradient-free methods with corresponding nonasymptotic convergence guarantees. The LME smoothing approach employed in this paper differs from the”uniform smoothing" technique introduced by [29]. Specifically, our method smooths the function \(\Phi\), which incorporates the maximum operator within an expectation-over-maximization objective, a scenario not covered by the framework of “uniform smoothing”. Nevertheless, Goldstein stationarity serves a convenient tool in our analysis.

1.4 Organization↩︎

The remainder of the paper is organized as follows. Section 2 provides a review of the related literature. In Section 3, we introduce the SSPG framework for solving the problem 1 . We first construct a smoothing function, then give our proposed algorithmic framework, and finally provide the convergence analysis, with proofs given in Appendix 9. Numerical results are presented in Section 4. We conclude the paper in Section 5.

2 Related Work↩︎

In this section, we briefly survey several optimization approaches for WDRO problems. Additionally, we review smoothing methods and the LME smoothing function, which serve as foundational tools for the proposed framework.

2.1 Existing Methods for Solving WDRO↩︎

As introduced in [7], a widely adopted strategy in the distributionally robust optimization (DRO) literature is to solve problem 7 . Standard minimax solvers are applicable when \(\mathbb{P}\) is an empirical distribution. However, for general distributions, practical algorithms typically require additional structural conditions. Examples include: (i) discretization of the support set \(\mathcal{Z}\) [30]–[33]; (ii) representation of \(\ell({\boldsymbol{\theta}},\cdot)\) as a finite maximum of concave functions under an \(\ell_1\) transport cost [23], [34]; (iii) log-loss models, such as robust logistic regression [35]–[37]; or (iv) strong concavity of \(\ell({\boldsymbol{\theta}},\cdot)-\lambda d(\cdot,{\mathbf{x}})\) [38], [39]. Without these structural assumptions, solving WDRO problems computationally becomes significantly challenging.

Recognizing the inherent challenges in addressing the WDRO problem, [40] investigate a regularized version wherein entropic smoothing yields a sampled approximation to the original objective. They demonstrate convergence of approximate gradients toward Clarke subgradients of the unregularized WDRO objective as the regularization parameter tends toward zero. This entropic-regularization perspective has also influenced robust learning software development, notably in [41], [42]. In their analyses, the authors explicitly distinguish between regularized WDRO problem and the original WDRO problem, underscoring that their methods do not fundamentally solve the original WDRO problem.

[43] innovatively proposes addressing a dual problem of Sinkhorn DRO (SDRO), which can be written as \[\label{eq:smoothg} \min_{{\boldsymbol{\theta}}\in \Theta, \lambda\in[0,B_{\lambda}]}g_{\mathrm{s}}({\boldsymbol{\theta}}, \lambda) := \lambda \delta^p+ \lambda\eta \mathbb{E}_{{\mathbf{x}}\sim \mathbb{P}}\left[\log \mathbb{E}_{{\mathbf{z}}\sim \zeta}\left[e^{(\ell({\boldsymbol{\theta}},{\mathbf{z}}) -\lambda d({\mathbf{x}},{\mathbf{z}})) / \lambda\eta}\right]\right]\tag{9}\] for some \(\eta>0\), where \(B_{\lambda}\) is the same as that in 7 , and \(\zeta\) is a distribution supported on \({\mathcal{Z}}\). This perspective connects naturally to regularized WDRO [40], since Sinkhorn distance is itself a Wasserstein-type discrepancy based on entropic regularization. The objective function in 9 can be viewed as a smoothing function of 7 by Definition 3 with \(\eta\) as the smoothing parameter. Assuming \(\ell(\cdot,{\mathbf{z}})\) is convex for all feasible \({\mathbf{z}}\), [43] develop a convergent triple-loop algorithm by skillfully combining a bisection method, a stochastic mirror descent method, and multilevel Monte-Carlo simulation. However, they fix \(\eta>0\). In addition, in their complexity result [43], they assume that any optimal solution \(({\boldsymbol{\theta}}^*,\lambda^*)\) of problem 9 satisfies \(\lambda^*\ge \overline{\lambda}\) for some positive scalar \(\overline{\lambda}\). This assumption circumvents significant computational challenges by preventing both the Lipschitz constant and the gradient Lipschitz constant from exploding if \(\lambda\) approaches 0 as the algorithm progresses, in which case, an exponential increase in the sampling complexity of their multilevel Monte-Carlo simulations would occur with respect to \(\frac{1}{\lambda}\) (see the proof in [43]). Hence, their algorithm does not solve the original WDRO problem.

In contrast, when applied to WDRO, our proposed method for 1 smooths the function in 7 by the LME function and solves \[\label{eq:smoothg2} \min_{{\boldsymbol{\theta}}\in \Theta, \lambda\in[0,B_{\lambda}]} \widetilde{g}({\boldsymbol{\theta}}, \lambda, \mu) := \lambda \delta^p + \mu \mathbb{E}_{{\mathbf{x}}\sim \mathbb{P}}\left[\log \mathbb{E}_{{\mathbf{z}}\sim \zeta}\left[e^{(\ell({\boldsymbol{\theta}}, {\mathbf{z}}) -\lambda d({\mathbf{x}}, {\mathbf{z}})) / \mu}\right]\right].\tag{10}\] In contrast to [43], our formulation allows the dual variable \(\lambda\) to approach zero. The formulation in 10 replaces \(\lambda\eta\) in 9 with \(\mu\). However, this subtle difference in the formulation leads to a fundamentally different method from the algorithm in [43]. First, [43] solves a dual problem of SDRO with a fixed \(\eta\). In contrast, we solve WDRO directly, requiring \(\mu \downarrow 0\) in our algorithmic framework. This distinction changes the limiting problem and requires controlling the smoothing bias as \(\mu\downarrow0\), while the gradient Lipschitz constants of the smoothed objectives typically grow as \(\mu\) decreases. Second, \(\lambda\) is a primal variable in problem 9 and appears in the denominator of the exponential term. This results in both the Lipschitz constant and the gradient Lipschitz constant of \(g_s\) with respect to \(\lambda\) becoming unbounded as \(\lambda\) approaches zero, which potentially yields slow convergence and high complexity. In contrast, the Lipschitz constant of \(\widetilde{g}\) with respect to \(\lambda\) is bounded, and the gradient Lipschitz constant of \(\widetilde{g}\) with respect to \(\lambda\) can be bounded by \(\frac{C}{\mu}\) for some positive scalar \(C\).

2.2 Smoothing Methods and the LME Smoothing Function↩︎

Smoothing methods provide a powerful and widely-used framework for solving nonsmooth optimization problems; we refer readers to foundational discussions and examples in [28], [44]–[48]. Typically, smoothing methods involve three essential steps: (i) constructing smooth approximations \(\{\widetilde{g}(\cdot,\mu)\}_{\mu>0}\) to the nonsmooth function \(g\), parameterized by a smoothing parameter \(\mu>0\); (ii) designing algorithms to solve the smoothed problems \(\min_{\mathbf{y}}\widetilde{g}({\mathbf{y}},\mu_k)\) approximately for a sequence \(\{\mu_k\}\); and (iii) analyzing the convergence behavior of solutions as \(\mu_k\downarrow 0\). A theoretical challenge occurs as \(\mu\downarrow 0\): the gradients \(\nabla\widetilde{g}(\cdot,\mu)\) are not Lipschitz continuous with constants bounded uniformly over \(\mu\in(0,1]\); in many smoothing schemes the Lipschitz constant scales as \(\Theta(1/\mu)\).

For deterministic nonsmooth optimization, the theoretical foundations, including epi-convergence, iteration complexity analysis, and gradient bounds, are well-established. In contrast, the theory for stochastic smoothing methods remains underdeveloped. Significant challenges in the stochastic scenario include the intricate coupling between smoothing-induced bias (controlled by \(\mu\)) and variance in stochastic gradient estimators, and the rigorous justification of interchangeability between expectation and differentiation, requiring careful considerations of measurability and integrability. Recently, [49] analyzed convergence properties of stochastic smoothing methods under convexity of \(\Psi\), and bounded variance of the stochastic gradient estimators for \(g(\cdot,\mu)\) for every fixed smoothing parameter \(\mu>0\). Their analysis, however, does not establish convergence to Clarke stationary points almost surely, a critical theoretical guarantee in smoothing methods. Without these restrictive conditions, our paper directly addresses this gap by providing such convergence guarantees (see Theorem 1).

We now discuss the LME smoothing function given in 3 . When the reference measure \(\zeta\) is uniform on a finite set \(\mathcal{Z}=\{{\mathbf{z}}_1,\ldots,{\mathbf{z}}_q\}\), the LME smoothing becomes \(\widetilde{\Phi}({\mathbf{y}},\mu;{\mathbf{x}}) = \mu\log\left(\frac{1}{q}\sum_{j=1}^q \exp(\Psi({\mathbf{y}},{\mathbf{z}}_j;{\mathbf{x}})/\mu)\right),\) which is the log-sum-exp smoothing of the pointwise maximum up to the additive constant \(-\mu\log q\). This well-known LSE function, also termed the softmax maximum [50] or the Neural Networks smoothing function [28], [51], has been extensively applied in finite minimax optimization contexts [49], [52]–[55], particularly when a differentiable approximation of the maximum operator is desired. Its gradient coincides exactly with the softmax function commonly employed in neural network models [50], and extensive studies have documented its robust theoretical properties and practical effectiveness.

While the LME function generalizes LSE to more general settings, including connected compact sets of positive Lebesgue measure, its theoretical properties and convergence analyses differ substantially. Unlike LSE, whose function value can be directly calculated, the LME value generally requires numerical integration or sampling, except in special cases where the exponential integral is available in closed form. Hence, convergence analyses of stochastic smoothing methods employing the LME smoothing require separate and more intricate mathematical treatments, which have not been previously developed in the literature.

3 SSPG: A Stochastic Smoothing Proximal Gradient Framework for Solving Problem 1↩︎

We focus on solving the primal problem \(\min_{{\mathbf{y}}} g({\mathbf{y}})\) in 1 . We first show a few important properties of the smoothing function \(\widetilde{\Phi}(\cdot,\mu;{\mathbf{x}})\) of \(\Phi(\cdot;{\mathbf{x}})\) for each \({\mathbf{x}}\in{\mathcal{X}}\) in Section 3.1. With access to a stochastic gradient estimator, we then introduce the SSPG method in Section 3.2 and establish its convergence results in Section 3.3.

3.1 The LME Smoothing Function↩︎

For every \({\mathbf{x}}\in\mathcal{X}\), we define the LME smoothing function by 3 . The reference probability measure \(\zeta\) is fixed throughout the algorithm: it is the uniform counting measure when \(\mathcal{Z}\) is finite, and the normalized Lebesgue measure when \(\mathcal{Z}\) has positive Lebesgue measure. We make the following assumption to ensure the smoothness of \(\widetilde{\Phi}(\cdot,\mu;{\mathbf{x}})\).

The following lemma shows that \(\widetilde{\Phi}(\cdot,\cdot;{\mathbf{x}})\) is a smoothing function of \(\Phi(\cdot;{\mathbf{x}})\) and establishes key properties of the smoothing function. In particular, it shows that \(\widetilde{\Phi}({\mathbf{y}},\mu;{\mathbf{x}})\) is nonincreasing with respect to \(\mu\). The proof of this lemma is provided in Appendix 9.1.

For any \(\mu_1,\mu_2>0\) and all \({\mathbf{x}}\in{\mathcal{X}}\), we establish upper bounds on \(|\widetilde{\Phi}({\mathbf{y}},\mu_1;{\mathbf{x}})-\widetilde{\Phi}({\mathbf{y}},\mu_2;{\mathbf{x}})|\) with respect to \(|\mu_1-\mu_2|\) in the following two lemmas. Their proofs are given in Appendices 9.2 and 9.3, respectively.

Based on the LME function in 3 , we introduce a (nonconvex) “smoothing” problem of 1 as follows: \[\label{eq:smootho} \min_{{\mathbf{y}}} \left\{\widetilde{g}({\mathbf{y}}, \mu) := \varphi({\mathbf{y}}) + \mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}} [\widetilde{\Phi}({\mathbf{y}}, \mu;{\mathbf{x}})]\right\}.\tag{11}\] It is important to note that problem 11 corresponds precisely to the subproblem solved at each iteration of our proposed SSPG method in Section 3.2, with the smoothing parameter \(\mu\) dynamically decreasing toward zero as the iterations proceed.

Definition 4 is particularly useful in smoothing methods as it quantifies stationarity for the smoothing problem \(\min_{{\mathbf{y}}}\widetilde{g}({\mathbf{y}},\mu)\). However, it does not directly characterize stationarity for the original problem \(\min_{{\mathbf{y}}}g({\mathbf{y}})\). Furthermore, providing nonasymptotic certification of Clarke stationarity remains generally infeasible. We therefore directly work with the \((\mu,\epsilon)\)-Goldstein stationarity notion introduced in Definition 2, and establish in Lemma 5 a conversion from scaled stationarity to Goldstein stationarity. The proof to Lemma 5 is given in Appendix 9.4.

For \(\mu\in(0,1]\), define the LME approximation gap \[\Delta_\mu := \sup_{{\mathbf{y}}\in\operatorname{dom}(\varphi)} \mathbb{E}_{{\mathbf{x}}\sim P} \left[ \Phi({\mathbf{y}};{\mathbf{x}})-\widetilde{\Phi}({\mathbf{y}},\mu;{\mathbf{x}}) \right].\] By Lemma 5(a, b), \(\widetilde{\Phi}({\mathbf{y}},\mu;{\mathbf{x}})\uparrow \Phi({\mathbf{y}};{\mathbf{x}})\) as \(\mu\downarrow0\), and hence \(\Delta_\mu\ge0\). Moreover, we have \[\label{eq:omega}\Delta_\mu \le \omega(\mu) := \begin{cases} 2\mu\log|{\mathcal{Z}}| , & \text{if } {\mathcal{Z}}\text{ is finite discrete,} \\[1.5mm] 2\mu\left[ l_\Psi + m_2\log\left(\dfrac{1+2D_{\mathcal{Z}}}{\mu}\right) \right], & \text{if } {\mathcal{Z}}\text{ is connected compact with diameter }D_{\mathcal{Z}}, \end{cases}\tag{12}\] where the finite case follows from Lemma 3 by taking \(\mu_1=\mu\) and \(\mu_2\downarrow 0\), and the connected compact case follows from Lemma 4 by taking \(\mu_1=\mu\) and letting \(\mu_2\downarrow0\).

We end this section with an almost-sure convergence result, whose proof is given in Appendix 9.5.

3.2 The SSPG method↩︎

In this subsection, we present a stochastic smoothing proximal gradient (SSPG) method for solving problem 1 . At iteration \(k\), the SSPG method approximately solves the smoothed subproblem 11 with smoothing parameter \(\mu=\mu_k\) via a single projected stochastic proximal gradient step. The smoothing parameter is then updated according to a predefined nonincreasing rule. Specifically, at iteration \(k\), we sample \(M_k\) independent and identically distributed points \(\{{\mathbf{x}}_{k_j}\}_{j=1}^{M_k}\) from the distribution \(\mathbb{P}\). Under standard bounded-variance assumptions, stochastic methods typically approximate the gradient \(\mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}}[\nabla_{{\mathbf{y}}}\widetilde{\Phi}({\mathbf{y}}^{(k)},\mu_k;{\mathbf{x}})]\) using the sample average \(\frac{1}{M_k}\sum_{j=1}^{M_k}\nabla_{{\mathbf{y}}}\widetilde{\Phi}({\mathbf{y}}^{(k)},\mu_k;{\mathbf{x}}_{k_j})\). However, as the proof of Lemma 2(b) reveals, exact evaluation of \(\nabla_{{\mathbf{y}}}\widetilde{\Phi}\) generally involves nested expectations and can thus be computationally prohibitive. To address this issue, we instead assume access to gradient estimators \(\mathcal{G}_{k_j}(\cdot,\cdot)\) approximating \(\nabla_{{\mathbf{y}}}\widetilde{\Phi}({\mathbf{y}}^{(k)},\mu_k;{\mathbf{x}}_{k_j})\) for all \(j=1,2,\ldots,M_k\). Specifically, we require \[\label{eq:gradinetbias2} \mathbb{E} \left[\|\mathcal{G}_{k_j}({\mathbf{y}}^{(k)},\mu_{k}) - \nabla_{{\mathbf{y}}} \widetilde{{\Phi}}( {\mathbf{y}}^{(k)},\mu_{k};{\mathbf{x}}_{k_j})\|^2\right]\leq \widehat{\epsilon}_k^2,\tag{13}\] for some scalar \(\widehat{\epsilon}_k\geq 0\) and all \(j=1,2,\ldots,M_k\). Specific constructions of these estimators are provided in Remarks [rem:exact]–[rem:exact2]. The numerical experiments presented in this paper adopt the techniques described in Remark [rem:exact]. The \(M_k\) gradient estimates are then aggregated into a stochastic gradient estimator \[\begin{align} \label{eq:gradinetg} \mathcal{G}({\mathbf{y}}^{(k)},\mu_{k}) = \frac{1}{M_k} \sum_{j=1}^{M_k}\mathcal{G}_{k_j}({\mathbf{y}}^{(k)},\mu_{k}) \end{align}\tag{14}\] of \(\mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}}[\nabla_{{\mathbf{y}}}\widetilde{\Phi}({\mathbf{y}}^{(k)},\mu_k;{\mathbf{x}})]\). We proceed with a proximal gradient update step, followed by updating the value of \(\mu\) according to a predefined nonincreasing rule. The SSPG method is summarized in Algorithm 1.

3.3 Convergence Results↩︎

The following lemma will be used for establishing the convergence results of Algorithm 1. Its proof is given in Appendix 9.6.

Now we are ready to present the main convergence rate result. The proof is given in Appendix 9.7.

There are multiple strategies for selecting \(\widehat{\epsilon}_k\) and \(\mu_k\). Typically, we set \(\widehat{\epsilon}_k\) in the same order as \(\epsilon\). For \(\mu_k\), it can be kept constant, decay over \(k\), or be updated based on the difference between the objective function values at two consecutive iterations [65]. Below, we present a corollary detailing specific choices for \(\widehat{\epsilon}_k\) and \(\mu_k\), along with an analysis of the computational complexity required by Algorithm 1. The proof is given in Appendix 9.8.

4 Numerical Experiments↩︎

We evaluate Algorithm 1 on four test problems related to WDRO and adversarially robustness: a newsvendor problem, a robust regression problem, an \(\infty\)-Wasserstein robust regression problem, and an adversarially robust image-classification problem. We compare SSPG with GDMax [66] and SDRO [43]. For the \(\infty\)-Wasserstein experiments, we additionally compare with FGSM [2] and PGD [4]. All experiments are implemented in Python; hardware information is summarized in Table 1. Our implementation of SDRO follows the publicly available implementation accompanying [43].

4.1 Newsvendor Problem↩︎

The newsvendor problem, which models the expected cost of a retailer under uncertain demand, takes the form of problem 4 . In this subsection, we consider solving problems 7 ,  9 and 10 using \[\label{eq:newsvendor} \ell(\theta, x) = v\theta - u\min(\theta,x)\quad\text{and}\quad d(x,z) = \frac{1}{2} (x - z)^{2},\tag{17}\] where \(\theta \in \mathbb{R}_+\) represents the inventory level, \(x \in \mathbb{R}\) denotes the demand, \(v = 5\) is the underage cost, and \(u = 7\) is the overage cost. To ensure that the inner maximization problem has a finite solution, we impose the constraint \(\lambda \geq 7\) on the target problems, as required in [67].

Problem parameters: We synthetically generate five different demand datasets, each consisting of \(n=100\) independent samples drawn from an exponential distribution with a rate parameter of 1. We set \(\delta = 1\) and \(p = 2\). For each dataset \(X_{\mathrm{train}}\), the empirical distribution \(\widehat{\mathbb{P}}_{n}\) is constructed from these \(n\) samples. The support set is defined as \(\mathcal{Z}= \{x:\min (X_{\mathrm{train}})\le x\le\max(X_{\mathrm{train}})\}\).

Algorithm parameters: For all methods, we use the same initialization \(\theta^{(0)}\sim \mathcal{U}(0, 1)\) on all five datasets, where \({\mathcal{U}}\) denotes the uniform distribution. In addition, for GDMax and SSPG, we use the same \(\lambda^{(0)} \sim \mathcal{U}(7, 15)\) across all five datasets. For SSPG, we set \(\mu_0 = \lambda^{(0)} \eta\), where \(\eta \sim \mathcal{U}(0.1, 1)\). Following [43], we use a grid search for SDRO to fine-tune the hyperparameters \(\lambda\) and \(\eta\) from the sets {7, 10, 15} and {0.1, 0.5, 1}, respectively. Each method is terminated after 1000 iterations.

Implementation of the compared methods at the \(k\)th iteration: For each method, we solve several inner maximization problems to obtain \(\{z_i^{(k+1)}\}_{i=1}^n\subset \mathbb{R}\). Specifically, for each \(i\in [n]\), we solve problem \(\max_{z \in {\mathcal{Z}}} \ell(\theta^{(k)}, z) - \lambda^{(k)} d(x_i, z)\) using the projected gradient ascent with fixed step size of \(10^{-2}\) for 20 iterations, starting from \(\mathrm{Proj}_{{\mathcal{Z}}}(x_i + 10^{-3}r)\) where \(r\) satisfies the normal distribution. After obtaining \(\{z_i^{(k+1)}\}_{i=1}^n\), GDMax performs projected gradient descent on \((\theta,\lambda)\) using \(\nabla_{(\theta,\lambda)}\frac{1}{n}\sum_{i=1}^n\left[\lambda^{(k)}\delta^p+\ell(\theta^{(k)},z_i^{(k+1)})-\lambda^{(k)} d(x_i,z_i^{(k+1)})\right]\) as a gradient and a fixed learning rate \(0.1\).

For SSPG and SDRO, for all \(i\in [n]\), we generate sets \(\Omega_i^k\) of size \(M=32\), containing samples near \(z_i^{(k+1)}\). For any \(\widehat{z}\in\Omega_i^k\), we set \(\widehat{z} = \mathrm{Proj}_{{\mathcal{Z}}}(z_i^{(k+1)} + r^{(k+1)})\) with \(r^{(k+1)}\) satisfying the normal distribution. SDRO then performs a projected gradient descent step on \(\theta\) with gradient \(\nabla_{\theta}{g}_s^k(\theta^{(k)},\lambda)\) and a fixed learning rate \(0.1\), where \[\begin{align} \label{eq:gradientgsk} {g}_s^k( \theta,\lambda) = \lambda \delta^{p} + \lambda\eta \frac{1}{n}\sum_{i=1}^n\left[ \log \left(\frac{1}{M} \sum_{\widehat{z}\in\Omega_i^k} \left[ e^{\frac{\ell(\theta, \widehat{z}) - \lambda d(x_i, \widehat{z})}{\lambda\eta}} \right] \right)\right]. \end{align}\tag{18}\] Our SSPG method conducts a projected gradient descent step on the primal variable \((\theta,\lambda)\) with \(\nabla_{(\theta,\lambda)}\widetilde{g}^k( \theta^{(k)},\lambda^{(k)} , \mu_k)\) and a fixed learning rate \(0.1\), where \[\begin{align} \label{eq:gradientgtildek} \widetilde{g}^k(\theta,\lambda,\mu) = \lambda \delta^{p} + \mu \frac{1}{n}\sum_{i=1}^n\left[ \log \left(\frac{1}{M} \sum_{\widehat{z}\in\Omega_i^k} \left[ e^{\frac{\ell(\theta, \widehat{z}) - \lambda d(x_i, \widehat{z})}{\mu}} \right] \right)\right]. \end{align}\tag{19}\] We then update \(\mu\) according to \[\label{eq:munews} \mu_{k+1} = \left\{ \begin{align} & \mu_k,\text{ if } \widetilde{g}^k( \theta^{(k+1)},\lambda^{(k+1)}, \mu_k) - \widetilde{g}^k( \theta^{(k)}, \lambda^{(k)},\mu_k)< -\mu_k^{2\sigma_2}, \\ &\max\{\sigma_1 \mu_k, 10^{-4}\mu_{0}\}, \,\,\, \text{otherwise}, \end{align} \right.\tag{20}\] where \(\sigma_{1} = 0.99\), \(\sigma_{2} = 0.5\).

Performance comparisons:

Figure 2 (a) reports the mean and standard deviation of \(g(\theta,\lambda)\) over the five datasets. In this experiment, SSPG attains lower final objective values than GDMax and SDRO. Figure 2 (b) reports the corresponding values of \(\lambda\). SSPG and GDMax move rapidly to the imposed lower bound \(\lambda=7\), while SDRO uses the fixed value \(\lambda=7\) selected by grid search.

Figure 2 (b) shows the mean and standard deviation of \(\lambda\). It reveals that SSPG and GDMax quickly converge to values of 7, while SDRO maintains a fixed \(\lambda = 7\). This choice results from a grid search indicating that the pair \(\lambda = 7, \eta=0.1\) yields the best performance among the tested configurations for SDRO.

4.2 Regression Problem↩︎

The distributionally robust regression problem aims to find a robust solution to the standard regression problem by minimizing the worst-case risk. In this subsection, we consider problems 7 , 9 , and 10 with \[\label{eq:ell-d-func} \ell({\boldsymbol{\theta}},{\mathbf{x}}) = (h_{{\boldsymbol{\theta}}}({\mathbf{a}}) - b)^{2}, \text{ and }d({\mathbf{x}},{\mathbf{z}}) = d(({\mathbf{a}}, b), (\overline{{\mathbf{a}}}, \overline{b})) = \frac{1}{2} \|{\mathbf{a}}- \overline{{\mathbf{a}}}\|_{2}^{2} + \infty |b - \overline{b}|,\tag{21}\] where \(h_{{\boldsymbol{\theta}}}: \mathbb{R}^{m_2-1} \to \mathbb{R}\) is a neural network parameterized by \({\boldsymbol{\theta}}\), \({\mathbf{x}}= ({\mathbf{a}}, b)\), \({\mathbf{z}}= (\overline{{\mathbf{a}}}, \overline{b})\) with \({\mathbf{a}}, \overline{{\mathbf{a}}} \in \mathbb{R}^{m_2-1}\) representing a vector of features and \(b, \overline{b} \in \mathbb{R}\) denoting a target. Specifically, we employ a neural network with a single hidden layer containing three neurons, using the ReLU activation function [68]. The \(\infty\) before \(|b - \overline{b}|\) means that there is no uncertainty in the target variable.

Datasets: We consider three real-world datasets, Space GA, BodyFat, and MG, from the LIBSVM repository. Each dataset, containing \(n\) data points, is randomly partitioned into training (80%) and testing (20%) sets using train_test_split from Scikit-learn, using the same random seeds across all methods. The resulting training and test sets are \(X_{\text{train}} \in \mathbb{R}^{\lfloor 0.8 \times n \rfloor \times m_2}\) and \(X_{\text{test}} \in \mathbb{R}^{\lceil 0.2 \times n \rceil \times m_2}\), respectively. These two sets are further normalized using StandardScaler. To ensure fair comparisons, we employ five distinct random seeds for data preparation, including splitting, and initialization, ensuring identical dataset configurations across all three methods.

Problem parameters: We set \(\delta = 10\) and \(p = 2\). For each dataset \(X_{\mathrm{train}}\), the empirical distribution \(\widehat{\mathbb{P}}_{n}\) is constructed from the training samples; the support set is given as \(\mathcal{Z} = \widetilde{\mathcal{Z}} \times \mathbb{R}\) with \(\widetilde{\mathcal{Z}}~=~\prod_{j=1}^{m_2-1} [ \min([X_{\text{train}}]_{\cdot,j}), \max([X_{\text{train}}]_{\cdot,j}) ]\).

Algorithm parameters: For all methods, we initialize \({\boldsymbol{\theta}}^{(0)}\) using PyTorch’s default parameter initialization under a fixed random seed. In addition, for GDMax and SSPG, we let \(\lambda^{(0)} \sim \mathcal{U}(1, 10)\). For SSPG, we set \(\mu_0 = \lambda^{(0)} \eta\), where \(\eta \sim \mathcal{U}(0.1, 1)\). Following [43], we utilize a grid search for SDRO to fine-tune the hyperparameters \(\lambda\) and \(\eta\) from the sets {1, 5,10} and {0.1, 0.5, 1}, respectively. For each comparison method, the training is terminated after 500 epochs.

Implementation of the compared methods at the \(k\)th iteration: For each method, we perform an inner maximization step to obtain \(\{{\mathbf{z}}_i^{(k+1)}= (\overline{{\mathbf{a}}}_i^{(k+1)}, \overline{b}_i^{(k+1)})\}_{i=1}^{\lfloor 0.8 \times n \rfloor} \subset \mathbb{R}^{m_2}\). Specifically, for each \(i\in[\lfloor 0.8 \times n \rfloor]\), we solve problem \(\max_{{\mathbf{z}}\in {\mathcal{Z}}} \ell({\boldsymbol{\theta}}^{(k)}, {\mathbf{z}}) - \lambda^{(k)} d({\mathbf{x}}_i, {\mathbf{z}})\) using projected gradient ascent with a fixed step size of \(10^{-2}\) for five iterations, starting from \({\mathbf{x}}_i + 10^{-3} {\mathbf{r}}\) where \({\mathbf{r}}\) satisfies the normal distribution. Notice that by the definition of \(d\) in 21 , we actually fix the last component of \({\mathbf{z}}_i\) as the label corresponding to \({\mathbf{x}}_i\).

After obtaining \(\{{\mathbf{z}}_i^{(k+1)}\}_{i=1}^{\lfloor 0.8 \times n \rfloor}\), GDMax updates \(({\boldsymbol{\theta}},\lambda)\) via a projected gradient descent step using the gradient \(\nabla_{({\boldsymbol{\theta}},\lambda)}\mathbb{E}_{{\mathbf{x}}\sim \widehat{\mathbb{P}}_n}\left[\lambda^{(k)}\delta^p+\ell({\boldsymbol{\theta}}^{(k)},{\mathbf{z}}^{(k+1)})-\lambda^{(k)} d({\mathbf{x}},{\mathbf{z}}^{(k+1)})\right]\) and a fixed learning rate \(\alpha>0\). For SSPG and SDRO, we generate, for each \(i\in[\lfloor 0.8 \times n \rfloor]\), a set \(\Omega_i^k\) of size \(M=32\), containing samples near \({\mathbf{z}}_i^{(k+1)}\). Specifically, for any \(\widehat{{\mathbf{z}}}\in\Omega_i^k\), we let \(\widehat{{\mathbf{z}}} = \mathrm{Proj}_{{\mathcal{Z}}}({\mathbf{z}}_i^{(k+1)} + {\mathbf{r}}^{(k+1)}),\) where \({\mathbf{r}}^{(k+1)})\) satisfies the normal distribution. SDRO then performs a projected gradient descent step on \({\boldsymbol{\theta}}\) with gradient \(\nabla_{{\boldsymbol{\theta}}}{g}_s^k({\boldsymbol{\theta}}^{(k)},\lambda)\) and a fixed learning rate \(\alpha\), where \({g}_s^k\) is given in 18 with \(\widehat{z}\) replaced by \(\widehat{{\mathbf{z}}}\). Our SSPG method similarly applies a projected gradient descent step on \(({\boldsymbol{\theta}},\lambda)\) with gradient \(\nabla_{({\boldsymbol{\theta}},\lambda)}\widetilde{g}^k( {\boldsymbol{\theta}}^{(k)},\lambda^{(k)} , \mu_k)\) and a fixed learning rate \(\alpha\), where \(\widetilde{g}^k\) is given in 19 with \(\widehat{z}\) replaced by \(\widehat{{\mathbf{z}}}\). We then update \(\mu\) by \[\label{eq:mu2} \begin{align} &\mu_{k+1} := \max\left\{10^{-4}\mu_{0},{\left(k+2\right)}^{-\frac{1}{3}}\mu_{0}\right\}. \end{align}\tag{22}\] For all compared methods, we perform a grid search to select the learning rate \(\alpha\) from the set \(\{10^{-1}, 5 \times 10^{-1}, 10^{-2}, 5 \times 10^{-2}, 10^{-3}\}\).

Performance comparisons: We measure model performance using the root mean square error (RMSE). Each model is evaluated on a modified version of the test set. Specifically, for each data point \({\mathbf{x}}=({\mathbf{a}}, b)\) in the test set, the feature vector \({\mathbf{a}}\) is perturbed according to \({\mathbf{a}}+\upsilon \boldsymbol{\omega}\|{\mathbf{a}}\|_2\), where \(\upsilon = 2\) and \(\boldsymbol{\omega} \sim [\text{Laplace}(0,1)]^{q}\) [68], where Laplace represents the Laplace distribution. For each candidate learning rate, we run 5 random seeds and determine the best learning rate \(\alpha\) for SSPG and GDMax based on the average RMSE. Additionally, for SDRO, we select the best-performing \(\alpha\) using the same criterion, while the best \(\lambda\) and \(\eta\) are determined per seed. Finally, we report the mean value and standard deviation of RMSE and training time across compared methods, as shown in Table 2.

From this table, we observe that SSPG achieves the lowest RMSE on two of the three datasets while also exhibiting smaller standard deviations, indicating superior error minimization and robust performance. SDRO, although it achieves the best performance on the MG dataset, requires significantly longer runtime, limiting its practical advantage. Overall, SSPG is an effective method for minimizing prediction errors and ensuring consistent performance, rendering it a good choice for robust regression tasks across diverse data environments.

4.3 Adversarially Robust Learning↩︎

We next consider the \(\infty\)-Wasserstein robust regression formulation discussed in Section 1.1. The problem can be expressed as \[\begin{align} \min_{{\boldsymbol{\theta}}} \frac{1}{n} \sum_{i=1}^{n} \max_{{\mathbf{z}}_i \in \mathcal{B}^{\infty}_{\hat{\epsilon}}({\mathbf{a}}_i)} \ell({\boldsymbol{\theta}}, ({\mathbf{z}}_i, b_{i}))~\text{ with }~ \ell({\boldsymbol{\theta}}, ({\mathbf{z}}_i, b_{i})):= \left(h_{{\boldsymbol{\theta}}}({\mathbf{z}}_i) - b_{i}\right)^{2} \label{eq:infinity95dro95objective} \end{align}\tag{23}\] where \(h_{{\boldsymbol{\theta}}}: \mathbb{R}^{m_2-1} \to \mathbb{R}\) is a neural network function parameterized by \({\boldsymbol{\theta}}\), \({\mathbf{a}}_{i} \in \mathbb{R}^{m_2-1}\) denotes a vector of features, \(b_{i} \in \mathbb{R}\) denotes the corresponding target, \(\mathcal{B}^{\infty}_{\hat{\epsilon}}({\mathbf{a}}_i)\) denotes the \(\infty\)-ball centered at \({\mathbf{a}}_i\) with radius \(\hat{\epsilon}\). The neural network architecture that we employ is a single hidden layer containing five neurons, using the ReLU activation function [68].

Datasets: The datasets are the same as those utilized in Section 4.2.

Algorithm parameters: For all methods, we initialize \({\boldsymbol{\theta}}^{(0)}\) using PyTorch’s default parameter initialization under a fixed random seed, and set \(\hat{\epsilon} = 0.25\). For SSPG, we set \(\mu_0 = \lambda^{(0)} \eta\), with \(\eta \sim \mathcal{U}(0.1, 1)\). The step size for the inner problem is \(\hat{\epsilon}\) for FGSM, while for GDMax and PGD, we use \(\frac{\hat{\epsilon}}{2}\). Following [43], we let \(\mu = \lambda \eta\) and perform a grid search for SDRO to fine-tune the hyperparameters \(\lambda\) and \(\eta\) from the sets {1, 5, 10} and {0.1, 0.5, 1} respectively. Additionally, for each method, we perform a grid search to select the optimal learning rate \(\alpha\), for the outer minimization problem, from the set \(\{10^{-1}, 5 \times 10^{-2}, 10^{-2}\},\) and the training is terminated after 250 epochs.

Implementation of the Compared Methods at the \(k\)th Iteration: For the inner maximization problem in 23 , each method first computes \(\{{\mathbf{z}}_i^{(k+1)}\}_{i=1}^{n} \subset \mathbb{R}^{m_2-1}\). Specifically, for each method, we initialize from points obtained by adding random noise uniformly sampled from the interval \([-\hat{\epsilon}, \hat{\epsilon}]\) to \({\mathbf{a}}_i\). For FGSM, the inner point is generated by one projected gradient step with stepsize \(\hat{\epsilon}\). For PGD, GDMax, SSPG, and SDRO, the inner maximization is approximated by 20 projected gradient-ascent steps with stepsize \(\hat{\epsilon}/2\), initialized from a uniformly perturbed point in \(\mathcal{B}^\infty_{\hat{\epsilon}}({\mathbf{a}}_i)\).

After obtaining \(\{{\mathbf{z}}_i^{(k+1)}\}_{i=1}^{n}\), GDMax updates \({\boldsymbol{\theta}}\) via a projected gradient descent step using the gradient \(\nabla_{{\boldsymbol{\theta}}}\mathbb{E}_{{\mathbf{x}}\sim \widehat{\mathbb{P}}_n}\left[\ell({\boldsymbol{\theta}}, ({\mathbf{z}}^{(k+1)}_i, b_{i}))\right]\) and a fixed learning rate \(\alpha>0\). For SSPG and SDRO, we generate, for each \(i\in[n]\), a set \(\Omega_i^k\) of size \(M=8\), containing samples near \({\mathbf{z}}_i^{(k+1)}\). Specifically, for any \(\widehat{{\mathbf{z}}}\in\Omega_i^k\), we let \(\widehat{{\mathbf{z}}} = \mathrm{Proj}_{\mathcal{B}^{\infty}_{\hat{\epsilon}}({\mathbf{a}}_i)}({\mathbf{z}}_i^{(k+1)} + {\mathbf{r}}^{(k+1)}),\) where \({\mathbf{r}}^{(k+1)})\) satisfies the normal distribution. SDRO then performs a projected gradient descent step on \({\boldsymbol{\theta}}\) with gradient \(\nabla_{{\boldsymbol{\theta}}}{g}_s^k({\boldsymbol{\theta}}^{(k)},\lambda)\) and a fixed learning rate \(\alpha\), where \({g}_s^k\) is given in 18 with \(\delta=0\), \(d\equiv 0\), and \(\ell(\theta,\hat{z})\) replaced by \(\ell({\boldsymbol{\theta}}, (\widehat{{\mathbf{z}}}, b_{i}))\). Our SSPG method similarly applies a projected gradient descent step on \({\boldsymbol{\theta}}\) with gradient \(\nabla_{{\boldsymbol{\theta}}}\widetilde{g}^k( {\boldsymbol{\theta}}^{(k)},0, \mu_k)\) and a fixed learning rate \(\alpha\), where \(\widetilde{g}^k\) is given in 19 with \(\ell(\theta,\hat{z})\) replaced by \(\ell({\boldsymbol{\theta}}, (\widehat{{\mathbf{z}}}, b_{i}))\). We then update \(\mu\) by 22 .

Performance comparisons: We conducted experiments on the \(\infty\)-Wasserstein DRO problem following the same experimental settings as those outlined for the numerical tests in Table 2. Then, we report the mean values and standard deviations of the MSE and training time for all compared methods in Table 3.

Table 3 shows that PGD attains the lowest mean MSE on all three datasets, while SSPG is close to PGD and substantially faster than SDRO. This is expected because PGD is tailored to the \(\ell_\infty\) perturbation model in 23 , whereas SSPG is designed for the broader expectation-over-maximization formulation.

4.4 Adversarially Robust Deep Learning Problem↩︎

We finally consider adversarially robust image classification, where the goal is to train a classifier that is robust to perturbations of the input images [4]. Given an image, we use a neural network prediction function \(h_{{\boldsymbol{\theta}}}: \mathbb{R}^{l \times w \times h} \to \mathbb{R}^{m_3}\), parameterized by \({\boldsymbol{\theta}}\), to predict the target class. The corresponding loss function and distance function in 7 , 9 , and 10 are given by \[\begin{align} \ell({\boldsymbol{\theta}},{\mathbf{x}}) = -\sum_{i=1}^{m_3} b_i \log \left( \frac{e^{\left[h_{{\boldsymbol{\theta}}}\left(\boldsymbol{a}\right)\right]_i}}{\sum_{j=1}^{m_3}e^{[h_{{\boldsymbol{\theta}}}(\boldsymbol{a})]_j}}\right), \text{ and } d(\mathbf{x}, {\mathbf{z}}) = \frac{1}{2} \|\mathbf{a} - \overline{{\mathbf{a}}}\|_{F}^{2} + \infty \|{\mathbf{b}}- \overline{{\mathbf{b}}}\|_1, \end{align}\] where \({\mathbf{x}}= ({\mathbf{a}}, {\mathbf{b}})\), \({\mathbf{z}}= (\overline{{\mathbf{a}}}, \overline{{\mathbf{b}}})\) with \({\mathbf{a}}, \overline{{\mathbf{a}}} \in {\mathbb{R}}^{l \times w \times h}\) representing an image and \({\mathbf{b}}, \overline{{\mathbf{b}}} \in \{0,1\}^{m_3}\) denoting their corresponding one-hot encoded label vectors. We set \(\delta = 10\), \(p = 2\), and \(\mathcal{Z} =[\min X_{\mathrm{train}},\max X_{\mathrm{train}}]^{l \times w \times h} \times \{0,1\}^{m_3}\).

Datasets and neural network architectures: To evaluate the efficacy of the compared methods, we use two benchmark datasets: Fashion-MNIST and CIFAR-10. These datasets are loaded using torchvision with the standard training/test split applied and the images are normalized using the usual mean-subtraction and standard-deviation scaling.

For Fashion-MNIST, we utilize a convolutional neural network (CNN) architecture consisting of three convolutional layers with 32, 64, and 128 filters, each using a \(3\times 3\) kernel, followed by ReLU activation and max pooling. The middle convolutional layer also employs dropout [68] with a probability of 0.3 and batch normalization [69]. The output of the convolutional layers is then passed through a fully connected layer with 512 hidden units, followed by ReLU activation, dropout (probability 0.25), and batch normalization.

For CIFAR-10, we adopt the All-CNN architecture [70], incorporating batch normalization after each ReLU activation in every convolutional layer.

Algorithm parameters: For all methods, we initialize \({\boldsymbol{\theta}}^{(0)}\) using PyTorch’s default parameter initialization under a fixed random seed. In addition, for GDMax and SSPG, we let \(\lambda^{(0)} \sim \mathcal{U}(1, 10)\). For SSPG, we set \(\mu_0 = \lambda^{(0)} \eta\), where \(\eta \sim \mathcal{U}(0.1, 1)\). We utilize a grid search for SDRO to fine-tune the hyperparameters \(\lambda\) and \(\eta\) from the sets {1, 10} and {0.1, 1}, respectively. For each comparison method, the training is terminated after 100 epochs.

Implementation of the compared methods at the \(k\)th iteration: For all methods, we first sample a mini-batch of size \(B=100\) from the training set, denoted as \(\{{\mathbf{x}}_{j_1^k},{\mathbf{x}}_{j_2^k},\ldots,{\mathbf{x}}_{j_B^k}\}\). Next, for each \(i~\in~[B]\), we perform an inner maximization step to obtain \(\{{\mathbf{z}}_{j_i}^{(k+1)}\}_{i=1}^B\). Specifically, we solve \(B\) problems of the form \(\max_{{\mathbf{z}}\in {\mathcal{Z}}} \ell({\boldsymbol{\theta}}^{(k)}, {\mathbf{z}}) - \lambda^{(k)} d({\mathbf{x}}_{j_i^k}, {\mathbf{z}})\) using projected gradient ascent with a fixed step size of \(10^{-2}\) for 15 iterations, starting from \({\mathbf{x}}_{j_i^k}+10^{-3}{\mathbf{r}}\) where \({\mathbf{r}}\) follows the standard normal distribution.

After obtaining \(\{{\mathbf{z}}_{j_i}^{(k+1)}\}_{i=1}^B\), GDMax updates the primal variable \(({\boldsymbol{\theta}},\lambda)\) via a projected gradient descent step with \(\nabla_{({\boldsymbol{\theta}},\lambda)} \frac{1}{B}\sum_{i=1}^B\left[\lambda^{(k)}\delta^p+\ell({\boldsymbol{\theta}}^{(k)},{\mathbf{z}}_{j_i^k}^{(k+1)})-\lambda^{(k)} d({\mathbf{x}}_{j_i^k},{\mathbf{z}}_{j_i^k}^{(k+1)})\right]\) and a learning rate \(\alpha_k>0\).

For SSPG and SDRO, we generate, for each \(i\in\{1,2,\ldots,B\}\), a set \(\Omega_{j_i}^k\) of size \(M=8\), containing samples near \({\mathbf{z}}_{j_i}^{(k+1)}\). Specifically, for any \(\widehat{{\mathbf{z}}}\in\Omega_{j_i}^k\), we form the immediate point \(\widehat{{\mathbf{z}}} = \mathrm{Proj}_{{\mathcal{Z}}}({\mathbf{z}}_{j_i}^{(k+1)} + {\mathbf{r}}^{(k+1)}),\) where \({\mathbf{r}}^{(k+1)})\) satisfies the normal distribution. For each \(j_i\), we retain only the samples that improve upon \({\mathbf{z}}_{j_i}^{(k+1)}\), defining the refined set as \[\begin{align} \overline{\Omega}_{j_i}^k:= \{{\mathbf{z}}_{j_i}^{(k+1)}\} \cup \left\{\widehat{{\mathbf{z}}}\in {\Omega}_{j_i}^k: \ell({\boldsymbol{\theta}}^{(k)}, \widehat{{\mathbf{z}}}) - \lambda^{(k)} d({\mathbf{x}}_{j_i^k}, \widehat{{\mathbf{z}}}) > \ell({\boldsymbol{\theta}}^{(k)}, {\mathbf{z}}_{j_i}^{(k+1)}) - \lambda^{(k)} d({\mathbf{x}}_{j_i^k}, {\mathbf{z}}_{j_i}^{(k+1)})\right\}. \end{align}\] SDRO then updates \({\boldsymbol{\theta}}\) via a projected gradient descent step with the gradient \(\nabla_{{\boldsymbol{\theta}}}{g}_s^k({\boldsymbol{\theta}}^{(k)},\lambda)\) and the learning rate \(\alpha_k\), where \[\begin{align} \label{eq:gradientgsk2} {g}_s^k( {\boldsymbol{\theta}},\lambda) = \lambda \delta^{p} + \lambda\eta \frac{1}{B}\sum_{i=1}^B\left[ \log \left(\frac{1}{|\overline{\Omega}_{j_i}^k|} \sum_{\widehat{{\mathbf{z}}}\in\overline{\Omega}_{j_i}^k} \left[ e^{\frac{\ell({\boldsymbol{\theta}}, \widehat{{\mathbf{z}}}) - \lambda d({\mathbf{x}}_{j_i^k}, \widehat{{\mathbf{z}}})}{\lambda\eta}} \right] \right)\right]. \end{align}\tag{24}\] Our SSPG method conducts a projected gradient descent step on the primal variable \(({\boldsymbol{\theta}},\lambda)\) using \(\nabla_{({\boldsymbol{\theta}},\lambda)}\widetilde{g}^k( {\boldsymbol{\theta}}^{(k)},\lambda^{(k)} , \mu_k)\) and the learning rate \(\alpha_k\), where \[\begin{align} \label{eq:gradientgtildek2} \widetilde{g}^k({\boldsymbol{\theta}},\lambda,\mu) = \lambda \delta^{p} + \mu \frac{1}{B}\sum_{i=1}^B\left[ \log \left(\frac{1}{|\overline{\Omega}_{j_i}^k|} \sum_{\widehat{{\mathbf{z}}}\in\overline{\Omega}_{j_i}^k} \left[ e^{\frac{\ell({\boldsymbol{\theta}}, \widehat{{\mathbf{z}}}) - \lambda d({\mathbf{x}}_{j_i}^k, \widehat{{\mathbf{z}}})}{\mu}} \right] \right)\right]. \end{align}\tag{25}\] We then update \(\mu\) by 22 . For all compared methods, we set \(\alpha_k= \alpha \gamma^{\lfloor k/20 \rfloor}\), where we perform a grid search to choose \(\alpha\) from the set \(\{1 \times 10^{-1}, 5 \times 10^{-1}, 5 \times 10^{-2}\}\), and \(\gamma\) is selected from \(\{0.5,0.9\}\).

Performance comparisons: We evaluate model performance using accuracy and overall training time. Each model is tested on a modified version of the test set, where each feature vector is perturbed in a manner similar to the distributionally robust regression setting. Specifically, for each data point \({\mathbf{x}}=({\mathbf{a}}, {\mathbf{b}})\) in the test set, we apply the perturbation \({\mathbf{a}}+\upsilon \boldsymbol{\omega}\|{\mathbf{a}}\|_2\), where \(\boldsymbol{\omega} \sim [\text{Laplace}(0,1)]^{q}\) [68], \(\upsilon = 2 \times 10^{-3}\) for Fashion-MNIST and \(\upsilon= 2 \times 10^{-4}\) for CIFAR-10.

To ensure a fair comparison, we use five distinct random seeds for initialization across all three methods. Given the large dataset sizes and high computational cost, we first sample 20% of the training set, ensuring that every class has the same number of samples, for hyperparameter tuning, optimizing \(\alpha\), \(\gamma\), \(\lambda\), and \(\eta\) from the specified choices in a manner similar to the one used for the regression problems. The best learning rate for SSPG and GDMax, as well as the hyperparameters for SDRO, are selected based on the highest accuracy. After selecting the parameters, we apply each method to the full training set. Finally, we aggregate the results across the five seeds and report in Table 4 the mean accuracy and training time, along with their standard deviations, for all methods.

Table 4 shows that SSPG attains the highest mean accuracy on Fashion-MNIST among the tested methods, while SDRO attains a slightly higher mean accuracy on CIFAR-10. SSPG is substantially faster than SDRO and has accuracy comparable to the best-performing method in both datasets. Overall, SSPG achieves an effective balance between accuracy and computational efficiency.

5 Conclusion↩︎

We study nonconvex minEmax problems in which the objective is an expectation of pointwise maxima. We introduce an LME smoothing of the random value function and developed SSPG, a stochastic smoothing proximal-gradient method for the resulting sequence of smoothed problems.

Our analysis provides two main guarantees. First, stationarity of the smoothed problems can be converted into Goldstein stationarity of the original nonsmooth problem through explicit smoothing-gap bounds. Second, SSPG attains an \(\epsilon\)-scaled stationary point in expectation with iteration complexity \(O(\epsilon^{-3})\) and sample complexity \(O(\epsilon^{-5})\), and it attains an \((O(\epsilon),O(\epsilon))\)-Goldstein stationary point with complexities \(\widetilde{O}(\epsilon^{-4})\) and \(\widetilde{O}(\epsilon^{-6})\), respectively. Under a refinement schedule with summable stationarity tolerances, any almost-sure cluster point of the generated approximate-stationary sequence is Clarke stationary, and hence directional stationary, for the original problem.

For WDRO, the derived dual-multiplier bound permits the dualized worst-case-risk problem to be treated on a compact feasible set. The numerical results on newsvendor, robust regression, \(\infty\)-Wasserstein/adversarially learning, and adversarially robust image classification demonstrate that SSPG is a practical first-order approach and is competitive with the tested baselines.

6 An Example Illustrating the Limitation of Value-Based Inner Accuracy↩︎

When solving \(\min_{{\mathbf{y}}}g({\mathbf{y}})\), subgradient methods may exhibit zigzagging behavior and may fail to approach stationary points, especially when the selected subgradients are unstable, [71]. In the expectation-over-maximization setting, a natural heuristic is to approximate a subgradient of \(g\) by \(\nabla_{{\mathbf{y}}}\Psi({\mathbf{y}},{\mathbf{z}}^\epsilon)\), where \({\mathbf{z}}^\epsilon\) is a near-optimal solution of the inner maximization problem. The following example shows that inner value accuracy alone does not control the distance from this vector to the Clarke subdifferential of the outer objective.

7 Boundedness of the Dual Multiplier of the Worst-case Risk Problem in WDRO↩︎

In this section, we present a lemma establishing an explicit bound \(B_{\lambda}\), ensuring that every optimal solution \(\lambda^*\) of problem 6 satisfies \(\lambda^*\in[0,B_{\lambda}]\).

Define the objective \[J(\lambda):=\lambda\delta^p+\mathbb{E}_{{\mathbf{x}}\sim \mathbb{P}} [\phi(\lambda;{\mathbf{x}})],\quad\text{with}\quad \phi(\lambda;{\mathbf{x}}):=\max_{{\mathbf{z}}\in\mathcal{Z}}\{\ell({\boldsymbol{\theta}},{\mathbf{z}})-\lambda\|{\mathbf{z}}-{\mathbf{x}}\|_p^p\}.\] Since each \(\phi(\lambda;{\mathbf{x}})\) is a pointwise supremum of affine functions in \(\lambda\), it is convex and non-increasing, hence \(J(\lambda)\) is convex.

7.0.0.1 Case of \(p>1\).

By Danskin’s theorem [72], \[\begin{align} \partial\phi(\lambda;{\mathbf{x}}) &= -\,\operatorname{co}\Big\{\,\|{\mathbf{z}}-{\mathbf{x}}\|_p^p \;:\;{\mathbf{z}}\in\arg\max_{{\mathbf{u}}\in{\mathcal{Z}}}\big(\ell({\boldsymbol{\theta}},{\mathbf{u}})-\lambda\|{\mathbf{u}}-{\mathbf{x}}\|_p^p\big)\Big\}. \label{eq:partiallb} \end{align}\tag{26}\] Let \({\mathbf{z}}(\lambda;{\mathbf{x}})\in\arg\max_{{\mathbf{z}}\in{\mathcal{Z}}}\{\ell({\boldsymbol{\theta}},{\mathbf{z}})-\lambda\|{\mathbf{z}}-{\mathbf{x}}\|_p^p\}\), and define \(t(\lambda;{\mathbf{x}})=\|{\mathbf{z}}(\lambda;{\mathbf{x}})-{\mathbf{x}}\|_p\).

Fix \({\mathbf{x}}\in{\mathcal{X}}\). We now prove \[\begin{align} \label{eq:partialphi} \partial\phi(\lambda;{\mathbf{x}})\;\subseteq\;\Big[-\Big(\tfrac{L\,C_{p,m_2}}{\lambda}\Big)^{\frac{p}{p-1}},\;0\Big], \text{ for all }\lambda>0. \end{align}\tag{27}\] If every active maximizer satisfies \({\mathbf{z}}={\mathbf{x}}\), then all active slopes in 26 are zero, and hence \(\partial q(\lambda;{\mathbf{x}})=\{0\}\). Thus, we conclude that \(\partial\phi(\lambda;{\mathbf{x}}) = \{0\}\), which then implies 27 . Otherwise, there exists \(z(\lambda;{\mathbf{x}})\in{\mathcal{Z}}\) such that \({\mathbf{z}}(\lambda;{\mathbf{x}})\neq{\mathbf{x}}\). We then consider the case when \(t(\lambda;{\mathbf{x}})=\|{\mathbf{z}}(\lambda;{\mathbf{x}})-{\mathbf{x}}\|_p>0\). Due to the fact that \({\mathbf{x}}\in{\mathcal{Z}}\), it follows directly from optimality of \({\mathbf{z}}(\lambda;{\mathbf{x}})\) that \(\ell({\boldsymbol{\theta}},{\mathbf{z}}(\lambda;{\mathbf{x}}))-\lambda t(\lambda;{\mathbf{x}})^p \;\ge\; \ell({\boldsymbol{\theta}},{\mathbf{x}}).\) Using \(L\)-Lipschitz continuity of \(\ell({\boldsymbol{\theta}},\cdot)\) and \(\|{\mathbf{v}}\|_2\le C_{p,m_2}\|{\mathbf{v}}\|_p\) for all \({\mathbf{v}}\in\mathbb{R}^{m_2}\), we arrive at \[\label{eq:basic-ti} \lambda\,t(\lambda;{\mathbf{x}})^p \;\le\; \ell({\boldsymbol{\theta}},{\mathbf{z}}(\lambda;{\mathbf{x}}))-\ell({\boldsymbol{\theta}},{\mathbf{x}}) \;\le\; L\,\|{\mathbf{z}}(\lambda;{\mathbf{x}})-{\mathbf{x}}\|_2 \;\le\; L\,C_{p,m_2}\,t(\lambda;{\mathbf{x}}).\tag{28}\] For \(\lambda>0\), dividing 28 by \(t(\lambda;{\mathbf{x}})\) yields \(t(\lambda;{\mathbf{x}})^{p-1}\;\le\;\frac{L\,C_{p,m_2}}{\lambda}\), and hence \(-t(\lambda;{\mathbf{x}})^p \ge- \Big(\frac{L\,C_{p,m_2}}{\lambda}\Big)^{\frac{p}{p-1}}.\) Since this bound holds for every maximizer, we have 27 .

By the optimality condition of problem 6 , any minimizer \(\lambda^*\) satisfies \[0\in \partial J(\lambda^*)+\mathcal{N}_{[0,\infty)}(\lambda^*).\] If \(\lambda^*>0\), then \(0\in\partial J(\lambda^*)=\delta^p+\partial \mathbb{E}_{{\mathbf{x}}\sim \mathbb{P}}[\phi(\lambda^*;{\mathbf{x}})]\subseteq \delta^p+\mathbb{E}_{{\mathbf{x}}\sim \mathbb{P}}[\partial\phi(\lambda^*;{\mathbf{x}})]\), where the second inclusion comes from [55]. Therefore, there exists a measurable selection \(s(\lambda^*;{\mathbf{x}})\in\partial\phi(\lambda^*;{\mathbf{x}})\) with \[0=\delta^p+\mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}}[s(\lambda^*;{\mathbf{x}})] \;\ge\;\delta^p-\Big(\tfrac{L\,C_{p,m_2}}{\lambda^*}\Big)^{\frac{p}{p-1}}.\] Hence \(\lambda^*\le L\,C_{p,m_2}\,\delta^{-(p-1)}\). If \(\lambda^*=0\), the bound is trivial.

7.0.0.2 Case of \(p=1\).

For any \({\mathbf{z}}\in{\mathcal{Z}}\), \[\ell({\boldsymbol{\theta}},{\mathbf{z}})-\lambda\|{\mathbf{z}}-{\mathbf{x}}\|_1 \le \ell({\boldsymbol{\theta}},{\mathbf{x}})+L\|{\mathbf{z}}-{\mathbf{x}}\|_2-\lambda\|{\mathbf{z}}-{\mathbf{x}}\|_1 \le \ell({\boldsymbol{\theta}},{\mathbf{x}})+(L-\lambda)\|{\mathbf{z}}-{\mathbf{x}}\|_1.\] Thus \(\phi(\lambda;{\mathbf{x}})\le \ell({\boldsymbol{\theta}},{\mathbf{x}})\) for all \(\lambda\ge L\). Since \({\mathbf{z}}={\mathbf{x}}\in{\mathcal{Z}}\), it holds \(\phi(\lambda;{\mathbf{x}})= \ell({\boldsymbol{\theta}},{\mathbf{x}})\) for \(\lambda\ge L\). Therefore \(J(\lambda)=\lambda\,\delta+\mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}}[\ell({\boldsymbol{\theta}},{\mathbf{x}})]\) is strictly increasing on \([L,\infty)\) (because \(\delta>0\)), so any minimizer satisfies \(\lambda^*\le L\) (equivalently, \(\lambda^*\le L\,C_{1,m_2}\,\delta^{0}\) with \(C_{1,m_2}=1\)).

The proof is then completed by combining the above two cases.

8 Methods to Generate a Stochastic Gradient Estimator↩︎

In this section, we detail two situations under which the gradient estimation condition 13 can be satisfied. As noted in Remark [rem:exact], if, for all \({\mathbf{x}}\in{\mathcal{X}}\), the function \(\Psi({\mathbf{y}}, \cdot; {\mathbf{x}})\) exhibits certain structures, we can efficiently generate a stochastic gradient estimator. For simplicity of notations, we fix \({\mathbf{x}}\in{\mathcal{X}}\) in this section, and define \(H({\mathbf{z}}):=\Psi({\mathbf{y}}, {\mathbf{z}}; {\mathbf{x}})\).

8.1 Computing a Stochastic Gradient Estimator with a Specific Loss Function and Support Set↩︎

As noted in Remark [rem:exact], if we can compute \(\mathbb{E}_{{\mathbf{z}}\sim \zeta}[e^{H({\mathbf{z}})/\mu}]\) for all \(\mu > 0\), we can generate a desired stochastic gradient estimator. We now describe several cases, in which this expectation can be computed. When \(H(\cdot)\) is a linear function and \(\zeta\) is the uniform distribution over its support set, we can calculate \(\mathbb{E}_{{\mathbf{z}}\sim \zeta}[e^{H({\mathbf{z}})/\mu}]\). This computation often reduces to evaluating the integrals of \(e^{{\mathbf{a}}^\top {\mathbf{z}}+ c}\) over \(\mathcal{Z}\). Notably, [73] demonstrates that for several structured sets \(\mathcal{Z}\), such exponential integrals admit efficient finite-dimensional formulas or algorithms. Specifically, the authors provide efficient methods for the following four cases:

1. \(\mathcal{Z}\) is a regular simplex, i.e., \(\mathcal{Z}=\left\{{\mathbf{z}}\in \mathbb{R}_{+}^{m_2}: \mathbf{1}_{m_2}^{\top}{\mathbf{z}}=1\right\}\);

2. \(\mathcal{Z} = [0,r]^{m_2}\) is a \(m_2\)-dimensional cube;

3. \(\mathcal{Z}=\operatorname{conic}\left\{{\mathbf{u}}_1, \ldots, {\mathbf{u}}_s\right\} \subset \mathbb{R}^{m_2}\) is a simple convex cone given as the conic hull of its extreme rays \({\mathbf{u}}_1, \ldots, {\mathbf{u}}_s \in \mathbb{R}^{m_2}\);

4. \(\mathcal{Z}\) is an intersection of several half-spaces.

8.2 A Sampling Method to Generate a Stochastic Gradient Estimator with a Connected Compact set \({\mathcal{Z}}\)↩︎

In this section, we consider the case where \({\mathcal{Z}}\) is a connected compact set, and introduce an algorithmic framework that outputs an approximate gradient \(\mathcal{G}_{k_j}({\mathbf{y}}^{(k)},\mu_k)\) satisfying condition 13 , i.e., \[\label{eq:gradinetbias2-restate} \mathbb{E}\!\big[\|\mathcal{G}_{k_j}({\mathbf{y}}^{(k)},\mu_k) -\nabla_{{\mathbf{y}}}\widetilde{\Phi}({\mathbf{y}}^{(k)},\mu_k;{\mathbf{x}}_{k_j})\|^2\big] \le \widehat\epsilon_k^2.\tag{29}\] The construction is used at a fixed outer iteration \(k\) and a fixed sample \({\mathbf{x}}_{k_j}\). For simplicity, in the remaining part of this section, unless with further specification, we denote and fix \[\label{eq:fix-x-y-mu} ({\mathbf{x}},{\mathbf{y}},\mu)=({\mathbf{x}}_{k_j},{\mathbf{y}}^{(k)},\mu_k).\tag{30}\] We will build a localized estimator for \(\nabla_{{\mathbf{y}}}\widetilde{\Phi}({\mathbf{y}},\mu;{\mathbf{x}})\) and choose its sample size and localization radius so that the mean-square error is at most \(\widehat\epsilon_k^2\).

The estimator analyzed here is a localized version of the LME Gibbs estimator. We define \[\label{eq:app-gap-definition} \begin{align} &F({\mathbf{z}}) = F_{{\mathbf{x}},{\mathbf{y}}}({\mathbf{z}}) := \Phi({\mathbf{y}};{\mathbf{x}})-\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}}) = \max_{{\mathbf{u}}\in{\mathcal{Z}}}\Psi({\mathbf{y}},{\mathbf{u}};{\mathbf{x}})-\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}})\ge0, \\ & w_\mu({\mathbf{z}}):=\exp\!\left(-\mu^{-1}{F({\mathbf{z}})} \right). \end{align}\tag{31}\] Then we have \(\exp\!\left(\mu^{-1}{\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}})} \right) = \exp\!\left(\mu^{-1}{\Phi({\mathbf{y}};{\mathbf{x}})} \right) w_\mu({\mathbf{z}}).\) For a localization radius \(r>0\), define the near-optimal level set \[\label{eq:app-local-set} {\mathcal{U}}_r={\mathcal{U}}_r({\mathbf{x}},{\mathbf{y}}):= \big\{{\mathbf{z}}\in{\mathcal{Z}}:F_{{\mathbf{x}},{\mathbf{y}}}({\mathbf{z}})\le r\big\},\tag{32}\] and let \(\nu_{r}\) be the conditional uniform law on \({\mathcal{U}}_r\): \[\label{eq:app-conditional-law} \nu_{r}(A):= \frac{\zeta(A\cap{\mathcal{U}}_r)}{\zeta({\mathcal{U}}_r)}, \quad \text{for any measurable set }A\subseteq{\mathcal{Z}}.\tag{33}\] Moreover, we define \[D_r=\int_{{\mathcal{U}}_r}e^{-F({\mathbf{z}})/\mu}\,\zeta(d{\mathbf{z}}), \text{ and } D_r^{\rm out}=\int_{{\mathcal{Z}}\setminus{\mathcal{U}}_r}e^{-F({\mathbf{z}})/\mu}\,\zeta(d{\mathbf{z}}).\] Notice that \(D_r\) is the contribution of the local region \({\mathcal{U}}_r\) to the Gibbs partition function, and \(D_r^{\rm out}\) is the contribution from the complement \({\mathcal{Z}}\setminus {\mathcal{U}}_r\).

With the above notations, we present the method to generate a stochastic gradient estimator in Algorithm 3. Equivalently, the output of Algorithm 3 equals \[\label{eq:app-local-estimator} G_{{\mathbf{x}},r}({\mathbf{y}},\mu) = \frac{\sum_{j=1}^M \exp\!\left(-\mu^{-1}F({\mathbf{z}}_j)\right) \nabla_{{\mathbf{y}}}\Psi({\mathbf{y}},{\mathbf{z}}_j;{\mathbf{x}})}{\sum_{j=1}^M \exp\!\left(-\mu^{-1}F({\mathbf{z}}_j)\right)} = \frac{\sum_{j=1}^M w_\mu({\mathbf{z}}_j)\nabla_{{\mathbf{y}}}\Psi({\mathbf{y}},{\mathbf{z}}_j;{\mathbf{x}})}{\sum_{j=1}^M w_\mu({\mathbf{z}}_j)}.\tag{34}\]

To guarantee that \(\mathcal{G}_{k_j}({\mathbf{y}}^{(k)},\mu_k):=G_{{\mathbf{x}},r}({\mathbf{y}},\mu)\) with notation given in 30 satisfies the condition in 29 , we need a flat-maximum condition assumed below. It is a positive-measure near-optimal-level-set condition. Related positive-measure optimality regions appear in random search through the essential optimum [74]; near-optimality dimension in optimistic optimization and \(X\)-armed bandits quantifies the size of near-optimal sets [75], [76]; in continuous evolutionary optimization it is explicitly described as a level set of positive Lebesgue measure [77]. Further discussions are provided in Remarks [rem:app-eb-vs-lower-growth]–[rem:app-flat-plateau-and-isolated-maximizers].

The next theorem gives the bound on the difference between \(G_{{\mathbf{x}},r}({\mathbf{y}},\mu)\) and \(\nabla_{{\mathbf{y}}}\widetilde{\Phi}({\mathbf{y}},\mu;{\mathbf{x}})\) in expectation.

9 Proofs of the Main Results↩︎

In this section we provide proofs of our results presented in Section 3.

9.1 Proof of Lemma 2↩︎

(a) We take \(\mu\downarrow 0\) in \(\mu \log \mathbb{E}_{{\mathbf{z}}\sim \zeta} [e^{\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}})/\mu}]\) to obtain \[\begin{align} \notag & \lim _{\mu \downarrow 0} \mu \log \mathbb{E}_{{\mathbf{z}}\sim \zeta} [e^{\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}})/\mu} ]= \lim _{\beta \rightarrow \infty} \frac{1}{\beta} \log \mathbb{E}_{{\mathbf{z}}\sim \zeta}[e^{\beta\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}}) }] \stackrel{(i)}{=}\lim_{\beta \rightarrow \infty} \nabla_\beta \log \mathbb{E}_{{\mathbf{z}}\sim \zeta}[e^{\beta\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}}) }] \\ = & {\lim _{\beta \rightarrow \infty} \frac{ \nabla_\beta \mathbb{E}_{{\mathbf{z}}\sim \zeta}[e^{\beta\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}}) } ]}{\mathbb{E}_{{\mathbf{z}}\sim \zeta}[e^{\beta\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}}) }] } \stackrel{(ii)}{=} \lim _{\beta \rightarrow \infty} \frac{ \mathbb{E}_{{\mathbf{z}}\sim \zeta}[\nabla_\beta e^{\beta\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}}) } ]}{\mathbb{E}_{{\mathbf{z}}\sim \zeta}[e^{\beta\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}}) }] } } = \lim _{\beta \rightarrow \infty} \frac{ \mathbb{E}_{{\mathbf{z}}\sim \zeta}[e^{\beta\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}}) } \Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}})]}{\mathbb{E}_{{\mathbf{z}}\sim \zeta}[e^{\beta\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}}) }] } \label{eq:limitu} \\\notag = &\lim _{\beta \rightarrow \infty} \frac{ \mathbb{E}_{{\mathbf{z}}\sim \zeta}[e^{\beta\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}}) -\beta \max_{{\mathbf{z}}\in \mathrm{supp}(\zeta)}\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}}) } \Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}})]}{\mathbb{E}_{{\mathbf{z}}\sim \zeta}[e^{\beta\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}}) -\beta \max_{{\mathbf{z}}\in \mathrm{supp}(\zeta)}\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}}) }] } \stackrel{(iii)}{=} \max_{{\mathbf{z}}\in \mathrm{supp}(\zeta)}\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}}) \stackrel{(iv)}{=}\max _{{\mathbf{z}}\in \mathcal{Z}}\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}}). \end{align}\tag{46}\] Here (i) comes from L’Hôpital’s rule [79]. When \(\mathcal{Z}\) is a finite set, (ii) holds directly; and when \(\mathcal{Z}\) is a connected compact set, (ii) holds by Leibniz integral rule [80] and the continuity of \(e^{\beta\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}}) } \Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}})\) with respect to \(\beta\). (iii) results from the fact that \({\mathbf{z}}\) is a continuous random variable, \[\begin{align} &\lim _{\beta \rightarrow \infty} e^{\beta\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}}) -\beta \max_{{\mathbf{z}}\in \mathrm{supp}(\zeta)}\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}})} = \left\{ \begin{aligned} 1, &\text{ if } {\mathbf{z}}\in\arg\max_{{\mathbf{z}}\in \mathrm{supp}(\zeta)}\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}}) \\ 0, & \text{ otherwise}. \end{aligned} \right., \text{ and}\\ &\lim _{\beta \rightarrow \infty} e^{\beta\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}}) -\beta \max_{{\mathbf{z}}\in \mathrm{supp}(\zeta)}\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}})} \Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}}) = \left\{ \begin{align} \Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}}), &\text{ if } {\mathbf{z}}\in\arg\max_{{\mathbf{z}}\in \mathrm{supp}(\zeta)}\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}}) \\ 0, \quad\quad& \text{ otherwise}. \end{align} \right., \end{align}\] In addition, (iv) follows from the definition of \(\zeta\).

Notice that \(\mu \log \mathbb{E}_{{\mathbf{z}}\sim \zeta} [e^{\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}})/\mu} ]\) is continuously differentiable with respect to \({\mathbf{y}}\) by Assumption 1. We then obtain that \(\widetilde{\Phi}(\cdot,\cdot;{\mathbf{x}})\) is a smoothing function of \(\Phi(\cdot;{\mathbf{x}})\) from Definition 3.

(b) The \({\mathbf{y}}\) and \(\mu\) partial gradients of \(\widetilde{\Phi}(\cdot,\cdot;{\mathbf{x}})\) are derived by direct calculation. For the \({\mathbf{y}}\)-partial gradient, we have \[\begin{align} \left\|\nabla_{{\mathbf{y}}} \widetilde{\Phi}({\mathbf{y}},\mu;{\mathbf{x}})\right\| &= \left\|\frac{\mathbb{E}_{{\mathbf{z}}\sim \zeta}[e^{\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}}) / \mu} \nabla_{{\mathbf{y}}}\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}})] }{\mathbb{E}_{{\mathbf{z}}\sim \zeta}[e^{\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}}) / \mu}] }\right\| \leq \frac{\mathbb{E}_{{\mathbf{z}}\sim \zeta}\left[\left\|e^{\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}}) / \mu} \nabla_{{\mathbf{y}}}\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}})\right\| \right] }{\left\|\mathbb{E}_{{\mathbf{z}}\sim \zeta}[e^{\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}}) / \mu}]\right\| }\\ &\leq \frac{\mathbb{E}_{{\mathbf{z}}\sim \zeta}\left[l_{\Psi}\left\|e^{\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}}) / \mu} \right\| \right] }{\left\|\mathbb{E}_{{\mathbf{z}}\sim \zeta}[e^{\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}}) / \mu}]\right\| }= l_{\Psi}. \end{align}\] For the \(\mu\)-partial gradient, it holds \(\lim_{ \mu \downarrow 0} \mu \nabla_\mu \widetilde{\Phi}(\mathbf{y}, \mu;{\mathbf{x}}) =0\) by 46 .

Let \(\sigma = \frac{\mu_2}{\mu_1}\leq 1\). We then have \[\begin{align} &\mu_1 \log \left(\mathbb{E}_{{\mathbf{z}}\sim \zeta} [e^{\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}})/\mu_1} ] \right) - \mu_2 \log \left(\mathbb{E}_{{\mathbf{z}}\sim \zeta} [e^{\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}})/\mu_2} ]\right) \\ =& \mu_1 \log \left(\mathbb{E}_{{\mathbf{z}}\sim \zeta} [e^{\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}})/\mu_1} ]\right) - \sigma\mu_1 \log \left(\mathbb{E}_{{\mathbf{z}}\sim \zeta} [e^{\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}})/(\sigma\mu_1)} ]\right) \\= & \mu_1 \log \frac{\mathbb{E}_{{\mathbf{z}}\sim \zeta} [e^{\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}})/\mu_1} ]}{\left(\mathbb{E}_{{\mathbf{z}}\sim \zeta} [e^{\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}})/(\sigma\mu_1)} ]\right)^\sigma} \leq 0, \end{align}\] where the inequality follows from Jensen’s inequality using \(\mathbb{E}\left[X^{1/\sigma}\right] \geq (\mathbb{E}[X])^{1/\sigma}\) with \(X=e^{\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}})/\mu_1}\). Thus \(\widetilde{\Phi}({\mathbf{y}},\mu_1;{\mathbf{x}})\leq \widetilde{\Phi}({\mathbf{y}},\mu_2;{\mathbf{x}})\) for any \(\mu_1\geq\mu_2>0\).

(c) Fix \({\mathbf{x}}\in{\mathcal{X}}\), \(\mu>0\), and \({\mathbf{y}}_1,{\mathbf{y}}_2\in\operatorname{dom}(\varphi)\). For each \({\mathbf{y}}\in\operatorname{dom}(\varphi)\), define the probability measure \(Q_{\mathbf{y}}\) on \(\mathcal{Z}\) by \[\mathrm{d}Q_{\mathbf{y}}({\mathbf{z}}) = \frac{\exp(\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}})/\mu)}{\mathbb{E}_{{\mathbf{z}}\sim\zeta}[\exp(\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}})/\mu)]}\,\zeta(d{\mathbf{z}}).\] Then \(\nabla_{\mathbf{y}}\widetilde{\Phi}({\mathbf{y}},\mu;{\mathbf{x}}) = \mathbb{E}_{{\mathbf{z}}\sim Q_{\mathbf{y}}}[\nabla_{\mathbf{y}}\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}})].\) Let \({\mathbf{d}}:={\mathbf{y}}_1-{\mathbf{y}}_2\), \(g_i({\mathbf{z}}):=\nabla_{\mathbf{y}}\Psi({\mathbf{y}}_i,{\mathbf{z}};{\mathbf{x}})\) for \(i=1,2\), and \(Q_i:=Q_{{\mathbf{y}}_i}\). We decompose \[\begin{align} &\|\nabla_{\mathbf{y}}\widetilde{\Phi}({\mathbf{y}}_1,\mu;{\mathbf{x}})-\nabla_{\mathbf{y}}\widetilde{\Phi}({\mathbf{y}}_2,\mu;{\mathbf{x}})\| \le \|\mathbb{E}_{Q_1}[g_1({\mathbf{z}})-g_2({\mathbf{z}})]\| + \|\mathbb{E}_{Q_1}[g_2({\mathbf{z}})]-\mathbb{E}_{Q_2}[g_2({\mathbf{z}})]\|. \end{align}\] By Assumption 3, the first term is bounded by \[\begin{align} \label{eq:bound1} \|\mathbb{E}_{Q_1}[g_1({\mathbf{z}})-g_2({\mathbf{z}})]\| \le L_\Psi\|{\mathbf{y}}_1-{\mathbf{y}}_2\|. \end{align}\tag{47}\] It remains to bound the second term. Let \({\mathbf{y}}_t:={\mathbf{y}}_2+t{\mathbf{d}}\) for \(t\in[0,1]\), and let \(Q_t:=Q_{{\mathbf{y}}_t}\). Define \(H(t):=\mathbb{E}_{{\mathbf{z}}\sim Q_t}[g_2({\mathbf{z}})].\) Since \(\|\nabla_{\mathbf{y}}\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}})\|\le l_\Psi\) by Assumption 3, differentiation under the expectation is justified by dominated convergence. By the standard covariance identity for normalized exponential tilts [81], for all \(t\in[0,1]\), we have \[H'(t) = \frac{1}{\mu} \left( \mathbb{E}_{Q_t}\left[g_2({\mathbf{z}})\left\langle \nabla_{\mathbf{y}}\Psi({\mathbf{y}}_t,{\mathbf{z}};{\mathbf{x}}),{\mathbf{d}}\right\rangle\right] - \mathbb{E}_{Q_t}[g_2({\mathbf{z}})]\, \mathbb{E}_{Q_t}\left[\left\langle \nabla_{\mathbf{y}}\Psi({\mathbf{y}}_t,{\mathbf{z}};{\mathbf{x}}),{\mathbf{d}}\right\rangle\right] \right).\] For any unit vector \({\mathbf{v}}\), the Cauchy–Schwarz inequality together with Assumption 2 gives \[\begin{align} |{\mathbf{v}}^\top H'(t)| &= \frac{1}{\mu} \left| \operatorname{Cov}_{Q_t} \left( {\mathbf{v}}^\top g_2({\mathbf{z}}), \left\langle \nabla_{\mathbf{y}}\Psi({\mathbf{y}}_t,{\mathbf{z}};{\mathbf{x}}),{\mathbf{d}}\right\rangle \right) \right| \\ &\le \frac{1}{\mu} \sqrt{ \operatorname{Var}_{Q_t}({\mathbf{v}}^\top g_2({\mathbf{z}})) \operatorname{Var}_{Q_t} \left( \left\langle \nabla_{\mathbf{y}}\Psi({\mathbf{y}}_t,{\mathbf{z}};{\mathbf{x}}),{\mathbf{d}}\right\rangle \right) } \\ &\le \frac{1}{\mu} { \|{\mathbf{v}}\|\|g_2({\mathbf{z}})\| \| \nabla_{\mathbf{y}}\Psi({\mathbf{y}}_t,{\mathbf{z}};{\mathbf{x}})\| \|{\mathbf{d}}\|} \le \frac{l_\Psi^2}{\mu}\|{\mathbf{d}}\|. \end{align}\] Taking the supremum over all unit vectors \({\mathbf{v}}\), we obtain \(\|H'(t)\| \le \frac{l_\Psi^2}{\mu}\|{\mathbf{y}}_1-{\mathbf{y}}_2\|.\) Therefore, \[\begin{align} \label{eq:bound2} \|\mathbb{E}_{Q_1}[g_2({\mathbf{z}})]-\mathbb{E}_{Q_2}[g_2({\mathbf{z}})]\| = \|H(1)-H(0)\| \le \int_0^1 \|H'(t)\|\,dt \le \frac{l_\Psi^2}{\mu}\|{\mathbf{y}}_1-{\mathbf{y}}_2\|. \end{align}\tag{48}\] Combining 47 and 48 yields \[\|\nabla_{\mathbf{y}}\widetilde{\Phi}({\mathbf{y}}_1,\mu;{\mathbf{x}})-\nabla_{\mathbf{y}}\widetilde{\Phi}({\mathbf{y}}_2,\mu;{\mathbf{x}})\| \le \left(L_\Psi+\frac{l_\Psi^2}{\mu}\right)\|{\mathbf{y}}_1-{\mathbf{y}}_2\|.\] The proof is then completed.

9.2 Proof of Lemma 3↩︎

Define \(\omega_t:\mathbb{R}^t \mapsto\mathbb{R}\) and \(\widetilde{b}_t:\mathbb{R}^t\times \mathbb{R}\mapsto\mathbb{R}\) by \[\omega_t({\mathbf{x}}):=\log \left( \sum_{i=1}^t e^{x_i} \right), \text{ and }\widetilde{b}_t({\mathbf{y}},\mu):= \mu \log \left(\sum_{i=1}^t e^{y_i / \mu}\right).\] We have \(\widetilde{b}_t({\mathbf{y}},\mu) = \mu \omega_t\left(\frac{{\mathbf{y}}}{\mu}\right)= \max _{{\mathbf{x}}\in D_t}\left\{\langle {\mathbf{x}}, {\mathbf{y}}\rangle -\mu \omega_t^*({\mathbf{x}})\right\},\) where \(D_t:=\{\mathbf{x} \in \mathbb{R}^t: \mathbf{x} \geq \mathbf{0}, \boldsymbol{1}_t^{\top}\mathbf{x}=1\}\), and the second equality follows from that the conjugate function of \(\omega\) over \(D_t\) is \(\omega_t^*({\mathbf{y}})=\sum_{i=1}^t y_i \log y_i\) (cf. [82]). Then \[\begin{align} \notag & \left|\widetilde{b}_t({\mathbf{y}},\mu_1)-\widetilde{b}_t({\mathbf{y}},\mu_2)\right| = \left| \max _{{\mathbf{x}}\in D_t }\left\{\langle {\mathbf{x}}, {\mathbf{y}}\rangle -\mu_2 \omega_t^*({\mathbf{x}})\right\}-\max _{{\mathbf{x}}\in D_t }\left\{\langle {\mathbf{x}}, {\mathbf{y}}\rangle -\mu_1 \omega_t^*({\mathbf{x}})\right\}\right| \\ \notag = & \max\bigg\{ \max _{{\mathbf{x}}\in D_t }\left\{\langle {\mathbf{x}}, {\mathbf{y}}\rangle -\mu_2 \omega_t^*({\mathbf{x}})\right\}-\max _{{\mathbf{x}}\in D_t }\left\{\langle {\mathbf{x}}, {\mathbf{y}}\rangle -\mu_1 \omega_t^*({\mathbf{x}})\right\}, \\ &\quad\quad\quad\quad\quad \max _{{\mathbf{x}}\in D_t }\left\{\langle {\mathbf{x}}, {\mathbf{y}}\rangle -\mu_1 \omega_t^*({\mathbf{x}})\right\} - \max _{{\mathbf{x}}\in D_t }\left\{\langle {\mathbf{x}}, {\mathbf{y}}\rangle -\mu_2 \omega_t^*({\mathbf{x}})\right\}\bigg\}\notag \\\notag \leq & \max\left\{ \max _{{\mathbf{x}}\in{D_t}} (\mu_1-\mu_2)\omega_t^*({\mathbf{x}}), \max _{{\mathbf{x}}\in{D_t}} (\mu_2-\mu_1)\omega_t^*({\mathbf{x}})\right\}\\ \notag = &\left(\mu_1-\mu_2\right) \max _{{\mathbf{x}}\in{D_t}}-\omega_t^*({\mathbf{x}})= \left(\mu_1-\mu_2\right) \left|\max _{{\mathbf{x}}\in{D_t}}[\langle\mathbf{0}, {\mathbf{x}}\rangle-\omega_t^*({\mathbf{x}})]\right| \\ = & \omega_t(\mathbf{0})\left(\mu_1-\mu_2\right) =\log(t)\left(\mu_1-\mu_2\right), \label{eq:boundbtilde} \end{align}\tag{49}\] where the fourth equality uses the fact that the conjugate function of \(\omega_t^*\) is \(\omega_t\) itself, and the inequality holds because for any continuous functions \(f_1, f_2: \mathbb{R}^t \rightarrow \mathbb{R}\), \[\max _{{\mathbf{u}}}\left\{f_1({\mathbf{u}})-f_2({\mathbf{u}})\right\}+\max _{{\mathbf{u}}}\left\{f_2({\mathbf{u}})\right\} \geq \max _{{\mathbf{u}}}\left\{f_1({\mathbf{u}})-f_2({\mathbf{u}})+f_2({\mathbf{u}})\right\}=\max _{{\mathbf{u}}}\left\{f_1({\mathbf{u}})\right\}.\] Since \(\mathcal{Z}\) is a finite discrete set, we let \(t=|\mathcal{Z}|\) and \(r_j =\Psi({\mathbf{y}},{\mathbf{z}}_j;{\mathbf{x}})\) for all \(j\in[t]\). We have \[\widetilde{\Phi}({\mathbf{y}},\mu;{\mathbf{x}})=\mu \log \mathbb{E}_{{\mathbf{z}}\sim \zeta} [e^{\Psi({\mathbf{y}},{\mathbf{z}};{\mathbf{x}})/\mu}] = \mu\log \frac{1}{t} + \mu \log \sum_{j=1}^t e^{\Psi({\mathbf{y}},{\mathbf{z}}_j;{\mathbf{x}})/\mu} =\widetilde{b}_t({\mathbf{r}},\mu)+\mu\log\frac{1}{t}.\] Thus it follows from 49 that \[\begin{align} |\widetilde{\Phi}({\mathbf{y}},\mu_1;{\mathbf{x}}) - \widetilde{\Phi}({\mathbf{y}},\mu_2;{\mathbf{x}})| \leq & |\widetilde{b}_t({\mathbf{r}},\mu_1) - \widetilde{b}_t({\mathbf{r}},\mu_2)| + \left|\log\frac{1}{t} \right| \left(\mu_1-\mu_2\right) \\ \leq &(\log(t)+\log({t}))(\mu_1-\mu_2)= 2\log(t)(\mu_1-\mu_2), \end{align}\] which indicates the desired result.

9.3 Proof of Lemma 4↩︎

Since \({\mathcal{Z}}\) is a connected compact set endowed with normalized Lebesgue measure, it admits diameter-bounded equal-measure partitions. Hence, there exists a constant \(C_{{\mathcal{Z}}}>0\), depending only on \({\mathcal{Z}}\) and \(m_2\), such that for every \(\rho\in(0,1]\), one can find a measurable partition \[{\mathcal{Z}}=\bigcup_{i=1}^{M_\rho}\mathcal{A}_i, \qquad \zeta(\mathcal{A}_i)=\frac{1}{M_\rho}, \qquad \operatorname{diam}(\mathcal{A}_i)\le \rho,\] with \[M_\rho\le \left(\frac{C_{{\mathcal{Z}}}}{\rho}\right)^{m_2}.\] Here \(C_{{\mathcal{Z}}}\) is a geometry-dependent constant; for instance, one may take \(C_{{\mathcal{Z}}}=1+2 D_{{\mathcal{Z}}}\) with \(D_{{\mathcal{Z}}}\) being the diameter of \({\mathcal{Z}}\). See, e.g., [83] for covering bounds of connected compact sets and [84] for diameter-bounded equal-measure partitions. In the subsequent lemma, we provide a more detailed and precise formulation of Lemma 4.

9.4 Proof of Lemma 5↩︎

9.5 Proof of Theorem 1↩︎

Before we give the proof of Theorem 1, we present a necessary lemma.

By Lemma 9, we know that any Clarke stationary point of problem 1 is also a directional stationary point.

(of Theorem 1) From the definition of \({\mathbf{y}}^{(k)}\), there exist \({\boldsymbol{\gamma}}_{{\mathbf{y}}}^{(k)} \in \partial \varphi ({\mathbf{y}}^{(k)}) \subset\mathbb{R}^{m_1}\) and \(0<\mu_k\leq \epsilon_k\) such that \[\label{eq:kkgamma} \mathbb{E}\left[\left(\mathrm{dist}\left(\mathbf{0}, \partial\widetilde{g}( {\mathbf{y}}^{(k)},\mu_k)\right)\right)^2\right]= \mathbb{E}\left[\left\| \mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}}[\nabla\widetilde{{\Phi}}( {\mathbf{y}}^{(k)},\mu_k;{\mathbf{x}})]+ {\boldsymbol{\gamma}}_{{\mathbf{y}}}^{{(k)}}\right\|^2\right] \leq \epsilon^2_k, \text{ for all }k\ge 0.\tag{52}\] Define \[\begin{align} &A_k(\delta): = \left\{\left\|\mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}}[\nabla\widetilde{{\Phi}}( {\mathbf{y}}^{(k)},\mu_k;{\mathbf{x}})]+ {\boldsymbol{\gamma}}_{{\mathbf{y}}}^{{(k)}}\right\|^2 \ge \delta \right\}, \text{ and }\\ &A^c_k(\delta): = \left\{\left\|\mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}}[\nabla\widetilde{{\Phi}}( {\mathbf{y}}^{(k)},\mu_k;{\mathbf{x}})]+ {\boldsymbol{\gamma}}_{{\mathbf{y}}}^{{(k)}}\right\|^2 < \delta \right\}. \end{align}\] By Markov’s inequality and 52 , we have for any \(\delta>0\), \[\begin{align} \sum_{k=0}^{\infty}\mathrm{Prob}\left( A_k(\delta) \right) \leq \sum_{k=0}^{\infty}\frac{1}{\delta} \mathbb{E}\left[\left\|\mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}}[\nabla\widetilde{{\Phi}}( {\mathbf{y}}^{(k)},\mu_k;{\mathbf{x}})]+ {\boldsymbol{\gamma}}_{{\mathbf{y}}}^{{(k)}}\right\|^2\right] \leq \frac{1}{\delta}\sum_{k=0}^{\infty}\epsilon^2_k<+\infty. \end{align}\] Setting \(\delta=\frac{1}{t}\) for \(t\in\mathbb{N}_{++}\), we apply the Borel-Cantelli Lemma [86] to the sequence of events \(\left(A_k(\delta): k\ge 0\right)\), and obtain \(\mathrm{Prob}(\limsup_{k\rightarrow\infty} A_k(\frac{1}{t}))=0.\) By the definition of the limit operator, \(\omega\) belongs to the event \[\Omega:=\left\{\omega: \lim_{k\rightarrow\infty}\left\|\mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}}[\nabla\widetilde{{\Phi}}( {\mathbf{y}}^{(k)}(\omega),\mu_k;{\mathbf{x}})]+ {\boldsymbol{\gamma}}_{{\mathbf{y}}}^{{(k)}}(\omega)\right\|=0\right\},\] if and only if \[\omega\in \bigcap_{t=1}^{\infty}\bigcup_{k=1}^{\infty}\bigcap_{i=k}^{\infty} A_i^c(\frac{1}{t})=\left(\bigcup_{t=1}^{\infty}\limsup_{k\rightarrow\infty}A_k(\frac{1}{t})\right)^c.\] It then follows \[\mathrm{Prob}\left( \Omega\right)=1-\mathrm{Prob}\left(\bigcup_{t=1}^{\infty}\limsup_{k\rightarrow\infty}A_k(\frac{1}{t})\right)\geq 1-\sum_{t=1}^{\infty}\mathrm{Prob}\left(\limsup_{k\rightarrow\infty}A_k(\frac{1}{t})\right)=1.\] Combining this with \(\mathrm{Prob}\left( \Omega\right)\leq 1\), we derive \(\mathrm{Prob}(\Omega)=1\). This implies ?? by the equation in 52 .

Define event \(\overline{\Omega}:= \{\omega: \lim_{k\rightarrow \infty}{\mathbf{y}}^{(j_k)}(\omega) = {\mathbf{y}}^{*}(\omega)\}\). According to our setting in this theorem, it holds \(\mathrm{Prob}(\overline{\Omega})=1\), and thus \(\mathrm{Prob}(\Omega \cap \overline{\Omega})=1\). Let \(\omega\) belong to the probability-one event on which \(\left\| \mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}} [\nabla_{\mathbf{y}}\widetilde{\Phi}({\mathbf{y}}^{(k)}(\omega),\mu_k;{\mathbf{x}})] +{\boldsymbol{\gamma}}_{\mathbf{y}}^{(k)}(\omega) \right\|\to0.\) Suppose \({\mathbf{y}}^{(j_k)}(\omega)\to {\mathbf{y}}^*(\omega)\). Define \({\mathbf{v}}_k(\omega):= \mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}} [\nabla_{\mathbf{y}}\widetilde{\Phi}({\mathbf{y}}^{(j_k)}(\omega),\mu_{j_k};{\mathbf{x}})].\) By Lemma 2, \(\|{\mathbf{v}}_k(\omega)\|\le l_\Psi\), hence \(\{{\mathbf{v}}_k(\omega)\}\) is bounded. Moreover, \({\boldsymbol{\gamma}}_{\mathbf{y}}^{(j_k)}(\omega)=-{\mathbf{v}}_k(\omega)+o(1)\) is also bounded. Passing to a further subsequence if necessary, we may assume \({\mathbf{v}}_k(\omega)\to {\mathbf{v}}^*(\omega)\) and \({\boldsymbol{\gamma}}_{\mathbf{y}}^{(j_k)}(\omega)\to{\boldsymbol{\gamma}}_{\mathbf{y}}^*(\omega)\). Next, the gradient-consistency result for the smoothing approximation, together with \(\mu_{j_k}\downarrow0\) and \({\mathbf{y}}^{(j_k)}\to{\mathbf{y}}^\star\), yields \[\operatorname{dist}\left( {\mathbf{v}}_k,\partial(g-\varphi)({\mathbf{y}}^\star) \right)\to0\] along the selected subsequence. Hence \({\mathbf{v}}^\star\in\partial(g-\varphi)({\mathbf{y}}^\star)\). Since \(g - \varphi\) is Clarke regular at \({\mathbf{y}}^*(\omega)\) by Lemma 9, and given that \(\lim_{k\to\infty} \mu_{j_k} = 0\), we invoke [55]. This allows us to conclude that \({\mathbf{v}}^*(\omega)\in \partial (g-\varphi)({\mathbf{y}}^*(\omega)).\) Since \({\boldsymbol{\gamma}}_{\mathbf{y}}^{(j_k)}(\omega)\in\partial\varphi({\mathbf{y}}^{(j_k)}(\omega))\) and \(\partial\varphi\) is outer semicontinuous [85], \({\boldsymbol{\gamma}}_{\mathbf{y}}^*(\omega)\in\partial\varphi({\mathbf{y}}^*(\omega)).\) Taking limits in \({\mathbf{v}}_k(\omega)+{\boldsymbol{\gamma}}_{\mathbf{y}}^{(j_k)}(\omega)\to0\) gives \(0={\mathbf{v}}^*(\omega)+{\boldsymbol{\gamma}}_{\mathbf{y}}^*(\omega)\in\partial g({\mathbf{y}}^*(\omega))\). Therefore \({\mathbf{y}}^*(\omega)\) is a Clarke stationary point of problem 1 . By Lemma 9, it is also a directional stationary point for all \(\omega\in \Omega\cap \overline{\Omega}\). Since \(\mathrm{Prob}(\Omega \cap \overline{\Omega}) = 1\), \({\mathbf{y}}^*\) is a directional stationary point of problem 1 almost surely.

9.6 Proof of Lemma 6↩︎

\((\)a\()\) From the updating rule [eq:updatex2] with \(\alpha_k=\frac{1}{L^{(k)}}\), it follows that \[\label{eq:KKT2} L^{(k)}\left({\mathbf{y}}^{(k+1)}-{\mathbf{y}}^{(k)}\right)+ \mathcal{G}({\mathbf{y}}^{(k)},\mu_{k}) + {\boldsymbol{\gamma}}_{{\mathbf{y}}}^{(k+1)} =\mathbf{0},\tag{53}\] for some \({\boldsymbol{\gamma}}_{{\mathbf{y}}}^{(k+1)}\in \partial \varphi({\mathbf{y}}^{(k+1)}).\) We derive \[\begin{align} \notag &\mathbb{E}\left[\widetilde{g}( {\mathbf{y}}^{(k+1)},\mu_k)- \widetilde{g}( {\mathbf{y}}^{(k)},\mu_k) \right]\\ \notag \leq& \mathbb{E}\left[\left\langle\mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}}[\nabla_{{\mathbf{y}}} \widetilde{{\Phi}}( {\mathbf{y}}^{(k)},\mu_k;{\mathbf{x}})], {\mathbf{y}}^{(k+1)}- {\mathbf{y}}^{(k)}\right\rangle\right]+\mathbb{E}\left[\frac{L^{(k)}}{2}\| {\mathbf{y}}^{(k+1)}- {\mathbf{y}}^{(k)}\|^2\right] \\ & +\mathbb{E}\left[ \varphi({\mathbf{y}}^{(k+1)}) - \varphi({\mathbf{y}}^{(k)})\right] \notag \stackrel{\eqref{eq:KKT2}}{=} \mathbb{E}\left[\left\langle -{\boldsymbol{\gamma}}_{{\mathbf{y}}}^{(k+1)}-L^{(k)}( {\mathbf{y}}^{(k+1)}- {\mathbf{y}}^{(k)}), {\mathbf{y}}^{(k+1)}- {\mathbf{y}}^{(k)}\right\rangle\right] \\ &+ \mathbb{E}\left[\frac{L^{(k)}}{2}\| {\mathbf{y}}^{(k+1)}- {\mathbf{y}}^{(k)}\|^2\right] \notag +\mathbb{E}\left[\left\langle\mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}}[\nabla_{{\mathbf{y}}} \widetilde{{\Phi}}( {\mathbf{y}}^{(k)},\mu_k;{\mathbf{x}})]-\mathcal{G}({\mathbf{y}}^{(k)},\mu_{k}), {\mathbf{y}}^{(k+1)}- {\mathbf{y}}^{(k)}\right\rangle\right] \\ \notag& + \mathbb{E}\left[\varphi({\mathbf{y}}^{(k+1)}) - \varphi({\mathbf{y}}^{(k)})\right]\\ \notag\leq & -\frac{L^{(k)}}{2}\mathbb{E} \left[\| {\mathbf{y}}^{(k+1)}- {\mathbf{y}}^{(k)}\|^2\right] + \mathbb{E}\left[-\left\langle {\boldsymbol{\gamma}}_{{\mathbf{y}}}^{{(k+1)}}, {\mathbf{y}}^{(k+1)}- {\mathbf{y}}^{(k)}\right\rangle + \varphi({\mathbf{y}}^{(k+1)}) - \varphi({\mathbf{y}}^{(k)})\right] \\\notag&\quad\quad\quad\quad\,\,+ \mathbb{E}\left[\frac{L^{(k)}}{4}\|{\mathbf{y}}^{(k+1)}- {\mathbf{y}}^{(k)}\|^2 + \frac{1}{L^{(k)}}\left\|\mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}} [\nabla_{{\mathbf{y}}}\widetilde{{\Phi}}_( {\mathbf{y}}^{(k)},\mu_k;{\mathbf{x}})]-\mathcal{G}({\mathbf{y}}^{(k)},\mu_{k})\right\|^2\right] \\ \leq &-\frac{L^{(k)}}{4}\mathbb{E}\left[\| {\mathbf{y}}^{(k+1)}- {\mathbf{y}}^{(k)}\|^2\right]+ \frac{1}{L^{(k)}}\mathbb{E}\left[\left\| \mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}}[\nabla_{{\mathbf{y}}}\widetilde{{\Phi}}( {\mathbf{y}}^{(k)},\mu_k;{\mathbf{x}})]-\mathcal{G}({\mathbf{y}}^{(k)},\mu_{k})\right\|^2\right], \label{eq:b2} \end{align}\tag{54}\] where the first inequality comes from the \(L^{(k)}\)-smoothness of \(\widetilde{\Phi}(\cdot,\mu;{\mathbf{x}})\) for each \({\mathbf{x}}\in{\mathcal{X}}\), the second inequality holds by Young inequality, and the last inequality results from the convexity of \(\varphi\) and \({\boldsymbol{\gamma}}_{{\mathbf{y}}}^{(k+1)}\in \partial \varphi({\mathbf{y}}^{k+1})\).

Next, we bound the second term in the right hand side of 54 , \[\begin{align} \notag & \mathbb{E}\left[\left\| \mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}}[\nabla_{{\mathbf{y}}}\widetilde{{\Phi}}( {\mathbf{y}}^{(k)},\mu_k;{\mathbf{x}})]-\mathcal{G}({\mathbf{y}}^{(k)},\mu_{k})\right\|^2\right] \\= \notag& \mathbb{E}\left[\left\| \mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}}[\nabla_{{\mathbf{y}}}\widetilde{{\Phi}}( {\mathbf{y}}^{(k)},\mu_k;{\mathbf{x}})]-\frac{1}{M_k}\sum^{M_k}_{j=1} \mathcal{G}_{k_j}({\mathbf{y}}^{(k)},\mu_{k})\right\|^2\right] \\ \notag\leq & 2\mathbb{E}\left[\left\| \frac{1}{M_k} \sum_{j=1}^{M_k}\nabla_{{\mathbf{y}}}\widetilde{{\Phi}}( {\mathbf{y}}^{(k)},\mu_k;{\mathbf{x}}_{k_j})-\frac{1}{M_k}\sum^{M_k}_{j=1} \mathcal{G}_{k_j}({\mathbf{y}}^{(k)},\mu_{k}) \right\|^2\right] \\ & \notag\quad\quad\quad\quad\quad + 2\mathbb{E}\left[\left\| \mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}}[\nabla_{{\mathbf{y}}}\widetilde{{\Phi}}( {\mathbf{y}}^{(k)},\mu_k;{\mathbf{x}})]-\frac{1}{M_k} \sum_{j=1}^{M_k}\nabla_{{\mathbf{y}}}\widetilde{{\Phi}} ( {\mathbf{y}}^{(k)},\mu_k;{\mathbf{x}}_{k_j}) \right\|^2\right] \\\notag \leq & 2\widehat{\epsilon}_k^2 + 2\mathbb{E}\left[\left\| \mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}}[\nabla_{{\mathbf{y}}}\widetilde{{\Phi}}( {\mathbf{y}}^{(k)},\mu_k;{\mathbf{x}})]-\frac{1}{M_k} \sum_{j=1}^{M_k}\nabla_{{\mathbf{y}}}\widetilde{{\Phi}} ( {\mathbf{y}}^{(k)},\mu_k;{\mathbf{x}}_{k_j}) \right\|^2\right] \\ = & {2\widehat{\epsilon}_k^2 + \frac{2}{M_k}\mathbb{E}\left[\left\| \mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}}[\nabla_{{\mathbf{y}}}\widetilde{{\Phi}}( {\mathbf{y}}^{(k)},\mu_k;{\mathbf{x}})]-\nabla_{{\mathbf{y}}}\widetilde{{\Phi}} ( {\mathbf{y}}^{(k)},\mu_k;{\mathbf{x}}_{k_1}) \right\|^2\right]} \leq 2\widehat{\epsilon}_k^2 + \frac{8}{M_k}l_{\Psi}^2\leq 4\widehat{\epsilon}_k^2, \label{eq:boundg2} \end{align}\tag{55}\] where the first inequality uses triangle inequality, the second one yields from 13 , the second equality holds because of the i.i.d. sampling, the third inequality uses \(\|\nabla_{{\mathbf{y}}} \widetilde{\Phi}({\mathbf{y}},\mu;{\mathbf{x}})\|\leq l_{\Psi}\) for all \({\mathbf{x}}\in{\mathcal{X}}\) from Lemma 2(b), and the last one follows by \(M_k\ge \lceil 4l_{\Psi}^2 \widehat{\epsilon}_k^{-2}\rceil\).

(b) We deduce that \[\begin{align} &\mathbb{E}\left[\left\|\mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}}[\nabla_{{\mathbf{y}}}\widetilde{{\Phi}}( {\mathbf{y}}^{(k+1)},\mu_k;{\mathbf{x}})]+ {\boldsymbol{\gamma}}_{{\mathbf{y}}}^{{(k+1)}}\right\|^2\right] \\ \leq & 2\mathbb{E}\left[\left\|\mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}}[\nabla_{{\mathbf{y}}}\widetilde{{\Phi}}( {\mathbf{y}}^{(k)},\mu_k;{\mathbf{x}})]+ {\boldsymbol{\gamma}}_{{\mathbf{y}}}^{{(k+1)}}\right\|^2\right] \\ & +2 \mathbb{E}\left[\left\|\mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}}[\nabla_{{\mathbf{y}}}\widetilde{{\Phi}}( {\mathbf{y}}^{(k)},\mu_k;{\mathbf{x}})]- \mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}}[\nabla_{{\mathbf{y}}}\widetilde{{\Phi}}( {\mathbf{y}}^{(k+1)},\mu_k;{\mathbf{x}})]\right\|^2\right] \\\leq & 2\mathbb{E}\left[\left\|\mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}}[\nabla_{{\mathbf{y}}}\widetilde{{\Phi}}( {\mathbf{y}}^{(k)},\mu_k;{\mathbf{x}})]+ {\boldsymbol{\gamma}}_{{\mathbf{y}}}^{{(k+1)}}\right\|^2\right] +2\mathbb{E}\left[(L^{(k)})^2\| {\mathbf{y}}^{(k+1)}- {\mathbf{y}}^{(k)}\|^2\right] \\ & \stackrel{\eqref{eq:KKT2}}{=} 2\mathbb{E} \left[\left\|-L^{(k)}( {\mathbf{y}}^{(k+1)}- {\mathbf{y}}^{(k)}) + \left(\mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}} [\nabla_{{\mathbf{y}}}\widetilde{{\Phi}}( {\mathbf{y}}^{(k)},\mu_k;{\mathbf{x}})]- \mathcal{G}({\mathbf{y}}^{(k)},\mu_{k})\right)\right\|^2\right]\\ &+2\mathbb{E}\left[(L^{(k)})^2\| {\mathbf{y}}^{(k+1)}- {\mathbf{y}}^{(k)}\|^2\right]\notag \\ \leq & \frac{5(L^{(k)})^2}{2} \mathbb{E}\left[ \| {\mathbf{y}}^{(k+1)}- {\mathbf{y}}^{(k)}\|^2\right] +10 \mathbb{E}\left[\left\|\mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}} [\nabla_{{\mathbf{y}}}\widetilde{{\Phi}}( {\mathbf{y}}^{(k)},\mu_k;{\mathbf{x}})]- \mathcal{G}({\mathbf{y}}^{(k)},\mu_{k})\right\|^2\right] \\ & \quad\quad\quad\quad\quad +2\mathbb{E}\left[(L^{(k)})^2\| {\mathbf{y}}^{(k+1)}- {\mathbf{y}}^{(k)}\|^2\right]\notag\\ =& \frac{9 (L^{(k)})^2}{2} \mathbb{E}\left[ \| {\mathbf{y}}^{(k+1)}- {\mathbf{y}}^{(k)}\|^2\right] +10 \mathbb{E}\left[\left\|\mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}}[\nabla_{{\mathbf{y}}}\widetilde{{\Phi}}( {\mathbf{y}}^{(k)},\mu_k;{\mathbf{x}})]- \mathcal{G}({\mathbf{y}}^{(k)},\mu_{k})\right\|^2\right] \\ \le& 18 L^{(k)} \left(\mathbb{E} \left[\widetilde{g}( {\mathbf{y}}^{(k)},\mu_k)- \widetilde{g}( {\mathbf{y}}^{(k+1)},\mu_k)\right] + \frac{4}{L^{(k)}}\widehat{\epsilon}_k^2 \right) +40 \widehat{\epsilon}_k^2 \\ =& 18 L^{(k)} \mathbb{E} \left[\widetilde{g}( {\mathbf{y}}^{(k)},\mu_k)- \widetilde{g}( {\mathbf{y}}^{(k+1)},\mu_k)\right] + {112}{}\widehat{\epsilon}_k^2, \end{align}\] where we use triangle inequality in the first inequality, Lemma 2(c) in the second one, and Young’s inequality in the third one. To obtain the last inequality, we have used \[\mathbb{E}[\|\mathbb{E}_{{\mathbf{x}}\sim\mathbb{P}}[\nabla_{{\mathbf{y}}}\widetilde{{\Phi}}( {\mathbf{y}}^{(k)},\mu_k;{\mathbf{x}})]-\mathcal{G}({\mathbf{y}}^{(k)},\mu_{k})\|^2] \le 4 \widehat{\epsilon}_k^2\] from 55 , and \(\frac{L^{(k)}}{4} \mathbb{E} \left[\| {\mathbf{y}}^{(k+1)}- {\mathbf{y}}^{(k)}\|^2\right]\leq -\mathbb{E}\left[\widetilde{g}( {\mathbf{y}}^{(k+1)},\mu_k)- \widetilde{g}( {\mathbf{y}}^{(k)},\mu_k)\right] + \frac{4}{ L^{(k)}}\widehat{\epsilon}_k^2\) from ?? . The proof is then completed.

9.7 Proof of Theorem 2↩︎

Summing up ?? for all \(k=k_1, k_1+1,\ldots,K-1\), and noticing \(L^{(k)}=C_2/\mu_k\), we have \[\begin{align} \notag &\sum_{k=k_1}^{K-1} \frac{\mu_k}{C_2}\mathbb{E}\left[\left(\mathrm{dist}\left(\mathbf{0}, \partial\widetilde{g}( {\mathbf{y}}^{(\tau+1)},\mu_\tau)\right)\right)^2\right] = \sum_{k=k_1}^{K-1} \frac{\mu_k}{C_2}\mathbb{E}\left[\left(\mathrm{dist}\left(\mathbf{0}, \partial\widetilde{g}( {\mathbf{y}}^{(k+1)},\mu_k)\right)\right)^2\right] \\ \leq & 18 \sum_{k=k_1}^{K-1} \mathbb{E} \left[\widetilde{g}( {\mathbf{y}}^{(k)},\mu_k)- \widetilde{g}( {\mathbf{y}}^{(k+1)},\mu_k)\right] + \sum_{k=k_1}^{K-1} \frac{112 \mu_k}{ C_2} \widehat{\epsilon}_k^2. \label{eq:sumgradient} \end{align}\tag{56}\] For the first term in the right hand side of 56 , we have \[\begin{align} \notag & \sum_{k=k_1}^{K-1} \mathbb{E} \left[\widetilde{g}( {\mathbf{y}}^{(k)},\mu_k)- \widetilde{g}( {\mathbf{y}}^{(k+1)},\mu_k)\right] \\ \notag= & \mathbb{E} \left[\widetilde{g}( {\mathbf{y}}^{(k_1)},\mu_{k_1})- \widetilde{g}( {\mathbf{y}}^{(K)},\mu_{K-1}) \right] +\sum_{k=k_1}^{K-2} \mathbb{E} \left[\widetilde{g}( {\mathbf{y}}^{(k+1)},\mu_{k+1})- \widetilde{g}( {\mathbf{y}}^{(k+1)},\mu_{k})\right] \\ \leq & \mathbb{E} \left[\widetilde{g}( {\mathbf{y}}^{(k_1)},\mu_{k_1})- \widetilde{g}( {\mathbf{y}}^{(K)},\mu_{K-1})\right] +\sum_{k=k_1}^{K-2} \omega\left(\mu_{k}- \mu_{k+1}\right), \label{eq:difffunc} \end{align}\tag{57}\] where the first inequality follows from Lemmas 3–4. Substituting 57 into 56 and multiplying each side by \(\frac{C_2}{\sum_{k=k_1}^{K-1}{\mu_k}}\) yields ?? .

9.8 Proof of Corollary 1↩︎

(a) We look at the three terms on the right-hand side of ?? . Since \(\mu_k=\epsilon\) and \(k_1=0\), we have \(\sum_{k=k_1}^{K-1}\mu_k=K\epsilon .\) The first term satisfies \[\frac{ 18C_2 \mathbb{E} \left[ \widetilde{g}({\mathbf{y}}^{(0)},\epsilon) - \widetilde{g}({\mathbf{y}}^{(K)},\epsilon) \right] }{ K\epsilon } \le \frac{ 18C_2 \left( \widetilde{g}({\mathbf{y}}^{(0)},\epsilon) - \min_{{\mathbf{y}}}\widetilde{g}({\mathbf{y}},\epsilon) \right) }{ K\epsilon } \le \frac{\epsilon^2}{2},\] where the last inequality follows from the definition of \(K\). The second term is zero because \(\mu_k=\epsilon\) for all \(k\in[K]\). The third term is \(\frac{ 112\sum_{k=0}^{K-1}\epsilon\widehat\epsilon_k^2 }{ K\epsilon } = 112\left(\frac{\epsilon}{16}\right)^2 \le \frac{\epsilon^2}{2}.\) Therefore, the right-hand side of ?? is at most \(\epsilon^2\). Hence \(\mathbb{E} \left[ \left( \operatorname{dist} \left( 0,\partial\widetilde{g}({\mathbf{y}}^{(\tau+1)},\mu_\tau) \right) \right)^2 \right] \le \epsilon^2 .\) Since \(\tau\) is sampled from \(\{0,1,\ldots,K-1\}\), there exists some \(0\le k<K\) satisfying the same bound.

(b) Since \(\mu_k=\mu_g^2\) and \(k_1=0\), we have \(\sum_{k=0}^{K-1}\mu_k=K\mu_g^2 .\) The first term on the right-hand side of ?? satisfies \[\frac{ 18C_2 \mathbb{E} \left[ \widetilde{g}({\mathbf{y}}^{(0)},\mu_g^2) - \widetilde{g}({\mathbf{y}}^{(K)},\mu_g^2) \right] }{ K\mu_g^2 } \le \frac{ 18C_2 \left( \widetilde{g}({\mathbf{y}}^{(0)},\mu_g^2) - \min_{{\mathbf{y}}}\widetilde{g}({\mathbf{y}},\mu_g^2) \right) }{ K\mu_g^2 } \le \frac{\epsilon_g^2}{2}.\] The second term is again zero. The third term is \(\frac{ 112\sum_{k=0}^{K-1}\mu_g^2\widehat\epsilon_k^2 }{ K\mu_g^2 } \le \frac{\epsilon_g^2}{2}.\) Therefore, Theorem 2 yields \(\mathbb{E} \left[ \left( \operatorname{dist} \left( 0,\partial\widetilde{g}({\mathbf{y}}^{(\tau+1)},\mu_g^2) \right) \right)^2 \right] \le\epsilon_g^2 .\)

By the proof of Lemma 5, and \(\Delta_{\mu_g^2}\le\omega(\mu_g^2)\), we obtain \[\begin{align} \left( \mathbb{E} \left[ \operatorname{dist} \left( 0,\partial^{\mu_g}g({\mathbf{y}}^{(\tau+1)}) \right)^2 \right] \right)^{1/2} & \le \left( 2\mathbb{E} \left[ \left( \operatorname{dist} \left( 0,\partial\widetilde{g}({\mathbf{y}}^{(\tau+1)},\mu_g^2) \right) \right)^2 \right] \right)^{1/2} + 2\sqrt{2L_\Psi\Delta_{\mu_g^2}} \\ & \le \sqrt{2}\epsilon_g + 2\sqrt{2L_\Psi\omega(\mu_g^2) }. \end{align}\] Thus the randomly sampled output \({\mathbf{y}}^{(\tau+1)}\) is a \(( \sqrt{\frac{2\omega(\mu_g^2)}{L_\Psi}},\sqrt{2}\epsilon_g + 2\sqrt{2L_\Psi\omega(\mu_g^2)})\)-Goldstein stationary point in expectation. Hence at least one iterate \({\mathbf{y}}^{(k+1)}\), \(0\le k<K\), satisfies the same bound. This proves part (b).

References↩︎