# Decentralized Non-Convex Learning with Linearly Coupled Constraints:
Algorithm Designs and Application to Vertical Learning Problem

Jiawei Zhang , Songyang Ge , Tsung-Hui Chang , and Zhi-Quan Luo Jiawei Zhang
is with Shenzhen Research Institute of Big Data, The Chinese University of
Hong Kong, Shenzhen, China 518172. E-mail<EMAIL_ADDRESS>Ge is with Shenzhen Research Institute of Big Data, The Chinese University of
Hong Kong, Shenzhen, China 518172. E-mail<EMAIL_ADDRESS>Tsung-
Hui Chang is the corresponding author. Address: Guangdong Provincial Key
Laboratory of Big Data Computing and the School of Science and Engineering,
The Chinese University of Hong Kong, Shenzhen, China 518172. E-mail:
<EMAIL_ADDRESS>Zhi-Quan Luo is with the Shenzhen Research Institute
of Big Data, The Chinese University of Hong Kong, Shenzhen, China 518172.
E-mail<EMAIL_ADDRESS>
# Decentralized Non-Convex Learning with Linearly Coupled Constraints:
Algorithm Designs and Application to Vertical Learning Problem††thanks:
Jiawei Zhang and Songyang Ge contributed equally. Part of the work was
presented in IEEE ICASSP 2020 [1].

Jiawei Zhang , Songyang Ge , Tsung-Hui Chang , and Zhi-Quan Luo Jiawei Zhang
is with Shenzhen Research Institute of Big Data, The Chinese University of
Hong Kong, Shenzhen, China 518172. E-mail<EMAIL_ADDRESS>Ge is with Shenzhen Research Institute of Big Data, The Chinese University of
Hong Kong, Shenzhen, China 518172. E-mail<EMAIL_ADDRESS>Tsung-
Hui Chang is the corresponding author. Address: Guangdong Provincial Key
Laboratory of Big Data Computing and the School of Science and Engineering,
The Chinese University of Hong Kong, Shenzhen, China 518172. E-mail:
<EMAIL_ADDRESS>Zhi-Quan Luo is with the Shenzhen Research Institute
of Big Data, The Chinese University of Hong Kong, Shenzhen, China 518172.
E-mail<EMAIL_ADDRESS>
###### Abstract

Motivated by the need for decentralized learning, this paper aims at designing
a distributed algorithm for solving nonconvex problems with general linear
constraints over a multi-agent network. In the considered problem, each agent
owns some local information and a local variable for jointly minimizing a cost
function, but local variables are coupled by linear constraints. Most of the
existing methods for such problems are only applicable for convex problems or
problems with specific linear constraints. There still lacks a distributed
algorithm for solving such problems with general linear constraints under the
nonconvex setting. To tackle this problem, we propose a new algorithm, called
_proximal dual consensus_ (PDC) algorithm, which combines a proximal technique
and a dual consensus method. We show that under certain conditions the
proposed PDC algorithm can generate an $\epsilon$-Karush-Kuhn-Tucker solution
in $\mathcal{O}(1/\epsilon)$ iterations, achieving the lower bound for
distributed non-convex problems up to a constant. Numerical results are
presented to demonstrate the good performance of the proposed algorithms for
solving two vertical learning problems in machine learning over a multi-agent
network.

##  1 Introduction

Recently, distributed optimization has attracted significant attention in
signal processing, control and machine learning societies [2, 3, 4] due to the
need of achieving decentralized control and decision making in sensor
networks, and large-scale data processing and learning. Since the distributed
agents/nodes process only local data and messages exchanged from direct
neighbors, distributed optimization is well suited for applications where the
data are acquired distributively [5] or where data privacy is of primary
concern [6].

In this paper we are interested in distributed optimization methods for
solving the following problem

$\displaystyle{\sf(P)}~{}~{}\min_{\begin{subarray}{c}{\bm{x}}_{i}\in\mathbb{R}^{n},\\\
i=1,\cdots,N\end{subarray}}~{}$
$\displaystyle\sum_{i=1}^{N}f_{i}({\bm{x}}_{i})~{}~{}\text{s.t.}~{}~{}\sum_{i=1}^{N}{\bm{B}}_{i}{\bm{x}}_{i}={\bm{q}},~{}~{}$
(1)

where $N$ denotes the number of agents,
$f_{i}:\mathbb{R}^{n}\rightarrow\mathbb{R}$ are smooth and possibly non-convex
functions, and ${\bm{B}}_{i}\in\mathbb{R}^{M\times n}$ and
${\bm{q}}\in\mathbb{R}^{M}$ are constants. The target is to solve (P) over a
multi-agent network with $N$ agents. It is assumed that each agent $i$ knows
only $f_{i}$, ${\bm{B}}_{i}$ and ${\bm{q}}$, and can only communicate with its
neighboring agents. The agent variables are coupled due to the linear
constraint $\sum_{i=1}^{N}{\bm{B}}_{i}{\bm{x}}_{i}={\bm{q}}$. Such problem
arises, for example, in machine learning with distributed features [7, 8],
power control in electrical power networks [9, 10], interference management in
wireless networks [11], and the network utility maximization problem [12].

###  1.1 Related Works

As a linearly constrained problem, problem (P) can be easily solved 111 Since
finding a global minimum of problem (P) is NP-hard in general, what we mean
“solved” here is to reach a Karush-Kuhn-Tucker (KKT) solution of problem (P).
under the centralized setting, where all information are gathered in a central
node. Centralized methods for solving linearly constrained problems are well
studied in the literature. For example, the penalty method [13] penalizes the
linear constraints in the objective function and minimizes the penalized
objective function instead. The augmented Lagrangian (AL) method [14] uses a
min-max approach to solve the resulting AL function. Problem (P) can also be
efficiently handled by the alternating direction method of multipliers (ADMM)
[15, 16, 17]. However, solving (P) in a distributed fashion over a multi-agent
network is challenging, due to the nonconvexity and coupled variables. In the
existing literature, most of the distributed methods for solving (P) assumed
convexity of $f_{i}$’s. For example, if (P) is convex, it is well known that
the dual decomposition method [18] is an efficient method to solve such
problem over a star network where all agents can communicate with a central
node. For a multi-agent network with general graph structures, convexity and
strong duality are the keys to most of the existing distributed methods for
solving (P); see [7, 19, 20, 21, 22, 23, 24, 25] and references therein.

Interestingly, when ${\bm{B}}=[{\bm{B}}_{1},\cdots,{\bm{B}}_{N}]$ is the
incidence matrix of the network graph and ${\bm{q}}$ is the zero vector, (P)
corresponds to a consensus formulation of the finite sum minimization problem
$\min_{{\bm{x}}}\sum_{i=1}^{N}f_{i}({\bm{x}})$, and there exist rich results
in the literature for both convex and nonconvex cases, see, for example, the
recent surveys [4, 26, 27]. Specifically, it is worth pointing out that the
recent work [28] has shown that, for solving the non-convex finite sum problem
in the distributed setting, the first-order algorithms have an iteration
complexity lower bound of $\Omega(1/\sqrt{\xi(\mathcal{G})}\times
1/\epsilon)$, where $\xi(\mathcal{G})$ is the spectral gap of the graph
Laplacian matrix and $\epsilon$ is the solution accuracy to a stationary
point. The authors of [28] also devised an optimal distributed algorithm that
can attain the lower bound up to a factor. However, since these decentralized
methods are designed for distributed finite sum problems, they cannot be used
to solve (P) which has a general linear equality constraint.

###  1.2 Contributions

In this paper, we propose a new distributed algorithm for handling (P) under
the general linearly coupled constraint and non-convex setting. The proposed
method, which we call the _proximal dual consensus_ (PDC) method [1]222
Compared with [1], this manuscript presents all the proof details and stronger
theoretical results. Furthermore, the current manuscript presents an inexact
update version of the proposed PDC algorithm as well as new experiment
results., is motivated by the proximal technique recently proposed in [17] and
the dual consensus method in [7, 19]. Specifically, a strongly convex
counterpart of (P) is obtained by introducing a quadratic proximal term
centered at an auxiliary primal variable. Then, the dual consensus approach is
applied to the proximal problem, resulting in a distributed algorithm for (P).
In the proposed algorithm, the agents can choose to solve their local
subproblem exactly or inexactly by simply performing a gradient descent.

While the proposed PDC algorithms are seemingly a simple combination of [17]
and [7], proving their convergence and convergence rate are far from easy as
several new error bounds are needed. Theoretically, we show that the proposed
PDC algorithms can converge to a Karush-Kuhn-Tucker (KKT) solution of (P), and
generate an $\epsilon$-optimal solution in $\mathcal{O}(1/\epsilon)$
iterations. This meets the iteration complexity lower bound (up to a constant)
of the first-order distributed algorithms for solving a nonconvex problem [29,
28]333 As mentioned, since the finite sum problem is a special instance of
problem (P), the iteration complexity lower bound in [28] also serves as a
lower bound for (P). . To evaluate the practical performance of the proposed
algorithms, we consider the less studied distributed vertical learning problem
[30, 31]. The presented numerical results show that the proposed algorithms
perform well in a non-convex regularized logistic regression problem and a
neural network based distributed classification problem under the vertical
learning setting. To the best of our knowledge, the proposed distributed
algorithms are the first to solve (P) and the distributed vertical learning
problem under the non-convex setting over the multi-agent network.

Synopsis: Section 2 presents two application examples of (P) about machine
learning from distributed features. The proposed PDC algorithms are developed
and analyzed in Section 3. The proof details of main theoretical results are
shown in Section 4. Numerical results are given in Section 5 and conclusions
are drawn in Section 6.

Notation:  ${\bf I}_{n}$ is denoted as the $n$ by $n$ identity matrix.
$\langle{\bf a},{\bm{b}}\rangle$ represents the inner product of vectors ${\bf
a}$ and ${\bm{b}}$, $\|{\bf a}\|$ is the Euclidean norm, and
$\|{\bm{a}}\|_{{\bf A}}^{2}$ represents ${\bm{a}}^{\top}{{\bf A}}{\bm{a}}$;
for a matrix ${\bf A}$, both $\|{\bf A}\|$ and $\sigma_{\max}({\bf A})$ denote
its matrix 2-norm (maximum singular value), and $\sigma_{\min}({\bf A})$
denotes the minimum non-zero singular value of ${\bf A}$; $\otimes$ denotes
the Kronecker product; $i\in[N]$ denotes $i\in\\{1,\ldots,N\\}$.

##  2 Application Examples

Figure 1:  Data partition for horizontal learning (Left) and vertical learning
(Right).

Problem (P) appears in many engineering applications. As mentioned in Section
1, one classical example is the consensus formulation of the finite sum
problem $\min_{{\bm{x}}}\sum_{i=1}^{N}f_{i}({\bm{x}})$. Specifically, consider
a multi-agent network with $N$ connected agents, and each agent $i$,
$i\in[N]$, has a local copy of ${\bm{x}}$, denoted as
${\bm{x}}_{i}\in\mathbb{R}^{n}$. Suppose that
$\tilde{}{\bm{A}}\in\mathbb{R}^{N\times N}$ is the incidence matrix of the
network graph. One can have the following consensus formulation of the finite
sum problem

$\displaystyle\min_{{\bm{x}}=[{\bm{x}}_{1}^{\top},\cdots,{\bm{x}}_{N}^{\top}]^{\top}}$
$\displaystyle~{}\sum_{i=1}^{N}f_{i}({\bm{x}}_{i})$ (2a) s.t.
$\displaystyle(\tilde{}{\bm{A}}\otimes{\bm{I}}_{n}){\bm{x}}=\mathbf{0},$ (2b)

where (2b) ensures all ${\bm{x}}_{i}$’s are the same to each other. In the
machine learning literature, the finite sum problem corresponds to the
conventional _horizontal learning_ problem where the data samples are
partitioned and distributively owned by the agents [4] (left of Fig. 1), and
such problem has been extensively studied [4, 26, 27]. Unlike (2), in this
section we discuss an important but less discussed learning problem: the
_vertical learning_ problem, where the agents own all data samples but part of
their features only (right of Fig. 1) [30, 31]. Such problems appear, for
example, in collaborative learning between multiple parties who hold different
aspects of information for the same set of samples. Using two examples below,
we show that the vertical learning problem can be expressed as the form of
(P).

Example 1 (Empirical risk minimization) Consider a standard empirical risk
minimization problem as follows

$\min_{{\bm{w}}}~{}\sum\limits_{k=1}^{M}\psi({\bm{b}}_{k}^{\top}{\bm{w}},{\bm{v}}_{k})+R({\bm{w}}),$
(3)

where $M$ is the number of data samples, ${\bm{b}}_{k}$ and ${\bm{v}}_{k}$ are
respectively the feature vector and (one-hot vector) label for the $k$-th
sample, ${\bm{w}}\in\mathbb{R}^{nN}$ is the parameter to be optimized, $\psi$
is the (possible non-convex) loss function, and $R(\cdot)$ is a regularization
function.

Let us consider to solve problem (3) over a multi-agent network with $N$
agents under the vertical learning setting. Each agent $i$ owns all $M$ data
samples but has part of the data features only. Specifically, assume that each
sample ${\bm{b}}_{k}$ is partitioned as
${\bm{b}}_{k}=[{\bm{b}}_{1,k}^{\top},\cdots,{\bm{b}}_{N,k}^{\top}]^{\top}$,
where ${\bm{b}}_{i,k}\in\mathbb{R}^{n}$, and each agent $i$ owns samples with
only partial features $\\{{\bm{b}}_{i,k}\\}_{k=1}^{M}$, for all $i\in[N]$. Let
the parameter vector ${\bm{w}}$ be partitioned in the same fashion, i.e.,
${\bm{w}}=[{\bm{w}}_{1}^{\top},\cdots,{\bm{w}}_{N}^{\top}]^{\top}$, and assume
the regularization term $R({\bm{w}})=\sum_{i=1}^{N}R_{i}({\bm{w}}_{i})$. So,
each agent $i$ under the vertical learning setting maintains only a subset of
the model parameters ${\bm{w}}_{i}$.

Then, problem (3) can be reformulated as

$\displaystyle\min_{{{\bm{w}}},{\bm{w}}_{0}}~{}$
$\displaystyle\sum\limits_{k=1}^{M}\psi\big{(}w_{0,k},{\bm{v}}_{k}\big{)}+\sum_{i=1}^{N}R_{i}({\bm{w}}_{i})$
(4a) $\displaystyle{\rm s.t.}~{}$
$\displaystyle\sum_{i=1}^{N}{\bm{b}}_{i,k}^{\top}{\bm{w}}_{i}-w_{0,k}=0,~{}k\in[M],$
(4b)

where ${\bm{w}}_{0}=[w_{0,1},\ldots,w_{0,M}]^{\top}$ is an introduced
variable. Further define
${{\bm{B}}}_{i}=[{\bm{b}}_{i,1},\cdots,{\bm{b}}_{i,M}]^{\top}$, $i\in[N]$, and
let $\Psi({\bm{w}}_{0})=\sum_{k=1}^{M}\psi(w_{0,k},{\bm{v}}_{k})$. Thus (4)
can be written as

$\displaystyle~{}~{}\min_{\begin{subarray}{c}{\bm{w}}_{1},\ldots,{\bm{w}}_{N},\\\
{\bm{w}}_{0}\end{subarray}}~{}$
$\displaystyle\Psi({\bm{w}}_{0})+\sum_{i=1}^{N}R_{i}({\bm{w}}_{i})$ (5a)
$\displaystyle~{}~{}\mathrm{s.t.}~{}$
$\displaystyle\sum_{i=1}^{N}{\bm{B}}_{i}{\bm{w}}_{i}-{\bm{w}}_{0}={\bm{0}}.$
(5b)

One can see that the above problem (5) is indeed an instance of (P). It
corresponds to the vertical learning problem where each agent $i$ knows
partial feature matrix ${\bm{B}}_{i}$, regularizer $R_{i}$ and partial model
parameter ${\bm{w}}_{i}$, for all $i\in[N]$; the variable ${\bm{w}}_{0}$ and
$\Psi$ can be handled by any one of the agents or an extra auxiliary agent.

Figure 2: A neural network structure used in problem (6).

Example 2 (Classification by neural network) Consider to train a neural
network (NN) for classification as shown in Fig. 2. The mapping of the NN is
expressed by the function $F_{{\bm{\theta}}}({\bm{b}}_{k}^{\top}{\bm{W}})$,
where ${\bm{W}}\in\mathbb{R}^{nN\times K}$ denotes the coefficients of the
first linear layer of the NN, $K$ is the dimension of the 1st layer output,
and ${\bm{\theta}}$ incorporates the parameters from the 1st layer output to
the output layer of the NN. The learning problem is given by

$\displaystyle\min_{\begin{subarray}{c}{\bm{W}},{\bm{\theta}}\end{subarray}}~{}\sum\limits_{k=1}^{M}L\big{(}F_{{\bm{\theta}}}({\bm{b}}_{k}^{\top}{\bm{W}}),{\bm{v}}_{k}\big{)},$
(6)

where $L$ is a loss function. By comparing (6) with (3), one can see that the
previous empirical risk minimization problem is equivalent to the NN problem
with only a single linear layer, $K=1$ and no ${\bm{\theta}}$.

Similarly, let us consider the vertical learning setting. Suppose that both
the data samples and model parameters are partitioned as
${\bm{b}}_{k}=[{\bm{b}}_{1,k}^{\top},\cdots,{\bm{b}}_{N,k}^{\top}]^{\top}$ for
all $k\in[M]$, and
${\bm{W}}=[{\bm{W}}_{1}^{\top},\ldots,{\bm{W}}_{N}^{\top}]^{\top}$. By
introducing

$\displaystyle{\bm{w}}_{0,k}=\sum_{i=1}^{N}{\bm{W}}_{i}^{\top}{\bm{b}}_{i,k}=\sum_{i=1}^{N}({\bm{b}}_{i,k}^{\top}\otimes{\bm{I}}_{K}){\bm{w}}_{i},$
(7)

where ${\bm{w}}_{i}={\rm vec}({\bm{W}}_{i})$ is obtained by stacking the
columns of ${\bm{W}}_{i}$, problem (6) can be written as

$\displaystyle\min_{\begin{subarray}{c}{\bm{w}}_{1},\cdots,{\bm{w}}_{N},\\\
({\bm{w}}_{0},{\bm{\theta}})\end{subarray}}~{}$
$\displaystyle\sum\limits_{k=1}^{M}L\big{(}F_{\bm{\theta}}({\bm{w}}_{0,k}^{\top}),{\bm{v}}_{k}\big{)}$
(8a) $\displaystyle{\rm s.t.}~{}$
$\displaystyle\sum_{i=1}^{N}{\bm{B}}_{i}{\bm{w}}_{i}-{\bm{w}}_{0}={\bm{0}},$
(8b)

where
$\textstyle{\bm{B}}_{i}=\textstyle[({\bm{b}}_{i,1}^{\top}\otimes{\bm{I}}_{K})^{\top},\cdots,({\bm{b}}_{i,M}^{\top}\otimes{\bm{I}}_{K})^{\top}]^{\top}$
and
${\bm{w}}_{0}=[{\bm{w}}_{0,1}^{\top},\ldots,{\bm{w}}_{0,M}^{\top}]^{\top}$.

As a result, the distributed optimization problem (P) includes both problems
(5) and (8) as special cases, and thereby any distributed algorithm that can
solve (P) can be applied for solving the above two vertical learning problems.
We should mention that the existing works for solving the vertical learning
problem in the distributed networks consider only convex problems or over a
simple network [30, 31]. In the next section, we present a new distributed
algorithm for non-convex problem (P) over general connected networks.

##  3 Proposed PDC Algorithm

The network models and assumptions are given in Section 3.1. To solve (P), we
develop the PDC algorithm in Section 3.2 and further give its convergence
analysis in Section 3.3. In Section 3.4, we further propose an inexact version
of PDC, referred as the IPDC.

###  3.1 Network Model and Assumptions

The multi-agent network is modeled as an undirected graph
$\mathcal{G}=(\mathcal{E},\mathcal{N})$, where $\mathcal{E}$ is the set of
edges and $\mathcal{N}$ is the set of agents. Each agent can only communicate
and exchange information with its neighbors, i.e., agents $i,j$ can
communicate with each other if and only if $(i,j)\in\mathcal{E}$.

Let us define some useful notations regarding the network graph. The
neighbouring set of agent $i$ is defined as
$\mathcal{N}_{i}=\\{j\in\mathcal{N}~{}|~{}(i,j)\in\mathcal{E}\\}$. Let
$\tilde{}{\bm{D}}={\rm
diag}\\{|\mathcal{N}_{1}|,\cdots,|\mathcal{N}_{N}|\\}\in\mathbb{R}^{N\times
N}$ be the degree matrix of $\mathcal{G}$. The agent-edge incidence matrix is
defined as $\tilde{}{\bm{A}}\in\mathbb{R}^{|\mathcal{E}|\times|\mathcal{N}|}$
which has $\tilde{}{\bm{A}}_{e,i}=1$ and $\tilde{}{\bm{A}}_{e,j}=-1$ for $j>i$
and $e=(i,j)\in\mathcal{E}$, and zeros otherwise.
$\tilde{}{\bm{L}}^{-}=\tilde{}{\bm{A}}^{\top}\tilde{}{\bm{A}}$ is the (signed)
Laplacian matrix of $\mathcal{G}$ . Then
$\tilde{}{\bm{L}}^{+}=2\tilde{}{\bm{D}}-\tilde{}{\bm{L}}^{-}$ is the signless
Laplacian matrix of $\mathcal{G}$. Both $\tilde{}{\bm{L}}^{-}$ and
$\tilde{}{\bm{L}}^{+}$ are positive semi-definite and describe the
connectivity of the network. For ease of later use, we also define

$\displaystyle{\bm{A}}=\tilde{}{\bm{A}}\otimes{\bm{I}}_{M}=[{\bm{A}}_{1},\ldots,{\bm{A}}_{N}],$
$\displaystyle{\bm{L}}^{+}=\tilde{}{\bm{L}}^{+}\otimes{\bm{I}}_{M}=\begin{bmatrix}({\bm{L}}_{1}^{+})^{\top},\cdots,({\bm{L}}_{N}^{+})^{\top}\end{bmatrix}^{\top},$
$\displaystyle{\bm{L}}^{-}=\tilde{}{\bm{L}}^{-}\otimes{\bm{I}}_{M}=\begin{bmatrix}({\bm{L}}_{1}^{-})^{\top},\cdots,({\bm{L}}_{N}^{-})^{\top}\end{bmatrix}^{\top},$

and let $\lambda_{\max}$ be the maximum eigenvalue of ${\bm{L}}^{+}$.

We make the following assumption for the network.

###### Assumption 1

The graph $\mathcal{G}$ is connected.

It means that for any $i,j\in\mathcal{N}$, there exists at least one path
connecting $i$ and $j$. Assumption 1 is commonly made in the literature of
distributed optimization.

For problem (P), we make the following assumption.

###### Assumption 2

Each $f_{i}$ is continuously differentiable and $\nabla f_{i}$ is Lipschitz
continuous with constant $\gamma^{+}>0$.

Under Assumption 2, there exists a (possibly negative) constant $\gamma^{-}$
such that

${\langle}\nabla f_{i}({\bm{x}}_{i})-\nabla
f_{i}({\bm{x}}_{i}^{\prime}),{\bm{x}}_{i}-{\bm{x}}_{i}^{\prime}{\rangle}\geq\gamma^{-}\|{\bm{x}}_{i}-{\bm{x}}_{i}^{\prime}\|^{2},$
(9)

for all ${\bm{x}}_{i},{\bm{x}}_{i}^{\prime}\in\mathbb{R}^{n}$. Note that if
$\gamma^{-}$ is negative, then $f_{i}$ is not a convex function.

###  3.2 Algorithm Development

The proposed PDC method is motivated by the proximal approach recently
presented in [17] for non-convex constrained optimization. In particular, like
[17], we consider a “proximal” counterpart of (P)

$\displaystyle\min_{{\bm{x}}_{i},i\in[N]}~{}$
$\displaystyle\sum_{i=1}^{N}f_{i}({\bm{x}}_{i})+\frac{p}{2}\sum_{i=1}^{N}\|{\bm{x}}_{i}-{\bm{z}}_{i}\|^{2}$
(10a) s.t. $\displaystyle\sum_{i=1}^{N}{\bm{B}}_{i}{\bm{x}}_{i}={\bm{q}},$
(10b)

where $p>-\gamma^{-}$ is a penalty parameter and
${\bm{z}}=({\bm{z}}_{1}^{\top},\ldots,{\bm{z}}_{N}^{\top})^{\top}$ is an
introduced proximal variable. Note that (10) is still not a convex problem for
the joint variable $({\bm{x}}^{\top},{\bm{z}}^{\top})^{\top}$, but it is
convex with respect to ${\bm{x}}$ when ${\bm{z}}$ is fixed. So, the proximal
problem (10) can be regarded as a convex approximation of problem (P) around
the given point ${\bm{z}}$.

Let us denote $({\bm{x}}({\bm{z}}),{\bm{y}}_{0}({\bm{z}}))$ as a primal-dual
optimal solution of (10) with fixed ${\bm{z}}$. One can verify that
$({\bm{x}}({\bm{z}}),{\bm{y}}_{0}({\bm{z}}))$ is a KKT solution of (P), i.e.,
satisfying

