ARock: an Algorithmic Framework for Asynchronous Parallel Coordinate Updates
Zhimin Peng
Yangyang Xu
Ming Yan
Wotao Yin
%\input{asyn_t_abstract Finding a fixed point to a nonexpansive operator, i.e., $x^*=Tx^*$, abstracts many problems in numerical linear algebra, optimization, and other areas of scientific computing. To solve fixed-point problems, we propose ARock, an algorithmic framework in which multiple agents (machines, processors, or cores) update $x$ in an asynchronous parallel fashion. Asynchrony is crucial to parallel computing since it reduces synchronization wait, relaxes communication bottleneck, and thus speeds up computing significantly. At each step of ARock, an agent updates a randomly selected coordinate $x_i$ based on possibly out-of-date information on $x$. The agents share $x$ through either global memory or communication. If writing $x_i$ is atomic, the agents can read and write $x$ without memory locks. Theoretically, we show that if the nonexpansive operator $T$ has a fixed point, then with probability one, ARock generates a sequence that converges to a fixed points of $T$. Our conditions on $T$ and step sizes are weaker than comparable work. Linear convergence is also obtained. We propose special cases of ARock for linear systems, convex optimization, machine learning, as well as distributed and decentralized consensus problems. Numerical experiments of solving sparse logistic regression problems are presented.
asyn_theory_MP.tex
%
\documentclass[smallextended,final]{svjour3}     % onecolumn (ditto)
\usepackage[usenames]{xcolor}
\usepackage[letterpaper,top=2in, bottom=1.5in, left=1in, right=1in]{geometry}

\smartqed % make Springer's \qed better placed
\newcommand\numberthis{\addtocounter{equation}{1}\tag{\theequation}} % add a number to equations between \begin{algin*} \end{align*}
%\usepackage{amsmath}

\input{macros}

% --- additional macros ---
\usepackage[pdftex]{graphicx}
\usepackage{caption}
\usepackage{multirow}
\usepackage{subcaption}
\usepackage{booktabs}
\usepackage[ruled,vlined]{algorithm2e}
\usepackage{amssymb}
\usepackage[backref,colorlinks=true, linkcolor=blue, citecolor=blue, urlcolor=blue, pdfborder={0 0 0}]{hyperref}

%%%%%%%%% suppress the %%%%%%%%%%%%%%%%%%%%%%
% ! Package amsmath Error: Multiple \label's: label 'important' will be lost.
% error message.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\makeatletter
\def\label@in@display{\gdef\df@label}
\makeatother


%% title
%\title{Parallel Asynchronous Methods: Theory, Algorithms and Applications}
\title{ARock: an Algorithmic Framework for Asynchronous Parallel Coordinate Updates}

%% author
\author{Zhimin Peng \and
        Yangyang Xu \and
        Ming Yan \and
        Wotao Yin}

\institute{Z. Peng \and W. Yin\at
              Department of Mathematics, University of California, Los Angeles, CA 90095, USA\\
              \email{zhimin.peng / wotaoyin@math.ucla.edu}\\
           Y. Xu\at Institute for Mathematics and its Application, University of Minnesota, Minneapolis, MN 55455, USA\\ \email{yangyang@ima.umn.edu} \\
            M. Yan\at Department of Computational Mathematics, Science and Engineering, Michigan State University, East Lansing, MI 48824, USA \\ \email{yanm@math.msu.edu}}
\date{\today}
\journalname{Report}  % Springer stuff


%% start the document


\begin{document}

\maketitle

\begin{abstract}
%\input{asyn_t_abstract}
Finding a fixed point to a nonexpansive operator, i.e., $x^*=Tx^*$, abstracts many problems in numerical linear algebra, optimization, and other areas of scientific computing. To solve fixed-point problems, we propose ARock, an algorithmic framework  in which multiple agents (machines, processors, or cores) update  $x$ in an  asynchronous parallel fashion. Asynchrony is crucial to parallel computing since it reduces  synchronization wait, relaxes communication bottleneck, and thus speeds up computing significantly. At each step of ARock, an  agent updates a randomly selected coordinate $x_i$ based on possibly out-of-date information on $x$. The agents share $x$  through either global memory or communication. If writing $x_i$ is atomic, the agents can read and write $x$ without memory locks.

Theoretically, we show that if the nonexpansive operator $T$ has a fixed point, then with probability one,  ARock generates a sequence that converges to a fixed points of $T$. Our conditions on $T$ and  step sizes are  weaker than comparable work. Linear convergence  is also obtained.

We propose special cases of ARock for linear systems, convex optimization, machine learning, as well as distributed and decentralized consensus problems. Numerical experiments of solving sparse logistic regression problems are presented.
\end{abstract}
\vspace{-4mm}
\tableofcontents

\thispagestyle{plain}
\markboth{Z. Peng, Y. Xu, M. Yan, and W. Yin}{ARock: Async-Parallel Coordinate Updates}

%\input{asyn_t_quick_intro}
\section{Introduction}
Technological advances in data gathering and storage have led to a rapid proliferation of big data in diverse areas such as climate studies, cosmology, medicine, the Internet, and engineering \cite{house2014big}. The data involved in many of these modern applications % is too large to be processed on a single workstation.
are large and grow quickly. Therefore, parallel computational approaches are needed. This paper introduces a new approach to asynchronous parallel computing with convergence guarantees.

In a synchronous(sync) parallel iterative algorithm, the  agents must wait for the slowest agent to finish an iteration before they can all proceed to the next one (Figure \ref{fig:async_vs_sync:a}). Hence, the slowest agent  may cripple the system. In contract,  the agents in an asynchronous(async) parallel iterative algorithm run continuously with little idling (Figure \ref{fig:async_vs_sync:b}). However, the iterations are disordered, and an agent may carry out an iteration without the newest information from other agents.

\begin{figure}[!h]
        \centering
       \begin{subfigure}[b]{0.4\textwidth}
                \includegraphics[width=\textwidth]{./figs/syn-simple}
                \caption{Sync-parallel computing\label{fig:async_vs_sync:a}}
        \end{subfigure}
        ~ %add desired spacing between images, e. g. ~, \quad, \qquad, \hfill etc.
          %(or a blank line to force the subfigure onto a new line)
        \begin{subfigure}[b]{0.4\textwidth}
                \includegraphics[width=\textwidth]{./figs/asyn-simple}
                \caption{Async-parallel computing\label{fig:async_vs_sync:b}}
        \end{subfigure}
        \vspace{-3mm}
        \caption{Sync-parallel computing (left) versus async-parallel computing (right).
}\label{fig:async_vs_sync}
\vspace{-3mm}
\end{figure}


Asynchrony has other advantages \cite{bertsekas1991some}: the system is more tolerant to computing faults and communication glitches; %since an asynchronous algorithm can stay running after some agents have failed;
it is also easy to incorporate new agents.


On the other hand, it is more difficult to analyze asynchronous algorithms and ensure their convergence. It becomes impossible to find a sequence of  iterates that one completely determines the next. Nonetheless, we let any update be a new iteration and propose an async-parallel algorithm (ARock) for the generic fixed-point iteration. It converges if  the fixed-point operator is nonexpansive (Def. \ref{def:quasi_lip}) and has a fixed point.

Let $\cH_1,\ldots,\cH_m$ be Hilbert spaces and  $\cH:=\cH_1\times\cdots\times \cH_m$ be their Cartesian product. For a \emph{nonexpansive operator}  $T:\cH\to\cH$, our problem is to
\begin{equation}\label{eqn:fix_point_set}
\text{find } x^* \in \cH \qquad \text{ such that }  \qquad x^*=Tx^* .
\end{equation}
{Finding a fixed point to $T$ is equivalent to finding a zero of $S\equiv I-T,$ denoted by $x^*$ such that $0= Sx^*$. Hereafter, we will use both $S$ and $T$ for convenience.}

Problem \eqref{eqn:fix_point_set} is widely applicable in linear and nonlinear equations, statistical regression, machine learning, convex optimization, and optimal control. %; see Section \ref{sec:application}.
 A generic framework for problem \eqref{eqn:fix_point_set} is the Krasnosel'ski\u i--Mann (KM) iteration~\cite{krasnosel1955two}:
\begin{equation}\label{eqn:km}
x^{k+1} = x^k + \alpha \, (Tx^k-x^k),\quad \mbox{or equivalently,}\quad x^{k+1} = x^k - \alpha  Sx^k,
\end{equation}
where $\alpha \in (0, 1)$ is the step size. % ($\alpha\ge 1$ is permitted \cite[Proposition 5.15]{bauschke2011convex} when $T$ is a $\alpha$-averaged operator).
If  $\Fix T$ --- the set of fixed points of $T$ (zeros of $S$) --- is nonempty, then the sequence $(x^k)_{k\ge 0}$ converges weakly\cut{\footnote{A sequence $(y^k)_{k\ge 0}\subset\cH$ converges weakly to a point $y$ if $\dotp{y^k,z}\to\dotp{y,z}$ for every $z\in\cH$. A sequence $(y^k)_{k\ge 0}\subset\cH$ converges strongly to a point $y$ if $\|y^k-y\|\rightarrow 0$. Strong convergence gives weak convergence, but not vice versa. If $\cH$ has a finite dimension, then they are equivalent. }} to a point in $\Fix T$ and $(Tx^k-x^k)_{k\ge 0}$ converges strongly to 0. The KM iteration generalizes algorithms in convex optimization, linear algebra, differential equations, and monotone inclusions.  Its special cases include the following iterations: alternating projection, gradient descent, projected gradient descent, proximal-point algorithm,
Forward-Backward Splitting (FBS) \cite{passty1979ergodic}, Douglas-Rachford Splitting (DRS)~\cite{lions1979splitting}, %relaxed Peaceman-Rachford splitting (PRS) \cite{lions1979splitting},
 a three-operator splitting~\cite{davis2015three}, and the Alternating Direction Method of Multipliers (ADMM)~\cite{lions1979splitting,glowinski1975approximation}. %which involve multiple nonexpansive operators.

In ARock, a set of  $p$ agents, $p\ge 1$, solve problem~\eqref{eqn:fix_point_set} by updating the coordinates $x_i\in\cH_i$, $i=1,\ldots,m$, in a  random and asynchronous fashion. Algorithm~\ref{alg:asyn_core} describes the framework. Its special forms for several applications are given in Section \ref{sec:application} below.

\begin{algorithm}[H]\label{alg:asyn_core}
%\SetAlgoLined
\SetKwInOut{Input}{Input}\SetKwInOut{Output}{output}
\Input{$x^0\in\cH$,  $K>0$, a distribution $(p_1,\ldots,p_m)>0$ with $\sum_{i=1}^mp_i=1$;}
global iteration counter $k\gets 0$\;
 \While{$k < K$, every agent asynchronously and continuously}{
  select $i_k\in\{1,\ldots,m\}$ with $\mathrm{Prob}(i_k=i)=p_i$\;%from $\{1, ..., m\}$ under a discrete positive distribution\;
  % $\hat{x}^k\gets$ reads the variable $x$\;
  % computes $\Delta x=-\frac{\eta_k}{mp_{i_k}} S_{i_k}\hat{x}^k$; ~(see \eqref{eqn:asyn_update}--\eqref{eqn:inconsist} for definitions)\;
  % adds $(\Delta x)_{i_k}$ to the component $x_{i_k}$ of the variable $x$\;
  perform an update to $x_{i_k}$ according to  \eqref{eqn:asyn_update}\;
  update the global counter $k \leftarrow k+1$\;
 }
 \caption{ARock: a framework for async-parallel coordinate updates}
\end{algorithm}

Whenever an agent updates a coordinate, the global iteration counter $k$ increases by one.
% In this paper, we introduces an async-parallel coordinate update framework for the KM iteration running on a set of $p$ agents (machines, processors, or cores), all of which can read the common variable $x= (x_1, ..., x_m)$ and update any coordinate $x_i\in\cH_i$ of $x$ in an uncoordinated fashion.  We assume that
The $k$th update is applied to $x_{i_k}\in\cH_{i_k}$, where $i_k\in \{1,\ldots,m\}$ is an independent random variable. Each coordinate update  has the form:
\begin{equation}
\label{eqn:asyn_update}
\textstyle x^{k+1} = x^k - \frac{\eta_k}{mp_{i_k}} \, S_{i_k}  \hat{x}^{k},
\end{equation}
where $\eta_k>0$ is a scalar whose range will be set later, $S_{i_k} {x} := (0, ..., 0, (Sx)_{i_k}, 0, ..., 0),$ and $mp_{i_k}$ is used to normalize nonuniform selection probabilities.
In the uniform case, namely,  $p_i \equiv \frac{1}{m}$ for all $i$, we have $mp_{i_k} \equiv 1$, which simplifies the update \eqref{eqn:asyn_update} to
\beq\label{eqn:simple_asyn_update}
x^{k+1} = x^k - \eta_k S_{i_k}  \hat{x}^{k}.
\eeq %The uncoordinated fashion means that while an agent computes $- \frac{\eta_k}{mp_{i_k}} S_{i_k}\hat{x}^{k}$, $x$ is usually updated by other agents. Therefore, $\hat{x}^k$, which denotes a copy of $x$, % that is used to evaluate $S_{i_k}\hat{x}^{k}$,
% is usually different from $x^k$.
Here, the point $\hat x^k$ is what an agent reads from global memory to its local cache and to which $S_{i_k}$ is applied, and $x^k$ denotes the state of $x$ in global memory just before the update \eqref{eqn:asyn_update} is applied. In a sync-parallel algorithm, we have $\hat{x}^k=x^k$, but  in ARock, due to possible  updates  to $x$ by other agents, $\hat x^k$ can be different from $x^k$. This is a key difference between  sync-parallel and async-parallel algorithms. In Subsection \ref{subset:aync_mem} below, we will establish the relationship between $\hat x^k$ and $x^k$ as
\begin{equation}\label{eqn:inconsist}
\textstyle\hat{x}^k= x^k+\sum_{d\in J(k)}(x^{d}-x^{d+1}),
\end{equation}
where $J(k)\subseteq \{k-1,\ldots,k-\tau\}$ and $\tau \in \mathbb{Z}^+$ is the maximum number of other updates  to $x$ during the computation of \eqref{eqn:asyn_update}. Equation \eqref{eqn:inconsist} has appeared in \cite{liu2014asynchronous}.

The update \eqref{eqn:asyn_update} is only computationally worthy if  $S_{i} x$ is much cheaper to compute than   $Sx$. Otherwise, it is more preferable to apply the full KM update \eqref{eqn:km}. In Section \ref{sec:application}, we will present several applications that have the favorable structures for ARock. The recent work \cite{peng2016coordinate} studies coordinate friendly structures more thoroughly.

%We propose the following algorithm to solve the problem \eqref{eqn:fix_point_set}:\\
The convergence of ARock (Algorithm~\ref{alg:asyn_core}) is stated in Theorems~\ref{thm:convergence} and \ref{thm:strongly_monotone}. Here we include a shortened version, leaving detailed bounds to the full theorems:
\begin{theorem}[Global and linear convergence]\label{thm:simple_convergence}
Let $T:\cH\to\cH$ be a nonexpansive operator that has a  fixed point. %$\Fix T$, which is the set of fixed points of $T$.
Let $(x^k)_{k\ge 0}$ be the sequence generated by Algorithm \ref{alg:asyn_core} with  properly bounded step sizes $\eta_k$. Then, with probability one, $(x^k)_{k\ge 0}$  converges weakly to a fixed point of $T$. {This  convergence becomes  strong if $\cH$ has a finite dimension.}

In addition, if $T$ is demicompact (see Definition~\ref{def:monotone} below), then with probability one, $(x^k)_{k\ge 0}$ converges strongly to a fixed point of $T$ .
 %to a random variable $x^*\in \Fix T$.

Furthermore, if $S\equiv I-T$ is quasi-strongly monotone (see Definition~\ref{def:quasi_lip} below), then $T$ has a unique fixed-point $x^*$,  $(x^k)_{k\ge 0}$  converges strongly to $x^*$ with probability one, and $\EE\|x^k-x^*\|^2$ converges to 0 at a linear rate.
\end{theorem}
\cut{\commrw{Do your examples demonstrate the importance of each of these three results (corresponding to three different conditions on the operator T)?}}

In the theorem, the weak convergence result only requires $T$ to be nonexpansive and has a fixed point. In addition, the computation requires: (a) bounded step sizes; (b) random coordinate selection; and (c) a finite maximal delay $\tau$. Assumption (a) is standard, and we will see the bound can be $O(1)$.  Assumption (b) is essential to both the analysis and the numerical performance of our algorithms. Assumption (c) is \emph{not} essential; an infinite delay with a light tail is allowed (but we leave it to future work). {The strong convergence result applies to all the examples in Section \ref{sec:application}, and the linear convergence result applies to Examples \ref{subsec:smooth} and \ref{subsec:nonsmooth}  when the corresponding operator $S$ is quasi-strongly monotone. Step sizes $\eta_k$ are discussed in Remarks \ref{rm:stepsize} and \ref{rmk:linear}.}

\subsection{On random coordinate selection}
ARock employs \emph{random coordinate selection}. This subsection discusses its advantages and disadvantages.

Its main disadvantage is that an agent cannot  caching
the data associated with a coordinate. The variable $x$ and its related data must be either stored  in  global memory or  passed through communication.\cut{An exception arises in the network setting where the coordinates are locally assigned to the agents and the coordinate updates are activated according to independent Poisson processes; then, the coordinate selection is effectively random and independent though no data needs to be moved.} A secondary disadvantage is that pseudo-random number generation takes  time, which becomes relatively significant if each coordinate update is cheap. (The network optimization examples  in Subsections~\ref{sec:dgd} and~\ref{subsec:async_admm_des} are exceptions, where  data  are naturally stored in a  distributed fashion and random coordinate assignments are  the  results of  Poisson processes.) \cut{\commrw{Add the explanation of Poisson Processes here}}

There are several  advantages of random coordinate selection. It realizes the user-specified update frequency $p_i$ for every component $x_i$, $i=1,\ldots,m$, even when different agents have different computing powers and different coordinate updates cost different amounts of computation. Therefore, random assignment ensures load balance.  The algorithm is also fault tolerant in the sense that if one or more agents fail, it will still converge to a fixed-point of $T$. In addition, it has been observed numerically on certain problems~\cite{chang2008coordinate} that random coordinate selection  accelerates convergence. \cut{Furthermore, in~\cite{xu2013block,xu2014globally}, the algorithms with random coordinate selection keep clear of low-quality local solutions in nonconvex optimization, which are difficult to avoid by fixed-order coordinate selection. \commrw{Is the utility of random coordinate updates in nonconvex optimization relevant for the argument that it's an interesting model to analyze in the convex case?}}



%\subsubsection{Sampling order vs updating order}
%Sampling order can be different from the actual updating order. However, if the delay is independent with the random choice of blocks, the actual block of coordinates will still be updated with the same probability. Without loss of generality, our algorithm assumes that the sampling order is the same as the updating order.



%\begin{theorem}[Linear convergence for quasi contractive operator]\label{thm:simple_linear}
%Let $T$ be a quasi contractive operator. Let $(x^k)_{k\geq0}$ be the sequence generated by Algorithm~\ref{alg:asyn_core} with proper bounded constant stepsize. Then
%$\EE\|x^{k}-x^*\|^2\le \|x^0-x^*\|^2 \cdot \rho^k$ with $\rho \in (0, 1)$.
%\end{theorem}


%\input{asyn_t_memory_access}
\subsection{Uncoordinated memory access}\label{subset:aync_mem}
In ARock, since multiple agents simultaneously read and update $x$ in global memory, $\hat x^k$ --- the result of  $x$ that is read from global memory by an agent to its local cache for computation --- may not equal $x^j$ for any $j\le k$, that is, $\hat x^k$ may never be consistent with a state of $x$ in global memory. This is known as \emph{inconsistent read}. In contrast, \emph{consistent read} means that $\hat x^k = x^j$ for some $j \le k$, i.e., $\hat x^k$ is consistent with a state of $x$ that existed in global memory.

We illustrate inconsistent read and consistent read in the following example, which is depicted in Figure \ref{fig:consist_vs_inconsistent}. Consider $x=[x_1,x_2, x_3, x_4]^T\in\RR^4$ and $x^0=[0, 0, 0, 0]^T$ initially, at time $t_0$. Suppose at time $t_1$,  agent 2 updates $x_1$ from 0 to 1, yielding  $x^1=[1, 0, 0, 0]^T$;  then, at time $t_2$, agent 3 updates $x_4$ from 0 to 2, further yielding $x^2=[1, 0, 0, 2]^T$. Suppose that agent 1 starts reading $x$ from the first component $x_1$ at $t_0$. For consistent read (Figure \ref{fig:consist}), agent 1 acquires a memory lock and only releases the lock after finishing reading all of $x_1$, $x_2$, $x_3$, and $x_4$. Therefore, agent 1 will read in $[0, 0, 0, 0]^T$. Inconsistent read, however, allows agent 1 to proceed without a memory lock: agent 1 starts reading $x_1$ at $t_0$ (Figure \ref{fig:inconsist}) and reaches the last component, $x_4$, after $t_2$; since $x_4$ is updated by agent 3 prior to it is read by agent 1, agent 1 has read $[0, 0, 0, 2]^T$, which is different from any of $x^0,x^1$, and $x^2$.
\begin{figure}[!h]
        \centering
       \begin{subfigure}[b]{0.47\textwidth}
                \includegraphics[width=\textwidth]{./figs/consistent_read}
                \caption{Consistent read. While agent 1 reads $x$ in memory, it acquires a global lock.}\label{fig:consist}
        \end{subfigure}
        ~ %add desired spacing between images, e. g. ~, \quad, \qquad, \hfill etc.
          %(or a blank line to force the subfigure onto a new line)
        \begin{subfigure}[b]{0.47\textwidth}
                \includegraphics[width=\textwidth]{./figs/inconsistent_read}
                \caption{Inconsistent read. No lock. Agent 1 reads $(0,0,0,2)^T$, a non-existing state of $x$.}\label{fig:inconsist}
        \end{subfigure}
        \vspace{-3mm}
        \caption{Consistent read versus inconsistent read: A demonstration.}\label{fig:consist_vs_inconsistent}
        \vspace{-3mm}
\end{figure}

Even with inconsistent read, each component is consistent under the \emph{atomic coordinate update} assumption, which will be defined below. Therefore, we can express what has been read in terms of the changes of individual coordinates. In the above example, the first change is $x^1_1-x^0_1=1$, which is added to $x_1$ just before time $t_1$ by agent 2, and the second change is $x^2_4 - x^1_4=2$, added to $x_4$ just before time $t_2$ by agent 3. The  inconsistent read by agent 1, which gives the result $[0, 0, 0, 2]^T$, % misses the first change but includes the second change. Therefore, $[0, 0, 0, 2]^T$
equals  $x^0 +0\times (x^1-x^0)+ 1\times (x^2 - x^1)$.

We have demonstrated that $\hat x^k$ can be  inconsistent, but each of its coordinates is consistent, that is, for each $i$, $\hat x^k_i$ is an ever-existed state of $x_i$ among $x_i^k,\ldots,x_i^{k-\tau}$. Suppose that  $\hat x^k_i=x_i^{\underline{d}}$, where $\underline{d}\in\{k,k-1,\ldots,k-\tau\}$. Therefore, $\hat x^k_i$ can be related to $x^k_i$ through the \emph{interim  changes} applied to $x_i$. Let $J_i(k)\subset\{k-1,\ldots,k-\tau\}$ be the index set of these interim changes. If $J_i(k)\not=\emptyset$, then $\underline{d}=\min\{d\in J_i(k)\}$; otherwise, $\underline{d}=k$. In addition, we have
 $\hat{x}^k_i =x_i^{\underline{d}}= x^k_i+\sum_{d\in J_i(k)}(x_i^{d}-x_i^{d+1})$.
Since the global counter $k$ is increased after each coordinate update, updates to $x_i$ and $x_j$, $i\not=j$, must occur at different $k$'s and thus $J_i(k)\cap J_j(k)=\emptyset,\,\forall i\neq j$. Therefore, by letting $J(k):=\cup_iJ_i(k)\subset\{k-1,\ldots,k-\tau\}$ and noticing $(x_i^{d}-x_i^{d+1})=0$ for $d\in J_j(k)$ where $i\not = j$, we have
$\textstyle \hat{x}^k_i = x^k_i+\sum_{d\in J(k)}(x_i^{d}-x_i^{d+1}), \forall i=1,\ldots,m$,
which is equivalent to~\eqref{eqn:inconsist}.
%
%
%Inconsistent read means that $\hat{x}^k$ is not necessary equals to $x^j$ for some $j\in\{k,k-1,\ldots,k-\tau\}$ but $\hat x^k$ can relate to $(x^k)_{k\geq 0}$ through equation \eqref{eqn:inconsist}.
%
%\emph{Consistent read} means $\hat x^k = x^j$ for some $j \in \{k, k-1, ..., k -\tau \}$, i.e., $\hat x^k$ is consistent with a state of $(x^k)_{k\geq 0}$ existed in memory. This can be implemented with a global memory lock.  This can be achieved as long as the coordinate update is atomic.
%
%We distinguish two different cases for the uncoordinated memory access including \emph{consistent read} and \emph{inconsistent read}.
%Since it takes time for an agent to complete a read--compute--update cycle,  it is possible that some coordinates $x_i$ are updated by other agents. Therefore, in equation \eqref{eqn:asyn_update}, the read copy $\hat{x}^{k}$ is generally different from the to-be-updated copy $x^k$. Under the following assumptions, the relationship between $\hat x^k$ and $x^k$ can be modeled.
Here, we have made two assumptions:
\begin{itemize}
\item \emph{atomic coordinate update}: a coordinate is not further broken to smaller components during an update; they are all updated at once. %Therefore, the components of each $x_i$ are consistently read;
\item \emph{bounded maximal delay} $\tau$: during any update cycle of an agent, $x$ in global memory is updated at most $\tau$ times by other agents.
\end{itemize}
When each coordinate is a single scalar, updating the scalar is a single atomic instruction on most modern hardware, so the first assumption naturally holds, and our algorithm is \emph{lock-free}. The case where a coordinate is a block that includes multiple scalars is discussed in the next subsection. % \ref{atomlock} below.

%In general, under the atomic coordinate update and bounded maximal delay assumptions, reading each  $x_i$ is  \emph{consistent}, that is, there must exist some $0\le \tau_k\le \tau$ so that  $\hat{x}_i^k=x_i^{k-\tau_k}$.%
%Therefore, $\hat{x}^k_i$ and $x^k_i$ can be related through the interim updates applied to $x_i$, specifically,
%\beq\label{atomcoor}\hat{x}^k_i = x^k_i+\sum_{d\in
% J_i(k)}(x_i^{d}-x_i^{d+1}),
%\eeq
%where $J_i(k)\subset\{k-1,\ldots,k-\tau_k\}$ denotes the set of iterations of the interim updates and $0\le \tau_k\le \tau$.
%
%
%\cut{However, the entire $x$ may be read \textit{inconsistently}, causing  $\hat{x}^k\not= x^j$ for all $j\in\{k,k-1,\ldots,k-\tau\}$. We illustrate \emph{inconsistent reading} with the following example: Consider $x=[x_1,x_2,\ldots,x_n]^T\in\RR^n$ and $x^0=[0,0,\ldots,0]^T$ initially. Suppose at $k=0$, an agent updates $x_1$ from 0 to 1, yielding  $x^1=[1,0,\ldots,0]^T$;  then, at $k=1$, another agent updates $x_n$ from 0 to 2, further yielding $x^2=[1,0,\ldots,2]^T$. Suppose a third agent started reading $x$ from the first component $x_1$ just before $k=0$ and reaching the last component $x_n$ after $k=2$; then this agent has actually read $[0,0,\ldots,2]^T$, which is different from $x^0,x^1,x^2$. From this example, we can see that inconsistent reading arises when a simultaneous update is missed.}
%
%Since each update only involves one coordinate in our algorithm, we have $J_i(k)\cap J_j(k)=\emptyset,\,\forall i\neq j$.  By letting $J(k):=\cup_iJ_i(k)\subset\{k-1,\ldots,k-\tau\}$, we obtain the equation \eqref{eqn:inconsist}. %we have $J_i(k)\cap J_j(k)=\emptyset$ if $i\not=j$. Therefore, letting $J(k):=\cup_iJ_i(k)\subset\{k-1,\ldots,k-\tau\}$, .
%Note that \eqref{eqn:inconsist} generalizes \emph{consistent read}, that is,  $\hat{x}^k= x^j$ for some $j\in\{k,k-1,\ldots,k-\tau\}$,  which  corresponds to $J(k)=\{k-1,\ldots,j\}$ in \eqref{eqn:inconsist}. Therefore, the variable history model \eqref{eqn:inconsist} permits both consistent reading and inconsistent reading.
%
%
%Note that under our assumptions,  a set of different agents can   write   to \emph{different} coordinates simultaneously, and they can also read any set of coordinates simultaneously as long as those coordinates are \emph{not} being updated.


\subsubsection{Block coordinate}\label{atomlock}
In the ``block coordinate" case (updating a block of several coordinates each time),  the atomic coordinate update assumption can be met by \emph{either} employing a per-coordinate memory  lock \emph{or} taking the following dual-memory approach: Store \emph{two} copies of each coordinate $x_i\in\cH_i$ in global memory, denoting them as $x_i^{(0)}$ and $x_i^{(1)}$; let a bit $\alpha_i\in\{0,1\}$ point to the active copy; an agent will only read $x_i$ from the active copy $x_i^{(\alpha_i)}$; before an agent updates the components of $x_i$, it  obtains a memory lock to the \emph{inactive} copy  $x_i^{(1-\alpha_i)}$ to prevent other agents from simultaneously updating it; then after it finishes updating $x_i^{(1-\alpha_i)}$, flip the bit $\alpha_i$ so that other agents will begin reading from the updated copy. This approach never blocks any read of $x_i$, yet it  eliminates inconsistency.

% \begin{remark}
% The bounded delay assumption is not an essential assumption. Our analysis can easily extend to unbounded delay case by defining a different metric. Due to page limit, we leave to future work.
% \end{remark}

% \subsubsection{Example of inconsistent read}\label{inconex}
%
%  In this paper, each coordinate corresponds to a Hilbert space $\cH_i$, which can for example be $\RR^{n_i}$ for some $n_i>0$; hence, a coordinate can involve multiple components. Such a coordinate is referred in the literature as ``block coordinate." When a coordinate has multiple components,

\subsection{Straightforward generalization}\label{sec:general}
Our async-parallel coordinate update scheme \eqref{eqn:asyn_update}  can be  generalized to (overlapping) block coordinate updates  after a  change to the step size.\cut{under the same convergence analysis. Such updates can be applied when an orthogonal matrix is a part of the problem data.} Specifically, the scheme \eqref{eqn:asyn_update} can be generalized to
\begin{equation}
\label{eqn:asyn_ortho_update}
\textstyle x^{k+1} = x^k - \frac{\eta_k}{n p_{i_k}} \, (U_{i_k} \circ S)\hat{x}^{k},
\end{equation}
where $U_{i_k}$ is randomly drawn from a set of operators $\{U_1,\ldots,U_n\}$ ($n \leq m$),  $U_i:\cH\to\cH$, following the probability $P(i_k = i) = p_{i}$, $i=1,\ldots, n$ ($p_i>0$, and $\sum_{i=1}^np_i=1$). The operators must satisfy
$\sum_{i=1}^n U_i = I_\cH$ and $\sum_{i=1}^n \|U_ix\|^2 \leq C\|x\|^2$ for some $C>0$.

Let $U_i: x\mapsto(0,\ldots,0,x_i,0,\ldots,0),~i=1,\ldots, m$, which has $C=1$; then~\eqref{eqn:asyn_ortho_update} reduces to~\eqref{eqn:asyn_update}. If  $\cH$ is endowed with a  metric  $M$ such that $\rho_1\|x\|^2\leq \|x\|_M^2\leq \rho_2\|x\|^2$ (e.g., the metric  in the Condat-V\~u primal-dual splitting \cite{condat2013primal,vu2013splitting}), then we have
\begin{align*}
\sum_{i=1}^m \|U_ix\|_M^2 \leq  \rho_2 \sum_{i=1}^m\|U_ix\|^2=\rho_2 \|x\|^2\leq {\frac{\rho_2}{\rho_1}}\|x\|_M^2.
\end{align*}
In general,  multiple coordinates can be updated in \eqref{eqn:asyn_ortho_update}. Consider linear $U_i:x\mapsto (a_{i1}x_1,\cdots,a_{im}x_{m})$, $i=1,\ldots, m$, where $\sum_{i=1}^na_{ij}=1$ for each $j$. Then, for $C:=\max\left\{\sum_{i=1}^n a_{i1}^2,\cdots,\sum_{i=1}^n a_{im}^2\right\}$, we have
\begin{align*}
\sum_{i=1}^n \|U_ix\|^2 = \sum_{i=1}^n\sum_{j=1}^m a_{ij}^2\|x_j\|^2=\sum_{j=1}^m\sum_{i=1}^n a_{ij}^2\|x_j\|^2\leq C\|x\|^2.
\end{align*}


\subsection{Special cases}
If there is only one agent ($p=1$), ARock (Algorithm \ref{alg:asyn_core}) reduces to  randomized coordinate update, which includes the special case of randomized coordinate descent  \cite{nesterov2012rcd} for  convex optimization. Sync-parallel coordinate update is another special case of ARock corresponding to $\hat{x}^k\equiv x^k$. In both cases, there is no delay, i.e., $\tau=0$ and $J(k)=\emptyset$. In addition, the step size $\eta_k$ can be more relaxed. In particular,  if $p_i=\frac{1}{m}$, $\forall i$, then we can let $\eta_k=\eta$, $\forall k$, for any $\eta<1$, or $\eta< 1/\alpha$ when $T$ is $\alpha$-averaged (see Definition~\ref{def:monotone} for the definition of an $\alpha$-averaged operator).

%Because of the orthogonality of $U^k$, the convergence analysis remains the same, and the same convergence result holds.

%\subsubsection{Non-uniform sampling vs. uniform sampling}
%
%\begin{equation}\label{eqn:asyn_ortho_update}
%x^{k+1} = x^k - \eta_k  \beta_{i_k}S_{i_k}\hat{x}^{k},
%\end{equation}
%such that $\beta_{i_k}p_{i_k}=1/m$.


%\input{asyn_t_review}
\begin{subsection}{Related work}\label{sec:rel-wk}

Chazan and Miranker~\cite{chazan1969chaotic} proposed the first async-parallel method  in 1969. The method was designed for solving linear systems. Later, async-parallel methods have been successful applied in many fields, e.g., linear systems~\cite{avron2014revisiting,bethune2014performance,FSS1997asyn-addSch,rosenfeld1969case}, nonlinear problems~\cite{BMR1997asyn-multisplit,baudet1978asynchronous}, differential equations~\cite{aharoni2000parallel,AAI1998implicit,Chau20081126,donzis2014asynchronous}, consensus problems~\cite{LMS1986asynchronous,leifang_information_2005}, and optimization~\cite{hsieh2015passcode,liu2014asynchronous,liu2013asynchronous,tai2002convergence,zhang2014asynchronous}. We review the theory for  async-parallel fixed-point iteration and its applications.



\cut{The asynchronous parallel computing method (also called \emph{chaotic relaxation}) appears to be first used by Rosenfeld~\cite{rosenfeld1969case} to solve linear equations arising in electrical network problem. Numerical simulations there show that using more processors can save total running time and sometimes achieve near-linear speedup. Chazan and Miranker~\cite{chazan1969chaotic} systematically analyzed a general asynchronous iterative method for solving linear systems. Assuming boundedness of the delay, it guarantees convergence to a solution if the spectral radius of the (componentwise) absolute of the iterating matrix is strictly less than \emph{one}. Since then, the theory and application has been studied by many authors.}

\textbf{General fixed point problems.} \emph{Totally async-parallel}\footnote{``Totally asynchronous" means  no upper bound on the delays; however, other conditions are required, for example: each coordinate must be updated infinitely many times. By default, ``asynchronous" in this paper assumes a finite maximum delay.} iterative methods for a  fixed-point problem go back as early as to  Baudet~\cite{baudet1978asynchronous}, where the operator was assumed to be \emph{P-contraction}.\footnote{An operator $T: \mathbb{R}^n \rightarrow \mathbb{R}^n$ is P-contraction if $|T(x) - T(y)| \leq P | x - y|$, component-wise, where $|x|$ denotes the vector with components $|x_i|, ~ i=1, ..., n$, and $P\in\RR^{n\times n}$ is a nonnegative matrix with a spectral radius strictly less than 1.
} Later, Bertsekas~\cite{bertsekas1983distributed} generalized the P-contraction assumption and showed convergence. \cut{Papers~\cite{Baz200591,el1998flexible} provide a convergence result for async-parallel iterations with simultaneous reading and writing of the same set of components. Strikwerda considered the unbounded but stochastic delays~\cite{Strikwerda2002125}.} Frommer and Szyld~\cite{frommer2000asynchronous} reviewed the theory and applications of totally async-parallel iterations prior to 2000. This review summarized convergence results under the conditions in~\cite{bertsekas1983distributed}.
%\commzp{ (Based on my discussion with Ming, this is not right.) If the conditions are satisfied, the operator is contractive, thus nonexpansive, under some norm for $\cH$~\cite[Exercise after Lemma 5.6.10]{Horn_matrix_2012}, then we can apply ARock to find the unique fixed point of the operator with a smaller stepsize and bounded delays (see Section~\ref{sec:general})}.
However, ARock can be applied to solve many more problems since our nonexpansive assumption, though not strictly weaker than P-contraction, is more pervasive.  %including those that  the conditions required by previous totally async-parallel methods.
As opposed to  totally asynchronous methods, Tseng, Bertsekas, and Tsitsiklis~\cite{bertsekas1989parallel,TB1990partially} assumed \emph{quasi-nonexpansiveness}\footnote{An operator $T: \cH \rightarrow \cH$ is quasi-nonexpansive if $\|Tx - x^*\| \leq \|x - x^*\|$, $\forall x \in \cH$, $x^*\in \Fix T$.} and proposed an {async-parallel} %\footnote{``Partially asynchronous" means that each agent performs an update at least once during any time interval of length $B$, and the information used by any agent is outdated by at most $B$ time units.}
method, converging under an additional assumption, which is   difficult to justify in general but can be established   for problems such as linear systems and strictly convex network flow problems \cite{bertsekas1989parallel,TB1990partially}. %, the additional assumption can be satisfied under certain conditions.   %In the rest of the paper, unless specified, ``asynchronous" refers ``partially asynchronous."

{The above  works assign coordinates  in a deterministic manner. Different from them, ARock is stochastic, works for nonexpansive operators, and is more applicable.}


\textbf{Linear, nonlinear, and differential equations.}
The first async-parallel method for solving linear equations was introduced by Chazan and Miranker in~\cite{chazan1969chaotic}. They proved that on solving linear systems, P-contraction was necessary and sufficient for convergence. The performance of the algorithm was studied by Iain et al.~\cite{bethune2014performance,rosenfeld1969case} on different High Performance Computing (HPC) architectures. %Frommer, Schwandt, and Szyld~\cite{FSS1997asyn-addSch} proposed a class of Schwarz methods.
Recently, Avron et al.~\cite{avron2014revisiting} revisited the async-parallel coordinate update and showed its linear convergence for solving positive-definite linear systems. %with a symmetric positive-definite matrix. %Linear convergence rate was derived.
Tarazi and Nabih~\cite{el1982some} extended the poineering work \cite{chazan1969chaotic} to solving nonlinear equations, and the async-parallel methods have also been applied \cut{Applications of async-parallel methods }for solving differential equations, \cut{can be found}e.g., in~\cite{aharoni2000parallel,AAI1998implicit,Chau20081126,donzis2014asynchronous}. Except for~\cite{avron2014revisiting}, all these methods are totally async-parallel with the P-contraction condition or its variants. \cut{except~\cite{avron2014revisiting} employs random selection to solve a special linear system of equations. }On solving a positive-definite linear system, \cite{avron2014revisiting} made assumptions similar to ours, and it obtained better linear convergence rate on that special problem. %However, our work applies to many more problems.

\textbf{Optimization.}
The first async-parallel coordinate update gradient-projection method was due to Bertsekas and Tsitsiklis~\cite{bertsekas1989parallel}. The method solves constrained optimization problems with a smooth objective and simple constraints. It was shown that the objective gradient  sequence converges to zero. Tseng~\cite{tseng1991rate-asyn} further analyzed the convergence rate and obtained local linear convergence based on the assumptions of isocost surface separation and a local Lipschitz error bound. Recently, Liu et al.~\cite{liu2013asynchronous} developed an async-parallel stochastic coordinate descent algorithm for minimizing convex smooth functions. Later, Liu and Wright~\cite{liu2014asynchronous} suggested an async-parallel stochastic proximal coordinate descent algorithm for minimizing convex composite objective functions. They established the convergence of the expected objective-error sequence for convex functions. Hsieh et al.~\cite{hsieh2015passcode} proposed an async-parallel dual coordinate descent method  for solving $\ell_2$ regularized empirical risk minimization problems. Other async-parallel approaches include asynchronous ADMM~\cite{hong2014distributed,wei2013on,zhang2014asynchronous,iutzeler2013asynchronous}. Among them,~\cite{wei2013on,iutzeler2013asynchronous} use an asynchronous clock, %assume that the nodes do not read and write the same data simultaneously,
and~\cite{hong2014distributed,zhang2014asynchronous} use a central node to update the dual variable; they do not deal with delay or inconsistency.  Async-parallel stochastic gradient descent methods have also been considered~in~\cite{nedic2001distributed,recht2011hogwild}.

Our framework differs from the recent surge of the aforementioned sync-parallel and async-parallel coordinate descent algorithms (e.g.,~\cite{peng2013parallel,kyrola2011parallel,liu2013asynchronous,liu2014asynchronous,hsieh2015passcode,richtarik2015parallel}). While they  apply to convex function minimization, ARock covers more cases (such as ADMM, primal-dual, and decentralized methods) and also provides sequence convergence. %aims at finding a fixed point to a nonexpansive operator, and its iterations decrease the expected distance to the fixed point set.
In Section \ref{sec:application}, we will show that some of the existing async-parallel coordinate descent algorithms are special cases of ARock, through relating their optimality conditions to nonexpansive operators. Another difference is that the convergence of ARock only requires a  nonexpansive operator with  a fixed point, whereas properties such as strong convexity, bounded feasible set, and bounded sequence, which  are  seen in some of the recent literature for async-parallel convex minimization, are unnecessary.


\textbf{Others.} Besides solving equations and optimization problems, there are also applications of  async-parallel algorithms to optimal control problems~\cite{LMS1986asynchronous}, network flow problems~\cite{ESMG1996asyn-flex}, and consensus problems of multi-agent systems~\cite{leifang_information_2005}.


\cut{The pioneering work for asynchronous iteration is due to Chazan and Miranker~\cite{chazan1969chaotic}, who proposed an asynchronous method for solving linear equations with convergence guarantee. The theory and application has since been studied by many authors. Noteworthy, Bertsekas and Tsitsiklis~\cite{bertsekas1989parallel} introduced an asynchronous parallel framework for general fixed point problems. More recently, Frommer and Szyld~\cite{frommer2000asynchronous} gave a review for the theory and application for asynchronous iterations. Baudet~\cite{baudet1978asynchronous} analyzed asynchronous iterations of contracting operators, so their results are only applicable to a very restive type of equations. Bertsekas and Tsitsiklis~\cite{bertsekas1989parallel} proved the convergence for asynchronous fixed point iterations of nonexpansive maps. However the convergence analysis is based on very restrictive conditions and in a hard to interpret manner.}
\end{subsection}


\begin{subsection}{Contributions}
%The techniques developed in this paper are different from those used in the existing  convergence analysis of async-parallel algorithms, which assumes contraction-type properties for the underlying operator. It is also different from many of the async-parallel optimization algorithms, which   monotonically reduce the expected objectives, since our algorithms are driven by a fixed-point iteration in which  the distance between the current point and the fixed-point set is stochastically monotone.
Our contributions and techniques are summarized below:
\begin{itemize}
\item ARock is the first async-parallel coordinate update framework for finding a fixed point to a nonexpansive operator. %, based on the Krasnosel'ski\u i--Mann (KM) iteration.
\item %We consider  a nonexpansive operator with a fixed point.
By introducing a new metric and establishing stochastic Fej\'er monotonicity, we show  that, with probability one,   ARock   converges to a point in the solution set; %strong convergence is obtained under slightly stronger conditions;
linear convergence is obtained for~\emph{quasi-strongly monotone} operators.
\item Based on ARock, we introduce an async-parallel algorithm for linear  systems, async-parallel ADMM algorithms for  distributed or decentralized  computing problems, as well as async-parallel operator-splitting algorithms for  nonsmooth minimization problems. Some problems are treated in they async-parallel fashion for the first time in history. The developed algorithms are \emph{not} straightforward modifications to their serial versions because their underlying nonexpansive operators must be identified before applying ARock.
\end{itemize}
\end{subsection}

%\input{asyn_t_notation}
\begin{subsection}{Notation, definitions, background of monotone operators}
Throughout this paper, $\cH$ denotes a separable Hilbert space equipped with the inner product $\langle \cdot, \cdot \rangle$ and norm $\| \cdot \|$, and $(\Omega, \cF, P)$ denotes the underlying probability space, where  $\Omega$, $\cF$, and $P$ are the sample space,  $\sigma$-algebra, and  probability measure, respectively. The map $x: (\Omega, \cF) \rightarrow (\cH, \cB)$, where  $\cB$ is the Borel $\sigma$-algebra, is an $\cH$-valued random variable. Let  $(x^k)_{k\ge 0}$ denote \emph{either} a sequence of deterministic points in $\cH$ \emph{or} a sequence of $\cH$-valued random variables, which will be clear from the context, and let $x_i\in\cH_i$ denote the $i$th coordinate of $x$. In addition, we let $\cX^k := \sigma (x^0, \hat{x}^1, x^1, ..., \hat{x}^k, x^k)$ denote the smallest $\sigma$-algebra generated by $x^0, \hat{x}^1, x^1, ..., \hat{x}^k, x^k$.  ``Almost surely'' is abbreviated as ``a.s.'', and the $n$ product space of $\cH$ is denoted by $\cH^n$.  We use $\to$ and $\rightharpoonup$ for strong convergence and weak convergence, respectively.



We define $\Fix T:= \{x\in \cH ~|~ Tx = x\}$ as the set of fixed points of operator $T$, and, in the product space,  we let $\vX^{*} :=\{(x^*, x^*, ..., x^*) ~|~ x^* \in \Fix T\}\subseteq \cH^{\tau+1}$.



%$\cB$ for Borel $\sigma$-algebra. The underlying probability space is denoted as $(\Omega, \cF, P)$, where $\Omega$, $\cF$ and $P$ represent the sample space, the $\sigma$-algebra, and the probability measure, respectively. The map $x: (\Omega, \cF) \rightarrow (\cH, \cB)$ is an $\cH$-valued random variable.


\begin{definition}\label{def:quasi_lip}
An operator $T: \cH \rightarrow \cH$ is \emph{$c$-Lipschitz}, where $c\ge 0$, if it satisfies $\|Tx - T y\| \leq c\|x - y\|$, $\forall x, y \in \cH$. In particular,
 $T$ is  \emph{nonexpansive} if $c\le 1$, and  \emph{contractive} if $c<1$.
\cut{
In addition, when $\Fix T\not=\emptyset$, $T$ is \emph{quasi $\bar{c}$-Lipschitz},  $\bar{c}\ge 0$, if it satisfies $\|Tx - y\| \leq \bar{c}\|x - y\|$, $\forall x \in \cH$ and $\forall y\in \Fix T$. $T$ is  \emph{quasi-nonexpansive} if $\bar{c}\le 1$, and \emph{quasi-contractive} if $\bar{c}<1$.}
\end{definition}

\begin{definition}\label{def:monotone} Consider an operator $T: \cH \rightarrow \cH$.
\begin{itemize}
\item $T$ is \emph{$\alpha$-averaged} with $\alpha \in (0, 1)$, if there is a nonexpansive operator $R:\cH \rightarrow \cH$ such that $T = (1 - \alpha) I_\cH + \alpha R$, where $I_\cH: \cH\rightarrow \cH$ is the identity operator.
\item $T$ is \emph{$\beta$-cocoercive} with $\beta>0$, if \cut{it satisfies }$\dotp{x - y, Tx - Ty} \geq \beta \|Tx - Ty\|^2,~\forall x, y \in \cH.$
\item $T$ is \emph{$\mu$-strongly monotone}, where $\mu>0$, if it satisfies $\dotp{x - y, Tx - Ty} \geq \mu\|x - y\|^2,~\forall x, y \in \cH.$ When the inequality holds for $\mu=0$, $T$ is \emph{monotone}.
\item $T$ is \emph{quasi-$\mu$-strongly monotone}, where $\mu>0$, if it satisfies $\dotp{x - y, Tx } \geq \mu\|x - y\|^2,~\forall x \in \cH,y\in\zer T:=\{y\in \cH\mid Ty=0\}$. When the inequality holds for $\mu=0$, $T$ is \emph{quasi-monotone}.
\item $T$ is \emph{demicompact~\cite{petryshyn1966construction} at $x\in\cH$} if for every bounded sequence $(x^k)_{k\geq0}$ in $\cH$ such that $Tx^k-x^k\to x$, there exists a strongly convergent subsequence.
\end{itemize}
\end{definition}
Averaged operators are nonexpansive. By the Cauchy-Schwarz inequality, a  {$\beta$-cocoercive} operator is $\frac{1}{\beta}$-Lipschitz; the converse is generally untrue, but true for the gradients of convex differentiable functions. Examples are given in the next section.\cut{(However, if $f$ is a convex differentiable function, by the Baillon-Haddad theorem, $\nabla f$ is  {$\beta$-cocoercive} if and only if it is $\frac{1}{\beta}$-Lipschitz.) It is worth noting that we can extend the definitions of (quasi) (strongly) monotone operators to \emph{set-valued operators} by replacing $Tx$ and $Ty$ in their definitions with any entries $p\in Tx$ and $q\in Ty$, respectively. When $T$ is monotone (possibly set-valued) and $\lambda >0$, its resolvent $(I+\lambda T)^{-1}$ is single-valued and $\frac{1}{2}$-averaged and thus also nonexpansive. If $T$ is strongly monotone, then$ (I+\lambda T)^{-1}$ is contractive.}

\cut{\begin{definition} Let $T:\cH\to\cH$ be a \emph{nonexpansive} or \emph{averaged} operator such that $\Fix T\not=\emptyset$. Let $\lambda\in(0,1]$. The Krasnosel'ski\u i--Mann (KM) iteration generates the sequence $(x^k)_{k\ge 0}$ by
$$ x^{k+1} = x^k + \lambda (Tx^k-x^k).$$
\end{definition}
The sequence $(x^k)_{k\ge 0}$ converges weakly to a point in $\Fix T$, and $(Tx^k-x^k)_{k\ge 0}$ converges strongly to 0.

Let us illustrate the above concepts using a proper closed convex function $f:\cH\to\RR\cup\{+\infty\}$. The following results are well known (cf.  textbook \cite{bauschke2011convex}). The subdifferential $\partial f$ is set-valued and~maximally monotone, and thus  $(I+\lambda\partial f)^{-1}$ is single-valued and  $\frac{1}{2}$-averaged. For any $x\in\cH$, $(I+\lambda\partial f)^{-1}x=\argmin_y f(y)+\frac{1}{2\lambda}\|y-x\|^2$. Therefore, the proximal-point algorithm is a  KM iteration.

If $f$ is differentiable and $\nabla f$ is $L$-Lipschitz, then $\nabla f$ is $\frac{1}{L}$-cocoercive and, by \cite[Prop. 4.33]{bauschke2011convex}, the gradient descent operator $(I-\gamma\nabla f)$ is $\frac{\gamma L}{2}$-averaged, where $\gamma$ satisfies $0<\gamma <\frac{2}{L}$ so that $0<\frac{\gamma L}{2}<1$. Therefore, the gradient descent iteration  is a KM iteration.

In addition, other special cases of the KM iteration are the forward-backward splitting (FBS) iteration \cite{passty1979ergodic}, Douglas-Rachford splitting (DRS) iteration \cite{lions1979splitting}, relaxed Peaceman-Rachford splitting (PRS) \cite{lions1979splitting}, and a three-operator splitting iteration \cite{davis2015three}, which involve multiple nonexpansive or averaged operators. They have numerous applications in many areas of computational mathematics.
}
\cut{If $f$ is $\mu$-strongly convex, then $\partial f$ is $\mu$-strongly monotone and, by \cite[Prop. 23.11]{bauschke2011convex}, $(I+\lambda\partial f)^{-1}$ is $\frac{1}{1+\lambda\mu}$-contractive. In this case, the proximal-point algorithm is a  Picard iteration (of a contractive operator). }

\end{subsection}

%\input{asyn_t_apps}
\begin{section}{Applications}\label{sec:application}
%Our async-parallel framework ARock generates simple numerical schemes %(oftentimes with closed form expressions) for a large number of problems that appear in scientific computing, signal processing, machine learning, sensor network, and other data science areas.
In this section, we provide some applications that are special cases of the fixed-point problem~\eqref{eqn:fix_point_set}. For each application, we identify its  nonexpansive operator $T$ (or the corresponding operator $S$) and implement the conditions in Theorem~\ref{thm:simple_convergence}. For simplicity, we use the uniform distribution, $p_1=\cdots=p_m=1/m$, and apply the simpler update~\eqref{eqn:simple_asyn_update} instead of~\eqref{eqn:asyn_update}. %All the algorithms in this section can take non-uniform distributions.

\begin{subsection}{Solving linear equations}\label{sec:linearequations}
%The problem is
Consider the linear system $Ax = b,$
%\begin{equation}\label{eqn:linear_eqn}
%Ax = b,
%\end{equation}
where $A \in \mathbb{R}^{m\times m}$ is a nonsingular matrix with nonzero diagonal entries. Let $A = D + R$, where $D$ and $R$ are the diagonal and off-diagonal parts of $A$, respectively. Let $M:=-D^{-1}R$ and $T(x):=M x + D^{-1} b$. Then the system $Ax = b$ is equivalent to the fixed-point problem $x = D^{-1} (b - R x)= :T(x)$, where $T$ is nonexpansive if the spectral norm $\|M\|_2$ satisfies $\|M\|_2\leq 1$.
\cut{Define the relaxed operator as
$$T_{\alpha}(x) = \big((1 - \alpha) I + \alpha \, T \big) x,$$
where $\alpha \in (0, 1)$. Note that $T_{\alpha}$ is $\alpha$-averaged, $x^{k+1} = T_{\alpha} (x^k)$ gives the relaxed Jacobi iteration.}The iteration $x^{k+1}= T(x^k)$ is widely known as the Jacobi algorithm. Let $S = I -T$. Each update  $S_{i_k}\hat{x}^k$ involves multiplying just the $i_k$th  row of $M$ to $x$ and adding the $i_k$th entry of $D^{-1} b$, so we arrive at the following algorithm.

\begin{algorithm}[H]\label{alg:asyn_linear_eqn}
%\SetAlgoLined
\SetKwInOut{Input}{Input}\SetKwInOut{Output}{output}
\Input{$x^0 \in \mathbb{R}^n$,  $K>0$.}
set the global iteration counter $k=0$\;
 \While{$k < K$, every agent asynchronously and continuously}{
  select $i_k\in\{1,\ldots, m\}$ uniformly at random\;%from $\{1, ..., m\}$ under a discrete positive distribution\;
  %$\hat{x}^k\gets$ read the variable $x$\;
  subtract $\frac{\eta_k}{a_{i_k i_k}} (\sum_{j} a_{i_k j} \, \hat x_j^k - b_{i_k})$ from  the component $x_{i_k}$ of the variable $x$\;

  %add $\Delta x$ to the component $x_{i_k}$ of the variable $x$\;
  update the global counter $k \leftarrow k+1$\;
 }
 \caption{ARock for linear equations}
\end{algorithm}

\begin{proposition}\label{prop:lip_to_sm}
{\cite[Example~22.5]{bauschke2011convex}} Suppose that $T$ is $c$-Lipschitz continuous with $c \in [0, 1)$. Then, $I - T$ is $(1 - c)$-strongly monotone.
\end{proposition}

Suppose $\|M\|_2<1$. Since $T$ is $\|M\|_2$-Lipschitz continuous, by Proposition~\ref{prop:lip_to_sm}, $S$ is $(1-\|M\|_2)$-strongly monotone. By Theorem~\ref{thm:strongly_monotone}, Algorithm~\ref{alg:asyn_linear_eqn} converges linearly. % when properly step sizes $\eta_k$ are chosen.

\cut{\begin{lemma}\label{lem:arock-lin}
Assume that $A$ is symmetric positive definite. Then, there exists a constant $\mu>0$ such that
\begin{equation}\label{eq:S-lin} \langle S(x)-S(y), x - y \rangle \ge \mu\|x-y\|^2,\quad\forall x,
\end{equation}
where $S(x)=x-(M x + D^{-1} b)$ is  defined in this subsection, i.e., $S$ is $\mu$-strongly monotone.
\end{lemma}
\begin{proof}
Note
\begin{align*}
S(x)=x-(Mx+D^{-1}b)=x-\big(-D^{-1} (A-D)x+D^{-1}b\big)=D^{-1}(Ax-b).
\end{align*}
Hence,~\eqref{eq:S-lin} reduces to
$$ \langle D^{-1}A (x - y), x - y \rangle \geq \mu \|x - y\|^2$$
which holds for $\mu=  \lambda_{\min} (D^{-1} A) $ with $\lambda_{\min} (D^{-1}A) $ being the smallest eigenvalue of $D^{-1} A$. This completes the proof.\hfill
\end{proof}}

\cut{\begin{align}\label{eqn:linear_eqn_solver}
x^{k+1} &= x^k - \eta_k \, U_k( S({\hat x}^{k}) ) \nonumber \\
             &= x^{k} - \eta_k \, U_k( \alpha (I - T)\hat x^{k} ) \\
             &=x^{k} - \eta_k \, \alpha \, U_k (\hat x^{k} - M \hat x^{k} - D^{-1} b) \nonumber.
\end{align}
Since $T_{\alpha}$ is $\alpha$-averaged, Theorem ~\ref{thm:convergence} guarantees the convergence of iteration~\eqref{eqn:linear_eqn_solver}.} \cut{When $A$ is symmetric, we know that $\rho(|A|) \leq \rho(A) = \|A\|_2$, so our convergence condition is weaker compared to the existing result by \cite{chazan1969chaotic} which requires $\rho(|M|) < 1$. In addition, when the smallest eigenvalue of $I - M$ is strictly positive, Theorem~\ref{thm:strongly_monotone} implies the linear convergence of our algorithm.} \cut{, then if $\lambda_{\min}>0$, then we have
$$\big\langle S(x) - S(y), x - y \big\rangle = \alpha \big\langle (I - M) (x-y), x - y \big\rangle \geq \alpha \,\lambda_{\min} \|x - y\|^2,$$
i.e., $S$ is strongly monotone with modulus $\alpha \, \lambda_{\min}$. Hence, if  $\lambda_{\min} > 0$, then the sequence $(x^k)$ generated by~\eqref{eqn:linear_eqn_solver} converges to the linear equation solution with linear rate.}
\end{subsection}


\cut {\begin{subsection}{Solving nonlinear equations}
Consider the following nonlinear system
\begin{equation}
\min_{x\in {\cH}} F(x) = 0,
\end{equation}
\end{subsection}

\begin{subsection}{Solving ordinary differential equations (ODE)}
Consider the following nonlinear system
\begin{equation}
\min_{x\in {\cH}} F(x) = 0,
\end{equation}

\end{subsection}
}

\begin{subsection}{Minimize convex smooth function}\label{subsec:smooth}
Consider the optimization problem
\begin{equation}\label{eq:probmin}
\Min_{x\in {\cH}} f(x),
\end{equation}
where $f$ is a closed proper convex differentiable function and $\nabla f$ is $L$-Lipschitz continuous, $L>0$. Let $S := \frac{2}{L} \, \nabla f$. \cut{ By the Baillon-Haddad theorem, $S$ is $\frac{1}{\alpha L}$-cocoercive\footnote{An operator $S:\cH\to\cH$ is \emph{$\beta$-cocoercive},  where $\beta>0$, if it satisfies $\dotp{x - y, Tx - Ty} \geq \beta \|Tx - Ty\|^2,~\forall x, y \in \cH.$}, and by Lemma~\ref{lemma:a-avg}, $T = I - S$ is $\frac{\alpha L}{2}$-averaged, where $\frac{\alpha L}{2}<1$.} As $f$ is  convex and differentiable, $x$ is a minimizer of $f$ if and only if $x$ is a zero of $ S$. Note that $S$ is $\frac{1}{2}$-cocoercive. By Lemma \ref{lemma:a-avg}, $T \equiv I-S$ is nonexpansive. Applying ARock, we have the following iteration:
\begin{equation}\label{eqn:cvx_itr}
x^{k+1} = x^k - \eta_k  S_{i_k} \hat{x}^{k},
\end{equation}
where $S_{i_k} x = \frac{2}{L}(0, ..., 0, \nabla_{i_k}f(x), 0, ..., 0)^T$. Note that $\nabla f$ needs a structure that makes it cheap to compute $\nabla_{i_k} f(\hat{x}^{k})$. % than to evaluate the entire vector $\nabla f(\hat x^k)$. %Otherwise, one would rather compute $\nabla f(\hat x^k)$ and perform the standard gradient update.
Let us give two such examples%where $\nabla_i f(x)$ is much easier to compute than  $\nabla f(x)$
: (i) quadratic programming: $f(x)= \frac{1}{2} x^T A x - b^T x$, where $\nabla f(x) = Ax - b$ and $\nabla_{i_k} f(\hat{x}^{k})$ only depends  on a part of $A$ and $b$; (ii) sum of sparsely supported functions: $f=\sum_{j=1}^N  f_j$ and $\nabla f = \sum_{j=1}^N \nabla f_j$, where each  $f_j$ depends on just a few variables.

Theorem~\ref{thm:convergence} below guarantees the convergence of $(x^k)_{k\geq0}$ if $\eta_k \in [\eta_{\min}, \frac{1}{ 2\tau/\sqrt{m} +1})$. In addition, If $f(x)$ is \emph{restricted strongly convex}, namely, for any $x \in \cH$ and $x^* \in X^*$, where $X^*$ is the  solution set to~\eqref{eq:probmin}, we have
$\langle x- x^*, \nabla f(x) \rangle \geq \mu \|x - x^*\|^2$
for some $\mu>0$, then $S$ is quasi-strongly monotone with modulus $\mu$. According to Theorem~\ref{thm:strongly_monotone}, iteration~\eqref{eqn:cvx_itr} converges at a linear rate if the step size meets the condition therein.

Our convergence and rates are given in term of  the distance to the solution set $X^*$. In comparison, the results in the work \cite{liu2013asynchronous} are given in terms of objective error under the assumption of a uniformly bounded $(x^k)_{k\ge 0}$. In addition, their step size decays like $O(\frac{1}{\tau \rho^\tau})$ for some $\rho > 1$ depending on $\tau$, and our $O(\frac{1}{\tau})$ is better. Under similar assumptions, Bertsekas and Tsitsiklis~\cite[Section 7.5]{bertsekas1989parallel} also describes an algorithm for~\eqref{eq:probmin} and proves  only subsequence convergence~\cite[Proposition 5.3]{bertsekas1989parallel}  in $\RR^n$. % is established.
\end{subsection}

\begin{subsection}{Decentralized consensus optimization}\label{sec:dgd}
Consider that $m$ agents in a connected  network   solve the  consensus  problem of minimizing $\sum_{i=1}^m f_i (x)$, where $x\in\RR^d$ is the shared variable and the  convex differentiable function $f_i$ is held privately by agent $i$. We assume that $\nabla f_i$ is $L_i$-Lipschitz continuous for all $i$. A decentralized gradient descent algorithm \cite{nedic2009distributed} can be developed based on the equivalent formulation
\begin{equation}\label{eqn:dect_problem_org}
\Min_{x_1, ..., x_m \in \mathbb{R}^d}  ~ \textstyle f(\vx) := \sum_{i=1}^m f_i (x_i),\quad \St~  W\vx = \vx,
\end{equation}
%\begin{align}\label{eqn:dect_problem_org}
%\Min_{x_1, ..., x_m \in \mathbb{R}^d}  ~& \textstyle f(\vx) := \sum_{i=1}^m f_i (x_i) \nonumber \\
%                \St~ & W\vx = \vx,
%\end{align}
where $\vx = (x_1, ..., x_m)^T\in \RR^{m\times d}$ and $W \in \mathbb{R}^{m\times m}$ is the so-called mixing matrix satisfying: $W\vx = \vx$ if and only if $x_1=\cdots=x_m$. For  $i\not=j$, if $w_{i,j} \neq 0$, then agent $i$ can communicate with agent $j$; otherwise they cannot. We assume that $W$ is  symmetric and doubly stochastic. Then, the decentralized consensus algorithm~\cite{nedic2009distributed} can be expressed as $\vx^{k+1}=W\vx^{k}-\gamma\nabla f(\vx^k)=\vx^k-\gamma(\nabla f(\vx^k)+\frac{1}{\gamma}(I-W)\vx^k)$, where $\nabla f(\vx)\in \RR^{m\times d}$ is a matrix with its $i$th row equal to $(\nabla f_i(x_i))^T$; see~\cite{yuan2013convergence}. The computation of $W\vx^k$ involves communication between agents, and $\nabla f_i(x_i)$ is independently computed by each agent $i$. The  iteration is equivalent to  the gradient descent iteration applied to $\min_{\vx}  \sum_{i=1}^m f_i (x_i) + \frac{1}{2\gamma} \vx^T(I - W)\vx$. To apply our algorithm, we let $S := \frac{2}{L} \nabla F= \frac{2}{L} (\nabla f+\frac{1}{\gamma}(I-W))$ with $L = \max_{i} L_i + (1 - \lambda_{\min}(W))/\gamma$, where $\lambda_{\min} (A)$ is the smallest eigenvalue of $W$. Computing $S_{i}\hat{\vx}^k$ reduces to computing  $\nabla f_i(\hat{x}_i^k)$ and the $i$th entry of $W\hat{\vx}^k$ or $\sum_{j}w_{i,j}\hat{x}^k_j$, which involves only $\hat{x}_i^k$ and $\hat x_j^k$ from the neighbors of agent $i$. Note that since each agent $i$ can store its own $x_i$ locally, we have $\hat{x}^k_i\equiv x^k_i$.

If the agents are $p$ independent Poisson processes and that each agent $i$ has activation rate $\lambda_i$, then the probability that agent $i$ activates before other agents is equal to $\frac{\lambda_i}{\sum_{i=1}^p \lambda_i}$~\cite{larson1981urban} and therefore our random sample scheme holds and ARock applies naturally. The algorithm is summarized as follows:

\begin{algorithm}[H]\label{alg:asyn_dgd}
\SetKwInOut{Input}{Input}\SetKwInOut{Output}{output}
\Input{Each agent $i$ sets $x_i^0\in\RR^d$, $K>0$.}
 \While{$k < K$}{
  when an agent $i$ is activated, $x^{k+1}_{i} = x^k_{i} -  \frac{\eta_k}{L} \, (\nabla f_{i}( x_{i}^k) + \frac{1}{\gamma} (x^k_{i}  -  \sum_j w_{i,j} \, \hat x_j^k) )$\;
  increase the global counter $k \leftarrow k+1$\;
 }
 \caption{ARock for decentralized optimization~\eqref{eqn:dect_problem_org}}
\end{algorithm}
%Compared to the synchronized version of the DGD algorithm \cite{yuan2013convergence}, Algorithm~\ref{alg:asyn_dgd} doesn't require a global lock, and each agent can compute and communicate at its own pace and collects the neighboring information without waiting for the neighboring agents finishing the computing.

\end{subsection}



\begin{subsection}{Minimize smooth $+$ nonsmooth functions}\label{subsec:nonsmooth}
Consider the problem
\begin{equation}\label{eqn:prob_fbs}
\Min_{x\in \cH} f(x) + g(x),
\end{equation}
where $f$ is closed proper convex  and $g$ is convex and $L$-Lipschitz differentiable with $L>0$. Problems in the form  of~\eqref{eqn:prob_fbs} arise in statistical regression, machine learning, and signal processing and include well-known problems such as the  support vector machine, regularized least-squares, and regularized logistic regression. For any $x \in \cH$ and scalar $\gamma\in (0, \frac{2}{L})$,   define the proximal operator $\prox_f: \cH \rightarrow \cH$ and the reflective-proximal operator $\refl_f: \cH \rightarrow \cH$ as
\begin{equation}\label{eqn:prox}
\prox_{\gamma f} (x) := \argmin_{y \in \cH} f(y) + \frac{1}{2\gamma} \|y - x\|^2 \quad \text{~and~} \quad\refl_{\gamma f} := 2\prox_{\gamma f} - I_{\cH},
\end{equation}
respectively, and define the following forward-backward operator $\TFBS := \prox_{\gamma f}  \circ (I - \gamma \nabla g).$
Because $\prox_{\gamma f}$ is $\frac{1}{2}$-averaged and $(I - \gamma \nabla g)$ is $\frac{\gamma L}{2}$-averaged,  $\TFBS$ is $\alpha$-averaged for
%$$\alpha_{FBS} := \frac{\max\{1, \gamma L\} }{ \max\{1/2, \gamma L / 2\} + 1},$$
$\alpha \in [\frac{2}{3}, 1)$ \cite[Propositions 4.32 and 4.33]{bauschke2011convex}.
 Define
$S := I - \TFBS = I - \prox_{\gamma f}  \circ (I - \gamma \nabla g)$.   When we apply Algorithm~\ref{alg:asyn_core} to $T=\TFBS$ to solve ~\eqref{eqn:prob_fbs}, and assume $f$ is separable in  all coordinates, that is, $f(x) = \sum_{i=1}^mf_i(x_i)$,  the update for the $i_k$th selected coordinate is
\begin{equation}\label{asyncFBS}
x_{i_k}^{k+1} = x_{i_k}^k - \eta_k  \,   \big( \hat{x}_{i_k}^k - \prox_{\gamma f_{i_k}}( \hat{x}_{i_k}^k - \gamma \nabla_{i_k} g (\hat{x}^{k}))\big),
\end{equation}
%The other coordinates remain the same.
Examples of separable functions include $\ell_1$ norm, $\ell_2$ norm square, the Huber function, and the indicator function of box constraints, i.e., ${\{x| a_i \leq x_i \leq b_i,~\forall i \}}$.  They all have simple $\prox$ maps.
If $\eta_k \in [\eta_{\min}, \frac{1}{2\tau/\sqrt{m} +1})$, then the convergence   is guaranteed by Theorem ~\ref{thm:convergence}. To show  linear convergence, we need to assume that $g(x)$ is strongly convex. Then, Proposition~\ref{prop:prox} below shows that $\prox_{\gamma f}  \circ (I - \gamma \nabla g)$ is a quasi-contractive operator, and  by Proposition~\ref{prop:lip_to_sm}, operator $I - \prox_{\gamma f}  \circ (I - \gamma \nabla g)$ is quasi-strongly monotone. Finally, linear convergence and its rate follow from Theorem~\ref{thm:strongly_monotone}.
\begin{proposition}\label{prop:prox}
Assume that $f$ is a closed proper convex function, and $g$ is $L$-Lipschitz differentiable and strongly convex with modulus $\mu > 0$. Let $\gamma \in (0, \frac{2}{L})$. Then, both $I - \gamma \nabla g$ and $\prox_{\gamma f}  \circ (I - \gamma \nabla g)$ are quasi-contractive operators.
\end{proposition}
\begin{proof}
We first show that $I - \gamma \nabla g$ is a quasi-contractive operator. Note
\begin{align*}
~ &\|(x - \gamma \nabla g (x)) - (x^* - \gamma \nabla g(x^*))\|^2 \\
= & \|x - x^*\|^2 - 2\gamma \langle x - x^*, \nabla g(x) - \nabla g(x^*) \rangle + \gamma^2 \|\nabla g(x) - \nabla g(x^*)\|^2 \\
\leq & \|x - x^*\|^2 - \gamma(2 - \gamma L)  \langle x - x^*, \nabla g(x) - \nabla g(x^*) \rangle \\
\leq & (1 - 2\gamma \mu + {\mu \gamma^2}L) \|x - x^*\|^2,
\end{align*}
where the first inequality follows from the  Baillon-Haddad theorem\footnote{Let $g$ be a convex differentiable function. Then, $\nabla g$ is $L$-Lipschitz if and only if it is $\frac{1}{L}$-cocoercive.} and the second one from the strong convexity of $g$.
Hence, \cut{we have
\begin{equation*}
\|(x - \gamma \nabla g (x)) - (x^* - \gamma \nabla g(x^*))\| \leq \sqrt{1 - 2\gamma \mu + {\mu \gamma^2}L} \, \|x - x^*\|,
\end{equation*}
which means that }$I - \gamma \nabla g$ is quasi-contractive if $0<\gamma<2/L$. Since $f$ is convex, $\prox_{\gamma f}$ is firmly nonexpansive, and thus we immediately have the quasi-contractiveness of $\prox_{\gamma f}  \circ (I - \gamma \nabla g)$ from that of $I-\gamma\nabla g$.
%\begin{align*}
%~&\| \prox_{\gamma f}  \circ (I - \gamma \nabla g) (x) - \prox_{\gamma f}  \circ (I - \gamma \nabla g)(x^*) \| \\
%\leq & \|(x - \gamma \nabla g (x)) - (x^* - \gamma \nabla g(x^*))\| \\
%\leq & \sqrt{1 - 2\gamma \mu + {\mu \gamma^2}L} \, \|x - x^*\|.
%\end{align*}
%Hence, $\prox_{\gamma f}  \circ (I - \gamma \nabla g)$ is a quasi-contractive operator. %, which implies that $S = I - \prox_{\gamma f}  \circ (I - \gamma \nabla g)$ is quasi-strongly monotone.
%\qed
\hfill\end{proof}




%\begin{remark}
%Async-parallel FBS requires that $\prox_{\gamma f}$ can be evaluated in a coordinate-wise manner. %% However, when $\prox_{\gamma f}$ is not separable, then evaluating $U_i(\prox_{\gamma f}  \circ (I - \gamma \nabla g)) (x)$ has the same complexity as evaluating $\prox_{\gamma f}  \circ (I - \gamma \nabla g) (x)$, which is not efficient. To avoid this, we can actually consider the backward forward splitting. A more detailed discussion can be found in our following paper \cite{???}.
%\end{remark}

\cut{\commwy{I think we can remove Remarks 1 and 2, and cite our forthcoming 2nd paper here.}
\begin{remark}
When $f = 0$, then~\eqref{eqn:fbs_itr} reduces to~\eqref{eqn:cvx_itr}. A special case is when $f(x)$ is an indicator function of a feasible set $S \subset \cH$. In this case $\prox_{\gamma f} = \Proj_{S}$, which is a projection to the feasible set $S$.
\end{remark}
\begin{remark}
The parallel asynchronous forward backward splitting requires the proximal operator $\prox_{\gamma f}$ can be component-wisely (or block-wisely) evaluated. Examples of those functions include $\ell_1$ norm, $\|x\|_{2,1}$ norm, huber function, $\|x\|_2^2$, and indicator function of box constraints $\cI_{\{x| a \leq x \leq b \} }(x)$. However, when $\prox_{\gamma f}$ is not separable, then evaluating $U_i(\prox_{\gamma f}  \circ (I - \gamma \nabla g)) (x)$ has the same complexity as evaluating $\prox_{\gamma f}  \circ (I - \gamma \nabla g) (x)$, which is not efficient. To avoid this, we can actually consider the backward forward splitting operator
\begin{equation*}
T_{FBS} :=   (I - \gamma \nabla g) \circ \prox_{\gamma f},
\end{equation*}
then we have the following iteration
\begin{equation}\label{eqn:bfs_itr}
z^{k+1} = z^k - \eta_k  \, U_k \Big( \hat z^k -  (I - \gamma \nabla g) \circ \prox_{\gamma f} (\hat{z}^{k}) \Big),
\end{equation}
then the solution can be recovered by $x^{k} = \prox_{\gamma f} (z^k)$. Applications of backward forward splitting doesn't require $\prox_{\gamma f}$ to be separable, however, it requires the proximal operator to be easy to compute, and its complexity is smaller than the evaluating $x_i - \gamma \nabla_i g (x)$, since each update requires evaluating $\prox_{\gamma f}$ once. Applications of backward splitting include $f(x)$ being the indicator function of some constraints such as $\{x ~|~ \|x\|_1 \leq 1\}$, $\{x ~|~ \|x\|_{2} \leq 1\}$, $\{ x ~|~ \sum_{i=1}^n x_i = 1, x_i \geq 0\}$. The proximal operators of those indicator functions can be easily evaluated with low computational complexity. If evaluating forward operator dominants the computational costs, we can expect linear speedup for the asynchronous algorithm.
\end{remark}
}



\end{subsection}

\begin{subsection}{Minimize nonsmooth $+$ nonsmooth functions}
Consider
\begin{equation}\label{fpg}
\Min_{x\in \cH} f(x) + g(x),
\end{equation}
where both $f(x)$ and $g(x)$ are closed proper convex and their $\prox$ maps are easy to compute. Define the Peaceman-Rachford \cite{lions1979splitting} operator:
$$\TPRS :=~\refl_{\gamma f}  \circ~\refl_{\gamma g}.$$
Since both $\refl_{\gamma f}$ and $\refl_{\gamma g}$ are nonexpansive, their composition $\TPRS$ is also nonexpansive.\cut{ Based on $T_{PRS}$ we define the following relaxed PRS operator
$$(T_{PRS})_{\alpha} := (1 - \alpha) \, I + \alpha \, T_{PRS},$$
which is a $\alpha$-averaged operator.} Let
$S :=  I - \TPRS.$
When applying ARock to $T=\TPRS$ to solve problem~\eqref{fpg}, the update~\eqref{eqn:asyn_ortho_update} reduces to:
\beq\label{aprs}z^{k+1} = z^k - \eta_k \, U_{i_k} \circ \big( I - ~\refl_{\gamma f}  \circ~\refl_{\gamma g}\big) \hat{z}^{k},
\eeq
where we use $z$ instead of $x$ since the limit $z^*$ of $(z^k)_{k\ge 0}$ is not a solution  to~\eqref{fpg}; instead, a solution  must be recovered via $x^*=\prox_{\gamma g}z^*$.
The convergence follows from Theorem ~\ref{thm:convergence} and that $\TPRS$ is nonexpansive. If either $f$ or $g$ is strongly convex, then $\TPRS$ is contractive and thus by Theorem~\ref{thm:strongly_monotone}, ARock converges linearly. Finer convergence rates follow from \cite{DavisYin2014,DavisYin2014c}. A naive implementation of~\eqref{aprs} is
\begin{subequations}\label{eqn:relaxed_prs}
\begin{align}
\hat x^{k} &= \prox_{\gamma g} (\hat z^{k}), \label{eqn:relaxed_prs-a} \\
\hat y^{k} & = \prox_{\gamma f} (2 \hat x^{k} - \hat z^{k}),  \label{eqn:relaxed_prs-b}\\
z^{k+1} & = z^k + 2\eta_k \, U_{i_k} (\hat y^{k} - \hat x^{k}),\label{eqn:relaxed_prs-c}
\end{align}
\end{subequations}
where $\hat x^{k}$ and $\hat y^{k}$ are intermediate variables. Note that the order in which the proximal operators are applied to  $f$ and $g$  affects both  $z^k$ \cite{YanYin2014} and whether coordinate-wise updates can be efficiently computed.
%The convergence of Algorithm~\ref{alg:asyn_core}  can be guaranteed by Theorem ~\ref{thm:convergence}.
Next, we present two special cases of~\eqref{fpg} in Subsections~\ref{sec:Feasibility} and~\ref{sec:async-admm} and discuss how to efficiently implement the update~\eqref{eqn:relaxed_prs}.

\begin{subsubsection}{Feasibility problem}\label{sec:Feasibility}
Suppose that $C_1, ..., C_m$ are closed convex subsets of $\cH$ with a nonempty intersection. The problem is to find a point in the intersection. Let $\cI_{C_i}$ be the indicator function of the set $C_i$, that is, $\cI_{C_i}(x)=0$ if $x\in C_i$ and $\infty$ otherwise. The feasibility problem can be formulated as the following
\begin{align}%\label{eqn:feasibility_problem}
\Min_{x=(x_1,\ldots,x_m)\in \cH^m} ~& \sum_{i=1}^m \cI_{C_i}(x_i) + \cI_{\{x_1 = \cdots =x_m\}}(x)  \nonumber.
\end{align}
Let $z^k=(z^k_1,\ldots,z^k_m)\in \cH^m$, $\hat z^k=(\hat z^k_1,\ldots,\hat z^k_m)\in \cH^m$, and $\hat{\bar{z}}^k\in\cH$. We can implement~\eqref{eqn:relaxed_prs} as follows (see Appendix \ref{sec:appendix} for the step-by-step derivation):
\begin{subequations}\label{eqn:relaxed_prs_set}
\begin{align}
\label{eqn:relaxed_prs_set:a}
{\hat {\bar z}}^k &=\textstyle \frac{1}{m} \sum_{i=1}^m \hat z_i^k,  \\
\label{eqn:relaxed_prs_set:b}
\hat{y}_{i_k}^k & = \textstyle\Proj_{C_{i_k}} \big(2 {\hat {\bar z}}^k - \hat z^{k}_{i_k}\big),\\
\label{eqn:relaxed_prs_set:c}
z_{i_k}^{k+1} & = \textstyle z_{i_k}^k + 2\eta_k  (\hat y_{i_k}^{k} - {\hat {\bar z}}^k).
\end{align}
\end{subequations}
The update~\eqref{eqn:relaxed_prs_set} can be  implemented as follows. Let global memory hold $z_1,\ldots,z_m$, as well as  $\bar{z}=\frac{1}{m}\sum_{i=1}^mz_i$. At the $k$th update, an agent independently generates a random number $i_k\in\{1,\ldots,m\}$,  then reads $z_{i_k}$ as $\hat z_{i_k}^k$ and $\bar{z}$ as ${\hat {\bar z}}^k$, and finally computes $\hat{y}_{i_k}$ and updates $z_{i_k}$ in global memory according to~\eqref{eqn:relaxed_prs_set}. Since $\bar{z}$ is maintained in global memory, the agent  updates $\bar{z}$ according to $\bar{z}^{k+1} = \bar{z}^k+\frac{1}{m}(z_{i_k}^{k+1} - z_{i_k}^{k})$. This implementation saves each agent from computing~\eqref{eqn:relaxed_prs_set:a} or reading all $z_1,\ldots,z_m$. Each agent only reads $z_{i_k}$ and $\bar{z}$, executes~\eqref{eqn:relaxed_prs_set:b}, and updates  $z_{i_k}$ \eqref{eqn:relaxed_prs_set:c} and $\bar{z}$. % To improve the efficiency, $\hat w_g^k$ should be updated whenever there is a update to $\hat z_i^k$.
\end{subsubsection}


\end{subsection}

\begin{subsection}{Async-parallel ADMM}\label{sec:async-admm}
This is another application of~\eqref{eqn:relaxed_prs}. Consider %the following problem
%\begin{align}
\beq \label{eqn:admm_form}
\Min_{x \in \cH_1,~ y \in \cH_2}~  f(x) + g(y)\quad\St~  Ax + By = b,
\eeq
%\end{align}
where $\cH_1$ and $\cH_2$ are Hilbert spaces, $A$ and $B$ are bounded linear operators. We apply the update~\eqref{eqn:relaxed_prs} to the Lagrange dual of~\eqref{eqn:admm_form} (see~\cite{gabay1983chapter} for the derivation):
\begin{equation}\label{eq:lag-dual}
\Min_{w \in \cG}\, d_f(w) + d_g(w),
\end{equation}
where $d_f(w) := f^* (A^* w)$, $d_g(w) := g^*(B^* w) - \langle w, b \rangle$, and $f^*$ and $g^*$ denote the convex conjugates of $f$ and $g$, respectively. The proximal maps induced by  $d_f$ and $d_g$ can be computed via solving subproblems that involve only the original terms in~\eqref{eqn:admm_form}:  $z^+=\prox_{\gamma d_f}(z)$ can be computed by (see Appendix \ref{sec:appendix} for the derivation)
\begin{equation}\label{dualf}
\begin{cases}
x^+ \in \arg\min_x f(x) - \langle z, Ax \rangle + \frac{\gamma}{2} \|Ax\|^2, \\
z^+  = z - \gamma Ax^+,
\end{cases}
\end{equation}
and $z^+=\prox_{\gamma d_g}(z)$  by
\begin{equation}\label{dualg}
\begin{cases}
y^+ \in \arg\min_y g(y) - \langle z, By - b\rangle + \frac{\gamma}{2} \|By - b\|^2, \\
z^+  = z - \gamma (By^+ - b).
\end{cases}
\end{equation}
Plugging~\eqref{dualf} and~\eqref{dualg} into (\ref{eqn:relaxed_prs}) yields the following naive implementation
\begin{subequations}\label{eqn:asyn_admm}
\begin{align}
\label{eqn:asyn_admm:a}
\hat y^{k } &\in \arg\min_y g(y) - \langle \hat z^{k}, By - b \rangle + \frac{\gamma}{2}\|By -b\|^2, \\
\label{eqn:asyn_admm:b}
\hat w_{g}^{k} &= \hat z^{k} - \gamma (B \hat y^{k} - b), \\
\label{eqn:asyn_admm:c}
\hat x^{k} &\in \arg\min_x f(x) - \langle 2 \hat w_{g}^{k} -\hat z^{k}, Ax \rangle + \frac{\gamma}{2} \|Ax\|^2,  \\
\label{eqn:asyn_admm:d}
\hat w_{f}^{k} &=2 \hat w_{g}^{k} -\hat z^{k} - \gamma  A \hat x^k,\\
\label{eqn:asyn_admm:e}
z_{i_k}^{k+1} &= z_{i_k}^{k} + \eta_k (\hat w_{f,i_k}^{k} - \hat w_{g,i_k}^{k}).
\end{align}
\end{subequations}
Note that $2\eta_k$ in~\eqref{eqn:relaxed_prs-c} becomes $\eta_k$ in~\eqref{eqn:asyn_admm:e} because ADMM is equivalent to the Douglas-Rachford operator, which is the average of the Peaceman-Rachford operator and the identity operator~\cite{lions1979splitting}.
Under favorable structures,~\eqref{eqn:asyn_admm} can be implemented efficiently.
%In practice, depends on the structure of matrices $A$ and $B$ and  the separability of $f$ and $g$, we shall apply different strategies for async-parallel computation.
For instance, when $A$ and $B$ are block diagonal matrices and $f,g$ are corresponding block separable functions,  steps~\eqref{eqn:asyn_admm:a}--\eqref{eqn:asyn_admm:d} reduce to independent computation for each $i$. Since only $\hat w_{f,i_k}^{k}$ and $\hat w_{g,i_k}^{k}$ are needed to update the main variable $z^k$,  we only need to compute~\eqref{eqn:asyn_admm:a}--\eqref{eqn:asyn_admm:d} for the $i_k$th block. This is exploited in distributed and decentralized ADMM in the next two subsections.
\end{subsection}

\begin{subsubsection}{Async-parallel ADMM for consensus optimization}
Consider the consensus optimization problem:
%\begin{align}
\beq\label{eqn:consensus_problem}
\Min_{x_i, y \in \cH} ~ \textstyle\sum_{i=1}^m f_i(x_i)\quad \St~  x_i - y = 0, \quad\forall i = 1, ..., m,
\eeq
%\end{align}
where $f_i(x_i)$ are proper close convex functions. Rewrite~\eqref{eqn:consensus_problem} to the ADMM form:
 \begin{align}\label{eqn:consensus_admm_form}
\Min_{x_i, y \in \cH} ~~& \textstyle\sum_{i=1}^m f_i(x_i)  + g(y)\nonumber \\
                        \St~ &
                        \begin{bmatrix} I_{\cH} & 0 & \cdots &0\\
                                0 & I_{\cH} &\cdots & 0\\
                                & & \ddots & \\
                                0 & 0 &\cdots & I_{\cH}
                        \end{bmatrix}
                        \begin{bmatrix} x_1 \\ x_2\\ \vdots \\x_m\\
                        \end{bmatrix}
                        -
                        \begin{bmatrix} I_{\cH} \\ I_{\cH}\\ \vdots \\ I_{\cH}\\
                        \end{bmatrix}
                        y = 0,
\end{align}
where $g = 0$. Now apply the async-parallel ADMM~\eqref{eqn:asyn_admm} to~\eqref{eqn:consensus_admm_form} with dual variables  $z_1, ..., z_m \in \cH$. In particular, the update~\eqref{eqn:asyn_admm:a}, \eqref{eqn:asyn_admm:b}, \eqref{eqn:asyn_admm:c}, \eqref{eqn:asyn_admm:d} reduce to
\begin{align}\label{eqn:async-par-final}
\hat{y}^k&=\textstyle\argmin_y \big\{ \sum_{i=1}^m \langle \hat z_i^k, y \rangle + \frac{\gamma m}{2} \|y\|^2 \big\}= -\frac{1}{\gamma m} \sum_{i=1}^m \hat z_i^k \cr
(\hat w_{d_g}^{k})_i &= \hat z_i^{k} + \gamma \, \hat y^k \\
\hat x_i^{k} &= \textstyle\argmin_{x_i} \big\{ f_i(x_i) - \langle 2 (\hat w_{d_g}^{k})_i -\hat z_i^{k}, x_i \rangle + \frac{\gamma}{2} \|x_i\|^2 \big\}, \cr
(\hat w_{d_f}^{k})_i &=2 (\hat w_{d_g}^{k})_i -\hat z_i^{k} - \gamma \, \hat x_i^k \nonumber
\end{align}
Therefore, we obtain the following async-parallel ADMM algorithm for  the  problem~\eqref{eqn:consensus_problem}. This algorithm applies to all the distributed applications  in~\cite{boyd2011distributed}. %can all be solved with our async-parallel ADMM model~\eqref{eqn:consensus_problem}.


\begin{algorithm}[H]\label{alg:asyn_consensus}
%\SetAlgoLined
\SetKwInOut{Input}{Input}\SetKwInOut{Output}{output}
\Input{set shared variables $y^0, z_i^0, ~\forall i$, and $K>0$.}
 \While{$k < K$ every agent asynchronously and continuously}{
  choose $i_k$ from $\{1, ..., m\}$ with equal probability\;
  %read ${z_{i_k}^k}$ and $y^k$\;
  evaluate $ (\hat w_{d_g}^{k})_{i_k}$, $\hat x_{i_k}^k$, and $(\hat w_{d_f}^{k})_{i_k}$ following~\eqref{eqn:async-par-final}\;
  update $z_{i_k}^{k+1} = z_{i_k}^{k} +  \eta_k \,  ((\hat w_{d_f}^{k})_{i_k} - (\hat w_{d_g}^{k})_{i_k})$\;
  update $ y^{k+1} =  y^k + \frac{1}{\gamma m} (z_{i_k}^k - z_{i_k}^{k+1})$\;
  update the global counter $k \leftarrow k+1$\;
 }
 \caption{ARock for consensus optimization}
\end{algorithm}
%We see that the async-parallel ADMM enables the objective and constraints being distributed across multiple agents.


\begin{subsubsection}{Async-parallel ADMM for decentralized optimization}\label{subsec:async_admm_des}
Let $V=\{1, ..., m\}$ be a set of agents and $E = \{(i, j) ~|~ \text{if agent $i$ connects to agent $j$}, i < j\}$ be the set of undirected links between the agents.  Consider the following decentralized consensus optimization problem on the graph $G= (V, E)$:
\begin{equation}\label{eqn:des_prob}
%\begin{aligned}
\Min_{x_1, ..., x_m \in \mathbb{R}^d}  ~ \textstyle f(x_1, \ldots, x_m) := \sum_{i=1}^m f_i (x_i),\quad
                \St~  x_i = x_j, ~ \forall (i,j) \in E,
%\end{aligned}
\end{equation}
%Similar to the setting in Section~\ref{sec:dgd}, we consider that $n$ agents form a connected communication network, in which they collaboratively solve the following decentralized consensus optimization problem
%\begin{equation}\label{eqn:des_prob}
%\Min_{x_1, ..., x_p} ~\sum_{i=1}^p f_i (x_i),
%\end{equation}
where $x_1, ..., x_m \in \mathbb{R}^d$ are the local variables and each agent can only communicate with its neighbors in $G$. By introducing the auxiliary variable $y_{ij}$ associated with each edge $(i, j) \in E$, the problem~\eqref{eqn:des_prob} can be reformulated as:
\begin{equation}\label{eqn:des_prob2}
%\begin{aligned}
\Min_{x_i, y_{ij}} ~\textstyle\sum_{i=1}^m f_i (x_i),\quad \St~  x_i = y_{ij}, ~x_j = y_{ij}, \quad \forall (i, j) \in E.
%\end{aligned}
\end{equation}
Define $x = (x_1, ..., x_m)^T$ and $y = (y_{ij})_{(i, j) \in E} \in \mathbb{R}^{|E|d}$ to rewrite~\eqref{eqn:des_prob2} as
\begin{equation}\label{eqn:des_prob3}
%\begin{aligned}
\Min_{x, y} ~ \textstyle\sum_{i=1}^m f_i (x_i),\quad \St ~ Ax + By = 0,
%\end{aligned}
\end{equation}
for proper matrices $A$ and $B$.
%\cut{where $A \in \mathbb{R}^{2|E|d \times md}$ and $B = [-I_{|E|d}; -I_{|E|d}]$ with $I_{|E|d}$ being the $|E|d \times |E|d$ identity matrix.}
Applying the async-parallel ADMM~\eqref{eqn:asyn_admm} to~\eqref{eqn:des_prob3} gives rise to the following simplified update:
Let $E(i)$ be the set of edges connected with agent $i$ and $|E(i)|$ be its cardinality. Let $L(i) = \{j ~|~ (j, i) \in E(i), j < i\}$ and $R(i) = \{j ~|~ (i, j) \in E(i), j > i\}$. To every pair of constraints $x_i = y_{ij}$ and $x_j = y_{ij}$, $(i,j)\in E$, we associate the dual variables $z_{ij,i}$ and $z_{ij,j}$, respectively. Whenever some agent $i$ is activated, it calculates
\begin{subequations}\label{eqn:asyn_admm_des}
\begin{align}
\label{eqn:asyn_admm_des_a}
 \hat x^{k}_{i}  = &\argmin_{x_i} f_i(x_i) +  \big(\sum_{l \in L(i)}\hat z^k_{li, l}  + \sum_{r\in R(i)}\hat z^k_{ir, r} \big) x_i +  \frac{\gamma}{2} |E(i)|\cdot \|x_i\|^2,  \\\label{eqn:asyn_admm_des_b}
 z_{li, i}^{k+1} = &z_{li, i}^k -   \eta_k (( \hat z_{li, i}^k + \hat z_{li, l})/2 +  \gamma \hat x^k_i), \quad\forall l \in L(i),\\
z_{ir, i}^{k+1} = &z_{ir, i}^k -  \eta_k (( \hat z_{ir, i}^k + \hat z_{ir, r})/2 +  \gamma \hat x^k_i), \quad\forall r \in R(i).
\label{eqn:asyn_admm_des_c}
\end{align}
\end{subequations}
%Equation~\eqref{eqn:asyn_admm_des} is truly distributed.
We present the algorithm based on
\eqref{eqn:asyn_admm_des} for problem~\eqref{eqn:des_prob} in Algorithm~\ref{alg:asyn_decentralize}.


%%
%%\cut{Specifically, update~\eqref{eqn:asyn_admm:a} gives $\hat y^{k } \in \arg\min_y   \langle -\hat z^{k}, By  \rangle + \frac{\gamma}{2}\|By\|^2$,
%%where $z = (z_{ij, i}, z_{ij, j})_{(i,j)\in E} \in \mathbb{R}^{2|E|d}$ is the dual variable with $z_{ij, s} (s \in \{i, j\})$  representing the dual variable associate with the constraint $x_s = y_{ij}$.  Since $B=[-I_{|E|d}; -I_{|E|d}]$, we have
%%\begin{equation}\label{eqn:ayn_admm_1}
%%\hat y^k_{ij} = \textstyle - \frac{1}{4 \gamma} (\hat z_{ij,i}^k+\hat z_{ij,j}^k), ~\forall (i, j) \in E
%%\end{equation}
%%Update~\eqref{eqn:asyn_admm:b} becomes
%%\begin{equation}\label{eqn:ayn_admm_2_b}
%%(\hat w_{g}^{k})_{ij, s} = \hat z^{k}_{ij, s} + \gamma \,  \hat y^{k}_{ij}, ~\forall (i, j) \in E, s \in \{i, j\}.
%%\end{equation}
%%Substituting~\eqref{eqn:ayn_admm_1} to~\eqref{eqn:ayn_admm_2_b} gives
%%\begin{equation}\label{eqn:ayn_admm_2b}
%%(\hat w_{g}^{k})_{ij, i} = \frac{3}{4} \hat z_{ij, i}^k - \frac{1}{4} \hat z_{ij, j}^k, \quad(\hat w_{g}^{k})_{ij, j} = \frac{3}{4} \hat z_{ij, j}^k - \frac{1}{4} \hat z_{ij, i}^k, ~ \forall (i, j) \in E.
%%\end{equation}
%%
%%
%%Let $E(i)$ be the set of edges connected with agent $i$, and $|E(i)|$ be the cardinality of $E(i)$. Let $L(i) = \{j ~|~ (j, i) \in E(i), j < i\}$ and $R(i) = \{j ~|~ (i, j) \in E(i), j > i\}$. The update~\eqref{eqn:asyn_admm:c} can be evaluated component-wisely as following
%%\begin{align}\label{eqn:des_x_update}
%%\hat x^{k}_{i} \in &\arg\min_{x_i}  f_i(x_i) + \sum_{e\in E(i)} (\hat z^k_{e, i} -2 (\hat w_g^k)_{e, i})  \, x_i + \frac{\gamma}{2} |E(i)|^2 \|x_i\|^2 \nonumber \\
%% ~ \in &\arg\min_{x_i} f_i(x_i) + \sum_{l \in L(i)}(\hat z^k_{li, i} - 2 (\hat w_g^k)_{li,i} ) x_i + \sum_{r\in R(i)}(\hat z^k_{ir, i} - 2 (\hat w_g^k)_{ir,i} ) x_i +  \frac{\gamma}{2} |E(i)|^2 \|x_i\|^2  \nonumber \\
%% ~ \in &\arg\min_{x_i} f_i(x_i) + \frac{1}{2} \sum_{l \in L(i)}(\hat z^k_{li, l} - \hat z^k_{li,i} ) x_i + \frac{1}{2} \sum_{r\in R(i)}(\hat z^k_{ir, r} - \hat z^k_{ir,i} ) x_i +  \frac{\gamma}{2} |E(i)|^2 \|x_i\|^2.
%%\end{align}
%%\cut{Update~\eqref{eqn:asyn_admm:d} can also be written component-wisely as
%%
%%\begin{equation}
%%(\hat w_{f}^{k})_{ij, s} = 2(\hat w_{g}^k)_{ij, s}-\hat z^{k}_{ij, s} - \gamma \, \hat x^k_s, ~\forall (i, j)\in E, s \in \{i, j\}.
%%\end{equation}
%%Finally, for~\eqref{eqn:asyn_admm:e}, we have
%%\begin{equation}
%%z_{i_k j_k, s_k}^{k+1} = z_{i_k j_k, s_k}^{k} + 2\alpha \eta_k \, (\hat w_{f}^{k})_{i_k j_k, s_k}.
%%\end{equation}
%%}
%%Combining~\eqref{eqn:asyn_admm:d} with~\eqref{eqn:asyn_admm:e}, we have
%%\begin{equation}\label{eqn:update_z}
%%z_{ij, i}^{k+1} = z_{ij, i}^k + 2 \alpha \eta_k (\frac{1}{2} \hat z_{ij, i}^k - \frac{1}{2} \hat z_{ij, j} - \gamma \hat x^k_i), \quad
%%z_{ij, j}^{k+1} = z_{ij, j}^k + 2 \alpha \eta_k (\frac{1}{2} \hat z_{ij, j}^k - \frac{1}{2} \hat z_{ij, i} - \gamma \hat x^k_j), \forall (i, j) \in E.
%%\end{equation}
%%}


\begin{algorithm}[H]\label{alg:asyn_decentralize}
\SetKwInOut{Input}{Input}\SetKwInOut{Output}{output}
\Input{Each agent $i$ sets the dual variables $z^0_{e, i}=0$ for  $e \in E(i)$, $K>0$.}
 \While{$k < K$, any activated agent $i$}{
 (previously received $\hat z_{li, l}^k$ from neighbors $l \in L(i)$ and $\hat z_{ir,r}^k$ from  $r \in R(i)$)\;
 update $\hat x_i^{k}$ according to~\eqref{eqn:asyn_admm_des_a}\;
 update $z_{li, i}^{k+1}$ and $z_{ir, i}^{k+1}$ according to~\eqref{eqn:asyn_admm_des_b} and~\eqref{eqn:asyn_admm_des_c}, respectively\;
 send $z_{li, i}^{k+1}$ to neighbors $l \in L(i)$ and  $z_{ir, i}^{k+1}$ to neighbors $r \in R(i)$\;
 }
 \caption{ARock for the decentralized problem~\eqref{eqn:des_prob2}}
\end{algorithm}
Algorithm~\ref{alg:asyn_decentralize} activates one agent at each iteration and updates all the dual variables associated with the agent. In this case, only one-sided communication is needed, for sending the updated dual variables in the last step. We allow this communication to be delayed in the sense that agent $i$'s neighbors may be activated and start their computation before receiving the latest dual variables from agent $i$.

Our algorithm is different from the asynchronous ADMM algorithm by Wei and Ozdaglar~\cite{wei2013on}. Their algorithm activates an edge and its two associated agents at each iteration and thus requires two-sided communication at each activation.
We can recover their algorithm as a special case by activating an edge $(i,j)\in E$ and its associated agents $i$ and $j$ at each iteration, updating  the  dual variables $z_{ij,i}$ and $z_{ij,j} $ associated with the edge, as well as computing the intermediate variables $x_i$, $x_j$, and $y_{ij}$. The updates are derived from~\eqref{eqn:des_prob3} with the orders of $x$ and $y$ swapped. Note that \cite{wei2013on} does not consider the situation that adjacent edges are  activated in a short period of time, which may cause overlapped computation and delay communication. Indeed, their algorithm corresponds to $\tau=0$ and the corresponding stepsize $\eta_k\equiv1$. Appendix \ref{sec:appendix2} presents the steps to derive the algorithms in this subsection.
\end{subsubsection}

\cut{
\begin{subsubsection}{The generalized lasso}
The generalized Lasso is the following model
\begin{equation}
\min_{x} \lambda \, \|Dx\|_1 + \frac{1}{N}  \sum_{i=1}^N \ell_{i}(x),
\end{equation}
where $N$ is the number of samples, $x \in \mathbb{R}^n$ is the feature vector, $\ell_i (x)$ is the loss function, such as square loss, logistic loss or hinge loss. $D$ is the penalty matrix specifying the desired structured sparsity pattern of $x$. Different setting of $D$ reduces to different models such as fused lasso, trend filtering, and graph-guided fused lasso. The asynchronous ADMM can be applied by letting
$$f_i(x) =  \frac{\lambda}{M} \|Dx\|_1 + \frac{1}{N}  \sum_{j= i N/M}^{(i+1)N/M} \ell_{j}(x),$$
where $M$ is the number of workers. We assign $N/M$ samples to each worker. In practice, since the subproblems do not have closed form solution, they can be solved by using some inexact ADMM method \cite{zhang2011unified}.
\end{subsubsection}
}


\end{subsubsection}

\cut{
\begin{subsection}{Asynchronous Proximal Jacobi ADMM}
Consider the following problem,
\begin{align}\label{eqn:jacobi_admm_problem}
\min_{x_i, x} & \sum_{i=1}^m g_i(x_i) + f(x) \nonumber \\
                        \text{s.t.}~~ & A_1 x_1 + A_2 x_2 + \cdots + A_m x_m = b,
\end{align}
where $x_i \in \cH_i$, $x = (x_1, ..., x_m)^T \in \cH$, $b \in \cH_l$, $f(x):\cH \rightarrow R$ is a convex differentiable function, $g_i(x_i): \cH_i \rightarrow R$ are closed proper convex functions and $A_i: \cH_i \rightarrow \cH_l$ are linear operators. Define the Lagrangian function
\begin{equation}
L(x, \lambda) = \sum_{i=1}^m f_i(x_i) + g(x)  - \langle x_{m+1}, \sum_{i=1}^m A_i x_i - b \rangle,
\end{equation}
where $x_{m+1}$ is the Lagrangian multiplier. The optimality condition for~\eqref{eqn:jacobi_admm_problem} is the following
\begin{align}\label{eqn:jacobi_admm_kkt}
                        \begin{bmatrix} 0 \\ 0\\ \vdots \\ 0 \\ 0 \\
                        \end{bmatrix}
                         \in
                        \begin{bmatrix}
                        \nabla_1 f(x) + \partial g_1(x_1) \\
                        \nabla_2 f(x) + \partial g_2(x_2) \\
                        \vdots\\
                        \nabla_m f(x) + \partial g_m(x_m) \\
                        0 \\
                        \end{bmatrix}
                        +
                        \underbrace{
                        \begin{bmatrix} 0 & 0 & \cdots &0 &-A_1^*\\
                                0 & 0 &\cdots &0& -A_2^*\\
                                \vdots& \vdots& \vdots & \vdots &\vdots\\
                                0 & 0 &\cdots &0& -A_m^* \\
                                A_1 & A_2 & \cdots &A_m& 0 \\
                        \end{bmatrix}
                        }_{Q}
                        \begin{bmatrix} x_1 \\ x_2\\ \vdots \\x_m\\ x_{m+1} \\
                        \end{bmatrix}
                        -
                        \begin{bmatrix} 0 \\ 0\\ \vdots \\ 0 \\ b \\
                        \end{bmatrix}.
\end{align}
We know that $x$ is the optimal solution of~\eqref{eqn:jacobi_admm_problem} if and only if $x_1, ..., x_m, x_{m+1}$ satisfies~\eqref{eqn:jacobi_admm_kkt}. Define $z := (x_1, ..., x_m, x_{m+1})$, $F(z) := f(z)$, $G(z) := \sum_{i=1}^m g_i (x_i)$ and $Q$ is defined as shown in~\eqref{eqn:jacobi_admm_kkt}, then finding $z$ satisfying~\eqref{eqn:jacobi_admm_kkt} is equivalent to the following monotone inclusion problem
\begin{equation}\label{eqn:jacobi_admm_monotone}
\text{find $z$ such that } \quad 0 \in \nabla F(z) + \partial G(z) + Qz.
\end{equation}
Introduce a symmetric positive definite matrix $W$, and apply the forward backward splitting operator to this problem, then the monotone inclusion problem~\eqref{eqn:jacobi_admm_kkt} is equivalent to the following fixed point problem
\begin{equation}\label{eqn:jacobi_admm_fixed_point}
\text{find $z$ such that } z = J_{W^{-1} (\partial G + Q)} \big(z - W^{-1} \nabla F(z) \big) ,
\end{equation}
where $J_{A} := (I_{\cH} + A)^{-1}$ is called the resolvent of operator $A$. Let $T z = J_{W^{-1} (\partial G + Q)} \big(z - W^{-1} \nabla F(z) \big)$, we know that $T$ is a $\alpha$-averaged operator in the norm $\|\cdot\|_W$. Then we can let $S = I - T$, and apply the parallel asynchronous algorithm with $U_k$ as defined in Remark~\ref{remark:pd}, we have the following iteration
\begin{equation}\label{eqn:asyn_jacobi_admm_fixed_point}
z^{k+1} = z^k - \eta_k \, U_k \Big( \hat z^k - J_{W^{-1} (\partial G + Q)} \big(\hat z^k - W^{-1} \nabla F(\hat z^k) \big) \Big).
\end{equation}
Now we show that iteration~\eqref{eqn:asyn_jacobi_admm_fixed_point} corresponds to an asynchronous version of the proximal Jacobi ADMM method. Let
\begin{equation}\label{eqn:def_w}
W =
        \begin{bmatrix}
        P - \sigma A^T A & 0 \\
        0         & \frac{1}{\sigma} I \\
        \end{bmatrix},
\end{equation}
where $P = diag(P_1, ..., P_m)$, $P_i: \cH_i \rightarrow \cH_i$ are easy to invert symmetric matrices and $A = (A_1, ..., A_m)$. For convergence purposes, we require $P \succ \sigma A^T A$ so that $W$ is symmetric positive definite. The reason for choosing such a $W$ is to decouple $x_1, ..., x_m$ and $\lambda$. Define
\begin{equation}
\hat{y}^k := J_{W^{-1} (\partial G + Q)} \big(\hat z^k - W^{-1} \nabla F(\hat z^k) \big),
\end{equation}
where $\hat y^k = (\hat y_1^k, ..., \hat y_m^k, \hat y_{m+1}^k)$,
then we have
\begin{equation}\label{eqn:update_y}
W \hat y^k = W \hat z^k - \nabla F(\hat z^k) -  \tilde \nabla G(\hat y^k) - Q \hat y^k,
\end{equation}
where $\tilde \nabla G(\hat y^k) \in \partial G(\hat y^k)$. Substituting~\eqref{eqn:def_w} to~\eqref{eqn:update_y} gives the following updates
\begin{equation}\label{eqn:update_y_i}
\begin{cases}
\hat p^k = \sum_{j=1}^m A_j \hat x_j^k  \\
\hat y_i^k = \argmin_{x_i}  g_i(x_i) + \big\langle \nabla_i f(\hat x^k), x_i - \hat x_i^k \big\rangle + \big \langle \sigma A_i^T (p^k - b - \frac{\hat x_{m+1}^k}{\sigma}), x_i \big\rangle + \frac{1}{2} \|x_i - \hat x_i^k\|_{P_i}^2, \forall i = 1, ..., m \\
\hat y_{m+1}^k = \hat x_{m+1}^k - \sigma \sum_{i=1}^m A_i \hat y_i^k
\end{cases}
\end{equation}
{\color{red} [$\hat y^{k}_{m+1}$ depends on the $y^k_{i}$, for $i=1, ..., m$]}
If we let $P_i$ be a diagonal matrix, assume $\nabla_i f$ can be evaluated without calculating $\nabla f$, and assume $g_i(x_i)$ is simple functions, then the first step has closed form solution. Substituting $\hat y^k$ to~\eqref{eqn:asyn_jacobi_admm_fixed_point}, we have the following updates
\begin{equation}
z^{k+1} = z^k - \eta_k \, U_k (\hat z^k - \hat y^k).
\end{equation}

To conclude, we have the following asynchronous proximal Jacobi ADMM algorithm for solving~\eqref{eqn:jacobi_admm_problem}.

\begin{algorithm}[H]\label{alg:asyn_jacobi_admm}
\SetAlgoLined
\SetKwInOut{Input}{Input}\SetKwInOut{Output}{output}
\Input{Initial points $x_i^0, \lambda^0, y_i^0, p^0 = \sum_{i=1}^m A_i x_i^0$ are shared variables, number of iterations $K$}
 \While{$k < K$}{
  choose $i_k$ from $\{1, ..., m, m+1\}$ with equal probability\;
  read ${\hat x_{i_k}^k}$ and $\hat p^k$ from global memory\;
  evaluate $\hat y_{i_k}^k$ according to~\eqref{eqn:update_y_i}\;
  update $x_{i_k}^{k+1} =x_{i_k}^k - \eta_k \, (\hat x_{i_k}^k - \hat y_{i_k}^k)$\;
  \If{$i_k\leq m$}{
  update $p^{k+1} = p^k + A_{i_k} (x_{i_k}^{k+1} - \hat x_{i_k}^k)$\;
  }
  $k \leftarrow k+1$\;
 }
 \caption{Asynchronous proximal Jacobi ADMM for consensus optimization}
\end{algorithm}
Solve linear programming problem in standard form
\begin{align}\label{eqn:jacobi_admm_problem}
\min_{x} ~&   c^T x \nonumber \\
                        \text{s.t.}~ & A x = b \\
                                            & x \geq 0. \nonumber
\end{align}

Solve the Basis pursuit problem
\begin{align}\label{eqn:jacobi_admm_problem}
\min_{x} ~&  \|x\|_1 \nonumber \\
                        \text{s.t.}~ & A x = b.
\end{align}

Consider a classification problem. Given a datasets $\{(a_1, b_1), ..., (a_m, b_m)\}$ where $(a_i, b_i), i=1,..., m, a_i \in \mathbb{R}^n, b_i \in \{-1, 1\}$ are instance-label pairs. We use the following dual SVM method to train a classifier,
\begin{align}\label{eqn:jacobi_admm_problem}
\min_{x} ~& \frac{1}{2} x^T A x - 1^T x \nonumber \\
                        \text{s.t.}~ & \sum_{i=1}^n b_i x_i = 0 \\
                                            & 0 \leq x \leq C, \nonumber
\end{align}
where $a_{ij} = b_i b_j k(a_i, a_j)$, where $k(\cdot, \cdot)$ is a kernel function. All of the previously listed problems can be solved by using Algorithm~\ref{alg:asyn_jacobi_admm}.
\end{subsection}
}
\end{section}


%\input{asyn_t_analysis}
\cut{\begin{section}{Assumptions}
To derive the convergence for Algorithm~\ref{alg:asyn_core}, we make the following two assumptions.

\begin{assumption}[Bounded delay] We assume that there is an upper bound $\tau$ on the delay, that is, $J(k)\subset \{k-1,k-2,\cdots,k-\tau\}$. % on the age of each component in $\hat{x}^k$. %Specifically, for consistent read, we assume $ 0 \leq j_k \leq \tau$. For inconsistent read, for all $j \in  J(k)$, we assume $ k - \min\{j ~|~ j \in  J(k) \} \leq \tau$.
\end{assumption}

%\begin{assumption}{($\alpha$-averaged operator)}
%We assume that operator $T: \cH \rightarrow \cH$ is an $\alpha$-averaged operator where $\alpha \in (0, 1)$.
%\end{assumption}

%\begin{assumption}{(Uniform sampling)}
%We assume that the block of coordinates $x_{i_k}$ at each iteration is sampled with equal probability.
%\end{assumption}

\begin{assumption}[Sampling order vs updating order]
We assume that the sampling order is the same as the updating order.
\end{assumption}

\begin{remark}
Sampling order can be different from the actual updating order. However, if the delay is independent with the random choice of blocks, the actual block of coordinates will still be updated with the same probability. Hence, this assumption does not lose generality. %Without loss of generality, the algorithm assumes that the sampling order is the same as the updating order.
\end{remark}

%\begin{remark}[Non-uniformly sampling]
%The block of coordinates $x_{i_k}$ at each iteration does not have to be sampled uniformly. The results in this paper apply to the more general case where the block of coordinate $x_{i_k}$ is sampled following a fixed strictly positive distribution. If the $i$th block of coordinates is sampled with a positive probability $p_i$, we can change $\tau_k$ such that $\tau_k p_{i_k}$ is fixed.
%\end{remark}

%\begin{subsection}{Notations.}
%
%The following is a list of notations used in this paper.
%\begin{itemize}
%\item $\cH$ -- a separable real Hilbert space;
%\item $\cH^n$ -- $n$ product space of $\cH$;
%\item  $\cB$ --- Borel $\sigma$-algebra of {$\cH$};
%\item  $(x^k)\subset\cH$ ---  a sequence of points in $\cH$;
%\item $x_i$ --- the $i$th block of $x \in \cH$, and $x = (x_1, ..., x_m)^T$;
%\item $X^*$ --- the set of fixed points of $T$;
%\item $\tilde{X}^{*} := \{(x^*, x^*, ..., x^*)^T ~|~ x^* \in X^*\}\subseteq \cH^{\tau+1}$. Each element is $\tau + 1$ copies of $x^*$;
%\item  $\mathscr{Z}(x^k)$ --- the set of weak sequential cluster points of the sequence $(x^k)$;
%\item $\cX^k := \sigma (x^0, x^1, ..., x^k)$ --- the smallest $\sigma$-algebra generated by $x^0, x^1, ..., x^k$;
%\item$\mathscr{X} := (\cX^n)$, a sequence of sigma-algebras;
%\item $\Omega$ ---   sample space;
%\item $\cF$ --- $\sigma$-algebra of $\Omega$;
%\item $\cF^k$ --- sub-sigma algebra of $\cF$;
%\item  $P$ ---  probability measure;
%\item $(\Omega, \cF, P)$  ---  underlying probability space;
%\item $x: (\Omega, \cF) \rightarrow (\cH, \cB)$ --- $\cH$-valued random variable;
%\item $\mathscr{F} = (\cF^k)$ --- a sequence of sub-sigma algebras of $\cF$ satisfying $\cF^k \subset \cF^{k+1}$, $\forall k$;
%\item $\ell_+ (\mathscr{F})$ --- the set of sequences of $[0, +\infty)$-valued random variables $(\xi_k)$ where $\xi^k$ is $\cF^k$ measurable;
%\item ``a.s.'' ---  ``almost surely'';
%\item $\ell_+^p (\mathscr{F}) := \big\{(\xi_k) \in \ell_+ (\mathscr{F})| \sum_{k} (\xi_k)^p < + \infty~ \text{a.s.} \big\}$;
%\item $\odot$ --- pointwise product.
%\end{itemize}
%\rev{Superscripts are iteration counters except for scalars such as $\xi_k$, subscripts are used.}
%\end{subsection}

\end{section}
}


\begin{section}{Convergence}
%This section provides convergence analysis of ARock (Algorithm~\ref{alg:asyn_core}).
We establish weak and strong convergence  in Subsection~\ref{sec:weak} and  linear convergence in Subsection~\ref{sec:linear}. Step size selection is also discussed.

\subsection{Almost sure convergence}\label{sec:weak}
%In this subsection, we show the almost sure weak and strong convergence of Algorithm~\ref{alg:asyn_core}.
\begin{assumption} Throughout the our analysis, we assume $p_{\min}:=\min_ip_i>0$ and
\begin{equation}\label{eq:independence}
\Prob(i_k=i\,|\,\cX^k)=\Prob(i_k=i)=p_i,\quad\forall i,k.
\end{equation} 
\end{assumption}
We let $|J(k)|$ be the number of elements in $J(k)$ (see Subsection~\ref{subset:aync_mem}).
Only for the purpose of analysis, we define the (never computed) full update at $k$th iteration:
\begin{align}\label{eqn:def_bar_x}
\bar{x}^{k+1} := x^k - \eta_k S\hat{x}^{k}.
\end{align}
Lemma~\ref{lemma:a-avg} below shows that $T$ is nonexpansive if and only if $S$ is ${1/2}$-cocoercive.
\begin{lemma}\label{lemma:a-avg}
Operator $T:\cH\to\cH$ is nonexpansive if and only if $S = I - T$ is {${1}/{2}$-cocoercive}, i.e., $\langle x-y,Sx- Sy\rangle \ge \frac{1}{2}\|Sx- Sy\|^2, {\forall~x,y\in\cH}$.
\end{lemma}
\begin{proof}
See textbook \cite[Proposition 4.33]{bauschke2011convex} for the proof of the ``if'' part, and the ``only if'' part, though missing there, follows by just reversing the proof.
\hfill\end{proof}

The lemma below develops an an upper bound for the expected distance between $x^{k+1}$ and any $x^*\in \Fix T$.
\begin{lemma}\label{lemma:fund}
Let $(x^k)_{k\geq 0}$ be the sequence generated by Algorithm \ref{alg:asyn_core}. Then for any $x^*\in \Fix T$ and $\gamma>0$ (to be optimized later), we have % for both consistent read and inconsistent read
\begin{align}\label{eqn:fund_inquality0}
\begin{aligned}
\textstyle\mathbb{E}\big(\|x^{k+1} - x^* \|^2 \,\big|\, \cX^k\big)  \leq & \textstyle\|x^{k} - x^* \|^2  +\frac{\gamma}{m}\sum_{d\in J(k)}\|x^d-x^{d+1}\|^2\\
&\textstyle+ \frac{1}{m}\left(\frac{{|J(k)|}}{\gamma}+\frac{1}{mp_{\min}}-\frac{1}{\eta_k}\right)\|x^k-\bar x^{k+1}\|^2.
\end{aligned}
\end{align}
%where $\mathbb{E}(\cdot\,|\,\cX^k)$ denotes conditional expectation on $\cX^k$.
\end{lemma}

\begin{proof}
Recall $\Prob(i_k=i)=p_i$. Then we have
\begin{equation}\label{eqn:equality_inconsistent}
\begin{array}{rcl}
&&\mathbb{E}\left(\|x^{k+1} - x^*\|^2\,|\,\cX^k\right)\\
&\overset{\eqref{eqn:asyn_update}}=&\mathbb{E}\left(\|x^{k}  - \textstyle\frac{\eta_k}{mp_{i_k}}S_{i_k} \hat{x}^{k}- x^* \|^2\,|\,\cX^k\right)\cr
&=& \|x^k - x^*\|^2+\mathbb{E}\left(\textstyle\frac{2\eta_k}{mp_{i_k}} \left\langle S_{i_k} \hat{x}^{k}, x^* - x^k \right\rangle + \textstyle\frac{\eta_k^2}{m^2p_{i_k}^2} \|S_{i_k}\hat{x}^{k}\|^2\,\big|\,\cX^k\right)\cr
&\overset{\eqref{eq:independence}}=&\|x^k - x^*\|^2+\textstyle\frac{2\eta_k}{m} \sum_{i=1}^m\left\langle S_i \hat{x}^{k}, x^* - x^k \right\rangle + \frac{\eta_k^2}{m^2}\sum_{i=1}^m\frac{1}{p_i}\|S_i\hat{x}^{k}\|^2\cr
&=&\|x^k - x^*\|^2+\textstyle\frac{2\eta_k}{m} \left\langle S \hat{x}^{k}, x^* - x^k \right\rangle +\frac{\eta_k^2}{m^2} \sum_{i=1}^m\frac{1}{p_i} \|S_i\hat{x}^{k}\|^2.
\end{array}
\end{equation}
%where in the third equality, we have used \eqref{eq:independence}. %holds because the probability of choosing $i$ is $p_i$.

Note that
\begin{equation}\label{term2}\sum_{i=1}^m\frac{1}{p_i} \|S_i\hat{x}^{k}\|^2\le\frac{1}{p_{\min}} \sum_{i=1}^m\|S_i\hat{x}^{k}\|^2
=\frac{1}{p_{\min}} \|S\hat{x}^{k}\|^2\overset{\eqref{eqn:def_bar_x}}{=}\frac{1}{\eta_k^2p_{\min}}\|x^k-\bar{x}^{k+1}\|^2,\end{equation}
%\begin{equation}\label{term2}
%\begin{aligned}
%\textstyle\sum_{i=1}^m\frac{1}{p_i} \|S_i\hat{x}^{k}\|^2&\le\textstyle\frac{1}{p_{\min}} \sum_{i=1}^m\|S_i\hat{x}^{k}\|^2\\
%&\textstyle=\frac{1}{p_{\min}} \|S\hat{x}^{k}\|^2\overset{\eqref{eqn:def_bar_x}}{=}\begin{array}{l}\frac{1}{\eta_k^2p_{\min}}\end{array}\|x^k-\bar{x}^{k+1}\|^2,
%\end{aligned}
%\end{equation}
and
\begin{align}\label{term1}
%\begin{aligned}
&\langle S \hat{x}^{k}, x^* - x^k \rangle\cr
\overset{\eqref{eqn:inconsist}}=&\textstyle\langle S \hat{x}^{k},x^* - \hat{x}^k + \sum_{d\in J(k)} (x^{d} - x^{d+1})\rangle\cr
\overset{\eqref{eqn:def_bar_x}}=&\textstyle\langle S \hat{x}^{k}, x^* - \hat{x}^k\rangle + \frac{1}{\eta_k}\sum_{d\in J(k)}\langle x^k-\bar{x}^{k+1}, x^{d} - x^{d+1}\rangle\cr
\le&\textstyle \langle S \hat{x}^{k}-Sx^*, x^* - \hat{x}^k\rangle+\frac{1}{2\eta_k}\sum_{d\in J(k)}\big(\frac{1}{\gamma}\|x^k-\bar{x}^{k+1}\|^2+ \gamma\|x^{d} - x^{d+1}\|^2\big)\\
\le&\textstyle -\frac{1}{2}\|S \hat{x}^{k}\|^2+\frac{1}{2\eta_k}\sum_{d\in J(k)}(\frac{1}{\gamma}\|x^k-\bar{x}^{k+1}\|^2+ \gamma\|x^{d} - x^{d+1}\|^2)\cr
\overset{\eqref{eqn:def_bar_x}}=&\textstyle-\frac{1}{2\eta_k^2}\|x^k-\bar{x}^{k+1}\|^2+\frac{|J(k)|}{2\gamma\eta_k}\|x^k-\bar{x}^{k+1}\|^2+\frac{\gamma}{2\eta_k}\sum_{d\in J(k)}\|x^{d} - x^{d+1}\|^2,\nonumber
%\end{aligned}
\end{align}
where the first inequality follows from the Young's inequality. Plugging~\eqref{term2} and~\eqref{term1} into~\eqref{eqn:equality_inconsistent} gives the desired result.\hfill\end{proof}

%\begin{align}
%\quad\|x^{k+1} - x^*\|^2 =& \|x^{k}  - \begin{array}{l}\frac{\eta_k}{mp_{i_k}}\end{array}S_{i_k} \hat{x}^{k}- x^* \|^2  \nonumber\\
% = &\|x^k - x^*\|^2 +  \begin{array}{l}\frac{2\eta_k}{mp_{i_k}}\end{array} \left\langle S_{i_k} \hat{x}^{k}, x^* - x^k \right\rangle + \begin{array}{l}\frac{\eta_k^2}{m^2p_{i_k}^2}\end{array} \|S_{i_k}\hat{x}^{k}\|^2\nonumber\\
% \overset{\eqref{eqn:inconsist}}=& \|x^k - x^*\|^2 + \begin{array}{l}\frac{2\eta_k}{mp_{i_k}}\end{array} \big\langle S_{i_k} \hat{x}^k, x^* - \hat{x}^k + \sum_{d\in J(k)} (x^{d} - x^{d+1})  \big\rangle + \begin{array}{l}\frac{\eta_k^2}{m^2p_{i_k}^2}\end{array} \|S_{i_k}\hat{x}^{k}\|^2\nonumber\\
% = &\|x^k - x^*\|^2 + \underbrace{ \begin{array}{l}\frac{2\eta_k}{mp_{i_k}}\end{array}\big\langle S_{i_k} \hat{x}^k, x^* - \hat{x}^k  \big\rangle }_{A} + \underbrace{\begin{array}{l}\frac{2\eta_k}{mp_{i_k}}\end{array}\sum_{d\in J(k)} \big\langle S_{i_k}\hat{x}^k, x^{d} - x^{d+1} \big\rangle}_{B} + \underbrace{\begin{array}{l}\frac{\eta_k^2}{m^2p_{i_k}^2}\end{array} \|S_{i_k}\hat{x}^{k}\|^2}_{C} \label{eqn:equality_inconsistent}
%\end{align}
%Now we consider the conditional expectation of parts $A$, $B$, and $C$ for given $\cX^k$. For part $A$, we have
%\begin{align*}
%\mathbb{E} (A | \cX^k)         &= \mathbb{E} \Big( 2\eta_k\big\langle \beta_{i_k}S_{i_k} \hat{x}^{k},  x^* - \hat{x}^{k}  \big\rangle | \cX^k\Big)= 2\eta_k\big\langle \mathbb{E} \big(\beta_{i_k}S_{i_k} \hat{x}^{k}| \cX^k \big),  x^* - \hat{x}^{k}  \big\rangle \\
%                        &= \frac{2\eta_k}{m}\big\langle {S} \hat{x}^{k},  x^* - \hat{x}^{k}  \big\rangle \\
%                        &= - \frac{2\eta_k}{m} \big\langle {S} \hat{x}^{k} - Sx^*,   \hat{x}^{k}  - x^* \big\rangle \\%\quad (Sx^*=0)\\
%                        &\leq - \frac{\eta_k}{\alpha  m} \| {S} \hat{x}^{k} - Sx^*\|^2 \quad \text{~( Lemma \ref{lemma:a-avg}) ~}\\
%                        &= -\frac{\eta_k}{\alpha  m} \|S\hat{x}^{k}\|^2 \overset{\eqref{eqn:def_bar_x}}= - \frac{1}{ \alpha  m  \eta_k} \|\bar{x}^{k+1} - x^k\|^2.
%\end{align*}
%The conditional expectation for part $B$ is
%\begin{align*}
%\mathbb{E} (B | \cX^k) &= \frac{2\eta_k}{m} \sum_{d\in J(k)} \big\langle {S}\hat{x}^{k}, x^{d} - x^{d+1} \big\rangle \quad (\text{replace } S_{i_k} \text{ by } S)\\
%                        &=  \frac{2\eta_k}{m} \sum_{d\in J(k)} \big\langle \frac{1}{\eta_k} (x^k - \bar{x}^{k+1}), x^{d} - x^{d+1} \big\rangle \\
%                        &\leq  \frac{1}{ m} \sum_{d\in J(k)} \big( \| x^k - \bar{x}^{k+1} \|^2 + \|{x}^{d} - {x}^{d+1}\|^2 \big) \\
%                        &\leq  \frac{1}{ m} \left(|J(k)| \| x^k - \bar{x}^{k+1} \|^2 + \sum_{d\in J(k)} \|{x}^{d} - {x}^{d+1}\|^2\right) \\
%                        &=  \frac{|J(k)|}{ m}  \| x^k - \bar{x}^{k+1} \|^2 +  \frac{1}{ m} \sum_{d\in J(k)}\|{x}^{d} - {x}^{d+1}\|^2.
%\end{align*}
%The conditional expectation for part $C$ is
%\begin{align*}
%\mathbb{E} (C | \cX^k) = \mathbb{E} \big(\eta_k^2\|\beta_{i_k}S_{i_k}\hat x^k\|^2|\cX^k\big) \leq {\eta_k^2\max_i\beta_i\over m}\|S\hat x^k\|^2 =     {1\over m^2p_{\min}}\| x^k - \bar{x}^{k+1} \|^2.
%\end{align*}
\cut{({\color{red}$\|x^k-\bar x^{k+1}\|^2$ can be replaced by the maximal component: $\max_i(x_i^k-\bar x_i^{k+1})^2$})}
%Substituting the above derived inequalities for part $A$ and $B$ to \eqref{eqn:equality_consistent} gives
%\begin{align}\label{eqn:expect_consistent}
%\mathbb{E} (\|x^{k+1} - x^*\|^2| \cX^k) & = \|x^k - x^*\|^2  + 2\eta_k \big( \mathbb{E} (A | \cX^k) + \mathbb{E} (B | \cX^k) \big) + \eta_k^2 \mathbb{E}(\|S_k(x^{k - j_k})\|^2 | \cX^k) \nonumber \\
 %&= \|x^k - x^*\|^2  + 2\eta_k \big (\mathbb{E} (A | \cX^k) + \mathbb{E} (B | \cX^k) \big)+ \frac{\eta_k^2}{m} \|{S}(x^{k - j_k})\|^2 \nonumber \\
 %&= \|x^k - x^*\|^2  + 2\eta_k \big( \mathbb{E} (A | \cX^k) + \mathbb{E} (B | \cX^k) \big) + \frac{1}{m} \|x^{k} - \bar{x}^{k+1}\|^2 \nonumber \\
 %&\leq \|x^k - x^*\|^2 + \frac{1}{m} (j_k+1 - \frac{1}{\alpha \eta_k}) \|x^k - \bar{x}^{k+1}\|^2 + \frac{1}{m}\sum_{t=0}^{\tau - 1}\|{x}^{k - \tau +t} - {x}^{k - \tau + t+1}\|^2
%\end{align}
%
%Now, we consider the inconsistent read case. As before, we consider the following equality
%\begin{align}\label{eqn:equality_inconsistent}
%\quad\|x^{k+1} - x^*\|^2 &= \|x^{k}  - \eta_k  S_k (\hat{x}^{k})- x^* \|^2 \nonumber \\
%& = \|x^k - x^*\|^2 \underbrace{- 2 \eta_k \big\langle S_k (\hat{x}^{k}), x^k - x^* \big\rangle}_{\text{cross term}} + \eta_k^2 \|S_k(\hat{x}^{k})\|^2.
%\end{align}
%Since we make the following assumption for inconsistent read
%$$x^k=\hat{x}^k+\sum_{d\in J(k)}(x^{d+1}-x^d),$$
%we can split the cross term into two parts as following
%\begin{align*}
%\text{cross term} &= - 2 \eta_k \big\langle S_k (\hat{x}^k), x^k - x^* \big\rangle \\
                %& = {2 \eta_k} \big\langle S_k (\hat{x}^k), x^* - \hat{x}^k + \sum_{d\in J(k)} (x^{d} - x^{d+1})  \big\rangle \\
                %&= 2 \eta_k \Big \{ \underbrace{ \big\langle S_k (\hat{x}^k), x^* - \hat{x}^k  \big\rangle }_{A} + \underbrace{\sum_{d\in J(k)} \big\langle S_k(\hat{x}^k), x^{d} - x^{d+1} \big\rangle}_{B} \Big \}.
%\end{align*}
%Similarly, by considering the conditional expectation for part $A$ given $x^0, ..., x^k$, we have
%\begin{align*}
%\mathbb{E} (A | \cX^k) &\leq  - \frac{1}{2 \alpha m \eta_k^2} \|\bar{x}^{k+1} - x^k\|^2,
%\end{align*}
%and the conditional expectation for part $B$ is
%\begin{align*}
%\mathbb{E} (B | \cX^k) &\leq  \frac{| J(k)|}{2 \eta_k m}  \| x^k - \bar{x}^{k+1} \|^2 +  \frac{1}{2 \eta_k m} \sum_{d\in J(k)}\|{x}^{d} - {x}^{d+1}\|^2.
%\end{align*}

%Substituting the above derived inequalities for parts $A$, $B$, and $C$ to \eqref{eqn:equality_inconsistent} gives
%\begin{align*}%\label{eqn:expect_inconsistent}
%\mathbb{E} (\|x^{k+1} - x^*\|^2| \cX^k) \leq \|x^k - x^*\|^2 + \frac{1}{m} (| J(k)|+{1\over mp_{\min}} - \frac{1}{\alpha \eta_k}) \|x^k - \bar{x}^{k+1}\|^2 + \frac{1}{m}\sum_{d\in J(k)}\|{x}^{d} - {x}^{d+1}\|^2.
%\end{align*}



%Later, we will define a metric for  the sequence $(x^k)_{k\geq0}$ of Algorithm~\ref{alg:asyn_core}.  follows a non-negative almost supermartingale and its convergence result as follows: %. Then we can apply the following convergence theorem for non-negative almost supermartingales.
%Before proving the convergence results, we state the following lemmas.
%The results provide weak convergence for nonexpansive operators, and establish the linear convergence rate for strongly monotone operators. Both convergence results depend strongly on the delay parameter $\tau$. In order to show the weak convergence, we need Lemma~\ref{thm:sup}.
We  need the following lemma on \emph{nonnegative almost supermartingales} \cite{robbins1985convergence}.
\begin{lemma}[{\cite[Theorem 1]{robbins1985convergence}}]\label{thm:sup}
Let $\mathscr{F} = (\cF^k)_{k\geq0}$ be a sequence of sub-sigma algebras of $\cF$ such that $\forall k \ge 0$, $\cF^k \subset \cF^{k+1}$. Define $\ell_+ (\mathscr{F})$ as the set of sequences of $[0, +\infty)$-valued random variables $(\xi_k)_{k\geq 0}$, where $\xi_k$ is $\cF^k$ measurable, and  $\ell_+^1 (\mathscr{F}) := \left\{(\xi_k)_{k\geq 0} \in \ell_+ (\mathscr{F})| \sum_{k} \xi_k < + \infty~ \text{a.s.} \right\}$. Let $(\alpha_k)_{k\geq 0}, (v_k)_{k\geq 0} \in \ell_{+}(\mathscr{F})$, and $(\eta_k)_{k\geq 0}, (\xi_k)_{k\geq 0} \in \ell_{+}^1(\mathscr{F})$ be such that
\begin{equation*}%\label{eqn:sup}
\mathbb{E} (\alpha_{k+1} | \cF^k) + v_k \leq (1 + \xi_k) \alpha_k + \eta_k.
\end{equation*}
Then $(v_k)_{k\ge0} \in \ell_{+}^1 (\mathscr{F})$ and $\alpha_k$ converges to a $[0, +\infty)$-valued random variable a.s..
\end{lemma}

Let $\cH^{\tau+1}=\prod_{i=0}^{\tau}\cH$ be a product space and $\langle\cdot\, |\,\cdot \rangle$ be the induced  inner product:
$$\left\langle (z^0,\ldots,z^{\tau})\,|\,(y^0,\ldots,y^{\tau})\right\rangle=\sum_{i=0}^{\tau}\langle z^i,y^i\rangle,\quad\forall (z^0,\ldots,z^{\tau}), (y^0,\ldots,y^{\tau})\in\cH^{\tau+1}.$$
%Define a $(\tau+1)\times(\tau+1)$ tri-diagonal matrix $M'$ by
%\begin{equation*}%\label{eqn:M_p}
%M':=\begin{bmatrix}1 & 0 & \cdots &0\\
%0 & 0 &\cdots & 0\\ \vdots &\vdots & \ddots & \vdots\\
%0 & 0 &\cdots & 0 \end{bmatrix} +\sqrt{p_{\min}}\begin{bmatrix}
%\tau & -\tau &  & \\
%-\tau & 2\tau-1 & 1-\tau & \\
% & 1-\tau & 2\tau-3 & 2-\tau  & \\
% & & \ddots & \ddots & \ddots &\\
% & & & -2 & 3  & -1 \\
% & & & &-1 & 1
%\end{bmatrix},
%\end{equation*}
Let $M'$ be a symmetric $(\tau+1)\times(\tau+1)$ tri-diagonal matrix with its main diagonal as $\sqrt{p_{\min}}[\frac{1}{\sqrt{p_{\min}}}+\tau,2\tau-1,2\tau-3,\ldots,1]$ and first off-diagonal as $-\sqrt{p_{\min}}[\tau,\tau-1,\ldots,1]$,
and let $M=M'\otimes I_\cH$. Here \cut{$I_\cH$ is the identity operator on $\cH$ and }$\otimes$ represents the Kronecker product. For a given $(y^0,\cdots,y^\tau)\in\cH^{\tau+1}$, $(z^0,\cdots,z^\tau)=M(y^0,\cdots,y^\tau)$ is given by:
\begin{align*}
&z^0=\textstyle y^0+\sqrt{p_{\min}}(y^0-y^1),\\
&z^i = \textstyle\sqrt{p_{\min}}\left((i-\tau-1)y^{i-1}+(2\tau-2i+1)y^i+(i-\tau)y^{i+1}\right),\text{ if } 1\le i\le \tau-1,\\
&z^{\tau}=\textstyle\sqrt{p_{\min}}(y^{\tau}-y^{\tau-1}).
\end{align*}
Then $M$ is a self-adjoint and positive definite linear operator since $M'$ is symmetric and positive definite, and we define $\dotp{\cdot\, |\, \cdot}_M=\dotp{\cdot\, |\, M\cdot}$ as the $M$-weighted inner product and $\|\cdot\|_M$ the induced norm. Let
\begin{equation*}%\label{eqn:def_x_k}
\vx^k=(x^k,x^{k-1},\ldots,x^{k-\tau})\in \cH^{\tau+1},~k\ge 0,\,\mbox{and}~ \vx^* =(x^*,x^*,\ldots,x^*)\in\vX^*\subseteq\cH^{\tau+1},
%\vx^k=\begin{bmatrix}x^k\\ x^{k-1}\\\vdots\\ x^{k-\tau}\end{bmatrix}\in \cH^{\tau+1},~k=0,1,\ldots,\quad\mbox{and}\quad \vx^* = \begin{bmatrix}x^*\\ x^*\\\vdots\\ x^*\end{bmatrix}\in \vX^* \subseteq\cH^{\tau+1},
\end{equation*}
where we set $x^{k}=x^{0}$ for $k<0$. With
\begin{equation}\label{eqn:xi}
\textstyle \xi_k(\vx^*) := \|\vx^k-\vx^*\|_M^2=\|x^{k} - x^*\|^2 + \sqrt{p_{\min}}\sum_{i=k-\tau}^{k-1} \left(i-(k-\tau)+1\right) \| x^{i} - x^{i+1}\|^2,%\quad\forall\,k\ge 0,
\end{equation}
we have the following fundamental inequality:


%Given a self-adjoint and positive definite linear operator $M: \cH^{\tau+1}\to \cH^{\tau+1}$, we let
%$\dotp{\cdot\, |\, \cdot}_M=\dotp{\cdot\, |\, M\cdot}$ be the weighted inner product and $\|\cdot\|_M$ the induced norm. For our purpose, we take $M: (z_1,\ldots,z_{\tau+1})\to (y_1,\ldots,y_{\tau+1})$ such that
%\begin{align*}
%&y_1=z_1+\frac{1}{m}(z_1-z_2),\\
%&y_i = \frac{1}{m}((i-\tau-2)z_{i-1}+(2\tau-2i+3)z_i+(i-\tau-1)z_{i+1}),\text{ if } 2\le i\le \tau,\\
%&y_{\tau+1}=\frac{1}{m}(z_{\tau+1}-z_\tau).
%\end{align*}
%Regarding $(z_1,\ldots,z_{\tau+1})$ as a column vector, we can write $M=M'\otimes I$, where $I$ is the identity operator on $\cH$, $\otimes$ represents the Kronecker product, and
%\begin{equation}\label{eqn:M_p}
%M'=\begin{bmatrix}1 & 0 & \cdots &0\\
%0 & 0 &\cdots & 0\\ \vdots &\vdots & \ddots & \vdots\\
%0 & 0 &\cdots & 0 \end{bmatrix}_{(\tau+1)\times (\tau+1)} +\frac{1}{m}\begin{bmatrix}
%\tau & -\tau &  & \\
%-\tau & 2\tau-1 & 1-\tau & \\
 %& 1-\tau & 2\tau-3 & 2-\tau  & \\
 %& & \ddots & \ddots & \ddots &\\
 %& & & -2 & 3  & -1 \\
 %& & & &-1 & 1
%\end{bmatrix}.
%\end{equation}
%%for $\tau > 0$. When $\tau=0$, we have $M' = [1]$.
%Notice that $M'$ is symmetric positive definite matrix and thus $M$ is a self-adjoint and positive definite operator.  Denote
%\begin{equation}\label{eqn:def_x_k}
%\vx^k=(x^k,x^{k-1},\ldots,x^{k-\tau})\in \cH^{\tau+1},~k\ge 0,\quad\mbox{and}\quad \vx^* =(x^*,x^*,\ldots,x^*)\in\vX^*\subseteq\cH^{\tau+1},
%%\vx^k=\begin{bmatrix}x^k\\ x^{k-1}\\\vdots\\ x^{k-\tau}\end{bmatrix}\in \cH^{\tau+1},~k=0,1,\ldots,\quad\mbox{and}\quad \vx^* = \begin{bmatrix}x^*\\ x^*\\\vdots\\ x^*\end{bmatrix}\in \vX^* \subseteq\cH^{\tau+1},
%\end{equation}
%where we regard $x^{k}=x^{0}$ for $k<0$.
%%\cut{\begin{align*}
%%\xi_{k+1} &:= \|x^{k+1} - x^*\|^2 + \frac{\tau}{n}\|x^k - x^{k+1}\|^2 + {{\tau -1} \over {n}} \|x^{k-1} - x^k\|^2 + \cdots + \frac{1}{n} \|x^{k - \tau +1} - x^{k-\tau +2}\|^2 \\
         %%& = \|x^{k+1} - x^*\|^2 + \sum_{i=1}^{\tau} \frac{i}{n} \| x^{k - \tau + i} - x^{k-\tau +i+1}\|^2.
%%\end{align*}}






%The analysis for the weak convergence follows the next steps. First, we introduce a special metric and show that Algorithm~\ref{alg:asyn_core} generates a sequence $(x^k)_{k\geq 0}$ satisfying~\eqref{eqn:sup} based on the given metric. Then we apply Lemma~\ref{thm:sup} to $(x^k)_{k\geq 0}$ with the special metric and obtain the intermediate results in Theorem~\ref{thm:convergence} (1-2). Finally, we give the weak convergence of $(x^k)_{k\geq 0}$ in Theorem~\ref{thm:convergence} (6).


%\begin{subsection}{Weighted metric and fundamental inequality}\label{sec:weak}




%\cut{This formulation can be written in a more compact form as we now derive. Since
%$$\| x^{k} - x^{k+1}\|^2= \| (x^{k}-x^*) - (x^{k+1}-x^*)\|^2=\|x^k-x^*\|^2 - 2\langle x^k-x^*, x^{k+1}-x^* \rangle+ \|x^{k+1}-x^*\|^2,$$
%we have
%\begin{align*}
%y^{k+1} & = \|x^{k+1} - x^*\|^2 + \sum_{i=1}^{\tau} \frac{i}{n} \left(\| x^{k - \tau + i} -x^*\|^2 -2\langle x^{k - \tau + i} -x^*, x^{k-\tau +i+1}-x^*\rangle+\| x^{k-\tau +i+1}-x^*\|^2\right)\\
%&= \|z^{k+1}-\tilde{x}^*\|_M^2,
%\end{align*}
%where
%$$z^k=\begin{bmatrix}x^k\\ x^{k-1}\\\vdots\\ x^{k-\tau}\end{bmatrix},\quad M=\begin{bmatrix}1 & 0 & \cdots &0\\
%0 & 0 &\cdots & 0\\ & & \cdots & \\
%0 & 0 &\cdots & 0 \end{bmatrix}_{(\tau+1)\times (\tau+1)} +\frac{1}{n}\begin{bmatrix}
%\tau & -\tau &  & \\
%-\tau & 2\tau-1 & 1-\tau & \\
 %& 1-\tau & 2\tau-3 & 2-\tau  & \\
 %& & \ddots & \ddots & \ddots &\\
 %& & & -2 & 3  & -1 \\
 %& & & &-1 & 1
%\end{bmatrix}_{(\tau+1)\times (\tau+1)},$$
%and $\tilde{x}^* = (x^*, x^*, ..., x^*)^T \in {\vX}^*$. Note that $M$ is a positive definite matrix.}
%Let $\tau_k$ be the number of delays at iteration $k$, and
%Before proving the fundamental inequality in Theorem~\ref{thm:fund_inquality}, we d
%Denote $\mathbb{E}\big( \cY \big| \cX \big)$ as the conditional expectation of $\cY$ for given $\cX$ and $\mathbb{E}\big(\cY \big)$ the expectation of $\cY$.
%Throughout the following analysis, we let $p_{\min}=\min_i p_i>0$ and
%\begin{align}\label{eqn:def_bar_x}
%\bar{x}^{k+1} := x^k - \eta_k S\hat{x}^{k}.
%\end{align}
%Then, we have the Lemma~\ref{lemma:fund}.


%\begin{remark}[Non-uniform sampling] For the non-uniform sampling, the conditional expectations for parts A and B are the same and the condition expectation for part C is
%\begin{align*}
%\mathbb{E} (C | \cX^k) = \mathbb{E} \big(\eta_k^2\beta^2_{i_k}\|S_{i_k}\hat x^k\|^2|\cX^k\big) \leq  {\eta_k^2\over m^2 \min_{i}p_i}\|S\hat x^k\|^2 =       {1\over m^2\min_i p_i}\| x^k - \bar{x}^{k+1} \|^2.
%\end{align*}
%\end{remark}


%By Lemma~\ref{lemma:fund}, we show the following theorem.
\begin{theorem}[Fundamental inequality]\label{thm:fund_inquality}
Let $(x^k)_{k\geq 0}$ be the sequence generated by ARock. Then for any $\vx^*\in\vX^*$, it holds that %the following fundamental inequality% for both consistent read and inconsistent read
\begin{equation}\label{eqn:fund_inquality}
\begin{aligned}
\textstyle\mathbb{E}\left(\xi_{k+1}(\vx^*) \,\big|\, \cX^k\right) + \frac{1}{m} \left(\frac{1}{ \eta_k} - \frac{2\tau}{m\sqrt{p_{\min}}} - \frac{1}{mp_{\min}} \right) \|\bar{x}^{k+1} - x^k \|^2
\leq  \xi_k(\vx^*) .
\end{aligned}
\end{equation}
\end{theorem}
\begin{proof}Let $\gamma=m\sqrt{p_{\min}}$.  Since $J(k)\subset\{k-1,\cdots,k-\tau\}$, then \eqref{eqn:fund_inquality0} indicates
\begin{equation}\label{fund-ineq1}
\begin{aligned}
\textstyle\mathbb{E}\big(\|x^{k+1} - x^* \|^2 \,\big|\, \cX^k\big)  \leq & \textstyle\|x^{k} - x^* \|^2  +\frac{1}{\sqrt{p_{\min}}}\sum_{i=k-\tau}^{k-1}\|x^i-x^{i+1}\|^2\\
&\textstyle+ \frac{1}{m}\left(\frac{{\tau}}{m\sqrt{p_{\min}}}+\frac{1}{mp_{\min}}-\frac{1}{\eta_k}\right)\|x^k-\bar x^{k+1}\|^2.
\end{aligned}
\end{equation}
From \eqref{eqn:asyn_update} and \eqref{eqn:def_bar_x}, it is easy to have $\EE(\|x^k-x^{k+1}\|^2|\cX^k)\le\frac{1}{m^2p_{\min}}\|x^k-\bar{x}^{k+1}\|^2$, which together with \eqref{fund-ineq1} implies \eqref{eqn:fund_inquality} by using the definition of $\xi_k(\vx^*)$.
\cut{
\begin{align*}
~&\mathbb{E} (\xi_{k+1}(\vx^*) | \cX^k)  \\
\overset{\eqref{eqn:xi}}= & \textstyle \mathbb{E} (\|x^{k+1} - x^*\|^2| \cX^k) + \gamma\sum_{i=k+1-\tau}^{k} \frac{i-(k-\tau)}{m} \mathbb{E} (\| x^{ i} - x^{i+1}\|^2 | \cX^k) \\
 %=&  \mathbb{E} (\|x^{k+1} - x^*\|^2| \cX^k) + \frac{\tau}{m} \mathbb{E} (\| x^{k } - x^{k+1}\|^2 | \cX^k) + \sum_{i=k+1-\tau}^{k-1} \frac{i-(k-\tau)}{m} \| x^{i} - x^{i+1}\|^2 \\
  \overset{\eqref{eqn:asyn_update}}= & \textstyle\mathbb{E} (\|x^{k+1} - x^*\|^2| \cX^k) + \frac{\gamma\tau}{m} \mathbb{E}(\frac{\eta_k^2}{m^2p_{i_k}^2}\| S_{i_k}\hat x^k\|^2|\cX^k) +  \gamma\sum_{i=k+1-\tau}^{k-1} \frac{i-(k-\tau)}{m} \| x^{i} - x^{i+1}\|^2\\
          \le & \textstyle\mathbb{E} (\|x^{k+1} - x^*\|^2| \cX^k) + \frac{\gamma\tau}{m^3p_{\min}} \| x^{k } - \bar{x}^{k+1}\|^2 +  \gamma\sum_{i=k+1-\tau}^{k-1} \frac{i-(k-\tau)}{m} \| x^{i} - x^{i+1}\|^2\\
  \overset{\eqref{eqn:fund_inquality0}}\leq &\textstyle\|x^k - x^*\|^2 + \frac{1}{m} \left(\frac{|J(k)|}{\gamma} + \frac{\gamma\tau}{m^2p_{\min}} + \frac{1}{mp_{\min}} - \frac{1}{\eta_k}\right) \|x^k - \bar{x}^{k+1}\|^2   \\
&\textstyle+\frac{\gamma}{m}\sum_{d\in J(k)}\|{x}^{d} - {x}^{d+1}\|^2 +  \gamma\sum_{i=k+1-\tau}^{k-1} \frac{i-(k-\tau)}{m} \| x^{i} - x^{i+1}\|^2\\
 \leq &\textstyle\|x^k - x^*\|^2 + \frac{1}{m} \left(\frac{\tau}{\gamma} + \frac{\gamma\tau}{m^2p_{\min}} + \frac{1}{ mp_{\min}} - \frac{1}{ \eta_k}\right) \|x^k - \bar{x}^{k+1}\|^2   \\
&\textstyle+\frac{\gamma}{m}\sum_{i=k-\tau}^{k - 1}\|{x}^{i} - {x}^{i+1}\|^2 +  \gamma\sum_{i=k+1-\tau}^{k-1} \frac{i-(k-\tau)}{m} \| x^{i} - x^{i+1}\|^2\\
%& =  \|x^k - x^*\|^2 + \frac{1}{m} (|J(k)| + \frac{\tau}{m} + {1\over mp_{\min}} - \frac{1}{\alpha \eta_k}) \|x^k - \bar{x}^{k+1}\|^2   \\
%~~~&+ \frac{1}{m} \|x^{k - \tau} - x^{k- \tau +1}\|^2 + \sum_{i=k+1-\tau}^{k-1} \frac{i-(k-\tau)+1}{m} \| x^{i} - x^{i+1}\|^2 \\
%& =  \|x^k - x^*\|^2 + \frac{1}{m} (|J(k)| + \frac{\tau}{m^2p_{\min}} + {1\over mp_{\min}} - \frac{1}{\alpha \eta_k}) \|x^k - \bar{x}^{k+1}\|^2  + \sum_{i=k-\tau}^{k-1} \frac{i-(k-\tau)+1}{m} \| x^{i} - x^{i+1}\|^2 \\
 \overset{\eqref{eqn:xi}}= &\textstyle\xi_k(\vx^*)  + \frac{1}{m} \left(\frac{2\tau}{m\sqrt{p_{\min}}} + \frac{1}{mp_{\min}} - \frac{1}{ \eta_k}\right) \|x^k - \bar{x}^{k+1}\|^2.
\end{align*}
The first inequality follows from the computation of the expectation and~\eqref{term2}, the second inequality holds because $J(k)\subset\{k-1,k-2,\cdots,k-\tau\}$, and the last equality uses $\gamma=m\sqrt{p_{\min}}$, which minimizes $\frac{\tau}{\gamma} + \frac{\gamma\tau}{m^2p_{\min}}$ over $\gamma>0$. Hence, \eqref{eqn:fund_inquality}\cut{the desired inequality} holds.}
\hfill\end{proof}
\begin{remark}[Stochastic Fej\'er monotonicity]\label{rmk:fejer}
From \eqref{eqn:fund_inquality}, if $0<\eta_k \leq \frac{mp_{\min}}{2\tau \sqrt{p_{\min}} + 1}$, then we have $\mathbb{E} (\|\vx^{k+1} - {\vx}^* \|^2_M | \cX^k)  \leq  \|\vx^{k} - {\vx}^* \|^2_M,\,\forall \vx^*\in \vX^*$.\cut{, namely, $(\vx^k)_{k\ge0}$ is stochastic Fej\'er monotone.}
\end{remark}
\begin{remark}\label{rm:stepsize}  Let us check our step size bound  $\frac{mp_{\min}}{2\tau\sqrt{p_{\min}}+1}$.  Consider the uniform case: $p_{\min}\equiv p_i\equiv\frac{1}{m}$. Then, the bound simplifies to $\frac{1}{1+2\tau/\sqrt{m}}$. If the max delay is no more than the square root of coordinates, i.e., $\tau = O(\sqrt{m})$, then the bound is  $O(1)$.  In general,  $\tau$ depends on several factors such as problem structure, system architecture, load balance, etc. If all updates and  agents are identical, then $\tau$ is  proportional to $p$, the number of agents. Hence, ARock takes an $O(1)$ step size for solving a problem with $m$ coordinates by $p=\sqrt{m}$ agents under balanced loads.
\end{remark}

The next lemma is a direct consequence of the invertibility of the metric  $M$. %be proved by observing that operator $M$ is invertible.
\begin{lemma}\label{lem:equiv-weak}
A sequence $(\vz^k)_{k\ge0}\subset\cH^{\tau+1}$ (weakly) converges to $\vz\in\cH^{\tau+1}$ under the metric $\dotp{\cdot\,|\,\cdot}$ if and only if it does so under the metric $\dotp{\cdot\,|\,\cdot}_M$.
\end{lemma}

In light of Lemma~\ref{lem:equiv-weak}, the metric of the inner product for weak convergence in the next lemma is not specified. The lemma and its proof are adapted from~\cite{combettes2014stochastic}. %, weak convergence is equivalent in terms of the original and weighted inner products.
%Based on Lemma~\ref{thm:sup}, Theorem~\ref{thm:fund_inquality}, and Lemma~\ref{lem:equiv-weak}, we show the following almost sure convergence result.% by following~\cite{combettes2014stochastic}.

\begin{lemma}\label{lemma:convergence}
Let $(x^k)_{k\geq0}\subset \cH$ be the sequence generated by ARock with $\eta_k \in [\eta_{\min}, \frac{cmp_{\min}}{2\tau \sqrt{p_{\min}}+1}]$ for any $\eta_{\min}>0$ and $0<c<1$. \cut{$\vx^k$ and  ${\vx}^*$ are defined in~\eqref{eqn:def_x_k},} Then we have:
\begin{enumerate}
\item[\emph{(i)}] $\sum_{k=0}^{\infty}\|x^k - \bar{x}^{k+1}\|^2 < \infty$ a.s..
\item[\emph{(ii)}] $x^k-x^{k+1}\rightarrow 0$ a.s. and $\hat x^k-x^{k+1}\rightarrow 0$ a.s..
\item[\emph{(iii)}] The sequence $(\vx^k)_{k\geq 0}\subset \cH^{\tau+1}$ is bounded a.s..
\item[\emph{(iv)}] There exists $\tilde{\Omega} \in \cF$ such that $P(\tilde{\Omega}) = 1$ and, for every $\omega \in \tilde{\Omega}$ and every ${\vx}^* \in \vX^*$, $(\|\vx^k(\omega) - {\vx}^*\|_M)_{k\geq 0}$ converges.
\item[\emph{(v)}] Let $\mathscr{Z}(\vx^k)$ be the set of weakly convergent cluster points of $(\vx^k)_{k\geq 0}$. Then, $\mathscr{Z}(\vx^k) \subseteq \vX^*$ a.s..
\end{enumerate}
\end{lemma}
\begin{proof}
(i): Note that $\inf_k\left(\frac{1}{ \eta_k} - \frac{2\tau}{m\sqrt{p_{\min}}} - \frac{1}{mp_{\min}}\right)>0$. Also note that, in (\ref{eqn:fund_inquality}), $\|\bar{x}^{k+1} - x^k \|^2= \|\eta_k S\hat{x}^k\|^2$ is $\cX^k$-measurable. Hence, applying Lemma \ref{thm:sup} with $\xi_k=\eta_k=0$ and $\alpha_k=\xi_k(\vx^*),\,\forall k,$ to (\ref{eqn:fund_inquality}) gives this result directly.

(ii) From~(i), we have $x^k-\bar x^{k+1}\rightarrow 0$ a.s.. Since $\|x^k-x^{k+1}\|\leq \frac{1}{mp_{\min}}\|x^k-\bar x^{k+1}\|$, we have $x^k-x^{k+1}\rightarrow 0$ a.s.. Then from~\eqref{eqn:inconsist}, we have $\hat x^k-x^k\rightarrow 0$ a.s..

(iii): From Lemma \ref{thm:sup}, we have that $(\|\vx^{k} - {\vx}^* \|_M^2)_{k\geq 0}$ converges a.s. and so does $(\|\vx^{k} - {\vx}^* \|_M)_{k\geq 0}$, i.e., $\lim_{k\rightarrow \infty} \|\vx^{k} - {\vx}^* \|_M = \gamma$ a.s., where $\gamma$ is a $[0, +\infty )$-valued random variable. Hence, $(\|\vx^{k} - {\vx}^* \|_M)_{k\geq 0}$ must be bounded a.s. and so is $(\vx^k)_{k\geq 0}$.

\cut{(iv): From (ii), we have that for any $\vx^*=(x^*,x^*,\cdots,x^*)\in \vX^*$, there exists $\Omega_{x^*} \in \cF$ such that $P(\Omega_{x^*}) = 1$ and for any $\omega \in \Omega_{x^*}$, $(\|\vx^{k}(\omega) - {\vx}^* \|_M)_{k\geq 0}$ converges. Since $\cH$ is separable, there is a countable set $\cC \subseteq \cH$ such that $\bar{\cC} = \Fix T$, where $\bar{\cC}$ denotes the closure of $\cC$. Let $\tilde{\Omega} = \cap_{u\in \cC} \Omega_{u}$ and $\Omega_u^c$ be the complement of $\Omega_u$. We have
$$\textstyle P(\tilde{\Omega}) = 1 - P(\tilde{\Omega}^c)        = 1 - P(\cup_{u\in \cC} \Omega_{u}^c)
                \geq 1 - \sum_{u\in \cC} P(\Omega_{u}^c)  = 1.$$
%\begin{align*}
%P(\tilde{\Omega}) &= 1 - P(\tilde{\Omega}^c)  \\
%               &= 1 - P(\cup_{x^*\in \cC} \Omega_{x^*}^c) \\
%               &\geq 1 - \sum_{z\in \cC} P(\Omega_{x^*}^c)  = 1.
%\end{align*}
Because $\bar{\cC} =\Fix T$, for a given $x^* \in \Fix T$, there must exist a sequence $(u^l)_{l\ge0} \subset \cC$ such that $\lim_{l\to\infty}u^l = x^*$. In addition, there are random variables $\gamma_l,\,l\ge 0$ such that
\begin{equation}\label{limgamma}
\textstyle \lim_{k\rightarrow \infty} \|\vx^{k}(\omega) - {\vu}^l \|_M = \gamma_l (\omega),\,\forall \omega\in\Omega_{u^l},\,\forall l\ge0,
\end{equation}
where ${\vu}^l = (u^l, u^l, ..., u^l)\in \vX^*$. By the triangle inequality, we have,
$$- \|{\vu}^l - {\vx}^* \|_M  \leq   \|\vx^{k}(\omega) - {\vx}^*\|_M -\|\vx^{k}(\omega) - {\vu}^l \|_M \leq \|{\vu}^l - {\vx}^* \|_M, \forall l, k,\,\forall \omega.$$
Therefore, for any $\omega \in \tilde{\Omega}$ and any $l\ge0$, it holds that
\begin{align*}
- \|{\vu}^l - {\vx}^* \|_M  &\textstyle \leq   \liminf_{k \rightarrow \infty} \big(\|\vx^{k}(\omega) - {\vx}^*\|_M -\|\vx^{k}(\omega) - {\vu}^l \|_M\big) \\
                                         &\textstyle  \overset{\eqref{limgamma}}=   \liminf_{k \rightarrow \infty} \|\vx^{k}(\omega) - {\vx}^*\|_M - \gamma_l(\omega) \\
                                         &\textstyle  \leq   \limsup_{k \rightarrow \infty} \|\vx^{k}(\omega) - {\vx}^*\|_M - \gamma_l(\omega) \\
                                         &\textstyle  = \limsup_{k \rightarrow \infty} \big(\|\vx^{k}(\omega) - {\vx}^*\|_M -  \|\vx^{k}(\omega) - {\vu}^l \|_M\big)\\
                                         & \textstyle \leq \|{\vu}^l - {\vx}^* \|_M.
\end{align*}
Let $l \rightarrow \infty$ in the above inequalites. Noting $\|\vu^l-\vx^*\|_M\to 0$ and $\|\vx^k(\omega)-\vx^*\|_M$ is independent of $l$, we have that $\gamma_\ell(\omega)$ converges to some $\gamma(\omega)$ and thus
$$\textstyle \lim_{k \rightarrow \infty} \|\vx^{k}(\omega) - {\vx}^*\|_M=\gamma(\omega),\,\forall \omega\in\tilde{\Omega}.$$
Since $\vx^*$ is arbitrary and $P(\tilde{\Omega})=1$, this completes the proof of (iv).
%$$\lim_{k \rightarrow \infty} \|\vx^{k}(\omega) - {\vx}^*\|_M = \lim_{l\rightarrow \infty} \lim_{k\rightarrow \infty} \|\vx^{k}(\omega) - {\vu}^l \|_M  =  \lim_{l\rightarrow \infty} \tau_l(\omega).$$
}
(iv): The proof follows directly from \cite[Proposition 2.3 (iii)]{combettes2014stochastic}. It is worth noting that $\tilde \Omega$ in the statement works for all ${\vx}^* \in \vX^*$, namely, $\tilde \Omega$ does not depend on ${\vx}^*$.

(v): By~(ii),  there exists $\hat{\Omega} \in \cF$ such that $P(\hat{\Omega} )=1$ and
\begin{align}\label{lim-xk}x^k(w)-x^{k+1}(w)\rightarrow 0,\quad\forall w\in\hat\Omega.\end{align}
%$$\textstyle \lim_{k\to\infty}\|x^k(\omega) - \bar{x}^{k+1}(\omega)\| \rightarrow 0,\,\forall \omega\in\hat{\Omega}.$$
For any $\omega\in\hat{\Omega}$, let $(\vx^{k_n}(\omega))_{n\ge0}$ be a weakly convergent subsequence of $(\vx^k(\omega))_{k\geq 0}$, i.e., $\vx^{k_n}(\omega) \rightharpoonup \vx$, where $\vx^{k_n}(\omega)  = (x^{k_n}(\omega), x^{k_n - 1}(\omega)..., x^{k_n - \tau}(\omega) )$ and $\vx = (u^0, ..., u^{\tau})$. Note that $\vx^{k_n}(\omega) \rightharpoonup \vx$ implies
$x^{k_n - j}(\omega) \rightharpoonup u^j,\, \forall j.$ %In addition, since
%$\|x^k(\omega) - {x}^{k+1}(\omega) \| \leq \frac{1}{mp_{\min}}\|x^k(\omega) - \bar{x}^{k+1}(\omega) \|,$ we have
%\begin{equation}\label{lim-xk}
%\textstyle \lim_{k\to\infty}\|x^k(\omega) - {x}^{k+1}(\omega)\|= 0.
%\end{equation}
Therefore, $u^i=u^j$, for any $i,j\in\{0,\cdots,\tau\}$ because $x^{k_n - i}(\omega) - x^{k_n - j}(\omega)\rightarrow 0$.



%For any $i, j \in \{0, ..., \tau\}$, we have  $\|u_i - u_j\|  \leq \|x^{k_n - i} - u_i\| + \|x^{k_n - i} - x^{k_n - j}\| + \|x^{k_n - j} - u_j\|.$ Let $n \rightarrow \infty$, we have $\|u_i - u_j\| = 0$.
%Then we have $\vx = (u_0, u_0, ..., u_0)^T$, and we will show $u_0 \in \Fix T$.

Furthermore, observing $\eta_k\ge\eta_{\min}>0$, we have
\begin{equation}\label{lim-xhat}
%\begin{aligned}
\lim_{n\rightarrow \infty} \hat x^{k_n} (\omega) - T \hat x^{k_n}(\omega)=\lim_{n\to\infty}S\hat x^{k_n}(\omega)
=\lim_{n\to\infty}\frac{1}{\eta_{k_n}}(x^{k_n}(\omega)-\bar{x}^{k_n+1}(\omega))=0.
%\end{aligned}
\end{equation}
%$$\forall \omega \in \hat{\Omega}, S\hat x^k(\omega)=\alpha (\hat x^k(\omega) - T \hat x^k(\omega)) \rightarrow 0,$$
%which implies that
%$$\lim_{n\rightarrow \infty} \hat x^{k_n} (\omega) - T \hat x^{k_n}(\omega) =0.$$
From the triangle inequality and the nonexpansiveness of $T$, it follows that
\begin{align*}
&\left\|x^{k_n}(\omega) - Tx^{k_n}(\omega)\right\| \\
=& \left\|x^{k_n}(\omega) - \hat x^{k_n} (\omega) + \hat x^{k_n}(\omega) - T\hat x^{k_n}(\omega) + T\hat x^{k_n}(\omega) - Tx^{k_n}(\omega)\right\| \\
 \leq& \left\|x^{k_n}(\omega) - \hat x^{k_n}(\omega)\right\| + \left\| \hat x^{k_n}(\omega) - T\hat x^{k_n}(\omega)\right\| + \left\| T\hat x^{k_n}(\omega) - Tx^{k_n}(\omega)\right\| \\
 \leq& 2 \left\|x^{k_n}(\omega) - \hat x^{k_n} (\omega)\right\| + \| \hat x^{k_n}(\omega) - T\hat x^{k_n}(\omega)\|\\
 \leq&\textstyle 2\sum_{d\in J(k_n)}\left\|x^{d}(\omega) -  x^{d+1} (\omega)\right\|+ \| \hat x^{k_n}(\omega) - T\hat x^{k_n}(\omega)\|.
\end{align*}
From~\eqref{lim-xk}, \eqref{lim-xhat}, and the above inequality, it follows
$\lim_{n\rightarrow \infty} x^{k_n}(\omega) - T x^{k_n}(\omega) = 0.$
Finally, the demiclosedness principle~\cite[Theorem 4.17]{bauschke2011convex} implies $u^0 \in \Fix T$. %, and this completes the proof of (v).
\hfill\end{proof}


\begin{theorem}\label{thm:convergence}
%Let $(x^k)_{k\geq0}$ be the sequence generated by Algorithm \ref{alg:asyn_core} with $\eta_k \in [\eta_{\min}, \frac{cm^2p_{\min}}{\alpha (\tau m^2p_{\min}+\tau+ m)}]$ for certain $\eta_{\min}>0$ and $0<c<1$. \cut{$\vx^k$ and  ${\vx}^*$ are defined in \eqref{eqn:def_x_k},}
Under the assumptions of Lemma~\ref{lemma:convergence}, the sequence $(\vx^k)_{k\geq 0}$ weakly converges to an $\vX^*$-valued random variable a.s.. In addition, if $T$ is demicompact at 0, $(\vx^k)_{k\geq 0}$ strongly converges to an $\vX^*$-valued random variable a.s.. %During the proof, we also the following statements.
%\begin{enumerate}
%\item[\emph{(i)}] $\sum_{k=0}^{\infty}\|x^k - \bar{x}^{k+1}\|^2 < \infty$ a.s.;
%\item[\emph{(ii)}] $(\vx^k)_{k\geq 0}$ is bounded a.s.;
%\item[\emph{(iii)}] There exists $\tilde{\Omega} \in \cF$ such that $P(\tilde{\Omega}) = 1$ and for every $\omega \in \tilde{\Omega}$ and every ${\vx}^* \in \vX^*$, $(\|\vx^k(\omega) - {\vx}^*\|_M)_{k\geq 0}$ converges.
%\item[\emph{(iv)}] Let $\mathscr{Z}(\vx^k)$ be the set of weakly convergent cluster points of the sequence $(\vx^k)_{k\geq 0}$, then $\mathscr{Z}(\vx^k) \subseteq \vX^*$ a.s..
%\end{enumerate}
\end{theorem}



%Hence, there exists $\hat{\Omega} \in \cF$ such that $P(\hat{\Omega}) = 1$ and $\forall \omega \in \hat{\Omega}$, $\mathscr{Z}(\vx^k(\omega)) \subset \vX^*$.
\begin{proof}
\cut{From Lemma~\ref{lemma:convergence}(i), there exists $\hat{\Omega}\in\cF$ such that $P(\hat{\Omega})=1$ and for every $w\in\hat\Omega$, we have $\sum_{k=0}^\infty \|x^k(w)-\bar x^{k+1}(w)\|^2<\infty$. Lemma~\ref{lemma:convergence}(iii) shows that there exists $\tilde\Omega\in\cF$ such that $P(\tilde\Omega)=1$ and for every $w\in\tilde\Omega$ and $\vx^*\in\vX^*$, $(\|\vx^k(w)-\vx^*\|_M)_{k\geq0}$ converges. Then $P(\tilde{\Omega}\cap\hat{\Omega} ) = 1$. We go to show that for any $\omega \in \tilde{\Omega} \cap \hat{\Omega}$, $\vx^k(\omega)$ weakly converges to a point in $\vX^*$.

For any $\omega \in \tilde{\Omega} \cap \hat{\Omega}$, it follows from Lemma~\ref{lemma:convergence}(iii) that $(\vx^k(\omega))_{k\ge0}$ is a bounded sequence and thus must have a weak cluster point~\cite[Lemma~2.37]{bauschke2011convex}. Assume $\vu(\omega)$ and $\vv(\omega)$ to be two weak cluster points of $(\vx^k(\omega))_{k\ge0}$. From Lemma~\ref{lemma:convergence}(iv), we have $\vu(\omega),\vv(\omega)\in \vX^*$. %Assume that there are two weakly convergent subsequences  converge to two points in $\vX^*$, that is $\vx^{k_n}(\omega) \rightharpoonup \vu(\omega)\in \vx^*$ and $\vx^{l_n}(\omega) \rightharpoonup \vv(\omega)\in\vX^*$. Then (4) implies that $(\|\vx^k(\omega) - \vu(\omega)\|_M)_{k\geq 0}$ and $(\|\vx^k(\omega) - \vv(\omega)\|_M)_{k\geq 0}$ converge. In turn, since
Note
\begin{equation*}
\begin{aligned}
&\big\langle \vx^k(\omega) ~|~ \vu(\omega) - \vv(\omega) \big\rangle_{M} \\
=&\textstyle \frac{1}{2} (\|\vx^k(\omega) - \vv(\omega)\|_M^2 - \|\vx^k(\omega)- \vu(\omega)\|_M^2 + \|\vu(\omega)\|_M^2 - \|\vv(\omega)\|_M^2),
\end{aligned}
\end{equation*}
and from Lemma~\ref{lemma:convergence}(iii), the right hand side of the above equality converges. Thus $\langle \vx^k(\omega) ~|~ \vu(\omega) - \vv(\omega) \rangle_{M}$ converges. Since $\vu(\omega)$ and $\vv(\omega)$ are cluster points of $(\vx^k(\omega))_{k\ge0}$, we have
$$\langle \vu(\omega) ~|~ \vu(\omega) - \vv(\omega) \rangle_{M} = \langle \vv(\omega) ~|~ \vu(\omega) - \vv(\omega) \rangle_{M},$$
which indicates
%the sequence $\big(\langle \vx^k(\omega) ~|~ \vu(\omega) - \vv(\omega) \rangle_{M} \big)_{k\geq 0}$ converges, say
%$$\langle \vx^k(\omega) ~|~ \vu(\omega) - \vv(\omega) \rangle_{M} \rightarrow \rho(\omega).$$
%Then $\langle \vu(\omega) ~|~ \vu(\omega) - \vv(\omega) \rangle_{M} = \langle \vu(\omega) ~|~ \vu(\omega) - \vv(\omega) \rangle_{M} =\rho(\omega)$, because $\vx^{k_n}(\omega) \rightharpoonup \vu(\omega)$ and $\vx^{l_n}(\omega) \rightharpoonup \vu(\omega)$.
%However, since $\vx^{k_n}(\omega) \rightharpoonup \vu(\omega)$, we have $\langle \vu(\omega) ~|~ \vu(\omega) - \vv(\omega) \rangle_{M} = \rho(\omega),$ similarly, we also have $\langle \vv(\omega) ~|~ \vu(\omega) - \vv(\omega) \rangle_{M} = \rho(\omega)$.
%Thus
%$$0 = \big\langle \vu(\omega) ~|~ \vu(\omega) - \vv(\omega) \big\rangle_{M} - \big\langle \vv(\omega) ~|~ \vu(\omega) - \vv(\omega) \big\rangle_{M}  = \|\vu(\omega) - \vv(\omega)\|_M^2,$$
%which implies that
$\vu(\omega) = \vv(\omega)$. Therefore, $(\vx^k(\omega))_{k\ge0}$ has a unique weak cluster point. This completes the proof of almost sure weak convergence.}
%The proof for a.s. weak convergence follows directly from~\cite[Proposition~2.3~(iv)]{combettes2014stochastic}.
{The proof for a.s. weak convergence follows from  Opial's Lemma~\cite{peypouquet_evolution_2009,opial_weak_1967} and Lemma~\ref{lemma:convergence} (iv)-(v).}
Next we assume that $T$ is demicompact at 0. From the proof of Lemma~\ref{lemma:convergence} (v), there is $\hat\Omega\in \cF$ such that $P(\hat\Omega)=1$ and, for any $w\in\hat\Omega$ and any weakly convergent subsequence of $(\vx^{k_n}(w))_{n\ge0}$, $\lim_{n\rightarrow \infty}x^{k_n}(w)-Tx^{k_n}(w)=0$. Since $T$ is demicompact, $(x^{k_n}(w))_{n\ge 0}$ has a strongly convergent subsequence, for which we still use $(x^{k_n}(w))_{n\ge 0}$. Hence, $x^{k_n}(w)\rightarrow \bar x(w)\in \Fix T$.   Lemma~\ref{lemma:convergence} (ii) yields $\vx^{k_n}(w)\rightarrow \bar\vx(w)\in \vX^*$. Then by Lemma~\ref{lemma:convergence} (iv), there is $\tilde\Omega\in\cF$ such that $P(\tilde\Omega)=1$ and, for every $w\in\tilde\Omega$ and every $\vx^*\in\vX^*$, $(\|\vx^k(w)-\vx^*\|_M)_{k\geq 0}$ converges. Thus, for any $w\in\hat\Omega\cap\tilde\Omega$, we have $\lim_{k\to\infty}\|\vx^k(w)-\bar\vx(w)\|_M=0$. Because $P(\hat\Omega\cap\tilde\Omega)=1$, we conclude that $(\vx^k)_{k\geq0}$ strongly converges to an $\vX^*$-valued random variable a.s..
\hfill\end{proof}

\cut{Next we show the almost sure strong convergence for uniformly quasi-monotone operators $S$. If $S$ is uniformly quasi-monotone, $\Fix T$ is a singleton. Indeed, let $x,y\in \Fix T$. Since $0\in Sx$ and $0\in Sy$, we have $0=\langle 0-0,x-y\rangle\geq \phi(\|x-y\|)$ and thus $\|x-y\|=0$ or $x=y$.

For any $w\in \hat\Omega\cup \tilde\Omega$, we have $x^k(w)-\bar x^{k+1}(w)\rightarrow 0$ because of Lemma~\ref{lemma:convergence}(i). Therefore we have $\hat x^k(w)-x^k(w)\rightarrow 0$ because of~\eqref{eqn:inconsist}. From the almost sure weak convergence, we have $x^k(w)\rightharpoonup  x^*\in \Fix T$. Thus $\hat x^k(w)\rightharpoonup x^*$ and
\begin{align*}
\textstyle\phi(\|\hat x^k(w)-x^*\|)\leq& \langle S\hat x^k(w)-Sx^*,\hat x^k(w)-x^*\rangle\\
=&\frac{1}{\eta_k}\langle x^k(w)-\bar x^{k+1}(w), \hat x^k(w)-x^*\rangle \rightarrow 0,
\end{align*}
or $\hat x^k(w)\rightarrow x^*$. Combining with $\hat x^k(w)-x^k(w)\rightarrow 0$, we have $x^k(w)\rightarrow x^*$.
%Since $P(\hat{\Omega} \cap \tilde{\Omega}) = 1$, we have $\vx^k \rightharpoonup {\vx}^*$ a.s..
\hfill\end{proof}}

\begin{remark} For the generalization in Section~\ref{sec:general}, we need to replace \eqref{term2} by
\begin{align*}
\textstyle\sum_{i=1}^m\frac{1}{p_i} \|U_i\circ S\hat{x}^{k}\|^2\le\frac{1}{p_{\min}} \sum_{i=1}^m\|U_i\circ S\hat{x}^{k}\|^2\leq \frac{C}{p_{\min}} \|S\hat{x}^{k}\|^2=\frac{C}{\eta_k^2p_{\min}}\|x^k-\bar{x}^{k+1}\|^2,
\end{align*}
%instead of~\eqref{term2}.
and update the step size condition to $\eta_k \in [\eta_{\min}, \frac{cmp_{\min}}{2\tau \sqrt{p_{\min}}+C}]$. Then the proofs of Theorem~\ref{thm:convergence} and Lemma~\ref{lemma:convergence} will go through and yield the same convergence result. % except the step size has to satisfy the condition $\eta_k \in [\eta_{\min}, \frac{cmp_{\min}}{2\tau \sqrt{p_{\min}}+C}]$.
\end{remark}



%\end{subsection}

% !TEX root = ./asyn_notes.tex

\subsection{Linear convergence\cut{ for quasi-strongly monotone operator $S$}}\label{sec:linear}
In this section, we establish  linear convergence  under the assumption that $S$ is quasi-strongly monotone.
%\begin{assumption}[Quasi $\mu$-strongly monotone]\label{assumption:strongly_monotone}
%Assume $S: \cH \rightarrow \cH$ to be quasi $\mu$-strongly monotone with $\mu >0$, i.e., $\forall x\in \cH$, and $\forall y\in \mathrm{zer}S$, we have
%\begin{equation}\label{eq:stgmon}
%\langle S(x), x-y\rangle\ge \mu\|x-y\|^2.
%\end{equation}
%\end{assumption}
We first present a key lemma.

\begin{lemma}\label{lemma:linear_induction}
Assume that the step size is fixed, i.e., $\eta_k  = \eta$, and satisfies
\begin{equation}\label{eta-cond1}
\textstyle 0< \eta \leq \underline{\eta}_1:=(1 - \frac{1}{\rho}) \frac{m \sqrt{p_{\min}}}{8} \frac{\rho^{1/2} - 1}{ \rho^{(\tau + 1)/2} - 1}
\end{equation}
for some $\rho>1$. Then we have, for all $k\geq 1$,
\begin{equation} \label{eqn:linear_rate}
\mathbb{E} \|\bar{x}^{k} - x^{k-1}\|^2 \leq \rho \mathbb{E} \|\bar{x}^{k+1} - x^k\|^2.
\end{equation}
\end{lemma}

\cut{
\begin{proof}
We prove \eqref{eqn:linear_rate} by induction. First, for any $k \geq 1$, we observe that
%\begin{equation}
\begin{align}
& \|\bar{x}^k - x^{k-1}\|^2 - \|\bar{x}^{k+1} - x^k\|^2 \cr
\leq \, &2 \| \bar{x}^k - x^{k-1} \| \| \bar{x}^{k+1} - x^k -  \bar{x}^k + x^{k-1}\| \quad \text{(by $\|a\|^2 - \|b\|^2 \leq 2 \|a\| \|b - a\|$)}\cr
=\, & 2 \| \bar{x}^k - x^{k-1} \| \| \eta S(\hat x^k) -  \eta S(\hat x^{k -1})\| \cr
\textstyle \leq \, & 4 \eta \| \bar{x}^k - x^{k-1} \| \| \hat x^{k} - \hat x^{k -1}\| \label{eq-mid}\\
 \leq \, & 4 \eta \| \bar{x}^k - x^{k-1} \| \big ( \| x^k  - \hat x^{k}\| + \|x^k - x^{k-1}\| + \| x^{k-1}- \hat x^{k -1}\| \big) \cr
= \, & \textstyle 4 \eta \| \bar{x}^k - x^{k-1} \| \left ( \| \sum_{d \in J(k)}(x^d - x^{d+1})\| + \|x^k - x^{k-1}\|\right.\cr
& \textstyle \hspace{3cm}+ \left.\|\sum_{d \in J(k-1)}(x^d - x^{d+1})\| \right) \cr
 \leq \, & \textstyle 4 \eta \| \bar{x}^k - x^{k-1} \| \left ( \sum_{d \in J(k)} \| x^d - x^{d+1}\| + \|x^k - x^{k-1}\| \right.\cr
 & \textstyle \hspace{3cm}+ \left.\sum_{d \in J(k-1)}\|x^d - x^{d+1}\| \right) \cr
\leq \, & \textstyle 4 \eta \| \bar{x}^k - x^{k-1} \|  \big( 2 \sum_{t=0}^{\tau } \| x^{k-t}  - x^{k - t -1}\| \big)\cr
= \, & \textstyle 8 \eta \sum_{t=0}^{\tau }  \| \bar{x}^k - x^{k-1} \| \| x^{k-t}  - x^{k - t -1}\| \label{eqn:linear_convergence_basic_step}
\end{align}
%\end{equation}

For the basic case,  we know that $\hat x^0 = x^0$, $\hat x^1 \in \{x^0, x^1\}$. Letting $k=1$ in \eqref{eq-mid} gives
\begin{equation}\label{eqn:linear_basic_case}
\begin{aligned}
       & \mathbb{E} \|\bar{x}^1 - x^{0}\|^2 -  \mathbb{E} \|\bar{x}^{2} - x^1\|^2 \\
\leq \, &  4\eta \mathbb{E}  \| \bar{x}^1 - x^{0} \|   \| x^1  - x^{0}\| \\
\leq \, &  \textstyle 2 \eta \big(\frac{1}{m \sqrt{p_{\min}}} \mathbb{E}  \| \bar{x}^1 - x^{0} \|^2  + m \sqrt{p_{\min}} \mathbb{E} \| x^1  - x^{0}\|^2 \big)\\
= \, &\textstyle  2 \eta \big(\frac{1}{m \sqrt{p_{\min}}} \mathbb{E}  \| \bar{x}^1 - x^{0} \|^2  + m \sqrt{p_{\min}} \sum_{i = 1}^m p_{i} \frac{\eta^2}{m^2 p_{i}^2} (S_i x^0)^2  \big)\\
\leq \, & \textstyle 2 \eta \big(\frac{1}{m \sqrt{p_{\min}}} \mathbb{E}  \| \bar{x}^1 - x^{0} \|^2  + \frac{1}{m \sqrt{p_{\min}}} \mathbb{E}  \| \bar{x}^1 - x^{0} \|^2  \big) \\
= \, & \textstyle \frac{4 \eta} {m \sqrt{p_{\min}}}  \mathbb{E}  \| \bar{x}^1 - x^{0} \|^2.
\end{aligned}
\end{equation}
Rearranging \eqref{eqn:linear_basic_case} gives $\mathbb{E} \|\bar{x}^1 - x^{0}\|^2 \leq \frac{1}{1 - \frac{4 \eta} {m \sqrt{p_{\min}}}} \mathbb{E} \| \bar{x}^2  - x^{1}\|^2.$ From \eqref{eta-cond1} and $\rho>1$, it holds that
$$0 < \eta \leq (1 - \frac{1}{\rho}) \frac{m \sqrt{p_{\min}}}{8} \frac{\rho^{1 / 2} - 1}{ \rho^{(\tau + 1) 2} - 1} \leq (1 - \frac{1}{\rho}) \frac{m\sqrt{p_{\min}}}{4}.$$ Hence, $\mathbb{E} \|\bar{x}^1 - x^{0}\|^2 \leq \rho \mathbb{E} \| \bar{x}^2  - x^{1}\|^2$.

For the induction step, we have
\begin{align*}
&\mathbb{E} \|\bar{x}^k - x^{k-1}\| \|x^{k-t} - x^{k-t-1}\| \\
\leq \, & \textstyle \frac{1}{2} \mathbb{E} \big\{ a \|x^{k-t} -x^{k-t-1}\|^2 + \frac{1}{a} \|\bar{x}^{k} - x^{k-1}\|^2 \big\} \\
\leq \, &  \textstyle \frac{1}{2} \mathbb{E} \big\{ \frac{a}{m^2 p_{\min}} \|\bar{x}^{k-t} -x^{k-t-1}\|^2 + \frac{1}{a} \|\bar{x}^{k} - x^{k-1}\|^2 \big\} \\
\leq \, &\textstyle \frac{1}{2}  \big\{ \frac{a \rho^t}{m^2 p_{\min}}  + \frac{1}{a} \big\} \mathbb{E} \|\bar{x}^{k} - x^{k-1}\|^2 \\
= \, & \textstyle \frac{ {\rho}^{t/2} }{m\sqrt{p_{\min}}}   \mathbb{E} \|\bar{x}^{k} - x^{k-1}\|^2. \qquad (\text{by letting } a = m \sqrt{p_{\min}}\rho^{ - { t/2}})\\
\end{align*}
Taking the expectation on both sides of the inequality \eqref{eqn:linear_convergence_basic_step} gives
\begin{align*}
& \textstyle \mathbb{E}\|\bar{x}^k - x^{k-1}\|^2 - \mathbb{E} \|\bar{x}^{k+1} - x^k\|^2 \\
\leq & \textstyle 8 \eta \sum_{t=0}^{\tau }  \mathbb{E} \| \bar{x}^k - x^{k-1} \| \| x^{k-t}  - x^{k - t -1}\| \\
\leq & \textstyle \frac{8 \eta}{m \sqrt{p_{\min}}} \sum_{t=0}^{\tau } \rho^{t / 2} \mathbb{E} \| \bar{x}^k - x^{k-1} \|^2 \\
\leq & \textstyle \frac{8 \eta}{m \sqrt{p_{\min}}} \frac{1 - \rho^{(\tau + 1 )/2}}{1 - \rho^{1 / 2}}\mathbb{E}  \| \bar{x}^k - x^{k-1} \|^2, %\\
%= &\textstyle \eta \theta \mathbb{E} \| \bar{x}^k - x^{k-1} \|^2,
\end{align*}
%where $\theta = \frac{8}{m \sqrt{p_{\min}}} \frac{1 - \rho^{\tau + 1 \over 2}}{1 - \rho^{1 \over 2}}$.
Rearranging the above inequality and using \eqref{eta-cond1} gives
 $$\mathbb{E} \|\bar{x}^{k} - x^{k-1}\|^2 \cut{\leq \frac{1}{1 - \eta \theta} \mathbb{E} \| \bar{x}^{k+1}  - x^{k}\|^2 }\leq \rho \mathbb{E} \| \bar{x}^{k+1}  - x^{k}\|^2. $$
 \hfill
\end{proof}
}

\begin{proof}
We prove \eqref{eqn:linear_rate} by induction. First, based on the inequality $\|a\|^2 - \|b\|^2 \leq 2 \|a\| \|b - a\|$ we observe that, for any $k \geq 1$,
\begin{align}\label{eq-mid}
\|\bar{x}^k - x^{k-1}\|^2 - \|\bar{x}^{k+1} - x^k\|^2 \leq \, &2 \| \bar{x}^k - x^{k-1} \| \| \bar{x}^{k+1} - x^k -  \bar{x}^k + x^{k-1}\| \cr
=\, & 2 \| \bar{x}^k - x^{k-1} \| \| \eta S(\hat x^k) -  \eta S(\hat x^{k -1})\| \cr
\textstyle \leq \, & 4 \eta \| \bar{x}^k - x^{k-1} \| \| \hat x^{k} - \hat x^{k -1}\|.
\end{align}
Applying the triangle inequality and \eqref{eqn:inconsist} yields
\begin{align}\label{eqn:linear_convergence_basic_step}
\| \hat x^{k} - \hat x^{k -1}\| \leq & \| x^k  - \hat x^{k}\| + \|x^k - x^{k-1}\| + \| x^{k-1}- \hat x^{k -1}\| \cr
 \leq \, &  \textstyle \underset{d \in J(k)}\sum \| x^d - x^{d+1}\| + \|x^k - x^{k-1}\| + \underset{d \in J(k-1)}\sum\|x^d - x^{d+1}\| \cr
\leq \, & \textstyle  2 \sum_{t=0}^{\tau } \| x^{k-t}  - x^{k - t -1}\|.
\end{align}
%\end{equation}
For the basic case,  we have $\hat x^0 = x^0$, $\hat x^1 \in \{x^0, x^1\}$. Letting $k=1$ in \eqref{eq-mid} gets us
\begin{align*}%\label{eqn:linear_basic_case}
%\begin{aligned}
        \mathbb{E} \|\bar{x}^1 - x^{0}\|^2 -  \mathbb{E} \|\bar{x}^{2} - x^1\|^2
\leq \, &  4\eta \mathbb{E}  \| \bar{x}^1 - x^{0} \|   \| x^1  - x^{0}\| \\
\leq \, &  \textstyle 2 \eta \big(\frac{1}{m \sqrt{p_{\min}}} \mathbb{E}  \| \bar{x}^1 - x^{0} \|^2  + m \sqrt{p_{\min}} \mathbb{E} \| x^1  - x^{0}\|^2 \big)\\
= \, &\textstyle  2 \eta \big(\frac{1}{m \sqrt{p_{\min}}} \mathbb{E}  \| \bar{x}^1 - x^{0} \|^2  + m \sqrt{p_{\min}} \sum_{i = 1}^m p_{i} \frac{\eta^2}{m^2 p_{i}^2} (S_i x^0)^2  \big)\\
\leq \, & \textstyle 2 \eta \big(\frac{1}{m \sqrt{p_{\min}}} \mathbb{E}  \| \bar{x}^1 - x^{0} \|^2  + \frac{1}{m \sqrt{p_{\min}}} \mathbb{E}  \| \bar{x}^1 - x^{0} \|^2  \big) \\
= \, & \textstyle \frac{4 \eta} {m \sqrt{p_{\min}}}  \mathbb{E}  \| \bar{x}^1 - x^{0} \|^2.
%\end{aligned}
\end{align*}
Rearranging the above inequality yields $\mathbb{E} \|\bar{x}^1 - x^{0}\|^2 \leq \frac{1}{1 - \frac{4 \eta} {m \sqrt{p_{\min}}}} \mathbb{E} \| \bar{x}^2  - x^{1}\|^2.$ By \eqref{eta-cond1} and $\rho>1$, it holds that
$\textstyle 0 < \eta \leq (1 - \frac{1}{\rho}) \frac{m \sqrt{p_{\min}}}{8} \frac{\rho^{1/2} - 1}{ \rho^{(\tau + 1)/2} - 1} \leq (1 - \frac{1}{\rho}) \frac{m\sqrt{p_{\min}}}{4}.$ Hence, $\mathbb{E} \|\bar{x}^1 - x^{0}\|^2 \leq \rho \mathbb{E} \| \bar{x}^2  - x^{1}\|^2$.

For the induction step, applying Young's inequality gives us
\begin{align*}
\mathbb{E} \|\bar{x}^k - x^{k-1}\| \|x^{k-t} - x^{k-t-1}\|
\leq \, & \textstyle \frac{1}{2} \mathbb{E} \big\{ a \|x^{k-t} -x^{k-t-1}\|^2 + \frac{1}{a} \|\bar{x}^{k} - x^{k-1}\|^2 \big\} \\
\leq \, &  \textstyle \frac{1}{2} \mathbb{E} \big\{ \frac{a}{m^2 p_{\min}} \|\bar{x}^{k-t} -x^{k-t-1}\|^2 + \frac{1}{a} \|\bar{x}^{k} - x^{k-1}\|^2 \big\} \\
\leq \, &\textstyle \frac{1}{2}  \big\{ \frac{a \rho^t}{m^2 p_{\min}}  + \frac{1}{a} \big\} \mathbb{E} \|\bar{x}^{k} - x^{k-1}\|^2 \\
= \, & \textstyle \frac{ {\rho}^{t/2} }{m\sqrt{p_{\min}}}   \mathbb{E} \|\bar{x}^{k} - x^{k-1}\|^2. \qquad (\text{letting } a = m \sqrt{p_{\min}}\rho^{ - { t/ 2}})
\end{align*}
Taking the expectation on \eqref{eqn:linear_convergence_basic_step} and combining it with \eqref{eq-mid} yield
\begin{align*}
& \textstyle \mathbb{E}\|\bar{x}^k - x^{k-1}\|^2 - \mathbb{E} \|\bar{x}^{k+1} - x^k\|^2
\leq  \textstyle 8 \eta \sum_{t=0}^{\tau }  \mathbb{E} \| \bar{x}^k - x^{k-1} \| \| x^{k-t}  - x^{k - t -1}\| \\
\leq & \textstyle \frac{8 \eta}{m \sqrt{p_{\min}}} \sum_{t=0}^{\tau } \rho^{t /2} \mathbb{E} \| \bar{x}^k - x^{k-1} \|^2
\leq  \textstyle \frac{8 \eta}{m \sqrt{p_{\min}}} \frac{1 - \rho^{(\tau + 1 )/2}}{1 - \rho^{1 /2}}\mathbb{E}  \| \bar{x}^k - x^{k-1} \|^2. %\\
%= &\textstyle \eta \theta \mathbb{E} \| \bar{x}^k - x^{k-1} \|^2,
\end{align*}
%where $\theta = \frac{8}{m \sqrt{p_{\min}}} \frac{1 - \rho^{\tau + 1 \over 2}}{1 - \rho^{1 \over 2}}$.
Finally, rearranging the above inequality and using \eqref{eta-cond1} lead to
 $\mathbb{E} \|\bar{x}^{k} - x^{k-1}\|^2 \cut{\leq \frac{1}{1 - \eta \theta} \mathbb{E} \| \bar{x}^{k+1}  - x^{k}\|^2 }\leq \rho \mathbb{E} \| \bar{x}^{k+1}  - x^{k}\|^2. $
 This completes the proof.
 \end{proof}

With this lemma, we are ready to derive the linear convergence rate of ARock.

%%%%%%%%%%%%%%%%% Theorem %%%%%%%%%%%%%%%%%%%%
\begin{theorem}[Linear convergence\cut{ for quasi-$\mu$-strongly monotone operator}]\label{thm:strongly_monotone}
Assume that $S$ is quasi-$\mu$-strongly monotone with $\mu >0$. Let $\beta\in (0,1)$ and $(x^k)_{k\geq0}$ be the sequence generated by ARock with a constant stepsize $\eta \in (0, \min \{\underline{\eta}_1, \underline{\eta}_2\}]$, where $\underline{\eta}_1$ is given in~\eqref{eta-cond1} and
\begin{equation}\label{eta-cond2}\textstyle
\underline{\eta}_2=\frac{-b + \sqrt{b^2 + 4(1-\beta)a}}{2a},\quad a = \frac{2\beta \mu \tau }{ m^2 p_{\min}} \frac{\rho(\rho^\tau - 1)}{\rho - 1},\quad b = \frac{1}{ m p_{\min}} + \frac{2}{ m}\sqrt{ \frac{\rho(\rho^\tau - 1)\tau}{(\rho - 1) p_{\min}}}.
\end{equation}
Then
\begin{equation}\label{linear-convg}
\textstyle\mathbb{E}\left(\|x^{k} - x^* \|^2\right)\le \left(1-\frac{\beta\mu\eta}{m}\right)^k \|x^{0} - x^* \|^2.
\end{equation}


\end{theorem}
\begin{proof}
Following the proof of Lemma \ref{lemma:fund} and starting from \eqref{term1}, we have
\begin{equation*}\label{thm:linear_cvg}
\begin{aligned}
&\langle S \hat{x}^{k}, x^* - x^k \rangle\\
%=&\textstyle\langle S \hat{x}^{k},x^* - \hat{x}^k + \sum_{d\in J(k)} (x^{d} - x^{d+1})\rangle\cr
%=&\textstyle\langle S \hat{x}^{k}, x^* - \hat{x}^k\rangle + \frac{1}{\eta}\sum_{d\in J(k)}\langle x^k-\bar{x}^{k+1}, x^{d} - x^{d+1}\rangle\cr
\le&\textstyle \langle S \hat{x}^{k}-Sx^*, x^* - \hat{x}^k\rangle+\frac{1}{2\eta}\sum_{d\in J(k)}\big(\frac{1}{\gamma}\|x^k-\bar{x}^{k+1}\|^2+ \gamma\|x^{d} - x^{d+1}\|^2\big)\cr
\le&\textstyle -{\beta\mu}\|\hat{x}^k-x^*\|^2-\frac{1-\beta}{2}\|S \hat{x}^{k}\|^2+\frac{1}{2\eta}\sum_{d\in J(k)}(\frac{1}{\gamma}\|x^k-\bar{x}^{k+1}\|^2+ \gamma\|x^{d} - x^{d+1}\|^2)\cr
=&\textstyle -{\beta\mu}\|{x}^k-x^*+\sum_{d\in J(k)} (x^{d} - x^{d+1})\|^2-\frac{1-\beta}{2\eta^2}\|x^k-\bar{x}^{k+1}\|^2 \\
 &\textstyle +\frac{|J(k)|}{2\gamma\eta}\|x^k-\bar{x}^{k+1}\|^2+\frac{\gamma}{2\eta}\sum_{d\in J(k)}\|x^{d} - x^{d+1}\|^2\cr
\le&\textstyle -\frac{\beta\mu}{2}\|{x}^k-x^*\|^2+{\beta\mu}\|\sum_{d\in J(k)} (x^{d} - x^{d+1})\|^2-\frac{1-\beta}{2\eta^2}\|x^k-\bar{x}^{k+1}\|^2 \\
& \textstyle+\frac{|J(k)|}{2\gamma\eta}\|x^k-\bar{x}^{k+1}\|^2+\frac{\gamma}{2\eta}\sum_{d\in J(k)}\|x^{d} - x^{d+1}\|^2\cr
\le & \textstyle -\frac{\beta\mu}{2}\|{x}^k-x^*\|^2+{\beta\mu}|J(k)|\sum_{d\in J(k)}\| x^{d} - x^{d+1}\|^2-\frac{1-\beta}{2\eta^2}\|x^k-\bar{x}^{k+1}\|^2 \\
&\textstyle +\frac{|J(k)|}{2\gamma\eta}\|x^k-\bar{x}^{k+1}\|^2+\frac{\gamma}{2\eta}\sum_{d\in J(k)}\|x^{d} - x^{d+1}\|^2,
\end{aligned}
\end{equation*}
where the second inequality holds because $S$ is $\frac{1}{2}$-cocoercive and also quasi-$\mu$-strongly monotone, and the last one comes from the Cauchy-Schwartz inequality.
Plugging the above inequality and \eqref{term2} into \eqref{eqn:equality_inconsistent} and noting $|J(k)|\subset\{k-\tau,\ldots,k-1\}$ gives
\begin{align*}
\begin{aligned}
\textstyle\mathbb{E}\big(\|x^{k+1} - x^* \|^2 \,\big|\, \cX^k\big)  \leq & \textstyle\big(1-\frac{\beta\mu\eta}{m}\big)\|x^{k} - x^* \|^2  +\frac{1}{ m}\left(2\beta\eta\mu\tau+\gamma\right)\sum_{d=k-\tau}^{k-1}\|x^d-x^{d+1}\|^2\\
&\textstyle+ \frac{1}{m}\left(\frac{{\tau}}{\gamma}+\frac{1}{ mp_{\min}}-\frac{1-\beta}{{\eta}}\right)\|x^k-\bar x^{k+1}\|^2.
\end{aligned}
\end{align*}
Taking expectation over both sides of the above inequality, noting $\EE\|x^d-x^{d+1}\|^2\le \frac{1}{m^2p_{\min}}\EE\|x^d-\bar{x}^{d+1}\|^2$, and using Lemma \ref{lemma:linear_induction}, we have
\begin{align*}
\begin{aligned}
~ &\mathbb{E}\big(\|x^{k+1} - x^* \|^2\big)  \\
%\leq & \textstyle \big(1-\frac{\beta\mu\eta}{m}\big)\EE\|x^{k} - x^* \|^2  +\frac{{1}}{ m^3 p_{\min}}\left(2\beta\eta\mu\tau+{\gamma}\right)\sum_{d\in J(k)}\EE\|x^d-\bar{x}^{d+1}\|^2\\
%&\textstyle + \frac{1}{m}\left(\frac{{\tau}}{\gamma}+\frac{1}{mp_{\min}}-\frac{1-\beta}{{\eta}}\right)\EE\|x^k-\bar x^{k+1}\|^2\cr
%\le&\textstyle\big(1-\frac{\beta\mu\eta}{m}\big)\EE\|x^{k} - x^* \|^2  +\frac{{1}}{m^3 p_{\min}}\left(2\beta\eta\mu\tau+{\gamma}\right)\sum_{d= k-\tau}^{k-1}\EE\|x^d-\bar{x}^{d+1}\|^2\\
%&\textstyle + \frac{1}{m}\left(\frac{{\tau}}{\gamma}+\frac{1}{mp_{\min}}-\frac{1-\beta}{{\eta}}\right)\EE\|x^k-\bar x^{k+1}\|^2\cr
\le& \textstyle \big(1-\frac{\beta\mu\eta}{m}\big)\EE\|x^{k} - x^* \|^2  +\frac{{1}}{ m^3 p_{\min}}\left(2\beta\eta\mu\tau+{\gamma}\right)\sum_{d= 1}^{\tau}\rho^d\EE\|x^k-\bar{x}^{k+1}\|^2\\
&\textstyle + \frac{1}{m}\left(\frac{{\tau}}{\gamma}+\frac{1}{mp_{\min}}-\frac{1-\beta}{{\eta}}\right)\EE\|x^k-\bar x^{k+1}\|^2\cr
=& \textstyle \big(1-\frac{\beta\mu\eta}{m}\big)\EE\|x^{k} - x^* \|^2  +\frac{{1}}{m^3 p_{\min}}\left(2\beta\eta\mu\tau+{\gamma}\right)\frac{\rho(\rho^\tau-1)}{\rho-1}\EE\|x^k-\bar{x}^{k+1}\|^2\\
&\textstyle + \frac{1}{m}\left(\frac{{\tau}}{\gamma}+\frac{1}{mp_{\min}}-\frac{1-\beta}{{\eta}}\right)\EE\|x^k-\bar x^{k+1}\|^2\\
= & \textstyle \big(1-\frac{\beta\mu\eta}{m}\big)\EE\|x^{k} - x^* \|^2 \cr
&\textstyle+\frac{1}{m}\left(\frac{2\beta\eta \mu\tau}{ m^2p_{\min}}\frac{\rho(\rho^\tau - 1)}{\rho - 1}+\frac{2}{m}\sqrt{\frac{\rho(\rho^\tau-1)\tau}{(\rho-1) p_{\min}}}+\frac{1}{ mp_{\min}}-\frac{1-\beta}{{\eta}}\right)\EE\|x^k-\bar x^{k+1}\|^2\cr
\le & \textstyle \big(1-\frac{\beta\mu\eta}{m}\big)\EE\|x^{k} - x^* \|^2,
\end{aligned}
\end{align*}
\cut{Letting $\gamma = m\sqrt{\frac{\tau (\rho - 1) p_{\min}}{ \rho(\rho^\tau - 1)}}$ in the above inequality, we have
\begin{align*}
&\mathbb{E}\big(\|x^{k+1} - x^* \|^2\big)  \cr
\le & \textstyle \big(1-\frac{\beta\mu\eta}{m}\big)\EE\|x^{k} - x^* \|^2 \cr
&\textstyle+\frac{1}{m}\left(\frac{2\beta\eta \mu\tau}{ m^2p_{\min}}\frac{\rho(\rho^\tau - 1)}{\rho - 1}+\frac{2}{m}\sqrt{\frac{\rho(\rho^\tau-1)\tau}{(\rho-1) p_{\min}}}+\frac{1}{ mp_{\min}}-\frac{1-\beta}{{\eta}}\right)\EE\|x^k-\bar x^{k+1}\|^2\cr
\le & \textstyle \big(1-\frac{\beta\mu\eta}{m}\big)\EE\|x^{k} - x^* \|^2,
% +{{\eta p_{\min}}\over m}\left(\mu\tau+{\gamma\over{\eta} }\right)\frac{\rho(\rho^\tau-1)}{\rho-1}\EE\|x^k-\bar{x}^{k+1}\|^2\\
%&\textstyle + {1\over m}\left({{\tau}\over \gamma}+{1\over mp_{\min}}-{1\over {2\eta}}\right)\EE\|x^k-\bar x^{k+1}\|^2.
\end{align*}}where we have let $\gamma = m\sqrt{\frac{\tau (\rho - 1) p_{\min}}{ \rho(\rho^\tau - 1)}}$ in the second equality, and the last inequality holds because of the choice of $\eta$. %is because $\frac{2\beta\eta \mu\tau}{ m^2p_{\min}}\frac{\rho(\rho^\tau - 1)}{\rho - 1}+\frac{2}{ m}\sqrt{\frac{\rho(\rho^\tau-1)\tau}{(\rho-1) p_{\min}}}+\frac{1}{ mp_{\min}}-\frac{1-\beta}{{\eta}}\le 0$ from the choice of $\eta$ and the quadratic formula.
%and $\eta \in (0, \frac{1}{2 (b + \sqrt{b^2 + a})}]$, where $a = 2 \mu \tau p_{\min} \frac{\rho^\tau - 1}{\rho - 1}$, $b = \frac{1}{2 m p_{\min}} + \sqrt{p_{\min}\tau \frac{\rho^\tau - 1}{\rho - 1}}$, then we have
%$${{2\eta p_{\min}}\over m}\left(2\mu\tau+{\gamma\over{2\eta} }\right)\frac{\rho^\tau-1}{\rho-1}+{1\over m}\left({{\tau}\over \gamma}+{1\over mp_{\min}}-{1\over {2\eta}}\right)\le 0,$$
%which implies
Therefore, \eqref{linear-convg} holds. %$$\mathbb{E}\big(\|x^{k+1} - x^* \|^2\big)\le \big(1-\frac{\mu\eta}{2m}\big)\EE\big(\|x^{k} - x^* \|^2 \big),$$
%and linear convergence is established.
\hfill\end{proof}

\begin{remark}\label{rmk:linear}
Assume $i_k$ is chosen uniformly at random, so $p_{\min}=\frac{1}{m}$. We consider the case when $m$ and $\tau$ are large. Let $\sqrt{\rho}=1+\frac{1}{\tau}$.  Then from the fact that $(1+\frac{1}{k})^k$ increasingly converges to the natural number $e$, we have from \eqref{eta-cond1} that $\underline{\eta}_1=O(\frac{\sqrt{m}}{\tau^2})$. In addition, note from \eqref{eta-cond2} that $a=O(b^2)=O(\frac{\tau^2}{m})$, and thus $\underline{\eta}_2=O(\frac{\sqrt{m}}{\tau})$. Therefore, if $\tau=O(m^{\frac{1}{4}})$, then the stepsize in Theorem \ref{thm:strongly_monotone} can be $\eta=O(1)$. Hence, linear speedup can be achieved.
\end{remark}

\end{section}


%\input{asyn_t_numerical}
\begin{section}{Experiments}\label{sec:numerical}
We illustrate the behavior of ARock for solving the $\ell_1$ regularized logistic regression problem. Our primary goal is to show the efficiency of the async-parallel implementation compared to the single-threaded implementation and the sync-parallel implementation.

Our experiments run on 1 to 32 threads on a machine with eight Quad-Core AMD Opteron\textsuperscript{TM} Processors (32 cores in total) and $64$ Gigabytes of RAM. All of the experiments were coded in C++ and OpenMP. We use the Eigen library\footnote{\url{http://eigen.tuxfamily.org}} for sparse matrix operations. Our codes as well as numerical results for other applications will be publicly available on the authors' website.

The running times and speedup ratios of  both sync-parallel and async-parallel algorithms are sensitive to a number of factors, such as the size of each coordinate update (granularity), sparsity  of the problem data,  compiler optimization flags,  and  operations that affect cache performance and memory access contention. In addition, since all  agents in the sync-parallel implementation  must wait for the last agent to finish an iteration, a large load imbalance will significantly degrade the performance. We do not have the space in this paper to present  numerical results under all variations of these cases. %However, our codes will be available online.
\cut{
\begin{subsection}{Linear equations}\label{sec:lin-eq}
We apply ARock presented in Subsection \ref{sec:linearequations} to solve linear equations $Ax=b$ on two synthetic datasets. In Dataset I, $A \in \mathbb{R}^{n\times n}$ is a banded matrix with a small bandwidth. Given the bandwidth $w$, we have  $a_{ij}=0$ whenever $|j-i|>w$. In Dataset II, $A\in \mathbb{R}^{n\times n}$ is a dense random symmetric matrix.
Both matrices $A$ are   diagonally dominant. They are summarized in Table \ref{tab:lq_data}. For each dataset, the entries of $x$ independently follow the standard Gaussian distribution, and  $b$ is generated as $Ax$.
\begin{table}[htbp]
\centering
 \begin{tabular}{lccc}
  \toprule
  Name  & Type & Size ($n$) & Bandwidth ($w$)\\
  \midrule
 Dataset I & Sparse & 1, 000, 000 & 5\\
  Dataset II & Dense & 5, 000 & N/A \\
  \bottomrule
 \end{tabular}
 \caption{\label{tab:lq_data}Linear equation datasets}
\end{table}

Throughout the experiments, $A$, $b$ and $x$ are shared variables. We set $\eta_k  = 1, \,\forall k$. %($\alpha \, \eta_k$) with length $1$.
Although the theory (Theorem~\ref{thm:convergence}) requires a smaller step size if more cores are used (usually leading to a larger delay $\tau$), we observe numerically that the unit step size can serve well for the algorithm on different numbers of cores. %remain the same even the number of cores increases.
The components of $x$ are  evenly partitioned into coordinates, and each coordinate contains $\frac{n}{\text{number of cores}}$ components of $x$, where the number of cores varies from  1 to 32. Each core is in charge of updating one pre-assigned coordinate\footnote{Deterministic assignments are applied since generating pseudo-random numbers is more expensive than updating a coordinate in this test. We tested random assignments and observed that they  give nearly the same epoch-versus-convergence performance but takes a much longer running time.}. We run both  sync-parallel Jacobi and ARock  algorithms to 100 epochs on Dataset I, and to 50 epochs on Dataset II, where an epoch is counted for every   $n$ updates to the components of $x$.

Figure~\ref{fig:lq_speedup} depicts how the size of residual, $\|Ax - b\|$, reduces over the wall-clock time. From the figure, we see that  both sync-parallel Jacobi and ARock show almost-linear speedup as the number of cores increases, but ARock takes much less time than sync-parallel Jacobi. To compare their  epoch convergence rates, we compute the size of residual  after each epoch (though they take different times to finish an epoch). Figure~\ref{fig:lq_residual} plots the sizes of residual at the end of different epochs.  Note that the  curves of sync-parallel Jacobi  running on different numbers of cores are nearly identical. Hence, only one curve is shown. The curves  of ARock running on different numbers of cores are also nearly identical, %Asynchronous Jacobi method with 32 cores has slightly worse residual curve than the others.
and they closely match with the reference curve of the standard Gauss-Seidel method, which is a sequential method. Therefore, ARock enjoys not only almost-linear speedup but also the Gauss-Seidel type of fast convergence, thanks to the asynchrony. This phenomenon was also  observed in the previous work~\cite{goos_numerical_1992,bertsekas1991some}.

\begin{figure}[!h]%\label{fig:lq_speedup}
        \centering
       \begin{subfigure}[b]{0.48\textwidth}
                \includegraphics[width=\textwidth]{./figs/linear_eqn/1_m_11_diag_time_vs_res_cropped}
                \caption{Dataset I}
        \end{subfigure}
        ~ %add desired spacing between images, e. g. ~, \quad, \qquad, \hfill etc.
          %(or a blank line to force the subfigure onto a new line)
        \begin{subfigure}[b]{0.48\textwidth}
                \includegraphics[width=\textwidth]{./figs/linear_eqn/dense_time_vs_residual}
                \caption{Dataset II}
        \end{subfigure}
        \caption{Residual $\|Ax - b\|$ versus  wall-clock time for sync-parallel Jacobi and ARock (async-parallel) running on 1--32 cores. All runs are terminated at 100 epochs on Dataset I and 50 epochs on Dataset II.}\label{fig:lq_speedup}
\end{figure}


\begin{figure}[!h]
        \centering
       \begin{subfigure}[b]{0.48\textwidth}
                \includegraphics[width=\textwidth]{./figs/linear_eqn/1_m_11_diag_residual_cropped}
                \caption{Dataset I}
        \end{subfigure}
        ~ %add desired spacing between images, e. g. ~, \quad, \qquad, \hfill etc.
          %(or a blank line to force the subfigure onto a new line)
        \begin{subfigure}[b]{0.48\textwidth}
                \includegraphics[width=\textwidth]{./figs/linear_eqn/ds_epoch_vs_res}
                \caption{Dataset II}
        \end{subfigure}
        \caption{Epoch-convergence (regardless of time) of sync-parallel Jacobi and ARock (async-parallel). The black residual curve of the standard Gauss-Seidel method is plotted for reference.}\label{fig:lq_residual}
\end{figure}




%\begin{remark}
%If the data matrices are extremely sparse, generating a random number can be more expensive than evaluating one Jacobi iteration. So instead of selecting coordinates randomly, for performance purpose, we can actually partition $x$ into equal size blocks, and each thread updates its own partition of $x$ in a cyclic fashion. This can also preserve better data locality.
%\end{remark}

\end{subsection}
}

\begin{subsection}{$\ell_1$ regularized logistic regression}
In this subsection, we apply ARock with the update \eqref{asyncFBS} to the $\ell_1$ regularized logistic regression problem:
\begin{equation}\label{eqn:log}
\Min_{x \in \mathbb{R}^n} \lambda \|x\|_1 + \frac{1}{N} \sum_{i=1}^N \log\big(1 + \exp(-b_i \cdot a_i^T x)\big),
\end{equation}
where $\{(a_i, b_i)\}_{i=1}^N$ is the set of sample-label pairs with $b_i \in \{1, -1\}$, $\lambda=0.0001$, and $n$ and $N$ represent the numbers of features and samples, respectively. This test uses the datasets\footnote{\url{http://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/}}: rcv1 and news20, which are summarized in Table \ref{tab:log_data}.

\begin{table}[htbp]
\centering
{\small \begin{tabular}{crrr}
\hline
  Name & \# samples & \# features & \# nonzeros in $\{a_1,\ldots,a_N\}$ \\
  \hline
 rcv1 & 20, 242 & 47, 236 & 1, 498, 952\\
  news20 & 19, 996 & 1, 355, 191 & 9, 097, 916\\
  \hline
 \end{tabular}}
 \caption{\label{tab:log_data}Two  datasets for  sparse logistic regression.}
 \vspace{-8mm}
\end{table}

\begin{figure}[!h]%\label{fig:lq_speedup}
        \centering
       \begin{subfigure}[b]{0.4\textwidth}
                \includegraphics[width=\textwidth]{./figs/logistic/rcv1}
                \caption{ dataset rcv1}
        \end{subfigure}
        ~ %add desired spacing between images, e. g. ~, \quad, \qquad, \hfill etc.
          %(or a blank line to force the subfigure onto a new line)
        \begin{subfigure}[b]{0.4\textwidth}
                \includegraphics[width=\textwidth]{./figs/logistic/news20}
                \caption{  dataset news20}
        \end{subfigure}
        \vspace{-3mm}
        \caption{The distribution of coordinate sparsity. Each dot represents the total number of nonzeros in the vectors $a_i$ that correspond to each coordinate. The large distribution in (b) is responsible  for the large load imbalance and thus the poor sync-parallel performance.}\label{fig:sparsity}
 \vspace{-3mm}
\end{figure}

We let each coordinate hold roughly 50 features. {Since the total number of features is not divisible by 50, some coordinates have 51 features.} We let each agent draw a coordinate uniformly at random at each iteration. We stop all the tests after 100 epochs since they have nearly identical progress per iteration. %The size of the block of coordinates is set to 50. In each iteration, each core randomly draws a block and updates it.
The step size is set to $\eta_k=0.9,\,\forall k$. Let $A = [a_1, \ldots, a_N]^T$ and $b = [b_1, ..., b_N]^T$. In global memory, we store $A,~b$, and $x$. We also store the product $Ax$ in global memory so that the forward step can be efficiently computed. Whenever a coordinate of $x$ gets updated, $Ax$ is immediately updated at a low cost. Note that if $Ax$ is \emph{not} stored in global memory, every coordinate update will have to compute $Ax$ from scratch, which involves the entire $x$ and will be very expensive.  %Note that the renewal of $Ax$ can be computed  efficiently, and it enables ARock to have the same per epoch computational complexity as that of the sequential FBS.

Table \ref{tab:log_time} gives the running times of  the sync-parallel and ARock (async-parallel) implementations on the two datasets. We can observe that ARock achieves almost-linear speedup, but sync-parallel scales very poorly as we explain below.

In the sync-parallel implementation,  all the running cores have to wait for the last core to finish an iteration, and therefore if a core has a large load, it slows down the iteration. Although every core is (randomly) assigned to roughly the same number of features (either 50 or 51 components of $x$) at each iteration, their  $a_i$'s have very different numbers of nonzeros (see Figure \ref{fig:sparsity} for the distribution), and the core with the largest number of nonzeros is the slowest (Sparse matrix computation is used for both datasets, which are very large.) As more cores are used,  despite that they altogether do more work at each iteration,  the per-iteration time reduces as the slowest core tends to be slower. The very large imbalance of load explains why the 32 cores only give speedup ratios of 4.0 and 1.3 in Table \ref{tab:log_time}.

On the other hand, being asynchronous, ARock does not suffer from the  load imbalance. Its performance grows nearly linear with the number of cores. In theory, a large load imbalance may cause a large \cut{maximum delay, }$\tau$, and thus a small $\eta_k$. However, the uniform $\eta_k=0.9$ works well in all the tests, possibly because the $a_i$'s are sparse.

Finally, we have observed that the progress toward solving \eqref{eqn:log} is mainly a function of the number of epochs and does not change appreciably  when the number of cores increases or between sync-parallel and async-parallel. Therefore, we always stop at 100 epochs.


\begin{table}[htbp]
\centering
{\small \begin{tabular}{|c|r|r|r|r|r|r|r|r|}
  \hline
  \multirow{3}{*}{\# cores} & \multicolumn{4}{|c|}{rcv1} & \multicolumn{4}{c|}{news20} \\
  \cline{2-9}
  & \multicolumn{2}{|c|}{Time (s)} &  \multicolumn{2}{c|}{Speedup} &  \multicolumn{2}{c|}{Time (s)} & \multicolumn{2}{c|}{Speedup}\\
  \cline{2-9}
  & async & sync &  async & sync &  async & sync &  async & sync \\
  \hline
   1 &   122.0 &  122.0 & 1.0   & 1.0 & 591.1   & 591.3 & 1.0   & 1.0\\
   2 &   63.4   &  104.1 & 1.9   & 1.2 & 304.2   & 590.1 & 1.9   & 1.0\\
   4 &   32.7   &  83.7   & 3.7   & 1.5 & 150.4   & 557.0 & 3.9   & 1.1\\
   8 &   16.8   &  63.4   & 7.3   & 1.9 & 78.3     & 525.1 & 7.5   & 1.1\\
   16 & 9.1     &  45.4   & 13.5 & 2.7 & 41.6     & 493.2 & 14.2 & 1.2\\
   32 & 4.9     &  30.3   & 24.6 & 4.0 & 22.6     & 455.2 & 26.1 & 1.3\\
  \hline
 \end{tabular}}
 \caption{\label{tab:log_time}Running times of  ARock (async-parallel) and sync-parallel FBS implementations for the $\ell_1$ regularized logistic regression on two datasets. Sync-parallel has a very poor speedup  due to the large distribution of coordinate sparsity (Figure \ref{fig:sparsity}) and thus the large load imbalance across cores.}
 \vspace{-4mm}
\end{table}
\end{subsection}
\end{section}


%\input{asyn_t_conclu}
\begin{section}{Conclusion}
We have proposed an async-parallel framework, ARock, for finding a fixed-point of a nonexpansive operator by coordinate updates. We establish the almost sure weak and strong convergence, linear convergence rate and almost-linear speedup of ARock under certain assumptions.
Preliminary numerical results on real data illustrate the high efficiency of the proposed framework compared to the traditional parallel (sync-parallel) algorithms.
\end{section}

\begin{section}{Acknowledgements}
We would like to thank Brent Edmunds for offering invaluable suggestions on the organization and writing of this paper. We would also like to thank Robert Hannah for coming up with the dual-memory approach.
The authors are grateful to Kun Yuan for helpful discussions on decentralized optimization. \end{section}


\bibliographystyle{spmpsci}
\bibliography{asyn}

%\input{asyn_t_appendix}
{\small
\appendix
\section{Derivation of certain updates}\label{sec:appendix}
We show in details how to obtain the updates in \eqref{eqn:relaxed_prs_set} and \eqref{eqn:asyn_admm}.

\subsection{Derivation of updates in \eqref{eqn:relaxed_prs_set}}
Let $x=(x_1,\ldots,x_m)\in \cH^m$, $$f(x):=\sum_{i=1}^m\cI_{C_i}(x_i),\qquad g(x):=\cI_{\{x_1=\cdots=x_m\}}(x),$$
where $g(x)$ equals $0$ if $x_1=\cdots=x_m$ and $\infty$ otherwise.
Then \eqref{eqn:relaxed_prs-a} reduces to
\begin{align*}
\hat{x}^k=&\argmin_{z\in\cH^m}\, g(z) +\frac{1}{2\gamma}\|z-\hat{z}^k\|^2=\argmin_{z\in\cH^m: z_1=\cdots=z_m} \|z-\hat{z}^k\|^2\\
=&\argmin_{z\in\cH^m: z_1=\cdots=z_m} \sum_{i=1}^m\|z_1-\hat{z}_i^k\|^2
=\left(\frac{1}{m}\sum_{i=1}^m\hat{z}_i^k,\ldots,\frac{1}{m}\sum_{i=1}^m\hat{z}_i^k\right)\in\cH^m,
\end{align*}
where the last equality is obtained by noting that $z_1=\frac{1}{m}\sum_{i=1}^m\hat{z}_i^k$ is the unique minimizer of $\sum_{i=1}^m\|z_1-\hat{z}_i^k\|^2$. Next, \eqref{eqn:relaxed_prs-b} reduces to
\begin{align*}
\hat{y}^k=\argmin_{z\in\cH^m} f(z) +\frac{1}{2\gamma}\|z-(2\hat{x}^k-\hat{z}^k)\|^2
=\argmin_{z: z_i\in C_i,\forall i} \sum_{i=1}^m\|z_i-(2\hat{x}_i^k-\hat{z}_i^k)\|^2.
\end{align*}
It is easy to see that $\hat{y}_i^k= \Proj_{C_i}(2\hat{x}_i^k-\hat{z}_i^k),\,\forall i$.

Since \eqref{eqn:relaxed_prs-c} only updates the $i_k$th coordinate of $z$, we only need $\hat{x}_{i_k}^k$ and $\hat{y}_{i_k}^k$, and thus in \eqref{eqn:relaxed_prs_set:a} and \eqref{eqn:relaxed_prs_set:b}, we only compute $\hat{x}_{i_k}^k$ and $\hat{y}_{i_k}^k$. Plugging the above $\hat{x}^k$ and $\hat{y}^k$ into \eqref{eqn:relaxed_prs-c} gives \eqref{eqn:relaxed_prs_set:c} directly.

\subsection{Derivation of \eqref{eqn:asyn_admm}}
We first show how to get \eqref{eq:lag-dual}. The Lagrangian of \eqref{eqn:admm_form} is
$L(x,y,w)= f(x)+g(y)-\langle w, Ax+By-b\rangle,$
and the Lagrange dual function is
\begin{align*}
d(w)=&\min_{x\in \cH_1,y\in\cH_2}L(x,y,w)\\
=&\big(\min_{x\in\cH_1} f(x)-\langle A^*w, x\rangle\big) +\big( \min_{y\in\cH_2} g(y)-\langle B^*w, y\rangle\big) +\langle w, b\rangle\\
=&-\big(\max_{x\in\cH_1} -f(x)+\langle A^*w, x\rangle\big) -\big( \max_{y\in\cH_2} -g(y)+\langle B^*w, y\rangle\big) +\langle w, b\rangle\\
=&-f^*(A^*w)-g^*(B^*w)+\langle w, b\rangle,
\end{align*}
where the last equality is from the definition of convex conjugate: $f^*(z)=\max_x \langle z, x\rangle-f(x).$ Hence, the dual problem is $\max_{w} d(w)$, which is equivalent to \eqref{eq:lag-dual}.

Secondly, we show why $z^+=\prox_{\gamma\cdot d_g}(z)$ is given by \eqref{dualg}. Note
\begin{align*}
\textstyle\min_s d_g(s)+\frac{1}{2\gamma}\|s-z\|^2
=&\textstyle\min_s g^*(B^*s)-\langle s, b\rangle +\frac{1}{2\gamma}\|s-z\|^2\\
=&\textstyle\min_s\max_y \langle B^*s, y\rangle - g(y)-\langle s, b\rangle +\frac{1}{2\gamma}\|s-z\|^2\\
=&\textstyle\max_y\min_s \langle B^*s, y\rangle - g(y)-\langle s, b\rangle +\frac{1}{2\gamma}\|s-z\|^2\\
=&\textstyle\max_y\min_s \langle s, By-b\rangle - g(y) +\frac{1}{2\gamma}\|s-z\|^2\\
=&\textstyle\max_y -g(y)+\langle z, By-b\rangle -\frac{\gamma}{2}\|By-b\|^2\\
=&\textstyle-\min_y g(y)-\langle z, By-b\rangle +\frac{\gamma}{2}\|By-b\|^2,
\end{align*}
where the fifth equality holds because $s^*=z-\gamma(By-b) =\arg\min_s\langle s, By-b\rangle  +\frac{1}{2\gamma}\|s-z\|^2.$
Hence, by the definition of the proximal operator and the above arguments, we have that $z^+=\prox_{\gamma\cdot d_g}(z)$ can be obtained from \eqref{dualg}. Then \eqref{dualf} is from \eqref{dualg} through replacing $g$ to $f$, $B$ to $A$, and $b$ to $0$.

Finally, it is straightforward to have \eqref{eqn:asyn_admm} by plugging \eqref{dualf} and \eqref{dualg} into \eqref{eqn:relaxed_prs}.

\section{Derivation of async-parallel ADMM for decentralized optimization}\label{sec:appendix2}
This section describes how to implement the updates \eqref{eqn:asyn_admm} for the model~\eqref{eqn:des_prob3}.

In~\eqref{eqn:des_prob3}, $g(y)$ and $b$ vanish and, corresponding to the two constraints $x_i = y_{ij}$ and $x_j = y_{ij}$, the two rows of matrices $A$ and $B$ are
$\begin{bmatrix}
     \cdots &  1 & \cdots  & \cdots  & \cdots  \\
     \cdots & \cdots & \cdots & 1 & \cdots
  \end{bmatrix}\quad
  \begin{bmatrix}
     \cdots  & -1 & \cdots  \\
     \cdots  & -1 & \cdots
  \end{bmatrix}, $
where $\cdots$ are zeros, the two coefficients 1 correspond to $x_i$ and $x_j$, and the two coefficients $-1$ correspond to $y_{ij}$.
%There is only one non-zero component ``$-1$" in each row of $B$.
Then, \eqref{eqn:asyn_admm:a} and~\eqref{eqn:asyn_admm:b} can be calculated as
\begin{align*}
\hat y^k_{li} &= (\hat z_{li,l}^k+\hat z_{li,i}^k)/(2\gamma)\quad\forall l\in L(i),\\
 (\hat w_g^k)_{li,i} & = (\hat z^k_{li,i}-\hat z^k_{li,r})/2\quad\forall l\in L(i),\\
\hat y^k_{ir} &= (\hat z_{ir,i}^k+\hat z_{ir,r}^k)/(2\gamma)\quad\forall r\in R(i),\\
 (\hat w_g^k)_{ir,i} & = (\hat z^k_{ir,i}-\hat z^k_{ir,r})/2\quad\forall r\in R(i).
\end{align*}
In addition, $\hat x^k_i$ can be obtained by solving~\eqref{eqn:asyn_admm_des_a}, and both $z_{li,i}^{k+1}$ and $z_{ir,i}^{k+1}$ can be updated from~\eqref{eqn:asyn_admm_des_b} and~\eqref{eqn:asyn_admm_des_c}.

Furthermore, as mentioned in Section~\ref{subsec:async_admm_des}, we can derive another version of async-parallel ADMM for decentralized optimization, which reduces to the algorithm in \cite{wei2013on}, by activating an edge $(i,j)\in E$ instead of an agent $i$ each time. In this version, the agents $i$ and $j$ associated with the edge $(i,j)$ must also be activated. Here we derive the update~\eqref{eqn:asyn_admm} for the model~\eqref{eqn:des_prob3} with the update order of $x$ and $y$ swapped. Following~\eqref{eqn:asyn_admm} we obtain the following steps whenever an edge $(i,j)\in E$ is activated:
\begin{subequations}
\begin{align*}
\hat x^{k}_{i} & = \arg\min_{x_i} f_i(x_i) -  \big(\sum_{l \in L(i)}\hat z^k_{li, i}  + \sum_{r\in R(i)}\hat z^k_{ir, i} \big) x_i +  \frac{\gamma}{2} |E(i)|\cdot \|x_i\|^2 \\
\hat x^{k}_{j}  &= \arg\min_{x_j} f_j(x_j) -  \big(\sum_{l \in L(j)}\hat z^k_{lj, j}  + \sum_{r\in R(j)}\hat z^k_{jr, j} \big) x_j +  \frac{\gamma}{2} |E(j)|\cdot \|x_j\|^2 \\
(\hat w_{f}^{k})_{ij,i} &= \hat z_{ij,i}^{k} - \gamma \hat x_i^{k} \\
(\hat w_{f}^{k})_{ij,j} &= \hat z_{ij,j}^{k} - \gamma \hat x_j^{k} \\
\hat y_{ij}^{k} &= \arg\min_{y_{ij}}  \langle 2 (\hat w_{f}^{k})_{ij,i} -\hat z_{ij,i}^{k}+ 2(\hat w_{f}^{k})_{ij,j} -\hat z_{ij,j}^{k}, y_{ij} \rangle + \frac{\gamma}{2} \|y_{ij}\|^2  \\
(\hat w_{g}^{k})_{ij,i} &=2 (\hat w_{f}^{k})_{ij,i} -\hat z_{ij,i}^{k} + \gamma  \hat y_{ij}^k\\
(\hat w_{g}^{k})_{ij,j} &=2 (\hat w_{f}^{k})_{ij,j} -\hat z_{ij,j}^{k} + \gamma  \hat y_{ij}^k\\
z_{ij,i}^{k+1} &= z_{ij,i}^{k} + \eta_k ((\hat w_{g}^{k})_{ij,i} - (\hat w_{f}^{k})_{ij,i})\\
z_{ij,j}^{k+1} &= z_{ij,j}^{k} + \eta_k ((\hat w_{g}^{k})_{ij,j} - (\hat w_{f}^{k})_{ij,j}).
\end{align*}
\end{subequations}
Every agent $i$ in the network maintains the dual variables $z_{li,i}$, $l\in L(i)$, and $z_{ir,i}$, $r\in R(i)$, and the variables $x,y,w$ are intermediate and do not need to be maintained between the activations. When an edge $(i,j)$ is activated, the agents $i$ and $j$ first compute their $\{\hat x_i^k,(\hat w_{f}^{k})_{ij,i}\}$ and $\{\hat x_j^k,(\hat w_{f}^{k})_{ij,j}\}$ independently and respectively, then they collaboratively compute $\hat y_{ij}^k$, and finally they update their own $z_{ij,i}^k$ and $z_{ij,j}^k$, respectively. We allow adjacent edges (which share agents) to be activated in a short period of time when their updates are possibly overlapped in time. When $\tau=0$, i.e., there is no simultaneous activation or overlap, it reduces to the algorithm in~\cite{wei2013on}.
}


\end{document}



diff.tex
Running: latexdiff   "../ARock_arxiv_v3/asyn_theory_MP.tex" "asyn_theory_MP.tex" > "asyn_theory_MP-diff.tex"

macros.tex
%% macros for commenting
\usepackage[normalem]{ulem} % to use \sout

%% math packages
% \usepackage{amsmath} % not needed for Springer
\usepackage{amssymb}
\usepackage{mathrsfs}
\usepackage{amssymb,amsmath,amsfonts}
\usepackage{mathtools} % Bonus

%\usepackage{subfigure}
%\usepackage{url}

%% a light-weight algorithm environment
\newtheorem{algo}{Algorithm}

%% highlighting and commenting
\newcommand{\outline}[1]{{\color{brown}#1}}
\newcommand{\rev}[1]{{\color{blue}#1}}
\newcommand{\remove}[1]{{\sout{#1}}}
\newcommand{\cut}[1]{{}}
\newcommand{\commmy}[1]{{\color{red}(#1 -- Ming)}} % comments
\newcommand{\commwy}[1]{{\color{red}(#1 -- Wotao)}} % comments
\newcommand{\commyx}[1]{{\color{red}(#1 -- Yangyang)}} % comments
\newcommand{\commzp}[1]{{\color{red}(#1 -- Zhimin)}} % comments

%% macros for letters

\newcommand{\va}{{\mathbf{a}}}
\newcommand{\vb}{{\mathbf{b}}}
\newcommand{\vc}{{\mathbf{c}}}
\newcommand{\vd}{{\mathbf{d}}}
\newcommand{\ve}{{\mathbf{e}}}
\newcommand{\vf}{{\mathbf{f}}}
\newcommand{\vg}{{\mathbf{g}}}
\newcommand{\vh}{{\mathbf{h}}}
\newcommand{\vi}{{\mathbf{i}}}
\newcommand{\vj}{{\mathbf{j}}}
\newcommand{\vk}{{\mathbf{k}}}
\newcommand{\vl}{{\mathbf{l}}}
\newcommand{\vm}{{\mathbf{m}}}
\newcommand{\vn}{{\mathbf{n}}}
\newcommand{\vo}{{\mathbf{o}}}
\newcommand{\vp}{{\mathbf{p}}}
\newcommand{\vq}{{\mathbf{q}}}
\newcommand{\vr}{{\mathbf{r}}}
\newcommand{\vs}{{\mathbf{s}}}
\newcommand{\vt}{{\mathbf{t}}}
\newcommand{\vu}{{\mathbf{u}}}
\newcommand{\vv}{{\mathbf{v}}}
\newcommand{\vw}{{\mathbf{w}}}
\newcommand{\vx}{{\mathbf{x}}}
\newcommand{\vy}{{\mathbf{y}}}
\newcommand{\vz}{{\mathbf{z}}}
%
%\newcommand{\ta}{{\tilde{a}}}
%\newcommand{\tb}{{\tilde{b}}}
%\newcommand{\tc}{{\tilde{c}}}
%\newcommand{\td}{{\tilde{d}}}
%\newcommand{\te}{{\tilde{e}}}
%\newcommand{\tf}{{\tilde{f}}}
%\newcommand{\tg}{{\tilde{g}}}
%\newcommand{\th}{{\tilde{h}}}
%\newcommand{\ti}{{\tilde{i}}}
%\newcommand{\tj}{{\tilde{j}}}
%\newcommand{\tk}{{\tilde{k}}}
%\newcommand{\tl}{{\tilde{l}}}
%\newcommand{\tm}{{\tilde{m}}}
%\newcommand{\tn}{{\tilde{n}}}
%\newcommand{\to}{{\tilde{o}}}
%\newcommand{\tp}{{\tilde{p}}}
%\newcommand{\tq}{{\tilde{q}}}
%\newcommand{\tr}{{\tilde{r}}}
%\newcommand{\ts}{{\tilde{s}}}
%\newcommand{\tt}{{\tilde{t}}}
%\newcommand{\tu}{{\tilde{u}}}
%\newcommand{\tv}{{\tilde{v}}}
%\newcommand{\tw}{{\tilde{w}}}
%\newcommand{\tx}{{\tilde{x}}}
%\newcommand{\ty}{{\tilde{y}}}
%\newcommand{\tz}{{\tilde{z}}}

\newcommand{\vA}{{\mathbf{A}}}
\newcommand{\vB}{{\mathbf{B}}}
\newcommand{\vC}{{\mathbf{C}}}
\newcommand{\vD}{{\mathbf{D}}}
\newcommand{\vE}{{\mathbf{E}}}
\newcommand{\vF}{{\mathbf{F}}}
\newcommand{\vG}{{\mathbf{G}}}
\newcommand{\vH}{{\mathbf{H}}}
\newcommand{\vI}{{\mathbf{I}}}
\newcommand{\vJ}{{\mathbf{J}}}
\newcommand{\vK}{{\mathbf{K}}}
\newcommand{\vL}{{\mathbf{L}}}
\newcommand{\vM}{{\mathbf{M}}}
\newcommand{\vN}{{\mathbf{N}}}
\newcommand{\vO}{{\mathbf{O}}}
\newcommand{\vP}{{\mathbf{P}}}
\newcommand{\vQ}{{\mathbf{Q}}}
\newcommand{\vR}{{\mathbf{R}}}
\newcommand{\vS}{{\mathbf{S}}}
\newcommand{\vT}{{\mathbf{T}}}
\newcommand{\vU}{{\mathbf{U}}}
\newcommand{\vV}{{\mathbf{V}}}
\newcommand{\vW}{{\mathbf{W}}}
\newcommand{\vX}{{\mathbf{X}}}
\newcommand{\vY}{{\mathbf{Y}}}
\newcommand{\vZ}{{\mathbf{Z}}}

\newcommand{\cA}{{\mathcal{A}}}
\newcommand{\cB}{{\mathcal{B}}}
\newcommand{\cC}{{\mathcal{C}}}
\newcommand{\cD}{{\mathcal{D}}}
\newcommand{\cE}{{\mathcal{E}}}
\newcommand{\cF}{{\mathcal{F}}}
\newcommand{\cG}{{\mathcal{G}}}
\newcommand{\cH}{{\mathcal{H}}}
\newcommand{\cI}{{\mathcal{I}}}
\newcommand{\cJ}{{\mathcal{J}}}
\newcommand{\cK}{{\mathcal{K}}}
\newcommand{\cL}{{\mathcal{L}}}
\newcommand{\cM}{{\mathcal{M}}}
\newcommand{\cN}{{\mathcal{N}}}
\newcommand{\cO}{{\mathcal{O}}}
\newcommand{\cP}{{\mathcal{P}}}
\newcommand{\cQ}{{\mathcal{Q}}}
\newcommand{\cR}{{\mathcal{R}}}
\newcommand{\cS}{{\mathcal{S}}}
\newcommand{\cT}{{\mathcal{T}}}
\newcommand{\cU}{{\mathcal{U}}}
\newcommand{\cV}{{\mathcal{V}}}
\newcommand{\cW}{{\mathcal{W}}}
\newcommand{\cX}{{\mathcal{X}}}
\newcommand{\cY}{{\mathcal{Y}}}
\newcommand{\cZ}{{\mathcal{Z}}}

\newcommand{\ri}{{\mathrm{i}}}
\newcommand{\rr}{{\mathrm{r}}}

\newcommand{\EE}{{\mathbb{E}}}


%% macros for math notions and operators

\newcommand{\RR}{\mathbb{R}}
\newcommand{\CC}{\mathbb{C}}
\newcommand{\ZZ}{\mathbb{Z}}
\renewcommand{\SS}{{\mathbb{S}}}
\newcommand{\SSp}{\mathbb{S}_{+}}
\newcommand{\SSpp}{\mathbb{S}_{++}}
\newcommand{\sign}{\mathrm{sign}}
\newcommand{\vzero}{\mathbf{0}}
\newcommand{\vone}{{\mathbf{1}}}

\newcommand{\st}{{\text{s.t.}}} % subject to
\newcommand{\St}{{\mathrm{subject~to}}} % subject to
\newcommand{\op}{{\mathrm{op}}} % subscript for operator norm
\newcommand{\opt}{{\mathrm{opt}}} % subscript for optimal solution
%\newcommand{\supp}{{\mathrm{supp}}} % support
\newcommand{\Prob}{{\mathrm{Prob}}} % probability
\newcommand{\Diag}{{\mathrm{Diag}}} % vector -> diagonal matrix
%\newcommand{\diag}{{\mathrm{diag}}} % matrix diagonal -> vector
\newcommand{\dom}{{\mathrm{dom}}} % domain
\newcommand{\range}{{\mathrm{range}}} % domain
%\newcommand{\grad}{{\nabla}}    % gradient
\newcommand{\tr}{{\mathrm{tr}}} % trace
\newcommand{\TV}{{\mathrm{TV}}} % total variation
\newcommand{\Proj}{{\mathrm{Proj}}}
\newcommand{\prj}{{\mathrm{prj}}}
\newcommand{\prox}{\mathbf{prox}}
\newcommand{\refl}{\mathbf{refl}}
\newcommand{\reflh}{\refl^{\bH}}
\newcommand{\proxh}{\prox^{\bH}}
\newcommand{\minimize}{\text{minimize}}
\newcommand{\bgamma}{\boldsymbol{\gamma}}
\newcommand{\bsigma}{\boldsymbol{\sigma}}
\newcommand{\bomega}{\boldsymbol{\omega}}
\newcommand{\blambda}{\boldsymbol{\lambda}}
\newcommand{\bH}{\vH}
\newcommand{\bbH}{\mathbb{H}}
\newcommand{\bB}{\boldsymbol{\cB}}
\newcommand{\Tau}{\mathrm{T}}
\newcommand{\tnabla}{\widetilde{\nabla}}
\newcommand{\TDRS}{T_{\mathrm{DRS}}}
\newcommand{\TPRS}{T_{\mathrm{PRS}}}
\newcommand{\TFBS}{T_{\mathrm{FBS}}}
\newcommand{\best}{\mathrm{best}}
\newcommand{\kbest}{k_{\best}}
\newcommand{\diff}{\mathrm{diff}}
\newcommand{\barx}{\bar{x}}
\newcommand{\xgbar}{\bar{x}_g}
\newcommand{\xfbar}{\bar{x}_f}
\newcommand{\hatxi}{\hat{\xi}}
\newcommand{\xg}{x_g}
\newcommand{\xf}{x_f}
\newcommand{\du}{\mathrm{d}u}
\newcommand{\dy}{\mathrm{d}y}
\newcommand{\kconvergence}{\stackrel{k \rightarrow \infty}{\rightarrow }}
\DeclareMathOperator{\shrink}{shrink} % shrinkage
\DeclareMathOperator*{\argmin}{arg\,min}
\DeclareMathOperator*{\argmax}{arg\,max}
\DeclareMathOperator*{\Min}{minimize}
\DeclareMathOperator*{\Max}{maximize}
\DeclareMathOperator*{\Fix}{Fix}
\DeclareMathOperator*{\zer}{zer}
\DeclareMathOperator*{\nablah}{\nabla^{\bH}}
\DeclareMathOperator*{\gra}{gra}
\DeclarePairedDelimiter{\dotpb}{\langle}{\rangle_{\bH}}
\DeclarePairedDelimiter{\dotpv}{\langle}{\rangle_{\vH}}
\DeclarePairedDelimiter{\dotp}{\langle}{\rangle}

%% macros for environments math equations

\newcommand{\MyFigure}[1]{../fig/#1}

\newcommand{\bc}{\begin{center}}
\newcommand{\ec}{\end{center}}

\newcommand{\bdm}{\begin{displaymath}}
\newcommand{\edm}{\end{displaymath}}

\newcommand{\beq}{\begin{equation}}
\newcommand{\eeq}{\end{equation}}

\newcommand{\bfl}{\begin{flushleft}}
\newcommand{\efl}{\end{flushleft}}

\newcommand{\bt}{\begin{tabbing}}
\newcommand{\et}{\end{tabbing}}

\newcommand{\beqn}{\begin{align}}
\newcommand{\eeqn}{\end{align}}

\newcommand{\beqs}{\begin{align*}} % no equation numbers
\newcommand{\eeqs}{\end{align*}}  % no equation numbers

%% macros for theorem-like environments

%\newtheorem{theorem}{Theorem}
\newtheorem{assumption}{Assumption}
%\newtheorem{condition}{Condition}
%\newtheorem{rul}{Rule}
%\newtheorem{definition}{Definition}
%\newtheorem{corollary}{Corollary}
%\newtheorem{remark}{Remark}
%\newtheorem{lemma}{Lemma}
%\newtheorem{proposition}{Proposition}
%\newtheorem{example}{Example}
%\newtheorem{proof}{Proof}
%
%\newtheorem{example}[remark]{Example}




newcommands.tex
%\newcommand{\ra}[1]{\renewcommand{\arraystretch}{#1}}
%\newcommand{\pitem}[1]{\item [\mbox{[#1]}]}
%\newcommand{\btheta}{\ensuremath{\boldsymbol{\theta}}}
%\newcommand{\E}[1]{\ensuremath{\mathbb{E}_{#1}}}
%\newcommand{\piL}{\ensuremath{\pi_\text{L}}}
%\newcommand{\piU}{\ensuremath{\pi_\text{U}}}
%\newcommand{\red}[1]{{\color{red}{#1}}}
%\newcommand{\blue}[1]{{\color{blue}{#1}}}
%\newcommand{\green}[1]{{\color{green}{#1}}}
%\newcommand{\mybox}[2][]{{\begin{center}\fbox{\parbox{0.9\textwidth}{\textbf{#1}\\ #2}}\end{center}}}
%%\newcommand{\incomplete}{\mybox{\red{\[\text{\bfseries Complete this section!}\]}}}
%\newcommand{\incomplete}{\mybox{\begin{center}\bfseries \red{Complete this section!}\end{center}}}
%\newcommand{\multibox}[2][0.9]{\begin{center}\fbox{\parbox{#1\textwidth}{#2}}\end{center}}
%\newcommand{\parstart}[2]{\lettrine[lraise=-0.1, loversize=0.2]{#1}{\ \MakeUppercase{#2}}}
%\newcommand{\matlab}{\textsc{Matlab}}
%\newcommand{\fig}[1]{\vref{#1}}
%\newcommand{\eq}[1]{\ref{#1}}
%\labelformat{equation}{(#1)}
%\labelformat{figure}{Figure~#1}
%\labelformat{table}{Table~#1}
%\labelformat{section}{section~#1}
%\labelformat{subsection}{section~#1}
%\labelformat{subsubsection}{section~#1}
%\labelformat{chapter}{Chapter~#1}
%\labelformat{theorem}{Theorem~#1}
%\labelformat{definition}{Definition~#1}
%\labelformat{proposition}{Proposition~#1}
%\labelformat{remark}{Remark~#1}
%\labelformat{lemma}{Lemma~#1}

\DeclareMathOperator*{\argmin}{\arg\min}

% begin of personal macros
\newcommand{\half}{{\textstyle \frac{1}{2}}}
\newcommand{\eps}{\varepsilon}
\newcommand{\myth}{\vartheta}
\newcommand{\myphi}{\varphi}

\newcommand{\IN}{\mathbb{N}}
\newcommand{\IZ}{\mathbb{Z}}
\newcommand{\IQ}{\mathbb{Q}}
\newcommand{\IR}{\mathbb{R}}
\newcommand{\IC}{\mathbb{C}}
\newcommand{\Real}[1]{\mathrm{Re}\left({#1}\right)}
\newcommand{\Imag}[1]{\mathrm{Im}\left({#1}\right)}

\newcommand{\norm}[2]{\|{#1}\|_{{}_{#2}}}
\newcommand{\abs}[1]{\left|{#1}\right|}
\newcommand{\ip}[2]{\left\langle {#1}, {#2} \right\rangle}
\newcommand{\der}[2]{\frac{\partial {#1}}{\partial {#2}}}
\newcommand{\dder}[2]{\frac{\partial^2 {#1}}{\partial {#2}^2}}

\newcommand{\nn}{\mathbf{n}}
\newcommand{\xx}{\mathbf{x}}
\newcommand{\uu}{\mathbf{u}}
\newcommand{\junk}[1]{{}}

%\newcommand{\E}[1]{\ensuremath{\mathbb{E}_{#1}}}
\newcommand{\mb}{\ensuremath{\overline{m}}}
\newcommand{\mt}{\ensuremath{\tilde{m}}}
\newcommand{\nb}{\ensuremath{\overline{n}}}
\newcommand{\nt}{\ensuremath{\tilde{n}}}
\newcommand{\hb}{\ensuremath{\overline{h}}}
\newcommand{\htilde}{\ensuremath{\tilde{h}}}
\newcommand{\qb}{\ensuremath{\overline{q}}}
\newcommand{\qt}{\ensuremath{\tilde{q}}}
\newcommand{\minf}{\ensuremath{m_\infty}}
\newcommand{\ninf}{\ensuremath{n_\infty}}
\newcommand{\hinf}{\ensuremath{h_\infty}}
\newcommand{\qinf}{\ensuremath{q_\infty}}
\newcommand{\taum}{\ensuremath{\tau_m}}
\newcommand{\taun}{\ensuremath{\tau_n}}
\newcommand{\tauh}{\ensuremath{\tau_h}}
\newcommand{\tauq}{\ensuremath{\tau_q}}

\newcommand{ \bl }{\ensuremath{\big (}}
\newcommand{ \br }{\ensuremath{\big )}}


\newcommand{\INa}{\ensuremath{I_{\text{Na}}}}
\newcommand{\VNa}{\ensuremath{V_{\text{Na}}}}
\newcommand{\gNa}{\ensuremath{g_\text{Na}}}
\newcommand{\gbNa}{\ensuremath{\overline{g}_\text{Na}}}
\newcommand{\gtNa}{\ensuremath{\tilde{g}_\text{Na}}}

\newcommand{\IK}{\ensuremath{I_{\text{K}}}}
\newcommand{\VK}{\ensuremath{V_{\text{K}}}}
\newcommand{\gK}{\ensuremath{g_\text{K}}}
\newcommand{\gbK}{\ensuremath{\overline{g}_\text{K}}}
\newcommand{\gtK}{\ensuremath{\tilde{g}_\text{K}}}

\newcommand{\ICl}{\ensuremath{I_{\text{Cl}}}}
\newcommand{\VCl}{\ensuremath{V_{\text{Cl}}}}
\newcommand{\gCl}{\ensuremath{g_\text{Cl}}}
\newcommand{\gbCl}{\ensuremath{\overline{g}_\text{Cl}}}
\newcommand{\gtCl}{\ensuremath{\tilde{g}_\text{Cl}}}

\newcommand{\Ih}{\ensuremath{I_h}}
\newcommand{\Vh}{\ensuremath{V_h}}
\newcommand{\gh}{\ensuremath{g_h}}
\newcommand{\gbh}{\ensuremath{\overline{g}_h}}
\newcommand{\gth}{\ensuremath{\tilde{g}_h}}

\newcommand{\Iss}{\ensuremath{I^{ss}}}
\newcommand{\Istim}{\ensuremath{I_\text{stim}}}
\newcommand{\Vr}{\ensuremath{V_r}}
\newcommand{\Iij}{\ensuremath{\mathtt{I}_i^j}}
\newcommand{\Vij}{\ensuremath{\mathtt{V}_i^j}}
\newcommand{\Ii}{\ensuremath{\mathtt{V}_i}}
\newcommand{\Vj}{\ensuremath{\mathtt{V}^j}}
\newcommand{\Vtt}{\ensuremath{\mathtt{V}}}
\newcommand{\Vt}{\ensuremath{\tilde{V}}}

\newcommand{\dx}{\ensuremath{\text{d}x}}
\newcommand{\dt}{\ensuremath{\text{d}t}}
\newcommand{\bB}{\ensuremath{\mathbf{B}}}
\newcommand{\bS}{\ensuremath{\mathbf{S}}}

\newcommand{\va}{{\mathbf{a}}}
\newcommand{\vb}{{\mathbf{b}}}
\newcommand{\vc}{{\mathbf{c}}}
\newcommand{\vd}{{\mathbf{d}}}
\newcommand{\ve}{{\mathbf{e}}}
\newcommand{\vf}{{\mathbf{f}}}
\newcommand{\vg}{{\mathbf{g}}}
\newcommand{\vh}{{\mathbf{h}}}
\newcommand{\vi}{{\mathbf{i}}}
\newcommand{\vj}{{\mathbf{j}}}
\newcommand{\vk}{{\mathbf{k}}}
\newcommand{\vl}{{\mathbf{l}}}
\newcommand{\vm}{{\mathbf{m}}}
\newcommand{\vn}{{\mathbf{n}}}
\newcommand{\vo}{{\mathbf{o}}}
\newcommand{\vp}{{\mathbf{p}}}
\newcommand{\vq}{{\mathbf{q}}}
\newcommand{\vr}{{\mathbf{r}}}
\newcommand{\vs}{{\mathbf{s}}}
\newcommand{\vt}{{\mathbf{t}}}
\newcommand{\vu}{{\mathbf{u}}}
\newcommand{\vv}{{\mathbf{v}}}
\newcommand{\vw}{{\mathbf{w}}}
\newcommand{\vx}{{\mathbf{x}}}
\newcommand{\vy}{{\mathbf{y}}}
\newcommand{\vz}{{\mathbf{z}}}

\newcommand{\vA}{{\mathbf{A}}}
\newcommand{\vB}{{\mathbf{B}}}
\newcommand{\vC}{{\mathbf{C}}}
\newcommand{\vD}{{\mathbf{D}}}
\newcommand{\vE}{{\mathbf{E}}}
\newcommand{\vF}{{\mathbf{F}}}
\newcommand{\vG}{{\mathbf{G}}}
\newcommand{\vH}{{\mathbf{H}}}
\newcommand{\vI}{{\mathbf{I}}}
\newcommand{\vJ}{{\mathbf{J}}}
\newcommand{\vK}{{\mathbf{K}}}
\newcommand{\vL}{{\mathbf{L}}}
\newcommand{\vM}{{\mathbf{M}}}
\newcommand{\vN}{{\mathbf{N}}}
\newcommand{\vO}{{\mathbf{O}}}
\newcommand{\vP}{{\mathbf{P}}}
\newcommand{\vQ}{{\mathbf{Q}}}
\newcommand{\vR}{{\mathbf{R}}}
\newcommand{\vS}{{\mathbf{S}}}
\newcommand{\vT}{{\mathbf{T}}}
\newcommand{\vU}{{\mathbf{U}}}
\newcommand{\vV}{{\mathbf{V}}}
\newcommand{\vW}{{\mathbf{W}}}
\newcommand{\vX}{{\mathbf{X}}}
\newcommand{\vY}{{\mathbf{Y}}}
\newcommand{\vZ}{{\mathbf{Z}}}

\newcommand{\cA}{{\mathcal{A}}}
\newcommand{\cB}{{\mathcal{B}}}
\newcommand{\cC}{{\mathcal{C}}}
\newcommand{\cD}{{\mathcal{D}}}
\newcommand{\cE}{{\mathcal{E}}}
\newcommand{\cF}{{\mathcal{F}}}
\newcommand{\cG}{{\mathcal{G}}}
\newcommand{\cH}{{\mathcal{H}}}
\newcommand{\cI}{{\mathcal{I}}}
\newcommand{\cJ}{{\mathcal{J}}}
\newcommand{\cK}{{\mathcal{K}}}
\newcommand{\cL}{{\mathcal{L}}}
\newcommand{\cM}{{\mathcal{M}}}
\newcommand{\cN}{{\mathcal{N}}}
\newcommand{\cO}{{\mathcal{O}}}
\newcommand{\cP}{{\mathcal{P}}}
\newcommand{\cQ}{{\mathcal{Q}}}
\newcommand{\cR}{{\mathcal{R}}}
\newcommand{\cS}{{\mathcal{S}}}
\newcommand{\cT}{{\mathcal{T}}}
\newcommand{\cU}{{\mathcal{U}}}
\newcommand{\cV}{{\mathcal{V}}}
\newcommand{\cW}{{\mathcal{W}}}
\newcommand{\cX}{{\mathcal{X}}}
\newcommand{\cY}{{\mathcal{Y}}}
\newcommand{\cZ}{{\mathcal{Z}}}

\newcommand{\ri}{{\mathrm{i}}}
\newcommand{\rr}{{\mathrm{r}}}

\newcommand{\op}{{\mathrm{op}}} % subscript for operator norm
\newcommand{\opt}{{\mathrm{opt}}} % subscript for optimal solution
\newcommand{\supp}{{\mathrm{supp}}} % support
\newcommand{\sgn}{{\mathrm{sgn}}} % support
\newcommand{\Prob}{{\mathrm{Prob}}} % probability
\newcommand{\Diag}{{\mathrm{Diag}}} % vector -> diagonal matrix
\newcommand{\diag}{{\mathrm{diag}}} % matrix diagonal -> vector
\newcommand{\dom}{{\mathrm{dom}}} % domain
\newcommand{\grad}{{\nabla}}    % gradient
\newcommand{\tr}{{\mathrm{tr}}} % trace
\newcommand{\TV}{{\mathrm{TV}}} % total variation
\newcommand{\Proj}{{\mathrm{Proj}}}
\newcommand{\prox}{{\mathbf{prox}}}
\newcommand{\refl}{{\mathbf{refl}}}
\DeclareMathOperator{\shrink}{shrink} % shrinkage

\newcommand{\fF}{{\mathcal{F}}}
\newcommand{\fR}{{\mathcal{R}}}
\newcommand{\fL}{{\mathcal{L}}}

\newcommand{\1}{\ensuremath{\mathds{1}}}		% characteristic function

\DeclareMathOperator*{\esssup}{ess\,sup}
\DeclareMathOperator*{\mysin}{mysin}
\DeclareMathOperator*{\myasin}{myasin}
\DeclareMathOperator*{\pr}{pr}

asyn_theory_MP.bbl
\begin{thebibliography}{10}
\providecommand{\url}[1]{{#1}}
\providecommand{\urlprefix}{URL }
\expandafter\ifx\csname urlstyle\endcsname\relax
  \providecommand{\doi}[1]{DOI~\discretionary{}{}{}#1}\else
  \providecommand{\doi}{DOI~\discretionary{}{}{}\begingroup
  \urlstyle{rm}\Url}\fi

\bibitem{aharoni2000parallel}
Aharoni, D., Barak, A.: Parallel iterative discontinuous galerkin
  finite-element methods.
\newblock In: Discontinuous Galerkin Methods, pp. 247--254. Springer (2000)

\bibitem{AAI1998implicit}
Amitai, D., Averbuch, A., Israeli, M., Itzikowitz, S.: Implicit-explicit
  parallel asynchronous solver of parabolic pdes.
\newblock SIAM Journal on Scientific Computing \textbf{19}(4), 1366--1404
  (1998)

\bibitem{avron2014revisiting}
Avron, H., Druinsky, A., Gupta, A.: Revisiting asynchronous linear solvers:
  Provable convergence rate through randomization.
\newblock In: Parallel and Distributed Processing Symposium, 2014 IEEE 28th
  International, pp. 198--207. IEEE (2014)

\bibitem{BMR1997asyn-multisplit}
Bahi, J., Miellou, J.C., Rhofir, K.: Asynchronous multisplitting methods for
  nonlinear fixed point problems.
\newblock Numerical Algorithms \textbf{15}(3-4), 315--345 (1997)

\bibitem{baudet1978asynchronous}
Baudet, G.M.: Asynchronous iterative methods for multiprocessors.
\newblock Journal of the ACM (JACM) \textbf{25}(2), 226--244 (1978)

\bibitem{bauschke2011convex}
Bauschke, H.H., Combettes, P.L.: Convex analysis and monotone operator theory
  in Hilbert spaces.
\newblock Springer Science \& Business Media (2011)

\bibitem{bertsekas1983distributed}
Bertsekas, D.P.: Distributed asynchronous computation of fixed points.
\newblock Mathematical Programming \textbf{27}(1), 107--120 (1983)

\bibitem{bertsekas1989parallel}
Bertsekas, D.P., Tsitsiklis, J.N.: Parallel and distributed computation:
  numerical methods, vol.~23.
\newblock Prentice hall Englewood Cliffs, NJ (1989)

\bibitem{bertsekas1991some}
Bertsekas, D.P., Tsitsiklis, J.N.: Some aspects of parallel and distributed
  iterative algorithms—a survey.
\newblock Automatica \textbf{27}(1), 3--21 (1991)

\bibitem{bethune2014performance}
Bethune, I., Bull, J.M., Dingle, N.J., Higham, N.J.: Performance analysis of
  asynchronous jacobi's method implemented in mpi, shmem and openmp.
\newblock International Journal of High Performance Computing Applications
  \textbf{28}(1), 97--111 (2014)

\bibitem{boyd2011distributed}
Boyd, S., Parikh, N., Chu, E., Peleato, B., Eckstein, J.: Distributed
  optimization and statistical learning via the alternating direction method of
  multipliers.
\newblock Foundations and Trends{\textregistered} in Machine Learning
  \textbf{3}(1), 1--122 (2011)

\bibitem{chang2008coordinate}
Chang, K.W., Hsieh, C.J., Lin, C.J.: Coordinate descent method for large-scale
  l2-loss linear support vector machines.
\newblock The Journal of Machine Learning Research \textbf{9}, 1369--1398
  (2008)

\bibitem{Chau20081126}
Chau, M., Spiteri, P., Guivarch, R., Boisson, H.: Parallel asynchronous
  iterations for the solution of a 3d continuous flow electrophoresis problem.
\newblock Computers \& Fluids \textbf{37}(9), 1126 -- 1137 (2008).
\newblock \doi{http://dx.doi.org/10.1016/j.compfluid.2007.06.006}.
\newblock
  \urlprefix\url{http://www.sciencedirect.com/science/article/pii/S0045793007001995}

\bibitem{chazan1969chaotic}
Chazan, D., Miranker, W.: Chaotic relaxation.
\newblock Linear algebra and its applications \textbf{2}(2), 199--222 (1969)

\bibitem{combettes2014stochastic}
Combettes, P.L., Pesquet, J.C.: Stochastic quasi-fej{\'e}r block-coordinate
  fixed point iterations with random sweeping.
\newblock arXiv preprint arXiv:1404.7536  (2014)

\bibitem{condat2013primal}
Condat, L.: A primal--dual splitting method for convex optimization involving
  lipschitzian, proximable and linear composite terms.
\newblock Journal of Optimization Theory and Applications \textbf{158}(2),
  460--479 (2013)

\bibitem{DavisYin2014}
Davis, D., Yin, W.: Convergence rate analysis of several splitting schemes.
\newblock {arXiv} preprint {arXiv}:1406.4834  (2014)

\bibitem{DavisYin2014c}
Davis, D., Yin, W.: Faster convergence rates of relaxed peaceman-rachford and
  {ADMM} under regularity assumptions.
\newblock {arXiv} preprint {arXiv}:1407.5210  (2014)

\bibitem{davis2015three}
Davis, D., Yin, W.: A three-operator splitting scheme and its optimization
  applications.
\newblock arXiv preprint arXiv:1504.01032  (2015)

\bibitem{donzis2014asynchronous}
Donzis, D.A., Aditya, K.: Asynchronous finite-difference schemes for partial
  differential equations.
\newblock Journal of Computational Physics \textbf{274}, 370--392 (2014)

\bibitem{ESMG1996asyn-flex}
El~Baz, D., Spiteri, P., Miellou, J.C., Gazen, D.: Asynchronous iterative
  algorithms with flexible communication for nonlinear network flow problems.
\newblock Journal of Parallel and Distributed Computing \textbf{38}(1), 1--15
  (1996)

\bibitem{el1982some}
El~Tarazi, M.N.: Some convergence results for asynchronous algorithms.
\newblock Numerische Mathematik \textbf{39}(3), 325--340 (1982)

\bibitem{leifang_information_2005}
Fang, L., Antsaklis, P.J.: Information consensus of asynchronous discrete-time
  multi-agent systems.
\newblock In: American Control Conference, 2005. Proceedings of the 2005, pp.
  1883--1888. IEEE (2005)

\bibitem{FSS1997asyn-addSch}
Frommer, A., Schwandt, H., Szyld, D.B.: Asynchronous weighted additive schwarz
  methods.
\newblock Electronic Transactions on Numerical Analysis \textbf{5}, 48--61
  (1997)

\bibitem{frommer2000asynchronous}
Frommer, A., Szyld, D.B.: On asynchronous iterations.
\newblock Journal of computational and applied mathematics \textbf{123}(1),
  201--216 (2000)

\bibitem{gabay1983chapter}
Gabay, D.: Chapter ix applications of the method of multipliers to variational
  inequalities.
\newblock Studies in mathematics and its applications \textbf{15}, 299--331
  (1983)

\bibitem{glowinski1975approximation}
Glowinski, R., Marroco, A.: Sur l'approximation, par éléments finis d'ordre
  un, et la résolution, par pénalisation-dualité d'une classe de
  problèmes de dirichlet non linéaires.
\newblock ESAIM: Mathematical Modelling and Numerical Analysis-Mod{\'e}lisation
  Math{\'e}matique et Analyse Num{\'e}rique \textbf{9}(R2), 41--76 (1975)

\bibitem{hong2014distributed}
Hong, M.: A distributed, asynchronous and incremental algorithm for nonconvex
  optimization: An admm based approach.
\newblock arXiv preprint arXiv:1412.6058  (2014)

\bibitem{house2014big}
House, W.: Big data: Seizing opportunities, preserving values (2014)

\bibitem{hsieh2015passcode}
Hsieh, C.J., Yu, H.F., Dhillon, I.S.: Passcode: Parallel asynchronous
  stochastic dual co-ordinate descent.
\newblock arXiv preprint arXiv:1504.01365  (2015)

\bibitem{iutzeler2013asynchronous}
Iutzeler, F., Bianchi, P., Ciblat, P., Hachem, W.: Asynchronous distributed
  optimization using a randomized alternating direction method of multipliers.
\newblock In: Decision and Control (CDC), 2013 IEEE 52nd Annual Conference on,
  pp. 3671--3676. IEEE (2013)

\bibitem{krasnosel1955two}
Krasnosel'skii, M.A.: Two remarks on the method of successive approximations.
\newblock Uspekhi Matematicheskikh Nauk \textbf{10}(1), 123--127 (1955)

\bibitem{kyrola2011parallel}
Kyrola, A., Bickson, D., Guestrin, C., Bradley, J.K.: Parallel coordinate
  descent for l1-regularized loss minimization.
\newblock In: Proceedings of the 28th International Conference on Machine
  Learning (ICML-11), pp. 321--328 (2011)

\bibitem{LMS1986asynchronous}
Lang, B., Miellou, J., Spiteri, P.: Asynchronous relaxation algorithms for
  optimal control problems.
\newblock Mathematics and computers in simulation \textbf{28}(3), 227--242
  (1986)

\bibitem{larson1981urban}
Larson, R.C., Odoni, A.R.: Urban operations research.
\newblock Monograph (1981)

\bibitem{lions1979splitting}
Lions, P.L., Mercier, B.: Splitting algorithms for the sum of two nonlinear
  operators.
\newblock SIAM Journal on Numerical Analysis \textbf{16}(6), 964--979 (1979)

\bibitem{liu2014asynchronous}
Liu, J., Wright, S.J.: Asynchronous stochastic coordinate descent: Parallelism
  and convergence properties.
\newblock SIAM Journal on Optimization \textbf{25}(1), 351--376 (2015)

\bibitem{liu2013asynchronous}
Liu, J., Wright, S.J., R{\'e}, C., Bittorf, V., Sridhar, S.: An asynchronous
  parallel stochastic coordinate descent algorithm.
\newblock Journal of Machine Learning Research \textbf{16}, 285--322 (2015)

\bibitem{nedic2001distributed}
Nedi{\'c}, A., Bertsekas, D.P., Borkar, V.S.: Distributed asynchronous
  incremental subgradient methods.
\newblock Studies in Computational Mathematics \textbf{8}, 381--407 (2001)

\bibitem{nedic2009distributed}
Nedic, A., Ozdaglar, A.: Distributed subgradient methods for multi-agent
  optimization.
\newblock Automatic Control, IEEE Transactions on \textbf{54}(1), 48--61 (2009)

\bibitem{nesterov2012rcd}
Nesterov, Y.: Efficiency of coordinate descent methods on huge-scale
  optimization problems.
\newblock SIAM Journal on Optimization \textbf{22}(2), 341--362 (2012)

\bibitem{opial_weak_1967}
{Opial}, Z.: Weak convergence of the sequence of successive approximations for
  nonexpansive mappings.
\newblock Bulletin of the American Mathematical Society \textbf{73}(4),
  591--598 (1967).
\newblock \doi{10.1090/S0002-9904-1967-11761-0}

\bibitem{passty1979ergodic}
Passty, G.B.: Ergodic convergence to a zero of the sum of monotone operators in
  hilbert space.
\newblock Journal of Mathematical Analysis and Applications \textbf{72}(2),
  383--390 (1979)

\bibitem{peng2016coordinate}
Peng, Z., Wu, T., Xu, Y., Yan, M., Yin, W.: Coordinate friendly structures,
  algorithms and applications.
\newblock arXiv preprint arXiv:1601.00863  (2016)

\bibitem{peng2013parallel}
Peng, Z., Yan, M., Yin, W.: Parallel and distributed sparse optimization.
\newblock In: Signals, Systems and Computers, 2013 Asilomar Conference on, pp.
  659--646. IEEE (2013)

\bibitem{petryshyn1966construction}
Petryshyn, W.: Construction of fixed points of demicompact mappings in hilbert
  space.
\newblock Journal of Mathematical Analysis and Applications \textbf{14}(2),
  276--284 (1966)

\bibitem{peypouquet_evolution_2009}
Peypouquet, J., Sorin, S.: Evolution equations for maximal monotone operators:
  Asymptotic analysis in continuous and discrete time.
\newblock Journal of Convex Analysis \textbf{17}, 1113--1163 (2010)

\bibitem{recht2011hogwild}
Recht, B., Re, C., Wright, S., Niu, F.: Hogwild: A lock-free approach to
  parallelizing stochastic gradient descent.
\newblock In: Advances in Neural Information Processing Systems, pp. 693--701
  (2011)

\bibitem{richtarik2015parallel}
Richt{\'a}rik, P., Tak{\'a}{\v{c}}, M.: Parallel coordinate descent methods for
  big data optimization.
\newblock Mathematical Programming, Series A pp. 1--52 (2015)

\bibitem{robbins1985convergence}
Robbins, H., Siegmund, D.: A convergence theorem for non negative almost
  supermartingales and some applications.
\newblock In: Herbert Robbins Selected Papers, pp. 111--135. Springer (1985)

\bibitem{rosenfeld1969case}
Rosenfeld, J.L.: A case study in programming for parallel-processors.
\newblock Communications of the ACM \textbf{12}(12), 645--655 (1969)

\bibitem{tai2002convergence}
Tai, X.C., Tseng, P.: Convergence rate analysis of an asynchronous space
  decomposition method for convex minimization.
\newblock Mathematics of Computation \textbf{71}(239), 1105--1135 (2002)

\bibitem{tseng1991rate-asyn}
Tseng, P.: On the rate of convergence of a partially asynchronous gradient
  projection algorithm.
\newblock SIAM Journal on Optimization \textbf{1}(4), 603--619 (1991)

\bibitem{TB1990partially}
Tseng, P., Bertsekas, D.P., Tsitsiklis, J.N.: Partially asynchronous, parallel
  algorithms for network flow and other problems.
\newblock SIAM Journal on Control and Optimization \textbf{28}(3), 678--710
  (1990)

\bibitem{vu2013splitting}
V{\~u}, B.C.: A splitting algorithm for dual monotone inclusions involving
  cocoercive operators.
\newblock Advances in Computational Mathematics \textbf{38}(3), 667--681 (2013)

\bibitem{wei2013on}
Wei, E., Ozdaglar, A.: On the o(1/k) convergence of asynchronous distributed
  alternating direction method of multipliers.
\newblock In: Global Conference on Signal and Information Processing
  (GlobalSIP), 2013 IEEE, pp. 551--554. IEEE (2013)

\bibitem{YanYin2014}
Yan, M., Yin, W.: Self equivalence of the alternating direction method of
  multipliers.
\newblock {arXiv} preprint {arXiv}:1407.7400  (2014)

\bibitem{yuan2013convergence}
Yuan, K., Ling, Q., Yin, W.: On the convergence of decentralized gradient
  descent.
\newblock arXiv preprint arXiv:1310.7063  (2013)

\bibitem{zhang2014asynchronous}
Zhang, R., Kwok, J.: Asynchronous distributed admm for consensus optimization.
\newblock In: Proceedings of the 31st International Conference on Machine
  Learning (ICML-14), pp. 1701--1709 (2014)

\end{thebibliography}
