\section{Introduction}\label{sec:introduction}
\setcounter{section}{1}

Distributed optimization over networks has been an appealing area of research in recent years. In this regime, a collection of agents (e.g., processing units, sensors) communicate their local information with their neighboring agents to cooperatively optimize a global objective.  It is through this cooperation that learning from massive datasets can be made possible.  Moreover,  the decentralized storage of the data over the network may allow for preserving the privacy of the agents.  It is for these reasons that distributed optimization has found a wide range of applications in wireless sensor networks, machine learning, and signal processing{~\cite{AngeliaReviewDO2018}}.  Despite the significant advances in the design and analysis of the optimization methods over networks, the existing models and algorithms are still less satisfactory in some regimes than those of centralized optimization. For example, there is still much left to be understood about how to tackle the presence of nonlinearity and uncertainty in the functional constraints, while requiring a low number of communications and enforcing weak assumptions on the network topology. The goal in this work is to tackle some of the shortcomings in distributed constrained  optimization through considering a new unifying mathematical framework described as follows. Consider a system with $m$ agents where the $i^{\text{th}}$ agent is associated with a component function $f_i:\mathbb{R}^n \to \mathbb{R}$ and a mapping $F_i: \mathbb{R}^n \to \mathbb{R}^{n}$. Our goal is to solve the following distributed constrained optimization problem
\begin{alignat}{2}\label{prob:initial_problem}\tag{P}
&{\text{minimize }} \quad \sum\nolimits_{i=1}^m f_i(x) \\
 &\text{subject to} \quad x \in \text{SOL}\left(X,\sum\nolimits_{i=1}^mF_i \right), \nonumber
\end{alignat}
\noindent  where $X \subseteq \mathbb{R}^n$ is a set and $\text{SOL}\left(X,\sum\nolimits_{i=1}^mF_i \right)$ denotes the solution set of the variational inequality $\text{VI}\left(X,\sum\nolimits_{i=1}^mF_i \right)$ defined as follows: $x \in X$ solves $ \text{VI}\left(X,\sum\nolimits_{i=1}^mF_i \right)$ if we have $(y-x)^T\sum\nolimits_{i=1}^mF_i(x)\geq 0$ for all $y \in X$. Problem \eqref{prob:initial_problem} represents a distributed optimization framework in the sense that the information about $f_i$ and $F_i$ is locally known to the $i^{\text{th}}$ agent, while the set $X$ is globally known to all the agents. In this work, we consider the case where the local functions $f_i$ are nondifferentiable and merely convex, and mappings $F_i$ are single-valued, continuous, and merely monotone. Model \eqref{prob:initial_problem} captures the canonical formulation of distributed optimization  
\begin{alignat}{2}\label{prob:canonical_DO_problem} 
&{\text{minimize }} \quad \sum\nolimits_{i=1}^m f_i(x) \\
 &\text{subject to} \quad x \in X, \nonumber
\end{alignat}
that has been extensively studied in the literature. Indeed, by choosing $F_i(x):= 0_n$ for all $i$, model \eqref{prob:initial_problem} is equivalent to model \eqref{prob:canonical_DO_problem}. In what follows, we show how the proposed model offers a framework for extending \eqref{prob:canonical_DO_problem} to capture more challenging distributed constrained optimization problems. 

{\noindent {\bf(i) Distributed optimization problems with complementarity constraints:}} Nonlinear complementarity problems (NCP) have been employed to formulate diverse applications in engineering and economics. The celebrated Wardrop's principle of equilibrium in traffic networks and also, the Walras's law of competitive equilibrium in economics are among important examples that can be represented using NCP (cf.~\cite{FerrisPang97}). Formally, NCP is defined as follows. Given a mapping $F:\mathbb{R}_+^n\to \mathbb{R}^n$, $x \in \mathbb{R}^n$ solves $\text{NCP}(F)$ if $0\leq x \perp F(x) \geq  0$, where $\perp$ denotes the perpendicularity operator between two vectors. It is known that $\text{NCP}(F)$ can be cast as $\text{VI}( \mathbb{R}^n_{+},F)$ (see Proposition 1.1.3 in~\cite{FacchineiPang2003}). In many applications where $F$ is merely monotone, $\text{NCP}(F)$ may admit multiple equilibria. In such cases, one may consider finding the best equilibrium with respect to a global metric $f:\mathbb{R}^n \to\mathbb{R}$. For example in traffic networks, the total travel time of the network users can be considered as the objective $f$. In fact, the problem of computing the best equilibrium of an NCP is important to be addressed particularly in the design of transportation networks where there is a need to estimate the efficiency of the equilibrium \cite{NisanBook2007,JohariThesis}. In this regime, the goal is to minimize $f(x)$ where $x$ solves $\text{NCP}(F)$. Motivated by applications in traffic equilibrium problems, stochastic variants of NCP have been considered more recently~\cite{CZF2009NCP, ZC2008NCP}. Consider a stochastic NCP given by 
\begin{align*}
x \geq 0, \quad \mathbb{E}\left[F(x,\xi(\omega)) \right] \geq 0, \quad x^T \mathbb{E}\left[F(x,\xi(\omega))\right]=0,
\end{align*}
where $\xi:\Omega \to \mathbb{R}^d$ is a random variable associated with the probability space $(\Omega, \mathcal{F},\mathbb{P})$ and $F:\mathbb{R}_+^n \times \Omega \to \mathbb{R}^n$ is a stochastic single-valued mapping. Let $\mathcal{S}_i$ denote a local index set of independent and identically distributed samples from the random variable $\xi$. Employing a sample average approximation scheme, one can consider a distributed NCP given by
\begin{align*}
&x \geq 0, \quad  \textstyle\sum\nolimits_{i=1}^m \textstyle\sum\nolimits_{\ell \in \mathcal{S}_i} F(x,\xi_\ell) \geq 0, \quad x^T \left( \textstyle\sum\nolimits_{i=1}^m \textstyle\sum\nolimits_{\ell \in \mathcal{S}_i} F(x,\xi_\ell)\right)=0.
\end{align*}
Let $f:\mathbb{R}_+^n \times \Omega \to \mathbb{R}^n$ denote a stochastic objective function that measures the performance of a given equilibrium at a realization of $\xi$. Then, the problem of distributed computation of the best equilibrium of the preceding NCP is formulated as
\begin{alignat}{2}\label{prob:dist_best_NCP} 
&{\text{minimize }} \quad \textstyle\sum\nolimits_{i=1}^m \textstyle\sum\nolimits_{\ell \in \mathcal{S}_i} f(x,\xi_\ell) \\
 &\text{subject to} \quad x \in \mbox{SOL}\left(\mathbb{R}^n_+,\textstyle\sum\nolimits_{i=1}^m \textstyle\sum\nolimits_{\ell \in \mathcal{S}_i} F(\bullet,\xi_\ell)\right). \notag
\end{alignat}
\noindent {Importantly,} the proposed model \eqref{prob:initial_problem} captures problem \eqref{prob:dist_best_NCP} by defining $X\triangleq \mathbb{R}^n_+$, $f_i(x) \triangleq\sum\nolimits_{\ell \in \mathcal{S}_i} f(x,\xi_\ell)$, and $F_i(x) \triangleq\sum\nolimits_{\ell \in \mathcal{S}_i} F(x,\xi_\ell) $. In section \ref{sec:numerical_implementation}, we present {preliminary numerical experiments} where we solve problem \eqref{prob:dist_best_NCP} for a given transportation network. 

{\noindent {\bf(ii) Distributed optimization with local nonlinear inequality and linear equality constraints:} Another class of problems that can be captured {by} problem \eqref{prob:initial_problem} is given as
	\begin{equation}\label{prob:subclass_constrained_opt}
		\begin{aligned}
			& {\text{minimize}} && \textstyle\sum\nolimits_{i=1}^{m}f_i(x) \\
			& \text{subject to}  &&g_{i,1}(x) \leq 0\ ,\  \ldots\ , \  g_{i,n_i}(x) \leq 0,&&  {\hbox{for } i \in [m]},\\
		&	&&  A_ix=b_i, &&  {\hbox{for }  i \in [m]}, \\
	
		&	&&  x \in X,
		\end{aligned}
	\end{equation}
	where agent $i$ is associated with function $f_i:\mathbb{R}^n \to\mathbb{R}$, functions $g_{i,j}:\mathbb{R}^n \to \mathbb{R}$ for $j \in {[n_i]}$, and parameters $A_i \in \mathbb{R}^{d_i\times n} $ and $b_i \in \mathbb{R}^{d_i}$.  {The notation $[m]$ is used to abbreviate $\{1,\ldots,m\}$.} The set $X$ is globally known to all the agents while $f_i(x)$, $g_{i,j}(x)$, $A_i$, and $b_i$ are  locally known to agent $i$. In the following, we show that problem \eqref{prob:subclass_constrained_opt} can be represented as model \eqref{prob:initial_problem}.}
	\begin{lemma}\label{Lem:nonlinear_const_DO_problem} {Consider problem \eqref{prob:subclass_constrained_opt}. Let functions $g_{i,j}(x)$ be continuously differentiable and convex  for all {$i \in [m]$} and {$j \in[n_i]$}. Assume that the feasible region of problem \eqref{prob:subclass_constrained_opt} is nonempty and the set $X$ is closed and convex. Then, problem \eqref{prob:subclass_constrained_opt} is equivalent to \eqref{prob:initial_problem} where we define $F_i:\mathbb{R}^n \to \mathbb{R}^n$ as $F_i(x) \triangleq A_i^T(A_ix-b_i)+\textstyle\sum_{j=1}^{n_i}\max\{0,g_{i,j}(x)\}\nabla g_{i,j}(x)$.}
	\end{lemma} 
	\begin{proof}
	See Appendix \ref{app:lem:nonlinear_const_DO_problem}.
\end{proof}

{\noindent {\bf(iii) Distributed optimization with coupling nonlinear equality constraints:}
Another class of problems that can be reformulated as \eqref{prob:canonical_DO_problem} is given as
	\begin{equation}\label{prob:subclass_equ_constrained}
		\begin{aligned}
			& {\text{minimize}} & & 	 				\textstyle\sum_{i=1}^mf_i(x) \\
			& \text{subject to}  & &  \textstyle\sum_{i=1}^mF_i(x) =0_n,
		\end{aligned}
	\end{equation}
	where agent $i$ is associated with a local mapping $F_i:\mathbb{R}^n\to \mathbb{R}^n$ and a local objective $f_i:\mathbb{R}^n \to \mathbb{R}$.  This model finds relevance to the multi-user optimization problems in network resource allocation applications where the objective and constraints are not separable by each user~\cite{KoshalNedichShanbhag2011,TUWA2021,CMWRG2017}. Note that the feasible solution set of problem \eqref{prob:subclass_equ_constrained} is equal to $\text{SOL}(\mathbb{R}^n,\sum_{i=1}^mF_i)$, implying that \eqref{prob:subclass_equ_constrained} is captured by \eqref{prob:initial_problem}.}
	
{\noindent {{\bf  Scope and literature review:}} }
\begin{table*}[t]
	\caption{{Comparison of {incremental gradient} schemes for solving {finite-sum} problems} }
	\centering
	\begin{footnotesize}
	
		\begin{tabular}{ l  l l  l l p{10cm} }
		
			\hline
			Reference  & Method & {Problem class} & Problem formulation & Convergence rate(s)  \\ 
		
			\cmidrule(r){1-5}
			\begin{minipage}{0.4cm} { \centering \cite{nedichThesis} }\end{minipage} &  Projected IG   &    $f_i \in C^{0,0}_{0, L}$   &   {$\min_{x\in X} \ \sum_{i = 1}^m f_i(x)$}  & $\mathcal{O}\left(\frac{1}{\sqrt{k}}\right)$   \\
		
			\cite{BlattHeroGauchman2007, GurbuzbalanOzdaglarParrilo2017}    & IAG &    $f_i \in C_{\mu,L}^{1,1}$   &   {$\min_{x\in \mathbb{R}^n}\sum_{i = 1}^m f_i(x)$}  & ${\cal O} \left(\rho^k\right)$ \\
		
			\cite{DefazioBachJulien2014}    & SAGA& $f_i \in C_{0,L}^{1,1}, \ C_{\mu,L}^{1,1}$   &   {$\min_{x\in X} \ \sum_{i = 1}^m f_i(x)$}   & ${\cal O} \left(\frac{1}{{k}}\right)$, ${\cal O} \left(\rho^k\right)$ \\
		
			\cite{VanliGurbuzbalanOzdaglar2018}  &  Proximal IAG  & { $f_i \in C_{\mu,L}^{1,1}$}  &   {$\min_{x\in X} \sum_{i = 1}^m f_i(x)$}  &${\cal O} \left(\rho^k\right)$  \\
		
			\cite{GurbuzbalanOzdaglarParrilo2019}   & \begin{minipage}{2cm} IG\end{minipage} &  { $f_i \in C_{\mu,L}^{2,1}$}   & \begin{minipage}{3.2cm} {$\min_{x\in \mathbb{R}^n}\sum_{i = 1}^m f_i(x)$} \end{minipage} & ${\cal O} \left(\frac{1}{{k}}\right)$,   ${\cal O} \left(\frac{1}{{k^2}}\right)$   \\
		
			\begin{minipage} {0.9cm} {\cite{Jalilzadeh2020}} \end{minipage}  & \begin{minipage}{2.5cm}   {Primal-Dual IG} \end{minipage} & \begin{minipage}{2.2cm}$f_i \in C^{0,0}_{0, L}$ \end{minipage} & \begin{minipage}{3.2cm}   $\min_{x\in X} \ \sum_{i = 1}^m f_i(x)$ \\   \ $ Ax - b \in -\mathcal{K}$\end{minipage} & \begin{minipage}{3.0cm}  {${\cal O} \left(\frac{1}{k}\right)$}\end{minipage}    \\
		
		
			\cmidrule(r){1-5}
		
			\begin{minipage} {0.9cm} \textbf{This}\\ \textbf{work} \end{minipage}  & \begin{minipage}{2.5cm}   pair-IG \end{minipage} & \begin{minipage}{2.2cm}$f_i \in {C^{0,0}_{0, L}},$\\ $F_i$ is monotone \end{minipage} & \begin{minipage}{3.2cm}   {$\min  \ \sum_{i = 1}^m f_i(x)$ \\  $ x \in \text{SOL}\left(X,\sum_{i=1}^mF_i \right)$} \end{minipage} & \begin{minipage}{3.7cm}  suboptimality: ${\cal O} \left( k^{b-0.5}\right)$ \\ infeasibility: ${\cal O} \left(k^{-b}\right)$\\ { { where $0<b<0.5$}}\end{minipage}    \\
			\hline
		\end{tabular}
	
	\end{footnotesize}
	\label{table:schemes_literature}
