Year: 2021

Abstract:
The large-scale linearly constrained nonsmooth nonconvex optimization finds wide applications in machine learning, including non-PSD Kernel SVM, linearly constrained Lasso with nonsmooth nonconvex penalty, etc. To tackle this class of optimization problems, we propose an efficient algorithm called Nonconvex Randomized Primal-Dual Coordinate (N-RPDC) method. At each iteration, this method only randomly selects a block of primal variables to update rather than updating all the variables, which is suitable for large-scale problems. We provide two types of convergence results for N-RPDC. We first show that any cluster point of the sequence of iterates generated by N-RPDC is almost surely (i.e., with probability 1) a stationary point. In addition, we also provide an almost sure asymptotic convergence rate of O(1/ √ k). Next, we establish the expected O(ε−2) iteration complexity of N-RPDC in order to drive a natural stationarity measure below ε in expectation. The fundamental aspect to establishing the aforementioned convergence results is a surrogate stationarity measure we discovered for analyzing N-RPDC. Finally, we conduct a set of experiments to show the efficacy of NRPDC.

text: The large-scale linearly constrained nonsmooth nonconvex optimization finds wide applications in machine learning, including non-PSD Kernel SVM, linearly constrained Lasso with nonsmooth nonconvex penalty, etc. To tackle this class of optimization problems, we propose an efficient algorithm called Nonconvex Randomized Primal-Dual Coordinate (N-RPDC) method. At each iteration, this method only randomly selects a block of primal variables to update rather than updating all the variables, which is suitable for large-scale problems. We provide two types of convergence results for N-RPDC. We first show that any cluster point of the sequence of iterates generated by N-RPDC is almost surely (i.e., with probability 1) a stationary point. In addition, we also provide an almost sure asymptotic convergence rate of O(1/ √ k). Next, we establish the expected O(ε−2) iteration complexity of N-RPDC in order to drive a natural stationarity measure below ε in expectation. The fundamental aspect to establishing the aforementioned convergence results is a surrogate stationarity measure we discovered for analyzing N-RPDC. Finally, we conduct a set of experiments to show the efficacy of NRPDC.