$\displaystyle\nabla
f_{i}({\bm{x}}_{i}({\bm{z}}_{i}))+{\bm{B}}_{i}^{\top}{\bm{y}}_{0}({\bm{z}})=\mathbf{0},\forall
i\in[N]$ (11a)
$\displaystyle\sum_{i=1}^{N}{\bm{B}}_{i}{\bm{x}}_{i}({\bm{z}}_{i})-{\bm{q}}=\mathbf{0},$
(11b)

if and only if

$\displaystyle{\bm{x}}({\bm{z}})={\bm{z}}.$ (12)

In [17], the authors proposed to handle (10) by using one step of primal-dual
gradient update in each iteration $r$, followed by updating
$\\{{\bm{z}}^{r}\\}$ by

$\displaystyle{\bm{z}}^{r+1}={\bm{z}}^{r}+\beta({\bm{x}}^{r+1}-{\bm{z}}^{r}),$
(13)

where $0<\beta<1$ is a positive stepsize. Notice that for $\beta=1$, the
algorithm may oscillate. It is shown that (13) with sufficiently small $\beta$
can stabilize the primal-dual sequence and the algorithm can achieve a KKT
solution of (P). However, the algorithm in [17] is a centralized algorithm and
cannot solve (P) over the multi-agent network.

To obtain a distributed algorithm, we consider employing the dual consensus
method in [7] to deal with the (strongly convex) proximal problem (10) when
${\bm{z}}$ is fixed. First, by the Lagrange duality theory [32], (10) is
equivalent to its dual problem when ${\bm{z}}$ is fixed

$\displaystyle\max_{{\bm{y}}_{0}\in\mathbb{R}^{M}}~{}\sum_{i=1}^{N}\bigg{(}\phi_{i}({\bm{y}}_{0},{\bm{z}}_{i})-{\bm{y}}_{0}^{\top}\frac{{\bm{q}}}{N}\bigg{)},$
(14)

where ${\bm{y}}_{0}$ is the dual variable associated with the constraint
(10b), and

$\displaystyle\phi_{i}({\bm{y}}_{0},{\bm{z}}_{i})=\min_{{\bm{x}}_{i}\in\mathbb{R}^{n}}f_{i}({\bm{x}}_{i})+\frac{p}{2}\|{\bm{x}}_{i}-{\bm{z}}_{i}\|^{2}+{\bm{y}}^{\top}_{0}{\bm{B}}_{i}{\bm{x}}_{i},$
(15)

for $i\in[N],$ are concave dual functions of ${\bm{y}}_{0}$. In order to solve
this dual problem in a decentralized manner, we let
${\bm{y}}_{i}\in\mathbb{R}^{M}$ be the local copy of ${\bm{y}}_{0}$ at agent
$i$. Then, $\phi_{i}$ can be written as

$\displaystyle\phi_{i}({\bm{y}}_{i},{\bm{z}}_{i})=\min_{{\bm{x}}_{i}\in\mathbb{R}^{n}}f_{i}({\bm{x}}_{i})+\frac{p}{2}\|{\bm{x}}_{i}-{\bm{z}}_{i}\|^{2}+{\bm{y}}^{\top}_{i}{\bm{B}}_{i}{\bm{x}}_{i},$
(16)

for $i\in[N]$. As the right hand side of (16) is a strongly convex problem
with respect to ${\bm{x}}_{i}$, $\phi_{i}({\bm{y}}_{i},{\bm{z}}_{i})$ is
smooth and concave with respect to ${\bm{y}}_{i}$ according to the Danskin’s
theorem [33], and $\nabla\phi_{i}({\bm{y}}_{i},{\bm{z}}_{i})$ is Lipschitz
continuous with constant $\frac{1}{p+\gamma^{-}}$ [34]. Next, similar to (2),
we consider the following consensus constrained counterpart

$\displaystyle\max_{{\bm{y}}=({\bm{y}}_{1}^{\top},\cdots,{\bm{y}}_{N}^{\top})^{\top}}~{}$
$\displaystyle\sum_{i=1}^{N}\bigg{(}\phi_{i}({\bm{y}}_{i},{\bm{z}}_{i})-{\bm{y}}_{i}^{\top}\frac{{\bm{q}}}{N}\bigg{)}$
(17a) $\displaystyle{\rm s.t.}~{}$ $\displaystyle{\bm{A}}{\bm{y}}={\bm{0}}.$
(17b)

Note that under Assumption 1, ${\bm{A}}{\bm{y}}={\bm{0}}$ is equivalent to
${\bm{y}}_{1}={\bm{y}}_{2}=\ldots={\bm{y}}_{N}$.

We further employ the inexact augmented Lagrangian (AL) method in [35] to
problem (17). In particular, let

$\displaystyle{\mathcal{L}}_{\rho}({\bm{y}},{\bm{\mu}};{\bm{z}})$
$\displaystyle=\sum_{i=1}^{N}\bigg{(}\phi_{i}({\bm{y}}_{i},{\bm{z}}_{i})-{\bm{y}}_{i}^{\top}\frac{{\bm{q}}}{N}\bigg{)}-{\bm{\mu}}^{\top}{\bm{A}}{\bm{y}}-\frac{\rho}{2}\|{\bm{A}}{\bm{y}}\|^{2},$
(18)

be the augmented Lagrangian function of (17), in which
${\bm{\mu}}\in\mathbb{R}^{M|\mathcal{E}|}$ is the dual variable associated
with (17b), and $\rho>0$ is a penalty parameter. Then, according to [35], at
each iteration $r$, we perform

$\displaystyle{\bm{\mu}}^{r+1}={\bm{\mu}}^{r}+\alpha{\bm{A}}{\bm{y}}^{r},$
(19a)
$\displaystyle{\bm{y}}^{r+1}=\arg\max_{{\bm{y}}}{\mathcal{L}}_{\rho}({\bm{y}},{\bm{\mu}}^{r+1};{\bm{z}}^{r})-\frac{\rho}{2}\|{\bm{y}}-{\bm{y}}^{r}\|_{{\bm{L}}^{+}}^{2},$
(19b)

where $\alpha>0$ is a step size parameter. The proximal term
$\frac{\rho}{2}\|{\bm{y}}-{\bm{y}}^{r}\|_{{\bm{L}}^{+}}^{2}$ is used to make
the objective function of (19b) decomposable across
${\bm{y}}_{1},\cdots,{\bm{y}}_{N}$ so as to obtain a fully distributed
algorithm, as we delineate below.

Updates of ${\bm{x}}^{r}$ and ${\bm{y}}^{r}$: According to (16) and (3.2), we
define

$\displaystyle G({\bm{x}},{\bm{y}},{\bm{\mu}};{\bm{z}})$
$\displaystyle\triangleq\sum_{i=1}^{N}\bigg{(}f_{i}({\bm{x}}_{i})+\frac{p}{2}\|{\bm{x}}_{i}-{\bm{z}}_{i}\|^{2}+{\bm{y}}_{i}^{T}{\bm{B}}_{i}{\bm{x}}_{i}\bigg{.}$
.
$\displaystyle~{}~{}~{}-{\bm{y}}_{i}^{\top}\frac{{\bm{q}}}{N}\bigg{)}-{\bm{\mu}}^{\top}{\bm{A}}{\bm{y}}-\frac{\rho}{2}\|{\bm{A}}{\bm{y}}\|^{2}.$
(20)

Then, by (19b), $({\bm{y}}^{r+1},{\bm{x}}^{r+1})$ is the solution of the
following max-min problem:

$\displaystyle\max_{{\bm{y}}}\min_{{\bm{x}}\in\mathbb{R}^{nN}}\bigg{\\{}G({\bm{x}},{\bm{y}},{\bm{\mu}}^{r+1};{\bm{z}}^{r})-\frac{\rho}{2}\|{\bm{y}}-{\bm{y}}^{r}\|^{2}_{{\bm{L}}^{+}}\bigg{\\}}.$
(21)

Similar to [7], based on the saddle point theory, the above max-min problem is
equivalent to the following min-max problem

$\displaystyle\min_{{\bm{x}}\in\mathbb{R}^{nN}}\max_{{\bm{y}}}\bigg{\\{}G({\bm{x}},{\bm{y}},{\bm{\mu}}^{r+1};{\bm{z}}^{r})-\frac{\rho}{2}\|{\bm{y}}-{\bm{y}}^{r}\|^{2}_{{\bm{L}}^{+}}\bigg{\\}}.$
(22)

When fixing ${\bm{x}}$, the inner maximization problem for ${\bm{y}}$ is a
strongly concave quadratic problem. Then, one can show that the inner solution
${\bm{y}}^{r+1}$ are decomposable and has a closed-form expression as

$\displaystyle{\bm{y}}_{i}^{r+1}$
$\displaystyle=\frac{1}{2\rho|\mathcal{N}_{i}|}\big{(}{\bm{B}}_{i}{\bm{x}}_{i}-\frac{{\bm{q}}}{N}-{\bm{A}}^{\top}_{i}{\bm{\mu}}^{r+1}+\rho{\bm{L}}^{+}_{i}{\bm{y}}^{r}\big{)},$
(23)

for all $i\in[N]$. By substituting (23) into (22), the outer solution
${\bm{x}}^{r+1}$ of (22) can be determined by

$\displaystyle{\bm{x}}_{i}^{r+1}=\arg\min_{{\bm{x}}_{i}\in\mathbb{R}^{n}}\bigg{\\{}f_{i}({\bm{x}}_{i})+\frac{p}{2}\|{\bm{x}}_{i}-{\bm{z}}_{i}^{r}\|^{2}$
$\displaystyle~{}~{}~{}+\frac{1}{4\rho|\mathcal{N}_{i}|}\|{\bm{B}}_{i}{\bm{x}}_{i}-{\bm{q}}/N-{\bm{A}}_{i}^{\top}{\bm{\mu}}^{r+1}+\rho{\bm{L}}_{i}^{+}{\bm{y}}^{r}\|^{2}\bigg{\\}},$
(24)

for all $i\in[N]$, which are convex problems given $p>-\gamma^{-}$. Notice
that

$\displaystyle{\bm{L}}_{i}^{+}{\bm{y}}^{r}\
=\sum_{j\in\mathcal{N}_{i}}({\bm{y}}_{i}^{r}+{\bm{y}}_{j}^{r}),$ (25)

and therefore both (23) and (3.2) can be implemented fully distributively
using only local data and messages from the neighbors.

Update of ${\bm{\mu}}^{r}$: Since in (23) and (3.2), one needs
${\bm{A}}_{i}^{\top}{\bm{\mu}}$ instead of ${\bm{\mu}}$. We define
${\bm{p}}_{i}^{r}={\bm{A}}_{i}^{\top}{\bm{\mu}}^{r}$, and according to (19a),
obtain

$\displaystyle{\bm{p}}^{r+1}_{i}$
$\displaystyle={\bm{p}}_{i}^{r}+\alpha{\bm{A}}_{i}^{\top}{\bm{A}}{\bm{y}}^{r}$
$\displaystyle={\bm{p}}_{i}^{r}+\alpha{\bm{L}}_{i}^{-}{\bm{y}}^{r}$
$\displaystyle={\bm{p}}_{i}^{r}+{\alpha}\textstyle\sum_{j\in\mathcal{N}_{i}}({\bm{y}}_{i}^{r}-{\bm{y}}_{j}^{r}),$
(26)

for $i\in[N]$.

By combining (23), (3.2), (3.2), and (13), we obtain the PDC algorithm (see
Algorithm 1). One can see that the algorithm is fully distributed since each
agent $i$ uses only local information and variables
$\\{{\bm{y}}_{j}\\}_{j\in\mathcal{N}_{i}}$ from its neighbors.

###  3.3 Convergence Analysis

We first define the approximate solution of problem (P) as follows.

###### Definition 1

A vector ${\bm{x}}\in\mathbb{R}^{Nn}$ is said to be an $\epsilon$-KKT solution
of problem ${\sf{(P)}}$ if
$\|\sum_{i=1}^{N}{\bm{B}}_{i}{\bm{x}}_{i}-{\bm{q}}\|^{2}\leq\epsilon$ and
there exists ${\bm{y}}\in\mathbb{R}^{M}$ such that $\sum_{i=1}^{N}\|\nabla
f_{i}({\bm{x}}_{i})+{\bm{B}}_{i}^{\top}{\bm{y}}\|^{2}\leq\epsilon$.

The main theoretical results for the proposed PDC method in Algorithm 1 are
given in the following theorem.

###### Theorem 1

Let $p>-\gamma^{-}$, $\rho>0$. Then for sufficiently small $\alpha<\rho$ and
$\beta<1$ (see (62) and (63), the following results hold for Algorithm 1:

  1. (a)

Every limit point of the iteration sequence
$\\{({\bm{z}}^{r},{\bm{y}}^{r})\\}$ is a KKT solution of problem (P);

  2. (b)

There exists a $t<r$ such that ${\bm{z}}^{t}$ is an $\mathcal{O}(1/r)$-KKT
solution of problem (P).

Algorithm 1  Proposed PDC method for solving (P)

1: Given parameters $p,\beta,\alpha,\rho$, and initial variables
${\bm{x}}_{i}^{0}={\bm{z}}_{i}^{0}\in\mathbb{R}^{n}$,
${\bm{y}}_{i}^{0}\in\mathbb{R}^{M}$ and ${\bm{p}}_{i}^{0}={\bm{0}}$ for
$i\in[N]$. Set $r=1.$

2: repeat

3: For all $i\in[N]$ (in parallel), send ${\bm{y}}_{i}^{r}$ to neighbour
$j\in\mathcal{N}_{i}$ and receive ${\bm{y}}_{j}^{r}$ from neighbour
$j\in\mathcal{N}_{i}$. Perform
$\displaystyle{\bm{p}}_{i}^{r+1}=\textstyle~{}{\bm{p}}_{i}^{r}+{\alpha}\sum_{j\in\mathcal{N}_{i}}({\bm{y}}_{i}^{r}-{\bm{y}}_{j}^{r}),$
(27)
$\displaystyle{\bm{x}}_{i}^{r+1}=\arg\min_{{\bm{x}}_{i}}\textstyle\bigg{\\{}f_{i}({\bm{x}}_{i})+\frac{p}{2}\|{\bm{x}}_{i}-{\bm{z}}_{i}^{r}\|^{2}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}+\textstyle\frac{1}{4\rho|\mathcal{N}_{i}|}\\!\big{\|}{\bm{B}}_{i}{\bm{x}}_{i}\\!-\\!\frac{{\bm{q}}}{N}\textstyle\\!-\\!{\bm{p}}_{i}^{r+1}\\!+\\!\rho{\bm{L}}^{+}_{i}{\bm{y}}^{r}\big{\|}^{2}\bigg{\\}},$
(28)
$\displaystyle\textstyle{\bm{y}}_{i}^{r+1}\\!=\\!\frac{1}{2\rho|\mathcal{N}_{i}|}\big{(}{\bm{B}}_{i}{\bm{x}}_{i}^{r+1}\\!-\\!\frac{{\bm{q}}}{N}\\!-\\!{\bm{p}}_{i}^{r+1}\\!+\rho{\bm{L}}^{+}_{i}{\bm{y}}^{r}\big{)},$
(29)
$\displaystyle\textstyle{\bm{z}}_{i}^{r+1}={\bm{z}}_{i}^{r}+\beta({\bm{x}}_{i}^{r+1}-{\bm{z}}_{i}^{r}),$
(30)

4: Set $r\leftarrow r+1.$

5: until a predefined stopping criterion is satisfied.

Proof: The proof relies on analyzing how the gap between ${\bm{z}}^{t}$ and
${\bm{x}}({\bm{z}}^{t})$ progresses with the iterations, and relating it with
the gap between ${\bm{z}}^{t}$ and a KKT solution of (P) (Lemma 9). The
details are given in Section 4. $\blacksquare$

Theorem 1 implies that, given constants $\rho$ and $p>-\gamma^{-}$, if step
size parameters $\alpha$ and $\beta$ are sufficiently small,
$\\{({\bm{z}}^{r},{\bm{y}}^{r})\\}$ can asymptotically converge to the set of
KKT solutions of problem (P). Moreover, Theorem 1(b) shows that an
$\epsilon$-KKT solution can be obtained in $\mathcal{O}({1}/{\epsilon})$
iterations. As mentioned in Section 1, since the finite sum problem is a
special instance of problem (P), the iteration complexity lower bound in [28]
also serves as a lower bound for (P). Therefore, Theorem 1(b) means that the
proposed PDC algorithm can attain the iteration complexity lower bound of
first-order distributed algorithms up to a constant.

While Theorem 1(b) gives only the complexity order, the following corollary
presents the dependency of the bound on the graph topology and data matrix.

###### Corollary 1

Assume that ${\bm{B}}=[{\bm{B}}_{1},\ldots,{\bm{B}}_{N}]$ is full row rank. By
choosing appropriate values for $\alpha,\beta$ and $\rho$, the convergence gap
$\|{\bm{z}}^{t}-{\bm{x}}({\bm{z}}^{t})\|^{2}$ in (4.4) can be further bounded
as

$\displaystyle\|{\bm{z}}^{t}-{\bm{x}}({\bm{z}}^{t})\|^{2}$
$\displaystyle\leq\textstyle\frac{\Phi^{0}-\underline{\Phi}}{r}\textstyle\left(D_{1}B_{\max}^{2}\lambda_{\max}\sqrt{D_{2}\theta_{1}^{4}+D_{3}\theta_{1}^{2}}\right.$
$\displaystyle~{}~{}~{}\left.+\sqrt{D_{4}+\textstyle
B_{\max}^{2}\lambda_{\max}^{2}(D_{5}\theta_{1}^{4}+D_{6}\theta_{1}^{2})\textstyle+\textstyle\frac{D_{7}+D_{8}B_{\max}^{2}\theta_{2}^{2}}{B_{\max}^{2}\lambda_{\max}^{2}(D_{9}\theta_{1}^{4}+D_{10}\theta_{1}^{2})}}\right.$
$\displaystyle~{}~{}~{}\left.\textstyle+\textstyle
D_{11}B_{\max}\theta_{1}\sqrt{D_{12}\theta_{3}^{2}+D_{13}}\times\max\bigg{\\{}D_{14},\sqrt{D_{15}\lambda_{\max}^{2}\theta_{1}^{2}+D_{16}}\bigg{\\}}\right)^{2},$
(31)

where $\theta_{1}$, $\theta_{2}$ and $\theta_{3}$, respectively defined in
Lemma 5, Lemma 6 and Lemma 7, can be bounded as

$\displaystyle\textstyle\theta_{1}=\theta_{2}\leq\textstyle(1+\zeta_{B})^{2}\frac{\sigma_{\max}({\bm{L}}^{-})}{\sigma_{\min}^{2}({\bm{L}}^{-})}+\frac{4(1+\zeta_{B})^{2}}{\sigma_{\max}({\bm{B}}_{\mathrm{diag}})},$
(32)
$\displaystyle\textstyle\theta_{3}\leq\textstyle\frac{\left(2\sqrt{\sigma_{\max}({\bm{L}}^{-})}+\sigma_{\max}({\bm{L}}^{-})\right)(3+\sqrt{\sigma_{\max}({\bm{L}}^{-})})^{2}(1+\zeta_{B})^{2}}{\sigma_{\min}({\bm{L}}^{-})}+\textstyle\frac{4(1+\zeta_{B})^{2}}{\sigma_{\max}({\bm{B}}_{\mathrm{diag}})}\left(3+\sqrt{\sigma_{\max}({\bm{L}}^{-})}\right)^{2},$
(33) $\displaystyle\textstyle
B_{\max}\textstyle\triangleq\textstyle\max_{i\in[N]}\|{\bm{B}}_{i}\|,~{}~{}~{}\textstyle\zeta_{B}\triangleq\textstyle\frac{2\sqrt{N}\cdot\sigma_{\max}({\bm{B}}_{\mathrm{diag}})}{\sigma_{\min}({\bm{B}})},$
(34)

$\lambda_{\max}$ is the largest eigenvalue of ${\bm{L}}^{+}$,
${\bm{B}}_{\mathrm{diag}}$ is a block diagonal matrix with
${\bm{B}}_{1},\ldots,{\bm{B}}_{N}$ being the diagonal blocks, and $D_{1}$ to
$D_{16}$ are constants which depend only on $\gamma^{-}$ and $\gamma^{+}$.

Proof: The bound (1) is obtained by further analyzing the dependency of the
three parameters $\theta_{1}$, $\theta_{2},$ $\theta_{3}$ on the graph
topology and data matrix. The details can be found in Appendix K.
$\blacksquare$

As seen from Corollary 1, both the graph topology (${\bm{L}}^{-}$ and
${\bm{L}}^{+}$) and data matrix ${\bm{B}}$ can have impacts on the algorithm
convergence. However, it is difficult to give a quantitative conclusion due to
its intricacy, except for some simple graphs. The following corollary proved
in [36] shows that the upper bound of (1) can be reduced with the network
connectivity for a cycle graph.

###### Corollary 2

Suppose that the network is a cycle graph. Then the upper bound in (1)
decreases when the network size $N$ is decreased.

Since for the cycle graph, the smaller network size implies better network
connectivity, Corollary 2 also implies the algorithm convergence improves when
the connectivity increases. The proof of Corollary 2 is shown in Appendix O.

###  3.4 Inexact PDC Algorithm

Instead of globally solving the subproblem (3), one may choose to solve it
inexactly by using one-step gradient descent only. Specifically, we replace
(3) by

$\displaystyle{\bm{x}}_{i}^{r+1}$
$\displaystyle={\bm{x}}_{i}^{r}-\zeta\bigg{[}\nabla
f_{i}({\bm{x}}_{i}^{r})+p({\bm{x}}_{i}^{r}-{\bm{z}}_{i}^{r})+\frac{{\bm{B}}_{i}^{\top}}{2\rho|\mathcal{N}_{i}|}\bigg{(}{\bm{B}}_{i}{\bm{x}}_{i}^{r}-\frac{{\bm{q}}}{N}-{\bm{p}}_{i}^{r+1}+\rho{\bm{L}}_{i}^{+}{\bm{y}}^{r}\bigg{)}\bigg{]},$
(35)

where $\zeta>0$ is a stepsize. Then we obtain the inexact PDC (IPDC) as
summarized in Algorithm 2.

The convergence results of Algorithm 2 are shown in the following theorem.

###### Theorem 2

Let $\rho>1/2$, $\beta$ satisfy (63), and $p$, $\zeta$ and $\alpha$ satisfy

$\displaystyle\textstyle
p>\max\bigg{\\{}\gamma^{+}-2\gamma^{-},\textstyle\frac{32(2\rho-1)B_{\max}^{2}}{5\lambda_{\max}^{2}\rho^{2}}-\gamma^{-},\textstyle
32B_{\max}^{2}-\gamma^{-},\textstyle\frac{2\gamma^{+}-7\gamma^{-}}{5}\bigg{\\}},$
$\displaystyle\textstyle\frac{1}{p+\gamma^{-}}<\zeta<\textstyle\min\bigg{\\{}\frac{5}{2(\gamma^{+}-\gamma^{-})},\textstyle\frac{2}{p+\gamma^{+}},\textstyle\frac{1}{32B_{\max}^{2}},\textstyle\bigg{.}\frac{5a_{2}\lambda_{\max}^{2}\rho^{2}}{64(2\rho-1)\theta_{1}^{2}B_{\max}^{2}}\bigg{\\}},$
$\displaystyle\textstyle\alpha<\textstyle\min\bigg{\\{}\frac{\rho}{5},~{}\textstyle\frac{2\rho-1}{16a_{1}\lambda_{\max}^{2}},~{}\textstyle\frac{\rho\zeta(p+\gamma^{-})[1+\rho^{2}(p+\gamma^{+})]}{8B_{\max}^{2}[\zeta(p+\gamma^{-})+3]^{2}}\bigg{\\}},$

where $a_{1}$ and $a_{2}$ are defined in (53) and (54), respectively. Then,
Algorithm 2 has the same convergence properties as Algorithm 1.

Theorem 2 can be proved following similar ways as for Theorem 1; details are
relegated to Appendix P. Theorem 2 shows that the IPDC algorithm enjoys the
same iteration complexity as the PDC algorithm. In practice, the inexact
algorithm may exhibit slower convergence speed than its exact counterpart in
terms of the iteration number, but has a substantially smaller computation
time; see [7, 19].

Algorithm 2  Proposed IPDC method for solving (P)

1: Given parameters $p,\beta,\alpha,\rho$ and stepsize $\zeta$, and initial
variables ${\bm{x}}_{i}^{0}={\bm{z}}_{i}^{0}\in\mathbb{R}^{n}$,
${\bm{y}}_{i}^{0}\in\mathbb{R}^{M}$ and ${\bm{p}}_{i}^{0}={\bm{0}}$ for
$i\in[N]$. Set $r=1.$

2: repeat

