Accelerated Stochastic Algorithms for Convex-Concave Saddle-Point Problems

TL;DR

This paper introduces new stochastic primal-dual algorithms with restart schemes that improve convergence rates for convex-concave saddle-point problems, especially under strong convexity and sub-Gaussian noise conditions.

Contribution

The paper presents the first stochastic restart scheme for strongly convex saddle-point problems and a new primal-dual hybrid gradient-based algorithm with optimal oracle complexity.

Findings

Restart scheme outperforms existing methods in oracle complexity.

New algorithm achieves state-of-the-art oracle complexity for non-strongly convex cases.

Method is nearly optimal under known lower bounds.

Abstract

We develop stochastic first-order primal-dual algorithms to solve a class of convex-concave saddle-point problems. When the saddle function is strongly convex in the primal variable, we develop the first stochastic restart scheme for this problem. When the gradient noises obey sub-Gaussian distributions, the oracle complexity of our restart scheme is strictly better than any of the existing methods, even in the deterministic case. Furthermore, for each problem parameter of interest, whenever the lower bound exists, the oracle complexity of our restart scheme is either optimal or nearly optimal (up to a log factor). The subroutine used in this scheme is itself a new stochastic algorithm developed for the problem where the saddle function is non-strongly convex in the primal variable. This new algorithm, which is based on the primal-dual hybrid gradient framework, achieves the state-of-the-art oracle complexity and may be of independent interest.

Tables

Table 1. Table 1: Comparison of oracle complexities of algorithms to obtain an ϵ italic-ϵ \epsilon -expected duality gap when μ > 0 𝜇 0 \mu>0 .

Algorithm	Problem Class	Oracle Complexity
PDHG-type [19]	SPP($L$, $L_{xx}$, $L_{yx}$, 0, 0, $\mu$)	$O\left(\frac{L+L_{xx}+L_{yx}}{\sqrt{\mu ϵ}}\right)$
Mirror-Prox-B [23]	SPP($L$, $L_{xx}$, $L_{yx}$, 0, 0, $\mu$)	$O\left(\frac{L+L_{xx}}{\mu }\log \left(\frac{1}{ϵ}\right)+\frac{L_{yx}}{\sqrt{\mu ϵ}}\right)$
Algorithm 2S	SPP($L$, $L_{xx}$, $L_{yx}$, $L_{yy}$, $\sigma$, $\mu$)	(9)

Table 2. Table 2: Comparison of oracle complexities of algorithms to obtain an ϵ italic-ϵ \epsilon -expected duality gap when μ = 0 𝜇 0 \mu=0 .

Algorithm\tnotexftn:HPE	Problem Class	Oracle Complexity
PDHG-type [19]	SPP($L$, $L_{xx}$, $L_{yx}$, $L_{yy}$, 0, 0)	$O\left(\frac{L}{ϵ}+\frac{L_{xx}+L_{yx}+L_{yy}}{ϵ}\right)$
MP [32]	SPP($L$, $L_{xx}$, $L_{yx}$, $L_{yy}$, 0, 0)	$O\left(\frac{L}{ϵ}+\frac{L_{xx}+L_{yx}+L_{yy}}{ϵ}\right)$
SMP [24]	SPP($L$, $L_{xx}$, $L_{yx}$, $L_{yy}$, $\sigma$, 0)	$O\left(\frac{L}{ϵ}+\frac{L_{xx}+L_{yx}+L_{yy}}{ϵ}+\frac{{\left({\sigma }_{x,f}+{\sigma }_{x,\Phi }\right)}^2+{\sigma }_{y,\Phi }^2}{{ϵ}^2}\right)$
SAMP [12]	SPP($L$, $L_{xx}$, $L_{yx}$, $L_{yy}$, $\sigma$, 0)	$O\left(\sqrt{\frac{L}{ϵ}}+\frac{L_{xx}+L_{yx}+L_{yy}}{ϵ}+\frac{{\left({\sigma }_{x,f}+{\sigma }_{x,\Phi }\right)}^2+{\sigma }_{y,\Phi }^2}{{ϵ}^2}\right)$
Algorithm 1	SPP($L$, $L_{xx}$, $L_{yx}$, $L_{yy}$, $\sigma$, 0)	$O\left(\sqrt{\frac{L}{ϵ}}+\frac{L_{xx}+L_{yx}+L_{yy}}{ϵ}+\frac{{\left({\sigma }_{x,f}+{\sigma }_{x,\Phi }\right)}^2+{\sigma }_{y,\Phi }^2}{{ϵ}^2}\right)$

Table 3. Table 3 : Comparison of the running times in seconds (mead ± plus-or-minus \pm std) obtained by SRS and SMP to reach ε = 0.001 𝜀 0.001 \varepsilon=0.001 over five runs.

$n$	SRS	SMP
1000	$29.05\pm 1.32$	$33.45\pm 0.95$
2000	$66.19\pm 1.52$	$73.64\pm 1.98$
3000	$79.69\pm 2.33$	$88.75\pm 2.67$
4000	$100.63\pm 2.85$	$112.80\pm 2.69$

Equations

$$\min_{x\in\mathcal{X}}\max_{y\in\mathcal{Y}}\big{[}S(x,y)\triangleq f(x)+g(x)+\Phi(x,y)-J(y)\big{]},$$

$${\Vert \nabla f(x)-\nabla f(x^{\prime })\Vert }_{{\mathbb{X}}^{\ast }}\le L{\Vert x-x^{\prime }\Vert }_{\mathbb{X}},\forall x,x^{\prime }\in \mathcal{X},$$

$$f(x)\ge f(x^{\prime })+⟨\nabla f(x^{\prime }),x-x^{\prime }⟩+\frac{\mu }{2}\Vert x-x^{\prime }{\Vert }^2.$$

$${\Vert {\nabla }_x\Phi (x,y)-{\nabla }_x\Phi (x^{\prime },y)\Vert }_{{\mathbb{X}}^{\ast }}\le L_{xx}{\Vert x-x^{\prime }\Vert }_{\mathbb{X}},$$

$${\Vert {\nabla }_x\Phi (x,y)-{\nabla }_x\Phi (x,y^{\prime })\Vert }_{{\mathbb{X}}^{\ast }}\le L_{yx}{\Vert y-y^{\prime }\Vert }_{\mathbb{Y}},$$

$${\Vert {\nabla }_y\Phi (x,y)-{\nabla }_y\Phi (x^{\prime },y)\Vert }_{{\mathbb{Y}}^{\ast }}\le L_{yx}{\Vert x-x^{\prime }\Vert }_{\mathbb{X}},$$

$${\Vert {\nabla }_y\Phi (x,y)-{\nabla }_y\Phi (x,y^{\prime })\Vert }_{{\mathbb{Y}}^{\ast }}\le L_{yy}{\Vert y-y^{\prime }\Vert }_{\mathbb{Y}}.$$

$$S(x^{\ast },y)\le S(x^{\ast },y^{\ast })\le S(x,y^{\ast }),\forall (x,y)\in \mathcal{X}\times \mathcal{Y}.$$

$$f(x)\triangleq {\mathbb{E}}_{\xi \sim P}[\tilde{f}(x,\xi )]\text{mbox}and\Phi (x,y)\triangleq {\mathbb{E}}_{\zeta \sim Q}[\tilde{\Phi }(x,y,\zeta )],\forall x\in \mathbb{X},y\in \mathbb{Y},$$

$$f(x)\triangleq \frac{1}{n}\sum _{i=1}^n\tilde{f}(x,{\xi }_i)\text{mbox}and\Phi (x,y)\triangleq \frac{1}{m}\sum _{i=1}^m\tilde{\Phi }(x,y,{\zeta }_i),\forall x\in \mathbb{X},y\in \mathbb{Y}.$$

$$O\left(\bigg{(}\sqrt{\frac{L}{\mu}}+\frac{L_{xx}}{\mu}\bigg{)}\log\left(\frac{1}{\epsilon}\right)+\frac{L_{yx}}{\sqrt{\mu\epsilon}}+\frac{L_{yy}}{\epsilon}+\left(\frac{(\sigma_{x,f}+\sigma_{x,\Phi})^{2}}{\mu\epsilon}+\frac{\sigma_{y,\Phi}^{2}}{\epsilon^{2}}\right)\log\left(\frac{\log(1/\epsilon)}{\varsigma}\right)\right).$$

$$O\left(\bigg{(}\sqrt{\frac{L}{\mu}}+\frac{L_{xx}}{\mu}\bigg{)}\log\left(\frac{1}{\epsilon}\right)+\frac{L_{yx}}{\sqrt{\mu\epsilon}}+\frac{L_{yy}}{\epsilon}+\left(\frac{(\sigma_{x,f}+\sigma_{x,\Phi})^{2}}{\mu\epsilon}+\frac{\sigma_{y,\Phi}^{2}}{\epsilon^{2}}\right)\log\left(\frac{1}{\epsilon}\right)\right).$$

$$O\left(\sqrt{\frac{L}{ϵ}}+\frac{L_{xx}+L_{yx}+L_{yy}}{ϵ}+\frac{({\sigma }_{x,f}+{\sigma }_{x,\Phi }{)}^2+{\sigma }_{y,\Phi }^2}{{ϵ}^2}\right),$$

$$O\left(\sqrt{\frac{L}{\epsilon}}+\frac{L_{xx}+L_{yx}+L_{yy}}{\epsilon}+\frac{(\sigma_{x,f}+\sigma_{x,\Phi})^{2}+\sigma_{y,\Phi}^{2}}{\epsilon^{2}}\log\bigg{(}\frac{1}{\varsigma}\bigg{)}\right).$$

$$D_{h_{\mathcal{U}}}(u,u^{\prime })\triangleq h_{\mathcal{U}}(u)-h_{\mathcal{U}}(u^{\prime })-⟨\nabla h_{\mathcal{U}}(u^{\prime }),u-u^{\prime }⟩\ge (1/2){\Vert u-u^{\prime }\Vert }^2,\forall u\in \mathcal{U},\forall u^{\prime }\in {\mathcal{U}}^o,$$

$${\Omega }_{h_{\mathcal{U}}}\triangleq {\sup }_{u\in \mathcal{U},u^{\prime }\in {\mathcal{U}}^o}{D}_{h_{\mathcal{U}}}(u,u^{\prime }).$$

$$u^{\prime}\mapsto u^{+}\triangleq\argmin_{u\in\mathcal{U}}\big{[}\varphi_{\lambda}(u)\triangleq\varphi(u)+\left\langle{u^{*}},{u}\right\rangle+\lambda^{-1}D_{h_{\mathcal{U}}}(u,u^{\prime})\big{]}.$$

$$u_i^{+}=\frac{\exp ({\xi }_i)}{{\sum }_{j=1}^n\exp ({\xi }_j)},\forall i\in [n],$$

$$h_{\mathcal{U}}^{\ast }(U^{\ast })$$

$$={\sup }_{\sigma (U)\in {\Delta }_n}⟨\sigma (U^{\ast }),\sigma (U)⟩-{\sum }_{i=1}^n{\sigma }_i(U)\ln {\sigma }_i(U)$$

$$=\ln \left({\sum }_{i=1}^n\exp ({\sigma }_i(U^{\ast }))\right).$$

$${\sigma }_i(U^{+})=\frac{\exp ({\sigma }_i(\Lambda ))}{{\sum }_{j=1}^n\exp ({\sigma }_j(\Lambda ))},\forall i\in [n].$$

$$\displaystyle(\mathbb{P}):\;\,\min_{x\in\mathcal{X}}\Big{[}\bar{S}(x)\triangleq\sup_{y\in\mathcal{Y}}S(x,y)\Big{]}\quad\quad\mbox{and}\quad\quad(\mathbb{D}):\;\,\max_{y\in\mathcal{Y}}\left[\underline{S}(x)\triangleq\inf_{x\in\mathcal{X}}S(x,y)\right].$$

$$G(x,y)$$

$$y^{t+1}:=\underset{y\in \mathcal{Y}}{\mathrm{argmin}}J(y)-⟨s^t,y-y^t⟩+{\alpha }_t^{-1}{D}_{h_{\mathcal{Y}}}(y,y^t)$$

$${\tilde{x}}^{t+1}:=(1-{\beta }_t){\overline{x}}^t+{\beta }_t{x}^t$$

$$x^{t+1}:=\underset{x\in \mathcal{X}}{\mathrm{argmin}}g(x)+⟨{\hat{\nabla }}_x\Phi (x^t,y^{t+1},{\zeta }_x^t)+\hat{\nabla }f({\tilde{x}}^{t+1},{\xi }^t),x-x^t⟩+{\tau }_t^{-1}{D}_{h_{\mathcal{X}}}(x,x^t)$$

$$s^{t+1}:=(1+{\theta }_{t+1}){\hat{\nabla }}_y\Phi (x^{t+1},y^{t+1},{\zeta }_y^{t+1})-{\theta }_{t+1}{\hat{\nabla }}_y\Phi (x^t,y^t,{\zeta }_y^t)$$

$${\overline{x}}^{t+1}:=(1-{\beta }_t){\overline{x}}^t+{\beta }_t{x}^{t+1}$$

$${\overline{y}}^{t+1}:=(1-{\beta }_t){\overline{y}}^t+{\beta }_t{y}^{t+1}$$

$$t:=t+1$$

$$D_{h_{\mathcal{Y}}}(y,y^{\prime })\ge (1/2)\Vert y-y^{\prime }{\Vert }^2,$$

$$D_{h_{\mathcal{X}}}(x,x^{\prime })\ge (1/2)\Vert x-x^{\prime }{\Vert }^2,$$

$${\Omega }_{h_{\mathcal{X}}}\triangleq {\sup }_{x\in \mathcal{X},x^{\prime }\in {\mathcal{X}}^o}{D}_{h_{\mathcal{X}}}(x,x^{\prime }),{\Omega }_{h_{\mathcal{Y}}}\triangleq {\sup }_{y\in \mathcal{Y},y^{\prime }\in {\mathcal{Y}}^o}{D}_{h_{\mathcal{Y}}}(y,y^{\prime }).$$

$${\delta }_{y,\Phi }^t$$

$${\delta }_{x,\Phi }^t$$

$${\delta }_{x,f}^t$$

Full text

$-$

1 Introduction.

Let $\mathbb{X}$ and $\mathbb{Y}$ be two finite-dimensional real normed spaces with dual spaces $\mathbb{X}^{*}$ and $\mathbb{Y}^{*}$ respectively. Consider the following saddle-point problem (SPP)

[TABLE]

where $\mathcal{X}\subseteq\mathbb{X}$ and $\mathcal{Y}\subseteq\mathbb{Y}$ are nonempty closed and convex sets, and the functions $f:\mathbb{X}\to\overline{\mathbb{R}}\triangleq(-\infty,+\infty]$, $g:\mathbb{X}\to\overline{\mathbb{R}}$ and $J:\mathbb{Y}\to\overline{\mathbb{R}}$ are convex, closed and proper (CCP). In addition, the function $\Phi:\mathbb{X}\times\mathbb{Y}\to[-\infty,+\infty]$ is convex-concave, i.e., $\Phi(\cdot,y)$ is convex for any $y\in\mathbb{Y}$ and $\Phi(x,\cdot)$ is concave for any $x\in\mathbb{X}$. We assume that $f$, $g$, $J$ and $\Phi$ satisfy the following regularity conditions:

• •

$f$ is $L$-smooth on $\mathcal{X}$ (where $L\geq 0$). Specifically, $f$ is differentiable on $\mathcal{X}^{\prime}\supseteq\mathcal{X}$, where $\mathcal{X}^{\prime}$ is an open set in $\mathbb{X}$, and its gradient $\nabla f:\mathbb{X}\to\mathbb{X}^{*}$ is $L$-Lipschitz on $\mathcal{X}$, i.e.,

[TABLE]

where $\left\|\cdot\right\|_{\mathbb{X}^{*}}$ and $\left\|\cdot\right\|_{\mathbb{X}}$ denote the norms on $\mathbb{X}^{*}$ and $\mathbb{X}$ respectively.

• •

$f$ is $\mu$-strongly convex on $\mathcal{X}$ (where $\mu\geq 0$), i.e., for any $x,x^{\prime}\in\mathcal{X}$,

[TABLE]

In this work, we will consider both cases where $\mu=0$ and $\mu>0$.

• •

$g$ and $J$ admit tractable Bregman proximal projections on $\mathcal{X}$ and $\mathcal{Y}$, respectively (see Section 2 for details). Define $\mathsf{dom}\,g\triangleq\{x\in\mathbb{X}:g(x)<+\infty\}$ and $\mathsf{dom}\,J\triangleq\{y\in\mathbb{Y}:J(y)<+\infty\}$. To make (1) well-posed, we assume that $\mathsf{dom}\,g\cap\mathcal{X}\neq\emptyset$ and $\mathsf{dom}\,J\cap\mathcal{Y}\neq\emptyset$.

• •

$\Phi$ is differentiable on $\mathcal{X}^{\prime}\times\mathcal{Y}^{\prime}$, where $\mathcal{Y}^{\prime}\supseteq\mathcal{Y}$ is an open set in $\mathbb{Y}$. For any $(x,y)\in\mathcal{X}\times\mathcal{Y}$, denote the gradient of $\Phi(\cdot,y)$ and $\Phi(x,\cdot)$ by $x\mapsto\nabla_{x}\Phi(x,y)$ and $y\mapsto\nabla_{y}\Phi(x,y)$ respectively. For all $x,x^{\prime}\in\mathcal{X}$ and $y,y^{\prime}\in\mathcal{Y}$, we assume that there exist constants $L_{xx},L_{yx},L_{yy}\geq 0$ such that

[TABLE]