heading: 1 INTRODUCTION
text: Many large scale problems arising in machine learning amounts to solving the following linearly constrained nonsmooth nonconvex optimization problem:
min x∈X F (x) = f(x) + g(x)
s.t Ax− b = 0, (P)
where A ∈ Rn×d and b ∈ Rn. Throughout this paper, we impose the following assumptions on problem (P): the function f : Rd → R is possibly nonconvex and continuously differentiable with its gradient ∇f being Lf -Lipschitz continuous, X is a convex and compact set of Rd, i.e., there exits a positive number M such that M = maxx,x′∈X ‖x− x′‖, and g : Rd → R is nonsmooth and nonconvex. More precisely, g is assumed to be lower semicontinuous (l.s.c) ρg-weakly convex and has bounded subgradients over X, i.e., there exists Lg > 0 such that ‖s‖ ≤ Lg,∀s ∈ ∂g(x), x ∈ X, where ∂g is the subdifferential of g (see (4) for definition). Let F ∗ be the optimal value of (P).
Recall that a function g is said to be ρg-weakly convex if g(·) + ρg2 ‖ · ‖ 2 is convex for some constant ρg ≥ 0 (Vial, 1983). It is worth mentioning that a wide class of nonsmooth nonconvex functions belong to the weakly convex class; see, e.g., (Vial, 1983; Davis & Drusvyatskiy, 2019) for more discussions on weak convexity.
We further assume g(x) = ∑N i=1 gi(xi) is block separable with respect to the space decomposition
X = X1 × X2 × · · · × XN , xi ∈ Xi ⊂ Rdi , and N∑ i=1 di = d. (1)
Let A = (A1, A2, . . . , AN ) ∈ Rn×d be the corresponding partition of A, where Ai ∈ Rn×di . When N = d, each Xi ⊂ R corresponds to a box constraint on the coordinate xi for i = 1, . . . , d.
Motivations. Our motivation for studying large scale problem (P) stems from the fact that this problem class finds wide applications in machine learning. To be more specific, we will present several practical examples below that give rise to (P).
Application 1: Non-PSD kernel support vector machine. Support vector machine (SVM) is a widely utilized technique for supervised learning (Boser et al., 1992; Cortes & Vapnik, 1995). In order to improve the interpretability and enhance the robustness, some non-PSD kernels are used in the kernalized SVM method, such as the sigmoid kernel (Lin & Lin, 2003), the jittering kernel (DeCoste & Schölkopf, 2002), the tangent distance kernel (Haasdonk & Keysers, 2002, August), to name a few. Formally, the kernalized SVM problem can be written as
min x∈[0,c]d
1 2x >Qx− 1>d x
s.t. y>x = 0, (2)
where Q is a d× d non-PSD matrix, c ∈ R is the upper bound of all variables, y ∈ {−1, 1}d is the vector of labels, and 1d is an d-dimensional vector of all 1s.
Application 2: Linearly constrained Lasso with nonsmooth nonconvex penalty. The Lasso is one of the most popular methods for variable selection. The standard Lasso uses `1-norm to select important variables. Some works propose to utilize nonsmooth nonconvex regularizers to further improve the performance of the Lasso; see, e.g., (Breheny & Huang, 2011; 2015; Rakotomamonjy et al., 2019). If additional prior information is available, the works (Deng et al., 2020, May; Gaines et al., 2018; James et al., 2013; Won et al., 2019) propose the constrained Lasso model. The linearly constrained LASSO with nonsmooth nonconvex penalty can be formulated as follows:
min x∈X
1 2‖Ax− b‖ 2 + d∑ i=1 φ(xi)
s.t. Bx− c = 0 (3)
where x = (x1, ..., xd) >, X ⊂ Rd, A ∈ Rn×d is the design matrix, b ∈ Rn is the response vector, B ∈ Rm×d and c ∈ Rm are given constraints, and φ is a nonsmooth nonconvex regularizer such as MCP (Zhang et al., 2010) or SCAD (Fan & Li, 2001) penalty.
There are plenty of other applications that give rise to (P), e.g., the robust M-estimators, distributed learning, etc. Due to the limitation of space, we will not expand them in details here.
Related works. Linearly constrained nonconvex optimization. Perhaps the most widely utilized class of algorithms for solving constrained optimization problems are the Primal-dual methods. Previous works on single-loop primal-dual methods mainly consider linearly contained convex optimization problems. The work (Boţ & Nguyen, 2020) considers a nonconvex instance with a special linear constraint Ax− z = 0, where both x and z are decision variables and the matrix A is assumed to be surjective. Under these assumptions, convergence results were established for a proximal alternating direction method of multipliers (ADMM). The recent works (Zhang & Luo, 2020a;b) provide convergence results of ADMM using gradient-based updates for linearly constrained smooth nonconvex optimization problems with a general linear constraint Ax = b, i.e., problem (P) with g ≡ 0. The main insights of their results are a construction of an proximally regularized auxiliary problem that dates back to (Bertsekas, 1979) and a properly constructed Lyapunov function. Then, they established descent property on the Lyapunov function based on the dual error bound condition (see also (Hong & Luo, 2017)), which leads to iteration complexity results. The method introduced in (Zhu et al., 2020) can be used to solve problem (P) with specific linear constraints. This method is based on the auxiliary problem principle of augmented Lagrangian (APP-AL) method (Cohen & Zhu, 1984; Zhao & Zhu, 2019; 2021), which can be seen as a forward-backward splitting method and applies a Jacobian updating strategy.1 However, both ADMM-type and APP-AL-type methods compute the full gradient of primal variables and update all primal variables at each iteration. Therefore, the computation complexity of one iteration of these two types of methods are expensive for large-scale problems (See Remark 2.1).
Randomized coordinate methods for unconstrained optimization. In the past decade, big data applications are ubiquitous in machine learning. The formulated optimization problems often involve very large datasets, and hence computing the function value or the gradient can be very expensive. These
1ADMM applies a Gauss-Seidel like minimization strategy, which is related to the Douglas-Rachford splitting method.
observations motivate the study of the randomized coordinate methods. The researches on this type of methods can be traced back to (Nesterov, 2012; 2014). In (Nesterov, 2012), Nesterov studied the iteration complexity of the randomized coordinate gradient (RCG) descent method for smooth convex optimization. Later, in (Nesterov, 2014), the same author analyzed the randomized coordinate subgradient method for a set of piece-wise linear nonsmooth convex optimization problems. The works (Richtárik & Takáč, 2014) and (Lu & Xiao, 2015) extends Nesterov’s results to convex additive composite optimization. Iteration complexity and linear convergence of RCG were studied in (Patrascu & Necoara, 2015) and (Karimi et al., 2016) for nonsmooth nonconvex optimization and smooth nonconvex optimization, respectively. Moreover, by employing an asynchronous point to evaluate the gradient in each iteration, The works (Liu & Wright, 2015) and (Liu et al., 2014) establish the iteration complexity results of asynchronous RCG for convex smooth optimization and additive composite optimization, respectively.
A Few results on randomized coordinate methods for linearly constrained optimization. Though the randomized coordinate-type methods for unconstrained optimization are extensively studied, there are only a few results for this type of methods concerning coupled linearly constrained optimization problems. The works (Gao et al., 2019; Wang et al., 2014) study the randomized primal-dual coordinate (RPDC) method and provides iteration complexity results for linearly constrained convex optimization. Better complexity results of RPDC were obtained for linearly constrained strongly convex optimization (Xu & Zhang, 2018). The recent works (Zhu & Zhao, 2020; Latafat et al., 2019; Fercoq & Bianchi, 2019; Alacaoglu et al., 2020) provide almost sure (i.e., with probability 1) asymptotic convergence results of RPDC-type methods for linearly constrained convex optimization. It is worth emphasizing that the works (Latafat et al., 2019; Fercoq & Bianchi, 2019; Alacaoglu et al., 2020) establish sequential convergence results (i.e., the convergence of the whole sequence of iterates).
However, the above mentioned works only considers either ADMM-type and APP-AL-type methods for linearly constrained smooth nonconvex optimization problem or the randomized coordinate methods for linearly constrained convex optimization. To the best of our knowledge, no previous result concerns convergence and iterate complexity of the randomized coordinate methods for linearly constrained nonsmooth nonconvex optimization, which will be the main focus of this paper.
Main contributions. In this paper, we aim to solve the large-scale linearly constrained nonsmooth nonconvex optimization problem (P). To tackle it, we propose the Nonconvex Randomized PrimalDual Coordinate (N-RPDC) method based on an auxiliary problem (see Section 2). Relying on the mild local uniform metric subregularity property (see Section 3.3), we provide two types of convergence results for N-RPDC (see Section 4). We first show that any cluster point of the sequence of iterates generated by N-RPDC is almost surely (i.e., with probability 1) a stationary point of (P). In addition, we also provide an almost sure asymptotic convergence rate of O(1/ √ k), where k represents iteration number. Next, we establish the expected iteration complexity of N-RPDC. Namely, N-RPDC needs at most O(ε−2) number of iterations in order to drive a surrogate stationarity measure below ε in expectation.
The fundamental aspect to our convergence analysis is a surrogate stationarity measure we discovered for analyzing N-RPDC (see Section 3.1), which is one of the main contributions of this work. The standard stationarity measure (i.e., checking KKT conditions) is not directly applicable here due to the stochastic nature of N-RPDC. Instead, we define the notion of reference point and use its distance to the iterate generated by N-RPDC as a surrogate stationarity measure. Then, such a surrogate stationarity measure is clarified by showing that it is an upper bound of the standard stationarity measure at the reference point up to a numerical constant (see Proposition 3.1 and Remark 3.1).
There are also some other interesting techniques that we utilized for establishing the above convergence results. For example, the utilization of the local uniform metric subregularity is crucial in order to deal with the nonsmooth nonconvex term g in problem (P), as it is far from obvious how to prove the previously used error bound condition if g is not null (see Section 3.3).
Notations. We use 〈·, ·〉 and ‖ · ‖ to denote the Euclidean inner product and the Euclidean norm, respectively. For a matrix A, its minimum eigenvalue is denoted by λmin(A). The spectral norm of a matrix A is denoted by ‖A‖. Let S be a subset of Rd. We use projS(z) to denote the orthogonal projector onto S and dist(z,S) := infx∈S ‖x − z‖ to denote the distance between x and S. When
S = ∅, we set dist(z,S) = +∞. IS(z) = {
0, z ∈ S +∞, z /∈ S represents the indicator function of the
set S. If S is a convex set, then the limiting normal cone to S is defined asNS(z) = {ξ : 〈ξ, ζ−z〉 ≤ 0,∀ζ ∈ S}.

heading: 2 AUXILIARY PROBLEM AND PROPOSED ALGORITHM
text: In this section, we will define an auxiliary problem of (P). An important feature of the auxiliary problem is that it has exactly the same set of KKT points to the original problem (P). Furthermore, the auxiliary problem has a better geometric structure, which makes it easier to solve. Finally, we will propose a primal-dual coordinate method based on the auxiliary problem.
Let us first present some preliminaries on the subdifferential and KKT conditions. Due to the fact that problem (P) is highly nonsmooth and nonconvex, our general purpose is to design an efficient algorithm for finding a stationary point rather than a globally optimal solution. Hence, we will define certain suitable stationarity measures in this subsection.
Recall that the function g in problem (P) is τ -weakly convex. By (Vial, 1983, Proposition 4.6), we have
∂g(x) = ∂h(x)− τx, (4) where h is the associated convex function such that g(x) = h(x)− τ2‖x||
2 (see (Vial, 1983, Proposition 4.3) for the guarantee of the existence of such a function h) and ∂h(x) is the usual convex subdifferential. Thus, the subdifferential of a weakly convex function is always well defined.
The Lagrangian of problem (P) can be written as L(x, p) = f(x) + g(x) + 〈p,Ax− b〉,
where p is the dual variable. Then, we have the following KKT conditions for problem (P): 0 ∈ ∂L(x, p) := ( ∇f(x) + ∂g(x) +NX(x) +A>p
Ax− b
) . (5)
A feasible point (x, p) ∈ X × Rn is called a stationary point of (P) if it satisfies the above KKT conditions, i.e., it satisfies dist (0, ∂L(x, p)) = 0.

heading: 2.1 THE AUXILIARY PROBLEM
text: Thanks to the construction used in (Bertsekas, 1979), we introduce an auxiliary problem to (P) in the following:
min x∈X,z∈Rd
F+(x, z) = f(x) + g(x) + σ2 ‖x− z‖ 2
s.t Ax− b = 0, (P+)
where σ > Lf + ρg is a regularization parameter.
We now characterize the set of stationary points of (P+). Towards that end, let w = (x, z, p) ∈ X× Rd × Rn. The Lagrangian of (P+) is given by
L+(w) = L+(x, z, p) = f(x) + g(x) + σ 2 ‖x− z‖2 + 〈p,Ax− b〉.
Thus, we have the following KKT conditions for problem (P+):
0 ∈ ∂L+(w) = ∂L+(x, z, p) :=  ∇f(x) + ∂g(x) + σ(x− z) +NX(x) +A>pσ(z − x) Ax− b  . (6) A point w = (x, z, p) ∈ X × Rd × Rn satisfying the above KKT conditions is called a stationary point of (P+). One crucial feature of the auxiliary problem is that it has exactly the same set of stationary points as that of the original problem (P), which is presented in the following lemma.
Lemma 2.1 (relation between stationary points of (P) and (P+)). The following two statements are equivalent: (a) (x∗, p∗) is a stationary point of problem (P); and (b) (x∗, x∗, p∗) is a stationary point of problem (P+).
Algorithm 1 N-RPDC: Nonconvex Randomized Primal-Dual Coordinate Method Initialization: set x0 ∈ X, z0 ∈ Rd, and p0 ∈ Null(A>); Step sizes η, αx, and αz .
1: for k = 0, 1, . . . do 2: Update p as pk+1 = pk + η(Axk − b). 3: Choose i(k) from {1, ..., N} uniformly at random; 4: Update x and z through coordinate steps:
xk+1 = arg min x∈X
〈 ∇i(k)f(xk) + σ(xk − zk)i(k), xi(k) 〉 + gi(k)(xi(k)) + 〈 pk+1 + γ(Axk − b), Ai(k)xi(k) 〉 + D(x, xk)
αx ,
zk+1 = arg min z∈Rd
〈 σ(zk − xk)i(k), zi(k) 〉 + 1
2αz ‖z − zk‖2.
5: end for
Proof. We first prove that part (a) implies part (b). Since (x∗, p∗) is a stationary point of (P), we have 0 ∈ ∂L(x∗, p∗). Upon plugging z = x∗ in the L+(x∗, z, p∗) shows that 0 ∈ ∂L+(x∗, x∗, p∗). The inverse direction is proved by a similar argument.
The auxiliary problem can be viewed as a quadratically regularized version of the original problem (P) with an additional variable z. Recall that f has Lf -Lipschitz gradient and g is ρg-weakly convex. We have the following benign properties of the auxiliary problem (P+), whose derivations can be found in Appendix A.2:
(1) f(x) + σ2 ‖x− z‖ 2 is continuously differentiable in (x, z) with Lipschitz continuous gradient;
(2) F+(x, z) = f(x) + g(x) + σ2 ‖x− z‖ 2 is weakly convex in (x, z) with parameter Lf + ρg;
(3) F+(x, z) is bi-strongly convex, i.e., F+(x, z) is strongly convex in x for every fixed z and F+(x, z) is strongly convex in z for every fixed x.

heading: 2.2 NONCONVEX RANDOMIZED PRIMAL-DUAL COORDINATE METHOD
text: As we established in the last subsection, problem (P+) has the same set of stationary points as that of the original problem (P) and the former has a series of benign properties. These observations motivate us to design algorithm for solving the original problem by targeting on (P+). Similar to standard primal-dual methods, our algorithm for solving problem (P+) builds on the following augmented Lagrangian of (P+):
L+γ (x, z, p) = f(x) + g(x) + σ 2 ‖x− z‖2 + 〈p,Ax− b〉+ γ 2 ‖Ax− b‖2. (7)
With the augmented Lagrangian, our algorithm is designed by performing an APP-AL update over (x, z, p) at each iteration. Furthermore, the updates for variables x and z are achieved through randomized coordinate steps. We depict our algorithmic procedures—i.e., the N-RPDC method— in Algorithm 1.
The proximal term D(x, xk) used in N-RPDC is the so-called Bregman distance function. Let βK be the strong convexity parameter of the Bregman distance function D(x, xk). It is safe for the readers to think D(x, xk) = 12‖x− x
k‖2 for now. We put the more general choices of D(x, xk) in Appendix A.3.
It is worth mentioning that the subproblem for updating xk+1 in N-RPDC has a closed form solution once the Bregman distance function D(x, xk) = 12‖x− x
k‖2 and gi(xi) is SCAD, MCP, quadratic, or `ν-norms with ν ∈ {1, 2,∞}. Remark 2.1 (computational complexity). It is claimed in (Nesterov, 2012; Wright, 2015) that a coordinate method can have a much lower computational complexity than its full counterpart. Let us take the second motivating application (i.e., nonsmooth nonconvex constrained Lasso, problem (3) in Application 2) in Section 1 as an illustrative example to show that the computational complexity of one iteration of N-RPDC with N = d can be much cheaper than the full primal-dual counterpart (ADMM-type or APP-type method). Note that in this case, each Xi is a box constraint on the coordi-
nate xi. Suppose we set D(x, xk) = 12‖x− x k‖2. Then, the update of wk+1 = (xk+1, zk+1, pk+1) consists of: 1) pk+1 = pk + ηsk. 2) xk+1i(k) = proxαx,IXi(k)+φ ( xki(k) − αx( ( Ai(k) )> rk + σ(xk −
zk)i(k) + ( Bi(k) )> (pk+1 + γsk)) ) and xk+1j 6=i(k) = x k j 6=i(k), where IXi(k) is the indicator function of Xi(k). 3) zk+1i(k) = z k i(k) − αzσ(z k − xk)i(k) with zk+1j 6=i(k) = z k j 6=i(k). We also update the residual vectors as rk+1 = rk+Ai(k)(x k+1 i(k)−x k i(k)) and s k+1 = sk+Bi(k)(x k+1 i(k)−x k i(k)). If we precompute r0 = Ax0 − b and s0 = Bx0 − c, and the proximal mapping proxαx,IXi(k)+φ admits a closed form solution (which is the case when φ is the SCAD or MCP penalty), then the above updates have a total complexity O(max{m,n}). Therefore, the complexity of N-RPDC is d× cheaper than its full counterpart (which has a computational complexity of O(dmax{m,n}) in one iteration).