3: For all $i\in[N]$ (in parallel), send ${\bm{y}}_{i}^{r}$ to neighbour
$j\in\mathcal{N}_{i}$ and receive ${\bm{y}}_{j}^{r}$ from neighbour $j$.
Perform
$\displaystyle{\bm{p}}_{i}^{r+1}=\textstyle~{}{\bm{p}}_{i}^{r}+{\alpha}\sum_{j\in\mathcal{N}_{i}}({\bm{y}}_{i}^{r}-{\bm{y}}_{j}^{r}),$
(36)
$\displaystyle{\bm{x}}_{i}^{r+1}=\textstyle{\bm{x}}_{i}^{r}-\zeta\bigg{[}\nabla
f_{i}({\bm{x}}_{i}^{r})+p({\bm{x}}_{i}^{r}-{\bm{z}}_{i}^{r})\bigg{.}$
$\displaystyle~{}~{}~{}\bigg{.}\textstyle+\frac{{\bm{B}}_{i}^{\top}}{2\rho|\mathcal{N}_{i}|}\bigg{(}{\bm{B}}_{i}{\bm{x}}_{i}^{r}-\frac{{\bm{q}}}{N}-{\bm{p}}_{i}^{r+1}+\rho{\bm{L}}_{i}^{+}{\bm{y}}^{r}\bigg{)}\bigg{]},$
(37)
$\displaystyle\textstyle{\bm{y}}_{i}^{r+1}\\!=\\!\frac{1}{2\rho|\mathcal{N}_{i}|}\big{(}{\bm{B}}_{i}{\bm{x}}_{i}^{r+1}\\!-\\!\frac{{\bm{q}}}{N}\\!-\\!{\bm{p}}_{i}^{r+1}\\!+\rho{\bm{L}}^{+}_{i}{\bm{y}}^{r}\big{)},$
(38)
$\displaystyle\textstyle{\bm{z}}_{i}^{r+1}=\textstyle{\bm{z}}_{i}^{r}+\beta({\bm{x}}_{i}^{r+1}-{\bm{z}}_{i}^{r}),$
(39)

4: Set $r\leftarrow r+1.$

5: until a predefined stopping criterion is satisfied.

##  4 Proof of Theorem 1

This section presents the proof of Theorem 1. The potential function used in
the proof is firstly developed, followed by showing some key error bounds and
perturbation bounds. Then, we prove the claims in the theorem.

###  4.1 The Potential Function

To prove the convergence, we need to build a proper potential function that
can monotonically descend with the iteration number. Since our proposed
algorithms are primal-dual proximal algorithms, the potential function should
involve both primal and dual information of (17). For ease of later use, we
define

$\displaystyle\mathcal{Y}({\bm{z}})=\arg\max_{{\bm{y}}=[{\bm{y}}_{1}^{\top},\ldots,{\bm{y}}_{N}^{\top}]^{\top}}~{}$
$\displaystyle\sum_{i=1}^{N}\bigg{(}\phi_{i}({\bm{y}}_{i},{\bm{z}}_{i})-{\bm{y}}_{i}^{\top}\frac{{\bm{q}}}{N}\bigg{)}$
$\displaystyle{\rm s.t.}~{}$ $\displaystyle{\bm{A}}{\bm{y}}={\bm{0}}$ (40)

as the set of optimal solutions of (17), which is non-empty since (17) is a
convex problem satisfying the Slater’s condition [37]. By the definition of
$\phi_{i}$ in (16), we define for $i\in[N]$

$\displaystyle\\!\\!\\!{\bm{x}}_{i}({\bm{y}}_{i},{\bm{z}}_{i})=\arg\min\limits_{{\bm{x}}_{i}}f_{i}({\bm{x}}_{i})+\frac{p}{2}\|{\bm{x}}_{i}-{\bm{z}}_{i}\|^{2}+{\bm{y}}_{i}^{\top}{\bm{B}}_{i}{\bm{x}}_{i},$
(41)

and let
${\bm{x}}({\bm{y}},{\bm{z}})=[{\bm{x}}_{1}^{\top}({\bm{y}}_{1},{\bm{z}}_{1}),\cdots,{\bm{x}}_{N}^{\top}({\bm{y}}_{N},{\bm{z}}_{N})]^{\top}$.
Besides, recalling ${\mathcal{L}}_{\rho}$ in (3.2), the set of optimal primal
variables of (17) is given by

$\displaystyle\mathcal{Y}({\bm{\mu}},{\bm{z}})=\arg\max_{{\bm{y}}}{\mathcal{L}}_{\rho}({\bm{y}},{\bm{\mu}};{\bm{z}}),$
(42)

which can be shown to be nonempty given (3.2) (see Appendix A), and the
associated objective value is denoted by $d({\bm{\mu}};{\bm{z}})$.

Next, we show some descent lemmas for $-{\mathcal{L}}_{\rho}$ and $d$.

###### Lemma 1

For the AL function in (3.2), we have

$\displaystyle-{\mathcal{L}}_{\rho}({\bm{y}}^{r+1},{\bm{\mu}}^{r+1};{\bm{z}}^{r+1})+{\mathcal{L}}_{\rho}({\bm{y}}^{r},{\bm{\mu}}^{r};{\bm{z}}^{r})$
$\displaystyle\leq\underbrace{\frac{1}{\alpha}\|{\bm{\mu}}^{r+1}-{\bm{\mu}}^{r}\|^{2}}_{\triangleq\sf(A1)}-\frac{1}{2}\|{\bm{y}}^{r+1}-{\bm{y}}^{r}\|^{2}_{\rho{\bm{L}}^{+}}-\frac{\rho}{2}\|{\bm{y}}^{r+1}-{\bm{y}}^{r}\|^{2}$
$\displaystyle\underbrace{-\frac{p}{2}\sum_{i=1}^{N}{\langle}{\bm{z}}_{i}^{r+1}-{\bm{z}}_{i}^{r},{\bm{z}}_{i}^{r+1}+{\bm{z}}_{i}^{r}-2{\bm{x}}_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r+1}){\rangle}}_{\triangleq\sf(A2)}.$
(43)

###### Lemma 2

For the dual function $d({\bm{\mu}};{\bm{z}})$, we have

$\displaystyle
d({\bm{\mu}}^{r+1};{\bm{z}}^{r+1})-d({\bm{\mu}}^{r};{\bm{z}}^{r})$
$\displaystyle\leq\underbrace{\frac{p}{2}\sum_{i=1}^{N}{\langle}{\bm{z}}_{i}^{r+1}-{\bm{z}}^{r},-2{\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r})+{\bm{z}}_{i}^{r+1}+{\bm{z}}_{i}^{r}{\rangle}}_{\triangleq\sf(B1)}\underbrace{-\alpha{\langle}{\bm{A}}{\bm{y}}^{r},{\bm{A}}\bar{}{\bm{y}}^{r}{\rangle}}_{\triangleq\sf(B2)},$
(44)

where
$\bar{}{\bm{y}}^{r+1}=[(\bar{}{\bm{y}}_{1}^{r+1})^{\top},\ldots,(\bar{}{\bm{y}}_{N}^{r+1})^{\top}]^{\top}\in\mathcal{Y}({\bm{\mu}}^{r+1},{\bm{z}}^{r+1})$
and $\bar{}{\bm{y}}^{r}\in\mathcal{Y}({\bm{\mu}}^{r+1},{\bm{z}}^{r})$ .

We also consider the descent property of $G$ in (3.2).

###### Lemma 3

Define function

$\displaystyle\tilde{G}_{r}({\bm{x}}^{r},{\bm{y}}^{r},{\bm{\mu}}^{r};{\bm{z}}^{r})\triangleq
G({\bm{x}}^{r},{\bm{y}}^{r},{\bm{\mu}}^{r};{\bm{z}}^{r})+\|{\bm{y}}^{r}\|^{2}_{{\bm{L}}^{-}}+\frac{1}{2}\|{\bm{y}}^{r}-{\bm{y}}^{r-1}\|^{2}_{\rho{\bm{L}}^{+}}.$
(45)

Then, we have

$\displaystyle\tilde{G}_{r+1}({\bm{x}}^{r+1},{\bm{y}}^{r+1},{\bm{\mu}}^{r+1};{\bm{z}}^{r+1})-\tilde{G}_{r}({\bm{x}}^{r},{\bm{y}}^{r},{\bm{\mu}}^{r};{\bm{z}}^{r})$
$\displaystyle\leq\underbrace{-\frac{1}{\alpha}\|{\bm{\mu}}^{r+1}-{\bm{\mu}}^{r}\|^{2}}_{\triangleq\sf(C1)}-\sum_{i=1}^{N}\frac{p+\gamma^{-}}{2}\|{\bm{x}}_{i}^{r+1}-{\bm{x}}_{i}^{r}\|^{2}$
$\displaystyle~{}~{}~{}-\frac{1}{2}\|({\bm{y}}^{r}-{\bm{y}}^{r-1})-({\bm{y}}^{r+1}-{\bm{y}}^{r})\|^{2}_{\rho{\bm{L}}^{+}}$
$\displaystyle~{}~{}~{}-\bigg{(}\frac{\rho}{2}-\frac{\alpha}{2}\bigg{)}\|{\bm{y}}^{r+1}-{\bm{y}}^{r}\|^{2}_{{\bm{L}}^{-}}$
$\displaystyle~{}~{}~{}+\|{\bm{y}}^{r+1}-{\bm{y}}^{r}\|^{2}_{\rho{\bm{L}}^{+}}+\frac{p}{2}\bigg{(}1-\frac{2}{\beta}\bigg{)}\|{\bm{z}}^{r+1}-{\bm{z}}^{r}\|^{2}.$
(46)

The above three lemmas can be proved by similar techniques in [17, 35], and
one may refer to Appendix B-D. Based on the above three lemmas, let us define
the potential function as

$\displaystyle\\!\\!\\!\\!\Phi^{r}\\!=\\!\tilde{G}_{r}({\bm{x}}^{r},{\bm{y}}^{r},{\bm{\mu}}^{r};{\bm{z}}^{r})\\!-\\!2{\mathcal{L}}_{\rho}({\bm{y}}^{r},{\bm{\mu}}^{r};{\bm{z}}^{r})\\!+\\!2d({\bm{\mu}}^{r};{\bm{z}}^{r}).$
(47)

By combining Lemmas 1, 2 and 3, we have

$\displaystyle\Phi^{r+1}-\Phi^{r}$ $\displaystyle\leq
2{\sf(A1)}+2{\sf(A2)}+2{\sf(B1)}+2{\sf(B2)}+{\sf(C1)}$
$\displaystyle~{}~{}~{}-\sum_{i=1}^{N}\frac{p+\gamma^{-}}{2}\|{\bm{x}}_{i}^{r+1}-{\bm{x}}_{i}^{r}\|^{2}$
$\displaystyle~{}~{}~{}-\frac{1}{2}\|({\bm{y}}^{r}-{\bm{y}}^{r-1})-({\bm{y}}^{r+1}-{\bm{y}}^{r})\|^{2}_{\rho{\bm{L}}^{+}}$
$\displaystyle~{}~{}~{}-\bigg{(}\frac{\rho}{2}-\frac{\alpha}{2}\bigg{)}\|{\bm{y}}^{r+1}-{\bm{y}}^{r}\|^{2}_{{\bm{L}}^{-}}+\frac{p}{2}\bigg{(}1-\frac{2}{\beta}\bigg{)}\|{\bm{z}}^{r+1}-{\bm{z}}^{r}\|^{2}$
$\displaystyle~{}~{}~{}-{\rho}\|{\bm{y}}^{r+1}-{\bm{y}}^{r}\|^{2}.$ (48)

###  4.2 Perturbation Bounds and Error Bounds

To show that $\Phi^{r}$ is a descent function, we will need some perturbation
bound and error bounds. The following lemma gives the upper bound on
${\bm{x}}({\bm{y}},{\bm{z}})$ when ${\bm{y}}$ and ${\bm{z}}$ are perturbed.

###### Lemma 4

Consider (41) for $i\in[N]$. With $p>-\gamma^{-}$, we have

$\displaystyle\|{\bm{x}}_{i}({\bm{y}}_{i};{\bm{z}}_{i})-{\bm{x}}_{i}({\bm{y}}_{i}^{\prime};{\bm{z}}_{i}^{\prime})\|\leq\sigma_{1}\|{\bm{z}}_{i}-{\bm{z}}_{i}^{\prime}\|+\sigma_{2}\|{\bm{y}}_{i}-{\bm{y}}_{i}^{\prime}\|,$
(49)
$\displaystyle\|{\bm{x}}(\bar{}{\bm{y}}^{r},{\bm{z}}^{r})-{\bm{x}}(\bar{}{\bm{y}}^{r+1},{\bm{z}}^{r+1})\|\leq{\sigma_{3}}\|{\bm{z}}^{r}-{\bm{z}}^{r+1}\|,$
(50)

where $\bar{}{\bm{y}}^{r}\in\mathcal{Y}({\bm{\mu}}^{r+1},{\bm{z}}^{r})$,
$\bar{}{\bm{y}}^{r+1}\in\mathcal{Y}({\bm{\mu}}^{r+1},{\bm{z}}^{r+1})$, and
$\sigma_{1}\textstyle\triangleq\textstyle\frac{p{(p+\gamma^{+}+1)}}{p+\gamma^{-}},$
$\sigma_{2}\textstyle\triangleq\textstyle\frac{B_{\max}{(p+\gamma^{+}+1)}}{p+\gamma^{-}}$
with $B_{\max}\textstyle\triangleq\textstyle\max_{i\in[N]}\|{\bm{B}}_{i}\|$,
and $\sigma_{3}\triangleq\textstyle\frac{p}{p+\gamma^{-}}$.

Proof: See Appendix E. $\blacksquare$

Next, we determine the error bound that characterizes the distance between
${\bm{y}}^{r+1}$ and the solution set $\mathcal{Y}({\bm{\mu}},{\bm{z}})$ of
(42).

###### Lemma 5

With $p>-\gamma^{-}$, there exists a
$\bar{}{\bm{y}}^{r}\in\mathcal{Y}({\bm{\mu}}^{r+1},{\bm{z}}^{r})$ of (42) such
that

$\displaystyle\|{\bm{y}}^{r+1}-\bar{}{\bm{y}}^{r}\|_{{\bm{L}}^{-}}^{2}\leq
a_{1}\|{\bm{L}}^{+}({\bm{y}}^{r+1}-{\bm{y}}^{r})\|^{2},$ (51)
$\displaystyle\|{\bm{y}}^{r+1}-\bar{}{\bm{y}}^{r}\|^{2}\leq
a_{2}\|{\bm{L}}^{+}({\bm{y}}^{r+1}-{\bm{y}}^{r})\|^{2},$ (52)

where

$\displaystyle a_{1}$
$\displaystyle\triangleq\frac{\frac{\theta_{1}^{2}}{2\rho}\bigg{(}\frac{1}{p+\gamma^{-}}+\rho^{2}(p+\gamma^{+})\bigg{)}^{2}+\rho}{\frac{1}{p+\gamma^{-}}+\rho^{2}(p+\gamma^{+})},$
(53) $\displaystyle a_{2}$
$\displaystyle\triangleq\frac{\theta_{1}^{4}}{\rho^{2}}\bigg{(}\frac{1}{p+\gamma^{-}}+\rho^{2}(p+\gamma^{+})\bigg{)}^{2}+2\theta_{1}^{2},$
(54)

and $\theta_{1}>0$ is a constant only depending on the problem instance.

Proof: See Appendix F. $\blacksquare$

The following lemma describes the change of
$\bar{}{\bm{y}}\in\mathcal{Y}({\bm{\mu}},{\bm{z}})$ when ${\bm{z}}$ is
perturbed.

###### Lemma 6

Consider (42) with $p>-\gamma^{-}$. For any
$\bar{}{\bm{y}}^{r}\in\mathcal{Y}({\bm{\mu}}^{r+1},{\bm{z}}^{r})$, there
exists a $\bar{}{\bm{y}}^{r+1}\in\mathcal{Y}({\bm{\mu}}^{r+1},{\bm{z}}^{r+1})$
such that

$\displaystyle\|\bar{}{\bm{y}}^{r}-\bar{}{\bm{y}}^{r+1}\|^{2}\leq{a_{3}}\|{\bm{z}}^{r}-{\bm{z}}^{r+1}\|^{2}.$
(55)

where
$a_{3}\triangleq\theta_{2}^{2}\bigg{(}(p+\gamma^{+})^{2}{\sigma_{3}}^{2}+\frac{B_{\max}^{2}{\sigma_{3}}^{2}}{\rho^{2}}\bigg{)}$
and $\theta_{2}>0$ only depends on the problem instance.

Proof: See Appendix G. $\blacksquare$

From above (52) and (55), we can obtain that for ${\bm{y}}^{r+1}$, there
exists a $\bar{}{\bm{y}}^{r+1}\in\mathcal{Y}({\bm{\mu}}^{r+1},{\bm{z}}^{r+1})$
so that

$\displaystyle\|{\bm{y}}^{r+1}-\bar{}{\bm{y}}^{r+1}\|^{2}$ $\displaystyle\leq
2a_{2}\|{\bm{L}}^{+}({\bm{y}}^{r+1}-{\bm{y}}^{r})\|^{2}+2a_{3}\|{\bm{z}}^{r+1}-{\bm{z}}^{r}\|^{2}.$
(56)

Next, we consider the error bound between (40) and (42) due to the violation
of the consensus constraint.

###### Lemma 7

For any $\bar{}{\bm{y}}^{r}\in\mathcal{Y}({\bm{\mu}}^{r+1},{\bm{z}}^{r})$,
there exists a ${\bm{y}}({\bm{z}}^{r})\in\mathcal{Y}({\bm{z}}^{r})$ such that

$\displaystyle\|{{\bm{y}}}({\bm{z}}^{r})-\bar{}{\bm{y}}^{r}\|\leq\sqrt{a_{4}}\|{\bm{A}}\bar{}{\bm{y}}^{r}\|,$
(57)

where
$a_{4}\textstyle\triangleq{\textstyle\frac{1}{2}\bigg{(}\frac{\theta_{3}^{2}}{p+\gamma^{-}}+\rho^{2}\theta_{3}^{2}(p+\gamma^{+})\bigg{)}^{2}+\frac{\theta_{3}^{2}}{\rho^{2}}}$
for some $\theta_{3}>0$ only depending on the problem instance.

Proof: See Appendix H. $\blacksquare$

According to Lemma 4, Lemma 5, and Lemma 7, we can bound the difference
between $({\bm{x}}^{r+1},{\bm{y}}^{r+1})$ generated by Algorithm 1 and the
primal-dual solution $({\bm{x}}({\bm{z}}^{r}),{\bm{y}}({\bm{z}}^{r}))$ of
problem (10) as follows.

###### Corollary 3

For ${\bm{y}}^{r+1}$, there exists a
$\bar{}{\bm{y}}^{r}\in\mathcal{Y}({\bm{\mu}}^{r+1},{\bm{z}}^{r})$ and
${\bm{y}}({\bm{z}}^{r})\in\mathcal{Y}({\bm{z}}^{r})$ such that

$\displaystyle\|{\bm{y}}^{r+1}-{\bm{y}}({\bm{z}}^{r})\|\\!\leq\\!\sqrt{a_{2}}\lambda_{\max}\|{\bm{y}}^{r}\\!-\\!{\bm{y}}^{r+1}\|\\!+\\!\sqrt{a_{4}}\|{\bm{A}}\bar{}{\bm{y}}^{r}\|,$
(58)
$\displaystyle\|{\bm{x}}^{r+1}-{\bm{x}}({\bm{z}}^{r})\|\leq\sigma_{2}\|{\bm{y}}^{r+1}-{\bm{y}}({\bm{z}}^{r})\|.$
(59)

Proof: To show (58), we have

$\displaystyle\|{\bm{y}}^{r+1}-{\bm{y}}({\bm{z}}^{r})\|$
$\displaystyle=\|{\bm{y}}^{r+1}-\bar{}{\bm{y}}^{r}+\bar{}{\bm{y}}^{r}-{\bm{y}}({\bm{z}}^{r})\|$
$\displaystyle\leq\|{\bm{y}}^{r+1}-\bar{}{\bm{y}}^{r}\|+\|\bar{}{\bm{y}}^{r}-{\bm{y}}({\bm{z}}^{r})\|.$
(60)

By substituting (52) in Lemma 5 and (57) in Lemma 7 into (P.3), one obtains
(58). To prove (59), we first note from (21), (3.2) and (41) that
${\bm{x}}^{r+1}={\bm{x}}({\bm{y}}^{r+1},{\bm{z}}^{r})$. Thus, by applying
(49), we obtain

$\displaystyle\|{\bm{x}}^{r+1}-{\bm{x}}({\bm{z}}^{r})\|$
$\displaystyle=\|{\bm{x}}({\bm{y}}^{r+1},{\bm{z}}^{r})-{\bm{x}}({\bm{y}}({\bm{z}}^{r}),{\bm{z}}^{r})\|$
$\displaystyle\leq\sigma_{2}\|{\bm{y}}^{r+1}-{\bm{y}}({\bm{z}}^{r})\|,$ (61)

which is (59). $\blacksquare$

With the above results, we are ready to prove Theorem 1.

###  4.3 Proof of Theorem 1(a)

We first show that the potential function $\Phi^{r}$ in (47) is non-increasing
and is lower bounded.

###### Lemma 8

Let $p>-\gamma^{-}$, $\rho>0$ and

$\displaystyle\alpha\leq\min\bigg{\\{}\frac{\rho}{5},\frac{\rho}{8a_{1}\lambda_{\max}^{2}}\bigg{\\}},$
(62)
$\displaystyle\beta<\bigg{(}\frac{1}{2}+\frac{20p\sigma_{2}^{2}a_{2}\lambda_{\max}^{2}}{\rho}+\frac{\rho(\sigma_{1}^{2}+2\sigma_{2}^{2}a_{3})}{10p\sigma_{2}^{2}a_{2}\lambda_{\max}^{2}}\bigg{)}^{-1}.$
(63)

(a) Then, we have

$\displaystyle\Phi^{r+1}-\Phi^{r}$
$\displaystyle\leq-\sum_{i=1}^{N}\frac{p+\gamma^{-}}{2}\|{\bm{x}}_{i}^{r+1}-{\bm{x}}_{i}^{r}\|^{2}-C_{1}\|{\bm{y}}^{r+1}-{\bm{y}}^{r}\|^{2}_{{\bm{L}}^{-}}$
$\displaystyle~{}~{}~{}-\frac{1}{2}\|({\bm{y}}^{r}-{\bm{y}}^{r-1})-({\bm{y}}^{r+1}-{\bm{y}}^{r})\|^{2}_{{\bm{L}}^{+}}-\alpha\|{\bm{A}}\bar{}{\bm{y}}^{r}\|^{2}$
$\displaystyle~{}~{}~{}-C_{2}\|{\bm{y}}^{r+1}-{\bm{y}}^{r}\|^{2}-C_{3}\|{\bm{z}}^{r+1}-{\bm{z}}^{r}\|^{2}.$
(64)

where

$\displaystyle\textstyle C_{1}$
$\displaystyle\textstyle\triangleq\frac{\rho-5\alpha}{2}\geq
0,~{}C_{2}\triangleq\rho-\bigg{(}2\alpha
a_{1}+\frac{{\rho}}{{5\lambda_{\max}^{2}}}\bigg{)}\lambda_{\max}^{2}>0,$
$\displaystyle\textstyle C_{3}$ $\displaystyle\textstyle\triangleq
p\bigg{(}-\frac{1}{2}+\frac{1}{\beta}-\frac{20p\sigma_{2}^{2}a_{2}\lambda_{\max}^{2}}{\rho}-\frac{\rho(\sigma_{1}^{2}+2\sigma_{2}^{2}a_{3})}{10p\sigma_{2}^{2}a_{2}\lambda_{\max}^{2}}\bigg{)}>0.$

(b) Moreover, $\Phi^{r}$ is lower bounded, i.e.,
$\sum_{r=0}^{T}\Phi^{r+1}>\underline{\Phi}$ for some constant
$\underline{\Phi}>-\infty$.

Proof: Firstly, consider ${\sf(A1)}$ in (1), ${\sf(B2)}$ in (2), $\sf{(C1)}$
in (3), and their combination:

$\displaystyle 2{\sf(A1)}+2{\sf(B2)}+{\sf(C1)}$
$\displaystyle=\frac{1}{\alpha}\|{\bm{\mu}}^{r+1}-{\bm{\mu}}^{r}\|^{2}-2\alpha{\langle}{\bm{A}}{\bm{y}}^{r},{\bm{A}}\bar{}{\bm{y}}^{r}{\rangle}$
$\displaystyle\overset{\mathrm{(i)}}{=}\alpha\|{\bm{A}}{\bm{y}}^{r}\|^{2}-2\alpha{\langle}{\bm{A}}{\bm{y}}^{r},{\bm{A}}\bar{}{\bm{y}}^{r}{\rangle}$
$\displaystyle\overset{\mathrm{(ii)}}{\leq}2\alpha\|{\bm{y}}^{r}-{\bm{y}}^{r+1}\|^{2}_{{\bm{L}}^{-}}+2\alpha\|{\bm{y}}^{r+1}-\bar{}{\bm{y}}^{r}\|_{{\bm{L}}^{-}}^{2}-\alpha\|{\bm{A}}\bar{}{\bm{y}}^{r}\|^{2}$
$\displaystyle\overset{\mathrm{(iii)}}{\leq}2\alpha\|{\bm{y}}^{r}-{\bm{y}}^{r+1}\|^{2}_{{\bm{L}}^{-}}+2\alpha
a_{1}\|{\bm{L}}^{+}({\bm{y}}^{r+1}-{\bm{y}}^{r})\|^{2}-\alpha\|{\bm{A}}\bar{}{\bm{y}}^{r}\|^{2},$
(65)

where $\mathrm{(i)}$ is due to (19a), $\mathrm{(ii)}$ is obtained by adding
and subtracting a term ${\bm{y}}^{r+1}$ followed by using inequality
$\|{\bm{a}}+{\bm{b}}\|^{2}\leq 2\|{\bm{a}}\|^{2}+2\|{\bm{b}}\|^{2}$, and
lastly $\mathrm{(iii)}$ is owing to (51) of Lemma 5.