Note that the gradient operator $(x,y)\mapsto[\nabla_{x}\Phi(x,y),-\nabla_{y}\Phi(x,y)]$ is $M$-Lipschitz on $\mathcal{X}\times\mathcal{Y}$, where $M\triangleq L_{xx}+2L_{yx}+L_{yy}$.

Based on the assumptions above, we aim to design optimal (or nearly optimal) first-order algorithms that find a saddle-point $(x^{*},y^{*})\in\mathcal{X}\times\mathcal{Y}$ of Problem (1), i.e., $(x^{*},y^{*})$ that satisfies

[TABLE]

For well-posedness, we assume that such a saddle-point exists (see Assumption 3.2 for conditions that guarantee the existence).

1.1 Stochastic first-order oracles.

Since we aim to solve (1) via first-order information, we need to properly set up the oracle model. For generality, we do not assume that the exact gradients of $f$, $\Phi(\cdot,y)$ and $\Phi(x,\cdot)$ can be obtained. Rather, we only assume that we have access to the unbiased estimators of $\nabla f$, $\nabla\Phi(\cdot,y)$ and $\nabla\Phi(x,\cdot)$ (a.k.a., stochastic gradients), which we denote by $\hat{\nabla}f$, $\hat{\nabla}\Phi(\cdot,y)$ and $\hat{\nabla}\Phi(x,\cdot)$, respectively. In addition, we assume that the gradient noise on $\nabla f$, i.e., $\hat{\nabla}f-{\nabla}f$, has bounded second moment and we denote this bound as $\sigma_{x,f}^{2}$. Similarly, we also assume that gradient noises on $\nabla\Phi(\cdot,y)$ and $\nabla\Phi(x,\cdot)$ have bounded second-moments and denote the bounds as $\sigma^{2}_{x,\Phi}$ and $\sigma^{2}_{y,\Phi}$, respectively. In some situations, we will further assume that the gradient noises have sub-Gaussian distributions. For a formal description of the oracle model and a precise statement of the aforementioned assumptions, readers are referred to Section 3 and Assumption 3.2 respectively.

Indeed, the oracles described above are standard in the literature on stochastic approximation, which dates back to Robbins and Monro [45] and since then, has become a canonical approach to solve stochastic programming (SP) problems. In the standard SP formulation, the smooth functions $f$ and $\Phi$ are typically represented as expectations (see e.g., Nemirovski et al. [33]), i.e.,

[TABLE]

where $\xi$ and $\zeta$ denote the random variables with distributions $P$ (supported on $\mathcal{M}$) and $Q$ (supported on $\mathcal{Z}$) respectively, and the functions $\widetilde{f}:\mathbb{X}\times\mathcal{M}\to\mathbb{R}$ and $\widetilde{\Phi}:\mathbb{X}\times\mathbb{Y}\times\mathcal{Z}\to\mathbb{R}$ are such that $f$ and $\Phi$ satisfy the convexity and smoothness assumptions above. In particular, if we take $P=n^{-1}\sum_{i=1}^{n}\delta_{\xi_{i}}$ and $Q=m^{-1}\sum_{i=1}^{m}\delta_{\zeta_{i}}$, where $\{\xi_{i}\}_{i=1}^{n}$ and $\{\zeta_{i}\}_{i=1}^{m}$ are deterministic points in $\mathcal{M}$ and $\mathcal{Z}$, respectively, and $\delta_{\xi_{i}}$ denotes the delta measure at $\xi_{i}$ (and the same for $\delta_{\zeta_{i}}$), then $f$ and $\Phi$ in (6) assume finite-sum forms, i.e.,

[TABLE]

In this case, we can construct the stochastic (first-order) oracle for $f$ by first sampling an index set $\mathcal{B}$ from $[n]\triangleq\{1,\ldots,n\}$ uniformly randomly, and then output the gradient of $f_{\mathcal{B}}\triangleq\left\lvert\mathcal{B}\right\rvert^{-1}\sum_{i\in\mathcal{B}}\widetilde{f}(x,\xi_{i})$. The stochastic oracle for $\Phi$ can also be constructed in the same way.

From this point on, we will refer to the class of problems in (1) as SPP($L$, $L_{xx}$, $L_{yx}$, $L_{yy}$, $\sigma$, $\mu$), where $\sigma\triangleq\max\{\sigma_{x,f},\sigma_{x,\Phi},\sigma_{y,\Phi}\}$ represents the collective stochasticity in the (stochastic) gradients of $f$, $\Phi(\cdot,y)$ and $\Phi(x,\cdot)$. If $\sigma=0$, then (1) corresponds to a deterministic optimization problem.

Oracle complexity. To obtain an $\epsilon$-duality gap in expectation or with high probability (cf. Section 3.4), we will focus on analyzing the oracle complexity in terms of its dependence on $L$, $L_{xx}$, $L_{yx}$, $L_{yy}$, $\sigma_{x,f}$, $\sigma_{x,\Phi}$, $\sigma_{y,\Phi}$, $\mu$ and $\epsilon$. In addition, for the algorithms in this work and almost all the works in the literature, the number of calls to each of the oracle described above (which returns $\nabla f$, $\nabla\Phi(\cdot,y)$ or $\nabla\Phi(x,\cdot)$ or their stochastic versions) is the same. Therefore, in our complexity analysis and comparison of complexities with other algorithms, we do not distinguish among the these oracles. Instead, the word “oracle complexity” refers to the complexity of each of them.

1.2 Related work.

The SPP in (1) has a wide range of applications across many fields, including statistics, machine learning, operations research and game theory. When $\Phi$ is bilinear, i.e., there exists a (bounded) linear operator $\mathsf{A}:\mathbb{X}\to\mathbb{Y}^{*}$ such that $\Phi(x,y)=\langle{\mathsf{A}x},{y}\rangle$ (where $\left\langle{\cdot},{\cdot}\right\rangle:\mathbb{Y}^{*}\times\mathbb{Y}\to\mathbb{R}$ denotes the duality pairing between $\mathbb{Y}^{*}$ and $\mathbb{Y}$), the applications of (1) can be found in numerous previous works, e.g., Juditsky and Nemirovski [22, 23], Chambolle and Pock [10] and Zhao et al. [54]. Beyond bilinear $\Phi$, there are also rich applications, including two-player convex-concave zero-sum game (Chen et al. [12]), convex optimization with functional constraints (Boyd and Vandenberghe [8]), and kernel matrix learning (Lanckriet et al. [29]). Other applications can also be found in Balamurugan and Bach [1].

Previous works on solving convex-concave SPPs generally fall into two categories, depending on whether the primal-dual coupling term $\Phi$ is bilinear. Our focus will be on the non-bilinear SPPs (i.e., the problem where $\Phi$ is non-bilinear). Before doing that, we briefly review the works for bilinear SPPs first.

1.2.1 Bilinear SPPs.

This class of problems is indeed a special case of (1), i.e., when $L_{xx}=L_{yy}=0$. In recent years, both deterministic (i.e., $\sigma=0$) and stochastic (i.e., $\sigma>0$) versions of this problem have been thoroughly studied, for both $\mu=0$ and $\mu>0$. For the deterministic problems, some well-known algorithms include Nesterov smoothing (a.k.a., excessive gap technique, Nesterov [36, 35]), primal-dual hybrid gradient (PDHG, Chambolle and Pock [9, 10]), hybrid proximal extragradient-type algorithm (HPE-type, He and Monteiro [20]) and primal-dual operator splitting (e.g., Condat [13], Vũ [48] and Davis [15]). In addition, to tackle the stochastic problems, stochastic versions of these algorithms have also been developed, e.g., Chen et al. [11], Zhao and Cevher [53] and Zhao et al. [54].

1.2.2 Non-bilinear SPPs.

We first consider the case where $\mu=0$ and $\sigma=0$ (i.e., no primal strong convexity and all oracles are deterministic). When $\Phi$ is possibly nonsmooth, algorithms based on primal-dual subgradient have been developed in several works, including Nedić and Ozdaglar [31], Nesterov [38] and Juditsky and Nemirovski [22]. However, these methods typically incur high oracle complexity, i.e., $O(\epsilon^{-2})$ (where $\epsilon$ denotes the desired accuracy for the duality gap). As a result, they are not competitive when $\Phi$ is smooth (cf. (4a) to (4d)). The smoothness of $\Phi$ has been exploited in many algorithms to achieve better complexity results. These methods include Mirror-Prox (Nemirovski [32]), HPE-type algorithm (Kolossoski and Monteiro [27]) and PDHG-type algorithm (Hamedani and Aybat [19]). In particular, the last two algorithms are the extensions of their counterparts for solving bilinear SPPs. When $\sigma>0$, stochastic extensions of Mirror-Prox have been developed in the literature. Some representative works include the stochastic Mirror-Prox (SMP) method (Juditsky et al. [24]) and the stochastic accelerated Mirror-Prox (SAMP) method (Chen et al. [12]).

Unlike the case where $\mu=0$, there exist very few works that have considered the case where $\mu>0$. When $L_{yy}=0$ and $\sigma=0$ (i.e., the function $\Phi(x,\cdot)$ is linear and all oracles are deterministic), Juditsky and Nemirovski [23] and Hamedani and Aybat [19] have proposed algorithms, which are based on Mirror-Prox and PDHG respectively, that achieve better complexity than their counterparts that are designed for $\mu=0$. However, no algorithms exist when $L_{yy}=0$ or $\sigma=0$. Therefore, in this work, we seek to resolve two important questions:

1. (I)

Can we improve the oracle complexities of these two algorithms when $L_{yy}=0$ and $\sigma=0$?** 2. (II)

Can we develop an algorithm that works for all the cases where $\mu>0$, $L_{yy}>0$ and $\sigma>0$?

Indeed, we will provide affirmative answers to both questions above, by developing a stochastic restart scheme that is not only able to deal all with the cases listed in (II), but also significantly improves the oracle complexities of the algorithms in [23] and [19], even in the case where $L_{yy}=0$ and $\sigma=0$.

1.3 Main contributions.

Our main contributions are summarized below.

1. (I)

First, when $\mu>0$ (i.e., $f$ is strongly convex), we develop a (multi-stage) stochastic restart scheme (i.e., Algorithm 2S) for solving (1). Since SPPs have different structures from convex optimization problems (COPs), our stochastic restart scheme is different from those for COPs (e.g., Ghadimi and Lan [18]). Specifically, we use a distance-based quantity as the restart criterion, instead of the objective error. In addition, instead of focusing on expectation, we analyze the stochasticity via the error probability, which is obtained using techniques from finite-state Markov chains. (For detailed discussions, we refer readers to Section 4.4.1.) In principle, the subroutine used in our restart scheme can be any stochastic algorithm developed for the SPP in (1) when $\mu=0$ (i.e., $f$ is non-strongly convex). However, different subroutines may result in different oracle complexities. Therefore, to achieve the desired oracle complexity, we develop a new stochastic algorithm (i.e., Algorithm 1) for the case where $\mu=0$, based on the PDHG framework. (For details, see contribution (II) below.)

When the gradient noises obey sub-Gaussian distributions, to achieve an $\epsilon$-duality gap with probability (with probability) at least $1-\varsigma$, our scheme has oracle complexity

[TABLE]

Note that even in the deterministic case (i.e., $\sigma=0$), this complexity is strictly better than any in the previous works (cf. Table 1). Based on the complexity in (8), under very mild conditions on the nonsmooth functions $g$ and $J$, our scheme obtains an $\epsilon$-expected duality gap with oracle complexity

[TABLE]

Regarding the optimality of the complexity in (9), the complexities of $L_{xx}$ and $L_{yx}$ match the lower bounds derived in Nemirovskii and Yudin [34] and Ouyang and Xu [41], respectively. Additionally, the complexities of $\sigma_{x,f}+\sigma_{x,\Phi}$ and $\sigma_{y,\Phi}$ nearly match (up to $\log(1/\epsilon)$ factor) the lower bounds derived in Raginsky and Rakhlin [43]. The complexities of $L_{xx}$ and $L_{yy}$ are the best-known, although their lower bounds are not known, to the best of our knowledge.

2. (II)

Second, when $\mu=0$ (i.e., $f$ is non-strongly convex), we develop a stochastic algorithm (i.e., Algorithm 1) that is used as a subroutine in our restart scheme (see above), and also may be of independent interest by itself. We develop this algorithm by extending the PDHG framework, which was originally developed for bilinear SPPs (cf. Section 1.2.1), to handle the non-bilinear case. In addition, we incorporate the stochastic acceleration technique (see e.g., Lan [28]) into our algorithm, which indeed enables us to obtain the optimal oracle complexity for the smooth function $f$. By judicious choices of the algorithm parameters, we are able to obtain an $\epsilon$-expected duality gap with the state-of-the-art oracle complexity (cf. Table 2)

[TABLE]

when all the gradient noises have bounded second moments (but the distributions are not necessarily sub-Gaussian). Previously, the complexity in (10) has been achieved by the SAMP algorithm introduced in Chen et al. [12]. However, since the algorithm in [12] is based on Mirror-Prox, it significantly differs from our algorithm (which is based on PDHG). Consequently, our algorithm provides a valuable alternative approach to achieve the complexity in (10).

Regarding the optimality of the complexity in (10), we notice that the complexities for $L$ and $L_{yx}$ match the lower bounds derived in Ouyang and Xu [41], and the complexities for $\sigma_{x,f}+\sigma_{x,\Phi}$ and $\sigma_{y,\Phi}$ match the lower bounds derived in Nemirovskii and Yudin [34]. Therefore, all of these complexities are optimal. In addition, the complexities of $L_{xx}$ and $L_{yy}$ are also the best-known, although no lower bounds have been derived in the literature. (Note that in this case, by taking $\Phi(x,y)=\phi^{\mathrm{P}}(x)+\langle{\mathsf{A}x},{y}\rangle-\phi^{\mathrm{D}}(y)$, where the convex functions $\phi^{\mathrm{P}}$ and $\phi^{\mathrm{D}}$ are $L_{xx}$- and $L_{yy}$-smooth respectively, we can obtain lower complexity bounds $O(\sqrt{L_{xx}/\epsilon})$ and $O(\sqrt{L_{yy}/\epsilon})$. However, this essentially returns to the bilinear case. Therefore, these lower bounds may not be tight for the non-bilinear case.)

For the analysis of Algorithm 1, we also derive the following large-deviation-type convergence result, which complements the convergence result in expectation (cf. (10)). Specifically, if all the gradient noises obey sub-Gaussian distributions, we can obtain an $\epsilon$-duality gap with probability at least $1-\varsigma$ with oracle complexity

[TABLE]

1.4 Notation.

Denote the set of natural numbers by $\mathbb{N}\triangleq\{1,2,\ldots\}$ and define $\mathbb{Z}_{+}\triangleq\mathbb{N}\cup\{0\}$. For any finite-dimensional real normed space $\mathbb{U}$, we denote its dual space by $\mathbb{U}^{*}$. We denote the norms on $\mathbb{U}$ and $\mathbb{U}^{*}$ by $\left\|\cdot\right\|$ and $\left\|\cdot\right\|_{*}$ respectively. In addition, denote the duality pairing between $\mathbb{U}^{*}$ and $\mathbb{U}$ by $\left\langle{\cdot},{\cdot}\right\rangle:\mathbb{U}^{*}\times\mathbb{U}\to\mathbb{R}$. For any CCP function $h:\mathbb{U}\to\overline{\mathbb{R}}$, define its domain as $\mathsf{dom}\,h\triangleq\{u\in\mathbb{U}:h(u)<+\infty\}$. In addition, for any nonempty set $\mathcal{U}$ in $\mathbb{U}$, denote its interior by $\mathsf{int}\,\mathcal{U}$, its boundary by $\mathsf{bd}\,\mathcal{U}$ and its diameter by $D_{\mathcal{U}}\triangleq\sup_{u,u^{\prime}\in\mathcal{U}}\left\|u-u^{\prime}\right\|$. We also denote the indicator function of $\mathcal{U}$ as $\iota_{\mathcal{U}}:\mathbb{U}\to\overline{\mathbb{R}}$, i.e., $\iota_{\mathcal{U}}(u)\triangleq 0$ if $u\in\mathcal{U}$ and $\iota_{\mathcal{U}}(u)\triangleq+\infty$ otherwise. Note that $\iota_{\mathcal{U}}$ is closed and convex if and only if $\mathcal{U}$ is closed and convex.

1.5 Organization.

The paper is organized as follows. We first introduce some preliminary material in Section 2. Then, we present Algorithm 1 developed for the case where $f$ is non-strongly convex (i.e., $\mu=0$), together with convergence results. The detailed proofs of these results are deferred to Section 6. In Section 4, we develop our stochastic restart scheme for the case where $f$ is strongly convex (i.e., $\mu>0$), by using a modified version of Algorithm 1 as a subroutine. We analyze the oracle complexity for this scheme and discuss some technical issues.

2 Preliminaries.

We first introduce the distance generating function and Bregman proximal projection, followed by the primal and dual functions associated with $S(\cdot,\cdot)$ in (1).

2.1 Distance generating function and Bregman proximal projection.

Let $\mathbb{U}$ and $h$ be given as in Section 1.4. We say that $h$ is essentially smooth if $h$ is continuously differentiable on $\mathsf{int}\,\mathsf{dom}\,h\neq\emptyset$ and for any $u\in\mathsf{bd}\,\mathsf{dom}\,h$ and any sequence $\{u^{k}\}_{k\in\mathbb{Z}^{+}}\subseteq\mathsf{int}\,\mathsf{dom}\,h$ such that $u^{k}\to u$, $\|\nabla h(u^{k})\|_{*}\to+\infty$. Let $\mathcal{U}$ be any nonempty closed and convex set in $\mathbb{U}$. We call $h_{\mathcal{U}}$ a distance generating function (DGF) on $\mathcal{U}$ if it is continuous on $\mathcal{U}$, essentially smooth, and

