STOCHASTIC PROJECTIVE SPLITTING: SOLVING SADDLE-POINT PROBLEMS WITH MULTIPLE REGULARIZERS
Year: 2021

Abstract:
We present a new, stochastic variant of the projective splitting (PS) family of algorithms for monotone inclusion problems. It can solve min-max and noncooperative game formulations arising in applications such as robust ML without the convergence issues associated with gradient descent-ascent, the current de facto standard approach in ML applications. Our proposal is the first version of PS able to use stochastic gradient oracles. It can solve min-max games while handling multiple constraints and nonsmooth regularizers via projection and proximal operators. Unlike other stochastic splitting methods that can solve such problems, our method does not rely on a product-space reformulation of the original problem. We prove almost-sure convergence of the iterates to the solution and a convergence rate for the expected residual. By working with monotone inclusions rather than variational inequalities, our analysis avoids the drawbacks of measuring convergence through the restricted gap function. We close with numerical experiments on a distributionally robust sparse logistic regression problem.

text: We present a new, stochastic variant of the projective splitting (PS) family of algorithms for monotone inclusion problems. It can solve min-max and noncooperative game formulations arising in applications such as robust ML without the convergence issues associated with gradient descent-ascent, the current de facto standard approach in ML applications. Our proposal is the first version of PS able to use stochastic gradient oracles. It can solve min-max games while handling multiple constraints and nonsmooth regularizers via projection and proximal operators. Unlike other stochastic splitting methods that can solve such problems, our method does not rely on a product-space reformulation of the original problem. We prove almost-sure convergence of the iterates to the solution and a convergence rate for the expected residual. By working with monotone inclusions rather than variational inequalities, our analysis avoids the drawbacks of measuring convergence through the restricted gap function. We close with numerical experiments on a distributionally robust sparse logistic regression problem.

heading: 1 INTRODUCTION
text: The most prominent application of optimization in ML is empirical risk minimization. However, inspired by the success of GANs (Goodfellow et al., 2014). , ML practitioners have developed more complicated min-max and adversarial optimization formulations (Yu et al., 2021; Kuhn et al., 2019; Shafieezadeh-Abadeh et al., 2015; Sinha et al., 2018; Lin et al., 2020; Namkoong & Duchi, 2016; Huang et al., 2017; Wadsworth et al., 2018; Zhang et al., 2018; Edwards & Storkey, 2015; Celis & Keswani, 2019). Solving these multi-player games leads to issues not seen when minimizing a single-player loss function. The competitive nature of a game leads to rotational dynamics that can cause intuitive gradient-based methods to fail to converge (Gidel et al., 2019; Daskalakis et al., 2018; Hsieh et al., 2020).
A mathematical framework underlying both convex optimization and saddle-point problems is the monotone inclusion problem; see Ryu & Boyd (2016) for an introduction. Methods developed for monotone inclusions will converge for convex-concave, games as they are explicitly designed to handle such problems’ governing dynamics. In recent years, monotone inclusion methods and theory have started to receive attention in the ML community (Diakonikolas, 2020; Liu et al., 2021; Ryu et al., 2020; Pathak & Wainwright, 2020), with a focus on monotone variational inequalities, which form a special case of monotone inclusions (Antonakopoulos et al., 2019; Gidel et al., 2019; Daskalakis et al., 2018; Hsieh et al., 2020; Mertikopoulos et al., 2019).
The most prevalent methods for solving min-max games in ML are variants of gradient descent-ascent (GDA). This method alternates between a gradient-descent step for the minimizing player and a gradient-ascent step for the maximizing player. Unfortunately, GDA requires additional assumptions to converge on convex-concave games, and it even fails for some simple 2D bilinear games (Gidel et al., 2019, Prop. 1). While there have been several approaches to modify either GDA (Chavdarova et al., 2021; Grnarova et al., 2021; Balduzzi et al., 2018) or the underlying game objective (Mescheder et al., 2018; Nagarajan & Kolter, 2017; Mescheder et al., 2017) to ensure convergence, this paper instead develops a method for solving monotone inclusions that can naturally handle game dynamics.
Our approach builds upon the recently proposed projective splitting (PS) method with forward steps (Johnstone & Eckstein, 2020b). PS is designed specifically for solving monotone inclusions, thus does not fall prey to the convergence issues that plague GDA, at least for convex-concave games. PS is within the general class of projective splitting methods invented by Eckstein & Svaiter (2008) and developed further in Eckstein & Svaiter (2009); Alotaibi et al. (2014); Combettes & Eckstein (2018); Eckstein (2017); Johnstone & Eckstein (2019; 2021; 2020a). These methods work by creating a separating hyperplane between the current iterate and the solution and then moving closer to the solution by projecting the current iterate onto this hyperplane (see Section 3 for an overview). Other than being able to natively handle game dynamics, the primary advantage of PS is that it fully splits problems involving an arbitrary number of regularizers and constraints. “Full splitting” means that the method can handle multiple regularizers and constraints through their respective individual proximal and projection operators, along with the smooth terms via gradients. What makes this useful is that many of the regularizers used in ML have proximal operators that are relatively easy to compute; see for example Parikh & Boyd (2013).
Despite these advantages, the preexisting PS framework has a significant drawback: it requires deterministic gradient oracles. This feature makes it impractical for application to large datasets for which stochastic oracles may be the only feasible option.
Contributions The primary contribution of this work is a new projective splitting algorithm that allows for a stochastic gradient oracle. We call the method stochastic projective splitting (SPS). Our method “fully splits” the monotone inclusion problem
Find z ∈ Rd s.t. 0 ∈ ∑n i=1Ai(z) +B(z), (1)
where B is monotone and L-Lipschitz and each Ai is maximal monotone and typically set valued, usually arising from a constraint or a nonsmooth regularizer in the underlying optimization problem or game; see for example Ryu & Boyd (2016) for definitions. For some example ML applications of (1), see Section 2 and Appendix A. Here, an algorithm that “fully splits” (1) means one whose computational steps each involve only the individual operators A1, . . . , An, B. Ours is the first method that can accomplish full splitting without a product-space reformulation that recasts (1) as a two-operator problem on a higher-dimensional space, a tactic whose disadvantages are discussed in Appendix F.7. Our method interrogates the Lipschitz operator B through a stochastic oracle. Previous methods splitting (1) have either required a deterministic oracle for B, or have made far more restrictive assumptions on the noise or the operators (Briceño-Arias & Combettes, 2011; Combettes & Pesquet, 2012; Malitsky & Tam, 2020; Bot et al., 2019; Van Dung & Vu, 2021) than we will require below. However, the stochastic methods of Alacaoglu et al. (2021) and Böhm et al. (2020), when combined with a product-space reformulation, can solve (1) when all the Ai are subdifferentials of convex functions; see Section 6.
When moving away from a deterministic gradient oracle in projective splitting, a key difficulty is that the generated hyperplanes do not guarantee separation between the solution and the current point. We solve this issue by relaxing the projection: we only update each iterate in the direction of the noisy projection and scale its movement by a decreasing stepsize that allows for control of the stochastic error. Using the framework of stochastic quasi-Fejér monotonicity (Combettes & Pesquet, 2015), we prove almost-sure convergence of the final iterate and do not require averaging of the iterates (Theorem 1, Section 5). We also provide a non-asymptotic convergence rate for the approximation residual (Theorem 2, Section 5).
A special case of SPS is the recently-developed Double Stepsize Extragradient Method (DSEG) (Hsieh et al., 2020). When n = 0 and therefore only B is present in (1), DSEG and SPS coincide. Thus, our method extends DSEG to allow for regularizers and constraints. Our analysis also provides a new interpretation for DSEG as a special case of projective splitting. Our nonasymptotic convergence rate for SPS also applies to DSEG under no additional assumptions. By contrast, the original convergence rate analysis for DSEG requires either strong monotonicity or an error bound.
We close with numerical experiments on a distributionally robust sparse logistic regression problem. This is a nonsmooth convex-concave min-max problem which can be converted to (1) with n = 2 set-valued operators. On this problems class, SPS compares well to the possible alternative splitting methods.
Non-monotone problems The work of Hsieh et al. (2020) included a local convergence analysis for DSEG applied to locally monotone problems. For min-max problems, if the objective is locally convex-concave at a solution and DSEG is initialized in close proximity, then for small enough stepsizes it converges to the solution with high probability. It is possible to extend this result to SPS, along with our convergence rate analysis. This result is beyond the scope of this work, but Appendix J provides a proof sketch.

heading: 2 BACKGROUND ON MONOTONE INCLUSIONS
text: Since they are so important to SPS, this section provides some background material regarding monotone inclusions, along with their connections to convex optimization, games, and ML. Appendix G discusses their connections to variational inequalities. For a more thorough treatment, we refer to Bauschke & Combettes (2017). See Appendix A for a longer discussion of the applications of monotone inclusions to ML along with several examples.
Fundamentals Let f : Rd → R ∪ {∞} be closed, convex, and proper (CCP). Recall that its subdifferential ∂f is given by ∂f(x) .= {g : f(y) ≥ f(x) + g>(y − x)}. The map ∂f has the property
u ∈ ∂f(x), v ∈ ∂f(y) =⇒ (u− v)>(x− y) ≥ 0, and any point-to-set map having this property is called a monotone operator. A monotone operator T is called maximal if no additional points can be included in the image T (x) of any x ∈ Rd without violating the above property (Bauschke & Combettes, 2017, Def. 20.20). Subgradient maps of CCP functions are maximal (Bauschke & Combettes, 2017, Thm. 20.25). A minimizer of f is any x∗ such that 0 ∈ ∂f(x∗). This is perhaps the simplest example of a monotone inclusion, the problem of finding x such that 0 ∈ T (x), where T is a monotone operator. If f is smooth, then ∂f(x) = {∇f(x)} for all x, and the monotone inclusion 0 ∈ ∂f(x) is equivalent to the first-order optimality condition 0 = ∇f(x). Under certain regularity conditions (Bauschke & Combettes, 2017, Cor. 16.5), minimizing a sum of CCP functions f1, . . . , fn is equivalent to solving the monotone inclusion formed from the sum of their subdifferentials:
x∗ ∈ arg min x∈Rd n∑ i=1 fi(x) ⇐⇒ 0 ∈ n∑ i=1 ∂fi(x ∗). (2)
As throughout this paper for all set addition operations, the summation on the right-hand side of (2) is the Minkowski sum ∑n i=1 Si = { ∑n i=1 si | si ∈ Si ∀ i ∈ 1..n}. For a convex set X , a constraint x ∈ C for some convex set C may be imposed by setting one of the fi to be the indicator function ιC , defined by ιC(x) = 0 for x ∈ C and ιC(x) = +∞ for x 6∈ C. Indicator functions of closed convex sets are CCP (Bauschke & Combettes, 2017, Ex. 1.25), and the subgradient map of ιC is also referred to as the normal cone map NC of C (Bauschke & Combettes, 2017, Def. 6.37). Multiple constraints may be imposed by including multiple indicator functions in (2).
ML applications The form (2) can be used to model ML problems with multiple constraints and/or nonsmooth regularizers, including sparse and overlapping group lasso (Jacob et al., 2009), sparse and low-rank matrix estimation problems (Richard et al., 2012), and rare feature selection (Yan & Bien, 2020); see Pedregosa & Gidel (2018) for an overview.
Games Consider a two-player noncooperative game in which each player tries to selfishly minimize its own loss, with each loss depending on the actions of both players. Typically, the goal is to find a Nash equilibrium, in which neither player can improve its loss by changing strategy:
x∗ ∈ arg min x∈Θ F (x, y∗) and y∗ ∈ arg min y∈Ω G(x∗, y). (3)
Assuming that the admissible strategy sets Θ ⊆ Rdx and Ω ⊆ Rdy are closed and convex and that F and G are differentiable, then writing the first-order necessary conditions for each optimization problem in (3) yields
0 ∈ [ ∇xF (x∗, y∗) ∇yG(x∗, y∗) ] + ( NΘ(x ∗)×NΩ(y∗) ) . (4)
IfG = −F , then (3) is a min-max game. If F is also convex in x and concave in y, thenB : (x, y) 7→ (∇xF (x, y),−∇yF (x, y))> is monotone1 on Rdx+dy (Rockafellar, 1970). In many applications, B is also Lipschitz continuous. In this situation, (4) is a monotone inclusion involving two operators B and NΘ×Ω, with B being Lipschitz. Using the simultaneous version of GDA on (3) is equivalent to applying the forward-backward method (FB) (Bauschke & Combettes, 2017, Thm. 26.14) to (4). However, convergence of FB requires that the operator B be cocoercive (Bauschke & Combettes, 2017, Def. 4.10), and not merely Lipschitz (Bauschke & Combettes, 2017, Thm. 26.14). Thus, simultaneous GDA fails to converge for (3) without additional assumptions; see Gidel et al. (2019, Prop. 1) for a simple counterexample.
Regularizers and further constraints may be imposed by adding more operators to (4). For example, if one wished to apply a (nonsmooth) convex regularizer r : Rdx → R ∪ {+∞} to the x variables and a similar regularizer d : Rdy → R ∪ {+∞} to the y variables, one would add the operator A2 : (x, y) 7→ ∂r(x)× ∂d(y) to the right-hand side of (4).
ML applications of games Distributionally robust supervised learning (DRSL) is an emerging framework for improving the stability and reliability of ML models in the face of distributional shifts (Yu et al., 2021; Kuhn et al., 2019; Shafieezadeh-Abadeh et al., 2015; Sinha et al., 2018; Lin et al., 2020; Namkoong & Duchi, 2016). Common approaches to DRSL formulate the problem as a min-max game between a learner selecting the model parameters and an adversary selecting a worst-case distribution subject to some ambiguity set around the observed empirical distribution. This min-max problem is often further reduced to either a finite-dimensional saddlepoint problem or a convex optimization problem.
DRSL is a source of games with multiple constraints/regularizers. One such formulation, based on Yu et al. (2021), is discussed in the experiments below. The work in Namkoong & Duchi (2016) uses an ambiguity set based on f -divergences, while Sinha et al. (2018) introduce a Lagrangian relaxation of the Wasserstein ball. When applied to models utilizing multiple regularizers (Jacob et al., 2009; Richard et al., 2012; Yan & Bien, 2020), both of these approaches lead to min-max problems with multiple regularizers.
Other applications of games in ML, although typically nonconvex, include generative adversarial networks (GANs) (Goodfellow et al., 2014; Arjovsky et al., 2017; Loizou et al., 2020; 2021; Mishchenko et al., 2020), fair classification (Wadsworth et al., 2018; Zhang et al., 2018; Edwards & Storkey, 2015; Celis & Keswani, 2019), and adversarial privacy (Huang et al., 2017).
Resolvents, proximal operators, and projections A fundamental computational primitive for solving monotone inclusions is the resolvent. The resolvent of a monotone operator A is defined to be JA . = (I + A)−1, where I is the identity operator and the inverse of any operator T is simply T−1 : x 7→ {y : Ty 3 x}. If A is maximal monotone, then for any ρ > 0, JρA is single valued, nonexpansive, and has domain equal to Rd (Bauschke & Combettes, 2017, Thm. 21.1 and Prop. 23.8). Resolvents generalize proximal operators of convex functions: the proximal operator of a CCP function f is
proxρf (t) . = arg min
x∈Rd
{ ρf(x) + (1/2)‖x− t‖2 } .
It is easily proved that proxρf = Jρ∂f . Like proximal operators, resolvents generalize projection onto convex sets: if f = ιC , then JρNC = proxρf = projC for any ρ > 0. In many ML applications, proximal operators, and hence resolvents, are relatively straightforward to compute. For examples, see Parikh & Boyd (2013, Sec. 6).
Operator splitting methods Operator splitting methods attempt to solve monotone inclusions such as (1) by a sequence of operations that each involve only one of the operators A1, . . . , An, B. Such methods are often presented in the context of convex optimization problems like (2), but typically apply more generally to monotone inclusions such as (1). In the specific context of (1), each iteration of such a method ideally handles each Ai via its resolvent and the Lipschitz operator B by explicit (not stochastic) evaluation. This is a feasible approach if the original problem can be decomposed in
1Sufficient conditions for the monotonicity of (4) in the case where G 6= −F are discussed in e.g. Scutari et al. (2014); Briceño-Arias & Combettes (2013).
such a way that the resolvents of each Ai are relatively inexpensive to compute, and full evaluations of B are possible. Although not discussed here, more general formulations in which matrices couple the arguments of the operators can broaden the applicability of operator splitting methods.