heading: 3 STATIONARITY MEASURE, LYAPUNOV FUNCTION, AND UNIFORM METRIC SUBREGULARITY
text: In this section, we will establish a series of important results, which serve as the foundations for our later convergence analysis.

heading: 3.1 THE STATIONARITY MEASURE
text: Consider the current iterate wk. If wk+1 is updated by the full primal-dual counterpart to N-RPDC (i.e., N = 1 in Algorithm 1), then the step length ‖wk − wk+1‖ is a natural stationarity measure due to the facts that dist(0, ∂L+(wk+1)) = O(‖wk − wk+1‖) with some step sizes hidden in the big-O and limk→+∞ ‖wk − wk+1‖ = 0 can be established by an easy argument based on an appropriate Lyapunov function. However, if wk+1 is updated by our N-RPDC with N 6= 1, which is a stochastic algorithm that only updates part of wk to wk+1, it is no longer obvious to show limk→+∞ ‖wk − wk+1‖ = 0. Therefore, due to the stochastic nature of N-RPDC, the step length ‖wk − wk+1‖ cannot play the role of stationarity measure. This observation motivates us to work with a surrogate stationarity measure. Towards that end, for a given point w = (x, z, p), we define the reference point of w as T (w) := (Tx(w), Tz(w), Tp(w)) and Tp(w) = p+ η(Ax− b); Tx(w) = arg min y∈X 〈∇f(x) + σ(x− z), y〉+ g(y) + 〈Tp(w) + γ(Ax− b), Ay〉+ D(y, x) αx ;
Tz(w) = arg min y∈Rd 〈σ(z − x), y〉+ 1 2αz ‖y − z‖2.
(8)
One possible way to understanding the reference point is to let wk = (xk, zk, pk), then T (wk) is nothing else other than the next iterate wk+1 if the full primal-dual method is utilized, i.e., N = 1 in N-RPDC. However, it is worthing mentioning in the case where N > 1, we do not need to compute the reference point T (wk) explicitly. Instead, N-RPDC only use a block coordinate for updating.
The next proposition provides us a valid stationary measure, which is crucial for our later convergence analysis.
Proposition 3.1 (surrogate stationarity measure). For all w = (x, z, p) ∈ X× Rd × Rn, we have dist ( 0, ∂L+(T (w)) ) ≤ c‖w − T (w)‖ =: Φ(w).
Here, T (w) is the reference point of w defined in (8) and c := √ 3 max { LK αx +Lf +2σ+‖A‖, 1αz + 2σ, γ‖A‖+1η } is a positive constant.
Remark 3.1 (interpretation of the surrogate stationarity measure). The mapping Φ(w) will serve as the surrogate stationarity measure for analyzing the convergence of N-RPDC. Particularly, it is clear that dist ( 0, ∂L+γ (T (w)) ) = 0 when Φ(w) = 0. In order to characterize the iteration complexity of N-RPDC, we focusing on achieving Φ(w) ≤ ε since it implies that w is ε-close to the reference point T (w) which is an ε-stationary point (here, T (w) is defined to be ε-stationary since
dist (0, ∂L+(T (w))) ≤ ε). These interpretations clarify our surrogate stationarity measure. Note that similar notion of surrogate stationarity measure dates back to Ekeland’s variational principle (Ekeland, 1974) and also appeared in recent advances of the analysis of the subgradient-type methods; see, e.g., (Davis & Drusvyatskiy, 2019; Li et al., 2021). It is interesting to note that the main difference between our approach and the existing ones lies in the choice of the reference point. Let us take the approach used in (Davis & Drusvyatskiy, 2019) for analyzing subgradient-type method as an example. The authors define the proximal mapping of the objective function as the reference point, while we utilize the update of the full primal-dual method (i.e., N = 1 in Algorithm 1) as the reference point. Such a main difference comes from different sources of difficulties. The main difficulty for the analysis of subgradient-type method in (Davis & Drusvyatskiy, 2019) comes from the fact that this algorithm is intrinsically not a descent method on the objective function, while our main difficulty is due to the stochastic nature of N-RPDC.

heading: 3.2 THE LYAPUNOV FUNCTION AND ONE-STEP ANALYSIS
text: Recall that the function F+ defined in (P+) is strongly convex in x for every fixed z (see our analysis in Section 2.1). Let us define
ν(z) = min x∈X, Ax−b=0 F+(x, z) and x(z) = arg min x∈X, Ax−b=0 F+(x, z). (9)
In addition, we also define ψγ(z, p) = min
x∈X L+γ (x, z, p) and x(z, p) = arg min x∈X L+γ (x, z, p), (10)
where L+γ is defined in (7). Next, mimicking the construction of the Lyapunov function used in (Zhang & Luo, 2020b), we define
Λ(w) = L+γ (w) + 2 (ν(z)− ψγ(z, p)) (11)
as the Lyapunov function for our convergence analysis. The first part of Λ(·) is augmented Lagrangian L+γ (w) of (P+). The second part is dual gap of (P+) with fixed z. Additionally, by the definition of ψγ(z, p), we have that L+γ (w) ≥ ψγ(z, p). By the strong duality, we have ν(z) = max
p ψγ(z, p). Therefore, we have the fact that Λ(w) ≥ ν(z) ≥ F ∗.
Standard analysis for primal-dual methods develops certain descent property on the Lyapunov function. Before that, we first present a preliminary one-step analysis for N-RPDC. Note that indices i(k), k = 0, 1, 2, . . . in N-RPDC are random variables. Thus, it generates a random output. We use
Fk := {i(0), i(1), . . . , i(k)} to denote the filtration generated by the random variables i(0), i(1), . . . , i(k). Clearly, we have Fk ⊂ Fk+1. We use Ei(k) and EFk to denote the expectations taken over the random variable i(k) and the filtration Fk, respectively. The following lemma provides estimations for ν(z) − ψγ(z, p) and L+γ (w) after one-step execution of N-RPDC.
Lemma 3.1 (one-step analysis of N-RPDC). Suppose 0 < αx ≤ βK/ [ Lf + 2σ + γ‖A‖2 + 5 ] and
0 < η < 1/ [ 2N‖A‖2 ( Lf+σ+γ‖A‖2+ LK αx
σ−Lf−ρg + 1
)2] . For all k ≥ 0, we have
Λ(wk)− Ei(k)Λ(wk+1)
≥ 2 N ‖xk − Tx(wk)‖2 + 1 N
[ 1
αz − 2σ − σ λ − 3σ
2
σ − Lf − ρg
] ‖zk − Tz(wk)‖2
+ η‖Ax(zk, Tp(wk))− b‖2 − σλ‖x(zk, Tp(wk))− x(zk)‖2.
(12)
where σ is defined in (P+) and λ is any positive constant.
Toward establishing certain decent property on the Lyapunov function in terms of the stationarity measure Φ(wk), we have to connect all the terms on the right-hand side of the inequality of Lemma 3.1 to Φ2(wk) = ‖wk−T (wk)‖2. Thus, it remains to upper bound ‖x(zk, Tp(wk))−x(zk)‖ by ‖Ax(zk, Tp(wk))− b‖, which needs the so-called uniform metric subregularity property.

heading: 3.3 UNIFORM METRIC SUBREGULARITY AND ITS CONSEQUENCE
text: Let us define the parameter set-valued mapping ∂Lzγ(x, p) = ( ∇f(x) + ∂g(x) + σ(x− z) +NX(x) +A>(p+ γ(Ax− b))
Ax− b
) (13)
as the subdifferential of L+γ with respect to (x, p) with fixed z. Clearly, we have ∂Lzγ(x(z, p), p) =( 0
Ax(z, p)− b
) due to the optimality of x(z, p) in (10). To proceed, we introduce the notion of
local uniform metric subregularity property for the parameter set-valued mapping. Definition 3.1 (local uniform metric subregularity, cf. (Kruger & Duy Cuong, 2021)). Let Hz(u) : U ⇒ V be a set-valued mapping. We call that Hz(u) satisfies the locally metric subregularity property with parameter (δ, κ) uniformly over all z ∈ Z if
dist ( u,H−1z (0) ) ≤ κ · dist (0,Hz(u)) , ∀z ∈ Z, (14)
for any u ∈ {u′ ∈ U : dist(0,Hz(u′)) ≤ δ}, where δ, κ > 0 are uniform constants.
This local error bound condition gives the following result. Lemma 3.2. Suppose the set-valued mapping ∂Lzγ(x, p) defined in (13) satisfies the local uniform metric subregularity property, then there is κ > 0 (independent of z) such that
‖x(z, p)− x(z)‖ ≤ κ‖Ax(z, p)− b‖, whenever ‖Ax(z, p)− b‖ ≤ δ.
Equipped with this lemma, we are able to upper bound ‖x(zk, Tp(wk)) − x(zk)‖ by ‖Ax(zk, Tp(wk)) − b‖ locally in order to establish the descent property on the Lyapunov function. We will discuss this descent property in the next section. Now, let us provide the conditions on problem (P) in order to ensure the local uniform metric subregularity assumed in Lemma 3.2. Lemma 3.3 (sufficient condition for local uniform metric subregularity). Consider problem (P). Suppose that∇f(x) is piecewise affine and ∂g(x) is polyhedral multifunction and X is a polyhedral set, then ∂Lzγ(x, p) defined in (13) satisfies the local uniform metric subregularity property. Remark 3.2 (generality of local uniform metric subregularity). This sufficient condition ensures that ∂Lzγ(x, p) must satisfy the local uniform metric subregularity property whenever f is quadratic (can be nonconvex, i.e., its Hessian is not required to be positive semidefinite), g can be written as ‖ · ‖1 − ρg2 ‖ · ‖
2, `2-norm, `∞-norm structure, and X is a polyhedron, which covers all the two motivating applications listed in Section 1. Remark 3.3 (comparison to the error bound condition utilized in (Zhang & Luo, 2020b)). Recall that the work (Zhang & Luo, 2020b) considers problem (P) with g ≡ 0. They use the so-called dual error bound condition (see (Hong & Luo, 2017)) instead of the local uniform metric subregularity property. However, we argue that it is far form obvious how to extend their technique to cover our nonsmooth nonconvex case, as it is not clear whether one can establish the dual error bound condition if g is not null. Remark 3.4 (further comments on the surrogate stationarity measure Φ(w) defined in Proposition 3.1). It is interesting that the sufficient condition in Lemma 3.3 also ensures that ∂L+(w) satisfies the local metric subregularity when Φ(w) ≤ ε with sufficient small ε (Robinson, 1981), i.e., dist(w,W) ≤ κ′ dist(0, ∂L+(w)) for some κ′ > 0, where W is the set of stationary points of (P+). Then, we have dist(w,W) ≤ (1/c + κ′)ε which follows from dist(T (w),W) ≤ κ′ε (due to Φ(w) ≤ ε) and the triangle inequality.