[TABLE]

where $\mathcal{U}^{o}\triangleq\mathcal{U}\cap\mathsf{int}\,\mathsf{dom}\,h_{\mathcal{U}}$ and $D_{h_{\mathcal{U}}}:\mathcal{U}\times\mathcal{U}^{o}\to\mathbb{R}$ is called the Bregman distance associated with $h_{\mathcal{U}}$. Based on $D_{h_{\mathcal{U}}}(\cdot,\cdot)$, we define the Bregman diameter of $\mathcal{U}$ under $h_{\mathcal{U}}$ as

[TABLE]

In addition, for any $u^{\prime}\in\mathcal{U}^{o}$ and CCP function $\varphi:\mathbb{U}\to\overline{\mathbb{R}}$, define the Bregman proximal projection (BPP) of $u^{\prime}$ on $\mathcal{U}$ under $\varphi$ (associated with the DGF $h_{\mathcal{U}}$, $u^{*}\in\mathbb{U}^{*}$ and $\lambda>0$) as

[TABLE]

Note that if $\inf_{u\in\mathcal{U}}\varphi(u)>-\infty$ and $\mathcal{U}^{o}\cap\mathsf{dom}\,\varphi\neq\emptyset$, then the minimization problem in (14) always has a unique solution in $\mathcal{U}^{o}\cap\mathsf{dom}\,\varphi$ (see Lemma 7.1 in Appendix 7). We say that the function $\varphi$ has a tractable BPP on $\mathcal{U}$ if there exists a DGF $h_{\mathcal{U}}$ on $\mathcal{U}$ such that the minimization problem in (14) has a unique and easily computable solution in $\mathcal{U}^{o}\cap\mathsf{dom}\,\varphi$, for any $u^{*}\in\mathbb{U}^{*}$ and $\lambda>0$.

We now provide some examples of the quadruple $(\mathbb{U},\varphi,\mathcal{U},h_{\mathcal{U}})$ that satisfy the assumptions above. Clearly, if $\mathbb{U}$ is a Hilbert space with inner product $\langle{\cdot},{\cdot}\rangle$ and its induced norm $\left\|\cdot\right\|$, and $\mathcal{U}=\mathbb{U}$, we can take $h_{\mathcal{U}}=(1/2)\left\|\cdot\right\|^{2}$, so that the minimization problem in (14) becomes the usual proximal minimization problem associated with the function $\varphi$. Due to this, we will state three examples below where $\mathbb{U}$ has non-Hilbertian geometry and $\varphi\equiv 0$ (or equivalently, $\varphi=\iota_{\mathcal{U}}$ and the minimization problem in (14) becomes unconstrained). These examples can be found in several references, e.g., Nesterov [36], Juditsky and Nemirovski [22].

1. (E1)

We let $\mathbb{U}=(\mathbb{R}^{n},\left\|\cdot\right\|_{1})$, where $\left\|u\right\|_{1}\triangleq\sum_{i=1}^{n}\left\lvert u_{i}\right\rvert$, $\mathcal{U}=\Delta_{n}\triangleq\{u\in\mathbb{R}^{n}:u\geq 0,\sum_{i=1}^{n}u_{i}=1\}$ and $h_{\mathcal{U}}(u)=\sum_{i=1}^{n}u_{i}\ln u_{i}$. It is well-known that $h_{\mathcal{U}}$ is 1-strongly convex on $\Delta_{n}$ with respect to $\left\|\cdot\right\|_{1}$ (see e.g., Beck and Teboulle [3]), so the inequality in (12) is satisfied. In addition, the minimization problem in (14) has closed-form solution, namely,

[TABLE]

where $\xi_{i}\triangleq 1+\ln u^{\prime}_{i}-\lambda u^{*}_{i}$, for any $i\in[n]$. 2. (E2)

This example is the matrix analogy of (E1). Let $\mathbb{U}=(\mathbb{S}^{n},\left\|\cdot\right\|_{*})$, where $\mathbb{S}^{n}$ denotes the space of real symmetric $n\times n$ matrices, $\left\|U\right\|_{*}\triangleq\left\|\sigma(U)\right\|_{1}$ and $\sigma(U)\triangleq(\sigma_{1}(U),\ldots,\sigma_{n}(U))\in\mathbb{R}^{n}$ denotes the vector of eigenvalues of $U$. Note that $\mathbb{U}^{*}=(\mathbb{S}^{n},\left\|\cdot\right\|_{*})$, where $\|U^{*}\|_{*}\triangleq\max_{i\in[n]}\left\lvert\sigma_{i}(U^{*})\right\rvert$. In addition, let $\mathcal{U}=\{U\in\mathbb{S}^{n}:U\succeq 0,\operatorname{tr}(U)=1\}=\{U\in\mathbb{S}^{n}:\sigma(U)\geq 0,\textstyle{\sum}_{i=1}^{n}\sigma_{i}(U)=1\}$, i.e., the set of positive semi-definite matrices with unit trace and $h_{\mathcal{U}}(U)=\sum_{i=1}^{n}\sigma_{i}(U)\ln\sigma_{i}(U)$. Consider

[TABLE]

In Nesterov [37, Section 2.2], it is shown that $h_{\mathcal{U}}^{*}$ is differentiable and $\nabla h_{\mathcal{U}}^{*}$ is $1$-Lipschitz with respect to $\left\|\cdot\right\|_{*}$ on $\mathbb{S}^{n}$. Therefore, we conclude that $h_{\mathcal{U}}$ is 1-strongly convex with respect to $\left\|\cdot\right\|$ on $\mathcal{U}$. Additionally, the solution of the minimization problem in (14) admits closed form. Specifically, let $\Lambda\triangleq\nabla h_{\mathcal{U}}(U^{\prime})-\lambda U^{*}\in\mathbb{S}^{n}$ (see Ben-Tal and Nemirovski [5, Section 6.2] for the expression of $\nabla h_{\mathcal{U}}(U^{\prime})$), whose eigen-decomposition is denoted by $P\mathsf{diag}\,(\sigma(\Lambda))P^{T}$ (where $P$ is an orthogonal matrix). Then $U^{+}=P\mathsf{diag}\,(\sigma(U^{+}))P^{T}$, where

[TABLE] 3. (E3)

Finally, consider $\mathbb{U}=(\mathbb{R}^{n},\left\|\cdot\right\|_{p})$, where $p\in(1,2]$ and $\left\|u\right\|_{p}\triangleq\left(\sum_{i=1}^{n}\left\lvert u_{i}\right\rvert^{p}\right)^{1/p}$. Consequently, $\mathbb{U}^{*}=(\mathbb{R}^{n},\left\|\cdot\right\|_{q})$, where $q\triangleq 1/(1-p^{-1})\in[2,+\infty)$. Let $\mathcal{U}=\mathbb{R}_{+}^{n}$, where $\mathbb{R}_{+}\triangleq[0,+\infty)$ and $h_{\mathcal{U}}(u)=(1/2)\left\|u\right\|_{p}^{2}$. From Ben-Tal et al. [4, Section 8], we know that $h_{\mathcal{U}}$ is $(p-1)$-strongly convex with respect to $\left\|\cdot\right\|_{p}$ on $\mathbb{R}^{n}$. In addition, the minimization problem in (14) can be reduced to a one-dimensional optimization problem (which is easy to solve). To see this, from the KKT conditions, if $u^{*}_{i}\geq 0$, then $u^{+}_{i}=0$. Define \mathcal{I}^{\scalebox{0.8}{-}}\triangleq\{i\in[n]:u^{*}_{i}<0\}. If \begin{abstract}\mathcal{I}^{\scalebox{0.8}{-}}\end{abstract}\leq 1, the claim is trivially true. Otherwise, fix any i\in\mathcal{I}^{\scalebox{0.8}{-}} and for any j\in\mathcal{I}^{\scalebox{0.8}{-}}\setminus\{i\}, we have that $u_{j}=u_{i}/(u^{*}_{i}/u^{*}_{j})^{q-1}$. Thus the minimization problem in (14) still only involves one variable $u_{i}$.

For many constraint sets $\mathcal{U}$, compared with Hilbertian geometry, working under non-Hilbertian geometry offers better dependence on the dimension of $\mathbb{U}$ in the oracle complexities of several first-order algorithms (e.g., the subgradient and gradient methods). For detailed discussions, see e.g., Nemirovski et al. [33] and Juditsky and Nemirovski [22].

2.2 Primal function, dual function and duality gap.

For the SPP (1), we define the associated primal and dual problems as

[TABLE]

From the definition in (5), we can easily prove the following: Given that a saddle-point $(x^{*},y^{*})$ exists in (1), both $(\mathbb{P})$ and $(\mathbb{D})$ have nonempty solution sets $\mathcal{P}^{*}$ and $\mathcal{D}^{*}$, respectively. Furthermore, $x^{*}\in\mathcal{P}^{*}$ and $y^{*}\in\mathcal{D}^{*}$ and $\bar{S}(x^{*})=S(x^{*},y^{*})=\underline{S}(y^{*})$. Based on the functions $\bar{S}$ and $\underline{S}$, we define the duality gap

[TABLE]

3 Convex $f$: algorithm and convergence results.

We first consider the case where $\mu=0$. We begin with describing our algorithm, followed by the assumptions needed to analyze it, and finally its convergence results. The detailed convergence analysis is deferred to Section 6.

3.1 Algorithm description.

The pseudo-code of our algorithm is shown in Algorithm 1. For input, we require two CCP functions $h_{\mathcal{Y}}:\mathbb{Y}\to\overline{\mathbb{R}}$ and $h_{\mathcal{X}}:\mathbb{X}\to\overline{\mathbb{R}}$ which are DGFs on $\mathcal{Y}$ and $\mathcal{X}$ respectively, i.e., they are essentially smooth on their respective domains and

[TABLE]

where $\mathcal{Y}^{o}\triangleq\mathcal{Y}\cap\mathsf{int}\,\mathsf{dom}\,h_{\mathcal{Y}}$ and $\mathcal{X}^{o}\triangleq\mathcal{X}\cap\mathsf{int}\,\mathsf{dom}\,h_{\mathcal{X}}$. In addition, $h_{\mathcal{Y}}$ and $h_{\mathcal{X}}$ are chosen such that the minimization problems in (18) and (20) have (unique) and easily computable solutions in $\mathcal{Y}^{o}\cap\mathsf{dom}\,J$ and $\mathcal{X}^{o}\cap\mathsf{dom}\,g$ respectively (cf. Section 2). In addition, according to (13), we define the Bregman diameters of $\mathcal{X}$ and $\mathcal{Y}$ under $h_{\mathcal{X}}$ and $h_{\mathcal{Y}}$ as $\Omega_{h_{\mathcal{X}}}$ and $\Omega_{h_{\mathcal{Y}}}$ respectively, i.e.,

[TABLE]

We next formally describe the stochastic first-order oracles. For each $t\in\mathbb{N}$, the oracles return $\hat{\nabla}_{y}\Phi(x^{t},y^{t},\zeta^{t}_{y})$, $\hat{\nabla}_{x}\Phi(x^{t},y^{t+1},\zeta^{t}_{x})$ and $\hat{\nabla}f(\widetilde{x}^{t+1},\xi^{t})$, which are the unbiased estimators of the gradients $\nabla_{y}\Phi(x^{t},y^{t})$, $\nabla_{x}\Phi(x^{t},y^{t+1})$ and $\nabla f(\widetilde{x}^{t})$, respectively, conditioned on the past information. Here $\zeta^{t}_{y}$, $\zeta^{t}_{x}$ and $\xi^{t}$ are the underlying random variables that generate the stochasticity. For analysis purposes, let us define a filtration $\{\mathcal{F}_{t}\}_{t\in\mathbb{Z}_{+}}$ based on the stochastic process $\{(\zeta^{t}_{y},\zeta^{t}_{x},\xi^{t})\}_{t\in\mathbb{N}}$. Specifically, we first define the nested sequence of sets of random variables $\{\Xi_{t}\}_{t\in\mathbb{Z}_{+}}$ such that $\Xi_{0}\triangleq\{0\}$ and for any $t\in\mathbb{N}$, $\Xi_{t}\triangleq\{(\zeta_{x}^{i},\zeta_{y}^{i},\xi^{i})\}_{i=1}^{t}$. Then, for any $t\in\mathbb{Z}_{+}$, define $\mathcal{F}_{t}$ to be the $\sigma$-algebra generated by $\Xi_{t}$, i.e., the minimal $\sigma$-algebra with respect to which $\Xi_{t}$ is measurable. In addition, for any $t\in\mathbb{N}$, we define the stochastic gradient “noises”

[TABLE]

We briefly describe the structure of Algorithm 1, which falls under the PDHG framework (first introduced by Chambolle and Pock [9]). Specifically, we perform the dual ascent, primal descent and relaxation steps in (18), (20) and (21), respectively. These three steps constitute the main features of PDHG. Besides these steps, we also perform an interpolation step in (19) to accelerate the convergence of the smooth function $f$. Finally, we perform the averaging steps in (22) and (23) to obtain the ergodic iterate sequence $\{(\overline{x}^{t},\overline{y}^{t})\}_{t\in\mathbb{N}}$.

To measure the progress of Algorithm 1, we adopt the duality gap (defined in (17)) and analyze the convergence rate of the sequence $\{G(\overline{x}^{t},\overline{y}^{t})\}_{t\in\mathbb{N}}$ (in expectation or with high probability). Using the definition in (5), we easily see that if $\mathbb{E}[G(\overline{x}^{t},\overline{y}^{t})]\to 0$ and $\mathbb{E}[(\overline{x}^{t},\overline{y}^{t})]\to(x^{\dagger},y^{\dagger})$, then $(x^{\dagger},y^{\dagger})$ must be a saddle-point of Problem (1).

3.2 Assumptions.

Before presenting our convergence results, we first describe assumptions on the constraint sets $\mathcal{X}$ and $\mathcal{Y}$, as well as the stochastic gradient noises $\delta^{t}_{y,\Phi}$, $\delta_{x,\Phi}^{t}$ and $\delta^{t}_{x,f}$.

{assumption}

1. (A)

The Bregman diameters $\Omega_{h_{\mathcal{X}}}$ and $\Omega_{h_{\mathcal{Y}}}$ in (27) are bounded. 2. (B)

The set $\mathcal{X}$ is bounded and the Bregman diameter $\Omega_{h_{\mathcal{Y}}}$ is bounded.

{assumption}

Define the conditional expectation $\mathbb{E}_{t}[\cdot]\triangleq\mathbb{E}[\cdot\,|\,\mathcal{F}_{t}]$, for any $t\in\mathbb{Z}_{+}$. For any $x\in\mathcal{X}$ and $y\in\mathcal{Y}$ and any $t\in\mathbb{N}$, there exist positive constants $\sigma_{y,\Phi}$, $\sigma_{x,\Phi}$ and $\sigma_{x,f}$ such that

1. (A)

(Unbiasedness)   $\mathbb{E}_{t-1}[\delta_{y,\Phi}^{t}]=0$, $\mathbb{E}_{t-1}[\delta_{x,\Phi}^{t}]=0$, $\mathbb{E}_{t-1}[\delta^{t}_{x,f}]=0$ a.s., 2. (B)

(Bounded variance) $\mathbb{E}_{t-1}[\|\delta_{y,\Phi}^{t}\|_{*}^{2}]\leq\sigma_{y,\Phi}^{2}$, $\mathbb{E}_{t-1}[\|\delta_{x,\Phi}^{t}\|_{*}^{2}]\leq\sigma_{x,\Phi}^{2}$, $\mathbb{E}_{t-1}[\|\delta^{t}_{x,f}\|_{*}^{2}]\leq\sigma_{x,f}^{2}$ a.s., 3. (C)

(Sub-Gaussian distributions)

$\mathbb{E}_{t-1}\left[\exp\left(\|\delta_{y,\Phi}^{t}\|_{*}^{2}/\sigma_{y,\Phi}^{2}\right)\right]\leq\exp(1)$, $\mathbb{E}_{t-1}\left[\exp\left(\|\delta_{x,\Phi}^{t}\|_{*}^{2}/\sigma_{x,\Phi}^{2}\right)\right]\leq\exp(1)$,

$\mathbb{E}_{t-1}\left[\exp\left(\|\delta^{t}_{x,f}\|_{*}^{2}/\sigma_{x,f}^{2}\right)\right]\leq\exp(1)$ a.s..

3.3 Remarks on assumptions.

We provide remarks on the Assumptions 3.2 and 3.2.

3.3.1 Remarks on Assumption 3.2.

First, by (26), we note that Assumption 3.2(A) implies Assumption 3.2(B). These two assumptions will be used in proving different convergence results of the duality gap in the sequel, and appear to be common in the literature (see e.g., Juditsky and Nemirovski [22, 23]). In addition, note that in many scenarios, Assumptions 3.2(A) and (B) are both equivalent to the boundedness of $\mathcal{X}$ and $\mathcal{Y}$. For example, if both $\mathbb{X}$ and $\mathbb{Y}$ are finite dimensional real Hilbert spaces (with inner products $\langle{\cdot},{\cdot}\rangle_{\mathbb{X}}$ and $\langle{\cdot},{\cdot}\rangle_{\mathbb{Y}}$ and their induced norms $\left\|\cdot\right\|_{\mathbb{X}}$ and $\left\|\cdot\right\|_{\mathbb{Y}}$ respectively), and we take $h_{\mathcal{X}}(x)=(1/2)\left\|x\right\|_{\mathbb{X}}^{2}$ and $h_{\mathcal{Y}}(y)=(1/2)\left\|y\right\|_{\mathbb{Y}}^{2}$, then

[TABLE]