\end{table*}
In addressing the proposed formulation \eqref{prob:initial_problem}, our focus in this paper lies in the development of an incremental gradient (IG) method.  IG methods are among popular avenues for addressing the classical distributed optimization model \eqref{prob:canonical_DO_problem} and they have received an increasing attention in recent years in addressing applications arising in sensor networks and machine learning~\cite{nedichThesis, wang2015,GurbuzbalanOzdaglarParrilo2017,GurbuzbalanOzdaglarParrilo2019}. In these schemes, utilizing the additive structure of the problem, the algorithm cycles through the data blocks and updates the local estimates of the optimal solution in a sequential manner \cite{BertsekasNLPBook2016}.  While the first variants of IG schemes find their roots in addressing neural networks as early as in the 1980s \cite{BertsekasIG11}, the complexity analysis of these schemes has been a trending research topic in the fields of control and machine learning in the past two decades. In addressing {the} constrained problems {with easy-to-project constraint sets}, the projected incremental gradient (P-IG) method and its subgradient variant were developed~\cite{NedBert2001}. Considering the smooth case, the P-IG scheme is described as follows. Given an initial point $x_{0,1} \in X$ where $X \subseteq \mathbb{R}^n$ denotes the constraint set, for each $k \geq 0$, consider the update rules given by
 \begin{align*}
 	&x_{k,i+1}: = \mathcal{P}_X\left(x_{k,i}-\gamma_{k}\nabla f_i\left(x_{k,i}\right)\right), \quad {\text{for }i \in [m]},\\
 	&x_{k+1,1}: = x_{k,m+1},
 \end{align*} 
 where ${x_{k,i} \in \mathbb{R}^n}$ denotes agent $i$'s local copy of the decision variables at iteration $k$, $\mathcal{P}$ denotes the Euclidean projection operator defined as $\mathcal{P}_X(z)\triangleq \text{argmin}_{x\in X}\|x-z\|_2$, and $\gamma_k>0$ denotes the stepsize parameter. To motivate our research, we provide an overview of the different variants of existing IG schemes and then, highlight some of the shortcomings of these methods in the constrained regime, in particular, in addressing VI constraints in \eqref{prob:initial_problem}. Recently, under the assumption of strong convexity and twice continuous differentiability of the objective function, and also, boundedness of the generated iterates, the standard IG method was proved to converge with the rate $\mathcal{O}(1/k)$ in the unconstrained case \cite{GurbuzbalanOzdaglarParrilo2019}. This {is an} improvement to the {previously} known rate of $\mathcal{O}(1/\sqrt{k})$ for the merely convex case. Accelerated variants of IG schemes with provable convergence speeds were {also developed}, including  the incremental aggregated gradient method (IAG) \cite{BlattHeroGauchman2007,GurbuzbalanOzdaglarParrilo2017}, SAG \cite{RouxSchmidtBach2012}, and SAGA \cite{DefazioBachJulien2014}. While addressing the merely convex case, SAGA {using} averaging achieves a sublinear convergence rate, assuming strong convexity and smoothness, this is improved for non-averaging variants of SAGA and IAG to a linear rate.  Table \ref{table:schemes_literature} presents a summary of the standard IG schemes in addressing unconstrained and constrained finite-sum problems. As evidenced, most of the past research efforts on the design and analysis of algorithms for distributed constrained optimization problems have focused on addressing easy-to-project sets or sets with {linear} functional inequalities. This has been done through employing duality theory, projection, or penalty methods (see \cite{ChangNedichScaglione2014,AybatHamedani2016,ScutariSun2019, NedichTatarenko2020}).  Also, a celebrated variant of the dual based schemes is the alternating direction method of multipliers (ADMM) (e.g., see \cite{MakhdoumiOzdaglar2017,TangDaoutidis2019}). Other related papers that have utilized duality theory in distributed constrained regimes include \cite{Bertsekas2015,AybatHamedani2016, HamedaniAybat2019}.  Despite the extensive work in the area of constrained optimization, no provably convergent iterative method exits in the literature that can be employed to solve distributed optimization problems with VI constraints. In fact, we are unaware of any IG methods with complexity guarantees that can be employed for addressing any of the three individual subclass problems (i), (ii),  and (iii) mentioned earlier. 



 \noindent {{\bf Contributions:} {Motivated by the aforementioned research gap, we make the following main contributions.}
 
 \noindent  \textit{(i) Complexity guarantees for addressing model \eqref{prob:initial_problem}}: The main contribution of this work lies in the development of a distributed iterative method equipped with agent-specific iteration complexity guarantees for solving distributed optimization problems with VI constraints of the form \eqref{prob:initial_problem}. To this end, employing a regularization-based relaxation technique, we propose a projected averaging iteratively regularized incremental gradient method (pair-IG) presented by Algorithm \ref{alg:IR-IG_avg}. In Theorem \ref{thm:rates}, under merely convexity of the global objective function and merely monotonicity of the global mapping, we derive new non-asymptotic suboptimality and infeasibility convergence rates for each agent's generated iterates. This implies a total iteration complexity of $\mathcal{O}\left((C_f+C_F)^4\epsilon^{-4}\right)$ for obtaining an $\epsilon$-approximate solution where $C_f$ and $C_F$ denote the bounds on the global objective function's subgradients and the global mapping over {a} compact convex set $X$, respectively. Iterative regularization (IR) has been recently employed as a constraint-relaxation technique in a class of bilevel optimization problems \cite{FarzadPushPull2020,FarzadOMS19} and also in regimes where the duality theory may not be directly applied \cite{FarzadMathProg17,KaushikYousefianSIOPT2021}. Of these, in our recent work \cite{KaushikYousefianSIOPT2021} we employed the IR technique to derive a provably convergent method for solving problem \eqref{prob:initial_problem} in a centralized framework, where the information of the objective function is globally known by the agents. Unlike in \cite{KaushikYousefianSIOPT2021}, here we assume {that} the agents have {access} only {to} local information about both the objective function and the mapping. {It is worth emphasizing that this lack of centralized access to information introduces a major challenge in both the design and the complexity analysis of the new algorithmic framework in addressing the distributed model \eqref{prob:initial_problem}.} Motivated by the need for distributed implementations, in a preliminary version of this work, in \cite{KaushikYousefianACC21}, we addressed a subclass of the distributed model \eqref{prob:initial_problem} formulated as  \eqref{prob:subclass_constrained_opt}. Extending \cite{KaushikYousefianACC21}, here we show that the convergence and rate analysis in \cite{KaushikYousefianACC21} can be extended to the much broader class of problems of the form \eqref{prob:initial_problem} without any degradation {of} the speed of the algorithm. In addressing hierarchical optimization problems, there have been other iterative methods proposed in works such as \cite{BeckSabach2014,SabachShtern2017,GarrigosRosascoVilla2018,XuViscosity2004}. We, however, note that all of these schemes can only address a subclass of \eqref{prob:initial_problem} in centralized regimes and under more restrictive assumptions such as strong convexity of the global objective function. We also note that compared to the existing methods that address VI problems (e.g., see \cite{Nem11, IusemJofreThompson2019,FarzadMathProg17,wang2015}), our work provides an avenue for addressing a {significantly} more general class of problems where VI is employed as a tool to characterize the constraints in distributed optimization. 
 
  \noindent  \textit{(ii) Distributed averaging scheme}: In pair-IG, we employ a distributed averaging scheme where agents can choose their initial averaged iterate arbitrarily and independent from each other. This relaxation in the proposed IG method appears to be novel, even for the classical IG schemes {in addressing \eqref{prob:canonical_DO_problem}}. 
  
  \noindent  \textit{(iii) Rate analysis in the solution space}: Motivated by the recent developments of iterative methods for MPECs \cite{CSY2021MPEC}, it is important to characterize the speed of the proposed scheme in the solution space. To this end, under strong convexity of the global objective function, in Theorem \ref{thm:part2} we derive agent-specific rate statements that compare the generated sequence of each agent with the so-called Tikhonov trajectory, that is defined as the trajectory of the unique solutions to a family of regularized optimization problems.
  
 \noindent  \textit{(iv) Preliminary numerical results}: To validate the theoretical results, we provide preliminary numerical experiments for computing the best equilibrium in a multi-agent traffic network problem. We also compare the performance of pair-IG with that of the existing IG methods in addressing constrained finite-sum problems. 
 }  
 
{\textbf{Outline:} The remainder of the  paper is organized as {follows}. In section \ref{sec:algo_outline}, we provide the algorithm outline, main assumptions, and preliminaries. Section \ref{sec:convergence_analysis} includes the convergence and rate analysis of the proposed scheme.  In section \ref{sec:convergence_analysis_strong_convexity} we  study the solution quality and derive bounds on {the} distance between the {iterates} from Algorithm \ref{alg:IR-IG_avg} and the Tikhonov trajectory. Section \ref{sec:numerical_implementation} contains numerical implementations, and concluding remarks follow in section \ref{sec:concludingrem}.}

{\textbf{Notation:} A function $f : \mathbb{R}^n \rightarrow \mathbb{R}$ is said to be in the class $C_{\mu,L}^{k,r}$ if $f$ is $\mu$-strongly convex in  $\mathbb{R}^n$, $k$ times continuously differentiable, and its $r^{\text{th}}$ derivative is Lipschitz continuous with constant $L$.  Throughout, a vector $x \in \mathbb{R}^n$ is assumed to be a column vector and $x^T$ denotes its transpose. For any vector $x\in \mathbb{R}^n$, we use $\|x\|$ to denote the $\ell_2$-norm. For a convex function $f:\mathbb{R}^n\rightarrow\mathbb{R}$ with the domain $\text{dom}(f)$ and  any $x\in\text{dom}(f)$, a vector $ \tilde{\nabla}f(x) \in\mathbb{R}^n$ is called a subgradient of $f$ at $x$ if $f(x) +\tilde{\nabla}f(x)^Tf(y-x)\leq f(y)$ for  all $y\in \text{dom}(f)$. We let $\partial f(x)$ denote the subdifferential set of  function $f$ at $x$.  The Euclidean  projection of vector $x$ onto set $X$ is denoted by $\mathcal{P}_X(x)$. We let $\mathbb{R}^n_{+}$ and $\mathbb{R}^n_{++}$ denote the set $\left\{x\in \mathbb{R}^n\ |\ x \geq 0\right \}$ and $\left\{x\in \mathbb{R}^n\ |\ x >0\right \}$, respectively. Given a set $S \subseteq \mathbb{R}^n$, we let $\mbox{int}(S)$ denote the interior of $S$. {Given an integer $m$, we let $[m]$ abbreviate the set $\{1,\dots, m\}$. }}

\section{Algorithm outline}\label{sec:algo_outline}

{In this section we present the main assumptions on problem \eqref{prob:initial_problem}, the outline of the proposed algorithm, and a few preliminary results that will be applied later in the rate analysis. Throughout the paper, we let $f(x)\triangleq \sum_{i=1}^m f_i(x)$ and $F(x)\triangleq \sum_{i=1}^m F_i(x)$ denote the global objective and global mapping in problem \eqref{prob:initial_problem}, respectively. 
\begin{assumption}[Properties of problem \eqref{prob:initial_problem}]\label{assum:initial_prob} Suppose the following conditions hold. \\
	\noindent(a) {Function $f_i:\mathbb{R}^n\rightarrow\mathbb{R}$ is real-valued and merely convex (possibly nondifferentiable) on its domain for all {$i \in [m]$}.}

	\noindent (b) Mapping $F_i:\mathbb{R}^n\rightarrow\mathbb{R}^n$ is real-valued, continuous, and merely monotone on its domain for all  {$i \in [m]$}. 
	
	\noindent (c) The set $X \subseteq \mbox{int}\left(\mbox{dom}(f)\cap \mbox{dom}(F)\right) $ is nonempty, convex, and compact.
\end{assumption}
\begin{remark}\label{rem:nonemptiness_SOL}\em
Under Assumption \ref{assum:initial_prob} we have the following immediate results. From Theorem 2.3.5 and Corollary 2.2.5 in \cite{FacchineiPang2003}, the set $\mbox{SOL}(X,F)$ is nonempty, convex, and compact. For all $i$, the nonemptiness of the subdifferential set $\partial f_i(x)$ for any $x \in \mbox{int}(\mbox{dom}(f_i))$ is implied from Theorem 3.14 in \cite{Beck2017}. Also, Theorem 3.16 in \cite{Beck2017} implies that $f_i$ has bounded subgradients over the compact set $X$. Further, mapping $F_i$ is bounded over the set $X$. 
\end{remark}
In view of compactness of the set $X$ and continuity of $f$, throughout the paper we let positive scalars $M_X< \infty$ and $M_f < \infty$ be defined as $M_X \triangleq \sup_{x \in X} \|x\|$ and $M_f \triangleq \sup_{x \in X} |f(x)|$, respectively.  We also let $f^*\in \mathbb{R}$ denote the optimal objective value of problem \eqref{prob:initial_problem}. In view of Remark \ref{rem:nonemptiness_SOL}, throughout we let scalars $C_F>0$ and $C_f>0$ be defined such that for all {$i\in [m]$} and for all $x \in X$ we have $ \|F_i(x)\| \leq \frac{C_F}{m}$, and $\left\|\tilde{\nabla}{ f_{i}\left( x\right)}\right\|\leq \frac{C_f}{m}$ for all $\tilde{\nabla}{ f_{i}\left( x\right)} \in \partial f_i(x)$. In the following, we comment on the Lipschitz continuity of the local and global objective functions.
\begin{remark}\label{rem:lipschitz_parameters}
\em Under Assumption \ref{assum:initial_prob} and from Theorem 3.61 in \cite{Beck2017}, function $f_i$ is Lipschitz continuous with the parameter $\frac{C_f}{m}$ over the set $X$, i.e., for all {$i \in [m]$} we have $|f_i(x)-f_i(y)|\leq \frac{C_f}{m}\|x-y\|$ for all $x,y \in X$. We also have $\|\tilde \nabla f(x) \| \leq C_f$ for all $x \in X$ and all $\tilde \nabla f(x) \in \partial f(x)$. This implies that $|f(x)-f(y)|\leq C_f\|x-y\|$ for all $x,y \in X$. 
\end{remark}}
We now present an overview of the proposed method given by Algorithm \ref{alg:IR-IG_avg}.  We use vector $x_{k,i}$ to denote the local copy of the global decision vector maintained by agent $i$ at iteration $k$.  At each iteration, agents update their iterates in a cyclic manner.  Each agent {$i\in[m]$} uses only its local information including the subgradient of the function $f_i$ and mapping $F_i$ and evaluates the regularized mapping $F_i + \eta_k \tilde \nabla f_i$ at $x_{k,i}$.  Here, $\gamma_{k}$ and $\eta_k$ denote the stepsize and the regularization parameter at iteration $k$, respectively.  Importantly, through employing an iterative regularization technique, we let both of these parameters be updated iteratively at suitable prescribed rates (cf.  Theorem \ref{thm:rates}).  Each agent computes and returns a weighted averaging iterative denoted by $\bar x_{k,i}$ where the weights are characterized in terms of the stepsize $\gamma_k$ and an arbitrary scalar $r \in [0,1)$. Notably,  this averaging technique is carried out in a distributed fashion in the sense that agents do not require to start from the same initialized averaging iterate. This is in contrast with the standard incremental gradient {schemes where} the averaging scheme is limited to a centralized initialization. 
\begin{algorithm}[t]
	\caption{projected averaging iteratively regularized Incremental subGradient (pair-IG)}\label{alg:IR-IG_avg}