heading: 4 CONVERGENCE AND ITERATION COMPLEXITY OF N-RPDC
text: Equipped with all the machineries developed in Section 3, we are now ready to establish the convergence properties of N-RPDC.

heading: 4.1 DESCENT PROPERTY ON THE LYAPNOV FUNCTION
text: Upon invoking Lemma 3.2 in Lemma 3.1, we can derive the descent property in expectation on the Lyapnov Function.
Lemma 4.1 (expected sufficient decrease for Λk). Suppose that Assumptions of Lemma 3.1 hold. If the set-valued mapping ∂Lzγ(x, p) defined in (13) satisfies the local uniform metric subregularity property (see Definition 3.1). Let λ = min{ δη2σM2+1 , η σκ2+1}, the parameter αz satisfy 0 < αz <
1/ [ 2σ + σλ + 3σ2 σ−Lf−ρg + 1 ] , and c1 = min { 1, N2η(σκ2+2)+1 } . Then, we have
Λ(wk)− Ei(k)Λ(wk+1) ≥ c1 N ‖wk − T (wk)‖2.

heading: 4.2 ALMOST SURE CONVERGENCE RESULTS AND EXPECTED ITERATION COMPLEXITY
text: In this subsection, we establish the almost sure (i.e., with probability 1) convergence results for N-RPDC.
Theorem 4.1 (almost sure convergence result). Under the setting of Lemma 4.1, then
(a) limk→+∞ Φ(wk) = 0 almost surely.
(b) The sequence {wk} is almost surely bounded. (c) Any cluster point of the sequence {wk} is almost surely a stationary point of (P).
In addition, we can also provide the almost surely asymptotic rate of convergence in the sense of limit inferior for our N-RPDC method. Theorem 4.2 (almost sure O(1/ √ k) convergence rate). Under the setting of Lemma 4.1, we have
lim inf k→+∞
√ k + 1 · Φ(wk) = 0 almost surely
i.e., the asymptotic O(1/ √ k) convergence rate holds almost surely.
The almost sure convergence rate established above is of asymptotic nature, which concerns the behavior of N-RPDC when k → ∞. Next, we also provide the iteration complexity results of N-RPDC in expectation.
Theorem 4.3 (expected iteration complexity). Under the setting of Lemma 4.1. The sequence {wk} is generated by N-RPDC. Then
min 0≤k≤t
EFtΦ(wk) ≤
√ Nc2(Λ(w0)− F ∗)
c1(t+ 1) .
Consequently, to achieve EFtΦ(wt̄) ≤ ε for some 0 ≤ t̄ ≤ t, N-RPDC needs at most t = Nc2(Λ(w0)−F∗)
c1ε2 number of iterations.
Remark 4.1. Note that the closely related work (Zhang & Luo, 2020b) provides a similar iteration complexity result (i.e., O(ε2)) if g ≡ 0. By contrast, our results can deal with the additional nonsmooth nonconvex term g. In addition, we also provide the almost sure convergence results, complementing the results in (Zhang & Luo, 2020b). Let us also emphasize that the main motivation of our algorithm lies in that N-RPDC is a randomized coordinate method, which is more suitable for modern large-scale problems.
Remark 4.2. The recent work (Zhang & Luo, 2020a) establishes a global dual error bound condition, which allows them to avoid the compactness assumption on X. Due to the existence of g, we utilize the local uniform metric subregularity rather than this dual error bound condition. To further relax the compactness assumption on X in problem (P), one possible direction is to establish a global version of the utilized local uniform metric subregularity. We leave this direction as a future work.

heading: 5 NUMERICAL EXPERIMENTS
text: Because of the limitation of space, we put the numerical experiments in Appendix A.1. We can observe that N-RPDC with larger N slightly outperforms that with smaller N in terms of convergence speed. We also compared our N-RPDC to the algorithm proposed in (Zhang & Luo, 2020b); see Figure 1. We can observe that N-RPDC slightly outperforms their algorithm. Note that the main motivation of N-RPDC is that it is more suitable for modern large-scale problems due to its reasonably fast convergence speed and low computational complexity at each iteration.

heading: A APPENDIX
text: The Appendix is organized as follows. Subsection A.1 discusses implementation details for NRPDC and presents numerical experiments on non-PSD kernel SVM and sum to zero constrained LASSO problem with SCAD penalty. Some properties of the auxiliary problem (P+) are given in subsection A.2. Some examples of Bregman distance function are given in A.3. Subsection A.4 shows the proof of Proposition 3.1. Subsection A.5 provides the proof of Lemma 3.1. Proof of Lemma 3.2 is given in subsection A.6. Subsection A.7 provides the proof of Lemma 3.3. Subsection A.8 shows the proof of Lemma 4.1. Subsection A.9 provides the proof of Theorem 4.1. Subsection A.10 gives the proof of Theorem 4.2. The proof of Theorem 4.3 is given in subsection A.11.

heading: A.1 NUMERICAL EXAMPLES
text: This subsection discusses experiments conducted using MATLAB R2020a on a personal computer with Intel Core i5-6200U CPU (2.40GHz) and 8.00 GB RAM.