Next, we detail the purposes of Assumption 3.2. Since either part (A) or part (B) implies the boundedness of the constraint sets $\mathcal{X}$ and $\mathcal{Y}$, combined with other structural assumptions stated in Section 1, we see that Problem (1) has at least one saddle-point, hence being well-posed. In addition, the compactness of $\mathcal{X}$ and $\mathcal{Y}$ ensures that the sequence $\{\mathbb{E}[(\overline{x}^{t},\overline{y}^{t})]\}_{t\in\mathbb{N}}$ in Algorithm 1 has at least one limit point in $\mathcal{X}\times\mathcal{Y}$. Hence if $\mathbb{E}[G(\overline{x}^{t},\overline{y}^{t})]\to 0$ (which will be shown in Theorem 3.2), then any limit point of $\{\mathbb{E}[(\overline{x}^{t},\overline{y}^{t})]\}_{t\in\mathbb{N}}$ is a saddle-point of Problem (1). Moreover, in the convergence analysis of Algorithm 1 (cf. Section 6), Assumption 3.2(A) is also needed, since it provides a uniform upper bound on the sequence $\{D_{h_{\mathcal{X}}}(x,x^{t})\}_{t\in\mathbb{N}}$ and $\{D_{h_{\mathcal{Y}}}(y,y^{t})\}_{t\in\mathbb{N}}$, for any $(x,y)\in\mathcal{X}\times\mathcal{Y}$.

3.3.2 Remarks on Assumption 3.2.

First, we explain the implications of Assumption 3.2. In Assumption 3.2, part (A) states that the stochastic noise process $\{(\delta_{y,\Phi}^{t},\delta_{x,\Phi}^{t},\delta_{x,f}^{t})\}_{t\in\mathbb{Z}_{+}}$ forms a (vector-valued) martingale difference sequence (MDS) with respect to the filtration $\{\mathcal{F}_{t}\}_{t\in\mathbb{Z}_{+}}$. Part (B) states that the (conditional) second-moment of each of $\{\delta_{y,\Phi}^{t}\}_{t\in\mathbb{N}}$, $\{\delta_{x,\Phi}^{t}\}_{t\in\mathbb{N}}$ and $\{\delta_{x,f}^{t}\}_{t\in\mathbb{N}}$ is uniformly bounded. Under this assumption, we are able to prove certain finite-time convergence results of $\{G(\overline{x}^{t},\overline{y}^{t})\}_{t\in\mathbb{N}}$, both in expectation and with high probability (cf. Theorem 3.2(A)). In particular, for any $\varsigma\in(0,1]$, $\Pr\{G(\overline{x}^{T},\overline{y}^{T})=O(\varsigma^{-1})\}\geq 1-\varsigma$, where $T$ denotes the total number of iterations of Algorithm 1. In this work, however, we are also interested in obtaining large-deviation-type convergence results of the sequence $\{G(\overline{x}^{t},\overline{y}^{t})\}_{t\in\mathbb{N}}$, i.e., $\Pr\{G(\overline{x}^{T},\overline{y}^{T})=O(\log(1/\varsigma))\}\geq 1-\varsigma$. To achieve this, we need to assume that the (conditional) distributions of these stochastic noises are “light-tailed”. Specifically, in part (C), we assume that $\delta_{y,\Phi}^{t}$, $\delta_{x,\Phi}^{t}$ and $\delta_{x,f}^{t}$ are (conditional) sub-Gaussian random vectors with variance proxies $\sigma_{y,\Phi}^{2}$, $\sigma_{x,\Phi}^{2}$ and $\sigma_{x,f}^{2}$, respectively (see e.g., Vershynin [47]). As we will show in Theorem 3.2(B), this assumption indeed enables us to prove the aforementioned large-deviation-type results.

Next, we provide some justifications for Assumption 3.2(C). We will focus on the assumption $\mathbb{E}_{t-1}\left[\exp\left(\|\delta^{t}_{x,f}\|_{*}^{2}/\sigma_{x,f}^{2}\right)\right]\leq\exp(1)$ a.s., as the other two assumptions can be similarly justified. Note that under the compactness assumption of $\mathcal{X}$, for any $\xi\in\Xi$, if $\hat{\nabla}f(\cdot,\xi)$ is continuous on $\mathcal{X}$, then $\overline{M}_{f}(\xi)\triangleq\sup_{x\in\mathcal{X}}\|\hat{\nabla}f(x,\xi)\|_{*}<+\infty$. Similarly, by the $L$-smoothness of $f$ on $\mathcal{X}$ (cf. Section 1), we also have $M_{f}\triangleq\sup_{x\in\mathcal{X}}\|{\nabla}f(x)\|_{*}<+\infty$. Therefore, by the definition of $\delta^{t}_{x,f}$ , we have for any $t\in\mathbb{N}$,

[TABLE]

If ${\sup}_{\xi\in\Xi}\;\overline{M}^{2}_{f}(\xi)<+\infty$, which happens if $\overline{M}^{2}_{f}(\cdot)$ is continuous on $\Xi$ and $\Xi$ is compact, or $\Xi$ simply has finite size (which implies the finite-sum form of $f$ as in (7)), then we can take $\sigma_{x,f}^{2}=2(\sup_{\xi\in\Xi}\;\overline{M}^{2}_{f}(\xi)+M^{2}_{f})$.

Note that the justification above, although being general, leverages the compactness of $\mathcal{X}$ and $\Xi$. In fact, for some specific examples of $f$, the compactness assumptions are not needed. Let $\ell:\mathbb{R}\to\mathbb{R}$ be a loss function with bounded derivative on $\mathbb{R}$, i.e., $M_{\ell}\triangleq\sup_{s\in\mathbb{R}}\left\lvert\ell^{\prime}(s)\right\rvert<+\infty$. Examples of $\ell$ include the Huber loss function or the quadratically smoothed hinge loss function (see Zhang [50]). Let $f(x)=\mathbb{E}[\ell(\langle{\xi},{x}\rangle)]$, where $\xi\in\mathbb{X}^{*}$ is a sub-Gaussian random vector, i.e., $\mathbb{E}[\xi]=0$ and $\mathbb{E}[\exp(\left\|\xi\right\|^{2}_{*}/\sigma_{\xi}^{2})]\leq\exp(1)$, for some $\sigma_{\xi}>0$. (For convenience, we write $\xi\sim\mathsf{subG}(\sigma_{\xi}^{2})$.) Note that the sub-Gaussian assumption on $\xi$ is standard in statistical learning theory (see e.g., Vershynin [47]). By the Leibniz integral rule, we have $\nabla f(x)=\mathbb{E}[\ell^{\prime}(\langle{\xi},{x}\rangle)\xi]$. At each $t\in\mathbb{N}$, let $\hat{\nabla}f(x,\xi^{t})=\ell^{\prime}(\langle{\xi^{t}},{x}\rangle)\xi^{t}$, where $\xi^{t}\sim\mathsf{subG}(\sigma_{\xi}^{2})$. Then

[TABLE]

where in (a) we use $\left\lvert\ell^{\prime}(s)\right\rvert\leq M_{\ell}$, for any $s\in\mathbb{R}$, in (b) we use $\mathbb{E}[\|\xi\|^{2}_{*}]\leq\sigma_{\xi}^{2}$, in (c) we use Jensen’s inequality, and in (d) we use $\mathbb{E}[\exp(\left\|\xi^{t}\right\|^{2}_{*}/\sigma_{\xi}^{2})]\leq\exp(1)$. Therefore, we can take $\sigma_{x,f}^{2}=4M_{\ell}^{2}\sigma_{\xi}^{2}$.

3.4 Convergence results.

Before presenting our main results, we first introduce an important proposition that will be used in our stochastic restart scheme (cf. Section 4). Its proof is deferred to Section 6.

Proposition 3.1

Let Assumptions 3.2(A) and 3.2(A) hold. Define $\theta_{0}=0$, $\beta_{0}=2$, $\alpha_{0}=\tau_{0}=1$. (Note that these parameters are merely for analysis purpose, and do not appear in Algorithm 1, which starts with $t=1$.) If there exists a nonnegative sequence $\{\gamma_{t}\}_{t\in\mathbb{Z}_{+}}$ that satisfies $\gamma_{0}=0$ and

[TABLE]

for any $t\in\mathbb{N}$, then for any $T\geq 3$, we have the following:

1. (A)

If Assumption 3.2(B)* also holds, then*

[TABLE]

In addition, for any $\varsigma\in(0,1]$, with probability at least $1-\varsigma$, we have

[TABLE] 2. (B)

Define the stochastic sequence $\{\hat{x}^{t}\}_{t\in\mathbb{N}}$ such that $\hat{x}^{1}=x^{1}$ and

[TABLE]

If Assumption 3.2(C)* also holds, then for any $\varsigma\in(0,1/6)$, with probability at least $1-6\varsigma$, we have*

[TABLE]

Furthermore,

[TABLE]

Based on Proposition 3.1, we obtain the main convergence results for Algorithm 1.

Theorem 3.2

Let Assumptions 3.2(A) and 3.2(A) hold. In Algorithm 1, for any $t\in\mathbb{N}$, choose

[TABLE]

where $\rho,\rho^{\prime}>0$ are constants independent of $(L,L_{xx},L_{yx},L_{yy},\sigma_{x,f},\sigma_{x,\Phi},\sigma_{y,\Phi},t)$.

1. (A)

If Assumption 3.2(B)* also holds, then for any $T\geq 3$, we have*

[TABLE]

In addition, for any $\varsigma\in(0,1]$, with probability at least $1-\varsigma$, we have

[TABLE] 2. (B)

If Assumption 3.2(C)* also holds, then for any $T\geq 3$ and $\varsigma\in(0,1/6]$, with probability at least $1-6\varsigma$, we have*

[TABLE]

Proof 3.3

Proof. We first verify the choices of the input sequences $\{\beta_{t}\}_{t\in\mathbb{N}}$, $\{\alpha_{t}\}_{t\in\mathbb{N}}$, $\{\tau_{t}\}_{t\in\mathbb{N}}$ and $\{\theta_{t}\}_{t\in\mathbb{N}}$ in Theorem 3.2 indeed satisfy the conditions required in Proposition 3.1. Indeed, based on these choices, we can choose $\gamma_{t}=t$, for any $t\in\mathbb{N}$. We only show the steps to verify the conditions in (34). First, since $\tau_{t}^{-1}\geq(4L+2(L_{yx}+L_{xx})t)/t$ and $\alpha_{t}\leq 1/(16L_{yx})$, we have

[TABLE]

Also, since $\alpha_{t}^{-1}\geq 16L_{yy}$ and $\theta_{t}\leq 1$, we have

[TABLE]

Next, we bound the summation terms appearing in (35) and (42). Specifically, by noting that $\alpha_{t}\leq\left(16\rho\sigma_{y,\Phi}\sqrt{t}\right)^{-1}$ and $\tau_{t}\leq(2\rho^{\prime}\left(\sigma_{x,f}+\sigma_{x,\Phi})\sqrt{t}\right)^{-1}$, we have

[TABLE]

We then substitute (43), (44), (49), (50) and (51) into (35), (38) and (42) to obtain (45), (46) and (47). \Halmos

Remark 3.4

Note that the parameter choices in Theorem 3.2 do not involve the Bregman diameters $\Omega_{h_{\mathcal{X}}}$ and $\Omega_{h_{\mathcal{Y}}}$. However, if they are known (or can be estimated), we can choose $\rho={1}/(4\sqrt{\Omega_{h_{\mathcal{Y}}}})$ and $\rho^{\prime}=1/\sqrt{\Omega_{h_{\mathcal{X}}}}$ to “optimize” the bound in (45), (46) and (47). In Section 4.4.4, we will provide some methods to estimate $\Omega_{h_{\mathcal{X}}}$ and $\Omega_{h_{\mathcal{Y}}}$.

Remark 3.5

Theorem 3.2 indicates that to find a (stochastic) primal-dual pair $(x,y)\in\mathcal{X}\times\mathcal{Y}$ such that $\mathbb{E}[G(x,y)]\leq\epsilon$, the oracle complexity of Algorithm 1 is

[TABLE]

In addition, to find $(x,y)\in\mathcal{X}\times\mathcal{Y}$ such that $\Pr\{G(x,y)\leq\epsilon\}\geq 1-\varsigma$, the oracle complexity of Algorithm 1 is

[TABLE]

under Assumption 3.2(B) and

[TABLE]

under Assumption 3.2(C). Comparing (53) and (54), we observe that the sub-Gaussian assumption on the gradient noises (cf. Assumption 3.2(C)) enables us to obtain a significantly better complexity result, in terms of its dependence on $\varsigma$.

Remark 3.6

Note that as it becomes customary in the literature (see e.g., Lan [28], Ghadimi and Lan [18]), the optimality of the dependence of $B_{\mathrm{E}}(T)$ (i.e., the bound of $\mathbb{E}[G(\overline{x}^{T},\overline{y}^{T})]$ in (45)) on the diameters $\Omega_{h_{\mathcal{X}}}$, $\Omega_{h_{\mathcal{Y}}}$, $D_{\mathcal{X}}$ and $D_{\mathcal{Y}}$ is not discussed. As a result, we do not include these diameters in (52) and (54). With that said, the dependence of $B_{\mathrm{E}}(T)$ on these diameters is indeed the same as those in the existing methods (cf. Table 2). Similar comments also apply to all the oracle complexity results in Section 4.

4 Strongly Convex $f$: restart scheme and complexity analysis.

We consider the case where $\mu>0$. We aim to develop restart schemes based on Algorithm 1 that significantly improve the oracle complexities in (52) and (54).

4.1 Algorithm 1 with rescaled DGF.

We first introduce a variant of Algorithm 1 (developed for $\mu\geq 0$) that will be used as a subroutine in our restart scheme. Fix any $x_{\mathrm{c}}\in\mathcal{X}^{o}$ and $R>0$, and define a new DGF on $\bar{\mathcal{X}}(x_{\mathrm{c}},R)\triangleq R\mathcal{X}+x_{\mathrm{c}}$:

[TABLE]

*Based on the definition in (55), we easily see that the corresponding Bregman distance *

[TABLE]

for any $x\in\mathcal{X}$ and $x^{\prime}\in\mathcal{X}^{o}$. Since from (26), we know that $D_{h_{\mathcal{X}}}(x,x^{\prime})\geq(1/2)\|x-x^{\prime}\|^{2}$, by using the expression of $D_{\widetilde{h}_{\bar{\mathcal{X}}(x_{\mathrm{c}},R)}}(\cdot,\cdot)$ in (56), we also have $D_{\widetilde{h}_{\bar{\mathcal{X}}(x_{\mathrm{c}},R)}}(x,x^{\prime})\geq(1/2)\|x-x^{\prime}\|^{2}$.

Define $\mathcal{B}(x_{\mathrm{c}},R)\triangleq\{x\in\mathbb{X}:\left\|x-x_{\mathrm{c}}\right\|\leq R\}$. If $\mathcal{B}(0,1)\subseteq\mathsf{dom}\,h_{\mathcal{X}}$, we then have

[TABLE]

In words, the quantity $\Omega^{\prime}_{h_{\mathcal{X}}}$ is the Bregman diameter of the unit ball $\mathcal{B}(0,1)$ under $h_{\mathcal{X}}$, which can be interpreted as the normalized Bregman diameter of $\mathcal{B}(x_{\mathrm{c}},R)$ under $\widetilde{h}_{\bar{\mathcal{X}}(x_{\mathrm{c}},R)}$. Note that the condition $\mathcal{B}(0,1)\subseteq\mathsf{dom}\,h_{\mathcal{X}}$ is satisfied when $\mathbb{X}$ is a Hilbert space and $h_{\mathcal{X}}=(1/2)\left\|\cdot\right\|^{2}$, in which case $\mathsf{dom}\,h_{\mathcal{X}}=\mathbb{X}$. For some other examples, we refer readers to Nesterov [36, Section 4]* and Nemirovski [32, Section 5]**.*

Our modified algorithm is shown in Algorithm 1R. Compared to Algorithm 1, it has two differences. First, in Algorithm 1R, we fix the total number of steps $T-1$ before the algorithm starts. This is because Algorithm 1R will be used as the subroutine in our restart scheme, where at each stage, we terminate it once the total number of iterations is reached. Second, we replace step (20) in Algorithm 1 with step (58). Compared with step (20), in step (58), we change the constraint set from $\mathcal{X}$ to $\mathcal{X}^{\prime}$ and the Bregman distance from $D_{h_{\mathcal{X}}}(\cdot,\cdot)$ to $D_{\widetilde{h}_{\bar{\mathcal{X}}({x^{1},R})}}(\cdot,\cdot)$. The motivations for these changes will become apparent in the subsequent analysis.

In addition, we remark that if step (20) has an easily computable solution, then it is reasonable to assume that step (58) does as well. Three cases that guarantee this are:

• •

$\mathbb{X}$* is a normed space, $g\equiv 0$ and $\mathcal{X}^{\prime}=\mathcal{X}=\mathbb{X}$,*

• •

$\mathbb{X}$* is a Hilbert space, $\mathcal{X}^{\prime}=\mathcal{X}$ and $h_{\mathcal{X}}=(1/2)\left\|\cdot\right\|^{2}$,*

• •

$\mathbb{X}$* is a Hilbert space, $g\equiv 0$, $\mathcal{X}^{\prime}$ has an easily computable orthogonal projection and $h_{\mathcal{X}}=(1/2)\left\|\cdot\right\|^{2}$.*

4.2 Deterministic restart scheme for strongly convex $f$.

For ease of exposition, we first develop our restart scheme in the case where $\sigma=0$, i.e., we can obtain the deterministic gradients of $f$, $\Phi(\cdot,y)$ and $\Phi(x,\cdot)$. (The restart scheme for the stochastic case, where $\sigma_{x,f},\sigma_{x,\Phi},\sigma_{y,\Phi}>0$, will be developed in Section 4.3.) We start with analyzing the convergence properties of Algorithm 1R.

