Stochastic subgradient projection methods for composite optimization with functional constraints

The large sum of functions in the objective function and/or the large number of constraints in most of the practical optimization applications, including machine learning and statistics (Vapnik,, 1998; Bhattacharyya et al.,, 2004), signal processing (Necoara,, 2021; Tibshirani,, 2011), computer science (Kundu et al.,, 2018), distributed control (Nedelcu et al.,, 2014), operations research and finance (Rockafellar and Uryasev,, 2000), create several theoretical and computational challenges. For example, these problems are becoming increasingly large in terms of both the number of variables and the size of training data and they are usually nonsmooth due to the use of regularizers, penalty terms and the presence of constraints. Due to these challenges, (sub)gradient-based methods are widely applied. In particular, proximal gradient algorithms (Necoara,, 2021; Rosasco et al.,, 2019) are natural in applications where the function to be minimized is in composite form, i.e. the sum of a smooth term and a nonsmooth regularizer (e.g., the indicator function of some simple constraints), as e.g. the empirical risk in machine learning. In these algorithms, at each iteration the proximal operator defined by the nonsmooth component is applied to the gradient descent step for the smooth data fitting component. In practice, it is important to consider situations where these operations cannot be performed exactly. For example, the case where the gradient, the proximal operator or the projection can be computed only up to an error have been considered in (Devolder et al.,, 2014; Nedelcu et al.,, 2014; Necoara and Patrascu,, 2018; Rasch and Chambolle,, 2020), while the situation when only stochastic estimates of these operators are available have been studied in (Rosasco et al.,, 2019; Hardt et al.,, 2016; Patrascu and Necoara,, 2018; Hermer et al.,, 2020). This latter setting, which is also of interest here, is very important in machine learning, where we have to minimize an expected objective function with or without constraints from random samples (Bhattacharyya et al.,, 2004), or in statistics, where we need to minimize a finite sum objective subject to functional constraints (Tibshirani,, 2011).

Previous work. A very popular approach for minimizing an expected or finite sum objective function is the stochastic gradient descent (SGD) (Robbins and Monro,, 1951; Nemirovski and Yudin,, 1983; Hardt et al.,, 2016) or the stochastic proximal point (SPP) algorithms (Moulines and Bach,, 2011; Nemirovski et al.,, 2009; Necoara,, 2021; Patrascu and Necoara,, 2018; Rosasco et al.,, 2019). In these studies sublinear convergence is derived for SGD or SPP with decreasing stepsizes under the assumptions that the objective function is smooth and (strongly) convex. For nonsmooth stochastic convex optimization one can recognize two main approaches. The first one uses stochastic variants of the subgradient method combined with different averaging techniques, see e.g. (Nemirovski et al.,, 2009; Yang and Lin,, 2016). The second line of research is based on stochastic variants of the proximal gradient method under the assumption that the expected objective function is in composite form (Duchi and Singer,, 2009; Necoara,, 2021; Rosasco et al.,, 2019). For both approaches, using decreasing stepsizes, sublinear convergence rates are derived under the assumptions that the objective function is (strongly) convex. Even though SGD and SPP are well-developed methodologies, they only apply to problems with simple constraints, requiring the whole feasible set to be projectable.

In spite of its wide applicability, the study on efficient solution methods for optimization problems with many constraints is still limited. The most prominent line of research in this area is the alternating projections, which focus on applying random projections for solving problems that involve intersection of a (infinite) number of sets. The case when the objective function is not present in the formulation, which corresponds to the convex feasibility problem, stochastic alternating projection algorithms were analyzed e.g., in (Bauschke and Borwein,, 1996), with linear convergence rates, provided that the sets satisfy some linear regularity condition. Stochastic forward-backward algorithms have been also applied to solve optimization problems with many constraints. However, the papers introducing those general algorithms focus on proving only asymptotic convergence results without rates, or they assume the number of constraints is finite, which is more restricted than our settings, see e.g. (Bianchi et al.,, 2019; Xu,, 2020; Wang et al.,, 2015). In the case where the number of constraints is finite and the objective function is deterministic, Nesterov’s smoothing framework is studied in (Tran-Dinh et al.,, 2018) under the setting of accelerated proximal gradient methods. Incremental subgradient methods or primal-dual approaches were also proposed for solving problems with finite intersection of simple sets through an exact penalty reformulation in  (Bertsekas,, 2011; Kundu et al.,, 2018).