Secondly, consider

$\displaystyle{\sf 2(A2)}+2{\sf(B1)}$
$\displaystyle=2p\sum_{i=1}^{N}{\langle}{\bm{z}}_{i}^{r+1}-{\bm{z}}_{i}^{r},{\bm{x}}_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r+1})-{\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r}){\rangle}$
$\displaystyle\leq
p\delta\|{\bm{z}}^{r+1}-{\bm{z}}^{r}\|^{2}+\frac{p}{\delta}\sum_{i=1}^{N}\|{\bm{x}}_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r+1})-{\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r})\|^{2},$
(66)

where the last inequality is obtained by Young’s inequality and $\delta>0$. By
further applying (49) of Lemma 4 to (66), we obtain

$\displaystyle{\sf 2(A2)}+2{\sf(B1)}$ $\displaystyle\leq
p\delta\|{\bm{z}}^{r+1}-{\bm{z}}^{r}\|^{2}+\sum_{i=1}^{N}\frac{2p\sigma_{1}^{2}}{\delta}\|{\bm{z}}_{i}^{r+1}-{\bm{z}}_{i}^{r}\|^{2}+\sum_{i=1}^{N}\frac{2p\sigma_{2}^{2}}{\delta}\|{\bm{y}}_{i}^{r+1}-\bar{}{\bm{y}}_{i}^{r+1}\|^{2}$
$\displaystyle\overset{\mathrm{(i)}}{\leq}\bigg{(}p\delta+\frac{2p\sigma_{1}^{2}}{\delta}\bigg{)}\|{\bm{z}}^{r+1}-{\bm{z}}^{r}\|^{2}+\frac{2p\sigma_{2}^{2}}{\delta}\bigg{[}2a_{2}\|{\bm{L}}^{+}({\bm{y}}^{r+1}-{\bm{y}}^{r})\|^{2}+2a_{3}\|{\bm{z}}^{r+1}-{\bm{z}}^{r}\|^{2}\bigg{]}$
$\displaystyle=\bigg{(}p\delta+\frac{2p\sigma_{1}^{2}}{\delta}+\frac{4p\sigma_{2}^{2}a_{3}}{\delta}\bigg{)}\|{\bm{z}}^{r+1}-{\bm{z}}^{r}\|^{2}+\frac{4p\sigma_{2}^{2}a_{2}}{\delta}\|{\bm{L}}^{+}({\bm{y}}^{r+1}-{\bm{y}}^{r})\|^{2},$
(67)

where $\mathrm{(i)}$ is obtained by (56).

By substituting (4.3) and (67) into (4.1), we thereby obtain

$\displaystyle\Phi^{r+1}-\Phi^{r}$
$\displaystyle=-\sum_{i=1}^{N}\frac{p+\gamma^{-}}{2}\|{\bm{x}}_{i}^{r+1}-{\bm{x}}_{i}^{r}\|^{2}-\underbrace{\frac{\rho-5\alpha}{2}}_{\triangleq\rm{C_{1}}}\|{\bm{y}}^{r}-{\bm{y}}^{r+1}\|^{2}_{{\bm{L}}^{-}}$
$\displaystyle~{}~{}~{}-\frac{1}{2}\|({\bm{y}}^{r}-{\bm{y}}^{r-1})-({\bm{y}}^{r+1}-{\bm{y}}^{r})\|^{2}_{\rho{\bm{L}}^{+}}-\alpha\|{\bm{A}}\bar{}{\bm{y}}^{r}\|^{2}$
$\displaystyle~{}~{}~{}-\underbrace{[\rho-(2\alpha
a_{1}+\frac{4p\sigma_{2}^{2}a_{2}}{\delta})\lambda_{\max}^{2}]}_{\triangleq\rm{C_{2}}}\|{\bm{y}}^{r+1}-{\bm{y}}^{r}\|^{2}$
$\displaystyle~{}~{}~{}-\underbrace{p(-\frac{1}{2}+\frac{1}{\beta}-\delta-\frac{2\sigma_{1}^{2}}{\delta}-\frac{4\sigma_{2}^{2}a_{3}}{\delta})}_{\triangleq\rm{C_{3}}}\|{\bm{z}}^{r+1}-{\bm{z}}^{r}\|^{2}.$
(68)

We see that if $\alpha\leq\rho/5$, then $C_{1}$ is nonnegative. In addition,
choose
$\delta=\frac{{20}p\sigma_{2}^{2}a_{2}\lambda_{\max}^{2}}{\rho}>\frac{16p\sigma_{2}^{2}a_{2}\lambda_{\max}^{2}}{\rho}$,
i.e., $\frac{4p\sigma_{2}^{2}a_{2}\lambda_{\max}^{2}}{\delta}<\frac{\rho}{4}$,
then as long as $\alpha\leq\frac{\rho}{8a_{1}\lambda_{\max}^{2}},$ we can have
$C_{2}\geq{\rho/2}>0$. Lastly, if $\beta$ satisfies (63), then $C_{3}$ is
positive. Thus, Lemma 8(a) is proved.

The proof for showing $\Phi^{r}$ is lower bounded (Lemma 8(b)) is referred to
Appendix I. $\blacksquare$

According to Lemma 8, we have

$\displaystyle\|{\bm{x}}^{r}-{\bm{x}}^{r+1}\|\rightarrow
0,~{}~{}\|{\bm{y}}^{r}-{\bm{y}}^{r+1}\|\rightarrow 0,$ (69)
$\displaystyle\|{{\bm{A}}}\bar{}{\bm{y}}^{r}\|\rightarrow
0,~{}~{}\|{\bm{z}}^{r}-{\bm{z}}^{r+1}\|\rightarrow 0.$ (70)

By (13) and $\|{\bm{z}}^{r}-{\bm{z}}^{r+1}\|\rightarrow 0$, one can have
$\|{\bm{x}}^{r+1}-{\bm{z}}^{r+1}\|\rightarrow 0$. Besides, combing (52) and
$\|{{\bm{A}}}\bar{}{\bm{y}}^{r}\|\rightarrow 0$, one then obtains
$\|{{\bm{A}}}{\bm{y}}^{r+1}\|\rightarrow 0$. Therefore, we obtain

$\displaystyle\|{\bm{x}}^{r+1}-{\bm{z}}^{r+1}\|\rightarrow
0,~{}~{}\|{{\bm{A}}}{\bm{y}}^{r+1}\|\rightarrow 0.$ (71)

By applying (69) and (71) to the KKT conditions of (19b), one then concludes
that every limit point of $\\{({\bm{z}}^{r},{\bm{y}}^{r})\\}$ is a KKT
solution of (P) $\blacksquare$

###  4.4 Proof of Theorem 1(b)

Since $\Phi$ is non-increasing and lower bounded by Lemma 8, one can show that
for $r>0$, there exists a $t\leq r$ such that

$\displaystyle\textstyle\Phi^{t}-\Phi^{t+1}\leq\frac{\Phi^{0}-\underline{\Phi}}{r}.$
(72)

Then, by (13), (58) and (59), we can bound

$\displaystyle\|{\bm{z}}^{t}-{\bm{x}}({\bm{z}}^{t})\|\leq\|{\bm{z}}^{t}-{\bm{x}}^{t+1}\|+\|{\bm{x}}^{t+1}-{\bm{x}}({\bm{z}}^{t})\|$
$\displaystyle\leq\frac{1}{\beta}\|{\bm{z}}^{t}-{\bm{z}}^{t+1}\|+\sigma_{2}\big{(}\sqrt{a_{2}}\lambda_{\max}\|{\bm{y}}^{t}\\!-\\!{\bm{y}}^{t+1}\|\\!+\sqrt{a_{4}}\|{\bm{A}}\bar{}{\bm{y}}^{t}\|\big{)}$
$\displaystyle\overset{\mathrm{(i)}}{\leq}{\textstyle\sqrt{\frac{\Phi^{0}-\underline{\Phi}}{r}}}\bigg{(}\frac{1}{\beta\sqrt{C_{3}}}+\frac{\sigma_{2}\sqrt{a_{2}}\lambda_{\max}}{\sqrt{C_{2}}}+\frac{\sigma_{2}\sqrt{a_{4}}}{\sqrt{\alpha}}\bigg{)}$
$\displaystyle\triangleq\tilde{C}\sqrt{\frac{1}{r}},$ (73)

where $\mathrm{(i)}$ is obtained by applying (64).

Next, we connect (4.4) with an KKT solution of problem (P) by the following
lemma, which is proved in Appendix J.

###### Lemma 9

Suppose that $\|{\bm{z}}-{\bm{x}}({\bm{z}})\|^{2}\leq\epsilon$ for some
${\bm{z}}$. Then, ${\bm{z}}$ is a $\kappa\epsilon$-KKT solution of problem
(P), where $\kappa=\max\\{2(p^{2}+(\gamma^{+})^{2},NB_{\max}^{2}\\}$ and
$B_{\max}\triangleq\max_{i\in[N]}\|{\bm{B}}_{i}\|$.

By the above lemma, we conclude that ${\bm{z}}^{t}$ satisfying (4.4) is a
$\big{(}\kappa\tilde{C}^{2}/r\big{)}$-KKT solution of problem (P).
$\blacksquare$

##  5 Numerical Results

In this section, we present two numerical examples to illustrate the empirical
performance of the proposed algorithms.

###  5.1 Distributed Logistic Regression

Following (3), we formulate a non-convex regularized logistic regression
problem [38] with

$\displaystyle\psi({\bm{b}}_{k}^{\top}{\bm{w}};v_{k})=\log{\left(1+e^{-v_{k}({\bm{b}}_{k}^{\top}{\bm{w}})}\right)},~{}\forall
k\in[M],$ (74a) $\displaystyle
R_{i}({\bm{w}}_{i})=\lambda\sum_{s=1}^{n}\frac{\xi w_{i,s}^{2}}{1+\xi
w_{i,s}^{2}},~{}\forall i\in[N],$ (74b)

where $v_{k}\in\\{\pm 1\\}$ are binary labels, and $\lambda,\xi>0$ are the
parameters.

(a) $\beta=0.1$, $\rho=0.01$, $p=0.01$, and various values of $\alpha$

(b) $\beta=0.1$, $p=0.01$, $\alpha=0.01$, and various values of $\rho$

Figure 3: Convergence performance of the PDC algorithm with various values of
$\alpha$ and $\rho$.

We consider the two images D24 and D68 of Brodatz data
set444http://www.ux.uis.no/tranden/broadatz.html, and extract $M/2$
overlapping patches with dimension $\sqrt{Nn}\times\sqrt{Nn}$ from these two
images, respectively. Each patch is vectorized into $1\times Nn$-vector, and
they are randomly shuffled and stacked into the data matrix ${\bm{B}}$. Here
we set $M=100$ and $n=100$, and assumed a network with $N=25$ nodes. The
generation of the network graph follows the same method as in [39]. For (74),
we set $\lambda=0.01$ and $\xi=0.5$.

(a) $\beta=0.1$, $\rho=0.01$, $\alpha=0.01$, and various values of $p$

(b) $p=0.01$, $\rho=0.01$, $\alpha=0.01$, and various values of $\beta$

Figure 4: Convergence performance of the PDC algorithm with various values of
$p$ and $\beta$.

We employ the proposed PDC algorithm (Algorithm 1) to handle the non-convex LR
problem in the form of (5). Our aim is to use this example to examine how the
algorithm parameters $\alpha$, $p$, $\beta$, and $\rho$ can affect the
convergence of proposed PDC algorithms. The fast iterative shrinkage
thresholding algorithm (FISTA) [40] is used to solve subproblem (3), and the
stopping condition of FISTA is that the normalized proximal gradient [40] is
smaller than $10^{-5}$. The entries of the initial variables ${\bm{x}}^{0}$
and ${\bm{y}}^{0}$ are randomly generated following a uniform distribution
over the interval $[-1,1]$. The experiments are performed 10 times, each with
a different initial (${\bm{x}}^{0}$, ${\bm{y}}^{0}$), and the averaged results
are presented. In particular, we assess the algorithm performance by
calculating the following two terms

$\displaystyle{\rm
Graident~{}residue}=\textstyle\frac{1}{Nn}\sum_{i=1}^{N}\|\nabla
f_{i}({\bm{x}}_{i}^{r})+{\bm{B}}_{i}^{\top}{\bm{y}}_{i}^{r}\|^{2},$ (75)
$\displaystyle{\rm
Infeasibility}=\textstyle\frac{1}{M}\|\sum_{i=1}^{N}{\bm{B}}_{i}{\bm{x}}_{i}^{r}-{\bm{q}}\|^{2},$
(76)

with respect to the iteration number (communication round) $r$. The simulation
results are displayed in Fig. 3, where we respectively examine the impacts of
the four parameters $\alpha$, $\rho$, $p$ and $\beta$ on the algorithm
convergence.

In Fig. 3(a), we set $\beta=0.1$, $\rho=0.01$, $p=0.01$, and
$\alpha\in\\{10^{-1},10^{-2},10^{-3},10^{-4}\\}$. One can see that a smaller
$\alpha$ may slow down the proposed PDC algorithm in terms of the gradient
residue (left figure), whereas slightly reducing the feasibility level (right
figure). It is observed that when $\alpha=10^{-1}>\rho$, the algorithm cannot
satisfy the linear constraint even though it achieves the smallest gradient
residue. Such observation is consistent with Theorem 1.

In Fig. 3(b), we set $\beta=0.1$, $p=0.01$, $\alpha=0.01$, and
$\rho\in\\{10^{-2},10^{-1},1,10,10^{2}\\}$. One can see from the two figures
that parameter $\rho$ does not significantly impact the convergence of
gradient residue, but can greatly affect the constraint feasibility. In
particular, a larger $\rho$ can cause larger infeasibility values.

In Fig. 4(a), we set $\beta=0.1$, $\rho=0.01$, $\alpha=0.01$, and
$p\in\\{10^{-2},10^{-1},1,10,10^{2}\\}$. One can see from the two figures that
parameter $p$ does not significantly impact the constraint feasibility,
whereas a larger $p$ can slow down the convergence of gradient residue. This
is expected since with a larger $p$, problem (10) is a more conservative
convex approximation, and thus it slows down the algorithm convergence.

Lastly, in Fig. 4(b), we set $p=0.01$, $\rho=0.01$, $\alpha=0.01$, and
$\beta\in\\{10^{-3},10^{-2},10^{-1},1\\}$. One can see that the step size
$\beta$ has a mild influence on the convergence of both gradient residue and
infeasibility. Interestingly, while Theorem 1 suggests that $\beta$ should be
less than one and a small number, it is found that the algorithm behaves well
with $\beta=1$.

In Appendix Q, we further present the convergence results for the IPDC
algorithm with respect to the algorithm parameters $\alpha$, $\rho$, $p$,
$\beta$ and $\zeta$.

###  5.2 Distributed Neural Network

(a) Training loss

(b) Testing accuracy

Figure 5: Convergence performance of the IPDC algorithm with $p=10$,
$\beta=0.01$, $\alpha=0.01$, and various values of $\rho$ and $\zeta$. Left
column: $\zeta=0.1$; Right column: $\rho=10$.

In this simulation, following (6), we consider the use of the NN to classify
handwritten digits in the MNIST dataset. The $60000$ training images are
divided into $12$ batches. The dimension of each image is $28\times 28$ which
is vectorized into a $1\times{784}$ vector ($nN=784$). These vectors of the
5000 samples are stacked as the data matrix ${\bm{B}}$. Then this data matrix
is column partitioned into $8$ equal parts to be distributed to $8$ agents
($N=8$). The network topology is generated in the same way as in [39].

This NN has 2 layers, with $K=30$ neurons in the hidden layer. The activation
function of the hidden layer and output layer are Rectified Linear Unit (ReLU)
and softmax function, respectively. The cross-entropy loss function [41] is
applied in the last layer. Then, the proposed IPDC algorithm are applied to
train the classification NN by solving problem (8). The experiment is
performed on the PyTorch platform.

For the IPDC algorithms, it is set that $p=10$, $\alpha=0.01$, $\beta=0.01$,
$\rho\in\\{10^{-1},1,10,100\\}$ and $\zeta\in\\{10^{-3},10^{-2},10^{-1}\\}$.
Moreover, the simulations are preformed with $3$ randomly generated initial
points, and their average results are presented in Fig. 5.

Fig. 5 displays the training loss and testing accuracy of the NN versus the
iteration number. From the left column of Fig. 5, one can observe that the
IPDC algorithm with a larger $\rho$ has a sharper decrease in the training
loss and increase in the testing accuracy in the first few iterations, but the
algorithm can eventually achieve similar training loss and testing accuracy
for the tested values of $\rho$. In the right column of Fig. 5, one can
observe that the step size $\zeta$ can influence the performance of the
proposed IPDC algorithm significantly. Though it is shown in Theorem 2 that a
sufficiently small $\zeta$ is required for the algorithm convergence of the
IPDC algorithm, it is found from the figure that with a larger $\zeta$ the
IPDC algorithm converges faster in both training loss and testing accuracy.
Lastly, it is interesting to see that, while our theoretical claims in Theorem
2 are under the smoothness assumption (Assumption 2), the simulation results
above demonstrate that the proposed algorithm can still work well for the non-
smooth ReLU neural network.

##  6 Conclusions

In this paper, two new distributed algorithms, i.e., the PDC and the IPDC
algorithms, for solving the non-convex linearly constrained problem ${\sf(P)}$
over the multi-agent network have been proposed. These algorithms can be
regarded as extensions of the dual consensus ADMM method in [7] to the
challenging non-convex settings, by leveraging the proximal technique in [17].
While the algorithm development is rather straightforward, it turns out that
proving the convergence and convergence rate of the proposed algorithms are
not. By developing some key error bounds and perturbation bounds, we have
shown that, under mild conditions on the algorithm parameters, the proposed
algorithms both can converge to a KKT solution of problem ${\sf(P)}$ with a
complexity order of $\mathcal{O}(1/\epsilon)$. Through experiments on a non-
convex logistic regression problem and an NN-based MNIST classification task,
we have demonstrated the good convergence behaviors of the proposed algorithms
with respect to various parameter settings.

To the best of our knowledge, the presented distributed algorithms and
convergence analyses are the first for the linearly constrained non-convex
problem ${\sf(P)}$. In the future, it is worthwhile to extend the current PDC
framework to mini-batch stochastic gradient descent (SGD) [42, 43], local SGD
[44], online learning [45], and to handle more complex learning problems with
heterogeneous data distribution [46], and hybrid partitioned data samples and
features [47]. Inspired by the simulation results in Section 5.2, it is also
meaningful to analyze convergence properties of the proposed algorithm for
ReLU neural networks [48, 49] or devise distributed algorithms for general
non-smooth problem ${\sf(P)}$. Finally, as we discuss after Corollary 1, it
needs more efforts to understand how graph topology and data matrix can affect
the algorithm convergence. Besides, it is also important to devise optimal
distributed schemes as in [28] for the general problem (P).

##  Appendix A Proof of nonempty property of
$\mathcal{Y}({\bm{\mu}},{\bm{z}})$

To prove that the set $\mathcal{Y}({\bm{\mu}},{\bm{z}})$ is not empty, we
first give a general lemma and its proof.

###### Lemma 10

Let $F(x,y)=f(x)+y^{\top}Ax+v^{\top}y-\frac{1}{2}y^{\top}Qy$, where $f(x)$ is
a strongly convex, smooth function, $Q$ is a positive semi-definite matrix and
$v\in\mathrm{Image}(Q)$. Then the saddle point of $F(x,y)$ exists, i.e., the
solution of the problem $P1$ exists

$\displaystyle P1:\min_{x}\max_{y}F(x,y).$ (A.1)

Proof of Lemma 10: First, by the basic property of quadratic optimization,
$\arg\max_{y}F(x,y)$ exists if and only if

$\displaystyle Ax\in\mathrm{Image}(Q).$ (A.2)

Let $\mathrm{rank}(Q)=r$. Without loss of generality, we assume
$Q=\begin{bmatrix}Q_{11}&Q_{12}\\\ Q_{12}&Q_{22}\end{bmatrix}$ and
$Q_{11}\in\mathbb{R}^{r\times r}$ is of full rank. Let
$P=\begin{bmatrix}Q_{11}\\\ Q_{12}\end{bmatrix}$ and
$Y=\\{y\in\mathbb{R}^{m}\mid y_{r+1}=\cdots=y_{m}=0\\}$.

Also let $\mathrm{Image}(Q)=\\{w\mid Bw=0\\}$, where the rows of
$B\in\mathbb{R}^{(m-r)\times m}$ are the basis for $\mathrm{Null}(Q)$. Define
$X=\\{x\in\mathbb{R}^{n}\mid BAx=0\\}$. Then, $x\in X$ implies
$x\in\mathrm{Image}(Q)$. For any $\tilde{x}\in X$, there exists some
$\tilde{y}\in Y$ such that

$\displaystyle\tilde{y}=\arg\max_{y\in\mathbb{R}^{m}}F(\tilde{x},y).$ (A.3)

In fact, the matrix $P$ is of rank $r$. Therefore,
$\mathrm{Image}(Q)=\mathrm{Image}(P)$. Notice that
$-A\tilde{x}-v\in\mathrm{Image}(Q)=\mathrm{Image}(P)$. It implies that there
exists some $u\in\mathbb{R}^{r}$ such that

$-A\tilde{x}-v=Pu.$

Hence, if $\tilde{y}=(u^{T},0,\cdots,0)^{T}\in Y\subseteq\mathbb{R}^{m}$, we
have

$\nabla_{y}F(\tilde{x},\tilde{y})=A\tilde{x}+v+Q\tilde{y}=0.$

Consequently, For any $\tilde{x}\in X$, there exists some $\tilde{y}\in Y$
such that

$\tilde{y}=\arg\max_{y\in Y}F(\tilde{x},y).$

Then consider the following constrained strongly convex-strongly concave min-
max problem:

$\displaystyle P2:\min_{x\in X}\max_{y\in Y}F(x,y).$ (A.4)

Note that $F(x,y)$ is strongly concave for $y$ in $Y$, thus $P2$ is a strongly
convex-strongly concave min-max problem. Since a strongly convex (concave)
function has a compact sublevel set, by [50, Prop. 3.6.9], there exists a
saddle point for problem $P2$. Suppose that $(\bar{x},\bar{y})$ is a saddle
point of P2. Since $\bar{x}\in X$, there exists $y\in Y$ such that,

$\bar{y}=\arg\max_{y\in
Y}F(\bar{x},y)=\arg\max_{y\in\mathbb{R}^{m}}F(\bar{x},y).$ (A.5)

Hence, $\nabla_{y}F(\bar{x},\bar{y})=0.$ We construct a saddle point of P1.
Since $\bar{x}=\arg\min_{x\in X}F(x,\bar{y})$, by the KKT conditions, there
exists some multiplier $\lambda\in\mathbb{R}^{m-r}$ such that

$\displaystyle\nabla f(\bar{x})+A^{T}\bar{y}+A^{T}B^{T}\lambda=0$
$\displaystyle BA\bar{x}=0.$

Let $x^{*}=\bar{x}$ and $y^{*}=\bar{y}+B^{T}\lambda$, we have

$\nabla_{x}F(x^{*},y^{*})=0$

and

$\displaystyle\nabla_{y}F(x^{*},y^{*})$ $\displaystyle=$ $\displaystyle-
Qy^{*}+v+Ax^{*}$ $\displaystyle=$
$\displaystyle(-Q\bar{y}+v+A\bar{x})+Q(\bar{y}-y^{*})$ $\displaystyle=$
$\displaystyle\nabla_{y}F(\bar{x},\bar{y})+Q(-B^{T}\lambda)$ $\displaystyle=$
$\displaystyle 0,$

where the last equality is because the rows of $B$ spans $\mathrm{Null}(Q)$.
We finish the proof. $\blacksquare$

Note that this lemma implies the nonempty property of the
$\mathcal{Y}({\bm{\mu}},{\bm{z}})$, where $f(x),y^{\top}Ax,Q,v$ correspond
$\textstyle\sum_{i=1}^{N}f_{i}({\bm{x}}_{i})+\frac{p}{2}\|{\bm{x}}_{i}-{\bm{z}}_{i}\|^{2}$,
$\sum_{i=1}^{N}{\bm{y}}_{i}^{\top}{\bm{B}}_{i}{\bm{x}}_{i}-{\bm{y}}^{\top}{\bm{q}}/N$,
$\rho{\bm{L}}^{-}$, ${\bm{A}}^{\top}{\bm{\mu}}$, respectively.

##  Appendix B Proof of Lemma 1

Notice that the change of the ${\mathcal{L}}_{\rho}$ function after one
iteration can be bounded as

$\displaystyle-{\mathcal{L}}_{\rho}({\bm{y}}^{r+1},{\bm{\mu}}^{r+1};{\bm{z}}^{r+1})+{\mathcal{L}}_{\rho}({\bm{y}}^{r},{\bm{\mu}}^{r};{\bm{z}}^{r})$
$\displaystyle\leq\bigg{(}-{\mathcal{L}}_{\rho}({\bm{y}}^{r},{\bm{\mu}}^{r+1};{\bm{z}}^{r})+{\mathcal{L}}_{\rho}({\bm{y}}^{r},{\bm{\mu}}^{r};{\bm{z}}^{r})\bigg{)}+\bigg{(}-{\mathcal{L}}_{\rho}({\bm{y}}^{r+1},{\bm{\mu}}^{r+1};{\bm{z}}^{r})+{\mathcal{L}}_{\rho}({\bm{y}}^{r},{\bm{\mu}}^{r+1};{\bm{z}}^{r})\bigg{)}$
$\displaystyle~{}~{}~{}+\bigg{(}-{\mathcal{L}}_{\rho}({\bm{y}}^{r+1},{\bm{\mu}}^{r+1};{\bm{z}}^{r+1})+{\mathcal{L}}_{\rho}({\bm{y}}^{r+1},{\bm{\mu}}^{r+1};{\bm{z}}^{r})\bigg{)}.$
(B.6)

Then, we will analyze this desired bound by following four steps below.

  1. (a)

Bound of
$-{\mathcal{L}}_{\rho}({\bm{y}}^{r},{\bm{\mu}}^{r+1};{\bm{z}}^{r})+{\mathcal{L}}_{\rho}({\bm{y}}^{r},{\bm{\mu}}^{r};{\bm{z}}^{r})$:

$\displaystyle-{\mathcal{L}}_{\rho}({\bm{y}}^{r},{\bm{\mu}}^{r+1};{\bm{z}}^{r})+{\mathcal{L}}_{\rho}({\bm{y}}^{r},{\bm{\mu}}^{r};{\bm{z}}^{r})={\langle}{\bm{\mu}}^{r+1}-{\bm{\mu}}^{r},{\bm{A}}{\bm{y}}^{r}{\rangle}=\frac{1}{\alpha}\|{\bm{\mu}}^{r+1}-{\bm{\mu}}^{r}\|^{2},$
(B.7)

where the last equality is obtained from (19a).

  2. (b)

Bound of
$-{\mathcal{L}}_{\rho}({\bm{y}}^{r+1},{\bm{\mu}}^{r+1};{\bm{z}}^{r})+{\mathcal{L}}_{\rho}({\bm{y}}^{r},{\bm{\mu}}^{r+1};{\bm{z}}^{r})$:
Since
${\mathcal{L}}_{\rho}({\bm{y}},{\bm{\mu}}^{r+1};{\bm{z}}^{r})-\frac{1}{2}\|{\bm{y}}-{\bm{y}}^{r}\|_{\rho{\bm{L}}^{+}}^{2}$
is strongly concave with modulus $\rho$, and

$\displaystyle\nabla_{y}{\mathcal{L}}_{\rho}({\bm{y}}^{r+1},{\bm{\mu}}^{r+1};{\bm{z}}^{r})-\rho{\bm{L}}^{+}({\bm{y}}^{r+1}-{\bm{y}}^{r})={\bm{0}},$
(B.8)

by the optimality condition of (19b), we have

$\displaystyle{\mathcal{L}}_{\rho}({\bm{y}}^{r},{\bm{\mu}}^{r+1};{\bm{z}}^{r})-\frac{1}{2}\|{\bm{y}}^{r}-{\bm{y}}^{r}\|_{\rho{\bm{L}}^{+}}^{2}\leq{\mathcal{L}}_{\rho}({\bm{y}}^{r+1},{\bm{\mu}}^{r+1};{\bm{z}}^{r})-\frac{1}{2}\|{\bm{y}}^{r+1}-{\bm{y}}^{r}\|_{\rho{\bm{L}}^{+}}^{2}-\frac{\rho}{2}\|{\bm{y}}^{r+1}-{\bm{y}}^{r}\|^{2}.$
(B.9)

  3. (c)

Bound of
$-{\mathcal{L}}_{\rho}({\bm{y}}^{r+1},{\bm{\mu}}^{r+1};{\bm{z}}^{r+1})+{\mathcal{L}}_{\rho}({\bm{y}}^{r+1},{\bm{\mu}}^{r+1};{\bm{z}}^{r})$:
Note that

$\displaystyle-{\mathcal{L}}_{\rho}({\bm{y}}^{r+1},{\bm{\mu}}^{r+1};{\bm{z}}^{r+1})+{\mathcal{L}}_{\rho}({\bm{y}}^{r+1},{\bm{\mu}}^{r+1};{\bm{z}}^{r})$
$\displaystyle=\sum_{i=1}^{N}\bigg{(}\phi_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r})-\phi_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r+1})\bigg{)}$
$\displaystyle=\sum_{i=1}^{N}\bigg{(}f_{i}({\bm{x}}_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r}))+\frac{p}{2}\|{\bm{x}}_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r})-{\bm{z}}_{i}^{r}\|^{2}+{\langle}{\bm{y}}_{i}^{r+1},{\bm{B}}_{i}{\bm{x}}_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r}){\rangle}\bigg{)}$
$\displaystyle~{}~{}~{}-\sum_{i=1}^{N}\bigg{(}f_{i}({\bm{x}}_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r+1}))+\frac{p}{2}\|{\bm{x}}_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r+1})-{\bm{z}}_{i}^{r+1}\|^{2}+{\langle}{\bm{y}}_{i}^{r+1},{\bm{B}}_{i}{\bm{x}}_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r+1}){\rangle}\bigg{)}.$
(B.10)