Proposition 4.1

Assume that $\sigma=0$, $\mathcal{B}(0,1)\subseteq\mathsf{dom}\,h_{\mathcal{X}}$ and $\Omega_{h_{\mathcal{Y}}}<+\infty$. In Algorithm 1R, fix

[TABLE]

If $x^{*}\in\mathcal{X}^{\prime}$, $R\geq 2\|x^{0}-x^{*}\|$, and we choose $\{\beta_{t}\}_{t=1}^{T}$ and $\{\theta_{t}\}_{t=1}^{T}$ as in (43), and $\alpha_{t}=\alpha$ and $\tau_{t}=t\tau$ for any $t\in[T]$, where

[TABLE]

then

[TABLE]

and \|\overline{x}^{T}-x^{*}\|\leq\sqrt{(2/\mu)B^{\rm det}_{R}(T)}\leq R/\big{(}2\sqrt{2}\big{)}.

Proof 4.2

*Proof. From the choices of $\{\beta_{t}\}_{t=1}^{T}$, $\{\alpha_{t}\}_{t=1}^{T}$, $\{\tau_{t}\}_{t=1}^{T}$ and $\{\theta_{t}\}_{t=1}^{T}$, we can easily verify that the conditions (32) to (34) in Proposition 3.1 continue to hold with $\gamma_{t}=t$, for any $t\in[T]$. In particular, $\gamma_{t}/\tau_{t}=\tau^{-1}$, for any $t\in[T]$ and $\gamma_{0}/\tau_{0}=0$. Therefore, by substituting the parameter choices in Proposition 4.1 and $x^{1}=\hat{x}^{1}=x^{0}$ into (40), we have *

[TABLE]

(Note that since in (58), the DGF $\widetilde{h}_{\bar{\mathcal{X}}(x^{1},R)}=\widetilde{h}_{\bar{\mathcal{X}}(x^{0},R)}$ is used, we need to replace the DGF $h_{\mathcal{X}}$ in (40) with $\widetilde{h}_{\bar{\mathcal{X}}(x^{0},R)}$ in order to obtain (62).) Next, recall from (16) that $\bar{S}(x)=f(x)+g(x)+\max_{y\in\mathcal{Y}}\Phi(x,y)-J(y)$ and define

[TABLE]

Note that since $f$ is $\mu$-strongly convex on $\mathcal{X}$, the same holds for $\bar{S}$. Based on $\bar{S}$ and $\underline{\widehat{S}}_{x^{0},R}$, we can define the $R$-restricted duality gap

[TABLE]

The second equality in (64) suggests us to take supremum over $x\in\mathcal{X}\cap\mathcal{B}(x^{0},R)$ and $y\in\mathcal{Y}$ on both sides of (62), and by using (57) and the value of $\eta$ in (60), we have

[TABLE]

On the other hand, by the $\mu$-strong convexity of $\bar{S}$, $x^{*}$ is the unique minimizer of $\bar{S}$ on $\mathcal{X}$, and

[TABLE]

where (a) follows from $\bar{S}(x^{*})=S(x^{*},y^{*})\geq S(x^{*},\overline{y}^{T})\geq{\min}_{x\in\mathcal{X}\bigcap\mathcal{B}(x^{0},R)}\;S(x,\overline{y}^{T})=\underline{\widehat{S}}_{x^{0},R}(\overline{y}^{T})$ as $x^{*}\in\mathcal{X}\cap\mathcal{B}(x^{0},R)$ (note that $\|x^{0}-x^{*}\|\leq R/2$ in Proposition 4.1).

Next, we aim to show that $\widehat{G}_{R}(\overline{x}^{T},\overline{y}^{T})=G(\overline{x}^{T},\overline{y}^{T})$. To start, suppose that

[TABLE]

By the input condition $\|x^{0}-x^{*}\|\leq R/2$, we have $\|x^{*}_{x^{0},R}(\overline{y}^{T})-x^{0}\|<R$. In other words, $x^{*}_{x^{0},R}(\overline{y}^{T})\in\mathcal{X}\cap\mathsf{int}\,\mathcal{B}(x^{0},R)$. Note that by its definition, there exists $d\in\partial_{x}S(x^{*}_{x^{0},R}(\overline{y}^{T}),\overline{y}^{T})$

[TABLE]

Since $x^{*}_{x^{0},R}(\overline{y}^{T})\in\mathsf{int}\,\mathcal{B}(x^{0},R)$, for any $x\in\mathcal{X}\setminus\mathcal{B}(x^{0},R)$, there exists $\lambda\in(0,1)$ such that $\overline{x}\triangleq\lambda x+(1-\lambda)x^{*}_{x^{0},R}(\overline{y}^{T})\in\mathcal{B}(x^{0},R)$. Moreover, $\overline{x}\in\mathcal{X}$ since $\mathcal{X}$ is convex. Thus $\overline{x}\in\mathcal{X}\cap\mathcal{B}(x^{0},R)$ and

[TABLE]

On the other hand, we have $\overline{x}-x^{*}_{x^{0},R}(\overline{y}^{T})=\lambda(x-x^{*}_{x^{0},R}(\overline{y}^{T}))$. Consequently,

[TABLE]

Combining (68) and (70), we have

[TABLE]

This indicates that $x^{*}_{x^{0},R}(\overline{y}^{T})=\argmin_{x\in\mathcal{X}}S(x,\overline{y}^{T})$ and hence $\widehat{G}_{R}(\overline{x}^{T},\overline{y}^{T})=G(\overline{x}^{T},\overline{y}^{T})$.

Now, it remains to show that (67) holds. First, since $S(\cdot,\overline{y}^{T})$ is $\mu$-strongly convex on $\mathcal{X}$, by the definition of $x^{*}_{x^{0},R}(\overline{y}^{T})$ and the fact that $x^{*}\in\mathcal{X}\cap\mathcal{B}(x^{0},R)$, we have

[TABLE]

On the other hand, $S(x^{*},\overline{y}^{T})\leq\max_{y\in\mathcal{Y}}S(x^{*},y)=\bar{S}(x^{*})\leq\bar{S}(\overline{x}^{T})$, since $x^{*}$ minimizes $\bar{S}$ on $\mathcal{X}$. Thus

[TABLE]

However, note that from (65) and the choice of $T$ in (59), we have $\widehat{G}_{R}(\overline{x}^{T},\overline{y}^{T})\leq\mu R^{2}/16$. We hence complete the proof. \Halmos

From Proposition 4.1, we observe that $\|x^{0}-x^{*}\|\leq R/2$ and \|\overline{x}^{T}-x^{*}\|\leq R/\big{(}2\sqrt{2}\big{)}. This suggests that if Algorithm 1R is used as the subroutine in a restart scheme, then at each stage, the radius $R$ can be reduced by a factor of $\sqrt{2}$, and accordingly, the bound of the duality gap (i.e., $\mu R^{2}/16$) can be halved. This observation naturally leads us to the restart scheme in Algorithm 2D, which comprises $K$ stages. At each stage $k$, using the output primal variable $\overline{x}_{k-1}^{T_{k-1}}$ from the previous stage as the input, we run Algorithm 1R for a sufficiently large number of iterations (i.e., $T_{k}$ iterations), so as to ensure the output primal variable $\overline{x}_{k}^{T_{k}}$ in the current stage satisfies that

[TABLE]

Theorem 4.3

Let Assumption 3.2(B) hold. In Algorithm 2D, for any desired accuracy $\epsilon\in(0,\mu U^{2}/4]$, we have $G(x_{K+1},y_{K+1})\leq\epsilon$, and the total number of oracle calls $C^{\rm det}_{\epsilon}$ satisfies that

[TABLE]

Before presenting the proof, from Theorem 4.3, we know that the oracle complexity for Algorithm 2D to obtain an $\epsilon$-duality gap is

[TABLE]

Proof 4.4

Proof. Note that $R_{k}=2^{(3-k)/2}U$, for any $k\in[K]$. Therefore,

[TABLE]

where (a) follows from (61) and (b) follows from K\geq\log_{2}\big{(}\mu U^{2}/(4\epsilon)\big{)}+1. In addition, we can also substitute the value of $R_{k}$ into (75) and obtain (76). \Halmos

Remark 4.5

By restricting $\epsilon\in(0,\mu U^{2}/4]$, we see that \max\big{\{}0,\log_{2}\big{(}\mu U^{2}/(4\epsilon)\big{)}\big{\}} simply becomes \log_{2}\big{(}\mu U^{2}/(4\epsilon)\big{)}. By doing so, we have indeed simplified the bound in (76) (as compared to the bound derived for $\epsilon>0$). On the other hand, note that if $\epsilon>\mu U^{2}/4$, then $K=1$, and Algorithm 2D effectively becomes Algorithm 1R. For the same reason, we will also focus on analyzing the regime $\epsilon\in(0,\mu U^{2}/4]$ in the stochastic restart scheme (see Theorem 4.9).

Remark 4.6

We observe that in Algorithm 2D, the choice of $T_{k}$ (cf. (75)) involves some large absolute constants, e.g., 4096 and 8192. These constants are mainly for theoretical purposes, and not meant to be tight. In practice, to avoid running Algorithm 1R for an excessively large number of iterations at each stage $k$ of Algorithm 2D, we can simply replace all the absolute constants in (75) (except 3) by a single constant $c>0$, and choose $c$ accordingly. The same remark also applies to the stochastic restart scheme (i.e., Algorithm 2S) in Section 4.3.2.

4.3 Stochastic restart scheme.

We now consider the general case where $\sigma_{x,f},\sigma_{x,\Phi},\sigma_{y,\Phi}>0$.

4.3.1 Intuition.

At the first attempt, one may try to combine the techniques used in the proofs of Theorem 3.2(A)* and Proposition 4.1, and then analyze the convergence of $\mathbb{E}[G(\overline{x}^{T},\overline{y}^{T})]$ in Algorithm 1R. However, a close inspection shows that such a combination does not work. Indeed, with this combination, in the proof of Proposition 4.1, although one can ensure that $\mathbb{E}[\widehat{G}_{R}(\overline{x}^{T},\overline{y}^{T})]\leq\mu R^{2}/16$ and hence \mathbb{E}\big{[}\|x^{*}_{x^{0},R}(\overline{y}^{T})-x^{*}\|\big{]}\leq R/\big{(}2\sqrt{2}\big{)} (cf. (73)) by choosing $T$ properly, these conditions cannot guarantee that $\mathbb{E}[\widehat{G}_{R}(\overline{x}^{T},\overline{y}^{T})]=\mathbb{E}[G(\overline{x}^{T},\overline{y}^{T})]$. This is because it is unclear how to ensure that $x^{*}_{x^{0},R}(\overline{y}^{T})\in\mathcal{X}\cap\mathsf{int}\,\mathcal{B}(x^{0},R)$ holds in the “expectation” sense, in contrast to the deterministic setting, where this condition holds a.s. With that said, when Assumption 3.2(C)** holds (i.e., the gradient noises are sub-Gaussian), it is possible to show that $\widehat{G}_{R}(\overline{x}^{T},\overline{y}^{T})\leq\mu R^{2}/16$ (and hence all the rest steps in the proof of Proposition 4.1, including $G(\overline{x}^{T},\overline{y}^{T})\leq\mu R^{2}/16$) holds with high probability (cf. Theorem 3.2(B)**). Moreover, if we “concatenate” this result in our restart scheme, we will still end up with a high-probability bound on the duality gap, as long as we keep the number of stages “reasonably” small.*

4.3.2 Algorithmic details.

Let us first analyze the convergence of Algorithm 1R in the stochastic setting.

Proposition 4.7

Assume that $\mathcal{B}(0,1)\subseteq\mathsf{dom}\,h_{\mathcal{X}}$, and let Assumptions 3.2(B), 3.2(A) and 3.2(C) hold. Fix any $\varsigma\in(0,1/6]$. In Algorithm 1R, choose $\mathcal{X}^{\prime}$ such that $x^{*}\in\mathcal{X}^{\prime}$ and $D_{\mathcal{X}^{\prime}}\leq R$, and choose

[TABLE]

If we choose $R\geq 2\|x^{0}-x^{*}\|$, $\{\beta_{t}\}_{t=1}^{T}$ and $\{\theta_{t}\}_{t=1}^{T}$ as in (43), and $\alpha_{t}=\alpha$ and $\tau_{t}=t\tau$ for any $t\in[T]$, where

[TABLE]

and $\eta=(4/R)\textstyle{\sqrt{\Omega_{h_{\mathcal{Y}}}/\Omega^{\prime}_{h_{\mathcal{X}}}}}$, then with probability at least $1-6\varsigma$,

[TABLE]

where $B^{\rm det}_{R}(T)$ is defined in (61) and

[TABLE]

Furthermore, \|\overline{x}^{T}-x^{*}\|\leq\sqrt{(2/\mu)(B^{\rm det}_{R}(T)+B^{\rm var}_{R}(T))}\leq R/\big{(}2\sqrt{2}\big{)} with probability at least $1-6\varsigma$.

Proof 4.8

Proof. The proof is adapted from that of Proposition 4.1, but with a careful stochastic argument. First, we substitute the values of $\{\beta_{t}\}_{t=1}^{T}$, $\{\alpha_{t}\}_{t=1}^{T}$, $\{\tau_{t}\}_{t=1}^{T}$ and $\{\theta_{t}\}_{t=1}^{T}$ into (40) and then take supremum over $x\in\mathcal{X}\cap\mathcal{B}(x^{0},R)$ and $y\in\mathcal{Y}$ on both sides of the resulting inequality. By noting that $D_{\mathcal{X}^{\prime}}\leq R$, we have that with probability at least $1-6\varsigma$,

[TABLE]

Denote the event in (83) as $\mathcal{A}_{T}$, on which we condition till the end of the proof. By the choice of $T$ in (79), we have with probability 1 that $\widehat{G}_{R}(\overline{x}^{T},\overline{y}^{T})\leq\mu R^{2}/16$. Consequently, from (66) and (73), we see that \|\overline{x}^{T}-x^{*}\|\leq R/\big{(}2\sqrt{2}\big{)} and \|x^{*}_{x^{0},R}(\overline{y}^{T})-x^{*}\|\leq R/\big{(}2\sqrt{2}\big{)} respectively with probability 1. As a result, we conclude that $x^{*}_{x^{0},R}(\overline{y}^{T})\in\mathcal{X}\cap\mathsf{int}\,\mathcal{B}(x^{0},R)$ and hence $\widehat{G}_{R}(\overline{x}^{T},\overline{y}^{T})=G(\overline{x}^{T},\overline{y}^{T})$. \Halmos

Based on Proposition 4.7, we present our stochastic restart scheme in Algorithm 2S. Compared to the deterministic restart scheme (i.e., Algorithm 2D), a notable difference is that at each stage $k$, the constraint set $\mathcal{X}^{\prime}$ becomes $\mathcal{X}\cap\mathcal{B}(x_{k},R_{k}/2)$ (as opposed to $\mathcal{X}$ as in Algorithm 2D). As will be shown in Theorem 4.9, this step enables us to obtain the nearly optimal oracle complexity on the primal noise term $\sigma_{x,f}+\sigma_{x,\Phi}$. Note that in general, the minimization problem in (58) with constraint set $\mathcal{X}^{\prime}=\mathcal{X}\cap\mathcal{B}(x_{k},R_{k}/2)$ is harder to solve, compared to the problem with constraint set $\mathcal{X}$. However, in many special cases, the problem with constraint set $\mathcal{X}^{\prime}$ can still be easily solved. For instance, let $\mathbb{X}$ be a Hilbert space and $h_{\mathcal{X}}=(1/2)\left\|\cdot\right\|^{2}$, so that $D_{\widetilde{h}_{\bar{\mathcal{X}}(x_{\mathrm{c}},R)}}(x,x^{\prime})=(1/2)\left\|x-x^{\prime}\right\|^{2}$. We can rewrite $\mathcal{B}(x_{k},R_{k}/2)$ as the functional constraint $\|x-x_{k}\|\leq R_{k}/2$, and then consider the Lagrangian form of this problem, i.e.,

[TABLE]

where $\pi^{t}\triangleq\hat{\nabla}_{x}\Phi(x^{t},y^{t+1},\zeta^{t}_{x})+\hat{\nabla}f(\widetilde{x}^{t+1},\xi^{t})$. We observe that for any fixed $\lambda\geq 0$, the minimization problem in (84) has the same form as that in (14), hence it has an easily computable solution. Furthermore, since $\lambda$ is a scalar, we can optimize it over $[0,+\infty)$ efficiently (e.g., via Newton’s method or bisection).

To analyze the convergence of Algorithm 2S, let us define the event

[TABLE]

where $\Omega$ denotes the sample space for the stochastic process $\{(x_{k},y_{k})\}_{k=1}^{K+1}$. For any $k\geq 2$, by conditioning on $\mathcal{E}_{k}$, from the proof of Proposition 4.7, we see that both $\|x_{k}-x^{*}\|\leq R_{k}/2$ (since $R_{k}=R_{k-1}/\sqrt{2}$) a.s. and

[TABLE]

Therefore, given $\mathcal{E}_{k}$, $\mathcal{X}_{k}$ satisfies all the requirements stated in Algorithm 1R and Proposition 4.7, i.e., $x^{*}\in\mathcal{X}_{k}\subseteq\mathcal{X}$ and $D_{\mathcal{X}_{k}}\leq R_{k}$. (Note that for $k=1$, we have $\mathcal{X}_{1}=\mathcal{X}$, thus these requirements are clearly satisfied.) As a result, Proposition 4.7 can be applied to the $k$-th stage of Algorithm 2S. Based on this observation, we derive the oracle complexity of Algorithm 2S below.

Theorem 4.9