The papers most related to our work are (Nedich,, 2011; Nedich and Necoara,, 2019), where subgradient methods with random feasibility steps are proposed for solving convex problems with deterministic objective and many functional constraints. However, the optimization problem, the algorithm and consequently the convergence analysis are different from the present paper. In particular, our algorithm is a stochastic proximal gradient extension of the algorithms proposed in (Nedich,, 2011; Nedich and Necoara,, 2019). Additionally, the stepsizes in (Nedich,, 2011; Nedich and Necoara,, 2019) are chosen decreasing, while in the present work for strongly like objective functions we derive insightful stepsize-switching rules which describe when one should switch from a constant to a decreasing stepsize regime. Furthermore, in (Nedich,, 2011) and (Nedich and Necoara,, 2019) sublinear convergence rates are established either for convex or strongly convex deterministic objective functions, respectively, while in this paper we prove (sub)linear rates under an expected composite objective function which is either convex or strongly convex. Moreover, (Nedich,, 2011; Nedich and Necoara,, 2019) present separately the convergence analysis for smooth and nonsmooth objective, while in this paper we present a unified convergence analysis covering both cases through the so-called stochastic bounded gradient condition. Hence, since we deal with stochastic composite objective functions, smooth or nonsmooth, convex or strongly convex, and since we consider a stochastic proximal gradient with new stepsize rules, our convergence analysis requires additional insights that differ from that of (Nedich,, 2011; Nedich and Necoara,, 2019).

In (Patrascu and Necoara,, 2018) a stochastic optimization problem with infinite intersection of sets is considered and stochastic proximal point steps are combined with alternating projections for solving it. However, in order to prove sublinear convergence rates, (Patrascu and Necoara,, 2018) requires strongly convex and smooth objective functions, while our results are valid also for convex non-smooth functions. Lastly, (Patrascu and Necoara,, 2018) assumes the projectability of individual sets, whereas in our case, the constraints might not be projectable. Finally, in all these studies the nonsmooth component is assumed to be proximal, i.e. one can easily compute its proximal operator. This assumption is restrictive, since in many applications the nonsmooth term is also expressed as expectation or finite sum of nonsmooth functions, which individually are proximal, but it is hard to compute the proximal operator for the whole nonsmooth component, as e.g., in support vector machine or generalized lasso problems (Villa et al.,, 2014; Rosasco et al.,, 2019). Moreover, all the convergence results from the previous papers are derived separately for the smooth and the nonsmooth stochastic optimization problems.

Contributions. In this paper we remove the previous drawbacks. We propose a stochastic subgradient algorithm for solving general convex optimization problems having the objective function expressed in terms of expectation operator over a sum of two terms subject to (possibly infinite number of) functional constraints. The only assumption we require is to have access to an unbiased estimate of the (sub)gradient and of the proximal operator of each of these two terms and to the subgradients of the constraints. To deal with such problems, we propose a stochastic subgradient method with random feasibility updates. At each iteration, the algorithm takes a stochastic subgradient step aimed at only minimizing the expected objective function, followed by a feasibility step for minimizing the feasibility violation of the observed random constraint achieved through Polyak’s subgradient iteration, see (Polyak,, 1969). In doing so, we can avoid the need for projections onto the whole set of constraints, which may be expensive computationally. The proposed algorithm is applicable to the situation where the whole objective function and/or constraint set are not known in advance, but they are rather learned in time through observations.