\begin{algorithmic}[0]
	 \State \textbf{input}: Agent $1$ arbitrarily chooses an initial vector $ x_{0,1} \in X$. Agent $i$ arbitrarily chooses $\bar x_{0,i} \in X$, for all $i \in [m]$.  Let $S_{0} := \gamma_{0}^r$ with an arbitrary $0\leq r <1$.
		\For {$k = 0,1, \dots, {N-1}$}
		\State Update $S_{k+1}:=S_{k}+\gamma_{k+1}^r $
		\For {$i = 1, \dots, {m}$}
		\vspace{-0.2in}
		\State \begin{align
 		&x_{k,{i+1}} :=  \mathcal{P}_X \left({x_{k,i}}-\gamma_k\left(F_i\left(x_{k,i}\right)+ {\eta_{k}} \tilde{\nabla} f_{i}\left( x_{k,i}\right)\right)\right)\label{eqn:alg_step_1}\\
		&{\bar{x}_{k+1,i}:=\left(\tfrac{S_k}{S_{k+1}}\right)\bar{x}_{k,i}  + \left(\tfrac{\gamma_{k+1}^r}{S_{k+1}}\right) x_{k,i+1} } \label{eqn:alg_step_2}
		\end{align}
	 \EndFor
		\State Set $x_{k+1,1}:= x_{k,{m+1}}$
	 \EndFor
 	\State \textbf{return}: $\bar{x}_{N,i}$ for all $i\in [m]$ 
	\end{algorithmic}
\end{algorithm}
Next we show that {for any $i  \in [m]$, $\bar{x}_{N,i}$} is indeed a well-defined weighted average of $\bar x_{0,i}$ and {the iterates} $x_{k-1,i+1}$ for {$k\in [N]$}.
{\begin{lemma}\label{lem:averaging_seq}
	Consider the sequence {$\{\bar x_{k,i}\}$} generated by agent {$i\in [m]$} in Algorithm \ref{alg:IR-IG_avg}. For {$k\in \{0,\ldots,N\}$}, let us define the weights $\lambda_{k,N} \triangleq \tfrac{\gamma_k^r}{\sum_{j = 0}^{N} \gamma_j^r }$. Then for all {$i \in [m] $} we have 
	\begin{align*}
	\bar{x}_{N,i} = \lambda_{0,N}\bar  x_{0,i} + \sum_{k=1}^{N}\lambda_{k,N}x_{k-1,i+1}.\end{align*}
	 Further, for a convex set $X$ we have $\bar{x}_{N,i} \in X$. 
\end{lemma}}
\begin{proof}
	We use induction on $N\geq 0$ to show the equation.  For $N = 0$, from $\lambda_{0,0} = 1$ we have $\bar{x}_{0,i} = \lambda_{0,0}\bar x_{0,i}$. Now,  assume that the {equation} holds for some $N\geq 0$. {This implies}
	\begin{align}\label{eqn:avg_1}
		\bar{x}_{N,i} &= \lambda_{0,N} \bar x_{0,i} +  \sum_{k=1}^{N}\lambda_{k,N}x_{k-1,i+1}=  \frac{\gamma_0^r \bar   x_{0,i} +  \sum_{k=1}^{N}\gamma_k^r x_{k-1,i+1}}{\sum_{j = 0}^{N} \gamma_j^r }.
	\end{align}
	Using equation \eqref{eqn:avg_1}, we now show that the hypothesis statement holds for any $N+1$.  From equation \eqref{eqn:alg_step_2} we have
		$\bar{x}_{N+1,i}=\left(\frac{S_N}{S_{N+1}}\right)\bar{x}_{N,i}  + \left(\frac{\gamma_{N+1}^r}{S_{N+1}}\right) x_{N, i+1}$.  Note that from {equation \eqref{eqn:alg_step_2}} in Algorithm \ref{alg:IR-IG_avg} we have $S_k=\sum_{t=0}^{k}\gamma_{t}^r$ for all $k\geq 0$.  From this and using equation \eqref{eqn:avg_1} we obtain
	\begin{align*}
	\bar{x}_{N+1,i}&=\left(\frac{\sum_{t=0}^{N}\gamma_{t}^r}{\sum_{t=0}^{N+1}\gamma_{t}^r}\right)\bar{x}_{N,i}  + \left(\frac{\gamma_{N+1}^r}{\sum_{t=0}^{N+1}\gamma_{t}^r}\right) x_{N,i+1} = \frac{\gamma_0^r \bar  x_{0,i} +  \sum_{k=1}^{N}\gamma_k^r x_{k-1,i+1} + \gamma_{N+1}^r x_{N,i+1}}{\sum_{t=0}^{N+1}\gamma_{t}^r}\\&= \frac{\gamma_0^r \bar  x_{0,i} +  \sum_{k=1}^{N+1}\gamma_k^r x_{k-1,i+1}}{\sum_{t=0}^{N+1}\gamma_{t}^r}.	
	\end{align*}	
From the definition of $\lambda_{k,N}$ we conclude that the hypothesis holds for $N+1$ and thus, the result holds for all $N\geq 0$.  To show the second part, note that from the {initialization in Algorithm \ref{alg:IR-IG_avg}} and the projection in equation \eqref{eqn:alg_step_1}, we have $\bar x_{0,i}, x_{k-1,i+1}\in X$ for all $i$ and $k\geq 1$. From the first part, $\bar{x}_{N,i}$ is a convex combination of $\bar x_{0,i}, x_{0,i+1}, \ldots,x_{N-1,i+1}$. Therefore, from the convexity of the set $X$ we conclude that $\bar{x}_{N,i} \in X$.
\end{proof}
For the ease of presentation throughout the analysis, we define a sequence $\{x_k\}$ as follows.
\begin{definition}\label{def:convenient_notation}
Consider Algorithm \ref{alg:IR-IG_avg}. Let the sequence $\{x_k\}$ be {defined as  
$x_k \triangleq x_{k-1, m+1}=x_{k,1},$ for all $k \geq 1,$ with $x_0\triangleq x_{0,1}.$}
\end{definition}
In the following result, we characterize the distance between the local variable of {any} arbitrary agent with that of the first and the last agent at any given iteration. This result will be utilized in the analysis. 
\begin{lemma}\label{lem:neighbor_agents}
Consider Algorithm \ref{alg:IR-IG_avg}. Let Assumption \ref{assum:initial_prob} hold. Then the following inequalities hold for all {$i\in[m]$ and $k\geq 0$}

\noindent (a) $\| x_k - x_{k,i} \|\leq \frac{(i-1)\gamma_{k}\left(C_F+\eta_{k}C_f\right)}{m}$.

\noindent (b) $	\left\| x_{k,i+1} -  x_{k+1} \right\| \leq  \frac{(m-i){\gamma_{k}}\left({C_F} + {\eta_{k}} C_f\right)}{m}$. 
\end{lemma}
\begin{proof}
\noindent (a) Let $k\geq 0$ be an arbitrary integer. We use induction on $i$ to show this result. From Definition \ref{def:convenient_notation},  for $i=1$ and $k\geq 0$ we have $\|x_k - x_{k,1}\| = 0$, implying that the result holds for $i=1$. Now suppose the hypothesis statement holds for some {$i\in [m]$}. We have
	\begin{align*}
		\| x_k - x_{k,i+1} \| &=  \left\| \mathcal{P}_X \left(x_{k}\right) - \mathcal{P}_X \left( x_{k,i} - \gamma_{k}\left(F_i(x_{k,i}) + \eta_{k}\tilde{\nabla} f_i\left(x_{k,i}\right)\right)  \right)  \right\|\nonumber\\   &\leq  \left\| x_{k} - x_{k,i} \right\| +  \gamma_{k}\left\| F_i(x_{k,i}) + \eta_{k}\tilde{\nabla} f_i\left(x_{k,i}\right)  \right\|\nonumber\\ 
		& \leq  \left\| x_{k} - x_{k,i} \right\| +  \frac{\gamma_{k}\left(C_F + \eta_k C_f\right)}{m} \leq  \frac{i\gamma_{k}\left(C_F + \eta_k C_f\right)}{m},	
	\end{align*}
where the first inequality is obtained from the nonexpansivity property of the projection. Therefore the hypothesis statement holds for any $i$ and the proof of part (a) is completed. 

\noindent (b) To show this result, we use downward induction on {$i\in[m]$}. Note that the relation trivially holds for the base case $i=m$. Suppose it holds for some {$i \in \{2,\ldots,m\}$}. We show that it holds for $i-1$ as well.  From Definition \ref{def:convenient_notation} we have
\begin{align*}
	\left\| x_{k,i} -  x_{k+1} \right\| &= \| x_{k,i} - x_{k,i+1} +x_{k,i+1} -  x_{k,m+1} \| \leq \| x_{k,i} - x_{k,i+1} \| + \|x_{k,i+1} -  x_{k,m+1}\|.
\end{align*}
{From equation \eqref{eqn:alg_step_1}, the hypothesis statement, and the nonexpansivity property of the projection, we obtain}
\begin{align*}
	\left\| x_{k,i} -  x_{k+1} \right\|
	&\leq \left\| \mathcal{P}_X( x_{k,i}) - \mathcal{P}_X \left({x_{k,i}}-{\gamma_{k}}\left(F_i\left(x_{k,i}\right)+ {\eta_{k}} \tilde{\nabla} f_{i}\left( x_{k,i}\right)\right)\right)\right\| \\
	&+ \frac{(m-i){\gamma_{k}}\left({C_F} + {\eta_{k}} C_f\right)}{m}
	\leq \frac{(m-i+1){\gamma_{k}}\left({C_F} + {\eta_{k}} C_f\right)}{m}.
\end{align*}
This completes the proof of part (b).
\end{proof}
We note that the generated agent-wise iterates $\bar x_{k,i}$ in Algorithm \ref{alg:IR-IG_avg}, as the scheme proceeds,  may not be solutions to $\mbox{VI}(X,F)$ and so, they may not necessarily be feasible to problem \eqref{prob:initial_problem}.  To quantify the infeasibility of these iterates, we employ a dual gap function (cf. Chapter 1 in \cite{FacchineiPang2003}) defined as follows. 

\begin{definition}[The dual gap function {\cite{MarcotteZhu1998}}]\label{def:gap}\normalfont
Consider a closed convex set $X\subseteq \mathbb{R}^n$ and  a continuous mapping $F:X\rightarrow\mathbb{R}^n$. The dual gap function at $x\in X$ is defined as $\mbox{GAP}(x) \triangleq \sup_{y\in X} F(y)^T(x-y)$.
\end{definition}
Note that under Assumption \ref{assum:initial_prob}, the dual gap function is well-defined. This is because for any $x\in X$,  we have $\mbox{GAP}(x)\geq 0$. Further, it is known that when the mapping $F$ is continuous and monotone, $x \in \text{SOL}(X,F)$ if and only if $\mathrm{GAP}(x)=0$ (cf.~\cite{Nem11}). We conclude this section by presenting the following result that will be utilized in the rate analysis. 
\begin{lemma}\label{lem:harmonic_series_bound} Let $\beta \in [0,  1)$ and $\Gamma\geq 1$ be given scalars and $K$ be an integer. Then for all $K\geq \left(2^{\frac{1}{1-\beta}}-1\right)\Gamma$, we have \begin{align*}
\frac{{(K+\Gamma)}^{1-\beta}}{2(1-\beta)} \leq {\sum_{k = 0}^{K}}(k+\Gamma)^{-\beta}\leq \frac{{(K+\Gamma)}^{1-\beta}}{1-\beta}.
\end{align*}
\end{lemma}

\begin{proof}
See Appendix \ref{app:lem:harmonic_series_bound}.
\end{proof}

\section{Rate and complexity analysis}\label{sec:convergence_analysis}
{In this section we present the convergence and rate analysis of the proposed method under Assumption \ref{assum:initial_prob}.  After obtaining a preliminary inequality in Lemma \ref{lem:pseudo_bound_sequence} in terms of the sequence generated by the last agent, in Lemma \ref{lem:consensus_error_bound} we derive inequalities that relate the global objective and the dual gap function at the iterate of other agents with those of the last agent. Utilizing these results, in Proposition \ref{prop:error_bound_aIR-IG} we obtain agent-specific bounds on the objective function value and the dual gap function. Consequently, in Theorem \ref{thm:rates} we derive convergence rate statements under suitably chosen sequences for the stepsize and the regularization parameter. }
{\begin{lemma}\label{lem:pseudo_bound_sequence}
	Consider Algorithm \ref{alg:IR-IG_avg}. Let Assumption \ref{rem:nonemptiness_SOL} hold. Let $\{\gamma_{k}\}$ and $\{\eta_{{k}}\}$ be nonincreasing and strictly positive sequences.  For any arbitrary $y\in X$, for all $k\geq 0$ we have
	\begin{align}\label{eq_36}
			{2\gamma_{k}^r\left(\eta_{k}\left(f(x_{k}) -f(y) \right)   +F(y)^T\left(x_{k}-y\right)\right)}&\leq \gamma_{k}^{r-1} \left\| {x_{k}}-y\right\|^2  -\gamma_{k}^{r-1}\left\| x_{k+1} -  y \right\|^2    \nonumber\\
			& +\gamma_{k}^{r+1}\left(C_F +\eta_{{k}}C_f\right)^2.
	\end{align}
\end{lemma}}
\begin{proof}
	Let $y\in X$ be an arbitrary vector and $k\geq 0$ be fixed. From the update rule \eqref{eqn:alg_step_1}, for {$i \in  [m]$} we have
	\begin{align*}
	&\left\| x_{k,{i+1}} -  y \right\|^2=\left\| \mathcal{P}_X\left({x_{k,i}}-\gamma_k\left( F_i\left(x_{k,i}\right)    + {\eta_{k}} \tilde{\nabla} f_{i}\left( x_{k,i}\right)\right)\right)- \mathcal{P}_X (y) \right\|^2.
	\end{align*}
	{Employing the nonexpansivity of the projection we have}
	\begin{align*}
	&\left\| x_{k,{i+1}} -  y \right\|^2  \leq \left\| {x_{k,i}}-\gamma_k\left(F_i\left(x_{k,i}\right)+ {\eta_{k}} \tilde{\nabla} f_{i}\left( x_{k,i}\right)\right)-  y \right\|^2\\
&	=  \left\| {x_{k,i}}-y\right\|^2  + \gamma_k^2\left\|F_i\left(x_{k,i}\right)+ {\eta_{k}} \tilde{\nabla} f_{i}\left( x_{k,i}\right) \right\|^2   -  2\gamma_{k}\left(F_i\left(x_{k,i}\right)+ {\eta_{k}} \tilde{\nabla} f_{i}\left( x_{k,i}\right)\right)^T\left(x_{k,{i}} -  y \right).
	\end{align*}
	From the triangle inequality and recalling the  bounds on $\tilde\nabla f_i(x)$ and $F_i(x)$, we obtain
	\begin{align}\label{eqn:lemma2_1}
		&\left\| x_{k,{i+1}} -  y \right\|^2 \leq  \left\| {x_{k,i}}-y\right\|^2 + \gamma_{k}^2\left(\frac{C_F + \eta_{{k}}C_f}{m}\right)^2+ 2\gamma_{k}\left(F_i\left(x_{k,i}\right)+ {\eta_{k}} \tilde{\nabla} f_{i}\left( x_{k,i}\right)\right)^T\left(y-x_{k,{i}}  \right).
	\end{align}
	{The last term in the preceding relation is  bounded as follows}
	\begin{align*}
2\gamma_{k}\left(F_i\left(x_{k,i}\right)+ {\eta_{k}} \tilde{\nabla} f_{i}\left( x_{k,i}\right)\right)^T\left(y-x_{k,{i}}  \right)  
&= 2\gamma_{k}F_i\left(x_{k,i}\right)^T\left(y-x_{k,{i}}  \right)+2\gamma_{k} {\eta_{k}} \tilde{\nabla} f_{i}\left( x_{k,i}\right)^T\left(y-x_{k,{i}}  \right)\\
& \leq 2\gamma_{k}F_i\left(y\right)^T\left(y-x_{k,{i}}  \right)+2\gamma_{k} {\eta_{k}}\left( f_{i}\left( y\right)- f_{i}\left(x_{k,i}\right)\right),
	\end{align*}
where the last inequality is implied from the monotonicity of  $F_i$ and convexity of $f_i$. {Combining with equation \eqref{eqn:lemma2_1} we have}
	\begin{align*}
	\left\| x_{k,{i+1}} -  y \right\|^2   &\leq \left\| {x_{k,i}}-y\right\|^2 + \gamma_{k}^2\left(\frac{C_F + \eta_{{k}}C_f}{m}\right)^2 + 2\gamma_{k}F_i\left(y\right)^T\left(y-x_{k,{i}}  \right)\\
	&+2\gamma_{k} {\eta_{k}}\left( f_{i}\left( y\right)- f_{i}\left(x_{k,i}\right)\right).
	\end{align*}
Adding and subtracting $2\gamma_{k}F_i\left(y\right)^Tx_{k}+2\gamma_{k} {\eta_{k}} f_{i}\left( x_{k}\right)$ we {get}
	\begin{align*}
	\left\| x_{k,{i+1}} -  y \right\|^2   
	&\leq  \left\| {x_{k,i}}-y\right\|^2 + \gamma_{k}^2\left(\frac{C_F + \eta_{{k}}C_f}{m}\right)^2+ 2\gamma_{k}F_i\left(y\right)^T\left(y-x_{k}  \right)+2\gamma_{k} {\eta_{k}}\left( f_{i}\left( y\right)- f_{i}\left( x_{k}\right)  \right)\\
	& + 2\gamma_{k}\left( \left|F_i\left(y\right)^T\left(x_{k}-x_{k,i}  \right)\right| + {\eta_{k}} \left| f_{i}\left( x_k\right)- f_{i}\left( x_{k,i}\right)\right|\right).
	\end{align*}	
	Using the Cauchy-Schwarz inequality and Remark \ref{rem:lipschitz_parameters} we obtain
	\begin{align*}
	\left\| x_{k,{i+1}} -  y \right\|^2 &\leq  \left\| {x_{k,i}}-y\right\|^2 + \gamma_{k}^2\left(\frac{C_F + \eta_{{k}}C_f}{m}\right)^2+ 2\gamma_{k}F_i\left(y\right)^T\left(y-x_{k}  \right)+2\gamma_{k} {\eta_{k}}\left( f_{i}\left( y\right)- f_{i}\left( x_{k}\right)  \right) \\
	&+ 2\gamma_{k}\left(\frac{C_F}{m}\left\|x_{k}-x_{k,i}  \right\|+\frac{ {\eta_{k}}C_f}{m}\left\|x_k-  x_{k,i}\right\|\right).
	\end{align*}
	{Summing over $i\in[m]$} and considering Definition  \ref{def:convenient_notation} we have
	\begin{align*}
	\left\| x_{k+1} -  y \right\|^2  &\leq  \left\| {x_{k}}-y\right\|^2 {+ \frac{\gamma_{k}^2\left(C_F + \eta_{{k}}C_f\right)^2}{m}}+ 2\gamma_{k}F\left(y\right)^T\left(y-x_{k}  \right)+2\gamma_{k} {\eta_{k}}\left( f\left( y\right)- f\left( x_{k}\right)  \right)\\
	&+ \frac{2\gamma_{k}\left(C_F+{\eta_{k}}C_f\right)}{m}\sum_{i = 1}^m \left\|x_{k}-x_{k,i}  \right\|{.}
	\end{align*}
From Lemma \ref{lem:neighbor_agents} we obtain
	\begin{align*}
		\left\| x_{k+1} -y \right\|^2& \leq\left\| {x_{k}}-y\right\|^2 + {\frac{\gamma_{k}^2\left(C_F + \eta_{{k}}C_f\right)^2}{m}}  +2\gamma_kF(y)^T\left(y-x_{k}\right) + 2\gamma_{k}\eta_{k}\left(f(y) - f(x_{k})\right)  \\
		&+ \frac{2 \gamma_{k}\left(C_F +\eta_{k}C_f\right)}{m}\sum_{i=1}^m \frac{(i-1)\gamma_{k}\left(C_F + \eta_k C_f\right)}{m}   \\
	    &=  \left\| {x_{k}}-y\right\|^2 + \gamma_{k}^2\left(C_F+\eta_{{k}}C_f\right)^2 +2\gamma_kF(y)^T\left(y-x_{k}\right)+  2\gamma_{k}\eta_{k}\left(f(y) - f(x_{k})\right).
		\end{align*}
Multiplying the both sides by $\gamma_{k}^{r-1}$ {we can obtain} the result.
\end{proof}
In the next result we provide inequalities that relate the objective function and the dual gap function at the generated averaged iterate of the last agent with that of any other agent, respectively. This result will be utilized in Proposition \ref{prop:error_bound_aIR-IG}.
\begin{lemma}\label{lem:consensus_error_bound}
Consider problem \eqref{prob:initial_problem} and the sequences $\{\bar{x}_{N,i}\}$ generated in Algorithm \ref{alg:IR-IG_avg} for {$i \in [m]$} for some $N\geq 1$. Let Assumption \ref{assum:initial_prob} hold {and let $\{\gamma_k\}$ and $\{\eta_k\}$} be strictly positive {and} nonincreasing sequences. Then for any {$i \in [m]$} we have 
\begin{subequations}
	\begin{align}
		&f(\bar x_{N,i}) - f(\bar x_{N,m}) \leq C_f\lambda_{0,N}\left\|\bar   {x}_{0,i} -\bar   {x}_{0,m}\right\|+ \frac{(m-i)C_f\left({C_F} + \eta_0 C_f\right)}{m}\sum_{{k=0}}^{N} \lambda_{k,N}\gamma_{k},\label{eqn:consensus_error_bound_f} \\
		&\mathrm{GAP}\left(\bar x_{N,i}\right) - \mathrm{GAP}(\bar x_{N,m})  \leq {C_F}\lambda_{0,N}\left\| \bar  {x}_{0,i} -\bar  {x}_{0,m}\right\|+ \frac{(m-i){C_F}\left({C_F} + \eta_0 C_f\right)}{m}\sum_{{k=0}}^{N} \lambda_{k,N}\gamma_{k},\label{eqn:consensus_error_bound_GAP}
	\end{align}
\end{subequations}
where $\lambda_{k,N} \triangleq \tfrac{\gamma_k^r}{\sum_{j = 0}^{N} \gamma_j^r }$ for {$k \in \{0,\ldots,N\}$}.
\end{lemma}
\begin{proof} {Note that the results are trivial when $m=1$. Throughout, we assume that $m\geq 2$. From the Lipschitz continuity of function $f$ from Remark \ref{rem:lipschitz_parameters} and invoking Lemma \ref{lem:averaging_seq},  we can write the following for all {$i \in [m]$}.
{\begin{align}\label{eqn:consensus_error_begin}
		&f(\bar x_{N,i}) - f(\bar x_{N,m}) \leq C_f\lambda_{0,N}\left\| \bar {x}_{0,i} -\bar  {x}_{0,m}\right\| +C_f \sum_{k=1}^{N} \lambda_{k,N} \left\| x_{k-1,i+1} -  x_{k-1,m+1} \right\|.
\end{align}}
Next, using Lemma \ref{lem:neighbor_agents}(b) for any $k\geq 1$ and {$i\in[m]$} we have
\begin{align}\label{eqn:induction_bound_consensus}
	\left\| x_{k-1,i+1} -  x_{k-1,m+1} \right\| \leq  \frac{(m-i){\gamma_{k-1}}\left({C_F} + {\eta_{k-1}} C_f\right)}{m}. 
	\end{align}
From \eqref{eqn:consensus_error_begin},  \eqref{eqn:induction_bound_consensus}, and the nonincreasing sequence $\{\eta_k\}$, we have
\begin{align*}
	f(\bar x_{N,i}) - f(\bar x_{N,m})  
	& \leq  \frac{(m-i)C_f\left({C_F} + {\eta_0} C_f\right)}{m}\sum_{k=1}^{N} \lambda_{k,N} {\gamma_{k-1}}  + C_f\lambda_{0,N}\left\|  \bar {x}_{0,i} -  \bar {x}_{0,m}\right\|.
\end{align*}
{Since $\{\gamma_k\}$ is nonincreasing  and $0\leq r <1$, we obtain}
\begin{align*}
\sum_{k=1}^{N} \lambda_{k,N} {\gamma_{k-1}} & \leq \frac{1}{\sum_{j=0}^N \gamma_j^r}\sum_{k=1}^{N}  {\gamma_{k-1}^{r+1}}  \leq \frac{1}{\sum_{j=0}^N \gamma_j^r}\sum_{k=0}^{N-1}  {\gamma_{k}^{r+1}}  \leq \frac{1}{\sum_{j=0}^N \gamma_j^r}\sum_{k=0}^{N}  {\gamma_{k}^{r+1}} =\sum_{k=0}^{N} \lambda_{k,N} {\gamma_{k}}.
\end{align*}
From the last two relations we obtain equation \eqref{eqn:consensus_error_bound_f}. Next we show \eqref{eqn:consensus_error_bound_GAP}. From Definition \ref{def:gap} we have
\begin{align*}
	&\mathrm{GAP}(\bar x_{N,i}) = \sup_{y \in X}F(y)^T\left(\bar x_{N,i} - y\right)  = \sup_{y \in X} F(y)^T\left(\bar x_{N,i} + \bar x_{N,m} -\bar x_{N,m} - y\right)	\\
	& \leq\sup_{y \in X} F(y)^T\left(\bar x_{N,i} - \bar x_{N,m}\right) + \sup_{y \in X}  F(y)^T\left( \bar x_{N,m} - y\right) \leq {C_F} \left\|\bar x_{N,i} - \bar x_{N,m}\right\| + \mathrm{GAP}\left( \bar x_{N,m} \right),
\end{align*}
Rearranging the terms we obtain $\mathrm{GAP}(\bar x_{N,i}) -  \mathrm{GAP}\left( \bar x_{N,m} \right) \leq {C_F} \left\|\bar x_{N,i} - \bar x_{N,m}\right\|$. The rest of the proof can be done in a similar fashion to the proof of \eqref{eqn:consensus_error_bound_f}.}
\end{proof}
{Next we construct agent-wise error bounds in terms of the {objective function value} and the  dual gap function at the averaged iterates generated in Algorithm \ref{alg:IR-IG_avg}.
\begin{proposition}[Agent-wise error bounds]\label{prop:error_bound_aIR-IG}
	Consider problem \eqref{prob:initial_problem} and  the averaged sequence $\{\bar{x}_{k,i}\}$ generated by agent $i$ in Algorithm \ref{alg:IR-IG_avg} for {$i \in [m]$}. Let Assumption \ref{assum:initial_prob} hold and $\{\gamma_{k}\}$ and $\{\eta_{k}\}$ be nonincreasing and strictly positive sequences. {Then we have    for  $i \in [m]$,  $N\geq 1$, and $r \in [0,1)$:}
\begingroup
\allowdisplaybreaks
\begin{align*}
(a)\ &f(\bar{x}_{N,i}) -f^*\leq\left({\sum_{k = 0}^{N} \gamma_{k}^r}\right)^{-1}\left(\frac{2M_X^2\gamma_{N}^{r-1}}{\eta_{N}} +  {\frac{\left({C_F}+\eta_{0}C_f\right)^2}{2}} \sum_{k = 0}^{N}\frac{\gamma_{k}^{r+1}}{\eta_k}  \right. \\ &  \left.  + \frac{(m-i)C_f\left(C_F + \eta_0 C_f\right)}{m}\sum_{k=0}^{N}\gamma_{k}^{r+1}+{\gamma_{0}^rf(\bar  x_{0,m}) - \gamma_{0}^rf(x_{0,1})}+{C_f\gamma_{{0}}^r\left\| \bar {x}_{0,i} -\bar  {x}_{0,m}\right\|} \right).\\
(b)\ &\mathrm{GAP}(\bar{x}_{N,i})  \leq\left(\sum_{k = 0}^{N} \gamma_{k}^r\right)^{-1}\left(2M_X^2  \gamma_{N}^{r-1} +2M_f\sum_{k = 0}^{N}\gamma_{k}^r\eta_{k}+ { \frac{\left(C_F+C_f\eta_{{0}}\right)^2}{2}\sum_{k=0}^{N}\gamma_{k}^{r+1}}\right.\\
&\left.  + \frac{(m-i)C_F\left(C_F + \eta_0 C_f\right)}{m}\sum_{k=0}^{N}\gamma_{k}^{r+1}+\gamma_0^rC_F\|\bar x_{0,m}-x_{0,1} \|+{C_F\gamma_{{0}}^r\left\| \bar {x}_{0,i} -\bar  {x}_{0,m}\right\|} \right).
\end{align*}
\endgroup
\end{proposition}}
\begin{proof}
		{(a) Let $x^* \in X$ denote an arbitrary optimal solution to problem \eqref{prob:initial_problem}. From feasibility of $x^*$ we have  $F(x^*)\left(x_k-x^*\right)\geq 0$. Substituting $y$ by $x^*$ in relation \eqref{eq_36} and using the preceding relation we have
		\begin{align*}
			2\gamma_{k}^r\eta_{k}\left(f(x_{k}) -f^* \right)  &\leq  \gamma_{k}^{r-1} \left\| {x_{k}}-x^*\right\|^2 -\gamma_{k}^{r-1}\left\| x_{k+1} -  x^* \right\|^2 +  {\gamma_{k}^{r+1}\left(C_F +\eta_{{k}}C_f\right)^2} .
		\end{align*}
	Dividing both sides by $2\eta_k$ we have
		\begin{align}\label{eq_38}
			 \gamma_{k}^r\left(f(x_{k}) -f^* \right) \leq    \frac{\gamma_{k}^{r-1}}{2\eta_{{k}}}\left( \left\| {x_{k}}-x^*\right\|^2  - \left\| x_{k+1} -  x^* \right\|^2\right)   +{ \frac{\gamma_{k}^{r+1}}{2\eta_k}\left({C_F}+\eta_{k}C_f\right)^2}.
		\end{align}
		Adding and subtracting {the term} $\frac{\gamma_{k-1}^{r-1}}{2\eta_{k-1}}\|x_k-x^*\|^2$ we have
		\begin{align}\label{eqn:prop_1_1}
			\gamma_{k}^r\left(f(x_{k}) -f^* \right) &\leq  \frac{\gamma_{k-1}^{r-1}}{2\eta_{k-1}}\|x_k-x^*\|^2\nonumber -\frac{\gamma_{k}^{r-1}}{2\eta_{{k}}}\left\| x_{k+1} -  x^* \right\|^2 \\ 
			& + \underbrace{\left(\frac{\gamma_{k}^{r-1}}{2\eta_{k}} - \frac{\gamma_{k-1}^{r-1}}{2\eta_{k-1}}\right)}_{\text{term 1}}\|x_k-x^*\|^2 + { \frac{\gamma_{k}^{r+1}}{2\eta_k}\left({C_F}+\eta_{k}C_f\right)^2}.
		\end{align}
		Recalling {the definition of scalar $M_X$} we have
		\begin{align}\label{eq:bound_X}
			\|x_k-x^*\|^2 \leq 2\|x_k\|^2 + 2\|x^*\|^2  \leq 4M_X^2. 
		\end{align}
		Taking into account $r < 1$,  the nonincreasing property of
		the sequences $\{\gamma_k\}$ and $\{\eta_k\}$, we have: term 1 $\geq 0 $. Using \eqref{eq:bound_X} and taking summation from \eqref{eqn:prop_1_1} over {$k \in [N]$},  we obtain
		\begin{align}\label{eq_39}
			\sum_{k = 1}^{N} \gamma_{k}^r\left(f(x_{k}) -f^* \right)&\leq  \frac{\gamma_{0}^{r-1}}{2\eta_{0}}\|x_1-x^*\|^2-\frac{\gamma_{N}^{r-1}}{2\eta_{N}}\|x_{N+1}-x^*\|^2\nonumber\\ 
			&+\left(\frac{\gamma_{N}^{r-1}}{2\eta_{N}}-\frac{\gamma_{0}^{r-1}}{2\eta_{0}}\right)4M_X^2+{\frac{\left({C_F}+\eta_{0}C_f\right)^2}{2} \sum_{k = 1}^{N}\frac{\gamma_{k}^{r+1}}{\eta_k}}.
		\end{align}
		Rewriting equation \eqref{eq_38} for $k = 0$ and then, adding and subtracting $f(\bar x_{0,m})$, we have
		\begin{align*}
		\gamma_{0}^r\left(f(\bar x_{0,m})  -f^* + f(x_{0}) - f(\bar x_{0,m}) \right)&\leq    \frac{\gamma_{0}^{r-1} \left\| {x_{0}}-x^*\right\|^2}{2\eta_{{0}}} -\frac{\gamma_{0}^{r-1}\left\| x_{1} -  x^* \right\|^2}{2\eta_{{0}}}   \\
		&  + {\left({C_F}+\eta_{0}C_f\right)^2 \frac{\gamma_{0}^{r+1}}{2\eta_0}}.
		\end{align*}
		Adding the preceding equation with \eqref{eq_39} we obtain
		\begin{align*}
		&\gamma_{0}^r\left(f(\bar x_{0,m})  -f^* \right)	+\sum_{k = 1}^{N} \gamma_{k}^r\left(f(x_{k}) -f^* \right)\leq  \frac{\gamma_{0}^{r-1} \left\| {x_{0}}-x^*\right\|^2}{2\eta_{{0}}} + 2M_X^2\left(\frac{\gamma_{N}^{r-1}}{\eta_{N}}-\frac{\gamma_{0}^{r-1}}{\eta_{{0}}}\right)\\
		& -\frac{\gamma_{N}^{r-1}}{2\eta_{N}} \left\| x_{N+1} -  x^* \right\|^2   +{\frac{\left({C_F}+\eta_{0}C_f\right)^2}{2} \sum_{k = 0}^{N}\frac{\gamma_{k}^{r+1}}{\eta_k}}+{\gamma_{0}^r\left(f(\bar x_{0,m})  -f(x_0) \right)}.
		\end{align*}
{From  \eqref{eq:bound_X}} and neglecting the nonpositive term we obtain
		\begin{align*}
			  \gamma_{0}^r\left(f(\bar x_{0,m})  -f^* \right)	+\sum_{k = 1}^{N} \gamma_{k}^r\left(f(x_{k}) -f^* \right)&\leq  \frac{2M_X^2\gamma_{N}^{r-1}}{\eta_{N}}   + {\frac{\left({C_F}+\eta_{0}C_f\right)^2}{2} \sum_{k =0}^{N}\frac{\gamma_{k}^{r+1}}{\eta_k}}\\
			  &+{\gamma_{0}^r\left(f(\bar x_{0,m})  -f(x_0) \right)}.
		\end{align*}
		Next, dividing both sides by  $\sum_{k = 0}^{N} \gamma_{k}^r$ we have
		\begin{align*
			\frac{{\gamma_{0}^r f(\bar x_{0,m}) }	+\sum_{k = 1}^{N} \gamma_{k}^rf(x_{k}) }{\sum_{k = 0}^{N} \gamma_{k}^r} -f^*
			& \leq\left(\sum_{k = 0}^{N} \gamma_{k}^r\right)^{-1}\left(\frac{2M_X^2\gamma_{N}^{r-1}}{\eta_{N}}+ \frac{\left({C_F}+\eta_{0}C_f\right)^2}{2} \sum_{k = 0}^{N}\frac{\gamma_{k}^{r+1}}{\eta_k}\right.\\
			&\left.  +{\gamma_{0}^r\left(f(\bar x_{0,m})  -f(x_0) \right)}\right).
		\end{align*}
Taking into account the convexity of  $f$ we have
\begin{align*}
f\left(\frac{\gamma_{0}^r\bar x_{0,m} + \sum_{k = 1}^{N} \gamma_{k}^rx_{k-1,m+1}}{\sum_{k = 0}^{N} \gamma_{k}^r}\right)  \le\frac{\gamma_{0}^rf(\bar x_{0,m}) + \sum_{k = 1}^{N} \gamma_{k}^rf(x_{k-1,m+1})}{\sum_{k = 0}^{N} \gamma_{k}^r}.
\end{align*} 
Invoking Lemma \ref{lem:averaging_seq}, from the preceding two relations we obtain
	\begin{align*}
		&f(\bar{x}_{N,m}) -f^* \leq   \left( \sum_{k = 0}^{N} \gamma_{k}^r\right)^{-1}\left(\frac{2M_X^2\gamma_{N}^{r-1}}{\eta_{N}}+  {\frac{\left({C_F}+\eta_{0}C_f\right)^2}{2}} \sum_{k = 0}^{N}\frac{\gamma_{k}^{r+1}}{\eta_k}   +{\gamma_{0}^rf( \bar x_{0,m}) - \gamma_{0}^rf(x_{0,1})}\right).
	\end{align*}
Adding equation \eqref{eqn:consensus_error_bound_f} with the preceding inequality we obtain the desired result. 
		
\noindent(b) From equation \eqref{eq_36}, for an arbitrary $y\in X$ we have
		\begin{align*}
		2\gamma_k^rF(y)^T\left(x_{k}-y\right)  & \leq   \gamma_{k}^{r-1} \left(\left\| {x_{k}}-y\right\|^2-\left\| x_{k+1} - y \right\|^2\right)\\
		& + 2\gamma_{k}^r\eta_{k}\left(f\left(y\right) - f(x_{k})  \right)+ { \gamma_{k}^{r+1}\left(C_F +\eta_{{k}}C_f\right)^2}.
		\end{align*}
From the triangle inequality and definition of $M_f$ we have $|f(y) - f(x_k)|\leq 2M_f$. We obtain:
		\begin{align}\label{eq39}
				2\gamma_k^rF(y)^T\left(x_{k}-y\right)  & \leq  \gamma_{k}^{r-1} \left(\left\| {x_{k}}-y\right\|^2  -\left\| x_{k+1} - y \right\|^2\right)  \nonumber\\
				& + 4\gamma_{k}^r\eta_{k}M_f + \gamma_{k}^{r+1}\left(C_F +\eta_{{k}}C_f\right)^2.
		\end{align}
Adding and subtracting $\gamma_{k-1}^{r-1}\|x_k-y\|^2$, we have:
		\begin{align}\label{eqn:prop1_partb_ineq1}
			 2\gamma_k^rF(y)^T\left(x_{k}-y\right)   &\leq  \gamma_{k-1}^{r-1} \left\| {x_{k}}-y\right\|^2-\gamma_{k}^{r-1}\left\| x_{k+1} -  y \right\|^2 +4\gamma_{k}^r\eta_{k}M_f\nonumber \\
			 & + \underbrace{\left(\gamma_{k}^{r-1} -\gamma_{k-1}^{r-1} \right) \left\| {x_{k}}-y\right\|^2}_{\text{term 2}} +  {\gamma_{k}^{r+1}\left(C_F +\eta_{{k}}C_f\right)^2}  .
		\end{align}
		Using the nonincreasing property of $\{\gamma_{k}\}$ and recalling $0\leq r<1$, we have $ \gamma_{k}^{r-1} -\gamma_{k-1}^{r-1} \geq 0 $. Thus, we can write: term 2 $\leq \left( \gamma_{k}^{r-1} -\gamma_{k-1}^{r-1} \right) 4M_X^2 $. Taking summation over {$k\in [N]$ in equation \eqref{eqn:prop1_partb_ineq1}} and dropping a nonpositive term we obtain
		\begin{align}\label{eq40}
			2\sum_{k = 1}^{N}\gamma_k^rF(y)^T\left(x_{k}-y\right) & \leq  \gamma_{0}^{r-1} \left\| {x_{1}}-y\right\|^2 +4M_f\sum_{k = 1}^{N}\gamma_{k}^r\eta_{k} + 4M_X^2 \left(\gamma_{N}^{r-1} -\gamma_{0}^{r-1} \right)    \nonumber \\   
			 &  + { \left(C_F+\eta_{{0}}C_f\right)^2\sum_{k = 1}^{N}\gamma_{k}^{r+1}}.
\end{align}
Writing equation \eqref{eq39} for $k=0$ and adding and subtracting $2\gamma_0^rF(y)^T\bar x_{0,m}$, we have
		\begin{align*}
			 2\gamma_0^rF(y)^T\left(\bar x_{0,m}-y + x_0 -\bar  x_{0,m}\right)    &\leq  4\gamma_{0}^r\eta_{0}M_f + \gamma_{0}^{r-1} \left(\left\| {x_{0}}-y\right\|^2   -\left\| x_{1} -  y \right\|^2\right) \\
			 &  +  {\gamma_{0}^{r+1}\left(C_F+\eta_{{0}}C_f\right)^2}\nonumber.
		\end{align*}
Adding the preceding relation with equation \eqref{eq40} we have
		\begin{align*}
			&2\gamma_0^rF(y)^T\left(\bar x_{0,m}-y \right)  + 2\sum_{k = 1}^{N}\gamma_k^rF(y)^T\left(x_{k}-y\right) 
			 \leq  4M_X^2 \left( \gamma_{N}^{r-1}- \gamma_{0}^{r-1} \right)   + \left(C_F+\eta_{{0}}C_f\right)^2\sum_{k=0}^{N}\gamma_{k}^{r+1}\\
			&+ \gamma_{0}^{r-1} \left\| {x_{0}}-y\right\|^2+4M_f\sum_{k = 0}^{N}\gamma_{k}^r\eta_{k}+{\underbrace{2\gamma_0^rF(y)^T\left(\bar x_{0,m}-x_0 \right)}_{\text{term 3}}} .
		\end{align*}
		Using the Cauchy-Schwarz inequality we have: term 3 $\leq 2\gamma_0^rC_F\|x_{0,m}-x_{0,1} \|$. We also have $\|x_0-y\|^2 \leq 4M_X^2$. Dividing the both sides of the preceding inequality by  $2\sum_{k = 0}^{N} \gamma_{k}^r$ and invoking Lemma \ref{lem:averaging_seq}, we have
		 \begin{align*}
		 	F(y)^T\left({\bar{x}_{N,m}}-y\right)& \leq  \left(\sum_{k = 0}^{N}\gamma_k^r\right)^{-1}\hspace{-0.25cm}\left(2M_X^2  \gamma_{N}^{r-1} +2M_f\sum_{k = 0}^{N}\gamma_{k}^r\eta_{k}  + { \frac{\left(C_F+\eta_{{0}}C_f\right)^2}{2}\sum_{k=0}^{N}\gamma_{k}^{r+1}}\right.\\
		 	& \left. +{\gamma_0^rC_F\|\bar x_{0,m}-x_{0,1} \|}\right).
		 \end{align*}
	Taking the supremum on both sides with respect to $y$ over the set $X$ and recalling Definition \ref{def:gap} we have
	 \begin{align*}
	  \mathrm{GAP}\left({\bar{x}_{N,m}}\right) 	&\leq  \left(\sum_{k = 0}^{N}\gamma_k^r\right)^{-1}\left( {\frac{\left(C_F+\eta_{{0}}C_f\right)^2}{2}\sum_{k=0}^{N}\gamma_{k}^{r+1}}+2M_f\sum_{k = 0}^{N}\gamma_{k}^r\eta_{k}+2M_X^2  \gamma_{N}^{r-1}\right.\\
		&\left.   +{\gamma_0^rC_F\|\bar x_{0,m}-x_{0,1} \|}\right).
	 \end{align*}
Adding equation \eqref{eqn:consensus_error_bound_GAP} with the preceding inequality we obtain the desired inequality.}
\end{proof}
{In the following we present the main result of this section. We provide non-asymptotic rate statements for each agent {$i \in [m]$} in terms of suboptimality measured by the global objective function, and infeasibility characterized by the dual gap function. We note that unlike the analysis of the standard incremental gradient schemes in the literature, here we provide these rate results for individual agents {$i \in [m]$}. 
\begin{theorem}[Agent-wise rate statements for Algorithm \ref{alg:IR-IG_avg}]\label{thm:rates}
Consider problem \eqref{prob:initial_problem}.  Let the averaged sequence $\{\bar{x}_{k,i}\}$ be generated by agent {$i \in[m]$} using Algorithm \ref{alg:IR-IG_avg}. Let Assumption \ref{assum:initial_prob} hold. Let the stepsize sequence $\{\gamma_{k}\}$ and the regularization sequence $\{\eta_{{k}}\}$ be updated using $\gamma_{k} := \frac{\gamma_{0}}{\sqrt{k+1}}$ and $\eta_{k} := \frac{\eta_{0}}{(k+1)^b}$, respectively, where $\gamma_{0}, \eta_{0} > 0$ and $0< b <0.5$. Then the following inequalities hold for all {$i\in[m]$}, all $N\geq 2^{\frac{2}{1-r}}-1$, and all $r \in [0,1)${:}
\begingroup
\allowdisplaybreaks
\begin{align}
(a)\  & f(\bar{x}_{N,i}) -f^* \leq\frac{2-r}{(N+1)^{0.5-b}}\left(\frac{2M_X^2}{\eta_0\gamma_0}+  \frac{\gamma_0\left(C_F+\eta_0C_f\right)^2}{\eta_0(1-r+2b)} \right.\nonumber\\
 & \left. + f( \bar x_{0,m}) - f(x_{0,1})  +C_f\left\|\bar  {x}_{0,i} -\bar  {x}_{0,m}\right\| + {\frac{{2(m-i)}\gamma_{0}C_f\left(C_F + \eta_0 C_f\right)}{m(1-r)}} \right).\label{eqn:rate_f}\\ 
(b)\   & \mathrm{GAP}(\bar{x}_{N,i}) \leq 
	 \frac{2-r}{(N+1)^{b}} \left(\frac{2M_X^2}{\gamma_0} +\frac{2M_f \eta_0}{1-0.5r-b} \right.\nonumber\\
	 & \left.+{C_F\|\bar x_{0,m}-x_{0,1} \|}  + C_F\left\| \bar {x}_{0,i} - \bar {x}_{0,m}\right\|+ \frac{\left(C_F+\eta_0C_f\right)^2\gamma_0}{1-r}+\frac{2{(m-i)}C_F\left(C_F + \eta_0 C_f\right) \gamma_0}{m(1-r)}\right).		
		\label{eqn:rate_gap}
\end{align}
\endgroup
\end{theorem}}
\begin{proof} 
{(a) Consider the inequality in Proposition \ref{prop:error_bound_aIR-IG}(a). Substituting $\gamma_k$ and $\eta_k$ by their update rules we obtain
\begin{align*}
	&f(\bar{x}_{N,i}) -f^* \leq \left(\sum_{k=0}^{N}\frac{\gamma_0^r}{(k+1)^{0.5r}}\right)^{-1}\left(\frac{2M_X^2(N+1)^{0.5(1-r)+b}}{\eta_0\gamma_0^{1-r}} \right.\nonumber\\
	&  \left.  + \frac{\left(C_F+\eta_0C_f\right)^2}{2}\sum_{k=0}^{N}\frac{\gamma_0^{1+r}}{\eta_0(k+1)^{0.5(1+r)-b}}+ {\frac{{(m-i)}C_f\left(C_F + \eta_0 C_f\right)}{m}\sum_{k=0}^{N}\frac{\gamma_{0}^{r+1}}{(k+1)^{0.5(1+r)}}}\right.\nonumber\\
	&\left.+\gamma_{0}^rf(\bar  x_{0,m}) - \gamma_{0}^rf(x_{0,1})+{C_f\gamma_{{0}}^r\left\| \bar {x}_{0,i} -\bar  {x}_{0,m}\right\|} \right).
\end{align*}
In the next step, to apply Lemma \ref{lem:harmonic_series_bound} we need to ensure that the conditions in that result are met. From $0\leq r<1$ and $0<b<0.5$, we have $0\leq 0.5r <1$, $0\leq 0.5(1+r)-b <1$, $0\leq 0.5r+b<1$, and $0 \leq 0.5(1+r)<1$.  Further, from $N\geq 2^{\frac{2}{1-r}}-1$, $0<b<0.5$, and $0\leq r<1$ we have that $N\geq \max\left\{2^{1/(1-0.5r)},2^{1/(1-0.5(1+r)+b)}, 2^{1/(1-0.5(1+r))}\right\}-1$. 

{\noindent Therefore, all the necessary conditions of Lemma \ref{lem:harmonic_series_bound} are met.}
 \begin{align*
 &	{f(\bar{x}_{N,i}) -f^*} \leq \left(\frac{\gamma_0^r(N+1)^{1-0.5r}}{2(1-0.5r)}\right)^{-1} \left(\frac{2M_X^2(N+1)^{0.5(1-r)+b}}{\eta_0\gamma_0^{1-r}} \right.\nonumber\\
 & \left.   +  \frac{\gamma_0^{1+r}\left(C_F+\eta_0C_f\right)^2(N+1)^{1-0.5(1+r)+b}}{2\eta_0(1-0.5(1+r)+b)} + {\frac{{(m-i)}C_f\left(C_F + \eta_0 C_f\right)\gamma_{0}^{r+1}(N+1)^{1-0.5(1+r)}}{m(1-0.5(1+r))}} \right.\nonumber\\
 & \left.+{\gamma_{0}^rf( \bar x_{0,m}) - \gamma_{0}^rf(x_{0,1})}+{C_f\gamma_{{0}}^r\left\|\bar  {x}_{0,i} - \bar {x}_{0,m}\right\|} \right).
 \end{align*}
From the preceding relation we obtain
  \begin{align*} 
 &	{f(\bar{x}_{N,i}) -f^*} \leq \left({2-r}\right)\left(\frac{2M_X^2}{\eta_0\gamma_0(N+1)^{0.5-b}}+  \frac{\gamma_0\left(C_F+\eta_0C_f\right)^2}{2\eta_0(1-0.5(1+r)+b)(N+1)^{0.5-b}}    \right.\nonumber\\
 & \left. + {\frac{{(m-i)}C_f\left(C_F + \eta_0 C_f\right)\gamma_{0}}{m(1-0.5(1+r))(N+1)^{0.5}}} +\frac{f( \bar x_{0,m}) - f(x_{0,1})+C_f\left\|\bar  {x}_{0,i} - \bar {x}_{0,m}\right\|}{(N+1)^{1-0.5r}} \right).
 \end{align*}
Factoring out $1/(N+1)^{0.5-b}$ we obtain
  \begin{align*} 
 &	{f(\bar{x}_{N,i}) -f^*} \leq \frac{2-r}{(N+1)^{0.5-b}}\left(\frac{2M_X^2}{\eta_0\gamma_0}+  \frac{\gamma_0\left(C_F+\eta_0C_f\right)^2}{2\eta_0(1-0.5(1+r)+b)} \right.\nonumber\\
 & \left.   + {\frac{{(m-i)}C_f\left(C_F + \eta_0 C_f\right)\gamma_{0}}{m(1-0.5(1+r))(N+1)^{b}}}+\frac{f(\bar  x_{0,m}) - f(x_{0,1})+C_f\left\| \bar {x}_{0,i} - \bar {x}_{0,m}\right\|}{(N+1)^{0.5-0.5r+b}} \right).
 \end{align*}
Note that from $b>0$ and $r<1$ we have $0.5-0.5r+b>0$.  Hence equation \eqref{eqn:rate_f} holds. 

\noindent (b) Consider the inequality in Proposition \ref{prop:error_bound_aIR-IG}(b). We have
		\begin{align*
			&\mathrm{GAP}(\bar{x}_{N,i}) \leq
			 \left(\sum_{k = 0}^{N}\gamma_k^r\right)^{-1} \left( 
			  2M_X^2  \gamma_{N}^{r-1}   +2M_f\sum_{k = 0}^{N}\gamma_{k}^r\eta_{k}	+{\gamma_0^rC_F\|\bar x_{0,m}-x_{0,1} \|} \right.\\
			  &\left.  +{\gamma_{{0}}^rC_F\left\|\bar  {x}_{0,i} - \bar {x}_{0,m}\right\| }+\left( \frac{\left(C_F+\eta_{0}C_f\right)^2}{2} + \frac{{(m-i)}C_F\left(C_F + \eta_0 C_f\right)}{m}\right)\sum_{k=0}^{N}\gamma_{k}^{r+1}\right).
		\end{align*}
Substituting $\{\gamma_k\}$ and $\{\eta_k\}$ by their update rules we obtain
		\begin{align*
			&\mathrm{GAP}(\bar{x}_{N,i}) \leq
			\left(\sum_{k=0}^{N}\frac{\gamma_0^r}{(k+1)^{0.5r}}\right)^{-1} \left(\frac{2M_X^2(N+1)^{0.5(1-r)}}{\gamma_0^{1-r}}   \right.\\
			&\left.+\sum_{k=0}^{N}\frac{2M_f \eta_0\gamma_0^r}{(k+1)^{0.5r+b}} +{\gamma_0^rC_F\left(\|\bar x_{0,m}-x_{0,1} \| +\left\|\bar  {x}_{0,i} -\bar  {x}_{0,m}\right\| \right)} \right.\\ 
&\left.+\left( \frac{\left(C_F+\eta_{0}C_f\right)^2}{2} + \frac{{(m-i)}C_F\left(C_F + \eta_0 C_f\right)}{m}\right)\sum_{k=0}^{N}\frac{\gamma_0^{r+1}}{(k+1)^{0.5(1+r)}}\right).
		\end{align*}
Utilizing the bounds in Lemma \ref{lem:harmonic_series_bound} we obtain
	\begin{align*
		&\mathrm{GAP}(\bar{x}_{N,i}) \leq
		\left(\frac{\gamma_0^r(N+1)^{1-0.5r}}{2(1-0.5r)}\right)^{-1} \left(\frac{2M_f \eta_0\gamma_0^r (N+1)^{1-0.5r-b}}{1-0.5r-b}  \right.\\
		&\left.+\frac{2M_X^2(N+1)^{0.5(1-r)}}{\gamma_0^{1-r}} +{\gamma_0^rC_F\left(\|\bar x_{0,m}-x_{0,1} \|  +\left\|\bar  {x}_{0,i} - \bar {x}_{0,m}\right\|\right)}\right.\\
		&\left.+ {\left( \frac{\left(C_F+\eta_{0}C_f\right)^2}{2} + \frac{{(m-i)}C_F\left(C_F + \eta_0 C_f\right)}{m}\right)\frac{\gamma_0^{r+1}(N+1)^{1-0.5(1+r)}}{(1-0.5(1+r))}}\right).
	\end{align*}	
Rearranging the terms we obtain
	\begin{align*} 
		&\mathrm{GAP}(\bar{x}_{N,i})\leq
		\left({2-r}\right)\left(\frac{2M_X^2}{\gamma_0(N+1)^{0.5}}   +\frac{2M_f \eta_0}{(1-0.5r-b)(N+1)^{b}}\right.\\
		&\left.+\frac{{C_F\|\bar x_{0,m}-x_{0,1} \|}  + C_F\left\| \bar {x}_{0,i} -\bar  {x}_{0,m}\right\|}{(N+1)^{1-0.5r}}\right.\\
		&\left.+ {\left( \frac{\left(C_F+\eta_{0}C_f\right)^2}{2} + \frac{{(m-i)}C_F\left(C_F + \eta_0 C_f\right)}{m}\right)\frac{\gamma_0}{(1-0.5(1+r))(N+1)^{0.5}}}\right).
	\end{align*}	
Note that from $b<0.5$ and $0\leq r<1$ we have $1-0.5r \geq b$.  Hence equation \eqref{eqn:rate_gap} holds. }
\end{proof} 
{\begin{remark}[Iteration complexity of Algorithm \ref{alg:IR-IG_avg}] \normalfont
	Consider the rate results presented by relations \eqref{eqn:rate_f} and \eqref{eqn:rate_gap}. Let us choose $r:=0$ and suppose $\gamma_{k} := \frac{(C_F+C_f)^{-1}}{\sqrt{k+1}}$ and $\eta_{k}:=\frac{1}{\sqrt[4]{k+1}}$ for $k\geq 0$. Let $\epsilon>0$ be an arbitrary small scalar such that $f\left(\bar{x}_{N_\epsilon,i}\right) -f^*+\mathrm{GAP}\left(\bar{x}_{N_\epsilon,i}\right)<\epsilon$ for all {$i\in [m]$}. Then,  we obtain the iteration complexity of $N_\epsilon=\mathcal{O}((C_F+C_f)^4\epsilon^{-4})$ for each agent.  Interestingly, this iteration complexity matches the complexity of the proposed method in our earlier work \cite{KaushikYousefianSIOPT2021} for addressing formulation \eqref{prob:initial_problem} in a centralized regime where the information of the objective function $f$ is globally known. Importantly, this indicates that there is no sacrifice in the iteration complexity in addressing the distributed formulation \eqref{prob:initial_problem}. Another {important} observation to make is that the iteration complexity of the proposed distributed method is independent of the number of agents $m$. 
\end{remark}}

\section{Rate analysis  in the solution space}\label{sec:convergence_analysis_strong_convexity}
{In this section we study the convergence rate properties of the proposed method in the solution space. To this end, we compare the sequences generated from Algorithm \ref{alg:IR-IG_avg} with the Tikhonov trajectory (formally introduced in Definition \ref{def:regularized_problem}). Throughout this section, we make the following additional assumption.
\begin{assumption}  \label{assum:initial_problem_3} Consider problem \eqref{prob:initial_problem}. For all  {$i \in [m]$} let {the} component function $f_i:\mathbb{R}^n\rightarrow {\mathbb{R}}$ be continuously differentiable and $\mu_{f_i}-$strongly convex over the set $X$.
\end{assumption}
Note that under this assumption, the equation \eqref{eqn:alg_step_1} in Algorithm \ref{alg:IR-IG_avg} can be written as
		\begin{align}\label{eqn:alg_2_step_1}
		x_{k,{i+1}} & :=  \mathcal{P}_X \left({x_{k,i}}-\gamma_k\left(F_i\left(x_{k,i}\right)+ {\eta_{k}} {\nabla} f_{i}\left( x_{k,i}\right)\right)\right).
	\end{align}
{Next, we comment on the strong convexity parameter of the global objective function.}
\begin{remark}\label{rem:strong_convex_parameter_f}\normalfont
Under Assumption \ref{assum:initial_problem_3}, we note that any function $f_i$ is strongly convex with a parameter $\mu_{\text{min}}\triangleq \min_{{i \in [m]}}\mu_{f_i}$. This also implies that the global function $f(x)\triangleq \sum_{i = 1}^m f_i(x)$ is $m \mu_{\text{min}}-$ strongly convex. Another implication is that under Assumption \ref{assum:initial_prob}(b), (c), and Assumption \ref{assum:initial_problem_3}, problem \eqref{prob:initial_problem} has a unique optimal solution. Throughout this section, we denote the {unique optimal} solution of \eqref{prob:initial_problem} by $x^*$.
\end{remark}
\subsection{Preliminaries}\label{subsec:preliminararies}
Here we provide some preliminary results that will be used {later}. We first introduce the notion of Tikhonov trajectory.
\begin{definition}[Tikhonov trajectory]\label{def:regularized_problem} \normalfont
Let $\{\eta_k\}\subseteq \mathbb{R}_{++}^n$ be a given sequence. Consider a class of regularized VI problems for $k\geq 0$ given by
 \begin{align}\label{prob:regularized_problem}
 	\text{VI}\left( X, \sum_{i = 1}^m \left(F_i + \eta_{k} {\nabla} f_i\right)\right).
 \end{align}
Let $x^*_{\eta_k}$ denote the unique solution to \eqref{prob:regularized_problem}. The sequence  $\{ x^*_{\eta_k}\}$ is called the Tikhonov trajectory associated with problem \eqref{prob:initial_problem}. 
\end{definition}
In the following, we establish the convergence of the Tikhonov trajectory to the unique optimal solution of problem \eqref{prob:initial_problem}. This result will be used later in Theorem \ref{thm:part2}.
\begin{lemma}[Properties of the Tikhonov trajectory]\label{lem:recursive_bound_reg_problem}
	Consider problem \eqref{prob:initial_problem} and Definition \ref{def:regularized_problem}. Let Assumption \ref{assum:initial_prob}(b), (c), and Assumption \ref{assum:initial_problem_3} hold. Then:\\

	\noindent (a) Let $\{\eta_k\}$ be a strictly positive sequence such that ${\text{lim}}_{k\rightarrow \infty}\eta_{{k}}=0$. Then, $\lim_{k\to\infty} x^*_{\eta_{{k}}}$ exists and is equal to $x^*$.
	
	\noindent (b) For any two nonnegative integers $k_1$ and $k_2$, we have $\left\| x^*_{\eta_{{k_2}}} - x^*_{\eta_{{k_1}}}  \right\|\leq \frac{C_f}{m \mu_{\text{min}}}\left|1-\frac{\eta_{{k_2}}}{\eta_{{k_1}}}\right|.$
\end{lemma}   
\begin{proof}
The proof can be done in a similar fashion to the proof of Lemma 4.5 in \cite{KaushikYousefianSIOPT2021}. 
\end{proof}
Next we obtain a recursive bound on an error metric that is characterized by the sequence $\{x_k\}$ in Definition \ref{def:convenient_notation} and the Tikhonov trajectory. This result will be utilized in Theorem \ref{thm:part2}.
\begin{lemma}
	Let  $\{x_k\}$ be given by Definition \ref{def:convenient_notation}. Let Assumption \ref{assum:initial_prob}(b), (c), and Assumption \ref{assum:initial_problem_3} hold. Suppose $\{\gamma_k\}$ and $\{\eta_k\}$ are strictly positive and nonincreasing such that $\gamma_0\eta_0 \mu_{\text{min}}\leq 0.5$. Then for any $k\geq 1$ we have
	\begin{align}\label{eqn:iterative_bound_0}
	&\left\| x_{k+1} - x^*_{\eta_k} \right\|^2 \leq 
 \left(1- \gamma_k\eta_k \mu_{\text{min}}\right)\left\| x_k-x^*_{\eta_{k-1}} \right\|^2 +\frac{1.5C_f^2}{m^2\gamma_k\eta_k\mu_{\text{min}}^3}\left(1-\frac{\eta_{k}}{\eta_{k-1}}\right)^2 +\gamma_k^2\left(C_F+\eta_kC_f\right)^2.
	\end{align}
\end{lemma}}
\begin{proof}
{From Algorithm \ref{alg:IR-IG_avg}, the nonexpansivity of the projection, and}   $x^*_{\eta_k} \in X$, for any {$i \in [m]$} and $k\geq 1$ we have
{\begin{align*} 
		\left\| x_{k,i+1} - x^*_{\eta_k} \right\|^2 &\leq  \left\| x_{k,i} - x^*_{\eta_k} \right\|^2 +\gamma_{k}^2\left\|F_i(x_{k,i}) +\eta_k\nabla f_i\left(x_{k,i}\right)\right\|^2 \\&-  2 \gamma_{k}\left(F_i(x_{k,i}) + \eta_k\nabla f_i\left(x_{k,i}\right)\right)^T\left( x_{k,i} - x^*_{\eta_k}\right).
	\end{align*}}
From the definition of $C_F$ and $C_f$ we have
	\begin{align*} 
		&\left\| x_{k,i+1} - x^*_{\eta_k} \right\|^2 
		\leq \left\| x_{k,i} - x^*_{\eta_k} \right\|^2+\gamma_k^2
		\left(\frac{C_F+\eta_kC_f}{m}\right)^2+   2 \gamma_{k}\left(F_i(x_{k,i}) + \eta_k\nabla f_i\left(x_{k,i}\right)\right)^T\left(x^*_{\eta_k}- x_{k,i}\right) .
	\end{align*}
From the strong monotonicity of  ${{\nabla}} f_i  $ and the monotonicity of $ F_i $ we can write
\begin{align*}
2 \gamma_{k}\left(F_i(x_{k,i}) + \eta_k\nabla f_i\left(x_{k,i}\right)\right)^T\left(x^*_{\eta_k}- x_{k,i}\right) 
& \leq 2\gamma_{k}F_i\left(x^*_{\eta_k}\right)^T\left(x^*_{\eta_k}- x_{k,i}\right) \\
&+{2\gamma_{k}\eta_k\left(\nabla f_i\left(x^*_{\eta_k}\right)^T\left(x^*_{\eta_k}- x_{k,i}\right) -  \mu_{\text{min}}\left\|x^*_{\eta_k}- x_{k,i}\right\|^2\right)}.
\end{align*}
From the preceding two relations we obtain
\begin{align*}
	\left\| x_{k,i+1} - x^*_{\eta_k} \right\|^2 
	&\leq  \left(1- 2\gamma_k\eta_k \mu_{\text{min}}\right)\left\| x_{k,i} - x^*_{\eta_k} \right\|^2+\gamma_k^2\left(\frac{C_F+\eta_kC_f}{m}\right)^2 +2\gamma_{k}F_i\left(x^*_{\eta_k}\right)^T\left(x^*_{\eta_k}- x_{k,i}\right)\nonumber \\
	& +2\gamma_{k}\eta_k\nabla f_i\left(x^*_{\eta_k}\right)^T\left(x^*_{\eta_k}- x_{k,i}\right) .
\end{align*}
Adding and subtracting $2\gamma_kF_i (x^*_{\eta_{{k}}} )^Tx_k + 2\gamma_k\eta_k\nabla f_i (x^*_{\eta_{{k}}} )^Tx_k$ in the previous relation we get
{\begin{align*}
	\left\| x_{k,i+1} - x^*_{\eta_k} \right\|^2 &\leq  \left(1- 2\gamma_k\eta_k \mu_{\text{min}}\right)\left\| x_{k,i} - x^*_{\eta_k} \right\|^2 +\gamma_k^2 \left(\frac{C_F+\eta_kC_f}{m}\right)^2 +2\gamma_{k}F_i\left(x^*_{\eta_k}\right)^T\left(x^*_{\eta_k}- x_{k}\right) \nonumber\\ 
	&+2\gamma_{k}\left(\eta_k\nabla f_i\left(x^*_{\eta_k}\right)^T\left(x^*_{\eta_k}- x_{k}\right) + \left|F_i\left(x^*_{\eta_k}\right)^T\left(x_{k}- x_{k,i}\right)\right|\right) \\
	& +2\gamma_{k}  \eta_k \left|\nabla f_i\left(x^*_{\eta_k}\right)^T\left(x_{k}- x_{k,i}\right) \right| .
\end{align*}}
Employing the Cauchy-Schwarz inequality we obtain
\begin{align*}
\left\| x_{k,i+1} - x^*_{\eta_k} \right\|^2 &\leq  \left(1- 2\gamma_k\eta_k \mu_{\text{min}}\right)\left\| x_{k,i} - x^*_{\eta_k} \right\|^2+\gamma_k^2 \left(\frac{C_F+\eta_kC_f}{m}\right)^2+2\gamma_{k}F_i\left(x^*_{\eta_k}\right)^T\left(x^*_{\eta_k}- x_{k}\right)\nonumber \\ 
& +2\gamma_{k}\eta_k\nabla f_i\left(x^*_{\eta_k}\right)^T\left(x^*_{\eta_k}- x_{k}\right)+\frac{2\gamma_{k}\left(C_F+\eta_{{k}}C_f\right)}{m}\left\|x_{k}- x_{k,i}\right\| .
\end{align*}
Next we take summations over {$i \in [m]$} from both sides. {Recall  $f(x)\triangleq\sum_{i = 1}^mf_i(x)$ and $F(x)\triangleq\sum_{i = 1}^mF_i(x)$.} Using  Definition \ref{def:convenient_notation} for $x_{k,1}$ and recalling  $1-2\gamma_k\eta_k\mu_{\text{min}}< 1$  we have
\begin{align}\label{eqn:iterative_bound_3}
	&\sum_{i = 1}^m\left\| x_{k,i+1} - x^*_{\eta_k} \right\|^2 \leq  \left(1- 2\gamma_k\eta_k \mu_{\text{min}}\right)\left\| x_{k} - x^*_{\eta_k} \right\|^2  +  \sum_{i = 2}^m\left\| x_{k,i} - x^*_{\eta_k} \right\|^2 +\gamma_k^2 
	\frac{\left(C_F+\eta_kC_f\right)^2}{m}\nonumber\\
	&+\frac{2\gamma_{k}\left(C_F+\eta_{{k}}C_f\right)}{m}\sum_{i = 1}^m \left\|x_{k}- x_{k,i}\right\|+2\gamma_{k}F\left(x^*_{\eta_k}\right)^T\left(x^*_{\eta_k}- x_{k}\right)  +2\gamma_{k}\eta_k\nabla f\left(x^*_{\eta_k}\right)^T\left(x^*_{\eta_k}- x_{k}\right).
\end{align}
{From Lemma \ref{lem:neighbor_agents}, $\| x_k - x_{k,i} \|\leq (i-1)\gamma_{k}\left({C_F}+\eta_{k}C_f\right)/m$} for all {$i\in[ m]$}. Invoking this relation and Definition \ref{def:convenient_notation} we obtain
\begin{align*}
	\left\| x_{k+1} - x^*_{\eta_k} \right\|^2 &\leq \left(1- 2\gamma_k\eta_k\mu_{\text{min}}\right)\left\| x_{k} - x^*_{\eta_k} \right\|^2   +\gamma_k^2\left(C_F+\eta_kC_f\right)^2 \\&+ 2\gamma_k\underbrace{\left(F\left(x^*_{\eta_k}\right)+\eta_k\nabla f\left(x^*_{\eta_k}\right)\right)^T\left(x^*_{\eta_k}- x_{k}\right)}_{\text{term 1}}.
\end{align*}
From Definition \ref{def:regularized_problem}, $x^*_{\eta_k}$ is the solution to problem \eqref{prob:regularized_problem}.  Recalling  $x_k \in X$,  we have: term 1 $\leq 0$. We obtain 
\begin{align}\label{eqn:iterative_bound_4}
	\left\| x_{k+1} - x^*_{\eta_k} \right\|^2 &\leq 
	 \left(1- 2\gamma_k\eta_k\mu_{\text{min}}\right) \left\| x_{k} - x^*_{\eta_k} \right\|^2 +\gamma_k^2\left(C_F+\eta_kC_f\right)^2.
\end{align}
{Next, consider the term $\| x_k-x^*_{\eta_k} \|^2$ as follows.}
{\begin{align*}
\left\| x_k-x^*_{\eta_k} \right\|^2 &= \left\| x_k-x^*_{\eta_{k-1}}\right\|^2+ \left\|x^*_{\eta_{k-1}}-x^*_{\eta_k} \right\|^2  + \underbrace{2\left(x_k-x^*_{\eta_{k-1}}\right)^T\left(x^*_{\eta_{k-1}}-x^*_{\eta_k}\right)}_{\text{term 2}}{.}
\end{align*}}
{Next, we  bound term 2 by recalling $2a^Tb\leq \|a\|^2/\alpha + \alpha\|b\|^2 $ where $a,b \in \mathbb{R}^n$ and $\alpha > 0$. For $\alpha := {1}/{\gamma_{k}\eta_k\mu_{\text{min}}}$, bounding term 2 in the preceding inequality we obtain}
\begin{align*}
	\left\| x_k-x^*_{\eta_k} \right\|^2 &=
	 \left(1+{\gamma_{k}\eta_k\mu_{\text{min}}}\right)\left\| x_k-x^*_{\eta_{k-1}} \right\|^2+\left(1+\frac{1}{\gamma_{k}\eta_k\mu_{\text{min}}}\right) \left\|x^*_{\eta_{k-1}}-x^*_{\eta_k} \right\|^2 .
\end{align*}
From Lemma \ref{lem:recursive_bound_reg_problem}(b) we obtain
\begin{align}\label{eqn:iterative_bound_5}
	\left\| x_k-x^*_{\eta_k} \right\|^2 &\leq \left(1+{\gamma_{k}\eta_k\mu_{\text{min}}}\right)\left\| x_k-x^*_{\eta_{k-1}} \right\|^2+\left(1+\frac{1}{\gamma_{k}\eta_k\mu_{\text{min}}}\right)\frac{C_f^2}{m^2\mu_{\text{min}}^2}\left(1-\frac{\eta_{k}}{\eta_{k-1}}\right)^2.
\end{align}
From equations \eqref{eqn:iterative_bound_4} and \eqref{eqn:iterative_bound_5} we have
\begin{align*}
	\left\| x_{k+1} - x^*_{\eta_k} \right\|^2 &\leq 
	\left(1- 2\gamma_k\eta_k\mu_{\text{min}}\right) \left(1+{\gamma_{k}\eta_k\mu_{\text{min}}}\right)\left\| x_k-x^*_{\eta_{k-1}} \right\|^2\\
	&+\left(1+\frac{1}{\gamma_{k}\eta_k\mu_{\text{min}}}\right)\frac{C_f^2}{m^2\mu_{\text{min}}^2}\left(1-\frac{\eta_{k}}{\eta_{k-1}}\right)^2 +\gamma_k^2\left(C_F+\eta_kC_f\right)^2.
\end{align*}
{Using $0<\gamma_k\eta_k\mu_{\text{min}}\leq 0.5$ we have the desired result.}
\end{proof}
\subsection{Convergence analysis}\label{subsec:convergence}
{In this section, our goal is to derive a {non-asymptotic} convergence rate statement that relates the generated sequences by Algorithm \ref{alg:IR-IG_avg} to the Tikhonov trajectory. We begin with providing a class of sequences for the stepsize and the regularization parameter and prove some properties for them that will be used in the analysis. 
\begin{definition}[{Stepsize and regularization parameter}]\label{def:sequences_tikh_section}
Let $\gamma_k:=\frac{\gamma}{(k+\Gamma)^a}$ and $\eta_k:=\frac{\eta}{(k+\Gamma)^b}$ for all $k\geq 0$ where $\gamma,\eta,\Gamma, a$ and $b$ are strictly positive scalars. Let $a>b$, $a+b<1$, and $3a+b<2$. Assume that $\Gamma \geq 1$ and it is sufficiently large such that $\Gamma^{a+b}\geq 2\gamma\eta\mu_{\text{min}}$ and $\Gamma^{1-a-b}\geq  \frac{4}{\gamma\eta\mu_{\text{min}}}$. 
\end{definition}
\begin{lemma}\label{lem:stepsize_reg_parameter_sequence_a} Consider Definition \ref{def:sequences_tikh_section}. The following results hold.
	
	\noindent(i) $\{\gamma_{k}\}$ and $\{\eta_{k}\}$ are strictly positive and nonincreasing such that $\gamma_0\eta_{{0}}\mu_{\text{min}} \leq  0.5$.
	
	\noindent(ii)  For all integers $k_1$ and $k_2$ such that $k_2\geq k_1\geq 0$ we have $1-\frac{\eta_{k_2}}{\eta_{k_1}}  \leq \frac{k_2-k_1}{k_2+\Gamma}$.
	
		\noindent(iii)  For all $k\geq 1$ we have $ \frac{1}{\gamma_k^3\eta_{k}} \left( 1- \frac{\eta_{{k}}}{\eta_{k-1}} \right)^2  \leq  \frac{1}{\gamma^3\eta\Gamma^{2-3a-b}}$.
	
	\noindent(iv) For all $k\geq 1$ we have $\frac{\gamma_{k-1}}{\eta_{k-1}} \leq \frac{\gamma_{k}}{\eta_{k}}(1+0.5 \gamma_{k}\eta_{k}\mu_{\text{min}})$.
\end{lemma}}
{\begin{proof}
	See Appendix \ref{app:lem:stepsize_reg_parameter_sequence_a}.
\end{proof}}
{The main contribution in this section is presented by Theorem \ref{thm:part2} where we derive agent-specific rates relating the sequences generated by Algorithm \ref{alg:IR-IG_avg} with the Tikhonov trajectory.
\begin{theorem}[Comparison with the Tikhonov trajectory]\label{thm:part2}
Consider problem \eqref{prob:initial_problem}. Let Assumption \ref{assum:initial_prob}(b), (c), and Assumption \ref{assum:initial_problem_3} hold. Consider $\{x_k\}$ and $\{x^*_{\eta_{{k}}}\}$ given in Definitions \ref{def:convenient_notation} and \ref{def:regularized_problem}, respectively. Let the stepsize sequence $\{\gamma_{k}\}$ and the regularization sequence $\{\eta_{{k}}\}$ be given by Definition \ref{def:sequences_tikh_section}. Then for all $ k\geq 0$ and all {$i\in [m]$} we have 
\begin{align*}
 \left\|x_{k+1,i}  - x^*_{\eta_k}\right\|^2 \leq \frac{2(i-1)^2\left(C_F+\eta_0C_f\right)^2\gamma^{2}}{m^2(k+\Gamma+1)^{2a}} +\frac{2\tau B_0\gamma}{\mu_{\text{min}}\eta(k+\Gamma)^{a-b}},
\end{align*}
where $\tau \triangleq \text{max}\left\{ \mu_{\text{min}}\eta\gamma^{-1}B_0^{-1}\Gamma^{a-b}\left\|x_1-x^*_{\eta_0}\right\|^2, 2 \right\}$ and $B_0 \triangleq\frac{1.5C_f^2}{m^2\mu_{\text{min}}^3\gamma^3\eta\Gamma^{2-3a-b}} + \left(C_F+\eta_0C_f\right)^2$.  
\end{theorem}}
\begin{proof}
{Consider  \eqref{eqn:iterative_bound_0}. From Lemma \ref{lem:stepsize_reg_parameter_sequence_a}, for  $k\geq 1$ we have}
\begin{align*}
&\left\| x_{k+1} - x^*_{\eta_k} \right\|^2 \leq \left(1 - \gamma_k \eta_{{k}}\mu_{\text{min}}\right)\left\| x_k - x^*_{\eta_{k-1}} \right\|^2+ \frac{1.5C_f^2\gamma_{k}^2}{m^2\mu_{\text{min}}^3\gamma^3\eta\Gamma^{2-3a-b}} + \gamma_k^2\left(C_F+\eta_0C_f\right)^2.
\end{align*}
Let us define the terms $ v_k {\triangleq} \left\| x_k-x^*_{\eta_{k-1}} \right\|^2$, $\alpha_{k} {\triangleq} \gamma_{k}\eta_{k}\mu_{\text{min}}$, and $\beta_k\triangleq B_0\gamma_k^2$ for $k \geq 1$. Therefore, for all $k \geq 1$ we have
\begin{align}\label{eqn:proposition2_1}
v_{k+1} \leq (1-\alpha_k)v_k + \beta_k.
\end{align}
From Lemma \ref{lem:stepsize_reg_parameter_sequence_a}(iii), for all $k\geq 1$ we have 
{\begin{align}\label{eqn:proposition2_2}
\frac{\beta_{k-1}}{\alpha_{k-1}} \leq  \frac{B_0\gamma_{k}}{\mu_{\text{min}}\eta_{k}}(1+0.5 \gamma_{k}\eta_{k}\mu_{\text{min}}) =\frac{\beta_{k}}{\alpha_{k}}(1+0.5 \alpha_{k}) .
\end{align}}
Next, we show that $v_{k+1} \leq \tau \frac{\beta_{k}}{\alpha_{k}}$ for all $k\geq 0$. We apply induction on $k\geq 0$. Note that this relation holds for $k:=0$ as an implication of the definition of $\tau$. Suppose $v_k \leq \tau\frac{\beta_{k-1}}{\alpha_{k-1}}$ holds for some $k \geq 1$. From \eqref{eqn:proposition2_1} we obtain $v_{k+1} \leq (1-\alpha_k)\frac{\beta_{k-1}}{\alpha_{k-1}} \tau + \beta_k$. Using the upper bound for the right-hand side given by \eqref{eqn:proposition2_2} we have
\begin{align*}
v_{k+1} &\leq \tau(1-\alpha_k)(1+0.5 \alpha_{k})\frac{\beta_{k}}{\alpha_{k}}  + \beta_k	 =  \tau(1-\alpha_{k} + 0.5 \alpha_k -0.5 \alpha_k^2)\frac{\beta_{k}}{\alpha_{k}}  + \beta_k
 \\&= \tau\frac{\beta_{k}}{\alpha_{k}}-\tau(1-0.5)\beta_{k} -0.5\tau\alpha_{k}\beta_{k} +\beta_{k}\leq\tau\frac{\beta_{k}}{\alpha_{k}} +  (1-0.5\tau)\beta_{k}.
\end{align*}
From the definition of $\tau$, $\tau \geq 2$ implying that $1-0.5\tau\leq 0$. This completes the proof of induction. Recall from Lemma \ref{lem:neighbor_agents} that we have $\| x_k - x_{k,i} \|\leq (i-1)\gamma_{k}\left(C_F+\eta_{k}C_f\right)/m$ for any {$i \in [m]$}.  {For all $k \geq 0$ and  $i\in[m]$ we have}
\begin{align*}
& \left\|x_{k+1,i}  - x^*_{\eta_k}\right\|^2 \leq 2 \left\| x_{k+1,i} -x_{k+1}\right\|^2+ 2\left\|x_{k+1} - x^*_{\eta_k}\right\|^2\\
 &\leq \frac{2(i-1)^2\left(C_F+\eta_0C_f\right)^2\gamma_{k+1}^2}{m^2} +\frac{2\tau B_0\gamma_{k}}{\mu_{\text{min}}\eta_k} =\frac{2(i-1)^2\left(C_F+\eta_0C_f\right)^2\gamma^2}{m^2(k+\Gamma+1)^{2a}} +\frac{2\tau B_0\gamma}{\mu_{\text{min}}\eta(k+\Gamma)^{a-b}}.
\end{align*}
Hence the proof is completed. 
\end{proof}


\section{Numerical results} \label{sec:numerical_implementation}
In this section we present the implementation results of Algorithm \ref{alg:IR-IG_avg} in addressing two of the motivating examples discussed earlier in section \ref{sec:introduction}. These include a traffic equilibrium problem and a soft-margin support vector classification problem. We also provide comparisons between the performance of the proposed method with that of other standard IG methods.


\noindent {\bf (i) A traffic equilibrium problem:} For an illustrative example, we consider the transportation network in \cite{CZF2009NCP}. We first describe the network and present the NCP formulation. Then, we implement Algorithm \ref{alg:IR-IG_avg} to solve model \eqref{prob:dist_best_NCP} and compute the best equilibrium.

\begin{wrapfigure}{r}{5.7cm}
\includegraphics[width=5.6cm]{Diagram1.pdf}
\caption{A transportation network with $2$ nodes and $5$ arcs}\label{fig:network}
\vspace*{-4mm}
\end{wrapfigure}
 Consider {a transportation network with the set of nodes $\{n_1,n_2\}$  and {the set of} directed arcs $ \{a_1,a_2,a_3,a_4,a_5\} $. As shown in Figure \ref{fig:network}, arcs $a_1, a_2$, and $a_3$ are directed from {node} $n_1$ to {node} $n_2$, and arcs $a_4$ and $a_5$ {are directed} in the reverse way. Note that $a_1$ and $a_4$ construct a two-way road. The same holds for $a_2$ and $a_5$. We let $d \triangleq [d_1,d_2]^T$ denote the expected travel demand vector where  $d_1$ and $d_2$ correspond to the demand from $n_1$ to $n_2$, and from $n_2$ to $n_1$, respectively. Let the vector $h\triangleq [h_1,\ldots,h_5]^T$ denote the traffic flow on the arcs. The travel cost on each arc is assumed to be a linear function in terms of $h$. More precisely, the travel cost on arc $i$ is equal to $[Ch+q]_i$ where we let the cost matrix $C\in \mathbb{R}^{5\times 5}$ and vector $q\in \mathbb{R}^5$ be given by
 
  { \scriptsize \begin{align}\label{eqn:C_q}
C := \begin{bmatrix}
	0.92& 0&0&5&0\\0&5.92&0&0&5\\0&0&10.92&0&0\\ 2&0&0&10.92&0 \\ 0&1&0&0&15.92
\end{bmatrix}, \qquad q :=  \begin{bmatrix} 1000 \\ 950 \\ 3000 \\ 1000 \\ 1300 \end{bmatrix}.
\end{align}}
\begin{table*}[t]
	\renewcommand\thetable{2}
	\setlength{\tabcolsep}{0pt}
	\centering
	\begin{tabular}{c | c  c  c c}
		{\footnotesize$(\gamma_0,\eta_0) \ $}& {\footnotesize  $(0.1,0.1)$} & {\footnotesize $(0.1,1)$} &  {\footnotesize $(1,0.1)$} & {\footnotesize $(1,1)$}  \\
		\hline\\
		\rotatebox[origin=c]{90}{{\footnotesize $\ln\left(\text{infeasibility}\right)$}}
		&
		\begin{minipage}{.22\textwidth}
			\includegraphics[scale=.125, angle=0]{Numerical_Expts_May_2021/Traffic_equilibrium_problem/1GAP_1.pdf}
		\end{minipage}
		&
		\begin{minipage}{.22\textwidth}
			\includegraphics[scale=.125, angle=0]{Numerical_Expts_May_2021/Traffic_equilibrium_problem/1GAP_2.pdf}
		\end{minipage}
		&
		\begin{minipage}{.22\textwidth}
			\includegraphics[scale=.125, angle=0]{Numerical_Expts_May_2021/Traffic_equilibrium_problem/1GAP_3.pdf}
		\end{minipage}
		&
		\begin{minipage}{.22\textwidth}
			\includegraphics[scale=.125, angle=0]{Numerical_Expts_May_2021/Traffic_equilibrium_problem/1GAP_4.pdf}
		\end{minipage}
		\\
		\hbox{}& & & \\
		\hline\\
		\rotatebox[origin=c]{90}{{\footnotesize$\ln\left(\text{objective}\right)$}}
		&
		\begin{minipage}{.22\textwidth}
			\includegraphics[scale=.125, angle=0]{Numerical_Expts_May_2021/Traffic_equilibrium_problem/1individual_objective_1.pdf}
		\end{minipage}
		&
		\begin{minipage}{.22\textwidth}
			\includegraphics[scale=.125, angle=0]{Numerical_Expts_May_2021/Traffic_equilibrium_problem/1individual_objective_2.pdf}
		\end{minipage}
		&
		\begin{minipage}{.22\textwidth}
			\includegraphics[scale=.125, angle=0]{Numerical_Expts_May_2021/Traffic_equilibrium_problem/1individual_objective_3.pdf}
		\end{minipage}
		&
		\begin{minipage}{.22\textwidth}
			\includegraphics[scale=.125, angle=0]{Numerical_Expts_May_2021/Traffic_equilibrium_problem/1individual_objective_4.pdf}
		\end{minipage}
	\end{tabular}
	\captionof{figure}{Performance of Algorithm \ref{alg:IR-IG_avg} in terms of infeasibility and the objective function value for finding the best equilibrium in the transportation network problem}\label{fig:TEP}
\end{table*} 
We note that the matrix $C$ is positive semidefinite. The diagonal values of $C$ are rounded by two decimal places for the ease of presentation. Intuitively speaking, the structure of $C$ implies that the cost of each arc in a two-way road depends on the flows on the both directions. Let $u\triangleq [u_1,u_2]^T$ denote the (unknown) vector of minimum travel costs between the origin-destination (OD) pairs, i.e., $u_1$ denotes the minimum travel cost from $n_1$ to $n_2$, and $u_2$ denotes the minimum travel cost from $n_2$ to $n_1$. The Wardrop user equilibrium principle represents the path choice behavior of the users based on the following rationale: (i) For any OD pair among all possible arcs, users tend to choose the {arc(s)} with the minimum cost. (ii) For any OD pair, the arc(s) that have the minimum cost will have positive flows and will have equal costs. (iii) For any OD pair, arcs with higher costs than the minimum value will have no flows. Mathematically the Wardrop's principle can be characterized as 
\begin{align}\label{eqn:wardrop}
	0 \leq C h + q - B^T u \perp h \geq 0, \quad 0 \leq Bh - d \perp u \geq 0,
\end{align} 
where $B \in \mathbb{R}^{2\times 5}$ denotes the (OD pair, arc)-incidence matrix given as {\scriptsize $B := \begin{bmatrix}1&1&1&0&0\\
0&0&0&1&1\end{bmatrix}$}. Throughout, we assume that the demand vector $d$ and the cost vector $q$ are subject to uncertainties. Let us define decision vector $x \in \mathbb{R}^7$, random variable $\xi \in \mathbb{R}^{10}$, and stochastic mapping $F(\bullet,\xi):\mathbb{R}^7\to\mathbb{R}^7$ as 
\begin{align*}
x\triangleq {\scriptsize \begin{bmatrix}h \\ u\end{bmatrix}}, \qquad \xi \triangleq {\scriptsize \begin{bmatrix}\tilde d \\ \tilde q\end{bmatrix}}, \qquad F(x,\xi)\triangleq{\scriptsize \begin{bmatrix}C &-B^T \\ B& 0\end{bmatrix}\begin{bmatrix}h \\ u\end{bmatrix} +\begin{bmatrix}\tilde q \\ -\tilde d\end{bmatrix}}.
\end{align*}
Then from section \ref{sec:introduction}, the Wardrop equation \eqref{eqn:wardrop} can be characterized as $\mbox{VI}\left(\mathbb{R}_+^7,\mathbb{E}[F(\bullet,\xi)]\right)$. Notably due to positive semidefinite property of $C$, the mapping $\mathbb{E}[F(\bullet,\xi)]$ is merely monotone. Consequently, the aforementioned VI may have multiple equilibria. Among them, we seek to find the best equilibrium with respect to a welfare function $f$ defined as the expected total travel time over the network by all users, i.e., $f(x)\triangleq \mathbb{E}\left[(Ch+\tilde q)^T\mathbf{1}_5\right]$ where $\mathbf{1}_5 {\in \mathbb{R}^n}$ denotes a vector with all unit elements. 
}


{\noindent \textit{Set-up}: For this experiment, we assume that $\tilde d_1\sim \mathcal{N}(210,10)$, $\tilde d_2\sim \mathcal{N}(120,10)$. Also for $i=1,\ldots,5$ we let $\tilde q_i$ be normally distributed with the mean equal to $q_i$ and the standard deviation of $300$, where the vector $q$ is given by \eqref{eqn:C_q}. Following the formulation \eqref{prob:dist_best_NCP} we generate $1000$ samples for each parameter and distribute the data equally among $10$ agents. We let $\gamma_k:=\frac{\gamma_0}{\sqrt{k+1}}$ and $\eta_k:=\frac{\eta_0}{(k+1)^{0.25}}$ and consider different values for the initial stepsize $\gamma_{0}$ and the initial regularization parameter $\eta_{0}$. The results are as shown in Figure \ref{fig:TEP}. We use standard averaging by assuming that $r=0$. Notably for quantifying the infeasibility, we consider the metric $\phi(x) \triangleq \left\|\max\{0,-x\}\right\|^2 + \left\|\max\{0,-F(x)\}\right\|^2 + \left|x^TF(x)\right|$, where $F(x)\triangleq \sum_{i=1}^mF_i(x)$ and $F_i(x)\triangleq \sum_{\ell \in \mathcal{S}_i}F(x,\xi_\ell)$. Note that $\phi(x) =0$ if and only if $0\leq x\perp F(x) \geq 0$. We choose this metric over the dual gap function employed earlier in the analysis because in this particular example, the dual gap function becomes infinity at some of the evaluations of the generated iterates. This is due to the unboundedness of the set $X:=\mathbb{R}^n_+$. Unlike the dual gap function, $\phi(x)$ stays bounded and is more suitable to plot.

{\noindent \textit{Insights}:} In Figure \ref{fig:TEP} we observe that in all four different settings the infeasibility metric decreases as the algorithm proceeds. This indeed implies that the generated iterates by the agents tend to satisfy the {NCP constraints} with an increasing accuracy. In terms of the suboptimality metric we observe that the each agent's objective value  becomes more and more stable over time. Intuitively this implies that the agents asymptotically reach to an equilibrium.  We should note that although the function $f$ is minimized, it is minimized only over the set of equilibria. The fact that the objective values in Figure \ref{fig:TEP} are not necessarily decreasing is mainly because of the impact of feasibility violation of the iterates with respect to the {NCP constraints} throughout the implementations. As evidenced, pair-IG performs with much robustness to the choice of the initial values of $\gamma_0$ and $\eta_0$. }

\begin{table*}[t]
	\renewcommand\thetable{3}
	\centering
	\begin{tabular}{c|c|c}
	
		$ {|S|} \backslash {n}$  & 50 & 100 \\
	
		\hline
		100&\begin{minipage}{7.5cm}\centering\includegraphics[width=0.5\textwidth]{f_11.pdf}\includegraphics[width=0.5\textwidth]{infeas_11.pdf}\end{minipage}
		&\begin{minipage}{7.5cm}\centering\includegraphics[width=0.5\textwidth]{f_12.pdf}\includegraphics[width=0.5\textwidth]{infeas_12.pdf}\end{minipage}
		\\\hline 200 &\begin{minipage}{7.5cm}\centering\includegraphics[width=0.5\textwidth]{f_21.pdf}\includegraphics[width=0.5\textwidth]{infeas_21.pdf}\end{minipage}
		&\begin{minipage}{7.5cm}\centering\includegraphics[width=0.5\textwidth]{f_22.pdf}\includegraphics[width=0.5\textwidth]{infeas_22.pdf}\end{minipage}
		\\\hline 500 &\begin{minipage}{7.5cm}\centering\includegraphics[width=0.5\textwidth]{f_31.pdf}\includegraphics[width=0.5\textwidth]{infeas_31.pdf}\end{minipage}
		&\begin{minipage}{7.5cm}\centering\includegraphics[width=0.5\textwidth]{f_32.pdf}\includegraphics[width=0.5\textwidth]{infeas_32.pdf}\end{minipage}\\\hline
	\end{tabular}
\captionof{figure}{Comparison of Algorithm \ref{alg:IR-IG_avg} with standard IG methods in solving an SVM model}
	\label{fig:SVM}
\end{table*}
\noindent {\bf (ii) Distributed support vector classification:} Consider a distributed soft-margin SVM problem described as follows.  Let a dataset be given as $\mathcal D \triangleq \{(u_j,v_j) \mid u_j\in \mathbb{R}^n, \ v_j \in \{-1,+1\}, \hbox{ for } j \in S\}$ where $u_j$ and $v_j$ denote the attributes and the binary label of the $j^{\text{th}}$ data point, respectively, and $S$ denotes the index set.  Let the data be distributed among $m$ agents by defining $\mathcal{D}_i $ such that $\cup_{i=1}^m \mathcal{D}_i =\mathcal D$. Let $S_i$ denote the index set corresponding to agent $i$ such that $\sum_{i=1}^{m}|S_i|=|S|$.  Then given $\lambda >0$, the distributed SVM problem is given as
	{\begin{equation}\label{prob:SVM}
			\begin{aligned}
	   & \underset{w,b,z}{\text{minimize}}&& \sum\nolimits_{i = 1}^m\left(\frac{1}{2m}\|w \|^2 + \frac{1}{\lambda}\sum_{j \in S_i} z_j\right) \\
	& \text{subject to}  &&   v_j(w^T{u}_j + b)   \geq 1 - z_j,  \  \hbox{for } j \in S_i \text{ and }i \in [m], \\
	&	&& z_j \geq 0,    \qquad \qquad \qquad \ \ \ \  \text{for } j \in S_i\hbox{ and  }i \in [m].
		\end{aligned}
	\end{equation}}
To solve problem \eqref{prob:SVM} using Algorithm \ref{alg:IR-IG_avg}, we apply Lemma \ref{Lem:nonlinear_const_DO_problem} by casting \eqref{prob:SVM} as model \eqref{prob:subclass_constrained_opt}. We also implement some of the well-known existing IG schemes, namely projected IG, proximal IAG, and SAGA.  Unlike Algorithm \ref{alg:IR-IG_avg}, these schemes require a projection step onto the constraint set.  To compute the projections we use the Gurobi-Python solver.

{\noindent \textit{Set-up}: We consider $20$ agents and assume that $\lambda := 10$. We let $\gamma_k:=\frac{1}{\sqrt{k+1}}$ and $\eta_k:=\frac{1}{(k+1)^{0.25}}$ in Algorithm \ref{alg:IR-IG_avg}.  We use identical initial stepsizes in all the methods. We provide the comparisons with respect to the runtime and report the performance of each scheme over $200$ seconds.  We use a synthetic dataset with different values for $n$ and $|S|$.  The suboptimality is characterized in terms of the global objective and the infeasibility is measured by quantifying the violation of constraints of problem \eqref{prob:SVM} aligned with ideas in Lemma \ref{Lem:nonlinear_const_DO_problem}.

\noindent \textit{Insights}: In Figure \ref{fig:SVM}, we observe that with an increase in the dimension of the solution space, i.e., $n$, or the size of the training dataset, i.e., $|S|$, the projection evaluations in the standard IG schemes take longer and consequently, the performance of the IG methods is deteriorated in large scale settings. However,  utilizing the reformulation in Lemma \ref{Lem:nonlinear_const_DO_problem}, the proposed method does not require any projection operations for addressing problem \eqref{prob:SVM}.  As such, the performance of Algorithm \ref{alg:IR-IG_avg} does not get affected severely with the increase in $n$ or $|S|$.  Note that in Figure \ref{fig:SVM}, the reason that the IG schemes do not show any updates for $ |S| = 200$ and $|S|=500$ beyond a time threshold is because of the interruption in their last
update when the method reaches the $200$ seconds time limit.}

\section{Conclusion}\label{sec:concludingrem}
{We {introduce} a new unifying formulation for distributed constrained optimization where the constraint set is characterized as the solution set of a merely monotone variational inequality problem.  We show that the proposed formulation can capture several important classes of multi-agent constrained optimization problems including those subject to complementarity, nonlinear equality, or equilibrium constraints. We develop an iteratively regularized incremental gradient method where at each iteration agents communicate over a cycle graph to update their iterates using their local information about the objective function and the mapping.  We derive new iteration complexity bounds in terms of the global objective function and a suitably defined infeasibility metric. To analyze the convergence in the solution space,  we also provide non-asymptotic agent-wise convergence rate statements that relate the iterate of each agent with that of the Tikhonov trajectory.  We validate the theoretical results on an illustrative transportation network problem and also, on a distributed support vector classification problem.  A potential, yet challenging, future direction to our research is to extend the current  algorithm to address distributed optimization problems with VI constraints requiring weaker assumptions on the network {and} with possibly asynchronous communications. }
 