heading: 3 THE PROJECTIVE SPLITTING FRAMEWORK
text: Before introducing our proposed method, we give a brief introduction to the projective splitting class of methods.
The extended solution set Projective splitting is a primal-dual framework and operates in an extended space of primal and dual variables. Rather than directly finding a solution to (1), we find a point in the extended solution set (or Kuhn-Tucker set)
S .= { (z, w1, . . . , wn+1) ∣∣∣ wi ∈ Ai(z)∀ i ∈ 1..n, wn+1 = B(z),∑n+1i=1 wi = 0} . (5)
Given p∗ = (z∗, w∗1 . . . , w ∗ n+1) ∈ S , it is straightforward to see that z∗ solves (1). Conversely, given a solution z∗ to (1), there must exist w∗1 , . . . , w ∗ n+1 such that (z ∗, w∗1 , . . . , w ∗ n+1) ∈ S. Suppose p∗ = (z∗, w∗1 . . . , w ∗ n+1) ∈ S. Since z∗ solves (1), z∗ is typically referred to as a primal solution. The vectors w∗1 , . . . , w ∗ n+1 solve a dual inclusion not described here, and are therefore called a dual solution. It can be shown that S is closed and convex; see for example Johnstone & Eckstein (2020b). We will assume throughout that a solution to (1) exists, therefore the set S is nonempty.
Separator-projection framework Projective splitting methods are instances of the general separator-projection algorithmic framework for locating a member of a closed convex set S within a linear space P . Each iteration k of algorithms drawn from this framework operates by finding a set Hk that separates the current iterate pk ∈ P from S, meaning that S is entirely in the set and pk typically is not. One then attempts to “move closer" to S by projecting the pk onto Hk. In the particular case of projective splitting applied to the problem (1) using (5), we select the space P to be
P .= { (z, w1, . . . , wn+1) ∈ R(n+2)d ∣∣∣ ∑n+1i=1 wi = 0} , (6)
and each separating set Hk to be the half space {p ∈ P | ϕk(p) ≤ 0} generated by an affine function ϕk : P → R. The general intention is to construct ϕk such that ϕk(pk) > 0, but ϕk(p∗) ≤ 0 for all p∗ ∈ S. The construction employed for ϕk in the case of (1) and (5) is of the form
ϕk(z, w1, . . . , wn+1) . = ∑n+1 i=1 〈z − xki , yki − wi〉 (7)
for some points (xki , y k i ) ∈ R2d, i ∈ 1..(n + 1), that must be carefully chosen (see below). Any function of the form (7) can be shown to be affine when restricted to P . As mentioned above, the standard separator-projection algorithm obtains its next iterate pk+1 by projecting pk onto Hk. This calculation involves the usual projection step for a half space, namely
pk+1 = pk − αk∇ϕk, where αk = ϕk(pk)/‖∇ϕk‖2, (8)
and the gradient ∇ϕk is computed relative to P , thus resulting in pk+1 ∈ P , i.e. ∇ϕk =(∑n+1 i=1 y k i , x k 1 − x̄k, . . . , xn+1 − x̄k ) where x̄k = 1n+1 ∑n+1 i=1 x k i .

heading: 4 PROPOSED METHOD
text: The proposed method is given in Algorithm 1 and called Stochastic Projective Splitting (SPS). Unlike prior versions of projective splitting, SPS does not employ the stepsize αk of (8) that places the next iterate exactly on the hyperplane given by ϕk(p) = 0. Instead, it simply moves in the direction−∇ϕk with a pre-defined stepsize {αk}. This fundamental change is required to deal with the stochastic noise on lines 6 and 8. This noise could lead to the usual choice of αk defined in (8) being unstable and difficult to analyze. In order to guarantee convergence, the parameters αk and ρk must be chosen to satisfy certain conditions given below. Note that the gradient is calculated with respect to the subspace P defined in (6); since the algorithm is initialized within P , it remains in P , within which ϕk is affine. Collectively, the updates on lines 9-10 are equivalent to pk+1 = pk − αk∇ϕk, where pk = (zk, wk1 , . . . , w k n+1).
Note that SPS does not explicitly evaluate ϕk, which is only used in the analysis, but it does keep track of (xki , y k i ) for i ∈ 1..(n+ 1). The algorithm’s memory requirements scale linearly with the number of nonsmooth operators n in the inclusion (1), with the simplest implementation storing (3n+5)d working-vector elements. This requirement can be reduced to (n+7)d through a technique discussed in Appendix H. In most applications, n will be small, for example 2 or 3.
Updating (xki , yki ) The variables (xki , yki ) are updated on lines 3-8 of Algorithm 1, in which ek and k are Rd-valued random variables defined on a probability space (Ω,F , P ). For B we use a new, noisy version of the two-forward-step procedure from Johnstone & Eckstein (2020b). For each Ai, i ∈ 1..n, we use the same resolvent step used in previous projective splitting papers, originating with (Eckstein & Svaiter, 2008). In the case k = ek = 0, the selection of the (xki , y k i ) is identical to that proposed by Johnstone & Eckstein (2020b), resulting in the hyperplane {p : ϕk(p) = 0} strictly separating pk from S. SPS achieves full splitting of (1): each Ai is processed separately using a resolvent and the Lipschitz term B is processed via a stochastic gradient oracle. When the Ai arise from regularizers or constraints, as discussed in Section 2, their resolvents can be readily computed so long as their respective proximal/projection operators have a convenient form.
Noise assumptions Let Fk . = σ(p1, . . . , pk) and Ek . = σ( k). The stochastic estimators for the gradients, rk and ykn+1, are assumed to be unbiased, that is, the noise terms have mean 0 conditioned on the past:
E[ k|Fk] = 0, E[ek|Fk] = 0 a.s. (9) We impose the following mild assumptions on the variance of the noise:
E [ ‖ k‖2|Fk ] ≤ N1 +N2‖B(zk)‖2 a.s. (10)
E [ ‖ek‖2|Fk, Ek ] ≤ N3 +N4‖B(xkn+1)‖2 a.s., (11)
where 0 ≤ N1, N2, N3, N4 <∞. We do not require ek and k to be independent of one another.
Stepsize choices The stepsizes ρk and αk are assumed to be deterministic. A constant stepsize choice which attains a non-asymptotic convergence rate will be considered in the next section (Theorem 2). The stepsize conditions we will impose to guarantee almost-sure convergence (Theorem 1) are ∑∞
k=1 αkρk =∞, ∑∞ k=1 α 2 k <∞, ∑∞ k=1 αkρ 2 k <∞, and ρk ≤ ρ < 1/L. (12)
For example, in the case L = 1, a particular choice which satisfies these constraints is
αk = k −0.5−p for 0 < p < 0.5, and ρk = k−0.5+t for p ≤ t < 0.5p+ 0.25.
For simplicity, the stepsizes τ used for the resolvent updates in lines 3-5 are fixed, but they could be allowed to vary with both i and k so long as they have finite positive lower and upper bounds.
Algorithm 1: Stochastic Projective Splitting (SPS) Input :p1 = (z1, w11, . . . , w1n+1) s.t. ∑n+1 i=1 w 1 i = 0, {αk, ρk}∞k=1, τ > 0
1 for k = 1, 2, . . . do 2 for i ∈ 1..n do 3 tki = z
k + τwki 4 xki = JτAi(t k i ) 5 yki = τ −1(tki − xki )
6 rk = B(zk) + k // k is unknown noise term 7 xkn+1 = z
k − ρk(rk − wkn+1) 8 ykn+1 = B(x k n+1) + e k // ek is unknown noise term
9 zk+1 = zk − αk ∑n+1 i=1 y k i
10 wk+1i = w k i − αk(xki − 1n+1 ∑n+1 i=1 x k i ) i ∈ 1..(n+ 1)

heading: 5 MAIN THEORETICAL RESULTS
text: Theorem 1. Suppose A1, . . . , An are maximal monotone, B is L-Lipschitz and monotone, and a solution to (1) exists. For Algorithm 1, suppose (9)-(12) hold. Then with probability one it holds that zk → z∗, where z∗ solves (1). Further, with probability one, xki → z∗ for i = 1, . . . , n.
Proof sketch Theorem 1 is proved in Appendix C, but we provide a brief sketch here. The proof begins by deriving a simple recursion inspired by the analysis of SGD (Robbins & Monro, 1951). Since pk+1 = pk − αk∇ϕk, a step of projective splitting can be viewed as GD applied to the affine hyperplane generator function ϕk. Thus, for any p∗ ∈ P ,
‖pk+1 − p∗‖2 = ‖pk − p∗‖2 − 2αk〈∇ϕk, pk − p∗〉+ α2k‖∇ϕk‖2
= ‖pk − p∗‖2 − 2αk(ϕk(pk)− ϕk(p∗)) + α2k‖∇ϕk‖2, (13) where in the second equation we have used that ϕk(p) is affine on P . The basic strategy is to show that, for any p∗ ∈ S,
E[‖∇ϕk‖2|Fk] ≤ C1‖pk − p∗‖2 + C2 a.s. for some C1, C2 > 0. This condition allows one to establish stochastic quasi-Fejér monotonicity (SQFM) (Combettes & Pesquet, 2015, Proposition 2.3) of the iterates to S. One consequence of SQFM is that with probability one there exists a subsequence vk such that ϕvk(p
vk) − ϕvk(p∗) converges to 0. Furthermore, roughly speaking, we show that ϕk(pk)− ϕk(p∗) provides an upper bound on the following “approximation residual" for SPS:
Gk . = ∑n i=1 ‖yki − wki ‖2 + ∑n i=1 ‖zk − xki ‖2 + ‖B(zk)− wkn+1‖2. (14)
Gk provides an approximation error for SPS, as formalized in the following lemma: Lemma 1. For SPS, pk = (zk, wk1 , . . . , wkn+1) ∈ S if and only if Gk = 0. Since yki ∈ Ai(xki ) for i ∈ 1..n, having Gk = 0 implies that zk = xki , wki = yki , and thus wki ∈ Ai(zk) for i ∈ 1..n. Since wkn+1 = B(zk) and ∑n+1 i=1 w k i = 0, it follows that z
k solves (1). The reverse direction is proved in Appendix D.
The quantity Gk generalizes the role played by the norm of the gradient in algorithms for smooth optimization. In particular, in the special case where n = 0 and B(z) = ∇f(z) for some smooth convex function f , one has Gk = ‖∇f(zk)‖2. Combining the properties of Gk with other results following from SQFM (such as boundedness) will allow us to derive almost-sure convergence of the iterates to a solution of (1).
Convergence rate We can also establish non-asymptotic convergence rates for the approximation residual Gk: Theorem 2. Fix the total iterations K ≥ 1 of Algorithm 1 and set
∀k = 1, . . . ,K : ρk = ρ . = min { K−1/4, 1/2L } and αk = Cfρ2 (15)
for some Cf > 0. Suppose (9)-(11) hold. Then (1/K) ∑K j=1 E[Gj ] = O(K−1/4)
where the constants are given (along with the proof) in Appendix E.
Theorem 2 implies that if we pick an iterate J uniformly at random from 1..K, then the expected value of GJ is O(K−1/4). As far as we know, this is the first convergence rate for a stochastic fullsplitting method solving (1) in the general discontinuous (i.e. set-valued) monotone inclusion case, and it is not clear whether it can be improved, either by a better analysis or a better method. Faster rates are certainly possible for deterministic methods under various continuity assumptions; Tseng’s method obtains O(K−1) rate (Monteiro & Svaiter, 2010) and the accelerated Halpern iteration under Lipschitz continuity obtains O(K−2) rate (Diakonikolas, 2020). While our rate may seem slow, it is worth remembering that (1) features n discontinuous operators Ai, so we expect rates at least as slow as nonsmooth convex optimization, but perhaps worse because (1) is far more general than convex optimization. For a different error metric, the restricted gap function, in the special case of variational inequalities, faster rates have been established in Juditsky et al. (2011) and Böhm et al. (2020). However, it is unclear how to relate the restricted gap function to Gk, so these rates may not be directly comparable to Theorem 2.