We present a general framework for the convergence analysis of this stochastic subgradient algorithm which is based on the assumptions that the objective function satisfies a stochastic bounded gradient condition, with or without strong convexity and the subgradients of the functional constraints are bounded. These assumptions include the most well-known classes of objective functions and of constraints analyzed in the literature: composition of a nonsmooth function and a smooth function, with or without strong convexity, and nonsmooth Lipschitz functions, respectively. Moreover, when the objective function satisfies the strong convexity property, we prove insightful stepsize-switching rules which describe when one should switch from a constant to a decreasing stepsize regime. Then, we prove sublinear convergence rates for the weighted averages of the iterates in terms of expected distance to the constraint set, as well as for the expected optimality of the function values/distance to the optimal set. Under some special conditions we also derive linear rates. Our convergence estimates are known to be optimal for this class of stochastic subgradient schemes for solving (nonsmooth) convex problems with functional constraints. Besides providing a general framework for the design and analysis of stochastic subgradient methods, in special cases, where complexity bounds are known for some particular problems, our convergence results recover the existing bounds. In particular, for problems without functional constraints we recover the complexity bounds from (Necoara,, 2021) and for linearly constrained least-squares problems we get similar convergence bounds as in (Leventhal and Lewis,, 2010).

Content. In Section 2 we present our optimization problem and the main assumptions. In Section 3 we design a stochastic subgradient projection algorithm and analyze its convergence properties, while is Section 4 we adapt this algorithm to constrained least-squares problems and derive linear convergence. Finally, in Section 5 detailed numerical simulations are provided that support the effectiveness of our method in real problems.

Recall that the proximal operator of the indicator function of a closed convex set becomes the projection. We consider additionally the following assumptions on the objective function and constraints:

From Jensen’s inequality, taking $x=x^{*}\in\mathcal{X}^{*}$ in (2), we get:

$$B^{2}\geq\mathbb{E}[\|\nabla F(x^{*},\zeta)\|^{2}]\geq\|\mathbb{E}[\nabla F(x^{*},\zeta)]\|^{2}=\|\nabla F(x^{*})\|^{2}\quad\forall x^{*}\in\mathcal{X}^{*}.$$

(3)

We also assume $F$ to satisfy a (strong) convexity condition:

Note that when $\mu=0$ relation (4) states that $F$ is convex on $\mathcal{Y}$. Finally, we assume boundedness on the subgradients of the functional constraints:

Note that our assumptions are quite general and cover the most well-known classes of functions analyzed in the literature. In particular, Assumptions 1 and 2, related to the objective function, covers the class of non-smooth Lipschitz functions and composition of a (potentially) non-smooth function and a smooth function, with or without strong convexity, as the following examples show.

Example 1 [Non-smooth (Lipschitz) functions satisfy Assumption 1]: Assume that the functions $f$ and $g$ have bounded (sub)gradients:

$$\|\nabla f(x,\zeta)\|\leq B_{f}\quad\text{and}\quad\|\nabla g(x,\zeta)\|\leq B_{g}\quad\forall x\in\mathcal{Y}.$$

Then, obviously Assumption 1 holds with $L=0\quad\text{and}\quad B^{2}=2B_{f}^{2}+2B_{g}^{2}.$

Example 2 [Smooth (Lipschitz gradient) functions satisfy Assumption 1]: Condition (2) includes the class of functions formed as a sum of two terms, $f(\cdot,\zeta)$ convex having Lipschitz continuous gradient and $g(\cdot,\zeta)$ convex having bounded subgradients, and $\mathcal{Y}$ bounded. Indeed, let us assume that the convex function $f(\cdot,\zeta)$ has Lipschitz continuous gradient, i.e. there exists $L_{f}(\zeta)>0$ such that:

$$\|\nabla f(x,\zeta)-\nabla f(\bar{x},\zeta)\|\leq L_{f}(\zeta)\|x-\bar{x}\|\quad\forall x\in\mathcal{Y}.$$

Using standard arguments (see Theorem 2.1.5 in (Nesterov,, 2018)), we have:

$$f(x,\zeta)-f(\bar{x},\zeta)\geq\langle\nabla f(\bar{x},\zeta),x-\bar{x}\rangle+\frac{1}{2L_{f}(\zeta)}\|\nabla f(x,\zeta)-\nabla f(\bar{x},\zeta)\|^{2}.$$

Since $g(\cdot,\zeta)$ is also convex, then adding $g(x,\zeta)-g(\bar{x},\zeta)\geq\langle\nabla g(\bar{x},\zeta),x-\bar{x}\rangle$ in the previous inequality, where $\nabla g(\bar{x},\zeta)\in\partial g(\bar{x},\zeta)$, we get:

$\displaystyle F(x,\zeta)-F(\bar{x},\zeta)$