Since ${\bm{x}}_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r})$ is a minimizer of
(3), we have the upper bound

$\displaystyle-{\mathcal{L}}_{\rho}({\bm{y}}^{r+1},{\bm{\mu}}^{r+1};{\bm{z}}^{r+1})+{\mathcal{L}}_{\rho}({\bm{y}}^{r+1},{\bm{\mu}}^{r+1};{\bm{z}}^{r})$
$\displaystyle\leq\sum_{i=1}^{N}\bigg{(}f_{i}({\bm{x}}_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r+1}))+\frac{p}{2}\|{\bm{x}}_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r+1})-{\bm{z}}_{i}^{r}\|^{2}+{\langle}{\bm{y}}_{i}^{r+1},{\bm{B}}_{i}{\bm{x}}_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r+1}){\rangle}\bigg{)}$
$\displaystyle~{}~{}~{}-\sum_{i=1}^{N}\bigg{(}f_{i}({\bm{x}}_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r+1}))+\frac{p}{2}\|{\bm{x}}_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r+1})-{\bm{z}}_{i}^{r+1}\|^{2}+{\langle}{\bm{y}}_{i}^{r+1},{\bm{B}}_{i}{\bm{x}}_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r+1}){\rangle}\bigg{)}$
$\displaystyle=\frac{p}{2}\sum_{i=1}^{N}\|{\bm{x}}_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r+1})-{\bm{z}}_{i}^{r}\|^{2}-\frac{p}{2}\sum_{i=1}^{N}\|{\bm{x}}_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r+1})-{\bm{z}}_{i}^{r+1}\|^{2}$
$\displaystyle=-\frac{p}{2}\sum_{i=1}^{N}{\langle}{\bm{z}}_{i}^{r+1}-{\bm{z}}_{i}^{r},{\bm{z}}_{i}^{r+1}+{\bm{z}}_{i}^{r}-2{\bm{x}}_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r+1}){\rangle}.$
(B.11)

Summing (B.7), (B.9) and (B.11) leads to the desired result. $\blacksquare$

##  Appendix C Proof of Lemma 2

For any $\bar{}{\bm{y}}^{r}\in\mathcal{Y}({\bm{\mu}}^{r+1},{\bm{z}}^{r})$ and
$\tilde{}{\bm{y}}^{r}\in\mathcal{Y}({\bm{\mu}}^{r},{\bm{z}}^{r})$,

$\displaystyle
d({\bm{\mu}}^{r+1};{\bm{z}}^{r})-d({\bm{\mu}}^{r};{\bm{z}}^{r})$
$\displaystyle=\bigg{[}\sum_{i=1}^{N}\bigg{(}\phi_{i}(\bar{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r})-{\langle}\bar{}{\bm{y}}_{i}^{r},{\bm{b}}/N{\rangle}\bigg{)}-{\langle}{\bm{\mu}}^{r+1},{\bm{A}}\bar{}{\bm{y}}^{r}{\rangle}-\frac{\rho}{2}\|{\bm{A}}\bar{}{\bm{y}}^{r}\|^{2}\bigg{]}$
$\displaystyle~{}~{}~{}-\bigg{[}\sum_{i=1}^{N}\bigg{(}\phi_{i}(\tilde{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r})-{\langle}\tilde{}{\bm{y}}_{i}^{r},{\bm{b}}/N{\rangle}\bigg{)}-{\langle}{\bm{\mu}}^{r},{\bm{A}}\tilde{}{\bm{y}}^{r}{\rangle}\bigg{.}-\frac{\rho}{2}\|{\bm{A}}\tilde{}{\bm{y}}^{r}\|^{2}\bigg{]}$
$\displaystyle\overset{\mathrm{(i)}}{\leq}\bigg{[}\sum_{i=1}^{N}\bigg{(}\phi_{i}(\bar{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r})-{\langle}\bar{}{\bm{y}}_{i}^{r},{\bm{b}}/N{\rangle}\bigg{)}-{\langle}{\bm{\mu}}^{r+1},{\bm{A}}\bar{}{\bm{y}}^{r}{\rangle}-\frac{\rho}{2}\|{\bm{A}}\bar{}{\bm{y}}^{r}\|^{2}\bigg{]}$
$\displaystyle~{}~{}~{}-\bigg{[}\sum_{i=1}^{N}\bigg{(}\phi_{i}(\bar{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r})-{\langle}\bar{}{\bm{y}}_{i}^{r},{\bm{b}}/N{\rangle}\bigg{)}-{\langle}{\bm{\mu}}^{r},{\bm{A}}\bar{}{\bm{y}}^{r}{\rangle}\bigg{.}-\frac{\rho}{2}\|{\bm{A}}\bar{}{\bm{y}}^{r}\|^{2}\bigg{]}$
$\displaystyle=-{\langle}{\bm{\mu}}^{r+1}-{\bm{\mu}}^{r},{\bm{A}}\bar{}{\bm{y}}^{r}{\rangle}$
$\displaystyle=-\alpha{\langle}{\bm{A}}{\bm{y}}^{r},{\bm{A}}\bar{}{\bm{y}}^{r}{\rangle},$
(C.12)

where $\mathrm{(i)}$ is due to the fact that $\tilde{}{\bm{y}}^{r}$ is a
maximizer of ${\mathcal{L}}_{\rho}({\bm{y}},{\bm{\mu}}^{r};{\bm{z}}^{r})$, and
the last equality is deduced by (19a).

On the other hand, we have

$\displaystyle
d({\bm{\mu}}^{r+1};{\bm{z}}^{r+1})-d({\bm{\mu}}^{r+1};{\bm{z}}^{r})$
$\displaystyle=\bigg{[}\sum_{i=1}^{N}\bigg{(}\phi_{i}(\bar{}{\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r+1})-{\langle}\bar{}{\bm{y}}_{i}^{r+1},{\bm{b}}/N{\rangle}\bigg{)}-{\langle}{\bm{\mu}}^{r+1},{\bm{A}}\bar{}{\bm{y}}^{r+1}{\rangle}-\frac{\rho}{2}\|{\bm{A}}\bar{}{\bm{y}}^{r+1}\|^{2}\bigg{]}$
$\displaystyle~{}~{}~{}-\bigg{[}\sum_{i=1}^{N}\bigg{(}\phi_{i}(\bar{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r})-{\langle}\bar{}{\bm{y}}_{i}^{r},{\bm{b}}/N{\rangle}\bigg{)}-{\langle}{\bm{\mu}}^{r+1},{\bm{A}}\bar{}{\bm{y}}^{r}{\rangle}-\frac{\rho}{2}\|{\bm{A}}\bar{}{\bm{y}}^{r}\|^{2}\bigg{]}$
$\displaystyle\overset{\mathrm{(i)}}{\leq}\bigg{[}\sum_{i=1}^{N}\bigg{(}\phi_{i}(\bar{}{\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r+1})-{\langle}\bar{}{\bm{y}}_{i}^{r+1},{\bm{b}}/N{\rangle}\bigg{)}-{\langle}{\bm{\mu}}^{r+1},{\bm{A}}\bar{}{\bm{y}}^{r+1}{\rangle}-\frac{\rho}{2}\|{\bm{A}}\bar{}{\bm{y}}^{r+1}\|^{2}\bigg{]}$
$\displaystyle~{}~{}~{}-\bigg{[}\sum_{i=1}^{N}\bigg{(}\phi_{i}(\bar{}{\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r})-{\langle}\bar{}{\bm{y}}_{i}^{r+1},{\bm{b}}/N{\rangle}\bigg{)}-{\langle}{\bm{\mu}}^{r+1},{\bm{A}}\bar{}{\bm{y}}^{r+1}{\rangle}-\frac{\rho}{2}\|{\bm{A}}\bar{}{\bm{y}}^{r+1}\|^{2}\bigg{]}$
$\displaystyle=\sum_{i=1}^{N}\bigg{(}\phi_{i}(\bar{}{\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r+1})-\phi_{i}(\bar{}{\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r})\bigg{)},$
(C.13)

where $\mathrm{(i)}$ is obtained from the fact that $\bar{}{\bm{y}}^{r}$ is a
maximizer of ${\mathcal{L}}_{\rho}({\bm{y}},{\bm{\mu}}^{r+1};{\bm{z}}^{r})$.
Then according to (16), we have

$\displaystyle
d({\bm{\mu}}^{r+1};{\bm{z}}^{r+1})-d({\bm{\mu}}^{r+1};{\bm{z}}^{r})$
$\displaystyle\leq\sum_{i=1}^{N}\bigg{(}f_{i}({\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r+1}))+\frac{p}{2}\|{\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r+1})-{\bm{z}}_{i}^{r+1}\|^{2}+{\langle}\bar{}{\bm{y}}_{i}^{r+1},{\bm{B}}_{i}{\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r+1}){\rangle}\bigg{)}$
$\displaystyle~{}~{}~{}-\sum_{i=1}^{N}\bigg{(}f_{i}({\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r}))+\frac{p}{2}\|{\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r})-{\bm{z}}_{i}^{r}\|^{2}+{\langle}\bar{}{\bm{y}}_{i}^{r+1},{\bm{B}}_{i}{\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r}){\rangle}\bigg{)}$
$\displaystyle\overset{\mathrm{(i)}}{\leq}\sum_{i=1}^{N}\bigg{(}f_{i}({\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r}))+\frac{p}{2}\|{\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r})-{\bm{z}}_{i}^{r+1}\|^{2}+{\langle}\bar{}{\bm{y}}_{i}^{r+1},{\bm{B}}_{i}{\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r}){\rangle}\bigg{)}$
$\displaystyle~{}~{}~{}-\sum_{i=1}^{N}\bigg{(}f_{i}({\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r}))+\frac{p}{2}\|{\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r})-{\bm{z}}_{i}^{r}\|^{2}+{\langle}\bar{}{\bm{y}}_{i}^{r+1},{\bm{B}}_{i}{\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r}){\rangle}\bigg{)}$
$\displaystyle=\frac{p}{2}\sum_{i=1}^{N}\bigg{(}\|{\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r})-{\bm{z}}_{i}^{r+1}\|^{2}-\|{\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r})-{\bm{z}}_{i}^{r}\|^{2}\bigg{)}$
$\displaystyle=\frac{p}{2}\sum_{i=1}^{N}{\langle}{\bm{z}}_{i}^{r+1}-{\bm{z}}^{r},-2{\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r})+{\bm{z}}_{i}^{r+1}+{\bm{z}}_{i}^{r}{\rangle},$
(C.14)

where $\mathrm{(i)}$ is due to the fact that
${\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r+1})$ is a minimizer of
$\phi_{i}(\bar{}{\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r+1})$.

$\blacksquare$

##  Appendix D Proof of Lemma 3

The proof of Lemma 3 is a direct corollary of the following two lemmas.

###### Lemma 11

We have

$\displaystyle
G({\bm{x}}^{r+1},{\bm{y}}^{r+1},{\bm{\mu}}^{r+1},{\bm{z}}^{r+1})-G({\bm{x}}^{r},{\bm{y}}^{r},{\bm{\mu}}^{r},{\bm{z}}^{r})$
$\displaystyle\leq-\frac{1}{\alpha}\|{\bm{\mu}}^{r+1}-{\bm{\mu}}^{r}\|^{2}-\sum_{i=1}^{N}\frac{p+\gamma^{-}}{2}\|{\bm{x}}_{i}^{r+1}-{\bm{x}}_{i}^{r}\|^{2}-\frac{\rho}{2}\|{\bm{y}}^{r+1}-{\bm{y}}^{r}\|^{2}_{{\bm{L}}^{-}}+\sum_{i=1}^{N}{\langle}{\bm{y}}_{i}^{r+1}-{\bm{y}}_{i}^{r},{\bm{B}}_{i}({\bm{x}}_{i}^{r}-{\bm{x}}_{i}^{r+1}){\rangle}$
$\displaystyle~{}~{}~{}+\sum_{i=1}^{N}(2\rho|\mathcal{N}_{i}|)\|{\bm{y}}_{i}^{r+1}-{\bm{y}}_{i}^{r}\|^{2}+\frac{p}{2}\bigg{(}1-\frac{2}{\beta}\bigg{)}\|{\bm{z}}^{r+1}-{\bm{z}}^{r}\|^{2}.$
(D.15)

###### Lemma 12

$\displaystyle\bigg{(}\frac{\alpha}{2}\|{\bm{y}}^{r+1}\|^{2}_{{\bm{L}}^{-}}+\frac{1}{2}\|{\bm{y}}^{r+1}-{\bm{y}}^{r}\|^{2}_{\rho{\bm{L}}^{+}}\bigg{)}-\bigg{(}\frac{\alpha}{2}\|{\bm{y}}^{r}\|^{2}_{{\bm{L}}^{-}}+\frac{1}{2}\|{\bm{y}}^{r}-{\bm{y}}^{r-1}\|^{2}_{\rho{\bm{L}}^{+}}\bigg{)}$
$\displaystyle\leq-\sum_{i=1}^{N}{\langle}{\bm{y}}_{i}^{r+1}-{\bm{y}}_{i}^{r},{\bm{B}}_{i}({\bm{x}}_{i}^{r}-{\bm{x}}_{i}^{r+1}){\rangle}-\frac{1}{2}\|({\bm{y}}^{r}-{\bm{y}}^{r-1})-({\bm{y}}^{r+1}-{\bm{y}}^{r})\|^{2}_{\rho{\bm{L}}^{+}}-\bigg{(}\rho-\frac{\alpha}{2}\bigg{)}\|{\bm{y}}^{r+1}-{\bm{y}}^{r}\|^{2}_{{\bm{L}}^{-}}.$
(D.16)

We then give the proofs of the above two lemmas.

Proof of Lemma 11: Notice that the change of the function $G$ after one
iteration can be bounded as

$\displaystyle
G({\bm{x}}^{r+1},{\bm{y}}^{r+1},{\bm{\mu}}^{r+1},{\bm{z}}^{r+1})-G({\bm{x}}^{r},{\bm{y}}^{r},{\bm{\mu}}^{r},{\bm{z}}^{r})$
$\displaystyle\leq\bigg{(}G({\bm{x}}^{r},{\bm{y}}^{r},{\bm{\mu}}^{r+1},{\bm{z}}^{r})-G({\bm{x}}^{r},{\bm{y}}^{r},{\bm{\mu}}^{r},{\bm{z}}^{r})\bigg{)}+\bigg{(}G({\bm{x}}^{r+1},{\bm{y}}^{r},{\bm{\mu}}^{r+1},{\bm{z}}^{r})-G({\bm{x}}^{r},{\bm{y}}^{r},{\bm{\mu}}^{r+1},{\bm{z}}^{r})\bigg{)}$
$\displaystyle~{}~{}~{}+\bigg{(}G({\bm{x}}^{r+1},{\bm{y}}^{r+1},{\bm{\mu}}^{r+1},{\bm{z}}^{r})-G({\bm{x}}^{r+1},{\bm{y}}^{r},{\bm{\mu}}^{r+1},{\bm{z}}^{r})\bigg{)}$
$\displaystyle~{}~{}~{}+\bigg{(}G({\bm{x}}^{r+1},{\bm{y}}^{r+1},{\bm{\mu}}^{r+1},{\bm{z}}^{r+1})-G({\bm{x}}^{r+1},{\bm{y}}^{r+1},{\bm{\mu}}^{r+1},{\bm{z}}^{r})\bigg{)}.$
(D.17)

Next, we will bound four terms in the right hand side (RHS) of the above
inequality respectively.

  1. (a)

Bound of
$G({\bm{x}}^{r},{\bm{y}}^{r},{\bm{\mu}}^{r+1},{\bm{z}}^{r})-G({\bm{x}}^{r},{\bm{y}}^{r},{\bm{\mu}}^{r},{\bm{z}}^{r})$:

$\displaystyle
G({\bm{x}}^{r},{\bm{y}}^{r},{\bm{\mu}}^{r+1},{\bm{z}}^{r})-G({\bm{x}}^{r},{\bm{y}}^{r},{\bm{\mu}}^{r},{\bm{z}}^{r})=-{\langle}{\bm{\mu}}^{r+1}-{\bm{\mu}}^{r},{\bm{A}}{\bm{y}}^{r}{\rangle}=-\frac{1}{\alpha}\|{\bm{\mu}}^{r+1}-{\bm{\mu}}^{r}\|^{2},$
(D.18)

where the last equality is due to (19a).

  2. (b)

Bound of
$G({\bm{x}}^{r+1},{\bm{y}}^{r},{\bm{\mu}}^{r+1},{\bm{z}}^{r})-G({\bm{x}}^{r},{\bm{y}}^{r},{\bm{\mu}}^{r+1},{\bm{z}}^{r})$:

Firstly, by the optimality condition of (16) , we have

$\displaystyle{\langle}\nabla
f_{i}({\bm{x}}_{i}^{r+1})+p({\bm{x}}_{i}^{r+1}-{\bm{z}}_{i}^{r})+{\bm{B}}_{i}^{\top}{\bm{y}}_{i}^{r+1},{\bm{x}}_{i}-{\bm{x}}_{i}^{r+1}{\rangle}\geq
0,~{}\forall~{}{\bm{x}}_{i}.$ (D.19)

Secondly, since
$f_{i}({\bm{x}}_{i})+\frac{p}{2}\|{\bm{x}}_{i}-{\bm{z}}_{i}^{r}\|^{2}$ is
strongly convex with modulus $p+\gamma^{-}$, we obtain

$\displaystyle
f_{i}({\bm{x}}_{i}^{r})+\frac{p}{2}\|{\bm{x}}_{i}^{r}-{\bm{z}}_{i}^{r}\|^{2}$
$\displaystyle\geq
f_{i}({\bm{x}}_{i}^{r+1})+\frac{p}{2}\|{\bm{x}}_{i}^{r+1}-{\bm{z}}_{i}^{r}\|^{2}+{\langle}\nabla
f_{i}({\bm{x}}_{i}^{r+1})+p({\bm{x}}_{i}^{r+1}-{\bm{z}}_{i}^{r}),{\bm{x}}_{i}^{r}-{\bm{x}}_{i}^{r+1}{\rangle}+\frac{p+\gamma^{-}}{2}\|{\bm{x}}_{i}^{r}-{\bm{x}}_{i}^{r+1}\|^{2}$
$\displaystyle\geq
f_{i}({\bm{x}}_{i}^{r+1})+\frac{p}{2}\|{\bm{x}}_{i}^{r+1}-{\bm{z}}_{i}^{r}\|^{2}+{\langle}-{\bm{B}}_{i}^{\top}{\bm{y}}_{i}^{r+1},{\bm{x}}_{i}^{r}-{\bm{x}}_{i}^{r+1}{\rangle}+\frac{p+\gamma^{-}}{2}\|{\bm{x}}_{i}^{r}-{\bm{x}}_{i}^{r+1}\|^{2},$
(D.20)

where the last inequality is due to (D.19).

Thus, we have the upper bound for the difference
$G({\bm{x}}^{r+1},{\bm{y}}^{r},{\bm{\mu}}^{r+1},{\bm{z}}^{r})-G({\bm{x}}^{r},{\bm{y}}^{r},{\bm{\mu}}^{r+1},{\bm{z}}^{r})$
as

$\displaystyle
G({\bm{x}}^{r+1},{\bm{y}}^{r},{\bm{\mu}}^{r+1},{\bm{z}}^{r})-G({\bm{x}}^{r},{\bm{y}}^{r},{\bm{\mu}}^{r+1},{\bm{z}}^{r})$
$\displaystyle=\sum_{i=1}^{N}\bigg{(}f_{i}({\bm{x}}_{i}^{r+1})+\frac{p}{2}\|{\bm{x}}_{i}^{r+1}-{\bm{z}}_{i}^{r}\|^{2}-f_{i}({\bm{x}}_{i}^{r})-\frac{p}{2}\|{\bm{x}}_{i}^{r}-{\bm{z}}_{i}^{r}\|^{2}\bigg{)}+\sum_{i=1}^{N}{\langle}{\bm{y}}_{i}^{r},{\bm{B}}_{i}({\bm{x}}_{i}^{r+1}-{\bm{x}}_{i}^{r}){\rangle}$
$\displaystyle\leq\sum_{i=1}^{N}{\langle}{\bm{y}}_{i}^{r+1}-{\bm{y}}_{i}^{r},{\bm{B}}_{i}({\bm{x}}_{i}^{r}-{\bm{x}}_{i}^{r+1}){\rangle}-\frac{p+\gamma^{-}}{2}\|{\bm{x}}_{i}^{r}-{\bm{x}}_{i}^{r+1}\|^{2},$
(D.21)