heading: 6 RELATED WORK
text: Arguably the three most popular classes of operator splitting algorithms are forward-backward splitting (FB) (Combettes & Pesquet, 2011), Douglas-Rachford splitting (DR) (Lions & Mercier, 1979), and Tseng’s method (Tseng, 2000). The extragradient method (EG) is similar to Tseng’s method, but has more projection steps per iteration and only applies to variational inequalities (Korpelevich, 1977; Nemirovski, 2004; Li et al., 2021). The popular Alternating Direction Method of Multipliers (ADMM), in its standard form, is a dual application of DR (Gabay, 1983). The three-operator splitting method (Davis & Yin, 2017) can only be applied to (1) if B is cocoercive rather than merely Lipchitz, and thus its usefulness is mostly limited to optimization applications and not games. FB, DR, and Tseng’s method apply to monotone inclusions involving two operators, with varying assumptions on one of the operators. It is possible to derive splitting methods for the more complicated inclusion (1), involving more than two operators, by applying an appropriate 2-operator splitting method such as Tseng’s method to a product-space reformulation (PSR) (Briceño-Arias & Combettes, 2011; Combettes & Pesquet, 2012) (for more on PSR, see Appendix F). The recently developed forward-reflected-backward (FRB) method (Malitsky & Tam, 2020) can be used in the same way. However, there are several disadvantages to using a PSR, as discussed in Appendix F.7.
By using a PSR, the stochastic methods of Alacaoglu et al. (2021) and Böhm et al. (2020) can be applied to (1) in the case that eachAi is a subdifferential. Both of these methods are analyzed in terms of the restricted gap function. This merit function has a drawback compared with our approximation residual in that it requires one to find a bound for the iterates. However, Alacaoglu et al. (2021) and Böhm et al. (2020) do not provide such a bound, meaning that their convergence rate results are somewhat incomplete. We discuss this issue in Appendix G.
Theoretical convergence of the method of Böhm et al. (2020) requires the use of averaging, since the final iterate does not converge for certain problems (Hsieh et al., 2020). Empirically, averaging tends to be slow and to destroy regularizer-induced structural properties such as sparsity or low matrix rank, so its utility is largely theoretical and it is usually avoided in practice. Furthermore, averaging loses even its theoretical benefits for nonconvex problems, so its use in such cases is rarer still. Another drawback of the analysis of Böhm et al. (2020) is that, unlike in SPS, the resolvent (proximal) stepsizes also need to vanish.
The method of Alacaoglu et al. (2021) applies variance reduction techniques to FRB. It only applies to finite-sum problems and requires the periodic computation of a full batch gradient, making it somewhat less flexible and scalable than our method. On the other hand, it has an accelerated ergodic rate for the restricted gap function in the variational inequality setting. We compare the empirical performance of SPS with Alacaoglu et al. (2021), Böhm et al. (2020), and several deterministic methods using PSR in the numerical experiments described in Section 7.
Additional related work is discussed in Appendix B.