$\displaystyle\geq\langle\nabla F(\bar{x},\zeta),x-\bar{x}\rangle+\frac{1}{2L_{f}(\zeta)}\|\nabla f(x,\zeta)-\nabla f(\bar{x},\zeta)\|^{2},$

where we used that $\nabla F(\bar{x},\zeta)=\nabla f(\bar{x},\zeta)+\nabla g(\bar{x},\zeta)\in\partial F(\bar{x},\zeta)$. Taking expectation w.r.t. $\zeta$ and assuming that the set $\mathcal{Y}$ is bounded, with the diameter $D$, then after using Cauchy-Schwartz inequality, we get (we also assume $L_{f}(\zeta)\leq L_{f}$ for all $\zeta$):

$\displaystyle F(x)-F^{*}$

$\displaystyle\geq\frac{1}{2L_{f}}\mathbb{E}[\|\nabla f(x,\zeta)-\nabla f(\bar{x},\zeta)\|^{2}]-D\|\nabla F(\bar{x})\|\quad\forall x\in\mathcal{Y}.$

Therefore, for any $\nabla g(x,\zeta)\in\partial g(x,\zeta)$, we have:

	$\displaystyle\mathbb{E}[\|\nabla F(x,\zeta)\|^{2}]$	$\displaystyle=\mathbb{E}[\|\nabla f(x,\zeta)-\nabla f(\bar{x},\zeta)+\nabla g(x,\zeta)+\nabla f(\bar{x},\zeta)\|^{2}]$
		$\displaystyle\leq 2\mathbb{E}[\|\nabla f(x,\zeta)-\nabla f(\bar{x},\zeta)\|^{2}]+2\mathbb{E}[\|\nabla g(x,\zeta)+\nabla f(\bar{x},\zeta)\|^{2}]$
		$\displaystyle\leq 4L_{f}(F(x)-F^{*})+4L_{f}D\|\nabla F(\bar{x})\|+2\mathbb{E}[\|\nabla g(x,\zeta)+\nabla f(\bar{x},\zeta)\|^{2}].$

Assuming now that the regularization functions $g$ have bounded subgradients on $\mathcal{Y}$, i.e., $\|\nabla g(x,\zeta)\|\leq B_{g}$, then the bounded gradient condition (2) holds on on $\mathcal{Y}$ with:

$$L=4L_{f}\quad\text{and}\quad B^{2}=4\left(B_{g}^{2}+\max_{\bar{x}\in\mathcal{X}^{*}}\left(\mathbb{E}[\|\nabla f(\bar{x},\zeta)\|^{2}]+DL_{f}\|\nabla F(\bar{x})\|\right)\right).$$

Note that $B^{2}$ is finite, since the optimal set $\mathcal{X}^{*}$ is compact (recall that in this example $\mathcal{Y}$ is assumed bounded).

Finally, we also impose some regularity condition for the constraints.

Note that this assumption holds e.g., when the index set $\Omega_{2}$ is arbitrary and the feasible set $\mathcal{X}$ has an interior point, see (Polyak,, 2001), or when the feasible set is polyhedral, see relation (24) below. However, Assumption (4) holds for more general sets, e.g., when a strengthened Slater condition holds for the collection of convex functional constraints, such as the generalized Robinson condition, as detailed in Corollary 3 of (Lewis and Pang,, 1998).

Here, $\alpha_{k}>0$ and $\beta>0$ are deterministic stepsizes and $(x)_{+}=\max\{0,x\}$. Note that $v_{k}$ represents a stochastic proximal (sub)gradient step (or a stochastic forward-backward step) at $x_{k}$ for the expected composite objective function $F(x)=\mathbb{E}[f(x,\zeta)+g(x,\zeta)]$. Note that the optimality step selects a random pair of functions $(f(x_{k},\zeta_{k}),g(x_{k},\zeta_{k}))$ from the composite objective function $F$ according to the probability distribution $\textbf{P}_{1}$, i.e., the index variable $\zeta_{k}$ with values in the set $\Omega_{1}$. Also, we note that the random feasibility step selects a random constraint $h(\cdot,\xi_{k})\leq 0$ from the collection of constraints set according to the probability distribution $\textbf{P}_{2}$, independently from $\zeta_{k}$, i.e. the index variable $\xi_{k}$ with values in the set $\Omega_{2}$. The vector $\nabla h(v_{k},\xi_{k})$ is chosen as:

$$\nabla h(v_{k},\xi_{k})=\begin{cases}\nabla h(v_{k},\xi_{k})\in\partial((h(v_{k},\xi_{k}))_{+})&\mbox{if }(h(v_{k},\xi_{k}))_{+}>0\\ s_{h}\neq 0&\mbox{if }(h(v_{k},\xi_{k}))_{+}=0,\end{cases}$$

where $s_{h}\in\mathbb{R}^{n}$ is any nonzero vector. If $(h(v_{k},\xi_{k}))_{+}=0$, then for any choice of nonzero $s_{h}$, we have $z_{k}=v_{k}$. Note that the feasibility step (7) has the special form of the Polyak’s subgradient iteration, see e.g., (Polyak,, 1969). Moreover, when $\beta=1$, $z_{k}$ is the projection of $v_{k}$ onto the hyperplane:

$$\mathcal{H}_{v_{k},\xi_{k}}=\{z:\;h(v_{k},\xi_{k})+\nabla h(v_{k},\xi_{k})^{T}(z-v_{k})\leq 0\},$$

that is $z_{k}=\Pi_{\mathcal{H}_{v_{k},\xi_{k}}}(v_{k})$ for $\beta=1$. Indeed, if $v_{k}\in\mathcal{H}_{v_{k},\xi_{k}}$, then $(h(v_{k},\xi_{k}))_{+}=0$ and the projection is the point itself, i.e., $z_{k}=v_{k}$. On the other hand, if $v_{k}\notin\mathcal{H}_{v_{k},\xi_{k}}$, then the projection of $v_{k}$ onto $\mathcal{H}_{v_{k},\xi_{k}}$ reduces to the projection onto the corresponding hyperplane:

$\displaystyle z_{k}=v_{k}-\frac{h(v_{k},\xi_{k})}{\|\nabla h(v_{k},\xi_{k})\|^{2}}\nabla h(v_{k},\xi_{k}).$

Combining these two cases, we finally get our update (7) . Note that when the feasible set of optimization problem (1) is described by (infinite) intersection of simple convex sets, see e.g. (Patrascu and Necoara,, 2018; Bianchi et al.,, 2019; Xu,, 2020; Wang et al.,, 2015):

$$\mathcal{X}=\mathcal{Y}\cap\left(\cap_{\xi\in\Omega_{2}}\mathcal{X}_{\xi}\right),$$

where each set $X_{\xi}$ admits an easy projection, then one can choose the following functional representation in problem (1):

$$h(x,\xi)=(h(x,\xi))_{+}={\rm dist}(x,\mathcal{X}_{\xi})\quad\forall\xi\in\Omega_{2}.$$

One can easily notice that this function is convex, nonsmooth and with bounded subgradients, since we have (Mordukhovich and Nam,, 2005):

$$\frac{x-\Pi_{\mathcal{X}_{\xi}}(x)}{\|x-\Pi_{\mathcal{X}_{\xi}}(x)\|}\in\partial h(x,\xi).$$

In this case, step (7) in SSP algorithm becomes a usual random projection step:

$$z_{k}=v_{k}-\beta(v_{k}-\Pi_{\mathcal{X}_{\xi_{k}}}(v_{k})).$$

Hence, our formulation is more general than (Patrascu and Necoara,, 2018; Bianchi et al.,, 2019; Xu,, 2020; Wang et al.,, 2015), as it allows to also deal with constraints that might not be projectable, but one can easily compute a subgradient of $h$. We mention also the possibility of performing iterations in parallel in steps (6) and (7) of SSP algorithm. However, here we do not consider this case. A thorough convergence analysis of minibatch iterations when the objective function is deterministic and strongly convex can be found in (Nedich and Necoara,, 2019). For the analysis of SSP algorithm, let us define the stochastic (sub)gradient mapping (for simplicity we omit its dependence on stepsize $\alpha$):

$$\mathcal{S}(x,\zeta)=\alpha^{-1}\left(x-\operatorname{prox}_{\alpha g(\cdot,\zeta)}(x-\alpha{\nabla}f(x,\zeta))\right).$$