Assume that $\mathcal{B}(0,1)\subseteq\mathsf{dom}\,h_{\mathcal{X}}$. Also, let Assumptions 3.2(B), 3.2(A) and (C) hold. In Algorithm 2S, for any $x_{0}\in\mathcal{X}^{o}$, desired accuracy $\epsilon\in(0,\mu U^{2}/4]$ and error probability $\nu\in(0,1]$, it holds that $G(x_{K+1},y_{K+1})\leq\epsilon$ with probability at least $1-\nu$. Furthermore, the total number of oracle calls $C^{\rm st}_{\epsilon}$ in Algorithm 2S satisfies that

[TABLE]

Before presenting the proof, from Theorem 4.9, we see that the oracle complexity of Algorithm 2S to obtain an $\epsilon$-duality gap with probability at least $1-\nu$ is

[TABLE]

Proof 4.10

Proof. For any $k=1,\ldots,K+1$, let $\mathbb{I}_{\mathcal{E}_{k}}$ denote the indicator function of the event $\mathcal{E}_{k}$ (cf. (85)). It is clear that $\{\mathbb{I}_{\mathcal{E}_{k}}\}_{k=1}^{K+1}$ forms a (finite-horizon) Markov chain, and therefore

[TABLE]

where (a) follows from Proposition 4.7, (b) follows from Bernoulli’s inequality and (c) follows from the choice of $\varsigma$ in Algorithm 2S. By the choice of $\{R_{k}\}_{k=1}^{K}$ in Algorithm 2S, we have

[TABLE]

Therefore, from (78), we know that $\mu R_{K}^{2}/16\leq\epsilon$. As a result,

[TABLE]

where (a) follows from (86). From the choice of $\varsigma$ and $K$, we have

[TABLE]

Now, we substitute both (91) and (93) into the choice of $T_{k}$ in (87), and then obtain the oracle complexity in (88). \Halmos

4.3.3 Complexity of convergence in expectation.

Based on the complexity results in Theorem 4.9, we aim to analyze the oracle complexity to obtain an $\varepsilon$-expected duality gap. To do so, we need one additional assumption on the nonsmooth functions $g$ and $J$.

{assumption}

*The nonsmooth functions $g$ and $J$ have closed domains. *

Remark 4.11

*Note that Assumption 4.3.3 is satisfied by many nonsmooth functions, e.g., an indicator function of a closed convex set or an absolutely homogeneous function (e.g., the norm function). *

To see the implication of Assumption 4.3.3, let us write down the explicit form of the primal and dual functions $\bar{S}$ and $\underline{S}$ in (16):

[TABLE]

In (94), by the closedness of $\mathsf{dom}\,J$ and the compactness of $\mathcal{Y}$ (cf. Assumption 3.2(B)), we see that $\mathsf{dom}\,J\cap\mathcal{Y}$ is compact. Furthermore, since $J$ is closed and convex, it is continuous on $\mathsf{dom}\,J\cap\mathcal{Y}$. Therefore, we can invoke Berge’s maximum theorem (see Theorem 8.1 in Appendix 8) to conclude that $x\mapsto{\max}_{y\in\mathcal{Y}\bigcap\mathsf{dom}\,J}\Phi(x,y)-J(y)$ is continuous on $\mathbb{X}$. Since $g$ is closed and convex, it is continuous on $\mathsf{dom}\,g\cap\mathcal{X}$, which is compact since $\mathsf{dom}\,g$ is closed and $\mathcal{X}$ is compact (cf. Assumption 3.2(B)). Therefore, we conclude that $\bar{S}$ is continuous on $\mathsf{dom}\,g\cap\mathcal{X}$. Similarly, we can show that $\underline{S}$ is continuous on $\mathsf{dom}\,J\cap\mathcal{Y}$. Hence, there exists a positive constant $\Gamma<+\infty$ such that

[TABLE]

Based on this observation, we can easily derive the following result.

Theorem 4.12

Assume that $\mathcal{B}(0,1)\subseteq\mathsf{dom}\,h_{\mathcal{X}}$. Also, let Assumptions 3.2(B), 3.2(A), 3.2(C), and 4.3.3 hold. In Algorithm 2S, for any $x_{0}\in\mathcal{X}^{o}$ and desired accuracy $\varepsilon\in(0,\mu U^{2}/2]$, if we choose K=\big{\lceil}\log_{2}\big{(}\mu U^{2}/(2\varepsilon)\big{)}\big{\rceil}+1 and $\nu=\min\{\varepsilon/(2\Gamma),1\}$, then $\mathbb{E}[G(x_{K+1},y_{K+1})]\leq\varepsilon$. Moreover, the oracle complexity of Algorithm 2S is

[TABLE]

Proof 4.13

Proof. Define the event $\mathcal{G}_{K,\varepsilon}\triangleq\{G(x_{K+1},y_{K+1})\leq\varepsilon/2\}$ and denote its complement as $\mathcal{G}_{K,\varepsilon}^{\mathrm{c}}$. From Theorem 4.9, we know that by the choice of $K$, $\Pr\{\mathcal{G}_{K,\varepsilon}\}\geq 1-\nu$, for any $\nu\in(0,1]$. Therefore,

[TABLE]

where (a) follows from (96) since $x_{K+1}\in\mathsf{dom}\,g\cap\mathcal{X}$ and $y_{K+1}\in\mathsf{dom}\,J\cap\mathcal{Y}$, and (b) follows from the choice of $\nu$ (which implies that $\nu\leq\varepsilon/(2\Gamma)$). To derive the oracle complexity in (97), we simply substitute $\epsilon=\varepsilon/2$ and $\nu=\varepsilon/(2\Gamma)$ into (89), and note that $O(\log(\log(1/\varepsilon)/\varepsilon))=O(\log(1/\varepsilon))$. \Halmos

4.4 Discussions.

We conclude this section by discussing some technical issues.

4.4.1 Restart schemes.

Note that our restart scheme in Algorithm 2S, which is developed for stochastic SPPs, greatly differs from the usual approaches for convex minimization problems (see e.g., Nesterov [39], Ghadimi and Lan [18] and in particular, Renegar and Grimmer [44]* for a comprehensive review). To be specific, consider $\min_{x\in\mathcal{X}}f(x)$, where $f$ is $L$-smooth and $\mu$-strongly convex on $\mathcal{X}$ (cf. Section 1) and denote the minimizer by $x^{*}$. Also, assume that we have access to stochastic gradients of $f$ that satisfy Assumption 3.2. Indeed, most of the existing restart schemes developed for this problem use the (primal) sub-optimality $f(x)-f(x^{*})$, either in expectation or with high probability, as the reduction criterion. To illustrate, consider the expected sub-optimality. These restart schemes require the subroutine to satisfy the following: for any starting point $\overline{x}^{1}\in\mathcal{X}$ and $\epsilon,\delta>0$, there exists $T\in\mathbb{N}$ such that*

[TABLE]

where $\overline{x}^{T}$ denotes $T$-th iterate produced by the subroutine. By the strong convexity of $f$, we have

[TABLE]

and therefore, we can establish an recursion between $\mathbb{E}[f(\overline{x}^{T})-f(x^{*})]$ and $\mathbb{E}[f(\overline{x}^{1})-f(x^{*})]$. (Indeed, by properly choosing $T$, this recursion becomes a contraction.) However, somewhat interestingly, this approach cannot be straightforwardly extended to solve stochastic SPPs, which would involve using the expected duality gap as the reduction criterion. There are at least two reasons behind this. First, by the definition of the duality gap in (17), we do not have the analog of (99), i.e., for any starting point $(\overline{x}^{1},\overline{y}^{1})\in\mathcal{X}^{o}\times\mathcal{Y}^{o}$ and $\epsilon,\delta>0$, there exists $T^{\prime}\in\mathbb{N}$ such that

[TABLE]

where recall that $(x^{*},y^{*})$ denotes a saddle-point of (1). Second, even if (101) holds, since we do not assume strong convexity on the dual function $\underline{S}$, there is no analog of (100) on the dual variable $\overline{y}^{1}$. Therefore, we still cannot establish a recursion between $\mathbb{E}[G(\overline{x}^{T^{\prime}},\overline{y}^{T^{\prime}})]$ and $\mathbb{E}[G(\overline{x}^{1},\overline{y}^{1})]$. As such, more sophisticated restart schemes need to be developed.

Inspired by Juditsky and Nemirovski [23]**, in the deterministic setting, we use the quantity $R$ as the reduction criterion (cf. Section 4.2). The quantity $R$ not only serves as an upper bound of the distance between the primal iterate $\overline{x}^{T}$ to $x^{*}$, but also of the duality gap $G(\overline{x}^{T},\overline{y}^{T})$ (cf. Proposition 4.1). Thus by sufficiently reducing $R$ (possibly over multiple stages), we can reduce the duality gap below any desired accuracy $\epsilon>0$. To achieve this, we modify Algorithm 1 into Algorithm 1R, whose convergence rate depends on $R$ and hence can be used as a subroutine in Algorithm 2D. Although our deterministic restart scheme is in the same spirit as the one in **[23]**, they are different in two important aspects. First, since our subroutine (i.e., Algorithm 1R) is based on PDHG, but the subroutine used in **[23]** is based on Mirror-Prox, our approach for developing the subroutine (cf. Section 4.1) is very different from that in **[23]**. As developing subroutines plays a crucial role in developing restart schemes, this difference clearly distinguishes our scheme from that in **[23]**. Second, in contrast to the approach in **[23]**, by using different DGFs (cf. (55)) and properly choosing stepsizes in our subroutine (cf. Propositions 4.1 and 4.7), our approach do not require defining new norms on $\mathbb{X}$ or $\mathbb{Y}$. Therefore, it is simpler and easier to understand.

In the stochastic setting, the problem becomes even more challenging. As discussed above (see also Section 4.3.1), the restart scheme based on the expected duality gap fails to work. Instead, we develop the restart scheme (i.e., Algorithm 2S) by leveraging the framework of finite-state Markov chain, which enables us to analyze the oracle complexity to obtain an $\epsilon$-duality gap with high probability. In addition, if the duality gap function $G(\cdot,\cdot)$ is bounded on $\mathcal{X}\times\mathcal{Y}$, we can also derive the oracle complexity to achieve an $\epsilon$-expected duality gap (cf. Section 4.3.3).

Shortly after the initial appearance of our paper on arXiv **[51], Yan et al. [49] also developed a stochastic restart scheme for a different class of SPP from the one studied in this work. Specifically, they assume $f\equiv 0$, $\mathbb{X}$ and $\mathbb{Y}$ are Hilbertian, and $\Phi(x,\cdot)$ is a linear function for any $x\in\mathcal{X}$. Also, instead of assuming $g$ and $J$ have easily computable BPPs on $\mathcal{X}$ and $\mathcal{Y}$, they assume a primal-dual stochastic sub-gradient oracle exists for $S(\cdot,\cdot)$. Moreover, such a stochastic sub-gradient is uniformly bounded on $\mathcal{X}\times\mathcal{Y}$ almost surely. As another important difference, they are interested in obtaining an $\epsilon$-primal sub-optimality gap, in expectation or with high probability, rather than an $\epsilon$-duality gap. Due to these differences, their restart scheme, as well as the subroutines therein, is different from ours. Interestingly, they consider the case where the primal function $\bar{S}$ (cf. (16)) satisfies the local error bound condition **[49, Definition 1]**, which might be worth exploring in our context. However, we defer this to future work.

4.4.2 Sub-Gaussian gradient noises.

In our stochastic restart scheme, we assume that the gradient noises $\{\delta_{y,\Phi}^{t}\}_{t\in\mathbb{N}}$, $\{\delta_{x,\Phi}^{t}\}_{t\in\mathbb{N}}$ and $\{\delta_{x,f}^{t}\}_{t\in\mathbb{N}}$ not only have uniformly bounded (conditional) second moments, but also follow sub-Gaussian distributions. The reason is that sub-Gaussianity leads to concentration, so that we have large-deviation-type results at each stage (cf. Proposition 4.7). However, these results fail to hold if the noises are “heavy-tailed”. Although we can still apply Markov’s inequality at each stage of Algorithm 2S (in lieu of concentration inequalities), the oracle complexity in (89) will have a significantly worse dependence on $\nu$, i.e., $O(1/\nu)$. Since in the proof of Theorem 4.12, we need to choose $\nu=\Theta(\epsilon)$, this worse dependence on $\nu$ leads to a worse dependence on $\epsilon$ in the complexity in (97). For example, to obtain an $\epsilon$-expected duality gap, the complexities corresponding to $(\sigma_{x,f}+\sigma_{x,\Phi})^{2}$ and $\sigma_{y,\Phi}^{2}$ degrade to $O((\sigma_{x,f}+\sigma_{x,\Phi})^{2}/(\mu\epsilon^{2}))$ and $O({\sigma_{y,\Phi}^{2}}/{\epsilon^{3}})$, respectively. Indeed, these complexities are worse than those in (52), which corresponds to the case where $\mu=0$. As future work, it would be interesting to develop new methodology that can effectively deal with the “heavy-tailed” noises.

4.4.3 Boundedness of $\mathcal{X}$.

We notice that the boundedness assumption on $\mathcal{X}$ is required by all the complexity results regrading both of the deterministic and stochastic restart schemes (i.e., Algorithms 2D and 2S). Since $f$ is $\mu$-strongly convex ($\mu>0$), in the following, we provide two ways to relax the assumption $D_{\mathcal{X}}<+\infty$ in Algorithm 2D, by assuming that $D_{\mathcal{Y}}<+\infty$. Note that this setting (i.e., $\mu>0$, $\mathcal{X}$ is unbounded, and $\mathcal{Y}$ is bounded) occurs in many applications. For example, consider the three-composite optimization problem (see e.g., Zhao et al. [54]**):

[TABLE]

where recall from Section 1.2 that $\mathsf{A}:\mathbb{X}\to\mathbb{Y}^{*}$ is a linear operator, and $h^{*}:\mathbb{Y}^{*}\to\mathbb{R}$ is $M_{h}$-Lipschitz on $\mathbb{Y}^{*}$ ($M_{h}<+\infty$). Instances of this problem include graph-guided fused Lasso (with full-column-rank data matrix) **[26]** and elastic-net regularization **[55]**. Note that the problem in (102) can be equivalently formulated into the SPP in (1), with $\Phi(x,y)=\langle{\mathsf{A}x},{y}\rangle$, and $\mathcal{Y}=\mathsf{dom}\,h$. In addition, since $M_{h}<+\infty$, we see that $\mathcal{Y}$ is indeed bounded.

To start with, let us observe that the assumption $D_{\mathcal{X}}<+\infty$ is only used in the input condition $U\geq D_{\mathcal{X}}$, as in general we can only use $D_{\mathcal{X}}$ to bound the initial distance $\|x^{1}-x^{*}\|$, i.e., $D_{\mathcal{X}}\geq\|x^{1}-x^{*}\|$. However, since $\bar{S}$ is $\mu$-strongly convex on $\mathcal{X}$, if we can find $\ell\in\mathbb{R}$ such that $\ell\leq\min_{x\in\mathcal{X}}\bar{S}(x)=\bar{S}(x^{*})$, then $\|x^{1}-x^{*}\|^{2}\leq({2}/{\mu})(\bar{S}(x^{1})-\bar{S}(x^{*}))\leq({2}/{\mu})(\bar{S}(x^{1})-\ell)$ and consequently,

[TABLE]

Thus, the assumption $D_{\mathcal{X}}<+\infty$ is no longer needed. We now provide a general approach to derive the lower bound $\ell$. Let us focus on the deterministic restart scheme (i.e., Algorithm 2D) first. Fix any $y\in\mathcal{Y}$. We then find an $\eta$-approximate optimal solution (where $\eta>0$) of the problem $\min_{x\in\mathcal{X}}S(x,y)$, denoted by $x^{*}_{\eta}(y)$, such that

[TABLE]

Since $f$ is $\mu$-strongly convex and $L$-smooth on $\mathcal{X}$ and $\Phi(\cdot,y)$ is $L_{xx}$-smooth on $\mathcal{X}$, from Lemma 9.1, to find $x^{*}_{\eta}(y)$, it suffices to find $\overline{x}\in\mathcal{X}$ such that its gradient mappings, i.e., $G_{\lambda}(\overline{x})\triangleq(\overline{x}-\overline{x}^{+})/\lambda$ and \overline{G}_{\lambda}(\overline{x})\triangleq\big{(}\nabla h_{\mathcal{X}}(\overline{x})-\nabla h_{\mathcal{X}}(\overline{x}^{+})\big{)}/\lambda, satisfy that

[TABLE]

where $0<\lambda\leq(L+L_{xx})^{-1}$ and

[TABLE]

and take $x^{*}_{\eta}(y)=\overline{x}^{+}$. Clearly, we can set the lower bound of $\bar{S}(x^{*})$, i.e., $\ell=S(x^{*}_{\eta}(y),y)-\eta$, since by the definition of $\bar{S}$ in (16),

[TABLE]

A potential issue of the approach above is that $\bar{S}(x^{1})$ in (103) may not be evaluated exactly. In this case, we need to find an upper bound $\upsilon\geq\bar{S}(x^{1})$. The way to find $\upsilon$ is similar to that for finding $\ell$. Specifically, we can find an $\eta$-approximate optimal solution of the problem $\max_{y\in\mathcal{Y}}\,\Phi(x^{1},y)-J(y)-(\eta/\Upsilon_{h_{\mathcal{Y}}})h_{\mathcal{Y}}(y)$, denoted by $y^{*}_{\eta}(x^{1})$, where \Upsilon_{h_{\mathcal{Y}}}=\max_{y\in\mathcal{Y}}\begin{abstract}h_\mathcal{Y}(y)\end{abstract}<+\infty. As a result,