heading: A.1.1 NON-PSD KERNEL SUPPORT VECTOR MACHINE PROBLEM
text: Consider the non-PSD kernel support vector machine problem,
min x∈[0,c]d
1 2x >Qx− 1>d x
s.t. y>x = 0 ,
where x ∈ Rd are the decision variables, and Q ∈ Rd×d is a matrix, possibly non-PSD. Let Q = (Q>1 , Q > 2 , · · · , Q>N )> ∈ Rd×d be an appropriate partition of matrix Q and Qi be an di × d matrix. L+γ of Non-PSD kernel SVM can be written as
L+γ (w) = L+γ (x, z, p) = 1
2 x>Qx− 1>d x+
σ 2 ‖x− z‖2 + 〈p, y>x〉+ γ 2 ‖y>x‖2.
Thanks of Proposition 3.3, for problem Non-PSD kernel SVM, ∂Lzγ is locally metric subregular in (x, z, p) uniformly in z at point 0. N-RPDC with D(x, xk) = 12‖x − x
k‖2 for non-PSD kernel support vector machine problem is
pk+1 ← pk + ηy>xk; Choose i(k) from {1, 2, . . . , N} with equal probability xk+1 ← min
x∈[0,c]d 〈Qi(k)xk + σ(xk − zk)i(k) − 1di(k) + yi(k)
( pk+1 + γy>xk ) , xi(k)〉
+ 1
2αx ‖x− xk‖2;
zk+1 ← min z∈Rd 〈σ(zk − xk)i(k), zi(k)〉+
1
2αz ‖z − zk‖2.
Thus, the primal subproblem of N-RPDC has the closed form xk+1i(k) = Π[0,c]di(k) [ xki(k) − αx ( Qi(k)x k + σ(xk − zk)i(k) − 1ni(k) + (pk + γy>xk)yi(k) )] ,
xk+1j 6=i(k) = x k j 6=i(k);
and { zk+1i(k) = z k i(k) − αzσ(z
k − xk)i(k), zk+1j 6=i(k) = z k j 6=i(k).
We used LIBSVM dataset heart scale and ionosphere in the experiment. Q was generated using the sigmoid kernel (Lin & Lin, 2003), and we selected c = 1. We compared two algorithms: nonconvex ADMM in (Zhang & Luo, 2020b) (nADMM), and N-RPDC from this paper on these two datasets.
For N-RPDC algorithm, we partitioned the variables N = 10, 90 blocks for heart scale and N = 10, 70 blocks for ionosphere.
In order to verify the theoretical results of this paper, we compute the sequence of reference point {T (wk)} along with the sequence {wk} generated by N-RPDC. In the real world applications of N-RPDC, the reference point sequence is not necessary to compute. In Figure 1, the left graph show the number of blocks and min0≤k≤t Φ(wk) with respect to epoch count for heart scale; the right graph show the number of blocks and min0≤k≤t Φ(wk) with respect to epoch count for ionosphere. For both datasets, min0≤k≤t Φ(wk) of each epoch of N-RPDC is the average value of 50 times result. Moreover, we can draw the conclusion that N-RPDC is slightly outperforms nADMM and larger N is slightly outperforms small N .

heading: A.1.2 SUM TO ZERO CONSTRAINED LASSO PROBLEM WITH SCAD PENALTY
text: Consider the sum to zero constrained LASSO problem with SCAD penalty:
(CLASSO-SCAD) min x∈[−1,1]d 1 2‖Ax− b‖ 2 + n∑ i=1 φ(xi)
s.t. 1>n x = 0; ,
where x ∈ Rd is the decision variables and A ∈ Rn×d is the regressors matrix. Let A = (A1, A2, · · · , AN ) ∈ Rn×d be an appropriate partition of matrix A and Ai be an n × di matrix. b ∈ Rn is the response vector. 1d is an d-dimensional vector of 1s. φ(·) is the SCAD penalty:
φ(|t|) =  |t| |t| ≤ 1 −t2+2θ|t|−1
2(θ−1) 1 < |t| ≤ θ 1+θ 2 |t| > θ. , t ∈ R
L+γ of (CLASSO-SCAD) can be written as
L+γ (w) = L+γ (x, z, p) = 1
2 ‖Ax− b‖2 + n∑ i=1 φ(xi) + σ 2 ‖x− z‖2 + 〈p,1>n x〉+ γ 2 ‖1>n x‖2.
Thanks of Proposition 3.3, for problem (CLASSO-SCAD), ∂Lzγ is locally metric subregular in (x, z, p) uniformly in z at point 0. We take D(x, xk) = 12‖x − x
k‖2 and construct N-RPDC algorithm for (CLASSO-SCAD) as following:
pk+1 ← pk + η(1>d xk), Choose i(k) from {1, 2, . . . , N} with equal probability xk+1 ← min
x∈[−1,1]d 〈A>i(k)(Ax k − b) + σ(xk − zk)i(k) + 1di(qk), xi(k)〉
+
∑i(k) j=1 dj∑
i=1+ ∑i(k)−1 j=1 dj
φ(xi) + 1
2αx ‖x− xk‖2;
zk+1 ← min z∈Rd 〈σ(zk − xk)i(k), zi(k)〉+
1
2αz ‖z − zk‖2;
where qk = pk+1 + γ(1>d x k). Thus, the primal subproblem of RPDC has the closed form xk+1i(k) =  Π [−1,1]di(k) {sign(ζ k) max(0, |ζk| − µ1di(k))} |ζk| ≤ 1 + µ Π [−1,1]di(k) { 1 θ−1−µ [(θ − 1)ζ k − µθsign(ζk)]} 1 + µ < |ζk| ≤ θ Π [−1,1]di(k) {ζ k} |ζk| > θ
xk+1j 6=i(k) = x k j 6=i(k);
with ζk = xki(k) − αx[A > i(k)(Ax k − b) + σ(xk − zk)i(k) + qk1ni(k) ] and µ = αx.{ zk+1i(k) = z k i(k) − αzσ(z
k − xk)i(k); zk+1j 6=i(k) = z k j 6=i(k).
In this experiment, the elements of A ∈ Rn×d are selected i.i.d. from a Gaussian N (0, 1) distribution. To construct a sparse true solution x∗ ∈ Rd, given the dimension d and sparsity s, we select s entries of x∗ at random to be nonzero and N (0, 1) normally distributed, and set the rest to zero. The measurement vector b ∈ Rn is obtained by b = Ax∗ + δ with the elements of the noise vector δ ∈ Rn are i.i.d. N (0, 0.001). We choose θ = 2.3 with n = 360, d = 1280, s = 8 for the first dataset and n = 720, d = 2560, s = 16 for the second dataset.
We partitioned the variables N = 10, 40, 80 blocks. Thus di = 128, 32, 16 for the first dataset and di = 256, 64, 32 for the second dataset.
In order to verify the theoretical results of this paper, we compute the sequence of reference point {T (wk)} along with the sequence {wk} generated by N-RPDC. In the real world applications of N-RPDC, the reference point sequence is not necessary to compute. In Figure 2, the left graph
show the number of blocks and min 0≤k≤t Φ(wk) with respect to epoch count for the first dataset; the
right graph show the number of blocks and min 0≤k≤t Φ(wk) with respect to epoch count for the second
dataset. For both datasets, min 0≤k≤t Φ(wk) for each epoch of N-RPDC is the average value of 50 times result. Additionally, we can draw the conclusion that larger N is slightly outperforms small N .

heading: A.2 SOME PROPERTIES OF THE AUXILIARY PROBLEM P+
text: The auxiliary problem can be viewed as a quadratically regularized version of the original problem (P) with an additional variable z. Recall that f has Lf -Lipschitz gradient and g is ρg-weakly convex. We have the following benign properties of the auxiliary problem (P+):
(1) f(x) + σ2 ‖x− z‖ 2 is continuously differentiable in (x, z) with Lipschitz continuous gradient;
(2) F+(x, z) = f(x) + g(x) + σ2 ‖x− z‖ 2 is weakly convex in (x, z) with parameter Lf + ρg;
(3) F+(x, z) is bi-strongly convex, i.e., F+(x, z) is strongly convex in x for every fixed z and F+(x, z) is strongly convex in z for every fixed x.
The first property is trivially true. To see the second property, we note that the Lf -Lipschitz gradient property of f implies that f(x) ≥ f(y) + 〈s, x − y〉 − Lf2 ‖x − y‖
2,∀x, y ∈ Rd (Nesterov, 2003, Lemma 1.2.3). This, together with (Vial, 1983, Proposition 4.8), establishes that f is Lf -weakly convex2. By (Vial, 1983, Proposition 4.1), we have that f + g is weakly convex with parameter Lf + ρg . We have shown the second property y recognizing that (x, z) 7→ σ2 ‖x − z‖
2 is a convex mapping. The third property is a direct consequence of the definition of weak convexity, which yields that f(x) + g(x) + σ2 ‖x− z‖ 2 is σ − (Lf + ρg) strongly convex in x.

heading: A.3 SOME EXAMPLES OF BREGMAN DISTANCE FUNCTION
text: The proximal term D(x, xk) used in N-RPDC is the so-called Bregman distance function, which is defined as D(x, xk) := K(x)−K(xk)− 〈 ∇K(xk), x− xk 〉 .
Here, K : Rd → R is the core function, which is βK-strongly convex and has LK-Lipschitz continuous gradient on Rd. One popular choice of the core function isK(x) = 12‖x‖
2, which yields D(x, xk) = 12‖x − x
k‖2. However, more generally, the core function can also be selected as the Q-quadratic proximal regularization term K(x) = 12‖x‖ 2 Q, where Q is a positive definite matrix. This choice covers some second-order methods as we can setQ to be the regularized Hessian matrix.
2Actually, this derivation verifies the fact that any continuously differentiable function with L-Lipschitz continuous gradient is L-weakly convex.

heading: A.4 PROOF OF PROPOSITION 3.1
text: Proof. The optimal condition of
Tx(w k) = arg min x∈X 〈∇f(xk) + σ(xk − zk), x〉+ g(x) + 〈Tp(wk) + γ(Axk − b), Ax〉+ D(x,x k) αx
and Tz(wk) = arg min z∈Rd 〈σ(zk − xk), z〉+ 12αz ‖z − z k‖2 yields that
0 ∈ ∇f(xk) + ∂g(Tx(wk)) + σ(xk − zk) +NX(Tx(wk)) +A>[Tp(wk) + γ(Axk − b)]
+ 1
αx [∇K(Tx(wk))−∇K(xk)];
0 = σ(zk − xk) + 1 αz (Tz(w k)− zk).
Let ξx = ∇f(Tx(wk))−∇f(xk) + σ(Tx(wk)− Tz(wk))− σ(xk − zk)− γA>(Axk − b)
− 1 αx [∇K(Tx(wk))−∇K(xk)]
= ∇f(Tx(wk))−∇f(xk) + σ(Tx(wk)− Tz(wk))− σ(xk − zk)− γ
η A>(Tp(w k)− pk)
− 1 αx [∇K(Tx(wk))−∇K(xk)]; (15)
ξz = σ(Tz(w k)− Tx(wk))− σ(zk − xk)−
1
αz (Tz(w
k)− zk); (16)
and
ξp = 1
η (Tp(w
k)− pk) +A(Tx(wk)− xk). (17)
Then we have that ξx ∈ ∇f(Tx(wk)) + ∂g(Tx(wk)) + σ(Tx(wk)− Tz(wk)) +NX(Tx(wk)) +A>Tp(wk); ξz = σ(Tz(w
k)− Tx(wk)); ξp = ATx(w k)− b;
By the expression (15-17) of ξ =
( ξx ξz ξp ) , αx, αz and η are given posi-
tive numbers and Assumptions of problem (P), there exists positive number c ≥√ 3 max { LK αx + Lf + 2σ + ‖A‖, 1αz + 2σ, γ‖A‖+1 η } and the following inequality holds
‖ξ‖ ≤ c‖wk − T (wk)‖. Combining ξ ∈ ∂L+(T (wk)) and it follow the desired statement.

heading: A.5 PROOF OF LEMMA 3.1
text: In order to provide the Lemma 3.1, we need the following two propositions as preparation.
Proposition A.1 (Lipschitz properties of x(z) and x(z, p)). Suppose σ > Lf + ρg . Then
(i) ‖x(z)− x(z′)‖ ≤ σσ−(Lf+ρg)‖z − z ′‖;
(ii) ‖x(z, p)− x(z′, p)‖ ≤ σσ−(Lf+ρg)‖z − z ′‖.
Proof. (i) Let φ(x, z) = f(x) + g(x) + σ2 ‖x− z‖ 2 + IX̃(x) with X̃ = {x ∈ X|Ax− b = 0}. Since σ > Lf + ρg , we have φ(x, z) is σ − (Lf + ρg)-strongly convex in x. Therefore φ (x(z), z′)− φ (x(z′), z′)
= [φ (x(z), z′)− φ (x(z), z)] + [φ (x(z), z)− φ (x(z′), z)] + [φ (x(z′), z)− φ (x(z′), z′)]
= σ
2 [‖x(z)− z′‖2 − ‖x(z)− z‖2] + [φ (x(z), z)− φ (x(z′), z)] + σ
2 [‖x(z′)− z‖2 − ‖x(z′)− z′‖2]
= [φ (x(z), z)− φ (x(z′), z)] + σ〈z′ − z, x(z′)− x(z)〉
≤ −σ − (Lf + ρg) 2 ‖x(z)− x(z′)‖2 + σ〈z′ − z, x(z′)− x(z)〉
Again using the strongly convex of φ(x, z) in x, we have that
φ (x(z), z′)− φ (x(z′), z′) ≥ σ − (Lf + ρg) 2 ‖x(z)− x(z′)‖2.
Therefore
‖x(z)− x(z′)‖ ≤ σ σ − (Lf + ρg) ‖z − z′‖.
(ii) Let ϕ(x, z, p) = L+γ (x, z, p) +IX(x). Since σ > Lf +ρg , function ϕ(x, z, p) also is σ− (Lf + ρg)-strongly convex in x. Therefore
ϕ (x(z, p), z′, p)− ϕ (x(z′, p), z′, p) = [ϕ (x(z, p), z′, p)− ϕ (x(z, p), z, p)] + [ϕ (x(z, p), z, p)− ϕ (x(z′, p), z, p)]
+[ϕ (x(z′, p), z, p)− ϕ (x(z′, p), z′, p)]
= σ
2 [‖x(z, p)− z′‖2 − ‖x(z, p)− z‖2] + [ϕ (x(z, p), z, p)− ϕ (x(z′, p), z, p)] + σ
2 [‖x(z′, p)− z‖2 − ‖x(z′, p)− z′‖2]
= [ϕ (x(z, p), z, p)− ϕ (x(z′, p), z, p)] + σ〈z′ − z, x(z′, p)− x(z, p)〉
≤ −σ − (Lf + ρg) 2 ‖x(z, p)− x(z′, p)‖2 + σ〈z′ − z, x(z′, p)− x(z, p)〉
Again using the strongly convex of ϕ(x, z, p) in x, we have that
ϕ (x(z, p), z′, p)− ϕ (x(z′, p), z′, p) ≥ σ − (Lf + ρg) 2 ‖x(z, p)− x(z′, p)‖2.
Therefore, ‖x(z, p)− x(z′, p)‖ ≤ σσ−(Lf+ρg)‖z − z ′‖.
Proposition A.2. Suppose σ > Lf + ρg . Let {wk} be generated by N-RPDC. For given wk, the primal output (xk+1, zk+1) are random variables, and the following assertion hold:
‖xk − x(zk, Tp(wk))‖ ≤
( Lf + σ + γ‖A‖2 + LKαx
σ − Lf − ρg + 1
) ‖xk − Tx(wk)‖. (18)
Proof. Since x(zk, Tp(wk)) = arg min x∈X L+γ (x, zk, Tp(wk)) and L+γ (x, z, p) is σ−Lf −ρg-strongly convex in x, then ∀ζ ∈ ∂g(Tx(wk)) +NX(Tx(wk)) we have that (σ − Lf − ρg)‖Tx(wk)− x(zk, Tp(wk))‖2
≤ 〈∇f(Tx(wk)) + σ(Tx(wk)− zk) + ζ, Tx(wk)− x(zk, Tp(wk))〉 +〈A>[Tp(wk) + γ(ATx(wk)− b)], Tx(wk)− x(zk, Tp(wk))〉.
It follows that (σ − Lf − ρg)‖Tx(wk)− x(zk, Tp(wk))‖ ≤ ‖∇f(Tx(wk)) + σ(Tx(wk)− zk) + ζ +A>[Tp(wk) + γ(ATx(wk)− b)]‖. Again using the optimal condition of
Tx(w k) = arg min x∈X 〈∇f(xk) + σ(xk − zk), x〉+ g(x) + 〈Tp(wk) + γ(Axk − b), Ax〉+ D(x,x k) αx ,
we have that 0 ∈ ∇f(xk) + ∂g(Tx(wk)) + σ(xk − zk) +NX(Tx(wk)) +A>[Tp(wk) + γ(Axk − b)]
+ 1
αx [∇K(Tx(wk))−∇K(xk)].
Therefore ξ ∈ ∇f(Tx(wk)) + σ(Tx(wk)− zk) + ∂g(Tx(wk)) +NX(Tx(wk))
+A>[Tp(w k) + γ(ATx(w k)− b)], with
ξ = ∇f(Tx(wk)) − ∇f(xk) + σ(Tx(wk) − xk) + γA>A(Tx(wk) − xk) − 1αx [∇K(Tx(w k)) − ∇K(xk)]. By Assumptions on (P), we obtain that
‖ξ‖ ≤ (Lf + σ + γ‖A‖2 + LK αx )‖xk − Tx(wk)‖,
and
‖Tx(wk)− x(zk, Tp(wk))‖ ≤ Lf + σ + γ‖A‖2 + LKαx
σ − Lf − ρg ‖xk − Tx(wk)‖.
Then by the triangle inequality, above inequality yields the desired result.
With these two propositions in hands, we can provide the technical proof of Lemma 3.1.

heading: Proof of Lemma 3.1:
text: Step 1: Estimate the term L+γ (wk)− Ei(k)L+γ (wk+1). By the dual update of N-RPDC, we have that
L+γ (xk, zk, pk)− L+γ (xk, zk, pk+1) = 〈pk − pk+1, Axk − b〉 = − 1
η ‖pk − pk+1‖2. (19)
Since xk+1 is the solution of subproblem of x in N-RPDC scheme, we have that
〈∇i(k)f(xk) + σ(xk − zk)i(k), (x− xk+1)i(k)〉+ gi(k)(xi(k))− gi(k)(xk+1i(k) )
+〈pk+1 + γ(Axk − b), Ai(k)(x− xk+1)i(k)〉+ 1 αx D(x, xk)− 1 αx D(xk+1, xk) ≥ 0 (20)
Take x = xk in (20), by the fact that gi(k)(xki(k))− gi(k)(x k+1 i(k) ) = g(x k)− g(xk+1)
and 〈∇i(k)f(xk) + σ(xk − zk)i(k) + (Ai(k))>[pk+1 + γ(Axk − b)], (xk − xk+1)i(k)〉
= 〈∇f(xk) + σ(xk − zk) +A>[pk+1 + γ(Axk − b)], xk − xk+1〉, we have that 〈∇f(xk) + σ(xk − zk), xk − xk+1〉+ g(xk)− g(xk+1) + 〈pk+1 + γ(Axk − b), A(xk − xk+1)〉
− 1 αx D(xk+1, xk) ≥ 0
It yields that
F+(xk, zk)− F+(xk+1, zk) + 〈pk+1 + γ(Axk − b), A(xk − xk+1)〉
≥ 1 αx D(xk+1, xk)− [f(xk+1)− f(xk)− 〈∇f(xk), xk+1 − xk〉]− σ 2 ‖xk − xk+1‖2 ≥ βK − αx(Lf + σ) 2αx ‖xk − xk+1‖2
Since 〈γ(Axk − b), A(xk − xk+1)〉 = γ2 ( ‖Axk − b‖2 − ‖Axk+1 − b‖2 + ‖A(xk − xk+1)‖2 ) , therefore
L+γ (xk, zk, pk+1)− L+γ (xk+1, zk, pk+1) ≥ βK − αx(Lf + σ + γ‖A‖2)
2αx ‖xk − xk+1‖2. (21)
Additionally, since zk+1 is the solution of subproblem of z in N-RPDC, we have that
〈σ(zk − xk)i(k), (z − zk+1)i(k)〉+ 1
αz 〈zk+1 − zk, z − zk+1〉 ≥ 0. (22)
Take z = zk in (22), by the fact that 〈σ(zk − xk)i(k), (zk − zk+1)i(k)〉 = 〈σ(zk − xk), zk − zk+1〉, we have that
〈σ(zk − xk), zk − zk+1〉 ≥ 1 αz ‖zk − zk+1‖2.
Since 〈σ(zk − xk), zk − zk+1〉 = 〈σ(zk − xk+1), zk − zk+1〉+ 〈σ(xk+1 − xk), zk − zk+1〉
= σ
2
( ‖xk+1 − zk‖2 − ‖xk+1 − zk+1‖2 + ‖zk − zk+1‖2 ) +〈σ(xk+1 − xk), zk − zk+1〉
≤ σ 2
( ‖xk+1 − zk‖2 − ‖xk+1 − zk+1‖2 + ‖zk − zk+1‖2 ) + σ
2
( ‖xk − xk+1‖2 + ‖zk − zk+1‖2 ) = σ
2
( ‖xk+1 − zk‖2 − ‖xk+1 − zk+1‖2 ) + σ‖zk − zk+1‖2
+ σ
2 ‖xk − xk+1‖2,
above inequality yields that
L+γ (xk+1, zk, pk+1)− L+γ (xk+1, zk+1, pk+1) = σ
2
( ‖xk+1 − zk‖2 − ‖xk+1 − zk+1‖2 ) ≥ ( 1
αz − σ)‖zk − zk+1‖2 − σ 2 ‖xk − xk+1‖2.
(23) By the combination of (19), (21) and (23), we obtain that
L+γ (xk, zk, pk)− L+γ (xk+1, zk+1, pk+1)
≥ βK αx − (Lf + 2σ + γ‖A‖2)
2 ‖xk − xk+1‖2 + ( 1 αz − σ)‖zk − zk+1‖2 − 1 η ‖pk − pk+1‖2.
Take expectation with respect to i(k) on both side of above inequality, we have that
L+γ (xk, zk, pk)− Ei(k)L+γ (xk+1, zk+1, pk+1)
≥ βK αx − (Lf + 2σ + γ‖A‖2)
2 Ei(k)‖xk − xk+1‖2 + (
1
αz − σ)Ei(k)‖zk − zk+1‖2
−1 η ‖pk − pk+1‖2.

heading: By the fact that Ei(k)‖xk−xk+1‖2 = 1N ‖x
text: k−Tx(wk)‖2, Ei(k)‖zk−zk+1‖2 = 1N ‖z k−Tz(wk)‖2 and pk+1 = Tp(wk), above inequality follows that
L+γ (xk, zk, pk)− Ei(k)L+γ (xk+1, zk+1, pk+1)
≥ βK αx − (Lf + 2σ + γ‖A‖2)
2N ‖xk − Tx(wk)‖2 +
1
N (
1
αz − σ)‖zk − Tz(wk)‖2
−1 η ‖pk − Tp(wk)‖2.
Step 2: Estimate the term ν(zk) − ν(zk+1). By the Danskin’s theorem, we have that ∇ν(z) = σ (z − x(z)). By statement (i) of Proposition A.1, we have that
‖∇ν(zk)−∇ν(zk+1)‖ ≤ (
σ2
σ − Lf − ρg + σ
) ‖zk − zk+1‖.
The gradient Lipschitz property of ν follows that ν(zk)− ν(zk+1) ≥ σ〈zk − x(zk), zk − zk+1〉 − (
σ2
2(σ − Lf − ρg) + σ 2
) ‖zk − zk+1‖2. (24)
Step 3: Estimate the term ψγ(zk+1, pk+1) − ψγ(zk, pk). Since x(zk, pk) = arg min
x∈X L+γ (x, zk, pk), we have
ψγ(z k+1, pk+1)− ψγ(zk, pk+1)
= L+γ (x(zk+1, pk+1), zk+1, pk+1)− L+γ (x(zk, pk+1), zk, pk+1) ≥ L+γ (x(zk+1, pk+1), zk+1, pk+1)− L+γ (x(zk+1, pk+1), zk, pk+1) = σ
2
( ‖x(zk+1, pk+1)− zk+1‖2 − ‖x(zk+1, pk+1)− zk‖2 ) ≥ σ〈zk − x(zk+1, pk+1), zk+1 − zk〉. (25)
By the concavity of ψγ(z, ·) and∇ψγ(z, p) = Ax(z, p)− b, we have ψγ(z k, pk+1)− ψγ(zk, pk) ≥ 〈pk+1 − pk, Ax(zk, pk+1)− b〉
= 1 2η ‖pk − pk+1‖2 + η 2 ‖Ax(zk, pk+1)− b‖2
− 1 2η ‖pk+1 − pk − η[Ax(zk, pk+1)− b]‖2.
By the dual update of N-RPDC, the above inequality yields that
ψγ(z k, pk+1)− ψγ(zk, pk)
≥ 1 2η ‖pk − pk+1‖2 + η 2 ‖Ax(zk, pk+1)− b‖2 − 1 2η ‖ηA[xk − x(zk, pk+1)]‖2
≥ 1 2η ‖pk − pk+1‖2 + η 2 ‖Ax(zk, pk+1)− b‖2 − η 2 ‖A‖2 · ‖xk − x(zk, pk+1)‖2. (26)
By Proposition A.2, (26) guarantees that
ψγ(z k, pk+1)− ψγ(zk, pk) ≥
1 2η ‖pk − pk+1‖2 + η 2 ‖Ax(zk, pk+1)− b‖2
−η 2 ‖A‖2
( Lf + σ + γ‖A‖2 + LKαx
σ − Lf − ρg + 1
)2 ‖xk − Tx(wk)‖2.
(27) By the combination of (25) and (27), we obtain that
ψγ(z k+1, pk+1)− ψγ(zk, pk)
≥ σ〈zk − x(zk+1, pk+1), zk+1 − zk〉+ 1 2η ‖pk − pk+1‖2
+ η 2 ‖Ax(zk, pk+1)− b‖2 − η 2 ‖A‖2
( Lf + σ + γ‖A‖2 + LKαx
σ − Lf − ρg + 1
)2 ‖xk − Tx(wk)‖2
(28)
Step 4: Estimate the term [ν(zk)− ψγ(zk, pk)]− Ei(k)[ν(zk+1)− ψγ(zk+1, pk+1)]. Summing statement (24) and (28), we obtain the variation of dual gap:
[ν(zk)− ψγ(zk, pk)]− [ν(zk+1)− ψγ(zk+1, pk+1)] ≥ σ〈x(zk+1, pk+1)− x(zk), zk − zk+1〉 − ( σ2
2(σ − Lf − ρg) + σ 2
) ‖zk − zk+1‖2
+ 1 2η ‖pk − pk+1‖2 + η 2 ‖Ax(zk, pk+1)− b‖2
−η 2 ‖A‖2
( Lf + σ + γ‖A‖2 + LKαx
σ − Lf − ρg + 1
)2 ‖xk − Tx(wk)‖2. (29)
By Proposition A.1, we have that
σ〈x(zk+1, pk+1)− x(zk), zk − zk+1〉 = σ[〈x(zk, pk+1)− x(zk), zk − zk+1〉+ 〈x(zk+1, pk+1)− x(zk, pk+1), zk − zk+1〉]
≥ − σ 2λ ‖zk − zk+1‖2 − σλ 2 ‖x(zk, pk+1)− x(zk)‖2
−σ‖zk − zk+1‖ · ‖x(zk, pk+1)− x(zk+1, pk+1)‖
≥ −( σ 2λ + σ2 σ − Lf − ρg )‖zk − zk+1‖2 − σλ 2 ‖x(zk, pk+1)− x(zk)‖2, (30)
with λ > 0 is any positive number. By the combination of (29) and (30), we have that
[ν(zk)− ψγ(zk, pk)]− [ν(zk+1)− ψγ(zk+1, pk+1)] ≥ − ( σ
2λ +
3σ2
2(σ − Lf − ρg) + σ 2
) ‖zk − zk+1‖2 − σλ
2 ‖x(zk, pk+1)− x(zk)‖2
+ 1 2η ‖pk − pk+1‖2 + η 2 ‖Ax(zk, pk+1)− b‖2
−η 2 ‖A‖2
( Lf + σ + γ‖A‖2 + LKαx
σ − Lf − ρg + 1
)2 ‖xk − Tx(wk)‖2. (31)
Take expectation with respect to i(k) on both side of above inequality, we have that
[ν(zk)− ψγ(zk, pk)]− Ei(k)[ν(zk+1)− ψγ(zk+1, pk+1)] ≥ − ( σ
2λ +
3σ2
2(σ − Lf − ρg) + σ 2
) Ei(k)‖zk − zk+1‖2 − σλ
2 ‖x(zk, pk+1)− x(zk)‖2
+ 1 2η ‖pk − pk+1‖2 + η 2 ‖Ax(zk, pk+1)− b‖2
−η 2 ‖A‖2
( Lf + σ + γ‖A‖2 + LKαx
σ − Lf − ρg + 1
)2 ‖xk − Tx(wk)‖2.

heading: By the fact that Ei(k)‖zk − zk+1‖2 = 1N ‖z
text: k − Tz(wk)‖2 and pk+1 = Tp(wk), above inequality follows that [ν(zk)− ψγ(zk, pk)]− Ei(k)[ν(zk+1)− ψγ(zk+1, pk+1)]
≥ − 1 N
( σ
2λ +
3σ2
2(σ − Lf − ρg) + σ 2
) ‖zk − Tz(wk)‖2 − σλ
2 ‖x(zk, Tp(wk))− x(zk)‖2
+ 1
2η ‖pk − Tp(wk)‖2 +
η 2 ‖Ax(zk, Tp(wk))− b‖2
−η 2 ‖A‖2
( Lf + σ + γ‖A‖2 + LKαx
σ − Lf − ρg + 1
)2 ‖xk − Tx(wk)‖2. (32)
Step 5: Estimate the term Λ(wk)−Ei(k)Λ(wk+1). By the combination of Step 1 and Step 4, we obtain that
Λ(wk)− Ei(k)Λ(wk+1)
≥  βKαx − (Lf + 2σ + γ‖A‖2) 2N − η‖A‖2 ( Lf + σ + γ‖A‖2 + LKαx σ − Lf − ρg + 1 )2 ‖xk − Tx(wk)‖2 + 1
N
[ 1
αz − 2σ − σ λ − 3σ
2
σ − Lf − ρg
] ‖zk − Tz(wk)‖2
+η‖Ax(zk, Tp(wk))− b‖2 − σλ‖x(zk, Tp(wk))− x(zk)‖2. (33)
Since 0 < αx ≤ βK/[Lf + 2σ + γ‖A‖2 + 5], then βKαx − (Lf + 2σ + γ‖A‖ 2) ≥ 5. Moreover,
0 < η ≤ 1/ [ 2N‖A‖2 ( Lf+σ+γ‖A‖2+ LK αx
σ−Lf−ρg + 1
)2] implies
Λ(wk)− Ei(k)Λ(wk+1)
≥ 2 N ‖xk − Tx(wk)‖2 + 1 N
[ 1
αz − 2σ − σ λ − 3σ
2
σ − Lf − ρg
] ‖zk − Tz(wk)‖2
+η‖Ax(zk, Tp(wk))− b‖2 − σλ‖x(zk, Tp(wk))− x(zk)‖2. (34)

heading: A.6 PROOF OF LEMMA 3.2
text: Proof. The result is directly by the fact that dist ( 0, ∂Lzγ (x(z, p), p) ) = ‖Ax(z, p)− b‖.

heading: A.7 PROOF OF LEMMA 3.3
text: Proof. The claim is provided by the results of Proposition 1 and Corollary of (Robinson, 1981).

heading: A.8 PROOF OF LEMMA 4.1
text: Proof. Taking λ = min{ δη2σM2+1 , η σκ2+1} and 0 < αz < 1/ [ 2σ + σλ + 3σ2 σ−Lf−ρg + 1 ]
in Lemma 3.1, we have that
Λ(wk)− Ei(k)Λ(wk+1)
≥ 2 N ‖xk − Tx(wk)‖2 + 1 N ‖zk − Tz(wk)‖2 + η‖Ax(zk, Tp(wk))− b‖2
−σmin{ δη 2σM2 + 1 , η σκ2 + 1 }‖x(zk, Tp(wk))− x(zk)‖2. (35)
For the last two terms of above inequality, we have two cases to discuss.
• Case 1: ‖Ax(zk, Tp(wk))− b‖ ≤ 2σλM 2 η ≤ δ.
Since ∂Lzγ(x, p) is locally metric subregular in (x, p) uniformly in z over Rd at point 0 with parameters κ, δ > 0, (35) yields that
Λ(wk)− Ei(k)Λ(wk+1) ≥ 2
N ‖xk − Tx(wk)‖2 +
1
N ‖zk − Tz(wk)‖2
+ η
σκ2 + 1 ‖Ax(zk, Tp(wk))− b‖2. (36)
• Case 2: ‖Ax(zk, Tp(wk))− b‖ > 2σλM 2
η .
Since M = max x,x′∈X
‖x− x′‖, we have that
η‖Ax(zk, Tp(wk))− b‖2 − σλ‖x(zk, Tp(wk))− x∗(zk)‖2
≥ η 2 ‖Ax(zk, Tp(wk))− b‖2 + η 2 · 2σλM
2
η − σλM2
= η
2 ‖Ax(zk, Tp(wk))− b‖2. (37)
For both case, (35) yields that
Λ(wk)− Ei(k)Λ(wk+1)
≥ 2 N ‖xk − Tx(wk)‖2 + 1 N ‖zk − Tz(wk)‖2 +
η
σκ2 + 2 ‖Ax(zk, Tp(wk))− b‖2. (38)
By Proposition A.2, above inequality yields that
Λ(wk)− Ei(k)Λ(wk+1)
≥ 1 N ‖xk − Tx(wk)‖2 +
1
N
( Lf+σ+γ‖A‖2+
LK αx
σ−Lf−ρg + 1 )2 ‖xk − x(zk, Tp(wk)‖2 + 1
N ‖zk − Tz(wk)‖2 +
η
σκ2 + 2 ‖Ax(zk, Tp(wk))− b‖2
≥ 1 N ‖xk − Tx(wk)‖2 +
1 N‖A‖2 ( Lf+σ+γ‖A‖2+ LK αx
σ−Lf−ρg + 1 )2 ‖A(xk − x(zk, Tp(wk))‖2 + 1
N ‖zk − Tz(wk)‖2 +
η
σκ2 + 2 ‖Ax(zk, Tp(wk))− b‖2
≥ 1 N ‖xk − Tx(wk)‖2 + 1 N ‖zk − Tz(wk)‖2
+ η
2(σκ2 + 2)
( 2‖A(xk − x(zk, Tp(wk))‖2 + 2‖Ax(zk, Tp(wk))− b‖2 ) ≥ 1
N ‖xk − Tx(wk)‖2 +
1
N ‖zk − Tz(wk)‖2 +
η
2(σκ2 + 2) ‖Axk − b‖2
= 1
N ‖xk − Tx(wk)‖2 +
1
N ‖zk − Tz(wk)‖2 +
1
2η(σκ2 + 2) ‖pk − Tp(wk)‖2
≥ c1 N ‖wk − T (wk)‖2, (39)
with c1 = min { 1, N2η(σκ2+2)+1 } .

heading: A.9 PROOF OF THEOREM 4.1
text: In order to provide the Theorem 4.1, we need the famous Robbins-Siegmund theorem.
Theorem A.1. (Robbins-Siegmund Theorem, Robbins & Siegmund (1971), Theorem 1.3.12 in Duflo (2013), Theorem 2.27 in Carpentier et al. (2015)) Let {Λk}k∈N, {µk}k∈N, {νk}k∈N and {ηk}k∈N be four positive sequences of real-valued random variables adapted to the filtration {ξk}k∈N. Assume that
EξkΛk+1 ≤ (1 + µk)Λk + νk − ηk, ∀k ∈ N, and that ∑
k∈N µk < +∞ and ∑ k∈N νk < +∞, a.s..
Then, the sequence {Λk}k∈N almost surely converges to a finite3 random variable Λ∞, and ∑ k∈N ηk < +∞, a.s..
Then we can derive the technical proof of this theorem.

heading: Proof of Theorem 4.1:
text: (a) By Lemma 4.1 and the fact that Λ(w) ≥ F ∗, we have that
Ei(k)[Λ(wk+1)− F ∗] ≤ [Λ(wk)− F ∗]− c1 N ‖wk − T (wk)‖2.
By the Robbins-Siegmund Theorem (Theorem A.1), we obtain that, lim k→+∞ Λ(wk) is almost
surely exists and +∞∑ k=0 ‖wk − T (wk)‖2 < +∞ a.s.. Therefore, we have limk→+∞ Φ(wk) = 0;
(b) By the compact of X, we have the sequences {xk} and {x(zk)} are bounded. By statement (a), we have lim
k→+∞ Λ(wk) is almost surely exists. It follows that the sequence
{Λ(wk)} is almost surely bounded. By the definition of Λ(wk), we also have that Λ(wk) ≥ ν(zk) ≥ F (x(zk)) ≥ F ∗. Therefore the sequence {ν(zk)} is almost surely bounded. Now we show the almost surely boundness of the sequence {zk} by contradiction. Suppose P { ‖zk‖ → +∞ } > 0. Since the sequence {x(zk)} is bounded, there exists a subsequence {zk′} such that x(zk′) → x̄. Therefore P { ν(zk ′ ) = F (x(zk
′ )) + σ2 ‖(x(z
k′))− zk′‖2 → +∞ } > 0, which follows a contradiction
with that {ν(zk′)} is almost surely bounded. Therefore, we obtain the almost surely boundness of the sequence {zk}. Next we propose the boundness of {pk}. By the optimal condition of problem of Tx(wk) = arg min
x∈X 〈∇f(xk) + σ(xk − zk), x〉+ g(x) + 〈Tp(wk) + γ(Axk − b), Ax〉+ D(x,x k) αx , we have
0 ∈ ∇f(xk) + σ(xk − zk) + ∂g(Tx(wk)) +NX(Tx(wk)) +A>[Tp(wk) + γ(Axk − b)]
+ 1
αx [∇K(Tx(wk))−∇K(xk)].
Since Tx(wk) = NEi(k)xk+1 − (N − 1)xk and the boundness of the sequence {xk}, we have the sequence {Tx(wk)} is bounded. By bounded of subgradient of g, the boundness of {xk}, {Tx(wk)} and the almost surely boundness of {zk}, it follows that ‖A>Tp(wk)‖ < +∞, ∀k > 0, a.s.. By the fact that pk+1 = Tp(wk), we have that
‖A>pk+1‖ < +∞, ∀k > 0, a.s. By update of p in N-RPDC, we have that
pk+1 = p0 + p̂k+1 when p̂k+1 = η k∑ j=1 (Axj − b) = ηA k∑ j=1 (xj − x∗) ∈ Im(A),
with x∗ be the optimal solution of (P). Since p0 ∈ Null(A>), we have A>p0 = 0 and ‖A>p̂k+1‖ < +∞, a.s.. Then we have two cases to discuss. • Case 1: rank(A) = n. Obviously, we have that
‖p̂k+1‖ ≤ ‖A >p̂k+1‖√
λmin(AA>) < +∞,
with λmin(AA>) be the smallest eigenvalue of matrix AA>. • Case 2: rank(A) = r < n. Without loss of generality, assume the first r rows of A
(denoted by Ar ∈ Rr×d) are linearly independent, we have
A =
( Ar
BAr
) = ( Ir×r B ) Ar,
3A random variable X is finite if P ({ω ∈ Ω|X(ω) = +∞}) = 0.
where Ir×r ∈ Rr×r is the indentity matrix and B ∈ R(n−r)×r. Let Q := Ir×r + B>B. It is easy to show that Q is symmetric and positive definite. By the fact that p̂k+1 = ηA ∑k j=1(x j − x∗), we have that
‖A>p̂k+1‖2 = ∥∥∥∥∥∥(Ar)> (Ir×r, B>) η ( Ir×r B ) Ar k∑ j=1 (xj − x∗) ∥∥∥∥∥∥ 2
= ∥∥∥∥∥∥(Ar)>QηAr k∑ j=1 (xj − x∗) ∥∥∥∥∥∥ 2 ≥ λmin ( Ar(Ar)>
) ∥∥∥∥∥∥QηAr k∑ j=1 (xj − x∗) ∥∥∥∥∥∥ 2
≥ λmin ( Ar(Ar)> ) λmin(Q >Q)‖ηAr k∑ j=1 (xj − x∗)‖2,
where λmin ( Ar(Ar)> ) is the smallest eigenvalue of positive semidefinite matrixAr(Ar)> and λmin(Q>Q) is the smallest eigenvalue of positive semidefinite matrixQ>Q. Using the fact that
‖p̂k+1‖2 = ‖ηA k∑ j=1 (xj − x∗)‖2
= ‖η ( Ir×r B ) Ar k∑ j=1 (xj − x∗)‖2 ≤ λmax (Q) ‖ηAr k∑ j=1 (xj − x∗)‖2,
where λmax(Q) is the largest eigenvalue of matrix Q. Therefore, we have that
‖p̂k+1‖2 ≤ λmax (Q) λmin (Ar(Ar)>)λmin(Q>Q) ‖A>p̂k+1‖2.
Therefore, for both cases, there exists d > 0 such that ‖p̂k+1‖ ≤ d‖A>p̂k+1‖ < +∞, a.s.
which yields the almost surely boundness of {pk} by the triangle inequality. Therefore, the sequence {wk} is almost surely bounded.
(c) By lim k→+∞ ‖wk − T (wk)‖ = 0 a.s. in statment (a), Lemma 2.1, Proposition 3.1, the almost
surely boundness of {wk}, and the closedness of ∂L+γ (·), we obtain the desired result.

heading: A.10 PROOF OF THEOREM 4.2
text: Proof. Suppose this statement does not hold, that is, P {
lim inf k→+∞
√ k + 1Φ(wk) ≥ δ } > 0,
for some d > 0. Then for k0 large enough, for all k ≥ k0 we have P {
Φ(wk) ≥ d√ k + 1
} > 0.
Therefore
P { +∞∑ k=k0 (Φ(wk))2 ≥ +∞∑ k=k0 (Φ(wk))2 ≥ +∞∑ k=k0 d2 k + 1 } > 0.
Since the fact that +∞∑ k=k0 d2 k+1 = +∞, we have that above inequality is contradicted with the fact that
+∞∑ k=0 (Φ(wk))2 < +∞ a.s. in the proof of Theorem 4.1. Therefore, we obtain the desired result.

heading: A.11 PROOF OF THEOREM 4.3
text: Proof. Recalling the inequality of Lemma 4.1 and the fact that Λ(wk) ≥ ν(zk) ≥ F ∗, we have
Ei(k)[Λ(wk+1)− F ∗] ≤ [Λ(wk)− F ∗]− c1 c2N (Φ(wk))2.
Taking expectation with respect to Ft, t > k for above inequality, we obtain that
EFt [Λ(wk+1)− F ∗] ≤ EFt [Λ(wk)− F ∗]− c1 c2N EFt(Φ(wk))2.
By the fact Λ(wk) ≥ ν(zk) ≥ F ∗, it follows
c1 c2N t∑ k=0 EFt(Φ(wk))2 ≤ Λ(w0)− F ∗.
By the Jensen’s inequality and the convexity of function h(x) = x2, x ∈ R, it follows that
c1 c2N t∑ k=0 ( EFtΦ(wk) )2 ≤ Λ(w0)− F ∗. By the fact that EFtΦ(wk) ≥ 0 and h(x) = x2, x ∈ R is monotonic increasing in x on [0,+∞), we have that
c1 c2N (t+ 1)
( min
0≤k≤t EFtΦ(wk)
)2 ≤ Λ(w0)− F ∗.
Here comes that
min 0≤k≤t
EFtΦ(wk) ≤
√ c2N c1
(Λ(w0)− F ∗) t+ 1 .

Review summaries:
This paper proposes a randomized primal-dual coordinate method for solving linearly constrained nonsmooth nonconvex optimization problems. At each iteration, this method only selects a primal variable to update randomly. The proposed algorithm can solve large-scale problems since the computational complexity in each iteration is very cheap compared with the full gradient methods. The authors prove that any cluster point of the sequence of iterates is almost surely a stationary point. Some experiments on the non-PSD kernel SVM problem and the linearly constrained Lasso problem show that it is practical.    

This paper develops a randomized primal-dual algorithm for solving linearly constrained problems with a nonsmooth nonconvex objective and convex constraint set. Specifically, the objective is composed of a smooth nonconvex term and a weakly convex and hence nonsmooth nonconvex term. The authors show that the limit points of the algorithm are stationary points. The authors then prove almost sure asymptotic rate of $1O(/\sqrt{k})$ and complexity of $O(\epsilon^{-2})$ for obtaining an $\epsilon$-stationary point, in expectation. Some preliminary empirical results are included. 

This paper studies a particular type of nonconvex nonsmooth problem. The authors propose to use the Bertsekas' convexification procedure on the nonconvex problem, and then use randomized coordinate method to solve the corresponding augmented Lagrangian. The paper also proposes a so-called surrogate stationarity measure to establish the convergence rate for the proposed method.

The paper proposes a randomized primal-dual method for solving the nonconvex nonsmooth optimization problem with linear constraints. The main analysis relies on the idea of using auxiliary problems and the Lyapunov function proposed in  (Zhang & Luo, 2020b).  The paper shows that any cluster point of the generated sequence is almost surely a stationary point of the nonconvex nonsmooth problem and also provides an almost sure asymptotic convergence rate together with the expected iteration complexity of the proposed method.