Then, it follows immediately that the stochastic proximal (sub)gradient step (6) can be written as:

$${\color[rgb]{0,0,0}u_{k}}=x_{k}-\alpha_{k}\mathcal{S}(x_{k},\zeta_{k}).$$

Moreover, from the optimality condition of the prox operator it follows that there exists $\nabla g({\color[rgb]{0,0,0}u_{k}},\zeta_{k})\in\partial g({\color[rgb]{0,0,0}u_{k}},\zeta_{k})$ such that:

$$\mathcal{S}(x_{k},\zeta_{k})=\nabla f(x_{k},\zeta_{k})+\nabla g({\color[rgb]{0,0,0}u_{k}},\zeta_{k}).$$

Let us also recall a basic property of the projection onto a closed convex set $\mathcal{X}\subseteq\mathbb{R}^{n}$, see e.g., (Nedich and Necoara,, 2019):

$\displaystyle\|\Pi_{\mathcal{X}}(v)-y\|^{2}\leq\|v-y\|^{2}-\|\Pi_{\mathcal{X}}(v)-v\|^{2}\qquad\forall v\in\mathbb{R}^{n}\;\text{and}\;y\in\mathcal{X}.$

(9)

Define also the filtration:

$$\mathcal{F}_{[k]}=\{\zeta_{0},\cdots,\zeta_{k},\;\;\xi_{0},\cdots,\xi_{k}\}.$$

The next lemma provides a key descent property for the sequence $v_{k}$ and for the proof we use as main tool the stochastic (sub)gradient mapping $\mathcal{S}(\cdot)$.

Proof  Recalling that $\bar{v}_{k}=\Pi_{\mathcal{X}^{*}}(v_{k})$ and $\bar{x}_{k}=\Pi_{\mathcal{X}^{*}}(x_{k})$, from the definition of ${\color[rgb]{0,0,0}u_{k}}$ and using (9) for $y=\bar{x}_{k}\in\mathcal{X^{*}}{\color[rgb]{0,0,0}\subseteq\mathcal{Y}}$ and $v=v_{k}$, we get:

	$\displaystyle\|v_{k}-\bar{v}_{k}\|^{2}$	$\displaystyle\overset{\eqref{proj_prop}}{\leq}\|v_{k}-\bar{x}_{k}\|^{2}{\color[rgb]{0,0,0}=\|\Pi_{\mathcal{Y}}(u_{k})-\Pi_{\mathcal{Y}}(\bar{x}_{k})\|^{2}\leq\|u_{k}-\bar{x}_{k}\|^{2}}=\|x_{k}-\bar{x}_{k}-\alpha_{k}\mathcal{S}(x_{k},\zeta_{k})\|^{2}$
		$\displaystyle=\|x_{k}-\bar{x}_{k}\|^{2}-2\alpha_{k}\langle\mathcal{S}(x_{k},\zeta_{k}),x_{k}-\bar{x}_{k}\rangle+\alpha_{k}^{2}\|\mathcal{S}(x_{k},\zeta_{k})\|^{2}$
		$\displaystyle=\|x_{k}-\bar{x}_{k}\|^{2}-2\alpha_{k}\langle\nabla f(x_{k},\zeta_{k})+\nabla g({\color[rgb]{0,0,0}u_{k}},\zeta_{k}),x_{k}-\bar{x}_{k}\rangle+\alpha_{k}^{2}\|\mathcal{S}(x_{k},\zeta_{k})\|^{2}.$		(11)

Now, we refine the second term. First, from convexity of $f$ we have:

$$\langle\nabla f(x_{k},\zeta_{k}),x_{k}-\bar{x}_{k}\rangle\geq f(x_{k},\zeta_{k})-f(\bar{x}_{k},\zeta_{k}).$$