where the last equality is obtained by (D.20).

  3. (c)

Bound of
$G({\bm{x}}^{r+1},{\bm{y}}^{r+1},{\bm{\mu}}^{r+1},{\bm{z}}^{r})-G({\bm{x}}^{r+1},{\bm{y}}^{r},{\bm{\mu}}^{r+1},{\bm{z}}^{r})$:

Since

$\displaystyle\frac{\rho}{2}\|{\bm{A}}{\bm{y}}^{r+1}\|^{2}$
$\displaystyle\geq\frac{\rho}{2}\|{\bm{A}}{\bm{y}}^{r}\|^{2}+\rho{\langle}{\bm{L}}^{-}{\bm{y}}^{r},{\bm{y}}^{r+1}-{\bm{y}}^{r}{\rangle}+\frac{\rho}{2}\|{\bm{y}}^{r+1}-{\bm{y}}^{r}\|^{2}_{{\bm{L}}^{-}},$
(D.22)

we can have

$\displaystyle
G({\bm{x}}^{r+1},{\bm{y}}^{r+1},{\bm{\mu}}^{r+1},{\bm{z}}^{r})-G({\bm{x}}^{r+1},{\bm{y}}^{r},{\bm{\mu}}^{r+1},{\bm{z}}^{r})$
$\displaystyle=\sum_{i=1}^{N}{\langle}{\bm{y}}_{i}^{r+1}-{\bm{y}}_{i}^{r},{\bm{B}}_{i}{\bm{x}}_{i}-\frac{{\bm{q}}}{N}{\rangle}-{\langle}{\bm{\mu}}^{r+1},{\bm{A}}({\bm{y}}^{r+1}-{\bm{y}}^{r}){\rangle}-\frac{\rho}{2}\|{\bm{A}}{\bm{y}}^{r+1}\|^{2}+\frac{\rho}{2}\|{\bm{A}}{\bm{y}}^{r}\|^{2}$
$\displaystyle=\sum_{i=1}^{N}{\langle}{\bm{y}}_{i}^{r+1}-{\bm{y}}_{i}^{r},{\bm{B}}_{i}{\bm{x}}_{i}-\frac{{\bm{q}}}{N}-{\bm{A}}_{i}^{\top}{\bm{\mu}}^{r+1}{\rangle}-\frac{\rho}{2}\|{\bm{A}}{\bm{y}}^{r+1}\|^{2}+\frac{\rho}{2}\|{\bm{A}}{\bm{y}}^{r}\|^{2}$
$\displaystyle\leq\sum_{i=1}^{N}{\langle}{\bm{y}}_{i}^{r+1}-{\bm{y}}_{i}^{r},{\bm{B}}_{i}{\bm{x}}_{i}-\frac{{\bm{q}}}{N}-{\bm{A}}_{i}^{\top}{\bm{\mu}}^{r+1}{\rangle}-\rho{\langle}{\bm{L}}^{-}{\bm{y}}^{r},{\bm{y}}^{r+1}-{\bm{y}}^{r}{\rangle}-\frac{\rho}{2}\|{\bm{y}}^{r+1}-{\bm{y}}^{r}\|^{2}_{{\bm{L}}^{-}}.$

Using the fact that ${\bm{L}}^{+}=2{\bm{D}}-{\bm{L}}^{-}$, we further have

$\displaystyle
G({\bm{x}}^{r+1},{\bm{y}}^{r+1},{\bm{\mu}}^{r+1},{\bm{z}}^{r})-G({\bm{x}}^{r+1},{\bm{y}}^{r},{\bm{\mu}}^{r+1},{\bm{z}}^{r})$
$\displaystyle\leq\sum_{i=1}^{N}{\langle}{\bm{y}}_{i}^{r+1}-{\bm{y}}_{i}^{r},{\bm{B}}_{i}{\bm{x}}_{i}-\frac{{\bm{q}}}{N}-{\bm{A}}_{i}^{\top}{\bm{\mu}}^{r+1}+\rho{\bm{L}}_{i}^{+}{\bm{y}}^{r}{\rangle}-2\rho{\langle}{\bm{D}}{\bm{y}}^{r},{\bm{y}}^{r+1}-{\bm{y}}^{r}{\rangle}-\frac{\rho}{2}\|{\bm{y}}^{r+1}-{\bm{y}}^{r}\|^{2}_{{\bm{L}}^{-}}$
$\displaystyle\overset{\mathrm{(i)}}{=}\sum_{i=1}^{N}{\langle}{\bm{y}}_{i}^{r+1}-{\bm{y}}_{i}^{r},2\rho|\mathcal{N}_{i}|{\bm{y}}_{i}^{r+1}{\rangle}-2\rho{\langle}{\bm{D}}{\bm{y}}^{r},{\bm{y}}^{r+1}-{\bm{y}}^{r}{\rangle}-\frac{\rho}{2}\|{\bm{y}}^{r+1}-{\bm{y}}^{r}\|^{2}_{{\bm{L}}^{-}}$
$\displaystyle=\sum_{i=1}^{N}(2\rho|\mathcal{N}_{i}|)\|{\bm{y}}_{i}^{r+1}-{\bm{y}}_{i}^{r}\|^{2}-\frac{\rho}{2}\|{\bm{y}}^{r+1}-{\bm{y}}^{r}\|^{2}_{{\bm{L}}^{-}},$
(D.23)

where $\mathrm{(i)}$ is obtained by (23).

  4. (d)

Bound of
$G({\bm{x}}^{r+1},{\bm{y}}^{r+1},{\bm{\mu}}^{r+1},{\bm{z}}^{r+1})-G({\bm{x}}^{r+1},{\bm{y}}^{r+1},{\bm{\mu}}^{r+1},{\bm{z}}^{r})$:
It can be obtained as

$\displaystyle
G({\bm{x}}^{r+1},{\bm{y}}^{r+1},{\bm{\mu}}^{r+1},{\bm{z}}^{r+1})-G({\bm{x}}^{r+1},{\bm{y}}^{r+1},{\bm{\mu}}^{r+1},{\bm{z}}^{r})$
$\displaystyle=\frac{p}{2}\|{\bm{x}}^{r+1}-{\bm{z}}^{r}\|^{2}-\frac{p}{2}\|{\bm{x}}^{r+1}-{\bm{z}}^{r}\|^{2}$
$\displaystyle=\frac{p}{2}{\langle}-{\bm{z}}^{r+1}+{\bm{z}}^{r},2({\bm{x}}^{r+1}-{\bm{z}}^{r})-{\bm{z}}^{r+1}+{\bm{z}}^{r}{\rangle}$
$\displaystyle\overset{\mathrm{(i)}}{=}\frac{p}{2}{\langle}-{\bm{z}}^{r+1}+{\bm{z}}^{r},\frac{2}{\beta}({\bm{z}}^{r+1}-{\bm{z}}^{r})-{\bm{z}}^{r+1}+{\bm{z}}^{r}{\rangle}$
$\displaystyle=\frac{p}{2}\bigg{(}1-\frac{2}{\beta}\bigg{)}\|{\bm{z}}^{r+1}-{\bm{z}}^{r}\|^{2},$
(D.24)

where $\mathrm{(i)}$ is due to the updating step (30) for ${\bm{z}}$ of
algorithm 1.

Summing (D.18), (D.21), (D.23), and (D.24) gives the desired bound in (11).
$\blacksquare$

Proof of Lemma 12: Firstly, the optimality condition of (19b) for iteration
$r$ and $r-1$ are given as

$\displaystyle\bigg{\langle}\nabla_{{\bm{y}}}\left({\mathcal{L}}_{\rho}({\bm{y}}^{r+1},{\bm{\mu}}^{r+1};{\bm{z}}^{r})-\frac{\rho}{2}\|{\bm{y}}^{r+1}-{\bm{y}}^{r}\|_{{\bm{L}}^{+}}^{2}\right),{\bm{y}}-{\bm{y}}^{r+1}\bigg{\rangle}\leq
0,\forall{\bm{y}},$ (D.25)
$\displaystyle\bigg{\langle}\nabla_{{\bm{y}}}\left({\mathcal{L}}_{\rho}({\bm{y}}^{r},{\bm{\mu}}^{r+1};{\bm{z}}^{r})-\frac{\rho}{2}\|{\bm{y}}^{r}-{\bm{y}}^{r-1}\|_{{\bm{L}}^{+}}^{2}\right),{\bm{y}}-{\bm{y}}^{r}\bigg{\rangle}\leq
0,\forall{\bm{y}}.$ (D.26)

According to Danskin’s theorem [14, Proposition B.22], we have

$\displaystyle\nabla_{{\bm{y}}_{i}}\phi_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r})={\bm{B}}_{i}{\bm{x}}_{i}^{r+1}.$
(D.27)

Hence, the gradient of ${\mathcal{L}}_{\rho}({\bm{y}},{\bm{\mu}};{\bm{z}})$ is
given by

$\nabla_{{\bm{y}}_{i}}{\mathcal{L}}_{\rho}({\bm{y}}^{r+1},{\bm{\mu}}^{r+1};{\bm{z}}^{r})={\bm{B}}_{i}{\bm{x}}_{i}^{r+1}-\frac{{\bm{q}}}{N}-{\bm{A}}_{i}^{\top}{\bm{\mu}}^{r+1}-\rho{\bm{L}}_{i}^{-}{\bm{y}}_{i}^{r+1}.$

Substituting the above equality into (D.25), we have

$\displaystyle{\langle}{\bm{B}}_{i}{\bm{x}}_{i}^{r+1}-{\bm{b}}/N-{\bm{A}}_{i}^{\top}{\bm{\mu}}^{r+1}-\rho{\bm{L}}_{i}^{-}{\bm{y}}_{i}^{r+1}-\rho{\bm{L}}_{i}^{+}({\bm{y}}_{i}^{r+1}-{\bm{y}}_{i}^{r}),{\bm{y}}_{i}-{\bm{y}}_{i}^{r+1}{\rangle}\leq
0,~{}\forall~{}{\bm{y}}_{i},$ (D.28)
$\displaystyle{\langle}{\bm{B}}_{i}{\bm{x}}_{i}^{r}-{\bm{b}}/N-{\bm{A}}_{i}^{\top}{\bm{\mu}}^{r}-\rho{\bm{L}}_{i}^{-}{\bm{y}}_{i}^{r}-\rho{\bm{L}}_{i}^{+}({\bm{y}}_{i}^{r}-{\bm{y}}_{i}^{r-1}),{\bm{y}}_{i}-{\bm{y}}_{i}^{r}{\rangle}\leq
0,~{}\forall~{}{\bm{y}}_{i}.$ (D.29)

Letting ${\bm{y}}_{i}={\bm{y}}_{i}^{r}$ and ${\bm{y}}_{i}={\bm{y}}_{i}^{r+1}$
in the first and second equations above respectively, and summing up these two
inequalities, we obtain

$\displaystyle{\langle}{\bm{B}}_{i}({\bm{x}}_{i}^{r}-{\bm{x}}_{i}^{r+1}),{\bm{y}}_{i}^{r+1}-{\bm{y}}_{i}^{r}{\rangle}-{\langle}{\bm{A}}_{i}^{\top}({\bm{\mu}}^{r}-{\bm{\mu}}^{r+1}),{\bm{y}}_{i}^{r+1}-{\bm{y}}_{i}^{r}{\rangle}-{\langle}\rho{\bm{L}}_{i}^{-}({\bm{y}}_{i}^{r}-{\bm{y}}_{i}^{r+1}),{\bm{y}}_{i}^{r+1}-{\bm{y}}_{i}^{r}{\rangle}$
$\displaystyle-{\langle}(\rho{\bm{L}}_{i}^{+})(({\bm{y}}_{i}^{r}-{\bm{y}}_{i}^{r-1})-({\bm{y}}_{i}^{r+1}-{\bm{y}}_{i}^{r})),~{}{\bm{y}}_{i}^{r+1}-{\bm{y}}_{i}^{r}{\rangle}\leq
0,$ (D.30)

for $\forall i\in[N]$. By summing the $N$ equations, one obtains

$\displaystyle\sum_{i=1}^{N}{\langle}{\bm{B}}_{i}({\bm{x}}_{i}^{r}-{\bm{x}}_{i}^{r+1}),{\bm{y}}_{i}^{r+1}-{\bm{y}}_{i}^{r}{\rangle}$
$\displaystyle\leq{\langle}{\bm{\mu}}^{r}-{\bm{\mu}}^{r+1},{\bm{A}}({\bm{y}}^{r+1}-{\bm{y}}^{r}){\rangle}-\rho\|{\bm{y}}^{r+1}-{\bm{y}}^{r}\|_{{\bm{L}}^{-}}^{2}$
$\displaystyle~{}~{}~{}+{\langle}({\bm{y}}^{r}-{\bm{y}}^{r-1})-({\bm{y}}^{r+1}-{\bm{y}}^{r}),\rho{\bm{L}}^{+}({\bm{y}}^{r+1}-{\bm{y}}^{r}){\rangle}.$
(D.31)

Notice that

$\displaystyle{\langle}{\bm{\mu}}^{r}-{\bm{\mu}}^{r+1},{\bm{A}}({\bm{y}}^{r+1}-{\bm{y}}^{r}){\rangle}$
$\displaystyle\overset{\mathrm{(i)}}{=}{\langle}\alpha{\bm{A}}{\bm{y}}^{r},{\bm{A}}({\bm{y}}^{r}-{\bm{y}}^{r+1}){\rangle}$
$\displaystyle=\frac{\alpha}{2}\|{\bm{y}}^{r+1}-{\bm{y}}^{r}\|_{{\bm{L}}^{-}}^{2}+\frac{\alpha}{2}\|{\bm{y}}^{r}\|_{{\bm{L}}^{-}}^{2}-\frac{\alpha}{2}\|{\bm{y}}^{r+1}\|_{{\bm{L}}^{-}}^{2},$
(D.32)

in which $\mathrm{(i)}$ is obtained from (19a). Moreover, we have

$\displaystyle{\langle}({\bm{y}}^{r}-{\bm{y}}^{r-1})-({\bm{y}}^{r+1}-{\bm{y}}^{r}),\rho{\bm{L}}^{+}({\bm{y}}^{r+1}-{\bm{y}}^{r}){\rangle}$
$\displaystyle=\frac{1}{2}\|{\bm{y}}^{r}-{\bm{y}}^{r-1}\|^{2}_{\rho{\bm{L}}^{+}}-\frac{1}{2}\|({\bm{y}}^{r}-{\bm{y}}^{r-1})-({\bm{y}}^{r+1}-{\bm{y}}^{r})\|^{2}_{\rho{\bm{L}}^{+}}-\frac{1}{2}\|{\bm{y}}^{r+1}-{\bm{y}}^{r}\|^{2}_{\rho{\bm{L}}^{+}}.$
(D.33)

By substituting (D.32) and (D.33) into (D.31), we obtain

$\displaystyle\sum_{i=1}^{N}{\langle}{\bm{y}}_{i}^{r+1}-{\bm{y}}_{i}^{r},{\bm{B}}_{i}({\bm{x}}_{i}^{r}-{\bm{x}}_{i}^{r+1}){\rangle}$
$\displaystyle\leq\frac{\alpha}{2}\|{\bm{y}}^{r}\|^{2}_{{\bm{L}}^{-}}-\frac{\alpha}{2}\|{\bm{y}}^{r+1}\|^{2}_{{\bm{L}}^{-}}+\frac{1}{2}\|{\bm{y}}^{r}-{\bm{y}}^{r-1}\|^{2}_{\rho{\bm{L}}^{+}}-\frac{1}{2}\|{\bm{y}}^{r+1}-{\bm{y}}^{r}\|^{2}_{\rho{\bm{L}}^{+}}-\frac{1}{2}\|({\bm{y}}^{r}-{\bm{y}}^{r-1})-({\bm{y}}^{r+1}-{\bm{y}}^{r})\|^{2}_{\rho{\bm{L}}^{+}}$
$\displaystyle~{}~{}~{}-\bigg{(}\rho-\frac{\alpha}{2}\bigg{)}\|{\bm{y}}^{r+1}-{\bm{y}}^{r}\|^{2}_{{\bm{L}}^{-}}.$
(D.34)

After rearranging the inequality above, we have

$\displaystyle\bigg{(}\frac{\alpha}{2}\|{\bm{y}}^{r+1}\|^{2}_{{\bm{L}}^{-}}+\frac{1}{2}\|{\bm{y}}^{r+1}-{\bm{y}}^{r}\|^{2}_{\rho{\bm{L}}^{+}}\bigg{)}-\bigg{(}\frac{\alpha}{2}\|{\bm{y}}^{r}\|^{2}_{{\bm{L}}^{-}}+\frac{1}{2}\|{\bm{y}}^{r}-{\bm{y}}^{r-1}\|^{2}_{\rho{\bm{L}}^{+}}\bigg{)}$
$\displaystyle\leq-\sum_{i=1}^{N}{\langle}{\bm{y}}_{i}^{r+1}-{\bm{y}}_{i}^{r},{\bm{B}}_{i}({\bm{x}}_{i}^{r}-{\bm{x}}_{i}^{r+1}){\rangle}-\frac{1}{2}\|({\bm{y}}^{r}-{\bm{y}}^{r-1})-({\bm{y}}^{r+1}-{\bm{y}}^{r})\|^{2}_{\rho{\bm{L}}^{+}}-\bigg{(}\rho-\frac{\alpha}{2}\bigg{)}\|{\bm{y}}^{r+1}-{\bm{y}}^{r}\|^{2}_{{\bm{L}}^{-}}.$
(D.35)

$\blacksquare$

##  Appendix E Proof of Lemma 4

###  E.1 Proof of (49) in Lemma 4

The proof of (49) is similar to [17, Lemma 3.6]. Recall from (41) that

$\displaystyle{\bm{x}}_{i}({\bm{y}}_{i},{\bm{z}}_{i})=\arg\min_{{\bm{x}}_{i}}\bigg{\\{}f_{i}({\bm{x}}_{i})+\frac{p}{2}\|{\bm{x}}_{i}-{\bm{z}}_{i}\|^{2}+{\bm{y}}_{i}^{\top}{\bm{B}}_{i}{\bm{x}}_{i}\bigg{\\}},$
(E.36)

which is a strongly convex problem with modulus $p+\gamma^{-}$. Thus,
according to the error bound for strongly convex problems in [51, Theorem
3.1], we have

$\displaystyle\|{\bm{x}}_{i}-{\bm{x}}_{i}({\bm{y}}_{i},{\bm{z}}_{i})\|\leq\frac{(p+\gamma^{+}+1)}{p+\gamma^{-}}\|\nabla
f_{i}({\bm{x}}_{i})+p({\bm{x}}_{i}-{\bm{z}}_{i})+{\bm{B}}_{i}^{\top}{\bm{y}}_{i}\|,$
(E.37)

for all ${\bm{x}}_{i}\in\mathbb{R}^{n}$. Since

$\displaystyle\nabla
f_{i}({\bm{x}}_{i}({\bm{y}}_{i}^{\prime},{\bm{z}}_{i}^{\prime}))+p({\bm{x}}_{i}({\bm{y}}_{i}^{\prime},{\bm{z}}_{i}^{\prime})-{\bm{z}}_{i}^{\prime})+{\bm{B}}_{i}^{\top}{\bm{y}}_{i}^{\prime}=0,$
(E.38)

by letting
${\bm{x}}_{i}={\bm{x}}_{i}({\bm{y}}_{i}^{\prime},{\bm{z}}_{i}^{\prime})$ in
(E.37), we obtain

$\displaystyle\|{\bm{x}}_{i}({\bm{y}}_{i}^{\prime},{\bm{z}}_{i}^{\prime})-{\bm{x}}_{i}({\bm{y}}_{i},{\bm{z}}_{i})\|\leq\frac{(p+\gamma^{+}+1)}{p+\gamma^{-}}\|\nabla
f_{i}({\bm{x}}_{i}({\bm{y}}_{i}^{\prime},{\bm{z}}_{i}^{\prime}))+p({\bm{x}}_{i}({\bm{y}}_{i}^{\prime},{\bm{z}}_{i}^{\prime})-{\bm{z}}_{i})+{\bm{B}}_{i}^{\top}{\bm{y}}_{i}\|.$
(E.39)

By substituting (E.38) into (E.39), we have

$\displaystyle\|{\bm{x}}_{i}({\bm{y}}_{i}^{\prime},{\bm{z}}_{i}^{\prime})-{\bm{x}}_{i}({\bm{y}}_{i},{\bm{z}}_{i})\|$
$\displaystyle\leq\frac{(p+\gamma^{+}+1)}{p+\gamma^{-}}\|p({\bm{z}}_{i}-{\bm{z}}_{i}^{\prime})+{\bm{B}}_{i}^{\top}({\bm{y}}_{i}-{\bm{y}}_{i}^{\prime})\|$
$\displaystyle\overset{\mathrm{(i)}}{\leq}\frac{p(p+\gamma^{+}+1)}{p+\gamma^{-}}\|{\bm{z}}_{i}-{\bm{z}}_{i}^{\prime}\|+\frac{\|{\bm{B}}_{i}\|(p+\gamma^{+}+1)}{p+\gamma^{-}}\|{\bm{y}}_{i}-{\bm{y}}_{i}^{\prime}\|$
$\displaystyle\leq\frac{p(p+\gamma^{+}+1)}{p+\gamma^{-}}\|{\bm{z}}_{i}-{\bm{z}}_{i}^{\prime}\|+\frac{B_{\max}(p+\gamma^{+}+1)}{p+\gamma^{-}}\|{\bm{y}}_{i}-{\bm{y}}_{i}^{\prime}\|,$
(E.40)

where $\mathrm{(i)}$ is obtained by the triangle inequality and Cauchy-
Schwartz inequality. $\blacksquare$

###  E.2 Proof of (50) in Lemma 4

To show (50), let

$\bar{h}({\bm{x}};{\bm{\mu}},{\bm{z}})=\max_{{\bm{y}}}\bigg{\\{}G({\bm{x}},{\bm{y}},{\bm{\mu}},{\bm{z}})-\frac{p}{2}\|{\bm{z}}\|^{2}\bigg{\\}}.$

By (41), ${\bm{x}}(\bar{}{\bm{y}},{\bm{z}})$ is the minimizer of the function
$\bar{h}({\bm{x}};{\bm{\mu}},{\bm{z}})$, where
$\bar{}{\bm{y}}\in\mathcal{Y}({\bm{\mu}},{\bm{z}})$. Then for any
$\bar{}{\bm{y}}^{r}\in\mathcal{Y}({\bm{\mu}}^{r+1},{\bm{z}}^{r})$ and
$\bar{}{\bm{y}}^{r+1}\in\mathcal{Y}({\bm{\mu}}^{r+1},{\bm{z}}^{r+1})$, we can
have
$\nabla\bar{h}({\bm{x}}(\bar{}{\bm{y}}^{r},{\bm{z}}^{r});{\bm{\mu}}^{r+1},{\bm{z}}^{r})={\bm{0}}$
and
$\nabla\bar{h}({\bm{x}}(\bar{}{\bm{y}}^{r+1},{\bm{z}}^{r+1});{\bm{\mu}}^{r+1},{\bm{z}}^{r+1})={\bm{0}}$.
Due to the fact that $\bar{h}$ is strongly convex with modulus
$(p+\gamma^{-})$, we have

$\displaystyle\bar{h}({\bm{x}}(\bar{}{\bm{y}}^{r+1},{\bm{z}}^{r+1});{\bm{\mu}}^{r+1},{\bm{z}}^{r})-\bar{h}({\bm{x}}(\bar{}{\bm{y}}^{r},{\bm{z}}^{r});{\bm{\mu}}^{r+1},{\bm{z}}^{r})\geq\frac{p+\gamma^{-}}{2}\|{\bm{x}}(\bar{}{\bm{y}}^{r},{\bm{z}}^{r})-{\bm{x}}(\bar{}{\bm{y}}^{r+1},{\bm{z}}^{r+1})\|^{2},$
(E.41)

and

$\displaystyle\bar{h}({\bm{x}}(\bar{}{\bm{y}}^{r},{\bm{z}}^{r});{\bm{\mu}}^{r+1},{\bm{z}}^{r+1})-\bar{h}({\bm{x}}(\bar{}{\bm{y}}^{r+1},{\bm{z}}^{r+1});{\bm{\mu}}^{r+1},{\bm{z}}^{r+1})\geq\frac{p+\gamma^{-}}{2}\|{\bm{x}}(\bar{}{\bm{y}}^{r},{\bm{z}}^{r})-{\bm{x}}(\bar{}{\bm{y}}^{r+1},{\bm{z}}^{r+1})\|^{2}.$
(E.42)

Then, by summing (E.41) and (E.42), we obtain

$\displaystyle\bar{h}({\bm{x}}(\bar{}{\bm{y}}^{r+1},{\bm{z}}^{r+1});{\bm{\mu}}^{r+1},{\bm{z}}^{r})-\bar{h}({\bm{x}}(\bar{}{\bm{y}}^{r},{\bm{z}}^{r});{\bm{\mu}}^{r+1},{\bm{z}}^{r})$
$\displaystyle+\bar{h}({\bm{x}}(\bar{}{\bm{y}}^{r},{\bm{z}}^{r});{\bm{\mu}}^{r+1},{\bm{z}}^{r+1})-\bar{h}({\bm{x}}(\bar{}{\bm{y}}^{r+1},{\bm{z}}^{r+1});{\bm{\mu}}^{r+1},{\bm{z}}^{r+1})$
$\displaystyle\geq(p+\gamma^{-})\|{\bm{x}}(\bar{}{\bm{y}}^{r},{\bm{z}}^{r})-{\bm{x}}(\bar{}{\bm{y}}^{r+1},{\bm{z}}^{r+1})\|^{2}.$
(E.43)