[TABLE]

To find $y^{*}_{\eta}(x^{1})$, following Lemma 9.1, it suffices to find $\overline{y}\in\mathcal{Y}$ such that

[TABLE]

where $0<\lambda\leq L_{yy}^{-1}$ (recall that $\Phi(x,\cdot)$ is $L_{yy}$-smooth on $\mathcal{Y}$) and

[TABLE]

and set $y^{*}_{\eta}(x^{1})=\overline{y}^{+}$.

The second approach aims to write the SPP in (1) into an equivalent SPP that has a bounded primal constraint set, and so it is independent of Algorithm 2D. For any $y\in\mathcal{Y}$, let us define

[TABLE]

From Lemma 10.1 (see Appendix 10), we know that $x^{*}(\cdot)$ is $(L_{yx}/\mu)$-Lipschitz on $\mathcal{Y}$. Thus, we have

[TABLE]

To bound $\|x^{*}(y)\|$, let us find $x^{*}_{\eta}(y)$ in the same way as in the first approach, and by definition, $x^{*}_{\eta}(y)$ satisfies (104). As a result,

[TABLE]

If we combine (112) and (113), then we have

[TABLE]

If we define $\widetilde{\mathcal{X}}\triangleq\{x^{*}(y):y\in\mathcal{Y}\}\neq\emptyset$, then from (114), we see that $\widetilde{\mathcal{X}}\subseteq\mathcal{X}\cap\mathcal{B}(0,\varkappa)$. As a result, the SPP in (1) can be equivalently written as

[TABLE]

Note that in this new SPP, the primal constraint set $\mathcal{X}\cap\mathcal{B}(0,\varkappa)$ is nonempty, convex and compact.

Although the two approaches above work well in Algorithm 2D, there is one major difficulty that prevents us from extending either of them to the stochastic case. Recall from Section 1.1 that we only have access to the stochastic gradients of $f$ and $\Phi(\cdot,\cdot)$. Therefore, in the first approach, the proximal point $\overline{x}^{+}$ defined in (106) becomes stochastic, hence the gradient mappings $G_{\lambda}(\overline{x})$ and $\overline{G}_{\lambda}(\overline{x})$, as well as $x^{*}_{\eta}(x)$ (since $x^{*}_{\eta}(x)=\overline{x}^{+}$). One of the most natural ways to handle such stochasticity is to let the event in (105) hold with high probability (say, no less than $1-\delta/4$), and perform the rest of the analysis by conditioning on this event. However, note that we only have information about the distributions of $\hat{\nabla}f(x)$ and $\hat{\nabla}\Phi(x,y)$. In addition, the mappings from $\hat{\nabla}f(x)+\hat{\nabla}\Phi(x,y)$ to $\overline{x}^{+}$, as well as from $\overline{x}^{+}$ to $\|G_{\lambda}(\overline{x})\|^{2}+\|\overline{G}_{\lambda}(\overline{x})\|_{*}^{2}$, are rather complicated. Therefore, it is unclear how to lower bound the probability of the event in (105). Certainly, there may exist more sophisticated approaches to relax the assumption $D_{\mathcal{X}}<+\infty$ in the stochastic case, and we leave the investigation of these approaches to future work. As a final remark, if both of the sample spaces ${\Xi}$ and $\mathcal{Z}$ have finite sizes (cf. (7)), we will still have access to the exact gradients $f$ and $\Phi(\cdot,\cdot)$. As a result, the aforementioned two approaches continue to work in this case.

4.4.4 Unknown problem parameters.

Note that our restart schemes (i.e., Algorithms 2D and 2S) require the knowledge of several problem parameters. These include the strong convexity parameter $\mu$, the smoothness parameters $L$, $L_{xx}$, $L_{xy}$ and $L_{yy}$, the diameters $\mathcal{D}_{\mathcal{X}}$ and $\Omega_{h_{\mathcal{Y}}}$, and the variance parameters $\sigma^{2}_{x,f}$, $\sigma^{2}_{x,\Phi}$ and $\sigma^{2}_{y,\Phi}$. In the following, we will discuss possible ways to estimate some of these parameters, by assuming the knowledge of the rest.

To begin with, let us focus on unknown $\mathcal{D}_{\mathcal{X}}$, or more precisely, unknown $\|x^{1}-x^{*}\|$. Note that in Section 4.4.3, we already provide two ways to estimate $\|x^{1}-x^{*}\|$, by assuming that $D_{\mathcal{Y}}<+\infty$. These methods apply to both the deterministic case and the stochastic case with finite sample spaces $\Xi$ and $\mathcal{Z}$ (i.e., the finite-sum case). In the future, it would be interesting to develop a procedure that starts with an initial guess of $\|x^{1}-x^{*}\|$, and based on the duality gap (and other information), to search for an upper bound on $\|x^{1}-x^{*}\|$.

Next, we discuss how to handle unknown $\mu$. Note that even in the deterministic case, it is non-trivial to develop a modified version of our restart scheme (i.e., Algorithm 2D), which automatically adapts to $\mu$ and has the same oracle complexities up as Algorithm 2D (up to a log factor). In the context of strongly convex (smooth and non-smooth) minimization, many adaptive restart schemes have been proposed to avoid the knowledge of $\mu$ (see e.g., Nesterov [39], Juditsky and Nesterov [25], Lin and Xiao [30]* and Renegar and Grimmer [44]**). Among them, the restart scheme in [25] is most relevant to ours, since it also uses the quantity $R$ as the reduction criterion (as opposed to the primal sub-optimality gap; cf. Section 4.4.1). The basic idea is to notice that $\mu$ resides in a bounded interval $(0,B]$, and partition this interval as $(0,\mu_{1}]\cup(\mu_{1},\mu_{2}]\cup\ldots\cup(\mu_{K-1},B]$, where $K$ denotes the total number of stages (cf. Algorithm 2D), and*

[TABLE]

Furthermore, given the total number of iterations $T$, number of iterations in each stage is the same and equal to $\lfloor{T/K}\rfloor$ (i.e., no longer exponentially increasing as in Algorithm 2D). By choosing $K$ properly (and independent of $\mu$), they obtain a convergence rate $O(\log T/T)$ of the primal sub-optimality gap. However, a close inspection reveals that their approach cannot be easily extended to Algorithm 2D. Let $k^{*}\in\mathcal{K}$ satisfy $\mu\in(\mu_{k^{*}},\mu_{k^{*}+1}]$. One of the major difficulties is that in their scheme, for $k\geq k^{*}+1$, we no longer have the condition in (74), which is crucial to establish $\widehat{G}_{R_{k}}(\overline{x}_{k+1},y_{k+1})=G(\overline{x}_{k+1},y_{k+1})$. Therefore, in this case, we cannot derive the convergence rate of the duality gap. Due to this, we leave developing $\mu$-free versions of Algorithms 2D and 2S as open questions. However, since these difficulties are mainly related to theoretical (complexity) analysis, in practice, we can still use the aforementioned restart schemes as heuristics.

Finally, we discuss how to handle the unknown smoothness parameters, i.e., $L$, $L_{xx}$, $L_{xy}$ and $L_{yy}$. For simplicity, let us focus on unknown $L$ (as the other three parameters can be handled similarly). Indeed, in the deterministic case, Nesterov [39, Section 3]* has proposed a backtracking method for estimating $L$. This method starts from a lower bound on $L$, and increase it by a constant factor until certain termination criterion is satisfied. This method can be straightforwardly integrated into the subroutine of Algorithm 2D, i.e., Algorithm 1R. However, in the stochastic case, this backtracking method fails to work properly. Specifically, since the termination criterion now involves random variables, we need it to hold in expectation or with high probability. However, both of which are hard to check, since we have little knowledge about the distribution of random variables involved. For this reason, we leave estimating smoothness parameters in the stochastic case to future work.*

5 Application and numerical experiments.

We consider the problem of two-player zero-sum game with non-linear payoff in Hien et al. [21]* (see also Chen et al. [12]**). The problem has the following form:*

[TABLE]

where

[TABLE]

There are several problem constants in (116) and (117), including $\varpi>0$, $c_{i}>0$, for all $i\in[n]$, $Q\in\mathbb{S}_{++}^{n}$ (i.e., $Q\in\mathbb{R}^{n\times n}$ is a symmetric and positive definite matrix), $N>0$, and $P>0$.

Indeed, the problem in (116) corresponds to the water-filling problem in information theory (see e.g., Cover and Thomas [14]**). Specifically, let Alice and Bob be two players, both of whom face $n$ Gaussian communication channels, each with noise power $c_{i}$, for all $i\in[n]$. Alice has noise power $N$, and wishes to minimize the total channel capacity, by allocating noise power $x_{i}$ to channel $i$, for each $i\in[n]$. Bob has transmission power $P$, and wishes to maximize the total channel capacity, by allocating transmission power $y_{i}$ to channel $i$, for each $i\in[n]$. The total channel capacity is given by $\sum_{i=1}^{n}\log\left(1+{y_{i}}/({c_{i}+x_{i}})\right)$. In addition, for Alice, injecting noise power to the channels incurs cost $(1/{2})\langle{x},{Qx}\rangle$ (where $Q$ denotes the cost matrix), and $\varpi$ reflects the trade-off between minimizing her cost and minimizing the total channel capacity. In fact, by finding a saddle point of (116) (which always exists by Sion’s minimax theorem **[46]**), we find a Nash equilibrium of this two-player game.

Let us complete setting up the problem in (116), by specifying the normed spaces $\mathbb{X}$ and $\mathbb{Y}$. Indeed, if we set $\mathbb{X}=(\mathbb{R}^{n},\left\|\cdot\right\|_{2})$ and $\mathbb{Y}=(\mathbb{R}^{n},\left\|\cdot\right\|_{1})$, then in (116), by letting

[TABLE]

the objective function $S(x,y)\triangleq f(x)+\Phi(x,y)$ satisfies all the assumptions stated in Section 1, where the constants $L$, $\mu$, $L_{xx}$, $L_{xy}$ and $L_{yy}$ are given by

[TABLE]

In (119), $\lambda_{\max}(Q)>0$ and $\lambda_{\min}(Q)>0$ denote the largest and smallest eigenvalues of $Q$, respectively, and $c_{\min}\triangleq\min_{i\in[n]}c_{i}>0$. (For the derivation of the constants $L_{xx}$, $L_{xy}$ and $L_{yy}$, see Appendix 11.) In addition, under this setting, it is natural to choose the DGFs $h_{\mathcal{X}}$ and $h_{\mathcal{Y}}$ to be

[TABLE]

In what follows, we will the stochastic Mirror-Prox algorithm, i.e., Algorithm 1 in Juditsky et al. [24]**, as the benchmark. For ease of reference, we will refer to our stochastic restart scheme (i.e., Algorithm 2S) and Algorithm 1 in **[24]** as SRS and SMP, respectively.

5.1 Implementation details.

We conducted four sets of experiments, corresponding to different values of $n\in\{1000,2000,3000,4000\}$. For simplicity, in all our experiments, we set $\varpi=0.1$ and $c_{i}=1$, for all $i\in[n]$, in (116). In addition, in (117), we set $P=N=1$, so $\mathcal{X}=\mathcal{Y}=\Delta_{n}$, i.e., the $(n-1)$-dimensional probability simplex. For each value of $n$, we generated the cost matrix $Q$ as follows. First we generated a matrix $\overline{Q}\in\mathbb{R}^{n\times n}$ such that its entries $\{\overline{Q}_{ij}\}_{i,j=1}^{n}$ are i.i.d. standard normal. (Therefore, with probability one, $\overline{Q}$ has full rank.) Then we let $Q=\overline{Q}^{T}\overline{Q}$, and consequently,

[TABLE]

where $\overline{Q}_{i}$ denotes the $i$-th column of $\overline{Q}^{T}$. From (121) and (118), we observe that both $f$ and $\Phi$ have finite-sum forms. Therefore, to generate stochastic estimators of the gradients $\nabla f$, $\nabla_{x}\Phi(\cdot,\cdot)$ and $\nabla_{y}\Phi(\cdot,\cdot)$, we uniformly randomly drew a mini-batch $\mathcal{B}\subseteq[n]$ without replacement, and set

[TABLE]

where $e_{i}$ denotes the $i$-th standard coordinate vector. Based on these estimators, we can set the variance parameters in Assumption 3.2(B)* as*

[TABLE]

(For the derivation of these parameters, see Appendix 11.) For simplicity, we set $\left\lvert\mathcal{B}\right\rvert=\lceil{n/2}\rceil$ in all our experiments.

For each value of $n$, after generating $Q\in\mathbb{S}_{++}^{n}$ and setting the problem parameters, we were able to run both algorithms SRS and SMP. Specifically, for both algorithms, we used the same initial primal-dual pair $(x_{0},y_{0})$, which was randomly chosen from $\Delta_{n}\times\Delta_{n}$. As stated in Remark 4.6, we replaced all the absolute constants in the definition of $T_{k}$ in Algorithm 2S by a constant $c=0.7$. For either SRS or SMP, we stopped the algorithm once the duality gap fell below $\varepsilon=0.001$, and recorded the running time accordingly. For each value of $n$, we repeated this process five times, and calculated the mean and standard deviation (std) of the running times for each algorithm. All the codes were implemented in Python 3.7.6 on a machine with 1.90 GHz CPU and 16 GB RAM, and the evaluation of the duality gap was done by the disciplined convex programming solvers in the CVXPY package **[16]**.

5.2 Numerical results.

The results are shown in Table 3. From this table, we can observe that our algorithm, i.e., SRS, outperforms SMP for different values of $n$. Specifically, to reach a duality gap that is no more than $0.001$, SRS uses less time than SMP on average. This is consistent with the superiority of the oracle complexity of SRS over SMP (cf. Tables 1 and 2).

6 Proof of Proposition 3.1.

Before proving Proposition 3.1, we first present the Bregman proximal inequality associated with the Bregman proximal projection in (14) and its corollary. The proof of this inequality can be found in many works, e.g., Ghadimi and Lan [17, Lemma 2]**.

Lemma 6.1

In (14), for any $u\in\mathcal{U}$, we have

[TABLE]

Corollary 6.2

In (14), for any $u\in\mathcal{U}$, we have

[TABLE]

Proof 6.3

Proof. First, by Young’s inequality, for any $u\in\mathbb{U}$ and $u^{*}\in\mathbb{U}^{*}$, we have

[TABLE]

Therefore, $\langle{u^{*}},{u^{\prime}-u^{+}}\rangle\leq(2\lambda)^{-1}\|u^{+}-u^{\prime}\|^{2}+(\lambda/2)\|u^{*}\|_{*}^{2}$. Add up this inequality with (125), we then complete the proof. \Halmos

Proof 6.4

Proof of Proposition 3.1. For convenience, for any $t\in\mathbb{Z}^{+}$ and $(x,y)\in\mathcal{X}\times\mathcal{Y}$, let us define

[TABLE]

The most crucial step in our proof is to establish the recursion between $\widetilde{G}(\overline{x}^{t+1},\overline{y}^{t+1};x,y)$ and $\widetilde{G}(\overline{x}^{t},\overline{y}^{t};x,y)$, which requires us to establish the recursion for each of the four terms in (128). By the convexity of $g$ and $J$ and convexity-concavity of $\Phi$, we easily see that

[TABLE]

To connect $f(\overline{x}^{t+1})$ and $f(\overline{x}^{t})$, we have

[TABLE]

where in (a) we use the descent lemma (e.g., Bertsekas [7]), resulted from the $L$-smoothness of $f$; in (b) we use the step (22) in Algorithm 1; in (c) we use the convexity of $f$ and that $\overline{x}^{t+1}-\widetilde{x}^{t+1}=\beta_{t}(x^{t+1}-x^{t})$ (resulted from the steps (19) and (22)); in (d) we again use the convexity of $f$.

Combining (129), (130), (131) and (135), we have the following recursion

[TABLE]

We now can apply the Bregman proximal inequality in (125) to the steps (19) and (22), and obtain bounds on $g(x^{t+1})-g(x)$ and $J(y^{t+1})-J(y)$, i.e.,

[TABLE]

Note that when $t=1$, (137) holds for any $\hat{\nabla}_{y}\Phi(x^{0},y^{0},\zeta_{y}^{0})\in\mathbb{Y}^{*}$ since $\theta_{1}=0$.

To bound $\Phi(x^{t+1},y)-\Phi(x,y^{t+1})$, we have

[TABLE]

where in (142), we use the concavity of $\Phi(x^{t+1},\cdot)$, the convexity of $\Phi(\cdot,y^{t+1})$ and the $L_{xx}$-smoothness of $\Phi(\cdot,y^{t+1})$, respectively.

If we multiply both sides of (136) by $\gamma_{t}\beta_{t}^{-1}$, and then substitute (138), (140) and (143) into the resulting inequality, we have

[TABLE]

Before we proceed, recall that $\gamma_{t}\theta_{t}=\gamma_{t-1}$ and $\gamma_{t}\beta_{t}^{-1}=\gamma_{t+1}(\beta_{t+1}^{-1}-1)$, for any $t\in\mathbb{N}$. This enables us to observe certain recursion patterns (e.g., on $\gamma_{t}(\beta_{t}^{-1}-1)\widetilde{G}(\overline{x}^{t},\overline{y}^{t};x,y)$) in (144). We now bound the terms (I), (II) and (III) in (144).

To bound (I), we make use of Young’s inequality (cf. (127)), i.e.,

[TABLE]

where in (a) we use the Lipschitz continuity of ${\nabla}_{y}\Phi(x^{t},\cdot)$ and ${\nabla}_{y}\Phi(\cdot,y^{t-1})$ respectively and in (b) we use the conditions that $\gamma_{t}\theta_{t}=\gamma_{t-1}$ and $\alpha_{t}\theta_{t}\leq\alpha_{t-1}$ for any $t\in\mathbb{N}$.