Then, from convexity of $g(\cdot,\zeta)$ and the definition of the gradient mapping $\mathcal{S}(\cdot,\zeta)$, we have:

	$\displaystyle\langle\nabla g({\color[rgb]{0,0,0}u_{k}},\zeta_{k}),x_{k}-\bar{x}_{k}\rangle=\langle\nabla g({\color[rgb]{0,0,0}u_{k}},\zeta_{k}),x_{k}-{\color[rgb]{0,0,0}u_{k}}\rangle+\langle\nabla g({\color[rgb]{0,0,0}u_{k}},\zeta_{k}),{\color[rgb]{0,0,0}u_{k}}-\bar{x}_{k}\rangle$
	$\displaystyle\geq\alpha_{k}\|\mathcal{S}(x_{k},\zeta_{k})\|^{2}-\alpha_{k}\langle\nabla f(x_{k},\zeta_{k}),\mathcal{S}(x_{k},\zeta_{k})\rangle+g({\color[rgb]{0,0,0}u_{k}},\zeta_{k})-g(\bar{x}_{k},\zeta_{k})$
	$\displaystyle\geq\alpha_{k}\|\mathcal{S}(x_{k},\zeta_{k})\|^{2}-\alpha_{k}\langle\nabla f(x_{k},\zeta_{k})+\nabla g(x_{k},\zeta_{k}),\mathcal{S}(x_{k},\zeta_{k})\rangle+g(x_{k},\zeta_{k})-g(\bar{x}_{k},\zeta_{k}).$

Replacing the previous two inequalities in (3), we obtain:

	$\displaystyle\|v_{k}-\bar{v}_{k}\|^{2}$	$\displaystyle\leq\|x_{k}-\bar{x}_{k}\|^{2}-2\alpha_{k}\left(f(x_{k},\zeta_{k})+g(x_{k},\zeta_{k})-f(\bar{x}_{k},\zeta_{k})-g(\bar{x}_{k},\zeta_{k})\right)$
		$\displaystyle\quad+2\alpha_{k}^{2}\langle\nabla f(x_{k},\zeta_{k})+\nabla g(x_{k},\zeta_{k}),\mathcal{S}(x_{k},\zeta_{k})\rangle-\alpha_{k}^{2}\|\mathcal{S}(x_{k},\zeta_{k})\|^{2}.$

Using that $2\langle u,v\rangle-\|v\|^{2}\leq\|u\|^{2}$ for all $v\in\mathbb{R}^{n}$, that $F(x_{k},\zeta_{k})=f(x_{k},\zeta_{k})+g(x_{k},\zeta_{k})$ and $\nabla F(x_{k},\zeta_{k})=\nabla f(x_{k},\zeta_{k})+\nabla g(x_{k},\zeta_{k})\in\partial F(x_{k},\zeta_{k})$, we further get:

$\displaystyle\|v_{k}-\bar{v}_{k}\|^{2}$

$\displaystyle\leq\|x_{k}-\bar{x}_{k}\|^{2}-2\alpha_{k}\left(F(x_{k},\zeta_{k})-F(\bar{x}_{k},\zeta_{k})\right)+\alpha_{k}^{2}\|\nabla F(x_{k},\zeta_{k})\|^{2}.$

Since $v_{k}$ depends on $\mathcal{F}_{[k-1]}\cup\{\zeta_{k}\}$, not on $\xi_{k}$, and the stepsize $\alpha_{k}$ does not depend on $(\zeta_{k},\xi_{k})$, then from basic properties of conditional expectation, we have:

	$\displaystyle\mathbb{E}_{\zeta_{k}}[\|v_{k}-\bar{v}_{k}\|^{2}|\mathcal{F}_{[k-1]}]$
	$\displaystyle\leq\|x_{k}-\bar{x}_{k}\|^{2}-2\alpha_{k}\mathbb{E}_{\zeta_{k}}[F(x_{k},\zeta_{k})-F(\bar{x}_{k},\zeta_{k})|\mathcal{F}_{[k-1]}]+\alpha_{k}^{2}\mathbb{E}_{\zeta_{k}}[\|\nabla F(x_{k},\zeta_{k})\|^{2}|\mathcal{F}_{[k-1]}]$
	$\displaystyle=\|x_{k}-\bar{x}_{k}\|^{2}-2\alpha_{k}\left(F(x_{k})-F(\bar{x}_{k})\right)+\alpha_{k}^{2}\mathbb{E}_{\zeta_{k}}[\|\nabla F(x_{k},\zeta_{k})\|^{2}|\mathcal{F}_{[k-1]}]$
	$\displaystyle\leq\|x_{k}-\bar{x}_{k}\|^{2}-2\alpha_{k}\left(F(x_{k})-F(\bar{x}_{k})\right)+\alpha_{k}^{2}\left(B^{2}+L\left(F(x_{k})-F(\bar{x}_{k})\right)\right)$
	$\displaystyle=\|x_{k}-\bar{x}_{k}\|^{2}-\alpha_{k}(2-\alpha_{k}L)\left(F(x_{k})-F(\bar{x}_{k})\right)+\alpha_{k}^{2}B^{2},$