On the other hand, according to the definition of $\bar{h}$, we have

$\displaystyle\bar{h}({\bm{x}}(\bar{}{\bm{y}}^{r},{\bm{z}}^{r});{\bm{\mu}}^{r+1},{\bm{z}}^{r})-\bar{h}({\bm{x}}(\bar{}{\bm{y}}^{r},{\bm{z}}^{r});{\bm{\mu}}^{r+1},{\bm{z}}^{r+1})$
$\displaystyle=\frac{p}{2}\|{\bm{x}}(\bar{}{\bm{y}}^{r},{\bm{z}}^{r})-{\bm{z}}^{r}\|^{2}-\frac{p}{2}\|{\bm{z}}^{r}\|^{2}-\bigg{(}\frac{p}{2}\|{\bm{x}}(\bar{}{\bm{y}}^{r},{\bm{z}}^{r})-{\bm{z}}^{r+1}\|^{2}-\frac{p}{2}\|{\bm{z}}^{r+1}\|^{2}\bigg{)}$
$\displaystyle=p{\bm{x}}(\bar{}{\bm{y}}^{r},{\bm{z}}^{r})^{\top}({\bm{z}}^{r+1}-{\bm{z}}^{r}),$
(E.44)

and

$\displaystyle\bar{h}({\bm{x}}(\bar{}{\bm{y}}^{r+1},{\bm{z}}^{r+1});{\bm{\mu}}^{r+1},{\bm{z}}^{r})-\bar{h}({\bm{x}}(\bar{}{\bm{y}}^{r+1},{\bm{z}}^{r+1});{\bm{\mu}}^{r+1},{\bm{z}}^{r+1})=p{\bm{x}}(\bar{}{\bm{y}}^{r+1},{\bm{z}}^{r+1})^{\top}({\bm{z}}^{r+1}-{\bm{z}}^{r}).$
(E.45)

By taking the difference between (E.2) and (E.45), and substituting it into
(E.2), we obtain

$\displaystyle(p+\gamma^{-})\|{\bm{x}}(\bar{}{\bm{y}}^{r},{\bm{z}}^{r})-{\bm{x}}(\bar{}{\bm{y}}^{r+1},{\bm{z}}^{r+1})\|^{2}$
$\displaystyle\leq
p({\bm{x}}(\bar{}{\bm{y}}^{r},{\bm{z}}^{r})-{\bm{x}}(\bar{}{\bm{y}}^{r+1},{\bm{z}}^{r+1}))^{\top}({\bm{z}}^{r}-{\bm{z}}^{r+1})$
$\displaystyle\leq
p\|{\bm{x}}(\bar{}{\bm{y}}^{r},{\bm{z}}^{r})-{\bm{x}}(\bar{}{\bm{y}}^{r+1},{\bm{z}}^{r+1})\|\|{\bm{z}}^{r}-{\bm{z}}^{r+1}\|.$
(E.46)

Dividing
$(p+\gamma^{-})\|{\bm{x}}(\bar{}{\bm{y}}^{r},{\bm{z}}^{r})-{\bm{x}}(\bar{}{\bm{y}}^{r+1},{\bm{z}}^{r+1})\|$
from both sides of the above inequality, we have

$\displaystyle\|{\bm{x}}(\bar{}{\bm{y}}^{r},{\bm{z}}^{r})-{\bm{x}}(\bar{}{\bm{y}}^{r+1},{\bm{z}}^{r+1})\|\leq\frac{p}{p+\gamma^{-}}\|{\bm{z}}^{r}-{\bm{z}}^{r+1}\|$
(E.47)

which is (50). $\blacksquare$

##  Appendix F Proof of Lemma 5

To prove the dual error bounds, we recap the Hoffman error bound [52, Lemma
3.2.3] given below.

###### Lemma 13

Let ${\bm{A}}\in\mathbb{R}^{m\times n}$, ${\bm{C}}\in\mathbb{R}^{k\times n}$,
${\bm{b}}\in\mathbb{R}^{m}$, ${\bm{d}}\in\mathbb{R}^{k}$, and the linear set

$\displaystyle\mathcal{S}=\\{{\bm{x}}|{\bm{A}}{\bm{x}}\leq{\bm{b}},{\bm{C}}{\bm{x}}={\bm{d}}\\}.$

The distance from a point $\bar{{\bm{x}}}\in\mathbb{R}^{n}$ to $\mathcal{S}$
is bounded by:

$\displaystyle\mathrm{dist}^{2}(\bar{{\bm{x}}},\mathcal{S})\leq\theta^{2}(\|({\bm{A}}\bar{{\bm{x}}}-{\bm{b}})_{+}\|^{2}+\|{\bm{C}}\bar{{\bm{x}}}-{\bm{d}}\|^{2}),$

where $(\cdot)_{+}=\max\\{\cdot,0\\}$ denotes the projection to the
nonnegative orthant and $\theta$ is a positive constant depending on
${\bm{A}}$ and ${\bm{C}}$ only.

To prove Lemma 5, let us define
$F_{i}({\bm{x}}_{i})=f_{i}({\bm{x}}_{i})+\frac{p}{2}\|{\bm{x}}_{i}-{\bm{z}}_{i}^{r}\|^{2}$
for all $i\in[N]$, and thus $\nabla F_{i}({\bm{x}}_{i})=\nabla
f_{i}({\bm{x}}_{i})+p({\bm{x}}_{i}-{\bm{z}}^{r}_{i})$. Note that function
$F_{i}$ is strongly convex with modulus $(p+\gamma^{-})$ and its gradient
$\nabla F_{i}$ is $(p+\gamma^{+})$-Lipschitz continuous. Then by[53, Eqn.
(2.1.8)]), we have

$\displaystyle\|\nabla F_{i}({\bm{x}}_{i})-\nabla
F_{i}({\bm{x}}^{\prime}_{i})\|^{2}\leq(p+\gamma^{+}){\langle}\nabla
F_{i}({\bm{x}}_{i})-\nabla
F_{i}({\bm{x}}_{i}^{\prime}),{\bm{x}}_{i}-{\bm{x}}_{i}^{\prime}{\rangle}.$
(F.48)

Besides, since function $\phi_{i}$ in (16) is concave over ${\bm{y}}$ and
smooth with Lipschitz constant $\textstyle\frac{1}{p+\gamma^{-}}$, we have

$\displaystyle\|\nabla\phi_{i}({\bm{x}}_{i})-\nabla\phi_{i}({\bm{x}}^{\prime}_{i})\|^{2}\leq-\frac{1}{p+\gamma^{-}}{\langle}\nabla\phi_{i}({\bm{x}}_{i})-\nabla\phi_{i}({\bm{x}}_{i}^{\prime}),{\bm{x}}_{i}-{\bm{x}}_{i}^{\prime}{\rangle}.$
(F.49)

Recall (42) and let

$\displaystyle\mathcal{Y}({\bm{\mu}}^{r+1},{\bm{z}}^{r})$
$\displaystyle=\arg\max_{{\bm{y}}}~{}\sum_{i=1}^{N}\bigg{(}\phi_{i}({\bm{y}}_{i},{\bm{z}}_{i}^{r})-\frac{1}{N}{\bm{y}}_{i}^{\top}{\bm{q}}\bigg{)}-{\langle}{\bm{\mu}}^{r+1},{\bm{A}}{\bm{y}}{\rangle}-\frac{\rho}{2}\|{\bm{A}}{\bm{y}}\|^{2}.$
(F.50)

For $\bar{}{\bm{y}}^{r}\in\mathcal{Y}({\bm{\mu}}^{r+1},{\bm{z}}^{r})$, we have
the following KKT conditions

$\displaystyle\begin{bmatrix}\nabla\phi_{1}(\bar{}{\bm{y}}_{1}^{r},{\bm{z}}_{1}^{r})-\frac{{\bm{q}}}{N}\\\
\vdots\\\
\nabla\phi_{N}(\bar{}{\bm{y}}_{N}^{r},{\bm{z}}_{N}^{r})-\frac{{\bm{q}}}{N}\end{bmatrix}-{\bm{A}}^{\top}{\bm{\mu}}^{r+1}-\rho{\bm{L}}^{-}\bar{}{\bm{y}}^{r}={\bm{0}},$
(F.51) $\displaystyle\nabla
F_{i}({\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r}))+{\bm{B}}_{i}^{\top}\bar{}{\bm{y}}_{i}^{r}={\bm{0}},\forall
i\in[N],$ (F.52)

where
$\nabla\phi_{i}(\bar{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r})={\bm{B}}_{i}{\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r})$
for all $i$ obtained by the Danskin’s Theorem [14, Proposition B.22].
Analogously, from (19b) the KKT conditions for ${\bm{y}}^{r+1}$ are

$\displaystyle\begin{bmatrix}\nabla\phi_{1}({\bm{y}}_{1}^{r+1},{\bm{z}}_{1}^{r})-\frac{{\bm{q}}}{N}\\\
\vdots\\\
\nabla\phi_{N}({\bm{y}}_{N}^{r+1},{\bm{z}}_{N}^{r})-\frac{{\bm{q}}}{N}\end{bmatrix}-{\bm{A}}^{\top}{\bm{\mu}}^{r+1}-\rho{\bm{L}}^{-}{\bm{y}}^{r+1}-\rho{\bm{L}}^{+}({\bm{y}}^{r+1}-{\bm{y}}^{r})={\bm{0}},$
(F.53) $\displaystyle\nabla
F_{i}({\bm{x}}_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r}))+{\bm{B}}_{i}^{\top}{\bm{y}}_{i}^{r+1}={\bm{0}},\forall
i\in[N],$ (F.54)

where
$\nabla\phi_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r})={\bm{B}}_{i}{\bm{x}}_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r})$
for all $i$.

Let $\bar{}{\bm{y}}^{r}$ be the projection point of ${\bm{y}}^{r+1}$ onto
$\mathcal{Y}({\bm{\mu}}^{r+1},{\bm{z}}^{r})$, i.e., ${\rm
dist}({\bm{y}}^{r+1},\mathcal{Y}({\bm{\mu}}^{r+1},{\bm{z}}^{r}))=\|{\bm{y}}^{r+1}-\bar{}{\bm{y}}^{r}\|$.
With fixed
$\\{{\bm{x}}_{i}\big{(}\bar{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r}\big{)}\\}_{i=1}^{N}$,
(F.51)-(F.52) forms a linear system of variable $\rho\bar{}{\bm{y}}^{r}$, and
any $\rho\bar{}{\bm{y}}^{r}$ satisfying (F.51)-(F.52) is a point in
$\mathcal{Y}({\bm{\mu}}^{r+1},{\bm{z}}^{r})$. By (F.53)-(F.54) and the Hoffman
bound in Lemma 13, we can bound ${\rm
dist}({\bm{y}}^{r+1},\mathcal{Y}({\bm{\mu}}^{r+1},{\bm{z}}^{r}))$ for some
constant $\theta_{1}>0$ as follows

$\displaystyle\rho^{2}{\rm
dist}^{2}({\bm{y}}^{r+1},\mathcal{Y}({\bm{\mu}}^{r+1},{\bm{z}}^{r}))$
$\displaystyle\leq{\theta_{1}^{2}}\sum_{i=1}^{N}\|{\bm{B}}_{i}({\bm{x}}_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r}))-{\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r}))\|^{2}+{\theta_{1}^{2}}\rho^{2}\sum_{i=1}^{N}\|\nabla
F_{i}({\bm{x}}_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r}))-\nabla
F_{i}({\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r}))\|^{2}$
$\displaystyle~{}~{}~{}+{\theta_{1}^{2}}\rho^{2}\|{\bm{L}}^{+}({\bm{y}}^{r+1}-{\bm{y}}^{r})\|^{2}$
$\displaystyle\overset{\mathrm{(i)}}{\leq}{-\frac{\theta_{1}^{2}}{p+\gamma^{-}}}\sum_{i=1}^{N}{\langle}{\bm{B}}_{i}{\bm{x}}_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r})-{\bm{B}}_{i}{\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r}),{\bm{y}}_{i}^{r+1}-\bar{}{\bm{y}}_{i}^{r}{\rangle}$
$\displaystyle~{}~{}~{}+\theta_{1}^{2}\rho^{2}(p+\gamma^{+})\sum_{i=1}^{N}\left\langle\nabla
F_{i}({\bm{x}}_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r}))\right.\left.-\nabla
F_{i}({\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r})),{\bm{x}}_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r})-{\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r})\right\rangle$
$\displaystyle~{}~{}~{}+\theta_{1}^{2}\rho^{2}\|{\bm{L}}^{+}({\bm{y}}^{r+1}-{\bm{y}}^{r})\|^{2}$
$\displaystyle\overset{\mathrm{(ii)}}{=}-\frac{\theta_{1}^{2}}{p+\gamma^{-}}\sum_{i=1}^{N}{\langle}{\bm{B}}_{i}{\bm{x}}_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r})-{\bm{B}}_{i}{\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r}),{\bm{y}}_{i}^{r+1}-\bar{}{\bm{y}}_{i}^{r}{\rangle}$
$\displaystyle~{}~{}~{}-\theta_{1}^{2}\rho^{2}(p+\gamma^{+})\sum_{i=1}^{N}{\langle}{\bm{B}}_{i}^{T}{\bm{y}}_{i}^{r+1}-{\bm{B}}_{i}^{T}\bar{}{\bm{y}}_{i}^{r},{\bm{x}}_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r})-{\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r}){\rangle}$
$\displaystyle~{}~{}~{}+\theta_{1}^{2}\rho^{2}\|{\bm{L}}^{+}({\bm{y}}^{r+1}-{\bm{y}}^{r})\|^{2},$
(F.55)

where the first term of $\mathrm{(i)}$ is due to (F.49), the third term of
$\mathrm{(i)}$ is directly obtained by (F.48), and $\mathrm{(ii)}$ is owing to
(F.52) and (F.54).

Rearranging the right hand side (RHS) of (F), we further have the upper bound
as

$\displaystyle\rho^{2}{\rm
dist}^{2}({\bm{y}}^{r+1},\mathcal{Y}({\bm{\mu}}^{r+1},{\bm{z}}^{r}))$
$\displaystyle\leq-\bigg{(}\frac{\theta_{1}^{2}}{p+\gamma^{-}}+\theta_{1}^{2}\rho^{2}(p+\gamma^{+})\bigg{)}\sum_{i=1}^{N}{\langle}{\bm{B}}_{i}{\bm{x}}_{i}({\bm{y}}_{i}^{r+1},{\bm{z}}_{i}^{r})-{\bm{B}}_{i}{\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r}),{\bm{y}}_{i}^{r+1}-\bar{}{\bm{y}}_{i}^{r}{\rangle}$
$\displaystyle~{}~{}~{}+\theta_{1}^{2}\rho^{2}\|{\bm{L}}^{+}({\bm{y}}^{r+1}-{\bm{y}}^{r})\|^{2}$
$\displaystyle\overset{\mathrm{(i)}}{=}-\bigg{(}\frac{\theta_{1}^{2}}{p+\gamma^{-}}+\theta_{1}^{2}\rho^{2}(p+\gamma^{+})\bigg{)}\sum_{i=1}^{N}{\langle}\nabla\phi_{i}({\bm{y}}^{r+1},{\bm{z}}^{r})-\nabla\phi_{i}(\bar{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r}),{\bm{y}}_{i}^{r+1}-\bar{}{\bm{y}}_{i}^{r}{\rangle}+\theta_{1}^{2}\rho^{2}\|{\bm{L}}^{+}({\bm{y}}^{r+1}-{\bm{y}}^{r})\|^{2}$
$\displaystyle\overset{\mathrm{(ii)}}{=}-\bigg{(}\frac{\theta_{1}^{2}}{p+\gamma^{-}}+\theta_{1}^{2}\rho^{2}(p+\gamma^{+})\bigg{)}{\langle}\nabla{\mathcal{L}}_{\rho}({\bm{y}}^{r+1},\mu^{r+1};{\bm{z}}^{r})-\nabla{\mathcal{L}}_{\rho}(\bar{}{\bm{y}}^{r},\mu^{r+1};{\bm{z}}^{r})+\rho{\bm{L}}^{-}({\bm{y}}^{r+1}-{\bm{y}}({\bm{\mu}}^{r+1},{\bm{z}}^{r})),$
$\displaystyle~{}~{}~{}~{}~{}~{}{\bm{y}}^{r+1}-\bar{}{\bm{y}}^{r}{\rangle}+\theta_{1}^{2}\rho^{2}\|{\bm{L}}^{+}({\bm{y}}^{r+1}-{\bm{y}}^{r})\|^{2}$
$\displaystyle=-\bigg{(}\frac{\theta_{1}^{2}}{p+\gamma^{-}}+\theta_{1}^{2}\rho^{2}(p+\gamma^{+})\bigg{)}{\langle}\nabla{\mathcal{L}}_{\rho}({\bm{y}}^{r+1},\mu^{r+1};{\bm{z}}^{r})-\nabla{\mathcal{L}}_{\rho}(\bar{}{\bm{y}}^{r},\mu^{r+1};{\bm{z}}^{r}),{\bm{y}}^{r+1}-\bar{}{\bm{y}}^{r}{\rangle}$
$\displaystyle~{}~{}~{}-\bigg{(}\frac{\theta_{1}^{2}}{p+\gamma^{-}}+\theta_{1}\rho^{2}(p+\gamma^{+})\bigg{)}\rho\|{\bm{y}}^{r+1}-\bar{}{\bm{y}}_{i}^{r}\|_{{\bm{L}}^{-}}^{2}+\theta_{1}^{2}\rho^{2}\|{\bm{L}}^{+}({\bm{y}}^{r+1}-{\bm{y}}^{r})\|^{2},$
(F.56)

where $\mathrm{(i)}$ is due to the definition of $\nabla\phi_{i}$ , and
$\mathrm{(ii)}$ is obtained by substituting the gradient of
${\mathcal{L}}_{\rho}$ in (3.2).

On the other hand, by the optimality conditions of (42) for
$\bar{}{\bm{y}}^{r}$ and $\bar{}{\bm{y}}^{r+1}$ respectively, we have

$\displaystyle\nabla_{\bm{y}}{\mathcal{L}}_{\rho}(\bar{}{\bm{y}}^{r},{\bm{\mu}}^{r+1};{\bm{z}}^{r})={\bm{0}},$
(F.57)
$\displaystyle\nabla_{\bm{y}}{\mathcal{L}}_{\rho}({\bm{y}}^{r+1},{\bm{\mu}}^{r+1};{\bm{z}}^{r})-\rho{\bm{L}}^{+}({\bm{y}}^{r+1}-{\bm{y}}^{r})={\bm{0}}.$
(F.58)

Thus,

$\displaystyle{\langle}{\bm{y}}^{r+1}-\bar{}{\bm{y}}^{r},\nabla_{\bm{y}}{\mathcal{L}}_{\rho}(\bar{}{\bm{y}}^{r},{\bm{\mu}}^{r+1};{\bm{z}}^{r})-\nabla_{\bm{y}}{\mathcal{L}}_{\rho}(\bar{}{\bm{y}}^{r};{\bm{z}}^{r}){\rangle}={\langle}{\bm{y}}^{r+1}-\bar{}{\bm{y}}^{r},-\rho{\bm{L}}^{+}({\bm{y}}^{r+1}-{\bm{y}}^{r}){\rangle}.$
(F.59)

Substituting (F.59) into (F.56), we have

$\displaystyle\rho^{2}\|{\bm{y}}^{r+1}-\bar{}{\bm{y}}^{r}\|^{2}$
$\displaystyle\leq-\bigg{(}\frac{\theta_{1}^{2}}{p+\gamma^{-}}+\theta_{1}^{2}\rho^{2}(p+\gamma^{+})\bigg{)}{\langle}{\bm{y}}^{r+1}-\bar{}{\bm{y}}^{r},-\rho{\bm{L}}^{+}({\bm{y}}^{r+1}-{\bm{y}}^{r}){\rangle}$
$\displaystyle~{}~{}~{}-\bigg{(}\frac{\theta_{1}^{2}}{p+\gamma^{-}}+\theta_{1}^{2}\rho^{2}(p+\gamma^{+})\bigg{)}\rho\|{\bm{y}}^{r+1}-\bar{}{\bm{y}}^{r}\|_{{\bm{L}}^{-}}^{2}+\theta_{1}^{2}\rho^{2}\|{\bm{L}}^{+}({\bm{y}}^{r+1}-{\bm{y}}^{r})\|^{2}$
$\displaystyle\leq\frac{\theta_{1}^{2}}{2\eta_{1}}\bigg{(}\frac{1}{p+\gamma^{-}}+\rho^{2}(p+\gamma^{+})\bigg{)}\|{\bm{y}}^{r+1}-\bar{}{\bm{y}}^{r}\|^{2}$
$\displaystyle~{}~{}~{}\\!+\\!\bigg{[}\frac{\eta_{1}}{2}\bigg{(}\frac{1}{p+\gamma^{-}}\\!+\\!\rho^{2}(p+\gamma^{+})\bigg{)}\\!+\\!1\bigg{]}\theta_{1}^{2}\rho^{2}\|{\bm{L}}^{+}({\bm{y}}^{r+1}-{\bm{y}}^{r})\|^{2}$
$\displaystyle~{}~{}~{}-\bigg{(}\frac{\theta_{1}^{2}}{p+\gamma^{-}}+\theta_{1}^{2}\rho^{2}(p+\gamma^{+})\bigg{)}\rho\|{\bm{y}}^{r+1}-\bar{}{\bm{y}}^{r}\|_{{\bm{L}}^{-}}^{2},$
(F.60)

where the last inequality is by Young’s inequality and $\eta_{1}>0$.

Choose
$\eta_{1}=\textstyle\frac{\theta_{1}^{2}}{\rho^{2}}\big{(}\frac{1}{p+\gamma^{-}}+\rho^{2}(p+\gamma^{+})\big{)}$,
then
$\frac{\theta_{1}^{2}}{2\eta_{1}}\big{(}\frac{1}{p+\gamma^{-}}+\rho^{2}(p+\gamma^{+})\big{)}=\frac{\rho^{2}}{2}<\rho^{2}$.
Then, we can obtain (51) from (F).

In addition, by substituting the chosen $\eta_{1}$ into (F), we further obtain
the desired result in (52). $\blacksquare$

##  Appendix G Proof of Lemma 6

For any $\bar{}{\bm{y}}^{r}\in\mathcal{Y}({\bm{\mu}}^{r+1},{\bm{z}}^{r})$, we
have the KKT condition in (F.51) and (F.52).

Analogously, recall from (42) that

$\displaystyle\mathcal{Y}({\bm{\mu}}^{r+1},{\bm{z}}^{r+1})$
$\displaystyle=\arg\max_{{\bm{y}}}~{}\sum_{i=1}^{N}\bigg{(}\phi_{i}({\bm{y}}_{i},{\bm{z}}_{i}^{r+1})-{\bm{y}}_{i}^{\top}\frac{{\bm{q}}}{N}\bigg{)}-{\langle}{\bm{\mu}}^{r+1},{\bm{A}}{\bm{y}}{\rangle}-\frac{\rho}{2}\|{\bm{A}}{\bm{y}}\|^{2}.$
(G.61)

For $\bar{}{\bm{y}}^{r+1}\in\mathcal{Y}({\bm{\mu}}^{r+1},{\bm{z}}^{r+1})$, we
have the following KKT condition

$\displaystyle\begin{bmatrix}{\bm{B}}_{1}{\bm{x}}_{1}(\bar{}{\bm{y}}_{1}^{r+1},{\bm{z}}_{1}^{r+1})-\frac{{\bm{q}}}{N}\\\
\vdots\\\
{\bm{B}}_{N}{\bm{x}}_{N}(\bar{}{\bm{y}}_{N}^{r+1},{\bm{z}}_{N}^{r+1})-\frac{{\bm{q}}}{N}\end{bmatrix}-{\bm{A}}^{\top}{\bm{\mu}}^{r+1}-\rho{\bm{L}}^{-}\bar{}{\bm{y}}^{r+1}={\bm{0}},$
(G.62) $\displaystyle\nabla
F_{i}({\bm{x}}_{i}(\bar{}{\bm{y}}^{r+1},{\bm{z}}_{i}^{r+1}))+{\bm{B}}_{i}^{\top}\bar{}{\bm{y}}^{r+1}={\bm{0}},\forall
i\in[N].$ (G.63)

With fixed
$\\{{\bm{x}}_{i}\big{(}\bar{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r}\big{)}\\}_{i=1}^{N}$,
(F.51) and (F.52) form a linear system of variable $\rho\bar{}{\bm{y}}^{r}$
and ${\bm{\mu}}^{r+1}$. Then, by applying the Hoffman bound in Lemma 13, the
distance from $\bar{}{\bm{y}}^{r+1}$ satisfying (G.62) and (G.63) to the set
$\mathcal{Y}({\bm{\mu}}^{r+1},{\bm{z}}^{r})$ can be bounded as