To bound (II), we need to use a technique introduced in Nemirovski et al. [33], which involves the sequence $\{\hat{x}^{t}\}_{t\in\mathbb{N}}$ defined in (39). Based on $\{\hat{x}^{t}\}_{t\in\mathbb{N}}$, we decompose (II) into three parts, i.e.,

[TABLE]

To see the benefit of doing this, note that in (II.B), both $\hat{x}^{t},x^{t}\in\mathcal{F}_{t-1}$, i.e., $\hat{x}^{t}$ and $x^{t}$ are measurable with respect to $\mathcal{F}_{t-1}$. Therefore by Assumption 3.2(A), $\{\langle{\delta_{x,\Phi}^{t}+\delta_{x,f}^{t}},{\hat{x}^{t}-x^{t}}\rangle\}_{t\in\mathbb{Z}_{+}}$ is an MDS adapted to $\{\mathcal{F}_{t}\}_{t\in\mathbb{Z}_{+}}$. Moreover, (II.A) and (II.C) can also be bounded using Corollary 6.2 and Young’s inequality respectively, i.e.,

[TABLE]

In summary, we have

[TABLE]

We can bound (III) and (IV) in a similar fashion. Indeed, define $\hat{y}^{0}\triangleq y^{0}$ and for any $t\in\mathbb{N}$,

[TABLE]

We then have

[TABLE]

where in the last inequality we use the fact that $\theta_{t}\in[0,1]$. In addition,

[TABLE]

where in the last inequality we use $\|a+b\|^{2}\leq 2(\|a\|^{2}+\|b\|^{2})$, for any $a,b\in\mathbb{Y}$ and $\gamma_{t-1}/\alpha_{t-1}\leq\gamma_{t}/\alpha_{t}$ (since $\alpha_{t}\theta_{t}\leq\alpha_{t-1}$ and $\gamma_{t}\theta_{t}=\gamma_{t-1}$).

Now we can substitute (145), (149), (151) and (152) into the recursion (144) to obtain

[TABLE]

We then sum up the inequality (153) over $t=1,\ldots,T-1$ to obtain

[TABLE]

where we have used the facts that $\gamma_{t-1}/\alpha_{t-1}\leq\gamma_{t}/\alpha_{t}$.

In addition, we can bound (V) in a similar fashion to bounding (I) (cf. (145)), i.e.,

[TABLE]

where in the last inequality we have use $\alpha_{T-1}\leq(2L_{yy})^{-1}$ and $(1/2)\|y-y^{T}\|^{2}\leq D_{h_{\mathcal{Y}}}(y,y^{T})$.

We now substitute (155) into (154) and simplify the resulting inequality by noting that

$L\beta_{t}+L_{xx}-({2\tau_{t}})^{-1}+{4\alpha_{t}L_{yx}^{2}}\leq 0$, $(1+\theta_{t})L_{yy}-({8\alpha_{t}})^{-1}\leq 0$. 2. 2)

$D_{h_{\mathcal{X}}}(x,x^{T}),D_{h_{\mathcal{X}}}(x,\hat{x}^{T})\leq\Omega_{h_{\mathcal{X}}}$*, $\gamma_{t-1}/\alpha_{t-1}\leq\gamma_{t}/\alpha_{t}$, *

$(2+\theta_{t-1})\gamma_{t-1}/\alpha_{t-1}\leq(2+\theta_{t})\gamma_{t}/\alpha_{t}$* (since $\theta_{t-1}\leq\theta_{t}$), $D_{h_{\mathcal{Y}}}(y,y^{T}),D_{h_{\mathcal{Y}}}(y,\hat{y}^{T})\leq\Omega_{h_{\mathcal{Y}}}$.*

As a result, we have

[TABLE]

In addition, since $\gamma_{t-1}/\tau_{t-1}\leq\gamma_{t}/\tau_{t}$, we have

[TABLE]

We obtain the proof of parts (A) and (B) by treating the inequality (156) in different ways. For part (A), we simply take expectation on both sides of (156); whereas for part (B), we need to use concentration inequalities and the Chernoff bound. The details are shown below.

Proof of Part (A).* For any $t\in\mathbb{N}$, since $x^{t},\hat{x}^{t},y^{t},\hat{y}^{t}\in\mathcal{F}_{t-1}$, both $\{\langle{\delta_{x,\Phi}^{t}+\delta_{x,f}^{t}},{\hat{x}^{t}-x^{t}}\rangle\}_{t\in\mathbb{N}}$ and $\{\langle{\delta^{t}_{y,\Phi}},{y^{t}-\hat{y}^{t}}\rangle\}_{t\in\mathbb{N}}$ are MDS adapted to $\{\mathcal{F}_{t}\}_{t\in\mathbb{Z}_{+}}$. Therefore,*

[TABLE]

In addition, by Assumption 3.2(B), we have

[TABLE]

We then take expectation on both sides of (156) and substitute (157), (158), (159) and (160) to the resulting inequality to obtain (35). In addition, by Markov’s inequality, for any $\varsigma\in(0,1]$,

[TABLE]

where the last step follows from (35). This proves (38).

Proof of Part (B).* We first present the Asuma-Hoeffding lemma for sub-Gaussian MDS.*

Lemma 6.5** **(Nemirovski et al. [33])

Let $\{\epsilon_{t}\}_{t\in\mathbb{N}}$ be a real-valued MDS adapted to a filtration $\{\mathcal{F}_{t}\}_{t\in\mathbb{Z}_{+}}$, such that for any $t\in\mathbb{N}$, $\mathbb{E}_{t-1}[\epsilon_{t}]=0$ and there exists a constant $d_{t}>0$ such that $\mathbb{E}_{t-1}[\epsilon_{t}^{2}/d_{t}^{2}]\leq\exp(1)$. Then for any $p>0$ and $T\in\mathbb{N}$,

[TABLE]

For convenience, define $C_{T}\triangleq\sum_{t=1}^{T-1}\gamma_{t}\alpha_{t}$ and $C^{\prime}_{T}\triangleq\sum_{t=1}^{T-1}\gamma_{t}\tau_{t}$. Then

[TABLE]

where in (a) we use Markov’s inequality, in (b) we use the convexity of $\exp(\cdot)$ and in (c) we use Assumption 3.2(C). Similarly, we can also show that for any $\widetilde{p}>0$,

[TABLE]

*Next, since $\{\gamma_{t}\langle{\delta_{x,f}^{t}},{x^{t}-\hat{x}^{t}}\rangle\}_{t\in\mathbb{N}}$, $\{\gamma_{t}\langle{\delta_{x,\Phi}^{t}},{x^{t}-\hat{x}^{t}}\rangle\}_{t\in\mathbb{N}}$ and $\{\gamma_{t-1}\langle{\delta_{y,\Phi}^{t}},{\hat{y}^{t}-y^{t}}\rangle\}_{t\in\mathbb{N}}$ are MDS, we can apply Lemma 6.5 to the last three terms in (156). Specifically, let us define *

[TABLE]

Then by Assumption 3.2(C), for any $t=1,\ldots,T-2$, we have

[TABLE]

and

[TABLE]

Similarly, for any $t=1,\ldots,T-1$, we have

[TABLE]

Thus by Lemma 6.5, for any $q,\widetilde{q}>0$, we have

[TABLE]

We then combine (156), (162), (163), (164), (170) and (171), and take $p=\widetilde{p}={\log(1/\varsigma)}$ and $q=\widetilde{q}=2\sqrt{\log(1/\varsigma)}$ to obtain (40). To obtain (42), we simply substitute (157) into (40). \Halmos

{APPENDICES}

7 Unique output of the Bregman proximal projection (BPP).

Lemma 7.1

For the minimization problem in (14), if $\varphi^{*}\triangleq\inf_{u\in\mathcal{U}}\varphi(u)>-\infty$ and $\mathcal{U}^{o}\cap\mathsf{dom}\,\varphi\neq\emptyset$, then a unique solution exists in $\mathcal{U}^{o}\cap\mathsf{dom}\,\varphi$.

Proof 7.2

Proof. By condition (12) and $\varphi^{*}>-\infty$, we see that $\varphi_{\lambda}$ is coercive on $\mathbb{U}$, i.e., $\lim_{\left\|u\right\|\to+\infty}\varphi_{\lambda}(u)=+\infty$. In addition, $\varphi_{\lambda}$ is CCP since both $\varphi$ and $h_{\mathcal{U}}$ are CCP. Choose any point $u\in\mathcal{U}\cap\mathsf{dom}\,\varphi_{\lambda}\neq\emptyset$ and any $\alpha\geq\varphi_{\lambda}(u)$. The closedness and coercivity of $\varphi_{\lambda}$ together imply that the sub-level set $\mathcal{S}_{\alpha}(\varphi_{\lambda})\triangleq\{u\in\mathbb{U}:\varphi_{\lambda}(u)\leq\alpha\}$ is compact. Since $\mathcal{U}$ is closed, $\mathcal{U}\cap\mathcal{S}_{\alpha}(\varphi_{\lambda})$ is compact and nonempty. This, together with the lower semi-continuity of $\varphi_{\lambda}$, implies that the solution set of Problem (14), denoted by $\mathcal{U}_{\rm opt}$, is nonempty and contained in $\mathcal{U}\cap\mathsf{dom}\,\varphi_{\lambda}$. Since $\mathsf{int}\,\mathsf{dom}\,h_{\mathcal{U}}\cap\mathsf{dom}\,\varphi\neq\emptyset$, by Moreau-Rockafellar theorem (see e.g., Peypouquet [42, Theorem 3.30]), $\partial\varphi_{\lambda}(u)=\partial\varphi(u)+\partial h_{\mathcal{U}}(u)+u^{*}-\nabla h(u^{\prime})$, for any $u\in\mathcal{U}$. Since $h_{\mathcal{U}}$ is essentially smooth, for any $u\in\mathsf{dom}\,h_{\mathcal{U}}\setminus\mathsf{int}\,\mathsf{dom}\,h_{\mathcal{U}}$, $\partial h_{\mathcal{U}}(u)=\emptyset$, and hence $\partial\varphi_{\lambda}(u)=\emptyset$. As a result, $\mathcal{U}_{\rm opt}\subseteq\mathcal{U}^{o}\cap\mathsf{dom}\,\varphi$. Finally, since $h_{\mathcal{U}}$ is 1-strongly convex on $\mathcal{U}^{o}$, so is $\varphi_{\lambda}$, and therefore, $\mathcal{U}_{\rm opt}$ must be a singleton. \Halmos

8 Berge’s maximum theorem.

The following theorem is adapted from Ok [40, Section E.3], which was originally proved by Berge [6].

Theorem 8.1

Let $\mathbb{U}$ and $\mathbb{V}$ be two metric spaces and $\Psi:\mathbb{U}\times\mathbb{V}\to\mathbb{R}$ be a continuous function on $\mathbb{U}\times\mathbb{V}$. If $\mathcal{V}\subseteq\mathbb{V}$ is nonempty and compact, then $\psi(u)\triangleq\sup_{v\in\mathcal{V}}\Psi(u,v)$ is continuous on $\mathbb{U}$.

9 Stopping criterion via gradient mappings.

Consider the BPP in (14) with $u^{*}=\nabla f(u^{\prime})$, where the function $f$ is $L$-smooth and $\mu$-strongly convex on $\mathcal{U}$. Accordingly, let us define two gradient mappings, i.e.,

[TABLE]

As we will show in the lemma below (i.e., Lemma 9.1), the norm of the two gradient mappings, i.e., $G_{\lambda}(u^{\prime})$ and $\overline{G}_{\lambda}(u^{\prime})$, can be used as the stopping criterion for the following optimization problem:

[TABLE]

Clearly, if $\mathsf{dom}\,\varphi\cap\mathcal{U}\neq\emptyset$, then (173) has a unique optimal solution, which we denote by $u_{\rm opt}\in\mathcal{U}$. Also, note that if $u_{\rm opt}\in\mathcal{U}^{o}$, then $G_{\lambda}(u_{\rm opt})=\overline{G}_{\lambda}(u_{\rm opt})=0$, since

[TABLE]

where in the last step, we use the uniqueness of the BPP (cf. Lemma 7.1).

Lemma 9.1

Consider the setting in Appendix 9. For any $\epsilon>0$ and $0<\lambda\leq 1/L$, we have

[TABLE]

Proof 9.2

Proof. By the first-order optimality condition of (14), we have

[TABLE]

where $\mathcal{N}_{\mathcal{U}}(u^{+})$ denotes the normal cone of $\mathcal{U}$ at $u^{+}\in\mathcal{U}$. Consequently, we have

[TABLE]

The $\mu$-strong convexity of $P$ on $\mathcal{U}$ indicates that

[TABLE]

On the other hand, we have

[TABLE]

where in the last step we use that $0<\lambda\leq 1/L$. Now, combining (177) and (178), we have

[TABLE]

We hence complete the proof. \Halmos

Remark 9.3

Two remarks are in order. First, note that both $G_{\lambda}(u^{\prime})$ and $\overline{G}_{\lambda}(u^{\prime})$ can be easily computed, since 1) we assume that the BPP in (14) has an easily computable solution (cf. Section 2.1), and 2) the DGF $h_{\mathcal{U}}$ typically has simple structure, so that $\nabla h_{\mathcal{U}}$ can be obtained easily. Second, if $\nabla h_{\mathcal{U}}$ is $\eta$-Lipschitz on $\mathcal{U}$ ($\eta\geq 1$), as in the Hilbertian case where $h_{\mathcal{U}}=(1/2)\|\cdot\|^{2}$, then $\|\overline{G}_{\lambda}(u^{\prime})\|_{*}\leq\eta\|G_{\lambda}(u^{\prime})\|$ and we need not compute $\nabla h_{\mathcal{U}}$.

10 Lipschitz continuity of $x^{*}(\cdot)$.

The following lemma is similar to Zhao [52, Lemma 3.4]**. For completeness, we provide it proof here.

Lemma 10.1

Let $x^{*}(\cdot):\mathcal{Y}\to\mathbb{X}$ be defined in (111). Then it is $(L_{yx}/\mu)$-Lipschitz on $\mathcal{Y}$.

Proof 10.2

Proof. By the $\mu$-strong convexity of $f$ (hence $\psi^{\mathrm{P}}(\cdot,y)$, for any $y\in\mathcal{Y}$), we have

[TABLE]

On the other hand,

[TABLE]

Combining (180) to (183), we have

[TABLE]

Note that (184) implies that $\left\|x^{*}(y^{\prime})-x^{*}(y)\right\|\leq(L_{yx}/\mu)\|y-y^{\prime}\|$, for both $x^{*}(y^{\prime})=x^{*}(y)$ and $x^{*}(y^{\prime})\neq x^{*}(y)$. We hence complete the proof. \Halmos

11 Technical proofs in Section 5.

We first derive the constants $L_{xx}$, $L_{xy}$ and $L_{yy}$ in (119). Let $i,j\in[n]$ and $\delta_{ij}$ be the Kronecker delta, i.e., $\delta_{ij}=1$ if $i=j$ and [math] otherwise. It is easy to see that

[TABLE]

In addition, we have that for any $i\in[n]$,

[TABLE]

where in (a) we use $c_{i}+x_{i}\geq c_{\min}$ and that the function $a\mapsto a^{-2}-(a+1)^{-2}$ is strictly decreasing on $(0,+\infty)$. Therefore, we can let $L_{xx}=c_{\min}^{-2}-(c_{\min}+1)^{-2}$ and $L_{xy}=L_{yy}=c_{\min}^{-2}$.

Next, we derive the variance parameters in (124). To do that, let us first present a well-known lemma on the variance of sampling without replacement from a finite population (see e.g., Banerjee [2]**).

Lemma 11.1

Given a finite set $\mathcal{Z}\triangleq\{z_{i}\}_{i=1}^{n}\subseteq\mathbb{R}^{d}$, uniformly randomly draw a subset $\mathcal{B}\subseteq[n]$ without replacement. Define the population mean $\overline{z}\triangleq n^{-1}\sum_{i=1}^{n}z_{i}$. Then

[TABLE]

Based on Lemma 11.1, to derive $\sigma_{x,f}^{2}$, it suffices to upper bound $\left\|z_{i}\right\|_{2}^{2}$ and lower bound $\left\|\overline{z}\right\|_{2}^{2}$, where $z_{i}=\overline{Q}_{i}\overline{Q}_{i}^{T}x$, and $x\in\Delta_{n}$. By Hölder’s inequality, it is clear that $\left\|z_{i}\right\|_{2}^{2}=\left\|\overline{Q}_{i}\right\|_{2}^{2}(\overline{Q}_{i}^{T}x)^{2}\leq\left\|\overline{Q}_{i}\right\|_{2}^{2}\left\|\overline{Q}_{i}\right\|_{\infty}^{2}$. In addition, since $\overline{z}=n^{-1}\sum_{i=1}^{n}\overline{Q}_{i}\overline{Q}_{i}^{T}x=n^{-1}Qx$, we have $\left\|\overline{z}\right\|_{2}^{2}=n^{-2}\left\|Qx\right\|_{2}^{2}\geq n^{-2}\left\|x\right\|_{2}^{2}\sigma^{2}_{\min}(Q)\geq n^{-3}\sigma_{\min}^{2}(Q)$. As a result, for any $(x,y)\in\Delta_{n}\times\Delta_{n}$,

[TABLE]

Using the same argument, we have

[TABLE]

Acknowledgments.

The author would like to thank Robert Freund for his constructive feedback and helpful discussions during the preparation of this manuscript. The author would also like to thank two anonymous referees for their helpful and insightful suggestions, which have greatly improved the exposition of the paper. The author’s research is supported by AFOSR Grant No. FA9550-19-1-0240.