heading: 7 EXPERIMENTS
text: We now present some numerical results on distributionally robust supervised learning (DRSL) problems. We follow the approach of Yu et al. (2021), which introduced a min-max formulation of Wasserstein DRSL. While other approaches reduce the problem to convex optimization, Yu et al. (2021) reduce it to a finite-dimensional min-max problem amenable to the use of stochastic methods on large datasets. However, unlike our proposed SPS method, the variance-reduced extragradient method that Yu et al. (2021) propose cannot handle multiple nonsmooth regularizers or constraints on the model parameters. Consequently, we consider distributionally robust sparse logistic regression (DRSLR), a problem class equivalent to that considered in Yu et al. (2021), but with an added `1 regularizer, a standard tool to induce sparsity. See the Appendix I for the full problem definition.
We compared our SPS method to several methods for solving DRSLR for a collection of real datasets from the LIBSVM repository (Chang & Lin, 2011). We implemented SPS with αk = Cdk−0.51 and ρk = Cdk
−0.25 and called it SPS-decay. We also implement SPS with the fixed stepsize given in (15) and called it SPS-fixed. We compared the method to deterministic projective splitting (Johnstone & Eckstein, 2020b) and the following methods based on PSR: Tseng’s method (Tseng, 2000; Combettes & Pesquet, 2012), the forward-reflected-backward (FRB) method (Malitsky & Tam, 2020), the
stochastic Tseng (S-Tseng) method of Böhm et al. (2020), and the variance-reduced stochastic FRB method (Alacaoglu et al., 2021), abbreviated FRB-VR. The S-Tseng and FRB-VR algorithms appear to be the only stochastic splitting methods other than SPS applicable to the tested problem class.
Figure 1 show results for three LIBSVM standard datasets: epsilon2 (m = 4 · 105, d = 2000), SUSY (Baldi et al., 2014; Dua & Graff, 2017) (m = 2 · 106, d = 18), and real-sim3 (m = 72,309, d = 20,958).
To measure the progress of the algorithms, we used the “approximation residual” Rk defined in Appendix F. As with Gk, having Rk = 0 implies that zk solves (1). We use Rk instead of Gk because it is also possible to compute essentially the same measure of convergence from the iterates of the other tested algorithms, establishing a fair comparison. Appendix F provides the details of the derivation of the residual measure for each algorithm, explores the relationship between Rk and Gk, and provides additional implementation details.
Figure 1 plots the approximation residual versus running time for all seven algorithms under consideration. The computations were performed using Python 3.8.3 and numpy on a 2019 MacBook Pro with a 2.4GHz 8-core Intel I9 processor and 32GB of RAM . Being a stochastic method, SPS-decay seems to outperform the deterministic methods at obtaining a medium-accuracy solution quickly. It also seems to outperform the stochastic PSR-based methods S-Tseng and FRB-VR.

heading: 8 CONCLUSIONS AND FUTURE WORK
text: We have developed and analyzed a stochastic splitting method that can handle min-max problems with multiple regularizers and constraints. Going forward, this development should make it possible to incorporate regularizers and constraints into adversarial formulations trained from large datasets.
Recent versions of deterministic projective splitting (Combettes & Eckstein, 2018; Johnstone & Eckstein, 2020b) allow for asynchronous and incremental operation, meaning that not all operators need to be activated at every iteration, with some calculations proceeding with stale inputs. Such characteristics make projective splitting well-suited to distributed implementations. Many of our SPS results may be extended to allow for these variations, but we leave those extensions to future work.

heading: A ML APPLICATIONS OF THE MONOTONE INCLUSION (1)
text: There are two main classes of applications of (1) in ML: optimization problems and saddle-point games.
Optimization Problems In this case the monotone inclusion arises from finding the zero of a sum of subgradients of convex functions, as discussed in Section 2. It is typical in ML to solve the empirical risk minimization problem
min x∈Rd
1
m m∑ j=1 fj(x) + n∑ i=1 ri(x) (16)
over a size-m dataset. Usually, the gradient of the loss function fj for each datapoint j is Lipschitz continuous. The terms ri may be regularizers used to reduce overfitting or encourage structural properties such as sparsity or low matrix rank. They also may represent constraints on the parameters such as nonnegativity or the being in the probability simplex. Crucially, these regularizers are rarely differentiable. The first-order necessary condition for the solution of (16) is
0 ∈ ∇f(x∗) + n∑ i=1 ∂ri(x ∗), (17)
where f(x) .= 1m ∑m j=1 fj(x), thus ∇f(x) . = 1m ∑m j=1∇fj(x). The inclusion (17) is a special case of (1), and our method may use the standard stochastic oracle for∇f(x), namely 1
|B| ∑ j∈B ∇fj(z)
which subsamples a randomly selected minibatch of datapoints B ∈ {1, . . . ,m}.
Games Consider the following nonsmooth Nash equilibrium problem
x∗ ∈ arg min x∈Rdx F (x, y∗) + n1∑ i=1 ri(x) and y∗ ∈ arg min y∈Rdy G(x∗, y) + n2∑ i=1 di(y). (18)
The terms ∑n1 i=1 ri(x) and ∑n2 i=1 di(y) once again represent regularizers and constraints on each player’s strategy. Note that min-max (saddle-point) problems correspond to having F (x, y) = −G(x, y). Under appropriate convexity conditions and constraint qualifications, the solutions of (18) correspond to the solutions of the following monotone inclusion in the form of (1):
0 ∈ [ ∇xF (x∗, y∗) ∇yG(x∗, y∗) ] + max{n1,n2}∑ i=1 ( ∂ri(x ∗)× ∂di(y∗) )
(19)
where for i > min{n1, n2} we include “dummy functions", either ri(x) = 0 when n1 < n2 or di(y) = 0 when n1 < n2. If the functions F and G arise as averages in the same we as f in (16), then our method may again use a stochastic oracle for them.
Distributionally-Robust ML One example application of (19) is distributionally-robust ML, as demonstrated in the numerical experiment in Section 7. The full problem statement is given in Appendix I.
Lagrangian Duality Another application of (19) is constrained optimization via Lagrangian duality. Consider
min x∈Rd
{ f(x) +
n∑ i=1 ri(x)
} s.t. hj(x) ≤ 0 j = 1, . . . , p.
As in (16), f is a loss function and the ri may represent regularizers and (“simple”) constraints; in addition, there are p functional constraints on the model parameters x. Introducing Lagrange multipliers γ ∈ Rp, the problem can be written as
min x∈Rd max γ∈Rp+ f(x) + n∑ i=1 ri(x) + p∑ j=1 γjhj(x)  . Under appropriate convexity conditions and constraint-qualifications, this reduces to the following inclusion in the form of (1):
0 ∈ [ ∇f(x) + ∑p j=1 γj∇hj(x)
−h(x)
] + n∑ i=1 ( ∂ri(x ∗)× {0} )
where h(x) = [h1(x), h2(x), . . . , hp(x)]>. For certain choices of h, such as linear or quadratic functions, the first term above is monotone and (locally) Lipschitz continuous (Alacaoglu et al., 2021).
Bilinear Games with Many Constraints Finally, consider the bilinear saddlepoint problem subject to multiple constraints:
min x∈Rd max y∈Rd
x>Dy s.t. x ∈ C1j j = 1, . . . , n1,
y ∈ C2j j = 1, . . . , n2.
Under some regularity conditions, this problem reduces to the inclusion
0 ∈ [ Dy∗ −D>x∗ ] + max{n1,n2}∑ j=1 ( NC1j (x ∗)×NC2j (y ∗) ) ,
where we introduce additional “dummy” sets C1j = Rd or C2j = Rd when n1 6= n2. The first term is linear and skew symmetric, and therefore can easily be shown to be Lipschitz continuous and monotone. If all the constraint sets are closed and convex, then the rest of the terms are maximal monotone, then the problem is of the form (1), meaning that projective splitting may be applied, possibly using a stochastic oracle for the first term.

heading: B ADDITIONAL RELATED WORK
text: The preprint by Bot et al. (2019) develops a stochastic version of Tseng’s method under the requirement that the noise variance converges to 0. In ML, this could be achieved with the use of perpetually increasing batch sizes, a strategy that is impractical in many scenarios. The stochastic version of FRB proposed by Van Dung & Vu (2021) has more practical noise requirements, but has stronger assumptions on the problem which are rarely satisfied in ML applications: either uniform/strong monotonicity or a bounded domain. The papers by Yurtsever et al. (2016) and Pedregosa et al. (2019) consider stochastic variants of three-operator splitting, but require B in (1) to be cocoercive, essentially restricting them to optimization problems.
There are several alternatives to the (stochastic) extragradient method that reduce the number of gradient evaluations per iteration from two to one (Hsieh et al., 2019; Malitsky & Tam, 2020; Gidel et al., 2019). However, these methods have more stringent stepsize limits, making it unclear a priori whether they will outperform two-step methods.
DSEG is a stochastic version of EG (Hsieh et al., 2020). The primary innovation of DSEG is using different stepsizes for the extrapolation and update steps, thereby resolving some of the convergence issues affecting stochastic EG. As noted earlier, DSEG is the special case of our SPS method in which n = 0, that is, no regularizers/constraints are present in the underlying game. The analysis in (Hsieh et al., 2020) also did not consider the fixed stepsize choice given in Theorem 2.
In the context of GANs, several methods have been developed based on a variational inequality/monotone inclusion approach (Gidel et al., 2019; Daskalakis et al., 2018; Hsieh et al., 2019; 2020; Böhm et al., 2020). Many of these papers point out that variational inequalities provide a principled framework for studying the GAN training problem and correcting some of the flaws in the standard method GDA.

heading: C PROOF OF THEOREM 1

heading: C.1 STOCHASTIC QUASI-FEJER MONOTONICITY
text: The key to the analysis is showing that the algorithm satisfies Stochastic Quasi-Fejer Monotonicity (Combettes & Pesquet, 2015). Lemma 2 ((Combettes & Pesquet, 2015), Proposition 2.3). Suppose pk is a sequence of Rd-valued random variables defined on a probability space (Ω,F , P ). Let Fk = σ(p1, . . . , pk). Let F be a nonempty, closed subset of Rd. Suppose that, for every p ∈ F , there exists χk(p) ≥ 0, ηk(p) ≥ 0, νk(p) ≥ 0 such that ∑∞ k=1 χ k(p) <∞, ∑∞ k=1 η k(p) <∞ and
(∀k ∈ N) E[‖pk+1 − p‖2|Fk] ≤ (1 + χk(p))‖pk − p‖2 − νk(p) + ηk(p).

heading: Then the following hold:
text: 1. (∀p ∈ F ) : ∑∞ k=1 ν k(p) <∞ a.s.
2. pk is bounded a.s. 3. There exists Ω̃ such that P [Ω̃] = 1 and { ‖pk(ω) − p‖ } converges for every ω ∈ Ω̃ and
p ∈ F .
C.2 IMPORTANT RECURSION FOR SPS
The following lemma summarizes the key recursion satisfied by Algorithm 1, to which we will apply Lemma 2. Recall that L is the Lipschitz constant of B. Lemma 3. For Algorithm 1, suppose (9)–(11) hold and
ρk ≤ ρ < 1/L. (20)

heading: Let
text: Tk . = τ
ρ n∑ i=1 ‖yki − wki ‖2 + 1 ρτ n∑ i=1 ‖zk − xki ‖2 + 2(1− ρL)‖B(zk)− wkn+1‖2
then for all p∗ ∈ S, with probability one E[‖pk+1 − p∗‖2|Fk] ≤ (1 + C1α2k + C3αkρ2k)‖pk − p∗‖2 − αkρkTk + C2α2k + C4αkρ2k (21)
where C1, . . . , C4 are nonegative constants defined in (33), (34), (48), and (49) below, respectively.
Note that Tk is a scaled version of the approximation residual Gk defined in (14).
We proceed to first prove Lemma 3 and then exploit the implications of Lemma 2. Referring to (10) and (11), let N .= maxj∈1..4Nj . To simplify the constants, we will use N in place of Nj for the noise variance bounds given in (10)-(11).

heading: C.3 UPPER BOUNDING THE GRADIENT
text: Throughout the analysis, we fix some p∗ = (z∗, w∗1 . . . , w ∗ n+1) ∈ S. All statements are with probability one (almost surely), but for brevity we will omit this unless it needs to be emphasized.
In this section, we derive appropriate upper bounds for ‖∇ϕk‖2 to use in (13). We begin with∇zϕk:
‖∇zϕk‖2 = ∥∥∥ n+1∑ i=1 yki ∥∥∥2 ≤ 2‖ykn+1‖2 + 2∥∥∥ n∑ i=1 yki ∥∥∥2 = 2∥∥B(xkn+1) + ek∥∥2 + 2∥∥∥ n∑ i=1 yki ∥∥∥2 ≤ 4‖B(xkn+1)‖2 + 2
∥∥∥ n∑ i=1 yki ∥∥∥2 + 4‖ek‖2. Now next take expectations with respect to Fk and Ek, and use the bound on the variance of the noise in (11), obtaining
E [ ‖∇zϕk‖2|Fk, Ek ] ≤ E [ 4‖B(xkn+1)‖2 + 2 ∥∥∥ n∑ i=1 yki ∥∥∥2 + 4‖ek‖2 ∣∣∣ Fk, Ek]
≤ 4(N + 1)‖B(xkn+1)‖2 + 2 ∥∥∥ n∑ i=1 yki ∥∥∥2 + 4N, where we have used that yki is Fk-measurable for i ∈ 1..n. Thus, taking expectations over Ek conditioned on Fk yields
E [ ‖∇zϕk‖2|Fk ] ≤ 4(N + 1)E[‖B(xkn+1)‖2|Fk] + 2 ∥∥∥ n∑ i=1 yki ∥∥∥2 + 4N. (22) We will now bound the two terms on the right side of (22).

heading: C.3.1 FIRST TERM IN (22)
text: First, note that
‖B(zk)‖2 = ‖B(zk)−B(z∗) +B(z∗))‖2
≤ 2‖B(zk)−B(z∗)‖2 + 2‖B(z∗)‖2
≤ 2L2‖zk − z∗‖2 + 2‖B(z∗)‖2
≤ 2L2‖pk − p∗‖2 + 2‖B(z∗)‖2. (23)
Now, returning to the first term on the right of (22), we have
‖B(xkn+1)‖2 = ‖B(zk) +B(xkn+1)−B(zk)‖2
≤ 2‖B(zk)‖2 + 2‖B(xkn+1)−B(zk)‖2
≤ 2‖B(zk)‖2 + 2L2‖xkn+1 − zk‖2 ≤ 4L2‖pk − p∗‖2 + 4‖B(z∗)‖2 + 2L2‖xkn+1 − zk‖2 (24)
where we have used (23) to obtain (24).
For the third term in (24), we have from the calculation on line 7 of the algorithm that
xkn+1 − zk = −ρk(rk − wkn+1) = −ρk(B(zk) + k − wkn+1),
and therefore
‖xkn+1 − zk‖2 = ρ2k‖B(zk) + k − wkn+1‖2
≤ ρ2‖B(zk) + k − wkn+1‖2
≤ 3ρ2(‖B(zk)‖2 + ‖ k‖2 + ‖wkn+1‖2).
We next take expectations conditioned on Fk and use the noise variance bound (10) to obtain E [ ‖xkn+1 − zk‖2 | Fk ] ≤ E [ 3ρ2 ( ‖B(zk)‖2 + ‖ k‖2 + ‖wkn+1‖2 ) | Fk ] ≤ 3ρ2 ( (N + 1)‖B(zk)‖2 + ‖wkn+1‖2 +N ) .
Therefore E [ ‖xkn+1 − zk‖2 | Fk ] ≤ 6ρ2 ( (N + 1)‖B(zk)‖2 + ‖wkn+1 − w∗n+1‖2 + ‖w∗n+1‖2 ) + 3ρ2N
= 6ρ2 ( 2(N + 1)L2‖pk − p∗‖2 + 2(N + 1)‖B(z∗)‖2
+ ‖wkn+1 − w∗n+1‖2 + ‖B(z∗)‖2 ) + 3ρ2N
≤ 6ρ2 ( 2(N + 1)L2‖pk − p∗‖2 + ‖wkn+1 − w∗n+1‖2 ) + 18ρ2(N + 1)‖B(z∗)‖2 + 3ρ2N
≤ 18ρ2(N + 1) ( (L2 + 1)‖pk − p∗‖2 + ‖B(z∗)‖2 ) + 3ρ2N (25)
where in the equality uses (23) and w∗n+1 = B(z ∗). Combining (24) and (25), we arrive at E [∥∥B(xkn+1)∥∥2 ∣∣∣Fk] ≤ 4L2[1 + 9ρ2(L2 + 1)(N + 1)]‖pk − p∗‖2
+ 4 ( 1 + 9ρ2L2(N + 1) ) ‖B(z∗)‖2 + 6ρ2L2N. (26)

heading: C.3.2 SECOND TERM IN (22)
text: For i ∈ 1..n, line 5 of the algorithm may be rearranged into yki = τ−1(zk − xki ) + wki , so∥∥∥ n∑ i=1 yki ∥∥∥2 = ∥∥∥ n∑ i=1 (τ−1(zk − xki ) + wki ) ∥∥∥2
≤ 2 ∥∥∥τ−1 n∑
i=1
(zk − xki ) ∥∥∥2 + 2∥∥∥ n∑
i=1
wki ∥∥∥2 ≤ 2nτ−2
n∑ i=1 ‖zk − xki ‖2 + 2 ∥∥∥ n∑ i=1 wki ∥∥∥2 ≤ 4n2τ−2‖zk − z∗‖2 + 4nτ−2
n∑ i=1 ‖z∗ − xki ‖2 + 4n n∑ i=1 ‖wki − w∗i ‖2 + 4 ∥∥∥ n∑ i=1 w∗i ∥∥∥2 ≤ 4n2(τ−2 + 1)‖pk − p∗‖2 + 4nτ−2
n∑ i=1 ‖z∗ − xki ‖2 + 4 ∥∥∥ n∑ i=1 w∗i ∥∥∥2. (27)
By the definition of the solution set S in (5), w∗i ∈ Ai(z∗), so z∗ + τw∗i ∈ (I + τAi)(z∗), and since the resolvent is single-valued (Bauschke & Combettes, 2017, Cor. 23.9) we therefore obtain
z∗ = (I + τAi) −1(I + τAi)(z ∗) = JτAi(z ∗ + τw∗i ).
From lines 3 and 4 of the algorithm, we also have xki = JτAi(z k + τwki ) for i ∈ 1..n. Thus, using the nonexpansiveness of the resolvent (Bauschke & Combettes, 2017, Def. 4.1 and Cor. 23.9), we have
n∑ i=1 ‖z∗ − xki ‖2 = n∑ i=1 ∥∥JτAi(zk + τwki )− JτAi(z∗ + τw∗i )∥∥2 ≤
n∑ i=1 ‖zk + τwki − z∗ − τw∗i ‖2
= n∑ i=1 ‖zk − z∗ + τ(wki − w∗i )‖2
≤ 2n‖zk − z∗‖2 + 2τ2 n∑ i=1 ‖wki − w∗i ‖2 ≤ 2(n+ τ2)‖pk − p∗‖2. (28)
Combining (27) and (28) yields∥∥∥ n∑ i=1 yki ∥∥∥2 ≤ 12n2τ−2(n+ τ2)‖pk − p∗‖2 + 4∥∥∥ n∑ i=1 w∗i ∥∥∥2. (29) Combining (26) and (29) with (22) yields
E [ ‖∇zϕk‖2 | Fk ] ≤ 24 [ (1 + 9ρ2)(L2 + 1)2(N + 1)2 + n2τ−2(n+ τ2) ] ‖pk − p∗‖2
+ 16(N + 1) ( 1 + 9ρ2L2(N + 1) ) ‖B(z∗)‖2 + 8 ∥∥∥ n∑ i=1 w∗i ∥∥∥2 + 24ρ2L2(N + 1)N + 4N. (30)

heading: C.3.3 DUAL GRADIENT NORM
text: Considering that ∇ϕk is taken with respect to the subspace P , the gradients with respect to the dual variables are — see for example Eckstein & Svaiter (2009) — for each i ∈ 1..(n+ 1),
‖∇wiϕk‖2 = ∥∥∥xki − 1n+ 1 n+1∑ j=1 xkj ∥∥∥2 = ∥∥∥ 1 n+ 1 n+1∑ j=1 (xki − xkj ) ∥∥∥2
≤ n+1∑ j=1 ‖xki − xkj ‖2 ≤ 2 n+1∑ j=1 ( ‖xki − zk‖2 + ‖zk − xkj ‖2
) Summing this inequality for i ∈ 1..(n+ 1) and collecting terms yields
n+1∑ i=1 ‖∇wiϕk‖2 ≤ 4(n+ 1) n+1∑ i=1 ‖xki − zk‖2,
so taking expectations conditioned on Fk produces
n+1∑ i=1 E[‖∇wiϕk‖2 | Fk] ≤ 4(n+ 1) n+1∑ i=1 E[‖xki − zk‖2 | Fk]
≤ 4(n+ 1)E[‖xkn+1 − zk‖2 | Fk] + 4(n+ 1) n∑ i=1 E[‖xki − zk‖2 | Fk]
≤ 4(n+ 1)E[‖xkn+1 − zk‖2 | Fk]
+ 8(n+ 1) n∑ i=1 E[‖xki − z∗‖2 | Fk] + 8(n+ 1)2‖zk − z∗‖2
≤ 4(n+ 1)E[‖xkn+1 − zk‖2|Fk]
+ 8(n+ 1) n∑ i=1 E[‖xki − z∗‖2|Fk] + 8(n+ 1)2‖pk − p∗‖2
≤ 8(n+ 1) [ 3n+ 2τ2 + 1 + 9ρ2(L2 + 1)(N + 1) ] ‖pk − p∗‖2
+ 72ρ2(n+ 1)(N + 1)‖B(z∗)‖2 + 12ρ2(n+ 1)N, (31)
where the final inequality employs (25) and (28).
All told, using (30) and (31) and simplifying the constants, one obtains
E[‖∇ϕk‖2 | Fk] = E[‖∇zϕk‖2 | Fk] + n+1∑ i=1 E[‖∇wiϕk‖2|Fk]
≤ C1‖pk − p∗‖2 + C2, (32)
where
C1 = 24(1 + 10ρ 2)(n+ 1)(L2 + 1)2(N + 1)2 + 8(n+ 1) ( 2τ2 + 6(n+ 1) + 1 + 3(n+ 1)2τ−2 ) (33)
and
C2 = 16(N + 1) [ 1 + 4ρ2(n+ 1) + 9ρ2L2(N + 1) ] ‖B(z∗)‖2 + 8‖ n∑ i=1 w∗i ‖2
+ 12ρ2N(2L2(N + 1) + n+ 1) + 4N. (34)
C.4 LOWER BOUND FOR ϕk-GAP
Recalling (13), that is,
‖pk+1 − p∗‖2 = ‖pk − p∗‖2 − 2αk(ϕk(pk)− ϕk(p∗)) + α2k‖∇ϕk‖2.
We may use the gradient bound from (32) to obtain
E[‖pk+1 − p∗‖2 | Fk] ≤ (1 + C1α2k)‖pk − p∗‖2 − 2αkE[ϕk(pk)− ϕk(p∗) | Fk] + C2α2k. (35)
We now focus on finding a lower bound for the term E[ϕk(pk) − ϕk(p∗) | Fk], which we call the “ϕk-gap”. Recall that for p = (z, w1, . . . , wn+1),
ϕk(p) = n+1∑ i=1 〈z − xki , yki − wi〉.
For each i ∈ 1..(n+1), define ϕi,k(p) . = 〈z−xki , yki −wi〉. We will call E[ϕi,k(pk)−ϕi,k(p∗) | Fk] the “ϕi,k-gap”. Note that ϕk(p) = ∑n+1 i=1 ϕi,k(p).
C.5 LOWER BOUND FOR ϕi,k-GAP OVER i ∈ 1..n
For i ∈ 1..n, we have from line 5 of the algorithm that
zk − xki = τ(yki − wki ).
Since ϕi,k(pk) = 〈zk − xki , yki − wki 〉, one may conclude that for i ∈ 1..n,
ϕi,k(p k) =
τ 2 ‖yki − wki ‖2 + 1 2τ ‖zk − xki ‖2.
On the other hand, for p∗ ∈ S and i ∈ 1..n, one also has −ϕi,k(p∗) = 〈z∗ − xki , w∗i − yki 〉 ≥ 0 (36)
by the monotonicity of Ai. Therefore, for i ∈ 1..n, it holds that
ϕi,k(p k)− ϕi,k(p∗) ≥
τ 2 ‖yki − wki ‖2 + 1 2τ ‖zk − xki ‖2,
and taking expectations conditioned on Fk leads to
E[ϕi,k(pk)− ϕi,k(p∗) | Fk] ≥ τ
2 ‖yki − wki ‖2 +
1
2τ ‖zk − xki ‖2 (37)
where we have used that xki and y k i are both Fk-measurable for i ∈ 1..n.
C.6 LOWER BOUND FOR ϕn+1,k-GAP
From lines 6-7 of the algorithm, we have
zk − xkn+1 = ρk(B(zk)− wkn+1 + k). Therefore,
ϕn+1,k(p k) = 〈zk − xkn+1, ykn+1 − wkn+1〉 (38)
= 〈zk − xkn+1, B(zk)− wkn+1〉+ 〈zk − xkn+1, ykn+1 −B(zk)〉 = ρk〈B(zk)− wkn+1 + k, B(zk)− wkn+1〉+ 〈zk − xkn+1, ykn+1 −B(zk)〉 = ρk‖B(zk)− wkn+1‖2 + 〈zk − xkn+1, ykn+1 −B(zk)〉+ ρk〈 k, B(zk)− wkn+1〉 (a) = ρk‖B(zk)− wkn+1‖2 + 〈zk − xkn+1, B(xkn+1)−B(zk)〉+ 〈zk − xkn+1, ek〉
+ ρk〈 k, B(zk)− wkn+1〉 ≥ ρk‖B(zk)− wkn+1‖2 − L‖zk − xkn+1‖2 + 〈zk − xkn+1, ek〉
+ ρk〈 k, B(zk)− wkn+1〉 = ρk‖B(zk)− wkn+1‖2 − L‖ρk(B(zk)− wkn+1 + k)‖2 + 〈zk − xkn+1, ek〉
+ ρk〈 k, B(zk)− wkn+1〉 = ρk‖B(zk)− wkn+1‖2 − ρ2kL‖B(zk)− wkn+1 + k‖2 + 〈zk − xkn+1, ek〉
+ ρk〈 k, B(zk)− wkn+1〉 = ρk(1− ρkL)‖B(zk)− wkn+1‖2 − ρ2kL‖ k‖2 + 〈zk − xkn+1, ek〉
+ ρk(1− 2ρkL)〈 k, B(zk)− wkn+1〉, (39) where equality (a) uses line 8 of the algorithm and the inequality employs the Cauchy-Schwartz inequality followed by Lipschitz continuity of B.
On the other hand,
−ϕn+1,k(p∗) = 〈z∗ − xkn+1, w∗n+1 − ykn+1〉 = 〈z∗ − xkn+1, B(z∗)−B(xki )〉+ 〈xkn+1 − z∗, ek〉 ≥ 〈xkn+1 − z∗, ek〉, (40)
where the second equality uses line 8 of the algorithm and the inequality follows from the monotonicity of B.
Combining (39) and (40) yields ϕn+1,k(p
k)− ϕn+1,k(p∗) ≥ ρk(1− ρkL)‖B(zk)− wkn+1‖2 + ρk(1− 2ρkL)〈 k, B(zk)− wkn+1〉 + 〈zk − xkn+1, ek〉+ 〈xkn+1 − z∗, ek〉 − ρ2kL‖ k‖2
= ρk(1− ρkL)‖B(zk)− wkn+1‖2 − ρ2kL‖ k‖2
+ ρk(1− 2ρkL)〈 k, B(zk)− wkn+1〉+ 〈zk − z∗, ek〉. (41)
Now, if we take expectations conditioned on Fk and use (9), we obtain E [ 〈zk − z∗, ek〉 ∣∣ Fk] = 〈zk − z∗,E[ek | Fk]〉 = 0. (42) Similarly, (9) also yields
E [ 〈 k, B(zk)− wkn+1〉 ∣∣ Fk] = 〈E[ k|Fk], B(zk)− wkn+1〉 = 0. (43) Thus, using (42) and (43) and taking expectations of (41) yields
E[ϕn+1,k(pk)− ϕn+1,k(p∗) | Fk] ≥ ρk(1− ρkL)‖B(zk)− wkn+1‖2 − ρ2kLE[‖ k‖2|Fk] ≥ ρk(1− ρL)‖B(zk)− wkn+1‖2 − ρ2kNL(1 + ‖B(zk)‖2),
(44) where in the second inequality we used (12) and the noise variance bound (10). Recall from (12) that 1− ρL > 0. Next, we remark that
‖B(zk)‖2 = ‖B(zk)−B(z∗) +B(z∗)‖2
≤ 2L2‖zk − z∗‖2 + 2‖B(z∗)‖2 ≤ 2L2‖pk − p∗‖2 + 2‖B(z∗)‖2. Substituting this inequality into (44) yields
E[ϕn+1,k(pk)− ϕn+1,k(p∗)|Fk] ≥ ρk(1− ρL)‖B(zk)− wkn+1‖2
− 2ρ2kNL3‖pk − p∗‖2 − ρ2kNL(1 + 2‖B(z∗)‖2). (45)
Finalizing the lower bound on the ϕk-gap Summing (37) over i ∈ 1..n and using (45) yields
E[ϕk(pk)− ϕk(p∗)|Fk] = n+1∑ i=1 E[ϕi,k(pk)− ϕi,k(p∗)|Fk]
≥ τ 2 n∑ i=1 ‖yki − wki ‖2 + 1 2τ n∑ i=1 ‖zk − xki ‖2
+ ρk(1− ρL)‖B(zk)− wkn+1‖2 − 2ρ2kNL3‖pk − p∗‖2 − ρ2kNL(1 + 2‖B(z∗)‖2). (46)

heading: C.7 ESTABLISHING STOCHASTIC QUASI-FEJER MONOTONICITY
text: Returning to (35), E[‖pk+1 − p∗‖2 | Fk] ≤ (1 + C1α2k)‖pk − p∗‖2 − 2αkE[ϕk(pk)− ϕk(p∗) | Fk] + C2α2k, we may now substitute (46) for the expectation on the right-hand side. First, define
Tk . = τ
ρ n∑ i=1 ‖yki − wki ‖2 + 1 ρτ n∑ i=1 ‖zk − xki ‖2 + 2(1− ρL)‖B(zk)− wkn+1‖2,
after which we may use (46) in (35) to yield E[‖pk+1 − p∗‖2 | Fk] ≤ (1 + C1α2k + C3αkρ2k)‖pk − p∗‖2 − αkρkTk + C2α2k + C4αkρ2k (47) where C1 and C2 are defined as before in (33) and (34) and C3 = 4NL 3 (48)
C4 = 2NL(1 + 2‖B(z∗)‖2). (49) This completes the proof of Lemma 3.

heading: C.8 A CONVERGENCE LEMMA
text: Before establishing almost-sure convergence, we need the following lemma to derive convergence of the iterates from convergence of Tk defined above. Note that a more elaborate result would be needed in an infinite-dimensional setting.
Lemma 4. For deterministic sequences zk ∈ R(n+1)d, {(wki ) n+1 i=1 } ∈ P , and {(xki , yki ) n+1 i=1 } ∈ R2(n+1)d, suppose that yki ∈ Ai(xki ) for i ∈ 1..n, ∑n+1 i=1 w k i = 0,
ξ1 n∑ i=1 ‖yki − wki ‖2 + ξ2 n∑ i=1 ‖zk − xki ‖2 + ξ3‖B(zk)− wkn+1‖2 → 0 (50)
for scalars ξ1, ξ2, ξ3 > 0, and pk . = (zk, wk1 , . . . , w k n+1)→ p̂ . = (ẑ, ŵ1, . . . , ŵn+1). Then p̂ ∈ S.
Proof. Fix any i ∈ {1, . . . , n}. Since ‖yki − wki ‖ → 0 by (50) and wki → ŵi, we also have yki → ŵi. Similarly, (50) also implies that ‖zk − xki ‖ → 0, so from zk → ẑ we also have xki → ẑ. Since yki ∈ Ai(xki ) and (xki , yki ) → (ẑ, ŵi), (Bauschke & Combettes, 2017, Prop. 20.37) implies ŵi ∈ Ai(ẑ). Since i was arbitrary, the preceding conclusions hold for i ∈ 1..n. Now, (50) also implies that ‖B(zk)− wkn+1‖ → 0. Therefore, since wkn+1 → ŵn+1, we also have B(zk) → ŵn+1. Much as before, since (zk, B(zk)) → (ẑ, ŵn+1), we may apply (Bauschke & Combettes, 2017, Prop. 20.37) to conclude that that ŵn+1 = B(ẑ).
Since the linear subspace P defined in (6) must be closed, the limit (ẑ, ŵ1, . . . , ŵn+1) of {(zk, wk1 , . . . , wkn+1)} ⊂ P must be in P , hence ∑n+1 i=1 ŵi = 0. Thus, the point p̂ = (ẑ, ŵ1, . . . , ŵn+1) satisfies ŵi ∈ Ai(ẑ) for i ∈ 1..n, ŵn+1 = B(ẑ), and∑n+1 i=1 ŵi = 0. These are the three conditions defining membership in S from (5), so p̂ ∈ S.

heading: C.9 FINISHING THE PROOF OF THEOREM 1
text: Given ∑ k α 2 k < ∞, and ∑ αkρ 2 k < ∞, (47) satisfies the conditions of Stochastic Quasi-Fejer Monotonicity as given in Lemma 2. By applying Lemma 2, we conclude that there exist Ω1,Ω2,Ω3 such that P [Ωi] = 1 for i = 1, 2, 3 and
1. for all v ∈ Ω1 ∞∑ k=1 αkρkTk(v) <∞, (51)
2. for all v ∈ Ω2, and p∗ ∈ S , ‖pk(v)− p∗‖ converges to a finite nonnegative random-variable, 3. for all v ∈ Ω3, pk(v) remains bounded.
Since ∑∞ k=1 αkρk =∞, (51) implies that for all v ∈ Ω1 there exists a subsequence qk(v) such that
Tqk(v) → 0. (52)
Let Ω′ = Ω1 ∩ Ω2 ∩ Ω3 and note that P [Ω′] = 1. Choose v ∈ Ω′. Since pk(v) remains bounded, so does pqk(v)(v) for qk(v) defined above in (52). Thus there exists a subsequence rk(v) ⊆ qk(v) and p̂(v) ∈ R(n+2)d such that prk(v)(v)→ p̂(v). But since Tqk(v) → 0, it also follows that Trk(v) → 0, that is,
τ
ρ n∑ i=1 ‖yrk(v)i (v)− w rk(v) i (v)‖ 2 + 1 ρτ n∑ i=1 ‖zrk(v)(v)− xrk(v)i (v)‖ 2
+ 2(1− ρL)‖B(zrk(v)(v))− wrk(v)n+1 (v)‖2 → 0.
We then have from Lemma 4 that p̂(v) ∈ S.
Since prk(v)(v)→ p̂(v), it follows that ‖prk(v)(v)− p̂(v)‖ → 0. But since p̂(v) ∈ S , ‖pk(v)− p̂(v)‖ converges by point 2 above. Thus
lim k→∞ ‖pk(v)− p̂(v)‖ = lim k→∞ ‖prk(v)(v)− p̂(v)‖ = 0.
Therefore pk(v)→ p̂(v) ∈ S. Thus there exists p̂ ∈ S such that pk → p̂ a.s., which completes the proof of Theorem 1.

heading: C.10 TWO ADDITIONAL RESULTS
text: In this section, we prove two additional useful results about SPS. First, that xki → ẑ (a.s.) for i = 1, . . . , n. Second, that Gk → 0 (a.s.). Note that
xki = JτAi(z k + τwki )
and since zk and wki convergence a.s., so does x k i . Consider the subsequence qk(v) such that (52) holds. Then zqk(v) − xqk(v)i → 0
thus x qk(v) i → ẑ.
Since xki converges to some limit (a.s.), that limit must be ẑ.
Recall that
Gk . = ∑n i=1 ‖yki − wki ‖2 + ∑n i=1 ‖zk − xki ‖2 + ‖B(zk)− wkn+1‖2.
We have shown that zk and xki share the same limit for i = 1, . . . , n (a.s.). Therefore z k − xki → 0 (a.s.). Since yki − wki = τ−1(zk − xki ),
it follows that yki − wki → 0 (a.s.) for i = 1, . . . , n. Therefore
Gk → ‖B(ẑ)− ŵn+1‖2.
But since (z, ŵ1, . . . , ŵn+1) ∈ S, ŵn+1 = B(ẑ) implying that Gk → 0 (a.s.).

heading: D PROOF OF LEMMA 1
text: If Gk = 0, then
∀i = 1, . . . , n : yki = wki and zk = xki . (53)
Since yki ∈ Ai(xki ) for i = 1, . . . , n, (53) implies that that
∀i ∈ 1..n : wki ∈ Ai(zk). (54)
Furthermore Gk = 0 also implies that wkn+1 = B(z k). Finally, since ∑n+1 i=1 w k i = 0, we have that
(zk, wk1 , . . . , w k n+1) ∈ S.
Conversely, suppose (zk, wk1 , . . . , w k n+1) ∈ S. The definition of S implies that B(zk) = wkn+1 and furthermore that wki ∈ Ai(zk) for i ∈ 1..n. For any i ∈ 1..n, considering line 3 of Algorithm 1, we may write tki = z
k + τwki ∈ (I + τAi)(zk), implying zk ∈ (I + τAi)−1(tki ). But since the resolvent JτAi = (I + τAi)
−1 is single-valued (Bauschke & Combettes, 2017, Prop. 23.8), we must have zk = (I + τAi)−1(tki ). Thus, by line 4, we have x k i = z k. We may also derive from line 5 that
yki = τ −1(tki − xki ) = τ−1(zk + τwki − zk) = wki .
Thus, since xki = z k and yki = w k i for i = 1, . . . , n and w k n+1 = B(z k), we have that Gk = 0.

heading: E PROOF OF THEOREM 2
text: In addition to the proof, we provide a more detailed statement of the theorem: Theorem 3. Fix the total iterations K ≥ 1 of Algorithm 1 and set
∀k = 1, . . . ,K : ρk = ρ . = min { K−1/4, 1
2L
} (55)
∀k = 1, . . . ,K : αk = α . = Cfρ 2 (56)
for some Cf > 0. Suppose (9)-(11) hold. Then for any p∗ ∈ S,
1
K K∑ j=1 E[Gj ] ≤ 8L3 exp (Cf (C1 + C3)) Cf min{τ, τ−1}K ( ‖p1 − p∗‖2 + CfC2 + C4 CfC1 + C3 ) for K < (2L)4 (57)
1
K K∑ j=1 E[Gj ] ≤ exp (Cf (C1 + C3)) Cf min{τ, τ−1}K1/4 ( ‖p1 − p∗‖2 + CfC2 + C4 CfC1 + C3 ) for K ≥ (2L)4. (58)
where Gk is the approximation residual defined in (14), and C1, C2, C3, C4 are the nonegative constants defined in (33), (34), (48), and (49), respectively. Therefore,
1
K K∑ j=1 E[Gj ] = O(K−1/4).
Proof. Fix αk = α and ρk = ρ, where α and ρ are the respective right-hand sides of (55)-(56). Lemma 3 implies that (21) so long as (9)-(11) hold and the stepsize ρ satisfies ρ < L−1. Since
ρ = min { K−1/4, 1
2L
} ≤ 1
2L ,
we conclude that (21) applies.
Rewriting (21) with αk = α and ρk = ρ, we have
E[‖pk+1 − p∗‖2 | Fk] ≤ (1 + C1α2 + C3αρ2)‖pk − p∗‖2 − αρTk + C2α2 + C4αρ2. Therefore, taking expectations over Fk, we have
E‖pk+1 − p∗‖2 ≤ (1 + C1α2 + C3αρ2)E‖pk − p∗‖2 − αρETk + C2α2 + C4αρ2. (59) Recall that
Tk . = τ
ρ n∑ i=1 ‖yki − wki ‖2 + 1 ρτ n∑ i=1 ‖zk − xki ‖2 + 2(1− ρL)‖B(zk)− wkn+1‖2,
where for the first two terms we have simply set ρ = ρ because the stepsize is constant. However, for the final term, we will still use an upper bound, ρ, on ρ. In the current setting, we know that ρ ≤ (1/2)L−1 and therefore we may set ρ = (1/2)L−1. Thus 1− ρL = 1/2, leading to
ρETk = τ n∑ i=1 E‖yki − wki ‖2 + τ−1 n∑ i=1 E‖zk − xki ‖2 + ρE‖B(zk)− wkn+1‖2.
Let
Uk . = E‖B(zk)− wkn+1‖2 Wk . = τ n∑ i=1 E‖yki − wki ‖2 + τ−1 n∑ i=1 E‖zk − xki ‖2,
so that
ρETk = ρUk +Wk,
and also let
Vk . = E‖pk − p∗‖2.
Using these definitions in (59) we write
Vk+1 ≤ (1 + C1α2 + C3αρ2)Vk − αρUk − αWk + C2α2 + C4αρ2. Therefore,
Vk+1 + αρUk + αWk ≤ (1 + C1α2 + C3αρ2)Vk + C2α2 + C4αρ2
⇐⇒ Vk+1 + αρ k∑ j=1 Uj + α k∑ j=1 Wj ≤ (1 + C1α2 + C3αρ2)Vk + αρ k−1∑ j=1 Uj + α k−1∑ j=1 Wj
+ C2α 2 + C4αρ 2
≤ (1 + C1α2 + C3αρ2) Vk + αρ k−1∑ j=1 Uj + α k−1∑ j=1 Wj  + C2α 2 + C4αρ 2,
where we have used that Uk,Wk ≥ 0. Letting
Rk = Vk + αρ k−1∑ j=1 Uj + α k−1∑ j=1 Wj ,
we then have
Rk+1 ≤ (1 + C1α2 + C3αρ2)Rk + C2α2 + C4αρ2, which implies
Rk+1 ≤ (1 + C1α2 + C3αρ2)kR1 + (C2α2 + C4αρ2) k∑ j=1 (1 + C1α 2 + C3αρ 2)k−j .
Now, k∑ j=1 (1 + C1α 2 + C3αρ 2)k−j = k−1∑ j=0 (1 + C1α 2 + C3αρ 2)j
= (1 + C1α
2 + C3αρ 2)k − 1
(1 + C1α2 + C3αρ2)− 1
= (1 + C1α
2 + C3αρ 2)k − 1
C1α2 + C3αρ2
≤ (1 + C1α 2 + C3αρ 2)k
C1α2 + C3αρ2 .
Therefore,
Rk+1 ≤ (1 + C1α2 + C3αρ2)k ( R1 + C2α 2 + C4αρ 2
C1α2 + C3αρ2
) .
Fix the number of iterations K ≥ 1. Now
ρ = min { K−1/4, 1
2L } ≤ 1 K1/4 ≤ 1.
Therefore,
αρ K∑ j=1 (Uj +Wj) ≤ αρ K∑ j=1 Uj + α K∑ j=1 Wj
≤ RK+1 ≤ (1 + C1α2 + C3αρ2)K ( R1 + C2α 2 + C4αρ 2
C1α2 + C3αρ2
) .
Dividing through by αρK, we obtain
1
K K∑ j=1 (Uj +Wj) ≤ (1 + C1α 2 + C3αρ 2)K αρK ( R1 + C2α 2 + C4αρ 2 C1α2 + C3αρ2 ) , (60)
and since α = Cfρ2, we also have
C2α 2 + C4αρ 2
C1α2 + C3αρ2 = CfC2 + C4 CfC1 + C3 .
Furthermore,
ρ ≤ K− 14 =⇒ α ≤ CfK− 1 2 .
Substituting these into (60) yields
1
K K∑ j=1 (Uj +Wj) ≤
( 1 +
Cf (CfC1+C3) K )K αρK ( R1 + CfC2 + C4 CfC1 + C3 )
≤ exp(Cf (CfC1 + C3)) αρK
( R1 +
CfC2 + C4 CfC1 + C3
) , (61)
where we have used that for any t ≥ 0, 1 + t/K ≤ et/K , so therefore (1 + t/K)K ≤ et.
The worst-case rates in terms of K occur when ρ = K−1/4 and α = CfK−1/2. This is the case when K ≥ (2L)4. Substituting these into the denominator yields, for K ≥ (2L)4, that
1
K K∑ j=1 (Uj +Wj) ≤ exp(Cf (CfC1 + C3)) CfK1/4 ( R1 + CfC2 + C4 CfC1 + C3 ) .
Thus, since Gk ≤ max{τ, τ−1} (Uk +Wk) , we obtain
1
K K∑ j=1 E[Gj ] ≤ max{τ, τ−1} exp (Cf (CfC1 + C3)) CfK1/4 ( ‖p1 − p∗‖2 + CfC2 + C4 CfC1 + C3 ) ,
which is (58).
When K < (2L)4, (57) can similarly be obtained by substituting ρ = (2L)−1 and α = Cf (2L)−2 into (61).

heading: F APPROXIMATION RESIDUALS
text: In this section we derive the approximation residual used to assess the performance of the algorithms in the numerical experiments. This residual relies on the following product-space reformulation of (1).

heading: F.1 PRODUCT-SPACE REFORMULATION AND RESIDUAL PRINCIPLE
text: Recall (1), the monotone inclusion we are solving:
Find z ∈ Rd : 0 ∈ n∑ i=1 Ai(z) +B(z).
In this section we demonstrate a “product-space" reformulation of (1) which allows us to rewrite it in a standard form involving just two operators, one maximal monotone and the other monotone and Lipschitz. This approach was pioneered in (Briceño-Arias & Combettes, 2011; Combettes & Pesquet, 2012). Along with allowing for a simple definition of an approximation residual as a measure of approximation error in solving (1), it allows one to apply operator splitting methods originally formulated for two operators to problems such as (1) for any finite n.
Observe that solving (1) is equivalent to
Find (w1, . . . , wn, z) ∈ R(n+1)d : wi ∈ Ai(z), i ∈ 1..n
0 ∈ n∑ i=1 wi +B(z).
This formulation resembles that of the extended solution set S used in projective spitting, as given in (5), except that it combines the final two conditions in the definition of S , and thus does not need the final dual variable wn+1. From the definition of the inverse of an operator, the above formulation is equivalent to
Find (w1, . . . , wn, z) ∈ R(n+1)d : 0 ∈ A−1i (wi)− z, i ∈ 1..n
0 ∈ n∑ i=1 wi +B(z).
These conditions are in turn equivalent to finding (w1, . . . , wn, z) ∈ R(n+1)d such that 0 ∈ A (w1, . . . , wn, z) + B(w1, . . . , wn, z), (62)
where A is the set-valued map
A (w1, . . . , wn, z) 7→ A−11 (w1)×A −1 2 (w2)× . . .×A−1n (wn)× {0} (63)
and B is the single-valued operator
B(w1, . . . , wn, z) 7→  0 · · · 0 −I ... . . . ...
... 0 · · · 0 −I I · · · I 0   w1 ... wn z +  0 ... 0 B(z)  . (64) It is easily established that B is maximal monotone and Lipschitz continuous, while A is maximal monotone. Letting T .= A + B, it follows from (Bauschke & Combettes, 2017, Proposition 20.23) that T is maximal monotone. Thus, we have reformulated (1) as the monotone inclusion 0 ∈ T (q) for q in the product space R(n+1)d. A vector z ∈ Rd solves (1) if and only if there exists (w1, . . . , wn) ∈ Rnd such that 0 ∈ T (q), where q = (w1, . . . , wn, z). For any pair (q, v) such that v ∈ T (q), ‖v‖2 represents an approximation residual for q in the sense that v = 0 implies q is a solution to (62). One may take ‖v‖2 as a measure of the error of q as an approximate solution to (62), and it can only be 0 if q is a solution. Given two approximate solutions q1 and q2 with certificates v1 ∈ T (q1) and v2 ∈ T (q2), we will treat q1 as a “better” approximate solution than q2 if ‖v1‖2 < ‖v2‖2. Doing so is somewhat analogous to the practice, common in optimization, of using the gradient ‖∇f(x)‖2 as a measure of quality of an approximate minimizer of some differentiable function f . However, note that since T (q1) is a set, there may exist elements of T (q1) with smaller norm than v1. Thus any given certificate v1 only corresponds to an upper bound on dist2(0,T (q1)).

heading: F.2 APPROXIMATION RESIDUAL FOR PROJECTIVE SPLITTING
text: In SPS (Algorithm 1), for i ∈ 1..n, the pairs (xki , yki ) are chosen so that yki ∈ Ai(xki ). This can be seen from the definition of the resolvent. Thus xki ∈ A −1 i (y k i ). Observe that
vk . =  xk1 − zk ... xkn − zk
B(zk) + ∑n i=1 y k i  ∈ T (yk1 , . . . , ykn, zk). (65) The approximation residual for SPS is thus
Rk . = ‖vk‖2 = n∑ i=1 ‖zk − xki ‖2 + ∥∥B(zk) + n∑ i=1 yki ∥∥2 (66)
which is an approximation residual for (yk1 , . . . , y k n, z k) in the sense defined above. We may relate Rk to the approximation residual Gk for SPS from Section 5 as follows:
Rk = n∑ i=1 ‖zk − xki ‖2 + ∥∥∥∥∥B(zk) + n∑ i=1 yki ∥∥∥∥∥ 2
= n∑ i=1 ‖zk − xki ‖2 + ∥∥∥∥∥B(zk) + n∑ i=1 yki − n+1∑ i=1 wki ∥∥∥∥∥ 2 ≤ n∑ i=1 ‖zk − xki ‖2 + 2‖B(zk)− wkn+1‖2 + 2 ∥∥∥∥∥ n∑ i=1 (yki − wki ) ∥∥∥∥∥ 2
≤ n∑ i=1 ‖zk − xki ‖2 + 2‖B(zk)− wkn+1‖2 + 2n n∑ i=1 ∥∥yki − wki ∥∥2 ≤ 2nGk
where in the second equality we have used the fact that ∑n+1 i=1 w k i = 0. Thus, Rk has the same convergence rate as Gk given in Theorem 2.
Note that while the certificate given in (65) focuses on the primal iterate zk, it may be changed to focus on any xki for i = 1, . . . , n, by using
vki . =  xk1 − xki ... xkn − xki
B(xki ) + ∑n i=1 y k i  ∈ T (yk1 , . . . , ykn, xki ). The approximation residual ‖vki ‖2 may also be shown to have the same rate as Gk by following similar derivations to those above for Rk.

heading: F.3 TSENG’S METHOD
text: Tseng’s method (Tseng, 2000) can be applied to (62), resulting in the following recursion with iterates qk, q̄k ∈ R(n+1)d:
q̄k = JαA (q k − αB(qk)) (67) qk+1 = q̄k + α ( B(qk)−B(q̄k) ) , (68)
where A and B are defined in (63) and (64). The resolvent of A may be readily computed from the resolvents of the Ai using Moreau’s identity (Bauschke & Combettes, 2017, Proposition 23.20).
Analogous to SPS, Tseng’s method has an approximation residual, which in this case is an element of T (q̄k). In particular, using the general properties of resolvent operators as applied to JαA , we have
1 α (qk − q̄k)−B(qk) ∈ A (q̄k).
Also, rearranging (68) produces 1
α (q̄k − qk+1) + B(qk) = B(q̄k).
Adding these two relations produces 1
α (qk − qk+1) ∈ A (q̄k) + B(q̄k) = T (q̄k)
Therefore,
RTsengk . =
1
α2 ‖qk − qk+1‖2
represents a measure of the approximation error for Tseng’s method equivalent to Rk defined in (66) for SPS.

heading: F.4 FRB
text: The forward-reflected-backward method (FRB) (Malitsky & Tam, 2020) is another method that may be applied to the splitting T = A + B for A and B as defined in (63) and (64). Doing so yields recursion
qk+1 = JαA ( qk − α ( 2B(qk)−B(qk−1) )) .
Following similar arguments to those for Tseng’s method, it can be shown that
vkFRB . =
1
α
( qk−1 − qk ) + B(qk) + B(qk−2)− 2B(qk−1) ∈ T (qk).
Thus, FRB admits the following approximation residual equivalent to Rk for SPS:
RFRBk . = ‖vkFRB‖2.
Finally, we remark that the stepsizes used in both the Tseng and FRB methods can be chosen via a linesearch procedure that we do not detail here.

heading: F.5 STOCHASTIC TSENG METHOD
text: The stochastic version of Tseng’s method of (Böhm et al., 2020) (S-Tseng) may be applied to the inclusion 0 ∈ A (q) + B(q), since the operator A may be written as a subdifferential. However, unlike the deterministic Tseng method, it does not produce a valid residual. Note also that S-Tseng outputs an ergodic sequence qkerg. To construct a residual for the ergodic sequence, we compute a deterministic step of Tseng’s method according to (67)-(68), starting at qkerg. That is, letting
q̄k = JαA (q k erg −B(qkerg))
qk+1 = q̄k + α(B(qkerg)−B(q̄k)), we can then compute essentially the same residual as in Section F.3,
RS-Tsengk . =
1
α2 ‖qkerg − qk+1‖2.
To construct the stochastic oracle for S-Tseng, we assumed B(z) = 1m ∑m i=1Bi(z). Then we used
B̃(w1, . . . , wn, z) 7→  0 · · · 0 −I ... . . . ...
... 0 · · · 0 −I I · · · I 0   w1 ... wn z + 
0 ... 0
1 |B| ∑ j∈BBj(z)  . (69) for some minibatch B ∈ {1, . . . ,m}.
F.6 VARIANCE-REDUCED FRB
The FRB-VR method of Alacaoglu et al. (2021) can also be applied to 0 ∈ A (q) + B(q), using the same stochastic oracle B̃ defined in (69). if we let the iterates of FRB-VR be (qk, pk), then line 4 of Algorithm 1 of Alacaoglu et al. (2021) can be written as
q̂k = qk − τ(B(pk) + B̃(qk)− B̃(pk)) (70) qk+1 = JτA (q̂ k). (71)
Once again, the method does not directly produce a residual, but one can be developed from the algorithm definition as follows: (71) yields τ−1(q̂k − qk+1) ∈ A (qk+1) and hence
τ−1(q̂k − qk+1) + B(qk+1) ∈ (A + B)(qk+1). Therefore we use the residual
RFRB-VRk = ‖τ−1(q̂k − qk+1) + B(qk+1)‖2.
Figure 1 plots Rk for SPS, R Tseng k for Tseng’s method, R FRB k for FRB, R S-Tseng k for S-Tseng, and RFRB-VRk for FRB-VR.

heading: F.7 BENEFITS AND DRAWBACKS OF THE PRODUCT SPACE REFORMULATION
text: The main benefit of the product space reformulation (PSR) is that it allows one to use familiar 2-operator splitting schemes for solving 0 ∈ A (q) + B(q) to solve the more complicated recursion (1). However, one drawback of this approach is that the operator B, defined in (64), combines a skew-symmetric consensus matrix with the Lipschitz operator B. Treating B as a single operator necessitates using a single stepsize for both of its constituent operators, but the B component will generally have a much larger Lipschitz constant than the skew part, necessitating a smaller stepsize than is ideal for the skew operator. This difficulty can be countered by using different stepsizes for the primal and dual components, but that strategy introduces additional tuning parameters. In other works, methods based on PSR have exhibited slower convergence than deterministic projective splitting methods (Johnstone & Eckstein, 2021; 2020b). However, in our experiments in Section 7, the performance is comparable.
G VARIATIONAL INEQUALITIES
For a mapping B : Rd → Rd and a closed and convex set C, the variational inequality problem (Harker & Pang, 1990) is to find z∗ ∈ C such that
B(z∗)>(z − z∗) ≥ 0,∀z ∈ C. (72) Consider the normal cone mapping discussed in Section 2 and defined as
NC(x) . = {g : g>(y − x) ≤ 0 ∀y ∈ C}
It is easily seen that (72) is equivalent to finding z∗ such that −B(z∗) ∈ NC(z∗). Hence, if B is monotone, (72) is equivalent to the monotone inclusion 0 ∈ B(z∗) +NC(z∗). (73) Thus, monotone variational inequalities are a special case of monotone inclusions with two operators, one of which is single-valued and the other is the normal cone map of the constraint set C. As a consequence, methods for monotone inclusions can be used to solve monotone variational inequality problems. The reverse, however, may not be true. For example, the analysis of the extragradient method (Korpelevich, 1977) relies on the second operatorNC in (73) being a normal cone, as opposed to a more general monotone operator. We are not aware of any direct extension of the extragradient method’s analysis allowing a more general resolvent to be used in place of the projection map corresponding to NC .
The Restricted Gap Function There is a disadvantage to pursuing convergence rates based on variational inequalities (as in Böhm et al. (2020) and Alacaoglu et al. (2021)) rather than monotone inclusions. Convergence rate analyses for variational inequalities focus on the gap function:
GC(z) . = sup z′∈C B(z′)>(z − z′). (74)
It can be shown that GC(z) ≥ 0 and GC(z) = 0 if and only if z solves (72). However, (74) is meaningless for most problems, since unless C is compact, GC(z) is typically equal to +∞ for any nonsolution (Diakonikolas, 2020). Thus researchers instead focus on the restricted gap function (Nesterov, 2007)
GC2(z) . = sup z′∈C2 B(z′)>(z − z′). (75)
where C2 is an arbitrary compact set. However, now the results are only meaningful over the set C2. Thus, C2 must be chosen large enough so that the iterates of the algorithm remain in the interior of C2 (Böhm et al., 2020). Further, the convergence rate bound depends on the diameter of C2. For some algorithms (Mokhtari et al., 2020) a valid set is provided which bounds the iterates. However Böhm et al. (2020) and Alacaoglu et al. (2021) do not provide one, although in principle it could be done so long as the ergodic sequence can be bounded almost-surely. Thus, the convergence rates depending on (75) in Böhm et al. (2020) and Alacaoglu et al. (2021) are somewhat incomplete in that they depend on unknown constants.
In contrast, rates based on the approximation residual in the monotone inclusion setting, including ours given in (57)–(58), completely avoid this pitfall. There is no need to select a compact set containing the algorithm’s iterates and the constants in our rates are all explicit or depend on standard quantities such as the initial distance to a solution.

heading: H MEMORY-SAVING TECHNIQUE FOR SPS
text: The variables tki , x k i , and y k i on lines 3-5 of SPS are stored in variables t, x and y. Another two variables x̄ and ȳ keep track of ∑n i=1 x k i and ∑n i=1 y k i . The dual variables are stored as wi for i ∈ 1..n and the primal variable as z. Once x = xki is computed, the ith dual variable wi can be partially updated as wi ← wi−αkx. Once all the operators have been processed, the update for each dual variable may be completed via wi ← wi + αk(n+ 1)−1x̄. Also, the primal update is computed as z ← z − αkȳ. During the calculation loop for the xki , yki , the terms in approximation residual Rk may also be accumulated one by one. The total total number of vector elements that must be stored is (n+ 7)d.

heading: I ADDITIONAL INFORMATION ABOUT THE NUMERICAL EXPERIMENTS
text: We solve the following convex-concave min-max problem:
min β∈Rd λ∈R max γ∈Rm
{ λ(δ − κ) + 1
m m∑ i=1 Ψ(〈x̂i, β〉) + 1 m m∑ i=1 γi(ŷi〈x̂i, β〉 − λκ) + c‖β‖1 } s.t. ‖β‖2 ≤ λ/(LΨ + 1) ‖γ‖∞ ≤ 1.
(76)
This model is identical to that of (Yu et al., 2021, Thm. 4.3) except for the addition of the `1 regularization term c‖β‖1, where c ≥ 0 is a given constant. The goal is to learn the model weights β from a training dataset of m feature vectors x̂i and corresponding labels ŷi. Rather than computing the expected loss over the training set, the formulation uses, for each β, the worst possible distribution within a Wasserstein-metric ball around the empirical distribution of the {(x̂i, ŷi)}, with the parameter δ ≥ 0 giving the diameter of the ball and the parameter κ ≥ 0 specifying the relative weighting of features and labels. The variables γ and λ parameterize the selection of this worst-case distribution in response to the model weights β. Finally, Ψ is the logistic loss kernel t 7→ log(et + e−t) and LΨ = 1 is the corresponding Lipschitz constant. In all the experiments, we set δ = κ = 1 and c = 10−3.
We now show how we converted this problem to the form (1) for our experiments. Let z be a shorthand for (λ, β, γ), and define
L(z) .= λ(δ − κ) + 1 m m∑ i=1 Ψ(〈x̂i, β〉) + 1 m m∑ i=1 γi(ŷi〈x̂i, β〉 − λκ).
The first-order necessary and sufficient conditions for the convex-concave saddlepoint problem in (76) are
0 ∈ B(z) +A1(z) +A2(z) (77)
where the vector field B(z) is defined as
B(z) . = [ ∇λ,βL(z) −∇γL(z) ] , (78)
with
∇λ,βL(z) =
[ δ − κ(1 + 1m ∑m i=1 γi)
1 m ∑m i=1 Ψ ′(〈x̂i, β〉)x̂i + 1m ∑m i=1 γiŷix̂i ] and
∇γL(z) =  1 m (ŷ1〈x̂1, β〉 − λκ)
... 1 m (ŷm〈x̂m, β〉 − λκ)  . It is readily confirmed that B defined in this manner is Lipschitz. The monotonicity of B follows from its being the generalized gradient of a convex-concave saddle function (Rockafellar, 1970).
For the set-valued operators, A1(z) corresponds to the constraints and A2(z) to the nonsmooth `1 regularizer, and are defined as
A1(z) . = NC1(λ, β)×NC2(γ),
where
C1 . = { (λ, β) : ‖β‖2 ≤ λ/(LΨ + 1) } and C2 . = {γ : ‖γ‖∞ ≤ 1},
and
A2(z) . = {01×1} × c∂‖β‖1 × {0m×1}.
Here, the notation 0p×1 denotes the p-dimensional vector of all zeros. C1 is a scaled version of the second-order cone, well known to be a closed convex set, while C2 is the unit ball of the `∞ norm, also closed and convex. Since A1 is a normal cone map of a closed convex set and A2 is the subgradient map of a closed proper convex function (the scaled 1-norm), both of these operators are maximal monotone and problem (77) is a special case of (1) for n = 2.
Stochastic oracle implementation The operator B : Rm+d+1 7→ Rm+d+1, defined in (78), can be written as
B(z) = 1
m m∑ i=1 Bi(z)
where
Bi(z) . =  δ − κ(1 + γi) Ψ′(〈x̂i, β〉)x̂i + γiŷix̂i 0(i−1)×1
−(ŷi〈x̂i, β〉 − λκ) 0(m−i)×1
 .