$\displaystyle\mathrm{dist}^{2}(\bar{}{\bm{y}}^{r},\mathcal{Y}({\bm{\mu}}^{r+1},{\bm{z}}^{r+1}))$
$\displaystyle\leq\theta_{2}^{2}(\|\nabla_{{\bm{x}}}F(\bar{}{\bm{y}}^{r},{\bm{z}}^{r}))-\nabla_{{\bm{x}}}F({\bm{x}}(\bar{}{\bm{y}}^{r+1},{\bm{z}}^{r+1}))\|^{2}+\frac{1}{\rho^{2}}\|\tilde{{\bm{B}}}({\bm{x}}(\bar{}{\bm{y}}^{r}-{\bm{x}}(\bar{}{\bm{y}}^{r+1},{\bm{z}}^{r+1}))\|^{2})$
$\displaystyle\leq\theta_{2}^{2}((p+\gamma^{+})^{2}\|{\bm{x}}(\bar{}{\bm{y}}^{r},{\bm{z}}^{r})-{\bm{x}}(\bar{}{\bm{y}}^{r+1},{\bm{z}}^{r+1})\|^{2}+\frac{B_{\max}^{2}}{\rho^{2}}\|{\bm{x}}(\bar{}{\bm{y}}^{r},{\bm{z}}^{r})-{\bm{x}}(\bar{}{\bm{y}}^{r+1},{\bm{z}}^{r+1})\|^{2})$
$\displaystyle\overset{\mathrm{(i)}}{\leq}\underbrace{\theta_{2}^{2}\bigg{(}(p+\gamma^{+})^{2}\sigma_{3}^{2}+\frac{B_{\max}^{2}\sigma_{3}^{2}}{\rho^{2}}\bigg{)}}_{\triangleq
a_{3}}\|{\bm{z}}^{r}-{\bm{z}}^{r+1}\|^{2},$ (G.64)

where
$\tilde{{\bm{B}}}\triangleq\mathrm{diag}\\{{\bm{B}}_{1},\cdots,{\bm{B}}_{N}\\}$,
and $\mathrm{(i)}$ is because of (50). $\blacksquare$

##  Appendix H Proof of lemma 7

For $\bar{}{\bm{y}}^{r}\in\mathcal{Y}({\bm{\mu}}^{r+1},{\bm{z}}^{r})$, recall
the KKT conditions in (F.51) and (F.52).

For ${{\bm{y}}}({\bm{z}}^{r})\in\mathcal{Y}({\bm{z}}^{r})$, recall from (40)
and we have its KKT conditions as

$\displaystyle\begin{bmatrix}{\bm{B}}_{1}{\bm{x}}_{1}({\bm{z}}^{r})-{\bm{q}}/N\\\
\vdots\\\
{\bm{B}}_{N}{\bm{x}}_{N}({\bm{z}}^{r})-{\bm{q}}/N\end{bmatrix}-{\bm{A}}^{\top}{\bm{\mu}}({\bm{z}}^{r})-\rho{\bm{L}}^{-}{{\bm{y}}}({\bm{z}}^{r})={\bm{0}},$
(H.65) $\displaystyle\nabla
f_{i}({\bm{x}}_{i}({\bm{z}}^{r}))+p({\bm{x}}_{i}({\bm{z}}^{r})-{\bm{z}}_{i}^{r})+{\bm{B}}_{i}^{\top}{{\bm{y}}}_{i}({\bm{z}}^{r})={\bm{0}},\forall
i\in[N],$ (H.66) $\displaystyle{\bm{A}}{{\bm{y}}}({\bm{z}}^{r})=\mathbf{0}.$
(H.67)

With fixed
$\\{{\bm{x}}_{i}\big{(}\bar{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r}\big{)}\\}_{i=1}^{N}$,
(H.65)-(H.67) forms a linear system of variable $\rho\bar{}{\bm{y}}^{r}$ and
${\bm{\mu}}^{r+1}$. Then, by using the Hoffman bound in Lemma 13, the distance
from ${\bm{y}}({\bm{z}}^{r})$ and ${\bm{\mu}}({\bm{z}}^{r})$ satisfying
(H.65)-(H.67) to the set $\mathcal{Y}({\bm{\mu}}^{r+1},{\bm{z}}^{r})$ can be
bounded as

$\displaystyle\rho^{2}\|\bar{}{\bm{y}}^{r}-{{\bm{y}}}({\bm{z}}^{r})\|^{2}+\|{\bm{\mu}}^{r+1}-{\bm{\mu}}({\bm{z}}^{r})\|^{2}$
$\displaystyle\leq{\theta_{3}^{2}\sum_{i=1}^{N}\|{\bm{B}}_{i}{\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r})-{\bm{B}}_{i}{\bm{x}}_{i}({\bm{z}}^{r})\|^{2}}$
$\displaystyle~{}~{}~{}+\rho^{2}\theta_{3}^{2}\sum_{i=1}^{N}\|\nabla
f_{i}({\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r}))-\nabla
f_{i}({\bm{x}}_{i}({\bm{z}}^{r}))+p[{\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r})-{\bm{x}}_{i}({\bm{z}}^{r})]\|^{2}+\theta_{3}^{2}\|{\bm{A}}\bar{}{\bm{y}}^{r}\|^{2}$
$\displaystyle\overset{\mathrm{(i)}}{\leq}-\bigg{(}\frac{\theta_{3}^{2}}{p+\gamma^{-}}\bigg{)}\sum_{i=1}^{N}{\langle}{\bm{B}}_{i}{\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r})-{\bm{B}}_{i}{\bm{x}}_{i}({\bm{z}}^{r})),\bar{}{\bm{y}}_{i}^{r}-{{\bm{y}}}_{i}({\bm{z}}^{r}){\rangle}$
$\displaystyle~{}~{}~{}+\rho^{2}\theta_{3}^{2}(p+\gamma^{+})\sum_{i=1}^{N}{\langle}\nabla
f_{i}({\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r}))-\nabla
f_{i}({\bm{x}}_{i}({\bm{z}}^{r}))+p[{\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r})-{\bm{x}}_{i}({\bm{z}}^{r})],{\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r})-{\bm{x}}_{i}({\bm{z}}^{r}){\rangle}$
$\displaystyle~{}~{}~{}+\theta_{3}^{2}\|{\bm{A}}\bar{}{\bm{y}}^{r}\|^{2}$
$\displaystyle\overset{\mathrm{(ii)}}{\leq}-\bigg{(}\frac{\theta_{3}^{2}}{p+\gamma^{-}}\bigg{)}\sum_{i=1}^{N}{\langle}{\bm{B}}_{i}{\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r})-{\bm{B}}_{i}{\bm{x}}_{i}({\bm{z}}^{r})),\bar{}{\bm{y}}_{i}^{r}-{{\bm{y}}}_{i}({\bm{z}}^{r}){\rangle}$
$\displaystyle~{}~{}~{}+\rho^{2}\theta_{3}^{2}(p+\gamma^{+})\sum_{i=1}^{N}{\langle}-{\bm{B}}_{i}^{\top}\bar{}{\bm{y}}_{i}^{r}+{\bm{B}}_{i}^{\top}{{\bm{y}}}_{i}({\bm{z}}^{r}),{\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r})-{\bm{x}}_{i}({\bm{z}}^{r}){\rangle}+\theta_{3}^{2}\|{\bm{A}}\bar{}{\bm{y}}^{r}\|^{2}$
$\displaystyle=-\bigg{(}\frac{\theta_{3}^{2}}{p+\gamma^{-}}+\rho^{2}\theta_{3}^{2}(p+\gamma^{+})\bigg{)}\sum_{i=1}^{N}{\langle}{\bm{B}}_{i}{\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r})-{\bm{B}}_{i}{\bm{x}}_{i}({\bm{z}}^{r})),\bar{}{\bm{y}}_{i}^{r}-{{\bm{y}}}_{i}({\bm{z}}^{r}){\rangle}+\theta_{3}^{2}\|{\bm{A}}\bar{}{\bm{y}}^{r}\|^{2}$
(H.68)

where $\mathrm{(i)}$ is obtained by (F.49), $\mathrm{(ii)}$ is due to (F.52)
and (H.66). Further, note that we have

$\displaystyle-\sum_{i=1}^{N}{\langle}{\bm{B}}_{i}{\bm{x}}_{i}(\bar{}{\bm{y}}_{i}^{r},{\bm{z}}_{i}^{r})-{\bm{B}}_{i}{\bm{x}}_{i}({\bm{z}}^{r})),\bar{}{\bm{y}}_{i}^{r}-{{\bm{y}}}_{i}({\bm{z}}^{r}){\rangle}$
$\displaystyle\overset{\mathrm{(i)}}{=}-{\langle}{\bm{A}}^{\top}({\bm{\mu}}^{r+1}-{\bm{\mu}}({\bm{z}}^{r}))+\rho{\bm{L}}^{-}(\bar{}{\bm{y}}^{r}-{{\bm{y}}}({\bm{z}}^{r})),\bar{}{\bm{y}}^{r}-{{\bm{y}}}({\bm{z}}^{r}){\rangle}$
$\displaystyle\overset{\mathrm{(ii)}}{=}-{\langle}{\bm{\mu}}^{r+1}-{\bm{\mu}}({\bm{z}}^{r}),{\bm{A}}\bar{}{\bm{y}}^{r}{\rangle}-\|\bar{}{\bm{y}}^{r}-{{\bm{y}}}({\bm{z}}^{r})\|_{\rho{\bm{L}}^{-}}^{2}$
$\displaystyle\leq-{\langle}{\bm{\mu}}^{r+1}-{\bm{\mu}}({\bm{z}}^{r}),{\bm{A}}\bar{}{\bm{y}}^{r}{\rangle},$
(H.69)

where $\mathrm{(i)}$ can be deduced by taking the difference between (F.51)
and (H.65), and $\mathrm{(ii)}$ is due to (H.67). By inserting (H) into (H),
we have

$\displaystyle\rho^{2}\|\bar{}{\bm{y}}^{r}-{{\bm{y}}}({\bm{z}}^{r})\|^{2}+\|{\bm{\mu}}^{r+1}-{\bm{\mu}}({\bm{z}}^{r})\|^{2}$
$\displaystyle\leq-\bigg{(}\frac{\theta_{3}^{2}}{p+\gamma^{-}}+\rho^{2}\theta_{3}^{2}(p+\gamma^{+})\bigg{)}{\langle}{\bm{A}}\bar{}{\bm{y}}^{r},{\bm{\mu}}^{r+1}-{\bm{\mu}}({\bm{z}}^{r}){\rangle}+\theta_{3}^{2}\|{\bm{A}}\bar{}{\bm{y}}^{r}\|^{2}.$
(H.70)

Rearranging the above inequality, we have

$\displaystyle\|\bar{}{\bm{y}}^{r}-{{\bm{y}}}({\bm{z}}^{r})\|^{2}$
$\displaystyle\leq\frac{\theta_{3}^{2}}{\rho^{2}}\|{\bm{A}}\bar{}{\bm{y}}^{r}\|^{2}-\frac{1}{\rho^{2}}\|{\bm{\mu}}^{r+1}-{\bm{\mu}}({\bm{z}}^{r})\|^{2}-\bigg{(}\frac{\theta_{3}^{2}}{\rho^{2}(p+\gamma^{-})}+\theta_{3}^{2}(p+\gamma^{+})\bigg{)}{\langle}{\bm{A}}\bar{}{\bm{y}}^{r},{\bm{\mu}}^{r+1}-{\bm{\mu}}({\bm{z}}^{r}){\rangle},$
$\displaystyle\overset{\mathrm{(i)}}{\leq}\frac{\theta_{3}^{2}}{\rho^{2}}\|{\bm{A}}\bar{}{\bm{y}}^{r}\|^{2}-\frac{1}{\rho^{2}}\|{\bm{\mu}}^{r+1}-{\bm{\mu}}({\bm{z}}^{r})\|^{2}+\bigg{(}\frac{\theta_{3}^{2}}{\rho^{2}(p+\gamma^{-})}+\theta_{3}^{2}(p+\gamma^{+})\bigg{)}\|{\bm{A}}\bar{}{\bm{y}}^{r}\|\|{\bm{\mu}}^{r+1}-{\bm{\mu}}({\bm{z}}^{r})\|,$
$\displaystyle\overset{\mathrm{(ii)}}{\leq}\frac{\theta_{3}^{2}}{\rho^{2}}\|{\bm{A}}\bar{}{\bm{y}}^{r}\|^{2}+\frac{\tau}{2}\bigg{(}\frac{\theta_{3}^{2}}{\rho^{2}(p+\gamma^{-})}+\theta_{3}^{2}(p+\gamma^{+})\bigg{)}\|{\bm{A}}\bar{}{\bm{y}}^{r}\|^{2}$
$\displaystyle~{}~{}~{}+\bigg{[}\frac{\bigg{(}\frac{\theta_{3}^{2}}{\rho^{2}(p+\gamma^{-})}+\theta_{3}^{2}(p+\gamma^{+})\bigg{)}}{2\tau}-\frac{1}{\rho^{2}}\bigg{]}\|{\bm{\mu}}^{r+1}-{\bm{\mu}}({\bm{z}}^{r})\|^{2},$
(H.71)

where $\mathrm{(i)}$ is obtained by Cauchy-Schwartz inequality, and
$\mathrm{(ii)}$ is due to Young’s inequality with $\tau>0$.

If we choose
$\textstyle\tau=\frac{\theta_{3}^{2}}{p+\gamma^{-}}+\rho^{2}\theta_{3}^{2}(p+\gamma^{+})$,
then
$\textstyle\frac{1}{2\tau}\bigg{(}\frac{\theta_{3}^{2}}{\rho^{2}(p+\gamma^{-})}+\theta_{3}^{2}(p+\gamma^{+})\bigg{)}-\frac{1}{\rho^{2}}<0$.
By substituting this choice of $\tau$ into (H.71), we obtain the desired (57).
$\blacksquare$

##  Appendix I Proof of Lemma 8(b)

According to Lemma 8(a), the upper bound for $\Phi^{r+1}-\Phi^{r}$ is non-
positive. Therefore we can have the conclusion that the potential function
$\Phi$ is non-increasing with $r$.

Recall from (3.2) that:

$\displaystyle
G({\bm{x}}^{r+1},{\bm{y}}^{r+1},{\bm{\mu}}^{r+1},{\bm{z}}^{r+1})$
$\displaystyle=\sum_{i=1}^{N}\bigg{(}f_{i}({\bm{x}}_{i}^{r+1})+\frac{p}{2}\|{\bm{x}}_{i}^{r+1}-{\bm{z}}_{i}^{r+1}\|^{2}+({\bm{y}}_{i}^{r+1})^{\top}{\bm{B}}_{i}{\bm{x}}_{i}^{r+1}-\frac{1}{N}({\bm{y}}_{i}^{r+1})^{\top}{\bm{b}}\bigg{)}-({\bm{\mu}}^{r+1})^{\top}{\bm{A}}{\bm{y}}^{r+1}-\frac{\rho}{2}\|{\bm{A}}{\bm{y}}^{r+1}\|^{2}.$
(I.72)

Note that,

$\displaystyle\sum_{i=1}^{N}\bigg{(}({\bm{y}}_{i}^{r+1})^{\top}{\bm{B}}_{i}{\bm{x}}_{i}^{r+1}-\frac{1}{N}({\bm{y}}_{i}^{r+1})^{\top}{\bm{b}}\bigg{)}-({\bm{\mu}}^{r+1})^{\top}{\bm{A}}{\bm{y}}^{r+1}$
$\displaystyle=\sum_{i=1}^{N}({\bm{y}}_{i}^{r+1})^{\top}({\bm{B}}_{i}{\bm{x}}_{i}^{r+1}-{\bm{b}}/N-{\bm{A}}_{i}^{\top}{\bm{\mu}}^{r+1})$
$\displaystyle\overset{\mathrm{(i)}}{=}\sum_{i=1}^{N}({\bm{y}}_{i}^{r+1})^{\top}(2\rho|\mathcal{N}_{i}|{\bm{y}}_{i}^{r+1}-\rho{\bm{L}}^{+}_{i}{\bm{y}}^{r})$
$\displaystyle=\rho({\bm{y}}^{r+1})^{\top}(2{\bm{D}}{\bm{y}}^{r+1}-{\bm{L}}^{+}{\bm{y}}^{r})$
$\displaystyle\overset{\mathrm{(ii)}}{=}\rho({\bm{y}}^{r+1})^{\top}({\bm{L}}^{+}({\bm{y}}^{r+1}-{\bm{y}}^{r})+{\bm{L}}^{-}{\bm{y}}^{r+1}).$
$\displaystyle=\rho\|{\bm{y}}^{r+1}\|^{2}_{{\bm{L}}^{-}}+\frac{\rho}{2}\|{\bm{y}}^{r+1}-{\bm{y}}^{r}\|^{2}_{{\bm{L}}^{+}}+\frac{\rho}{2}\|{\bm{y}}^{r+1}\|^{2}_{{\bm{L}}^{+}}-\frac{\rho}{2}\|{\bm{y}}^{r}\|^{2}_{{\bm{L}}^{+}},$
(I.73)

where $\mathrm{(i)}$ is obtained by (23) and $\mathrm{(ii)}$ is due to
${\bm{L}}^{+}=2{\bm{D}}-{\bm{L}}^{-}$.

Thus, recall $\tilde{G}$, then combine the results in (3.2) and (I). We have

$\displaystyle\tilde{G}^{r+1}$ $\displaystyle\triangleq
G({\bm{x}}^{r+1},{\bm{y}}^{r+1},{\bm{\mu}}^{r+1},{\bm{z}}^{r+1})+\frac{\alpha}{2}\|{\bm{y}}^{r+1}\|^{2}_{{\bm{L}}^{-}}+\frac{\rho}{2}\|{\bm{y}}^{r+1}-{\bm{y}}^{r}\|^{2}_{{\bm{L}}^{+}}$
$\displaystyle=\sum_{i=1}^{N}f_{i}({\bm{x}}_{i}^{r+1})+\frac{p}{2}\|{\bm{x}}^{r+1}-{\bm{z}}^{r+1}\|^{2}+\bigg{(}\frac{\alpha}{2}+\frac{\rho}{2}\bigg{)}\|{\bm{y}}^{r+1}\|^{2}_{{\bm{L}}^{-}}$
$\displaystyle~{}~{}~{}+{\rho}\|{\bm{y}}^{r+1}-{\bm{y}}^{r}\|^{2}_{{\bm{L}}^{+}}+\frac{\rho}{2}\|{\bm{y}}^{r+1}\|^{2}_{{\bm{L}}^{+}}-\frac{\rho}{2}\|{\bm{y}}^{r}\|^{2}_{{\bm{L}}^{+}},$
(I.74)

which implies that

$\displaystyle\sum_{r=0}^{T}(\tilde{G}^{r+1}-\underline{f})>-\frac{\rho}{2}\|{\bm{y}}^{0}\|_{{\bm{L}}^{+}}^{2}>-\infty,$
(I.75)

for any $T$. Due to the weak duality,
$d({\bm{\mu}}^{r};{\bm{z}}^{r})\geq{\mathcal{L}}_{\rho}({\bm{y}}^{r},{\bm{\mu}}^{r};{\bm{z}}^{r})$.
We further obtain
$\sum_{r=0}^{T}(\Phi^{r+1}-\underline{f})>-\frac{\rho}{2}\|{\bm{y}}^{0}\|^{2}_{{\bm{L}}^{+}}>-\infty.$
By this fact and the conclusion of (a) that $\Phi^{r+1}$ is non-increasing,
$\Phi^{r+1}$ must be lower bounded, which is denoted as $\underline{\Phi}$.
$\blacksquare$

##  Appendix J Proof of lemma 9

Let $({\bm{x}}({\bm{z}}),{\bm{y}}_{0}({\bm{z}}))$ be the primal-dual solution
of problem (10). They satisfy the following KKT conditions

$\displaystyle\nabla_{{\bm{x}}}f_{i}({\bm{x}}_{i}({\bm{z}}_{i}))+p({\bm{x}}_{i}({\bm{z}}_{i})-{\bm{z}}_{i})+{\bm{B}}_{i}^{\top}{\bm{y}}_{0}({\bm{z}})=0,\forall
i\in[N],$ (J.76a)
$\displaystyle\sum_{i=1}^{N}{\bm{B}}_{i}{\bm{x}}_{i}({\bm{z}}_{i})={\bm{q}}.$
(J.76b)

Suppose that $\|{\bm{z}}-{\bm{x}}({\bm{z}})\|^{2}\leq\epsilon$ for some
${\bm{z}}$. We analyze the gap between ${\bm{z}}$ and a KKT solution of
problem (P) (see (11)). Firstly, consider the following chain of bounds

$\displaystyle\sum_{i=1}^{N}\left\|\nabla_{{\bm{x}}}f_{i}({\bm{z}}_{i})+{\bm{B}}_{i}^{\top}{\bm{y}}_{0}({\bm{z}})\right\|^{2}$
$\displaystyle\leq\sum_{i=1}^{N}\left\|\nabla_{{\bm{x}}}f_{i}({\bm{z}}_{i})\\!-\\!\nabla_{{\bm{x}}}f_{i}({\bm{x}}_{i}({\bm{z}}_{i}))\\!+\\!\nabla_{{\bm{x}}}f_{i}({\bm{x}}_{i}({\bm{z}}_{i}))\\!+\\!{\bm{B}}_{i}^{\top}{\bm{y}}_{0}({\bm{z}})\right\|^{2}$
$\displaystyle\overset{{\mathrm{(i)}}}{=}\sum_{i=1}^{N}\left\|\nabla_{{\bm{x}}}f_{i}({\bm{z}}_{i})-\nabla_{{\bm{x}}}f_{i}({\bm{x}}_{i}({\bm{z}}_{i}))+p({\bm{z}}_{i}-{\bm{x}}_{i}({\bm{z}}_{i}))\right\|^{2}$
$\displaystyle{\leq}\sum_{i=1}^{N}2\left\|\nabla_{{\bm{x}}}f_{i}({\bm{z}}_{i})-\nabla_{{\bm{x}}}f_{i}({\bm{x}}_{i}({\bm{z}}_{i}))\right\|^{2}$
$\displaystyle~{}~{}~{}~{}+\sum_{i=1}^{N}2\left\|p({\bm{z}}_{i}-{\bm{x}}_{i}({\bm{z}}_{i}))\right\|^{2}$
$\displaystyle\overset{{\mathrm{(ii)}}}{\leq}2(p^{2}+(\gamma^{+})^{2})\sum_{i=1}^{N}\left\|{\bm{z}}_{i}-{\bm{x}}_{i}({\bm{z}}_{i})\right\|^{2}$
$\displaystyle\leq 2(p^{2}+(\gamma^{+})^{2})\epsilon,$ (J.77)

where $\mathrm{(i)}$ is by (J.76a), and $\mathrm{(ii)}$ is by Assumption 2.

Secondly, we have

$\displaystyle\left\|\sum_{i=1}^{N}{\bm{B}}_{i}{\bm{z}}_{i}-{\bm{q}}\right\|^{2}$
$\displaystyle=\left\|\sum_{i=1}^{N}{\bm{B}}_{i}{\bm{z}}_{i}-\sum_{i=1}^{N}{\bm{B}}_{i}{\bm{x}}_{i}({\bm{z}}_{i})+\sum_{i=1}^{N}{\bm{B}}_{i}{\bm{x}}_{i}({\bm{z}}_{i})-{\bm{q}}\right\|^{2}$
$\displaystyle\overset{{\mathrm{(i)}}}{=}\left\|\sum_{i=1}^{N}{\bm{B}}_{i}({\bm{z}}_{i}-{\bm{x}}_{i}({\bm{z}}_{i}))\right\|^{2}$
$\displaystyle\overset{{\mathrm{(ii)}}}{\leq}NB_{\max}^{2}\sum_{i=1}^{N}\|{\bm{z}}_{i}-{\bm{x}}_{i}({\bm{z}}_{i})\|^{2}$
$\displaystyle\leq NB_{\max}^{2}\epsilon,$ (J.78)

where ${\mathrm{(i)}}$ is by(J.76b) and ${\mathrm{(ii)}}$ is by
$B_{\max}\triangleq\max_{i\in[N]}\|{\bm{B}}_{i}\|$. Thus, we can conclude that
${\bm{z}}$ is a $\kappa\epsilon$-KKT solution of problem (P) with

$\displaystyle\kappa\triangleq\max\\{2(p^{2}+(\gamma^{+})^{2},NB_{\max}^{2}\\}.$

$\blacksquare$

##  Appendix K Proof of Corollary 1

To determine the dependence of the graph topology in the convergence rate, we
make the analysis by 3 steps, including (1) simplifying the upper bound of the
term $\|{\bm{z}}^{t}-{\bm{x}}({\bm{z}}^{t})\|$, (2) deriving the bound for
$\theta_{1},\theta_{2},\theta_{3}$ which appear in Lemmas 5, 6 and 7, and (3)
determining the dependence of the bound on the topology and data matrix.

  * (1)

Simplify the upper bound of $\|{\bm{z}}^{t}-{\bm{x}}({\bm{z}}^{t})\|$: Recall
in (4.4) that

$\displaystyle\|{\bm{z}}^{t}-{\bm{x}}({\bm{z}}^{t})\|$
$\displaystyle\leq\|{\bm{z}}^{t}-{\bm{x}}^{t+1}\|+\|{\bm{x}}^{t+1}-{\bm{x}}({\bm{z}}^{t})\|$