where in the last inequality we used $x_{k}\in\mathcal{Y}$ and the stochastic bounded gradient inequality from Assumption 1. Now, taking expectation w.r.t. $\mathcal{F}_{[k-1]}$ we get:

$\displaystyle\mathbb{E}[\|v_{k}-\bar{v}_{k}\|^{2}]\leq\mathbb{E}[\|x_{k}-\bar{x}_{k}\|^{2}]-\alpha_{k}(2-\alpha_{k}L)\mathbb{E}[(F(x_{k})-F(\bar{x}_{k}))]+\alpha_{k}^{2}B^{2}.$

which concludes our statement.  

Next lemma gives a relation between $x_{k}$ and $v_{k-1}$ (see also (Nedich and Necoara,, 2019)).

Proof  Consider any $y\in\mathcal{Y}$ such that $(h(y,\xi_{k-1}))_{+}=0$. Then, using the nonexpansive property of the projection and the definition of $z_{k-1}$, we have:

	$\displaystyle\|x_{k}-y\|^{2}$	$\displaystyle=\|\Pi_{\mathcal{Y}}(z_{k-1})-y\|^{2}\leq\|z_{k-1}-y\|^{2}$
		$\displaystyle=\|v_{k-1}-y-\beta\frac{(h(v_{k-1},\xi_{k-1}))_{+}}{\|\nabla h(v_{k-1},\xi_{k-1})\|^{2}}\nabla h(v_{k-1},\xi_{k-1})\|^{2}$
		$\displaystyle=\|v_{k-1}-y\|^{2}+\beta^{2}\frac{(h(v_{k-1},\xi_{k-1}))_{+}^{2}}{\|\nabla h(v_{k-1},\xi_{k-1})\|^{2}}$
		$\displaystyle\quad-2\beta\frac{(h(v_{k-1},\xi_{k-1}))_{+}}{\|\nabla h(v_{k-1},\xi_{k-1})\|^{2}}\langle v_{k-1}-y,\nabla h(v_{k-1},\xi_{k-1})\rangle$
		$\displaystyle\leq\|v_{k-1}-y\|^{2}+\beta^{2}\frac{(h(v_{k-1},\xi_{k-1}))_{+}^{2}}{\|\nabla h(v_{k-1},\xi_{k-1})\|^{2}}-2\beta\frac{(h(v_{k-1},\xi_{k-1}))_{+}}{\|\nabla h(v_{k-1},\xi_{k-1})\|^{2}}(h(v_{k-1},\xi_{k-1}))_{+},$

where the last inequality follows from convexity of $(h(\cdot,\xi))_{+}$ and our assumption that $(h(y,\xi_{k-1}))_{+}=0$. After rearranging the terms and using that $v_{k-1}\in\mathcal{Y}$, we get:

	$\displaystyle\|x_{k}-y\|^{2}$	$\displaystyle\leq\|v_{k-1}-y\|^{2}-\beta(2-\beta)\left[\frac{(h(v_{k-1},\xi_{k-1}))_{+}^{2}}{\|\nabla h(v_{k-1},\xi_{k-1})\|^{2}}\right]$
		$\displaystyle\overset{\text{Assumption}\,\ref{assumption3}}{\leq}\|v_{k-1}-y\|^{2}-\beta(2-\beta)\left[\frac{(h(v_{k-1},\xi_{k-1}))_{+}^{2}}{B_{h}^{2}}\right],$

which concludes our statement.  
In the next sections, based on the previous two lemmas, we derive convergence rates for SSP algorithm depending on the (convexity) properties of $f(\cdot,\zeta)$.