In our SPS experiments, the stochastic oracle for B is simply B̃(z) = 1|B| ∑ i∈BBi(z) for some minibatch B ⊆ {1, . . . ,m}. We used a batchsize of 100.
Resolvent computations The resolvent of A1 is readily constructed from the projection maps of the simple sets C1 and C2, while the resolvent A2 involves the proximal operator of the `1 norm. Specifically,
JρA1(z) =
[ projC1(λ, β)
projC2(γ)
] and JρA2(z) =  01×1proxρc‖·‖1(β) 0m×1  . The constraint C1 is a scaled second-order cone and C2 is the `∞ ball, both of which have closed-form projections. The proximal operator of the `1 norm is the well-known soft-thresholding operator (Parikh & Boyd, 2013, Section 6.5.2). Therefore all resolvents in the formulation may be computed quickly and accurately.
SPS stepsize choices For the stepsize in SPS, we ordinarily require ρk ≤ ρ < 1/L for the global Lipschitz constant L of B. However, since the global Lipschitz constant may be pessimistic, better performance can often be achieved by experimenting with larger stepsizes. If divergence is observed, then the stepsize can be decreased. This type of strategy is common for SGD and similar stochastic methods. Thus, for SPS-decay we set αk = Cdk−0.51 and ρk = Cdk−0.25, and performed a grid search to select the best Cd from {0.1, 0.5, 1, 5, 10}, arriving at Cd = 1 for epsilon and SUSY, and Cd = 0.5 for real-sim. For SPS-fixed we used ρ = K−1/4 and α = Cfρ2, and performed a grid search to select Cf over {0.1, 0.5, 1, 5, 10}, arriving at Cf = 1 for epsilon and real-sim, and Cf = 5 for SUSY. The total number of iterations for SPS-fixed was chosen as follows: For the epsilon dataset, we used K = 5000, for SUSY we used K = 200, and for real-sim we used K = 1000.
Parameter choices for the other algorithms All methods are initialized at the same random point. For Tseng’s method, we used the backtracking linesearch variant with an initial stepsize of 1, θ = 0.8, and a stepsize reduction factor of 0.7. For FRB, we used the backtracking linesearch variant with the same settings as for Tseng’s method. For deterministic PS, we used a fixed stepsize of 0.9/L. For the stochastic Tseng’s method of Böhm et al. (2020), the stepsize αk must satisfy: ∑∞ k=1 αk =∞
and ∑∞ k=1 α 2 k < ∞. So we set αk = Ck−d and perform a grid search over {C, d} in the range [10−4, 10]× [0.51, 1], checking 5× 5 values to find the best setting for each of the three problems. The selected values are in Table 1.
The work of Böhm et al. (2020) also introduced FBFp, a stochastic version of Tseng’s method that reuses a previously-computed gradient and therefore only needs one additional gradient calculation per iteration. In our experiments, the performance of the two methods was about the same, so we only report the performance of stoch. Tseng’s method.
For variance-reduced FRB, the main parameter is the probability p. We hand-tuned p,arriving at p = 0.01 for all problems. We set the stepsize to its maximum allowed value of
τ = 1− √
1− p 2L .
Plots versus Epoch Figure 2 plots the performance of each method versus epoch (i.e. data pass). This shows an even more dramatic benefit for the stochastic methods than the plots versus time, since at each iteration the stochastic methods only need to process small amounts of data, whereas deterministic methods must process all of it. We believe these benefits do not fully manifest themselves in the plots versus time due to overheads in each iteration of the stochastic methods, multithreading providing a boost for the deterministic methods, memory access patterns, and other practical considerations.
Fraction of Nonzero Entries versus Running time Figure 3 plots the fraction of nonzero entries in the iterates of each method versus running time. For each method, we used output of proxc‖·‖1 . We observe that our methods produce sparse intermediate iterates for two of the three problems. This is one of the benefits of proximal splitting algorithms in general, including our method. For the other problem, SUSY, no method produces sparse iterates, suggesting that c should be increased if sparse solutions are desired.

heading: J LOCAL CONVERGENCE ON NON-MONOTONE PROBLEMS
text: The work by Hsieh et al. (2020) provides a local convergence analysis for DSEG applied to locally monotone problems. Recall that DSEG is equivalent to the special case of SPS for which n = 0. While extending this result to the more general setting of SPS is beyond the scope of this manuscript, we next provide a preliminary sketch of how the analysis of Hsieh et al. (2020) might be generalized to our setting. We leave a formal proof to future work.
Sketch of assumptions and main result The first assumption needed is the existence of an isolated solution p∗ = (z∗, w∗1 , . . . , w ∗ n+1) ∈ S. We then require that there exists a ball Br(z∗), centered at z∗, throughout which the operator B is “well-behaved”, meaning that it satisfies monotonicity and Lipschitz continuity. In addition, we need each Ai, for i ∈ 1..n, to be maximal monotone within this ball. Outside of the ball, the operators do not need to be monotone or Lipschitz.
Following (Hsieh et al., 2020, Assumption 2′), the noise variance assumptions are slightly stronger than in the monotone case. In particular, we require that E[‖ k‖q|Fk] ≤ Nq and E[‖ek‖q|Fk] ≤ Nq for some q > 2. As before, the noise must be zero-mean. Finally, the stepsize requirements are also slightly stronger than (12), having the added assumption that ∑∞ k=1 ρ q k <∞.
With these assumptions, the goal is to show that, so long as the initial point p1 is sufficiently close to p∗, then with high probability pk converges to p∗.
Proof strategy The initial strategy is to develop the following recursion, satisfied by SPS, that does not (yet) utilize local monotonicity or Lipschitz continuity:
‖pk+1 − p∗‖2 ≤ (1 + c1α2k)‖pk − p∗‖2 − c2αkρk(T ′k + lk + rk)− c3αk(r′k + qk) + c1α 2 k ( ‖ek‖2 + ‖ k‖2 + c4 ) + c5αkq ′ k (79)
for appropriate constants c1 . . . c5 ≥ 0. In this inequality, we use
T ′k . = τ
ρ n∑ i=1 ‖yki − wki ‖2 + 1 ρτ n∑ i=1 ‖zk − xki ‖2,
lk . = n∑ i=1 〈z∗ − xki , w∗i − yki 〉+ 〈 z∗ − xkn+1, w∗i −B(xkn+1) 〉 , rk . = 〈 k, B(x̃k)− wkn+1〉, r′k . = 〈zk − z∗, ek〉, qk . = (ρ−1k − d/2)‖x̃ k − zk‖2 − ‖x̃k − zk‖‖B(x̃k)−B(zk)‖
q′k . = ρk‖ k‖‖Bxkn+1 −Bx̃k‖+
1
2d ‖Bx̃kn+1 −Bxkn+1‖2,
where
x̃k . = zk − ρk ( B(zk)− wkn+1 ) d . =
1− ρL 1 + ρ/2 , (80)
with L being the local Lipschitz constant of B on Br(z∗). The iterate x̃k is the analog of the iterate X̃t+1/2 used in Hsieh et al. (2020).
The recursion (79) is derived by once again starting from (13) and following the arguments leading to (35), but this time not taking conditional expectations. In particular, the upper bounds on ‖∇zϕk‖2 and ‖∇wiϕk‖2 contribute the terms c1α2k ( ‖ek‖2 + ‖ k‖2 + c4) and c1α2k‖pk − p∗‖2. For i ∈ 1..n, the “ϕi,k-gap" term, ϕi,k(pk)− ϕi,k(p∗), is dealt with in a similar manner to Section C.5, but this time not using monotonicity as in (36). This contributes T ′k and the first term in lk. Finally, as we sketch below, the “ϕn+1,k-gap" term contributes rk, r′k, qk, q ′ k, and the last term in lk.
For the “ϕn+1,k-gap”, that is, ϕn+1,k(pk) − ϕn+1,k(p∗), we have to depart from the analysis in Section C.6 and use an alternative argument involving x̃k. We now provide some details of this argument: in the following, we use Bz as shorthand for B(z) for any vector z ∈ Rd. We begin the analysis with
ϕn+1,k(p k) = 〈zk − xkn+1, ykn+1 − wkn+1〉
= 〈zk − xkn+1, Bxkn+1 − wkn+1〉+ 〈zk − xkn+1, ek〉︸ ︷︷ ︸ part of r′k . (81)
The final term will combine with the term 〈xkn+1 − z∗, ek〉 coming from
−ϕn+1,k(p∗) = 〈z∗ − xkn+1, w∗n+1 − ykn+1〉 = 〈z∗ − xkn+1, w∗n+1 −Bxkn+1〉+ 〈xkn+1 − z∗, ekn+1〉 (82)
to yield r′k above. Equation (82) also yields the second term in lk. Using that x̃ k − xkn+1 = ρk k, we rewrite the first term in (81) as
〈zk − xkn+1, Bxkn+1 − wkn+1〉 = 〈zk − x̃k, Bxkn+1 − wkn+1〉+ 〈x̃k − xkn+1, Bxkn+1 − wkn+1〉 = 〈zk − x̃k, Bxkn+1 − wkn+1〉+ ρk〈 k, Bxkn+1 − wkn+1〉 = 〈zk − x̃k, Bxkn+1 − wkn+1〉+ ρk〈 k, Bxkn+1 −Bx̃k〉 (83)
+ ρk 〈 k, Bx̃k − wkn+1〉︸ ︷︷ ︸ rk .
Next, the terms in (83) admit the lower bound
〈zk − x̃k, Bxkn+1 − wkn+1〉+ ρk〈 k, Bxkn+1 −Bx̃k〉 ≥ 〈zk − x̃k, Bxkn+1 − wkn+1〉 − ρk‖ k‖‖Bxkn+1 −Bx̃k‖︸ ︷︷ ︸
first part of q′k
.
Considering the first term on right-hand side of this bound, we also have
〈zk − x̃k, Bxkn+1 − wkn+1〉 = 〈zk − x̃k, Bx̃k − wkn+1〉+ 〈zk − x̃k, Bxkn+1 −Bx̃k〉
≥ 〈zk − x̃k, Bx̃k − wkn+1〉 − d 2 ‖zk − x̃k‖2 − 1
2d ‖Bx̃k −Bxkn+1‖2︸ ︷︷ ︸
second part of q′k
for any d > 0, using Young’s inequality. Finally, for the first two terms of the right-hand side of the above relation, we may write
〈zk − x̃k, Bx̃k − wkn+1〉 − d
2 ‖zk − x̃k‖2
= 〈zk − x̃k, Bzk − wkn+1〉+ 〈zk − x̃k, Bx̃k −Bzk〉 − d
2 ‖zk − x̃k‖2
≥ (ρ−1k − d/2)‖z k − x̃k‖2 − ‖zk − x̃k‖‖Bx̃k −Bzk‖︸ ︷︷ ︸
qk
,
where in the final inequality we use the Cauchy-Schwartz inequality and substitute Bzk − wkn+1 = ρ−1k (z
k − x̃k), from the definition of x̃k in (80). We have now accounted for all the terms appearing in (79).
The recursion (79) is analogous to equation (F.7) on page 24 of Hsieh et al. (2020) and provides the starting point for the local convergence analysis. The next step would be to derive an analog of Theorem F.1. of Hsieh et al. (2020) using (79). The following translation to the notation of Theorem F.1. could be used (note that Hsieh et al. (2020) uses t for iteration counter):
Dk = ‖pk − p∗‖2, ζk = c2αkρk(T ′ k + lk) + c3αkqk,
ξk = −c2αkρkrk − c3αkr′k, χk = c1α 2 k ( ‖ek‖2 + ‖ k‖2 + ‖pk − p∗‖2 + c4 ) + c5αkq ′ k,
and the event Eρ∞ is translated to Eρ∞ = { xkn+1 ∈ Br(z∗), x̃k ∈ Bρr(z∗), pk ∈ Bρr(p∗) for all k = 1, 2, . . . } .
An analog of Theorem 2 of Hsieh et al. (2020) could then be developed based on this result.

Review summaries:
The paper focuses on the stochastic variant of the projective splitting (PS) algorithm. With a specific focus on monotone inclusion problems, the authors propose a novel separable algorithm featured by the ability to handle multiple constraints and non-smooth regularizers. Compared with similar approaches on variational inequality which is a special case of monotone inclusions, this paper uses a more direct error metric than the restricted gap function. Moreover, although with a slower convergence rate, this paper is the first discussion under the general discontinuous monotone inclusion case.