\section{Auxiliary Results}
\label{app:first}

We collect here in one place all the auxiliary results in the existing literature that will
be used in our proofs subsequently. 
The first one is a well-known characterization of convex sets and the projection operator given in~\cite{Nes04}:

\begin{lemma}\label{lem:aux1}\quad
	Let $\mathcal{\debug X}$ be a compact and convex subset of $\mathbb{R}^d$. Then for any $\debug x \in \mathcal{\debug X}, \debug y \in \mathbb{R}^d$:
	\begin{equation}
	\inner{\mirror(\debug y) - \debug x}{\mirror(\debug y) -\debug y}
	\leq 0.
	\end{equation}
	
\end{lemma}

The second one is an $L_{p}$-bounded martingale convergence theorem given in~\cite{HH80}:

\begin{lemma} \quad
	\label{thm:mg_convergence2}
	Let $S_n$ be a martingale adapted to the filtration $\mathcal{S}_n$.
	If for some $p \ge 1$, $\sup_{n \ge 0} \mathbf{E}[|S_n|^p] < \infty$, then $S_n$ converges almost surely to a limiting random variable $S_{\infty}$ with $\mathbf{E}[|S_{\infty}|^p] < \infty$.
\end{lemma}

\begin{remark}\quad
	Note that $\mathbf{E}[|S_{\infty}|^p] < \infty \text{ for } p \ge 1$ obviously implies $S_{\infty}$ is finite almost surely.
\end{remark}

The third one is the classical envelope theorem (see~\cite{carter2001foundations}).

\begin{lemma}\label{lem:env}\quad
	Let $f: \mathbb{R}^n \times \mathbb{R}^m \rightarrow \mathbb{R}$ be a continuously differentiable function.
	Let $U$ be a compact set and consider the problem 
	\begin{equation}\label{eq:env_opt}
	\max_{x \in U} f(x, \theta).
	\end{equation}
	Let $x^*: \mathcal{O} \rightarrow \mathbb{R}^m$ be a continuous function defined on an open set $\mathcal{O} \subset \mathbb{R}^m$ such that for each $\theta \in \mathcal{O}$,
	$x^*(\theta)$ solves the problem in Equation~\ref{eq:env_opt}.
	Define $V: \mathbb{R}^m \rightarrow \mathbb{R}$ where $V(\theta) = f(x^*(\theta), \theta)$.
	Then $V(\theta)$ is differentiable on $\mathcal{O}$ and:
	\begin{equation}
	\grad V(\theta) = \grad f(x^*(\theta), \theta).
	\end{equation}
\end{lemma}

The fourth one is an elementary sequence result (see~\cite{bertsekas1995nonlinear}).
\begin{lemma}\label{lem:A3}\quad
	Let $a_n, b_n$ be two non-negative sequences such that 
	$\sum_{n=1}^\infty a_n =\infty , \sum_{n=1}^\infty 
	a_nb_n < \infty$. 
	If there exists a real number $K > 0$ such that $|b_{n+1} - b_n| \leq Ka_n$. Then, $\lim_{n \rightarrow \infty} b_n =0$.
\end{lemma}	

The fifth one is law of large numbers for martingales given in~\cite{HH80}:
\begin{lemma}
	\label{thm:mg_convergence1}\quad
	Let $M_{\debug n} = \sum_{\debug k=0}^{\debug n} d_{\debug k}$ be a martingale adapted to $(\mathcal{\debug F}_{\debug n})_{\debug n=0}^{\infty}$ and let $(u_{\debug n})_{\debug n=0}^{\infty}$ be a nondecreasing sequence of positive numbers with $\lim_{\debug n\to\infty} u_{\debug n} = \infty$.
	If $\sum_{\debug n=0}^{\infty} u_{\debug n}^{-p} \exof{\abs{d_{\debug k}}^{p}\given\mathcal{\debug F}_{\debug n}} <\infty$ for some $p\in[1,2]$ \textup(a.s.\textup)\xspace, then:
	\begin{equation}
	\label{eq:LLN}
	\lim_{\debug n\to\infty} \frac{M_{\debug n}}{u_{\debug n}} = 0
	\quad
	\textup(a.s.\textup)\xspace
	\end{equation}
\end{lemma}

Finally, we recall the standard notion of semiflow. 
\begin{definition}\quad
	A semiflow $\debug P$ on a metric space $(M, d)$ is a continuous map $\debug P: \mathbb{R}_+ \times M \to M$:
	$$(t, x) \to \debug P_t(x),$$
	such that the semi-group properties hold: $\debug P_0$ = identity, $\debug P_{t+s} = \debug P_t \circ \debug P_s$ for all $(t,s) \in \mathbb{R}_+ \times \mathbb{R}_+$.
\end{definition}

\begin{remark}\label{rem:semi_flow}\quad
	A standard way to induce a semiflow is via an \ac{ODE}.
	Specifically, if $F: \mathbb{R}^m \to \mathbb{R}^m$ is a continuous function and if the following \ac{ODE} has a unique solution trajectory for each initial point $\tilde{x} \in \mathbb{R}^m$:
	\begin{eqnarray*}
		\frac{dx}{dt}=&F(x),\\
		x(0)&=\tilde{x},
	\end{eqnarray*}
	then $\debug P_t(\tilde{x})$ defined by the solution trajectory $x(t) \in \mathbb{R}^m$ as follows is a semiflow: $\debug P_t(\tilde{x}) \triangleq x(t)$ with $x(0) = \tilde{x}$. We say $\debug P$ defined in this way is the semiflow induced by the corresponding \ac{ODE}.
\end{remark}

\section{Technical Proofs}

\subsection{General Nonconvex Objectives}

\subsubsection{Proof of Lemma~\ref{lem:tail}}
{\em Proof:}
We start by rewriting the delayed gradient update $\debug X_{\debug n+1} = \debug X_{\debug n} - \debug\alpha_{\debug n+1} \nabla \debug F(\debug X_{s(\debug n)},\debug\omega_{\debug n+1})$ in Algorithm~\ref{alg:DASGD_un} in two forms as follows, both of which will be used subsequently:
\begin{equation}
\begin{split}
&\debug X_{\debug n+1} = \debug X_{\debug n} - \debug\alpha_{\debug n+1}\Bigg(\grad \debug f(\debug X_{s(\debug n)}) + \grad \debug F(\debug X_{s(\debug n)},\debug\omega_{\debug n+1}) -  \grad \debug f(\debug X_{s(\debug n)}) \Bigg)\\
& =\debug X_{\debug n} - \debug\alpha_{\debug n+1}\Bigg( \grad \debug f(\debug X_{\debug n}) + \grad \debug f(\debug X_{s(\debug n)}) - \grad \debug f(\debug X_{\debug n})  + \grad \debug F(\debug X_{s(\debug n)},\debug\omega_{\debug n+1}) -  \grad \debug f(\debug X_{s(\debug n)}) \Bigg).
\end{split}
\end{equation}

Denoted $B \triangleq \sup_{\debug x \in \mathcal{\debug X}} \exof{\|\nabla \debug F(x,\debug\omega)\|_2^2}$ per Assumption~\ref{asm:basic}.
Since $\grad f(x) = \exof{\grad \debug F(\debug x, \debug\omega)}$ is Lipschitz per Assumption~\ref{asm:basic}, letting $L$ be the Lipschitz constant, we have:
\begin{equation}\label{eq:Lip}
\begin{split}
&\debug f(\debug X_{\debug n+1}) - \debug f(\debug X_{\debug n}) 
\leq  \langle  \grad \debug f(\debug X_{\debug n}),  \debug X_{\debug n+1} - \debug X_{\debug n}\rangle + \frac{L}{2}\|\debug X_{\debug n+1} - \debug X_{\debug n} \|_2^2\\
& = - \debug\alpha_{\debug n+1}\langle  \grad \debug f(\debug X_{\debug n}),   \nabla \debug F(\debug X_{s(\debug n)},\debug\omega_{\debug n+1})\rangle + \frac{L}{2}\|\debug\alpha_{\debug n+1} \nabla \debug F(\debug X_{s(\debug n)},\debug\omega_{\debug n+1}) \|_2^2\\
& = - \debug\alpha_{\debug n+1}\langle  \grad \debug f(\debug X_{\debug n}),  \grad \debug f(\debug X_{\debug n}) + \grad \debug f(\debug X_{s(\debug n)}) - \grad \debug f(\debug X_{\debug n})  + \grad \debug F(\debug X_{s(\debug n)},\debug\omega_{\debug n+1}) -  \grad \debug f(\debug X_{s(\debug n)}) \rangle \\
&+ \frac{L}{2}\debug\alpha_{\debug n+1}^2\|  \grad \debug f(\debug X_{s(\debug n)})  + \grad \debug F(\debug X_{s(\debug n)},\debug\omega_{\debug n+1}) -  \grad \debug f(\debug X_{s(\debug n)})  \|_2^2\\
& = - \debug\alpha_{\debug n+1} \|\grad \debug f(\debug X_{\debug n})\|_2^2 -  \debug\alpha_{\debug n+1}\langle  \grad \debug f(\debug X_{\debug n}),   \grad \debug f(\debug X_{s(\debug n)}) - \grad \debug f(\debug X_{\debug n})\rangle
	\notag\\
	&-  \debug\alpha_{\debug n+1}\langle  \grad \debug f(\debug X_{\debug n}), \grad \debug F(\debug X_{s(\debug n)},\debug\omega_{\debug n+1}) -  \grad \debug f(\debug X_{s(\debug n)}) \rangle \\
&+ \frac{L}{2}\debug\alpha_{\debug n+1}^2\|  \grad \debug f(\debug X_{s(\debug n)})  + \grad \debug F(\debug X_{s(\debug n)},\debug\omega_{\debug n+1}) -  \grad \debug f(\debug X_{s(\debug n)})  \|_2^2\\
& \leq - \debug\alpha_{\debug n+1} \|\grad \debug f(\debug X_{\debug n})\|_2^2 +  \debug\alpha_{\debug n+1}\| \grad \debug f(\debug X_{\debug n})\|_2  \| \grad \debug f(\debug X_{s(\debug n)}) - \grad \debug f(\debug X_{\debug n})\|_2 
	\notag\\
	&- \debug\alpha_{\debug n+1}\langle  \grad \debug f(\debug X_{\debug n}), \grad \debug F(\debug X_{s(\debug n)},\debug\omega_{\debug n+1}) -  \grad \debug f(\debug X_{s(\debug n)}) \rangle \\
&+ \frac{L}{2}\debug\alpha_{\debug n+1}^2 \Big\{2 \|\grad \debug f(\debug X_{s(\debug n)})\|_2^2  + 2\|\grad \debug F(\debug X_{s(\debug n)},\debug\omega_{\debug n+1}) -  \grad \debug f(\debug X_{s(\debug n)})  \|_2^2\Big\}\\
& \leq - \debug\alpha_{\debug n+1} \|\grad \debug f(\debug X_{\debug n})\|_2^2 +  \sqrt{B}\debug\alpha_{\debug n+1}\| \grad \debug f(\debug X_{s(\debug n)}) - \grad \debug f(\debug X_{\debug n})\|_2
	\notag\\
	&-  \debug\alpha_{\debug n+1}\langle  \grad \debug f(\debug X_{\debug n}), \grad \debug F(\debug X_{s(\debug n)},\debug\omega_{\debug n+1}) -  \grad \debug f(\debug X_{s(\debug n)}) \rangle \\
& \leq L\debug\alpha_{\debug n+1}^2 \Big\{ \sqrt{B}  + \|\grad \debug F(\debug X_{s(\debug n)},\debug\omega_{\debug n+1}) -  \grad \debug f(\debug X_{s(\debug n)})  \|_2^2\Big\},\\
\end{split}
\end{equation}
where in the last inequality follows because by Jensen's inequality, we have:
$$\| \grad \debug f(\debug X_{\debug n})\|_2^2   \leq \sup_{\debug x \in \mathcal{\debug X}} \|\exof{\nabla \debug F(x,\debug\omega)}\|_2^2  \leq \sup_{\debug x \in \mathcal{\debug X}} \exof{\|\nabla \debug F(x,\debug\omega)\|_2^2} \leq B.$$ 

Denote the filtration generated by $\debug X_0, \debug X_1, \dots, \debug X_{\debug n}$ to be
$\mathcal{F}_{\debug n}$. We take the expectation of both sides of Equation~\eqref{eq:Lip} and obtain:
\begin{equation}\label{eq:Lip_ex}
\begin{split}
&\exof{\debug f(\debug X_{\debug n+1}) - \debug f(\debug X_{\debug n})}
\leq - \debug\alpha_{\debug n+1} \exof{\|\grad \debug f(\debug X_{\debug n})\|_2^2} +  \sqrt{B}\debug\alpha_{\debug n+1} \exof{\| \grad \debug f(\debug X_{s(\debug n)}) - \grad \debug f(\debug X_{\debug n})\|_2} \\
& -  \debug\alpha_{\debug n+1} \exof{\langle  \grad \debug f(\debug X_{\debug n}), \grad \debug F(\debug X_{s(\debug n)},\debug\omega_{\debug n+1}) -  \grad \debug f(\debug X_{s(\debug n)}) \rangle} +  L\sqrt{B}\debug\alpha_{\debug n+1}^2
	\notag\\
	&+ L\debug\alpha_{\debug n+1}^2\exof{\|\grad \debug F(\debug X_{s(\debug n)},\debug\omega_{\debug n+1}) -  \grad \debug f(\debug X_{s(\debug n)})  \|_2^2}\\
& = - \debug\alpha_{\debug n+1} \exof{\|\grad \debug f(\debug X_{\debug n})\|_2^2} +  \sqrt{B}\debug\alpha_{\debug n+1} \exof{\| \grad \debug f(\debug X_{s(\debug n)}) - \grad \debug f(\debug X_{\debug n})\|_2}
	\notag\\
	&-  \debug\alpha_{\debug n+1} 
\ex \Bigg\{\ex \Big\{\langle  \grad \debug f(\debug X_{\debug n}), \grad \debug F(\debug X_{s(\debug n)},\debug\omega_{\debug n+1}) -  \grad \debug f(\debug X_{s(\debug n)}) \rangle \Big | \mathcal{F}_{\debug n} \Big\} \Bigg\}\\
&+  L\sqrt{B}\debug\alpha_{\debug n+1}^2   
+ L\debug\alpha_{\debug n+1}^2\exof{\|\grad \debug F(\debug X_{s(\debug n)},\debug\omega_{\debug n+1}) -  \grad \debug f(\debug X_{s(\debug n)})  \|_2^2}\\
& =- \debug\alpha_{\debug n+1} \exof{\|\grad \debug f(\debug X_{\debug n})\|_2^2} +  \sqrt{B}\debug\alpha_{\debug n+1} \exof{\| \grad \debug f(\debug X_{s(\debug n)}) - \grad \debug f(\debug X_{\debug n})\|_2}
	\notag\\
	&-  \debug\alpha_{\debug n+1} 
\ex \Bigg\{\langle  \grad \debug f(\debug X_{\debug n}), \ex \Big\{\grad \debug F(\debug X_{s(\debug n)},\debug\omega_{\debug n+1}) -  \grad \debug f(\debug X_{s(\debug n)})  \Big | \mathcal{F}_{\debug n}\Big\}\rangle  \Bigg\}\\
&+  L\sqrt{B}\debug\alpha_{\debug n+1}^2   
+ L\debug\alpha_{\debug n+1}^2\exof{\|\grad \debug F(\debug X_{s(\debug n)},\debug\omega_{\debug n+1}) -  \grad \debug f(\debug X_{s(\debug n)})  \|_2^2}\\
& =- \debug\alpha_{\debug n+1} \exof{\|\grad \debug f(\debug X_{\debug n})\|_2^2} +  \sqrt{B}\debug\alpha_{\debug n+1} \exof{\| \grad \debug f(\debug X_{s(\debug n)}) - \grad \debug f(\debug X_{\debug n})\|_2}  
+  L\sqrt{B}\debug\alpha_{\debug n+1}^2
	\notag\\
	&+ L\debug\alpha_{\debug n+1}^2\exof{\|\grad \debug F(\debug X_{s(\debug n)},\debug\omega_{\debug n+1}) -  \grad \debug f(\debug X_{s(\debug n)})  \|_2^2}\\
& \leq - \debug\alpha_{\debug n+1} \exof{\|\grad \debug f(\debug X_{\debug n})\|_2^2} +  \sqrt{B}\debug\alpha_{\debug n+1} \exof{\| \grad \debug f(\debug X_{s(\debug n)}) - \grad \debug f(\debug X_{\debug n})\|_2}  
+  L\sqrt{B}\debug\alpha_{\debug n+1}^2   
+ 4LB\debug\alpha_{\debug n+1}^2,\\
\end{split}
\end{equation}
where the third equality follows from $\ex \Big\{\grad \debug F(\debug X_{s(\debug n)},\debug\omega_{\debug n+1}) -  \grad \debug f(\debug X_{s(\debug n)})  \Big | \mathcal{F}_{\debug n}\Big \} = 0$, since $\debug\omega_{\debug n +1}$ is independent of $\mathcal{F}_{\debug n}$, and
the last inequality follows from 
$\exof{\|\grad \debug F(\debug X_{s(\debug n)},\debug\omega_{\debug n+1}) -  \grad \debug f(\debug X_{s(\debug n)})  \|_2^2} \leq 
\sup_{\debug x \in \mathcal{\debug X}} \exof{\|\grad \debug F(\debug x,\debug\omega_{\debug n+1}) -  \grad \debug f(\debug x)  \|_2^2}  \leq 4B$.
Since $\grad \debug f$ is $L$-Liptchiz continuous, we have:
\begin{flalign}
\label{eq:expansion}
&\|\v(\debug X_{s(\debug n)}) - \v(\debug X_\debug n)\|_2 \leq L \|\debug X_{s(\debug n)} - \debug X_\debug n\|_2
	\notag\\
	&= L \Big\| \debug X_{s(\debug n)} - \debug X_{{s(\debug n)}+1} + \debug X_{{s(\debug n)}+1} - \debug X_{{s(\debug n)}+2} + \dots +  \debug X_{\debug n-1} - \debug X_{\debug n}  \Big\|_2
\notag\\
&= L\Big\|\sum_{r = s(\debug n)}^{\debug n-1} \Big\{\debug X_{r} - \debug X_{r+1} \Big\}\Big\|_2
= L\Big\|\sum_{r = s(\debug n)}^{\debug n-1} \debug\alpha_{r+1} \nabla \debug F(\debug X_{s(r)},\debug\omega_{r+1}) \Big\|_2
\notag \leq L \sum_{r = s(\debug n)}^{\debug n-1} \debug\alpha_{r+1} \Big\|\nabla \debug F(\debug X_{s(r)},\debug\omega_{r+1})\Big\|_2
\notag.\\
\end{flalign}
Taking the expectation of both sides of the above equation then yields:
\begin{flalign}
\label{eq:expansion1}
&\exof{\|\v(\debug X_{s(\debug n)}) - \v(\debug X_\debug n)\|_2}
\leq L \sum_{r = s(\debug n)}^{\debug n-1} \debug\alpha_{r+1} \exof{\|\nabla \debug F(\debug X_{s(r)},\debug\omega_{r+1})\|_2}
	\notag\\
	&\leq L \sum_{r = s(\debug n)}^{\debug n-1} \debug\alpha_{r+1} \sup_{\debug x \in \mathcal{\debug X}} \exof{\|\nabla \debug F(x,\debug\omega_{r+1})\|_2} \\
& \leq L\bar{B} \sum_{r = s(\debug n)}^{\debug n-1} \debug\alpha_{r+1},
\end{flalign}
where $\bar{B}  \triangleq \sup_{\debug x \in \mathcal{\debug X}} \exof{\|\nabla \debug F(x,\debug\omega_{r+1})\|_2} < \infty$ per Remark~\ref{rem:basic}.
Combining Equation~\eqref{eq:expansion1} with Equation~\eqref{eq:Lip_ex} then yields:
\begin{equation}\label{eq:Lip_ex2}
\begin{split}
&\exof{\debug f(\debug X_{\debug n+1})} - \exof{\debug f(\debug X_{\debug n})} \leq - \debug\alpha_{\debug n+1} \exof{\|\grad \debug f(\debug X_{\debug n})\|_2^2} +  L\bar{B}\sqrt{B}  \debug\alpha_{\debug n+1}\sum_{r = s(\debug n)}^{\debug n-1} \debug\alpha_{r+1}
+  L\sqrt{B}\debug\alpha_{\debug n+1}^2   
+ 4LB\debug\alpha_{\debug n+1}^2.\\
\end{split}
\end{equation}
Telescoping Equation~\eqref{eq:Lip_ex2} then yields:
\begin{equation}\label{eq:Lip_ex3}
\begin{split}
& -\infty < \inf_{\debug x \in \mathcal{\debug X}} \debug f(\debug x)- \debug f(\debug X_{0}) 
\le
\exof{\debug f(\debug X_{T+1})} - \debug f(\debug X_{0}) =
\sum_{\debug n=0}^T\Big\{\exof{\debug f(\debug X_{\debug n+1})} - \exof{\debug f(\debug X_{\debug n})} \Big\}\\
&\leq - \sum_{\debug n=0}^T\debug\alpha_{\debug n+1} \exof{\|\grad \debug f(\debug X_{\debug n})\|_2^2} +  L\bar{B}\sqrt{B} \sum_{\debug n=0}^T \Big(\debug\alpha_{\debug n+1}\sum_{r = s(\debug n)}^{\debug n-1} \debug\alpha_{r+1}\Big)
+  (L\sqrt{B}+ 4LB) \sum_{\debug n = 0}^T\debug\alpha_{\debug n+1}^2.\\
\end{split}
\end{equation}
Taking $T \to \infty$, the above equation yields:
\begin{equation}\label{eq:Lip_ex3}
\begin{split}
& -\infty < 
- \sum_{\debug n=0}^{\infty}\debug\alpha_{\debug n+1} \exof{\|\grad \debug f(\debug X_{\debug n})\|_2^2} +  L\bar{B}\sqrt{B} \sum_{\debug n=0}^{\infty} \Big(\debug\alpha_{\debug n+1}\sum_{r = s(\debug n)}^{\debug n-1} \debug\alpha_{r+1}\Big)
+  (L\sqrt{B}+ 4LB) \sum_{\debug n = 0}^{\infty}\debug\alpha_{\debug n+1}^2.\\
\end{split}
\end{equation}
Note that in all cases in Assumption~\ref{asm:delays}, we have $ \sum_{\debug n = 0}^{\infty}\debug\alpha_{\debug n+1}^2 < \infty$.
We next proceed to bound $\sum_{\debug n=0}^{\infty} \Big(\debug\alpha_{\debug n+1}\sum_{r = s(\debug n)}^{\debug n-1} \debug\alpha_{r+1}\Big)$ and show that
it is finite in each of the cases in Assumption~\ref{asm:delays}.


\begin{enumerate}
	\item In the bounded delays case, since $\sup_{\debug n} d_{\debug n} \leq D$, it follows that $s(n) + D \ge n$ and hence:
	\begin{equation}
	\begin{split}
	&\sum_{\debug n=0}^{\infty} \Big(\debug\alpha_{\debug n+1}\sum_{r = s(\debug n)}^{\debug n-1} \debug\alpha_{r+1}\Big) \leq 
	\sum_{\debug n=0}^{\infty} \Big(\debug\alpha_{\debug n+1}\sum_{r = \debug n - 1 -D}^{\debug n} \debug\alpha_{r+1}\Big) \leq 
	(D+1)\sum_{\debug n=0}^{\infty} \Big(\debug\alpha_{\debug n+1}  \max_{r \in \{\debug n - 1 -D, \dots, \debug n\}   } \debug\alpha_{r+1}\Big) \\
	&\le
	(D+1)\sum_{\debug n=0}^{\infty} \Big( \max_{r \in \{\debug n - 1 -D, \dots, \debug n\}   } \debug\alpha_{r+1} \cdot \max_{r \in \{\debug n - 1 -D, \dots, \debug n\}   } \debug\alpha_{r+1}\Big) = (D+1)\sum_{\debug n=0}^{\infty} \Big(\max_{r \in \{\debug n - 1 -D, \dots, \debug n\}   } \debug\alpha_{r+1}^2\Big)
	\notag\\
	&\leq  (D+1)\sum_{\debug n=0}^{\infty} \Big(\sum_{r =\debug n - 1 -D}^{\debug n} \debug\alpha_{r+1}^2\Big)\\
	& = (D+1)^2 \sum_{\debug n = 0}^{\infty}\debug\alpha_{\debug n+1}^2 < \infty,
	\end{split}
	\end{equation}
	where all the terms $\alpha_r$ are defined to be $0$ when $r$ drops below $0$.
	\item In the sublinearly growing delays case, since $d_{\debug n} = O(\debug n^p)$, it follows that $s(n) + K s^p(\debug n) \ge \debug n$ for some positive number $K$, which further implies that $s(n) + K s(\debug n) \ge s(n) + K s^p(\debug n) \ge \debug n$, thereby leading to
	$s(n) \ge \frac{n}{K+1}$.
	Consequently, we have:
	\begin{equation}
	\begin{split}
	&\sum_{\debug n=0}^{\infty} \Big(\debug\alpha_{\debug n+1}\sum_{r = s(\debug n)}^{\debug n-1} \debug\alpha_{r+1}\Big) \leq 
	\sum_{\debug n=0}^{\infty} \Big(\frac{1}{\debug n} \sum_{r = s(\debug n) }^{\debug n} \frac{1}{r}\Big) \leq 
	\sum_{\debug n=0}^{\infty} \Big(\frac{1}{\debug n}   \frac{Ks^p(\debug n)}{s(\debug n)}\Big)= K\sum_{\debug n=0}^{\infty} \frac{1}{\debug n}  s^{p-1}(\debug n) \leq K\sum_{\debug n=0}^{\infty} \frac{1}{\debug n}   (\frac{n}{K+1})^{p-1}\\
	& \leq  (K+1)^{-p}\sum_{\debug n=0}^{\infty} \debug n^{p-2}
	< \infty.
	\end{split}
	\end{equation}
	\item In the linearly growing delays case, since $d_{\debug n} = O(\debug n)$, it follows that $s(n) + K s(\debug n) \ge \debug n$ for some positive number $K$ and hence again
	$s(n) \ge \frac{n}{K+1}$.
	Consequently, we have:
\begin{flalign}
\sum_{\debug n=0}^{\infty} \Big(\debug\alpha_{\debug n+1}\sum_{r = s(\debug n)}^{\debug n-1} \debug\alpha_{r+1}\Big)
	&\leq  \sum_{\debug n=0}^{\infty} \Big(\frac{1}{\debug n \log \debug n} \sum_{r = s(\debug n) }^{\debug n} \frac{1}{r \log r}\Big)
	\leq \sum_{\debug n=0}^{\infty} \Big(\frac{1}{\debug n \log \debug n} \sum_{r = s(\debug n) }^{s(\debug n) + Ks(\debug n)} \frac{1}{r \log r}\Big)
	\notag\\
	&\leq  \sum_{\debug n=0}^{\infty} \Big(\frac{1}{\debug n \log \debug n}     \int_ {s(n)}^{s(\debug n) + K s(\debug n)} \frac{1}{r\log r} dr \Big)
	= \sum_{\debug n=0}^{\infty} \Big(\frac{1}{\debug n \log \debug n} \log \frac{\log (s(\debug n) + K s(\debug n))}{\log s(\debug n)} \Big)
	\notag\\
	&= \sum_{\debug n=0}^{\infty} \Big(\frac{1}{\debug n \log \debug n} \log \frac{\log (K+1) + \log s(\debug n) }{\log s(\debug n)}            \Big)
	= \sum_{\debug n=0}^{\infty} \Big(\frac{1}{\debug n \log \debug n} \log (1 + \frac{\log (K+1)}{\log s(\debug n)}  )          \Big)
	\notag\\
	& \leq \sum_{\debug n=0}^{\infty} \Big(\frac{1}{\debug n \log \debug n} \frac{\log (K+1)}{\log s(\debug n)} \Big)
	\leq \sum_{\debug n=0}^{\infty} \Big(\frac{1}{\debug n \log \debug n} \frac{\log (K+1)}{\log \debug n - \log (K+1)}
	\notag\\
	& \sim \overline{K}    \sum_{\debug n=0}^{\infty} \frac{1}{\debug n (\log \debug n)^2}     
	< \infty.
\end{flalign}
\item In the polynomially growing delays case, since $d_{\debug n} = O(\debug n^q)$, it follows that $s(n) + K s^q(\debug n) \ge \debug n$ for some positive number $K$. Note that in this case, $s(n) = \Omega (n^{\frac{1}{2q}})$, because otherwise, $s(n) + K s^q(\debug n) =  O (n^{\frac{1}{2q}} + K n^{\frac{1}{2}}) = o(n)$, which is a contradiction.
	Consequently, we have:
	\begin{flalign}
	\sum_{\debug n=0}^{\infty} \Big(\debug\alpha_{\debug n+1}\sum_{r = s(\debug n)}^{\debug n-1} \debug\alpha_{r+1}\Big)
	&\leq  \sum_{\debug n=0}^{\infty} \Big(\frac{1}{\debug n \log \debug n \log\log \debug n} \sum_{r = s(\debug n) }^{\debug n} \frac{1}{r \log r \log\log r}\Big)
	\notag\\
	&\leq \sum_{\debug n=0}^{\infty} \Big(\frac{1}{\debug n \log \debug n \log \log \debug n} \sum_{r = s(\debug n) }^{s(\debug n) + Ks^q(\debug n)} \frac{1}{r \log r \log\log r}\Big)
	\notag\\
	&\leq  \sum_{\debug n=0}^{\infty} \Big(\frac{1}{\debug n \log \debug n \log\log\debug n}  \int_ {s(n)}^{s(\debug n) + K s^q(\debug n)} \frac{1}{r\log r \log\log r} dr \Big)
	\notag\\
	&\leq \sum_{\debug n=0}^{\infty} \Big(\frac{1}{\debug n \log \debug n  \log\log\debug n} \log \frac{\log((K+1)\log s(\debug n) + a\log s(\debug n))}{\log \log s(\debug n)}     \Big)
	\notag\\
	& \leq \sum_{\debug n=0}^{\infty} \Big(\frac{1}{\debug n \log \debug n  \log\log\debug n} \log (1 + \frac{\log(K+1+a) }{\log \log s(\debug n)}  )   \Big)
	\notag\\
	&\leq \sum_{\debug n=0}^{\infty} \Big(\frac{1}{\debug n \log \debug n  \log\log\debug n} \frac{\log(K+1+a) }{\log \log s(\debug n)}    \Big)
	\notag\\
	& \leq \sum_{\debug n=0}^{\infty} \Big(\frac{1}{\debug n \log \debug n  \log\log\debug n} \frac{\log(K+1+a) }{\log \log \frac{n^{\frac{1}{2q}}}{K^\prime} }   \Big)
	\notag\\
	&\sim \overline{K}    \sum_{\debug n=0}^{\infty} \frac{1}{\debug n \log\debug n (\log\log \debug n)^2}     
	< \infty.
	\end{flalign}
\end{enumerate}

Consequently, in each of the above 4 cases, we have $\sum_{\debug n=0}^{\infty} \Big(\debug\alpha_{\debug n+1}\sum_{r = s(\debug n)}^{\debug n-1} \debug\alpha_{r+1}\Big) < \infty$. Therefore, Equation~\eqref{eq:Lip_ex3} yields
\begin{equation}\label{eq:Lip_ex4}
\begin{split}
& -\infty < 
- \sum_{\debug n=0}^{\infty}\debug\alpha_{\debug n+1} \exof{\|\grad \debug f(\debug X_{\debug n})\|_2^2} +  L\bar{B}\sqrt{B} \sum_{\debug n=0}^{\infty} \Big(\debug\alpha_{\debug n+1}\sum_{r = s(\debug n)}^{\debug n-1} \debug\alpha_{r+1}\Big)
+  (L\sqrt{B}+ 4LB) \sum_{\debug n = 0}^{\infty}\debug\alpha_{\debug n+1}^2 \\
&\leq - \sum_{\debug n=0}^{\infty}\debug\alpha_{\debug n+1} \exof{\|\grad \debug f(\debug X_{\debug n})\|_2^2} + \overline{C},\\
\end{split}
\end{equation}
for some finite constant $\overline{C}$. Reversing the above inequality immediately yields the result:
$$\sum_{\debug n=0}^{\infty}\debug\alpha_{\debug n+1} \exof{\|\grad \debug f(\debug X_{\debug n})\|_2^2} < \infty.$$
\hfill$\blacksquare$

\subsubsection{Proof of Lemma~\ref{lem:diff}}

{\em Proof:}
We first recall a useful fact: for any two vectors $\mathbf{a}, \mathbf{b}$ and any finite-dimensional vector norm $\|\cdot\|$, 
\begin{equation}\label{eq:help1}
\Big|(\|a\| + \|b\|)(\|a\|-\|b\|)\Big|\leq \|a+b\|\|a-b\|.
\end{equation}
Using this fact, we can expect bound the term in question as follows:
\begin{align*}
& \bigg|\mathbb{E}[\|\nabla \debug f(\debug X_{\debug n+1})\|_2^2] -
\mathbb{E}[\|\nabla \debug f(\debug X_{\debug n})\|_2^2]\bigg|  = \Bigg| 
\mathbb{E} \Bigg[
\bigg(\|\nabla \debug f(\debug X_{\debug n+1})\|_2 +
\|\nabla \debug f(\debug X_{\debug n})\|_2 \bigg)   
\bigg(\|\nabla \debug f(\debug X_{\debug n+1})\|_2 -
\|\nabla \debug f(\debug X_{\debug n})\|_2 \bigg)   
\Bigg] \Bigg| \\
& \leq  
\mathbb{E} \Bigg[
\bigg|\|\nabla \debug f(\debug X_{\debug n+1})\|_2 +
\|\nabla \debug f(\debug X_{\debug n})\|_2 \bigg| \cdot   
\bigg|\|\nabla \debug f(\debug X_{\debug n+1})\|_2 -
\|\nabla \debug f(\debug X_{\debug n})\|_2 \bigg|   
\Bigg] \\
& \leq 
\mathbb{E} \Bigg[
\bigg\|\nabla \debug f(\debug X_{\debug n+1}) +
\nabla \debug f(\debug X_{\debug n})\bigg\|_2  \cdot   
\bigg\|\nabla \debug f(\debug X_{\debug n+1}) -
\nabla \debug f(\debug X_{\debug n})\bigg\|_2    
\Bigg] 
\leq \mathbb{E} \Bigg[
2\sup_{\debug x \in \mathcal{\debug X}}\|\nabla \debug f(\debug x) \|_2 \cdot   
\bigg\|\nabla \debug f(\debug X_{\debug n+1}) -
\nabla \debug f(\debug X_{\debug n})\bigg\|_2    
\Bigg]\\
& \leq 2\sqrt{B}\mathbb{E} \Bigg[
\bigg\|\nabla \debug f(\debug X_{\debug n+1}) -
\nabla \debug f(\debug X_{\debug n})\bigg\|_2    
\Bigg] 
\leq 2\sqrt{B}\mathbb{E} \Bigg[
L\bigg\|\debug X_{\debug n+1} -
\debug X_{\debug n}\bigg\|_2    
\Bigg] 
= 2L\sqrt{B}\mathbb{E} \Bigg[
\bigg\|\debug X_{\debug n+1} -
\debug X_{\debug n}\bigg\|_2    
\Bigg] \\
&= 2L\sqrt{B}\mathbb{E} \Bigg[
\bigg\|   \debug\alpha_{\debug n+1} \nabla \debug F(\debug X_{s(\debug n)},\debug\omega_{\debug n+1})         \bigg\|_2    
\Bigg] \le
2L\sqrt{B}\debug\alpha_{\debug n+1}\sup_{\debug x \in \mathcal{\debug X}}\mathbb{E} \Bigg[
\bigg\|    \nabla \debug F(\debug x,\debug\omega_{\debug n+1})         \bigg\|_2    
\Bigg]
	\notag\\
	&= 2L\sqrt{B}\debug\alpha_{\debug n+1}\sup_{\debug x \in \mathcal{\debug X}}\mathbb{E} \Bigg[
\sqrt{\bigg\|   \nabla \debug F(\debug x,\debug\omega_{\debug n+1})         \bigg\|_2 ^2   }
\Bigg]  \\
& \leq 2L\sqrt{B}\debug\alpha_{\debug n+1}\sup_{\debug x \in \mathcal{\debug X}} 
\sqrt{\mathbb{E}\Bigg[\bigg\|   \nabla \debug F(\debug x,\debug\omega_{\debug n+1})         \bigg\|_2 ^2  \Bigg] }
	\notag\\
	&\leq 2L\sqrt{B}\debug\alpha_{\debug n+1}\sqrt{B} =2LB\debug\alpha_{\debug n+1},
\end{align*}
where the first inequality is an application of Jensen's inequality,
the second inequality follows from Equation~\eqref{eq:help1}
the fifth inequality follows from Lipschitz continuity and the second-to-last inequality follows from another application of Jensen's inequality and that the squre root function is concave.
\hfill$\blacksquare$

\subsection{Variationally Coherent Objectives}
We define the following constants that will be handy later:
\begin{enumerate}
	\item $C_ 1 = \sup_{\debug x \in \mathcal{\debug X}}\exof{\dnorm{\nabla\debug F(\debug x;\debug\omega)}}$.
	\item $C_ 2 = \sup_{\debug x \in \mathcal{\debug X}}\exof{\dnorm{\nabla\debug F(\debug x;\debug\omega)}^{2}}$.
	\item  $C_3$ is the Lipschitz constant for $\grad\debug f(\debug x) (= \exof{\nabla\debug F(x;\debug\omega)}):$
	\begin{equation}
	\|\grad\debug f(\debug x) - \grad\debug f(\debug x^{\prime})\|_2 \leq C_3 \|\debug x - \debug x^{\prime}\|_2.
	\end{equation}


	\item $C_4 = \sup_{x, x^{\prime} \in \mathcal{\debug X}} \|x - x^{\prime}\|_2.$
\end{enumerate}




\subsubsection{Proof of Lemma~\ref{lem:energy}}

	Per the definition of the energy function, we have:
	\begin{equation}
	\begin{split}
	&\fench_{\debug x^{\ast}}(\debug y) - \|\mirror(\debug y) - \debug x^{\ast}\|_2^2 = \|\debug x^{\ast}\|^2_2 -  \|\mirror(\debug y)\|_2^2 + 
	2\langle\debug y, \mirror(\debug y) - \debug x^{\ast}\rangle - \Big\{\|\mirror(\debug y) \|_2^2 - 2\langle \mirror(\debug y), \debug x^{\ast} \rangle + \|\debug x^{\ast}\|_2^2 \Big\} \\
	&= -2\|\mirror(\debug y) \|_2^2 + 2\langle \debug y - \mirror(\debug y), \mirror(\debug y) - \debug x^{\ast}\rangle + 
	2\langle \mirror(\debug y), \mirror(\debug y) - \debug x^{\ast}\rangle +
	2\langle \mirror(\debug y), \debug x^{\ast} \rangle \\
	& = 2\langle \debug y - \mirror(\debug y), \mirror(\debug y) - \debug x^{\ast}\rangle \ge 0,
	\end{split}
	\end{equation}
	where the last inequality follows from Lemma~\ref{lem:aux1}.
	Consequently, $\fench_{\debug x^{\ast}}(\debug y) \ge \|\mirror(\debug y) - \debug x^{\ast}\|_2^2 \ge 0$.
	Taking the infimum over $\feas^{\ast}$ yields $\fench(\debug y) \ge 0$.
	
	Next, we establish that $\fench(\debug y) = 0$ if and only if $\mirror(\debug y) \in \feas^{\ast}$. The if part is already established in Remark~\ref{rem:prop1}.
	It suffices to show that $\fench(\debug y) = 0$ implies $\mirror(\debug y) \in \feas^{\ast}$.
	To see this, observe that if $\fench_{\debug x^{\ast}}(\debug y) = 0$, we must have $\|\mirror(\debug y) - \debug x^{\ast}\|_2^2 = 0$, therefore implying
	$\mirror(\debug y) = \debug x^{\ast}.$ Since $\fench_{\debug x^{\ast}}(\debug y)$ is a continuous function of $\debug x^{\ast}$ for each fixed $\debug y$, and since $\feas^{\ast}$ is a compact set, $\inf_{\debug x^{\ast} \in \feas^{\ast}} \fench_{\debug x^{\ast}}(\debug y)$ must be achieved by
	some $z^* \in \feas^{\ast}$. Namely $\fench(\debug y) = \fench_{z^*}(\debug y)$.
	Consequently, by the preceding observation, $\fench(\debug y) = 0$ implies $\mirror(\debug y) = z^* \in \feas^{\ast}$.
	
	For the second statement, suppose on the contrary $E(\debug y_{\debug n}) \to 0$ but
	$\mirror(\debug y_{\debug n})$ does not converge to $\feas^{\ast}$. Then there must exist  a subsequence $\debug n_k$ such that $\mirror(\debug y_{\debug n_k})$ is bounded away from $\feas^{\ast}$. 
	Denoting by $\mathcal{N}(\feas^{\ast}, \epsilon)$ the $\epsilon$-open ball around $\feas^{\ast}$ (i.e. $\mathcal{N}(\feas^{\ast}, \epsilon) \triangleq \{x\in\mathbb{R}^{\debug d} \mid \text{dist}(x, \feas^{\ast}) < \epsilon\}$).
	Then for the subsequence $\debug n_k$, there must exist some positive $\epsilon$ such that $\mirror(\debug y_{\debug n})$ remains in $\mathcal{\debug X} \cap \mathcal{N}^c(\feas^{\ast}, \epsilon)$.
	Since $\mathcal{\debug X} \cap \mathcal{N}^c(\feas^{\ast}, \epsilon)$ is an intersection of two closed sets, it is itself closed; it is also bounded since $\mathcal{\debug X}$ is bounded: hence it is compact. Further, since for each fixed $\debug x^{\ast}$, $E_{\debug x^{\ast}}(\debug y)$ is a continuous function 
	of $\mirror(\debug y)$, $E_{\debug x^{\ast}}(\debug y)$ must achieve the minimum value on the compact set $\mathcal{\debug X} \cap \mathcal{N}^c(\feas^{\ast}, \epsilon)$, where the minimum value $a_{\debug x^{\ast}}$ is positive per the first statement:
	$$E_{\debug x^{\ast}}(\debug y) \ge a_{\debug x^{\ast}} > 0, \forall \mirror(\debug y) \in \mathcal{\debug X} \cap \mathcal{N}^c(\feas^{\ast}, \epsilon), \forall \debug x^{\ast} \in \feas^{\ast}$$
	 Finally, since $E_{\debug x^{\ast}}(\debug y)$ is continuous in $\debug x^{\ast}$, it must achieve the minimum value (over $\debug x^{\ast}$) on the compact set $\feas^{\ast}$, where the minimum value $a$ (corresponding to some $E_{\debug x^{\ast}}(\debug y)$) must again be positive:
	 $$E(\debug y) = \inf_{\debug x^{\ast} \in \feas^{\ast}} E_{\debug x^{\ast}}(\debug y) \ge \inf_{\debug x^{\ast} \in \feas^{\ast}} a_{\debug x^{\ast}} = a > 0.$$
	
	Consequently, $E(\debug y_{\debug n_k}) \ge a > 0$ since $\mirror(\debug y_{\debug n_k}) \in \mathcal{\debug X} \cap \mathcal{N}^c(\feas^{\ast}, \epsilon), \forall k$.
	However, on this subsequence, the energy function still converges to $0$ by assumption: $E(\debug y_{\debug n_k}) \to 0$, which immediately yields a contradiction.
	The claim is hence established.
\hfill $\blacksquare$




\subsubsection{Proof of Lemma~\ref{lem:energy_prop}}
	We first prove the first claim. By expanding it, we have:
	\begin{equation}\label{eq:iden}
	\begin{split}
	&\|\mirror(\debug y) - \hat{\debug y}\|_2^2 - \|\mirror(\hat{\debug y}) -\hat{\debug y}\|_2^2 =
	\|\mirror(\debug y) - \debug y + \debug y - \hat{\debug y}\|_2^2 - \|\mirror(\hat{\debug y}) -\hat{\debug y}\|_2^2  \\
	& = \|\debug y - \hat{\debug y}\|_2^2 + \|\mirror(\debug y) - \debug y\|_2^2 + 2\langle \mirror(\debug y) - \debug y , \debug y - \hat{\debug y} \rangle- \|\mirror(\hat{\debug y}) -\hat{\debug y}\|_2^2 \\
	& = \|\debug y - \hat{\debug y}\|_2^2 -\Big\{ \|\mirror(\hat{\debug y}) -\hat{\debug y}\|_2^2  - \|\mirror(\debug y) - \debug y\|_2^2 - 2\langle \debug y - \mirror(\debug y)  ,  \hat{\debug y} -\debug y  \rangle\Big\}. \\
	\end{split}
	\end{equation}
	Now define the function $f(\debug x, \debug y) = \|\debug x- \debug y\|_2^2$. It follows easily that the solution to the problem $\max_{\debug x \in \mathcal{\debug X}} \|\debug x- \debug y\|_2^2$  is $x^*(\debug y) = \mirror(\debug y)$.
	Consequently, by Lemma~\ref{lem:env}, $V(\debug y) = f(x^*(\debug y), \debug y)$ is a differential function in $\debug y$ and its derivative can be computed explicitly as follows:
	\begin{equation}\label{eq:gradient}
	\grad V(\debug y) =  \grad f(\debug x^{\ast} (\debug y), \debug y) = 2 (\debug y - x^*(\debug y)) = 2(\debug y - \mirror(\debug y)).
	\end{equation}
	
	Futher, since for each $\debug x \in \mathcal{\debug X}$, $f(\debug x, \debug y)$ is a convex function in $\debug y$, and taking the maximum preserves convexity, we have $V(\debug y)$ is also a convex function in $\debug y$.
	This means that $$V(\hat{\debug y}) - V(\debug y) - \langle \grad V(\debug y), \hat{\debug y} -\debug y\rangle \ge 0.$$
	By Equation~\eqref{eq:gradient} and that $V(\debug y) = f(x^*(\debug y), \debug y) = \|\mirror(\debug y) -\debug x\|_2^2$, the above equation becomes:
	$$ \|\mirror(\hat{\debug y}) -\hat{\debug y}\|_2^2  - \|\mirror(\debug y) - \debug y\|_2^2 - 2\langle \debug y - \mirror(\debug y)  ,  \hat{\debug y} -\debug y \rangle \ge 0.$$
	Consequently, Equation~\eqref{eq:env_opt} then immediately yields:
	$$\|\mirror(\debug y) - \hat{\debug y}\|_2^2 - \|\mirror(\hat{\debug y}) -\hat{\debug y}\|_2^2 \le
	\|\debug y-\hat{\debug y}\|_2^2.$$
	
	We now prove the second part. Expanding using the definition of the Lyapunov function (and skipping some tedious algebra in between), we have:
	\begin{flalign}
	\label{eq:energy_change}
	\fench_{\debug x^{\ast}}(\debug y + \Delta) - \fench_{\debug x^{\ast}}(\debug y)
	&= \|\debug x^{\ast}\|^2_2 -  \|\mirror(\debug y + \Delta)\|_2^2 + 
	2\langle\debug y + \Delta, \mirror(\debug y + \Delta) - \debug x^{\ast}\rangle
	\notag\\
	&- \Big\{ \|\debug x^{\ast}\|^2_2 -  \|\mirror(\debug y)\|_2^2 + 
	2\langle\debug y, \mirror(\debug y) - \debug x^{\ast}\rangle  \Big\}
	\notag\\
	& = \|\mirror(\debug y)\|_2^2  - \|\mirror(\debug y + \Delta)\|_2^2  - 2\langle\debug y, \mirror(\debug y) - \debug x^{\ast}\rangle  +  2\langle \debug y + \Delta, \mirror(\debug y + \Delta) - \debug x^{\ast} \rangle
	\notag\\
	& = \|\mirror(\debug y)\|_2^2  - \|\mirror(\debug y + \Delta)\|_2^2 +
	2\langle y, \mirror(\debug y + \Delta) - \mirror(\debug y) \rangle
	\notag\\
	&\hfill
	+2\langle  \Delta, \mirror(\debug y) - \debug x^{\ast} + \mirror(\debug y + \Delta) -\mirror(\debug y)  \rangle
	\notag\\
	&= 2\langle  \Delta, \mirror(\debug y) - \debug x^{\ast}  \rangle +
	2\langle \debug y + \Delta, \mirror(\debug y + \Delta) - \mirror(\debug y) \rangle
	+ \|\mirror(\debug y)\|_2^2  - \|\mirror(\debug y + \Delta)\|_2^2
	\notag\\
	& =  2\langle  \Delta,  \mirror(\debug y)- \debug x^{\ast}  \rangle +
	\|\mirror(\debug y) - (\debug y + \Delta)\|_2^2 - \|\mirror(\debug y + \Delta) - (\debug y + \Delta)\|_2^2
	\notag\\
	& \leq 2\langle  \Delta,  \mirror(\debug y) - \debug x^{\ast}  \rangle +
	\|\Delta\|_2^2
	\end{flalign}
where the last equality follows from completing the squares and the last inequality follows from the first part of the lemma.
\hfill $\blacksquare$

\subsubsection{Proof of Proposition~\ref{prop:recur_d}}
	We provide the details for all the steps.
	\begin{enumerate}
		\item 
		Defining $b_{\debug n} = \v(x_{s(n)}) - \v(x_\debug n)$,
		we can rewrite the gradient update in \ac{DAGD} as:
		\begin{flalign}
		y_{\debug n+1}
		&= y_{\debug n} - \alpha_{\debug n+1} \v(x_{s(n)})
		\notag\\
		&= y_{\debug n} - \alpha_{\debug n+1} \v(x_\debug n)  - \alpha_{\debug n+1} \{\v(x_{s(n)}) - \v(x_\debug n)\}
		\notag\\
		&= y_{\debug n} - \alpha_{\debug n+1} (\v(x_\debug n) + b_{\debug n}).
		\end{flalign}
		Recall here once again that  $s(\debug n)$ denotes the previous iteration count whose gradient becomes available only at the current iteration $\debug n$.
	To establish the claim, we start by expanding $b_{\debug n}$ as follows:
		\begin{flalign}
		\|b_{\debug n}\|_2 
		&= \|\v(\debug x_{s(\debug n)}) - \v(\debug x_\debug n)\|_2 \leq C_3 \|\debug x_{s(\debug n)} - \debug x_\debug n\|_2
		\notag\\
		& = C_3 \|\mirror(\debug y_{s(\debug n)}) - \mirror(\debug y_\debug n)\|_2  \leq  C_3 \|\debug y_{s(\debug n)} - \debug y_\debug n\|_2
		\notag\\
		& \leq  C_3 \Big\{ \|\debug y_{s(\debug n)} - \debug y_{{s(\debug n)}+1}\|_2 + \|\debug y_{{s(\debug n)}+1} - \debug y_{{s(\debug n)}+2}\|_2 + \dots +  \|\debug y_{\debug n-1} - \debug y_{\debug n}\|_2  \Big\}
		\notag\\
		&= C_3 \sum_{r = s(\debug n)}^{\debug n-1} \|\debug\alpha_{r+1} \v(\debug x_{s(r)}) \|_2 \leq  C_3\sup_{\debug x \in \mathcal{\debug X}} \|\v(\debug x)\|_2 \sum_{r = s(n)}^{n-1} \debug\alpha_{r+1} = C_3V_{\max} \sum_{r = s(n)}^{n-1} \debug\alpha_{r+1}.
		\end{flalign}
		We now consider two cases, depending on whether the delays are bounded or not.
		\begin{enumerate}
			\item
			If $\{\debug\alpha_{\debug n}\}_{\debug n = 1}^{\infty}$ and $d_{\debug n} \leq D, \forall \debug n$ satisfy \cref{asm:delays}, then $d_{s(\debug n)} = \debug n - s(\debug n) \leq D$. Consequently, 
			\begin{equation}
			0
			< \sum_{r = s(n)}^{n-1} \debug\alpha_{r+1} = \sum_{r = s(n)+1}^{n} \debug\alpha_{r}\leq \sum_{r = n-D}^{n} \debug\alpha_{r}
			\leq D \max_{r \in \{n-D, \dots, n\}} \debug\alpha_{r} \to 0
			\quad\text{ as $n\to\infty$}, 
			\end{equation}
			where the limit approaching $0$ follows from $\lim_{\debug n \to \infty} \debug\alpha_{\debug n} = 0$, which itself is a consequence of \cref{asm:delays}. This implies $\lim_{\debug n \to \infty} C_3V_{\max} \sum_{r = s(n)}^{n-1} \debug\alpha_{r+1} = 0$ and consequently, $\lim_{\debug n \to \infty} \|b_{\debug n}\|_2 = 0$.
			
			\item 
			We consider each of the three conditions in turn.
			
			When $\debug\alpha_{\debug n-1} = \frac{1}{\debug n}$  and $d_{\debug n}= o(\debug n)$,
			we have $\debug n - s(\debug n) \leq K o(s(\debug n))$ for some universal constant $K > 0$, which means $\debug n \leq s(\debug n) + K o(s(\debug n))$. 
			Consequently, we have:
			\begin{flalign}
			0
			&< \sum_{r = s(n)}^{n-1} \debug\alpha_{r+1} = \sum_{r = s(n) }^{s(\debug n) + K o(s(\debug n))} \debug\alpha_{r}
			\leq \int_ {s(n)}^{s(\debug n) + K o(s(\debug n))} \frac{1}{r} dr
			\notag\\
			&= \log (s(\debug n) + K o(s(\debug n))) - \log s(\debug n)
			= \log \frac{Ko(s(\debug n)) + s(\debug n)}{s(\debug n)} \to \log 1 = 0
			\quad \text{as $\debug n\to\infty$},
			\end{flalign}
			where the last limit follows from $s(\debug n) \to \infty$ as $\debug n \to \infty$ because
			$\debug n \leq s(\debug n) + K o(s(\debug n))$. 
			
			When $\debug\alpha_{\debug n-1} = \frac{1}{\debug n \log \debug n }$  and $d_{\debug n}= O(\debug n)$,
			it is easy to verify (by integration) that this particular choice of sequence satisfies $\sum_{\debug n=\debug 1}^{\infty} \debug\alpha_{\debug n}^{2}
			<\infty,
			\sum_{\debug n=\debug 1}^{\infty} \debug\alpha_{\debug n}
			= \infty.$
			Since $d_{\debug n}= O(\debug n)$, we have $\debug n - s(\debug n) \leq K s(\debug n)$ for some universal constant $K > 0$, which means $\debug n \leq s(\debug n) + K s(\debug n)$. 
			Consequently, we have:
			\begin{flalign}
			0
			&< \sum_{r = s(n)}^{n-1} \debug\alpha_{r+1} = \sum_{r = s(n) }^{s(\debug n) + K s(\debug n)} \debug\alpha_{r}
			\leq \int_ {s(n)}^{s(\debug n) + K s(\debug n)} \frac{1}{r\log r} dr
			\notag\\
			&= \log \frac{\log (s(\debug n) + K s(\debug n))}{\log s(\debug n)}
			= \log \frac{\log (K+1) + \log s(\debug n) }{\log s(\debug n)} \to 0
			\quad \text{as $\debug n\to\infty$},
			\end{flalign}
			where the last limit follows from $s(\debug n) \to \infty$ as $\debug n \to \infty$ because
			$\debug n \leq s(\debug n) + K s(\debug n)$. 
			
			When $\debug\alpha_{\debug n-1} = \frac{1}{\debug n \log \debug n \log\log\debug n}$  and $d_{\debug n}= O(\debug n^a), a>1$,
			it is again easy to verify (by integration) that this particular choice of sequence satisfies \cref{asm:delays}.
			Since $d_{\debug n}= O(\debug n^a)$, we have $\debug n - s(\debug n) \leq K s(\debug n)^a$ for some universal constant $K > 0$, which means $\debug n \leq s(\debug n) + K s(\debug n)^a$. 
			Consequently, we have:
			\begin{flalign}
			0 < \sum_{r = s(n)}^{n-1} \debug\alpha_{r+1} = \sum_{r = s(n) }^{s(\debug n) + K s(\debug n)^a} \debug\alpha_{r}
			&\leq \int_ {s(n)}^{s(\debug n) + K s(\debug n)^a} \frac{1}{r\log r \log\log r}dr
			\notag\\
			&= \log \frac{\log\log (s(\debug n) + K s(\debug n)^a)}{\log\log s(\debug n)}
			\notag\\
			&\leq  \log \frac{\log\log (s(\debug n)^a + K s(\debug n)^a)}{\log\log s(\debug n)}
			\notag\\
			&< \log \frac{\log((K+1)\log s(\debug n) + a\log s(\debug n))}{\log \log s(\debug n)}
			\notag\\
			&= \log \frac{\log(K+1+a) + \log\log s(\debug n)}{\log \log s(\debug n)} \to 0
			\quad
			\text{as $\debug n \to \infty$},
			\end{flalign}
			where the last limit follows from the fact that $s(\debug n) \to \infty$ as $\debug n \to \infty$ (again, because
			$\debug n \leq s(\debug n) + K s(\debug n)^a$). 
			
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
			
		\end{enumerate}
		\item 
		With the defintion of $b_{\debug n}$, \ac{DAGD} can be written as:
		\begin{equation}
		\begin{aligned}
		\debug x_{\debug n}
		&= \mirror(\debug y_{\debug n}),
		\\
		\debug y_{\debug n+1}
		&= y_{\debug n} - \alpha_{\debug n+1} (\v(x_\debug n) + b_{\debug n}).
		\end{aligned}
		\end{equation}
To prove the claim, we start by fixing an arbitrary $\debug x^{\ast} \in \feas^{\ast}$ and applying Lemma~\ref{lem:energy_prop} to bound the energy change in a single gradient update as follows:
		\begin{equation}\label{eq:energy_change}
		\begin{split}
		\fench_{\debug x^*}(\debug y_{\debug n+1}) - \fench_{\debug x^*}(\debug y_{\debug n})
		&\leq 2\langle \debug y_{\debug n+1} - \debug y_{\debug n}, \mirror(\debug y_{\debug n}) -\debug x^{\ast} \rangle
		 + 
		 \|\debug y_{\debug n}-\debug y_{\debug n+1}\|_2^2\\
		 &= -2\langle  \alpha_{\debug n+1} (\v(x_\debug n) + b_{\debug n}), \debug x_{\debug n} - \debug x^*  \rangle +
		\|\debug y_{\debug n}-\debug y_{\debug n+1}\|_2^2. \\
		\end{split}
		\end{equation}

		
		Now telescoping the above inequality yields:
		\begin{flalign}
		\label{eq:telescoping}
		\fench_{\debug x^*}(\debug y_{\debug n+1}) - \fench_{\debug x^*}(\debug y_0)
		&= \sum_{r=0}^{\debug n} \{\fench_{\debug x^*}( \debug y_{r+1}) - \fench_{\debug x^*}(\debug y_{r})\}
		\notag\\
		& \leq \sum_{r=0}^{\debug n} \{-2 \alpha_{r+1} \langle  \v(x_r) + b_r, \debug x_{r} - \debug x^*  \rangle +
		\alpha_{r+1}^2\| \v(x_r) + b_{r}\|_2^2 \}
		\notag\\
		&\leq -2 \sum_{r=0}^{\debug n} \alpha_{r+1}\{ \langle  \v(x_r), \debug x_{r} - \debug x^*  \rangle  - \|b_r\|_2 \|\debug x_{r} - \debug x^*  \|_2 \}+
		2\sum_{r=0}^{\debug n}\alpha_{r+1}^2\{\| \v(x_r)\|_2^2  + \|b_{r}\|_2^2 \}
		\notag\\
		& \leq -2 \sum_{r=0}^{\debug n} \alpha_{r+1}\{ \langle  \v(x_r), \debug x_{r} - \debug x^*  \rangle  - C_4\|b_r\|_2  \}+
		2\sum_{r=0}^{\debug n}\alpha_{r+1}^2\{C_2  + B \},
		\end{flalign}
		where the last inequality follows from the fact that $b_{\debug n}$'s must be bounded
		(since $\lim_{\debug n \to \infty} \|b_{\debug n}\|_2 = 0$) and hence let $B \triangleq  \sup_{\debug n} \|b_{\debug n}\|_2$.
		By \cref{asm:delays}, we have $2\sum_{r=0}^{\debug n}\alpha_{r+1}^2\{C_2  + B \} = \overline{B} < \infty$. Now fix any positive number $\epsilon$. Assume for contradiction purposes
		$\debug x_{\debug n}$ only enters $\mathcal{N}(\debug x^*, \epsilon)$ a finite number of times and let $t_1$ be the last time 
		this occurs. This means that for all $ r > t_1$, $\debug x_r $ is outside the open set $\mathcal{N}(\feas^{\ast}, \epsilon)$.
		Therefore, since a continuous function always achieves its minimum on a compact set, we have:
		$\langle  \v(x_r), \debug x_{r} - \debug x^*  \rangle \ge \min_{\debug x \in \mathcal{\debug X} -  \mathcal{N}(\feas^{\ast}, \epsilon)}\langle  \v(x), \debug x - \debug x^*  \rangle \triangleq a > 0, \forall r > t_1$ (note that $a$ depends on $\epsilon$).
		Further, since $b_r \to 0$ as $r \to \infty$, pick $t_2$ such that $\|b_r\|_2 < \frac{a}{2C_4}, \forall r \ge t_2$. Denoting $t = \max (t_1, t_2)$, we 
		continue the chain of inequalities in Equation~\eqref{eq:telescoping} below:
		\begin{equation}\label{eq:telescoping2}
		\begin{split}
		- \fench_{\debug x^*}(\debug y_0)
		\notag\\
		&\leq \fench_{\debug x^*}(\debug y_{\debug n+1}) - \fench_{\debug x^*}(\debug y_0)  \le
		-2 \sum_{r=0}^{t} \alpha_{r+1}\{ \langle  \v(x_r), \debug x_{r} - \debug x^*  \rangle  - C_4\|b_r\|_2  \} \\
		&  - 2 \sum_{r=t+1}^{\debug n} \alpha_{r+1}\{ \langle  \v(x_r), \debug x_{r} - \debug x^*  \rangle  - C_4\|b_r\|_2  \} +
		2\sum_{r=0}^{\debug n}\alpha_{r+1}^2\{C_2  + B \} \\
		& \leq  -2 \sum_{r=0}^{t} \alpha_{r+1}\{ \langle  \v(x_r), \debug x_{r} - \debug x^*  \rangle  - C_4\|b_r\|_2  \} 
		-2 \sum_{r=t+1}^{\debug n} \alpha_{r+1}\{ a - C_4\|b_r\|_2  \} + \overline{B} \\
		& \leq 2 C_4\sum_{r=0}^{t} \alpha_{r+1} |b_r\|_2   + \overline{B} - 2 \sum_{r=t+1}^{\debug n} \alpha_{r+1}\{ a - \frac{a}{2}  \}  \\
		& = \overline{\overline{B}} -  a\sum_{r=t+1}^{\debug n} \alpha_{r+1}  
		\to -\infty, \text{ as } \debug n \to \infty
		\end{split}
		\end{equation}
		where the first inequality follows from the energy function always being positive (\cref{lem:energy}), the second-to-last inequality follows from variational coherence and the limit on the last line follows from
		\cref{asm:delays} and that
		$\overline{\overline{B}} \triangleq  2 C_4\sum_{r=0}^{t} \alpha_{r+1} |b_r\|_2   + \overline{B}$ is just some finite constant.  
		This yields an immediate contradiction and the claim is therefore established.
	\end{enumerate}
\hfill$\blacksquare$


\subsubsection{Proof of Theorem~\ref{thm:convergence}}


	Fix a given $\delta > 0$.
	Since $\debug\alpha_{\debug n} \to 0, b_\debug n \to 0$ as $\debug n \to \infty$, for any $a > 0$, we can pick an $N$ large enough (depending on $\delta$ and $a$) such that $\forall n \ge N$, the following three statements all hold:
	\begin{equation}\label{eq:threshold}
	\begin{split}
	& 2BC_4\debug\alpha_{\debug n+1}  +  2\debug\alpha_{\debug n+1}^2(C_2 + B^2) \leq \frac{\delta}{2},\\
	& C_4 \| b_\debug n\|_2 \leq \frac{a}{2},\\
	&  \debug\alpha_{\debug n+1}(C_2 + B^2) < \frac{a}{2} .
	\end{split}
	\end{equation}
	
	
	We show that under either of the following two (exhaustive) possibilities, if $\fench(\debug y_{\debug n})$ is less than $\delta$,
	$\fench(\debug y_{\debug n+1})$ is less than $\delta$ as well, where $n \ge N$.
	\begin{enumerate}
		\item 
		Case 1: $\fench(\debug y_{\debug n}) < \frac{\delta}{2}$.
		\item
		Case 2: $\frac{\delta}{2} \leq \fench(\debug y_{\debug n}) < \delta$.
	\end{enumerate}
	
	Under Case 1, it follows from Equation~\eqref{eq:energy_change}:
	\begin{equation}
	\begin{split}
	&\fench_{\debug x^*}(\debug y_{\debug n+1}) -
	\fench_{\debug x^*}(\debug y_{\debug n}) \leq -2\debug\alpha_{\debug n+1} \langle \v(\debug x_\debug n) + b_\debug n, x_\debug n - x^*\rangle +  \debug\alpha_{\debug n+1}^2\|\v(\debug x_\debug n) + b_\debug n \|_2^2 \\
	& \leq -2\debug\alpha_{\debug n+1} \langle b_\debug n, x_\debug n - x^*\rangle +  \debug\alpha_{\debug n+1}^2\|\v(\debug x_\debug n) + b_\debug n \|_2^2 \\
	& \leq 2\debug\alpha_{\debug n+1} \| b_\debug n\|_2 \|x_\debug n - x^*\|_2+  2\debug\alpha_{\debug n+1}^2(C_2 + B^2) \\
	&\leq 
	2BC_4\debug\alpha_{\debug n+1}  +  2\debug\alpha_{\debug n+1}^2(C_2 + B^2) 


	 \leq \frac{\delta}{2},
	\end{split}
	\end{equation}
	where the second inequality follows from variational coherence.
	Taking the infimum of $\debug x^{\ast}$ over $\feas^{\ast}$ then yields:
	 $\fench(\debug y_{\debug n+1}) -
	 \fench(\debug y_{\debug n}) \leq \frac{\delta}{2}.$
	This then implies that $\fench{(\debug x^*, \debug y_{\debug n+1})} \leq 
	\fench(\debug x^*, \debug y_{\debug n}) +  \frac{\delta}{2} < \delta$.

	
	Under Case 2, \cref{eq:energy_change} readily yields:
	\begin{equation}
	\begin{split}
	& \fench_{\debug x^*}(\debug y_{\debug n+1}) -
	\fench_{\debug x^*}(\debug y_{\debug n}) \leq -2\debug\alpha_{\debug n+1} \langle \v(\debug x_\debug n) + b_\debug n, x_\debug n - x^*\rangle +  \debug\alpha_{\debug n+1}^2\|\v(\debug x_\debug n) + b_\debug n \|_2^2 \\
	& = -2\debug\alpha_{\debug n+1} \langle \v(\debug x_\debug n), x_\debug n - x^*\rangle -2\debug\alpha_{\debug n+1} \langle b_\debug n, x_\debug n - x^*\rangle +  \debug\alpha_{\debug n+1}^2\|\v(\debug x_\debug n) + b_\debug n \|_2^2 \\
	& \leq -2\debug\alpha_{\debug n+1} a + 2\debug\alpha_{\debug n+1} \| b_\debug n\|_2 \|x_\debug n - x^*\|_2+  2\debug\alpha_{\debug n+1}^2(C_2 + B^2) \\
	& \leq -2\debug\alpha_{\debug n+1} \Big\{a - C_4 \| b_\debug n\|_2  - \debug\alpha_{\debug n+1}(C_2 + B^2)\Big\} \\
	& \leq -2\debug\alpha_{\debug n+1} \Big\{a - \frac{a}{2}  - \debug\alpha_{\debug n+1}(C_2 + B^2)\Big\} \\
	& = -2\debug\alpha_{\debug n+1} \Big\{ \frac{a}{2}  - \debug\alpha_{\debug n+1}(C_2 + B^2)\Big\} \\
	& < 0,
	\end{split}
	\end{equation}
	where the second inequality follows from $\langle \v(\debug x_\debug n), x_\debug n - x^*\rangle \ge a$ under Case 2\footnote{Here $a$ is a constant that only depends on $\delta$: if $\fench{(\debug x^*, y)} \ge \frac{\delta}{2}$, then by part 2 of \cref{lem:energy}, $\mirror(y)$ must be outside an $\epsilon$-neighborhood of $\debug x^*$, for some $\epsilon > 0$. On this neighborhood, the strictly positive continuous function $\langle \v(\debug x), x - x^*\rangle$ must achieve a minimum value $a > 0$.}. Taking the infimum of $\debug x^{\ast}$ over $\feas^{\ast}$ then yields:
	$\fench(\debug y_{\debug n+1}) \le
	\fench(\debug y_{\debug n}) .$
	This then implies that $\fench{( \debug y_{\debug n+1})} \leq 
	\fench(\debug y_{\debug n}) <  \frac{\delta}{2}$.
	
	Consequently, putting the above two cases together therefore yields $\fench{(\debug y_{\debug n+1})} < \delta$.
\hfill{$\blacksquare$}

\subsubsection{Proof of Proposition~\ref{prop:recur_s}}

Using the definitions introduced in the main text, we rewrite the gradient update in \ac{DASGD} as:
		\begin{equation}\label{eq:DSLGD_abbre}
		\debug Y_{\debug n+1} = \debug Y_{\debug n} - \debug\alpha_{\debug n+1}\{\v(\debug X_\debug n) + B_{\debug n} + U_{\debug n+1}\}.
		\end{equation} 
		
		\begin{enumerate}
			\item
		To see that 
		$\sum_{r=0}^{\debug n} U_{\debug n+1}$ is a martingale adapted to $\debug Y_0, \debug Y_1 \dots, \debug Y_{\debug n+1}$,
		first note that, by defintion, $B_{\debug n}$ is adapted to $\debug Y_0, \debug Y_1 \dots, \debug Y_{\debug n}$ (since $\debug X_{\debug n}$ is a deterministic function of $\debug Y_{\debug n}$) and $\debug Y_{\debug n+1}, \debug Y_{\debug n}, B_{\debug n}$ together determine $U_{\debug n+1}$.
		We then check that their first moments are bounded: 
		\begin{equation}
		\begin{split}
		&\exof{\|\sum_{r=0}^ \debug n \|U_{r+1} \|_2} \leq   \sum_{r=0}^\debug n \exof{ \|U_{r+1} \|_2} = \sum_{r=0}^\debug n \exof{\| \nabla \debug F(\debug X_{s(r)},\debug\omega_{r+1}) -  \v(\debug X_{s(r)})\|_2}\\
		& \leq \sum_{r=0}^\debug n \Big\{\exof{\| \nabla \debug F(\debug X_{s(r)},\debug\omega_{r+1})\|_2} +  \exof{\|\v(\debug X_{s(r)})\|_2}\Big\} \\
		&\leq 
		\sum_{r=0}^\debug n \Big\{\sup_{\debug x \in \mathcal{\debug X}} \exof{\| \nabla \debug F(\debug x,\debug\omega)\|_2} +  \sup_{\debug x \in \mathcal{\debug X}}\|\v(\debug x)\|_2\Big\}  \\
		&= \sum_{r=0}^\debug n \Big\{\sup_{\debug x \in \mathcal{\debug X}} \exof{\| \nabla \debug F(\debug x,\debug\omega)\|_2} +  \sup_{\debug x \in \mathcal{\debug X}}\|\exof{\nabla \debug F(\debug x, \debug\omega)}\|_2\Big\}\\
		& \leq \sum_{r=0}^\debug n 2 \sup_{\debug x \in \mathcal{\debug X}} \exof{\| \nabla \debug F(\debug x,\debug\omega)\|_2} = \sum_{r=0}^\debug n 2 C_1
		= 2(n+1) C_1 < \infty,
		\end{split}
		\end{equation}
		where the last inequality follows from Jensen's inequality (since $\|\cdot\|_2$ is a convex function).
		Finally, the martingale property holds because:
		$\exof{\sum_{r=0}^ \debug n U_{r+1} \mid \debug Y_1, \dots, \debug Y_{\debug n+1}} = 
		\exof{\nabla \debug F(\debug X_{s(\debug n)},\debug\omega_{\debug n+1}) -  \v(\debug X_{s(\debug n)}) \mid \debug Y_1, \dots, \debug Y_{\debug n}}+ \sum_{r=0}^ {\debug n-1} U_{r+1}   = \sum_{r=0}^ {\debug n-1} U_{r+1}.$, since $\debug\omega_{\debug n+1}$ is \textit{iid}.
		Therefore, $\sum_{r=0}^{\debug n} U_{r+1}$ is a martingale, and $U_{\debug n+1}$ is a martingale difference sequence adapted to $\debug Y_0, \debug Y_1 \dots, \debug Y_{\debug n+1}$.
		
		
		Next, we show that $\lim_{\debug n \to \infty} \|B_{\debug n}\|_2 = 0, { a.s.}.$ By definition, we can expand $B_{\debug n}$ as follows:
		\begin{flalign}
		\label{eq:expansion}
		\|B_{\debug n}\|_2
		&= \|\v(\debug X_{s(\debug n)}) - \v(\debug X_\debug n)\|_2 \leq C_3 \|\debug X_{s(\debug n)} - \debug X_\debug n\|_2
		= C_3 \|\mirror(\debug Y_{s(\debug n)}) - \mirror(\debug Y_\debug n)\|_2
		\notag\\
		&\leq  C_3 \|\debug Y_{s(\debug n)} - \debug Y_\debug n\|_2
		= C_3 \Big\| \debug Y_{s(\debug n)} - \debug Y_{{s(\debug n)}+1} + \debug Y_{{s(\debug n)}+1} - \debug Y_{{s(\debug n)}+2} + \dots +  \debug Y_{\debug n-1} - \debug Y_{\debug n}  \Big\|_2
		\notag\\
		&= C_3\Big\| \sum_{r = s(\debug n)}^{\debug n-1} \Big\{\debug Y_{r} - \debug Y_{r+1} \Big\}\Big\|_2
		= C_3\Big\|\sum_{r = s(\debug n)}^{\debug n-1} \debug\alpha_{r+1} \nabla \debug F(\debug X_{s(r)},\debug\omega_{r+1}) \Big\|_2
		\notag\\
		&= C_3\Big\|\sum_{r = s(\debug n)}^{\debug n-1} \debug\alpha_{r+1} \Big\{\v(\debug X_{s(r)}) + \nabla \debug F(\debug X_{s(r)},\debug\omega_{r+1}) -  \v(\debug X_{s(r)})\Big\}\Big\|_2
		\notag\\
	
		&= C_3\Big\|\sum_{r = s(\debug n)}^{\debug n-1} \debug\alpha_{r+1} \v(\debug X_{s(r)}) +\sum_{r = s(\debug n)}^{\debug n-1} \debug\alpha_{r+1} U_{r+1}\Big\|_2
		\notag\\
		&\leq  C_3 \sum_{r = s(\debug n)}^{\debug n-1} \debug\alpha_{r+1} \|\v(\debug X_{s(r)})\|_2 +  C_3\Big\| \sum_{r = s(\debug n)}^{\debug n-1}\debug\alpha_{r+1}U_{r+1}\Big\|_2
		\notag\\
	
		&\leq C_3C_1 \sum_{r = s(\debug n)}^{\debug n-1} \debug\alpha_{r+1} + C_3 \|\sum_{r = s(\debug n)}^{\debug n-1} \debug\alpha_{r+1} U_{r+1}\|_2
		\notag\\
		& = C_3C_1 \sum_{r = s(\debug n)}^{\debug n-1} \debug\alpha_{r+1} + C_3 \|\sum_{r = 0}^{\debug n-1} \debug\alpha_{r+1} U_{r+1} -\sum_{r = 0}^{s(\debug n)-1} \debug\alpha_{r+1} U_{r+1}\|_2,
	
	
	
		\end{flalign}
		where the first inequality follows from $\grad\debug f$ being Liptichz-continuous (Assumption~\ref{assump:0.6}) and the second inequality follows from $\mirror$ is a non-expansive map.
	
		
		By the same analysis as in the deterministic case, 
		the first part of the last line of Equation~\eqref{eq:expansion} converges to $0$ (under each one of the conditions on step-size and delays in \cref{asm:delays}):
		\begin{equation}\label{eq:limit1}
		\lim_{\debug n \to \infty} C_3C_1 \sum_{r = s(n)}^{n-1} \debug\alpha_{r+1} = 0.
		\end{equation}
		
		
		
		We then analyze the limit of $\|\sum_{r = 0}^{\debug n-1} \debug\alpha_{r+1} U_{r+1} -\sum_{r = 0}^{s(\debug n)-1} \debug\alpha_{r+1} U_{r+1}\|_2$. 
		Define:
		$$ M_{\debug n} = \sum_{r = 0}^{\debug n-1} \debug\alpha_{r+1} U_{r+1}.$$
		Since $U_{r+1}$'s are martingale differences, $M_{\debug n}$ is a martingale.
		Further, in each of the three conditions, $\sum_{\debug n=\debug 1}^{\infty} \debug\alpha_{\debug n}^{2}
		<\infty$. This implies that $M_{\debug n}$ is an $L_2$-bounded martingale because:
		\begin{equation}\label{eq:reason}
		\begin{split}
		& \sup_{\debug n} \exof{\|M_{\debug n}\|_2^2} = \sup_{\debug n} \exof{\Big\|\sum_{r = 0}^{\debug n-1} \debug\alpha_{r+1} U_{r+1}\Big\|_2^2}
		= \sup_{\debug n} \exof{\langle \sum_{r = 0}^{\debug n-1} \debug\alpha_{r+1} U_{r+1}, \sum_{r = 0}^{\debug n-1} \debug\alpha_{r+1} U_{r+1}\rangle}\\
		&= \sup_{\debug n} \exof{\sum_{i,j} \langle  \debug\alpha_{i+1} U_{i+1},  \debug\alpha_{j+1} U_{j+1}\rangle} = \sup_{\debug n} \sum_{r =0}^{\debug n-1} \exof{\langle  \debug\alpha_{r+1} U_{r+1},  \debug\alpha_{r+1} U_{r+1}\rangle} \\
		& =  \sup_{\debug n} \sum_{r =0}^{\debug n-1} \debug\alpha_{r+1}^2\exof{\| U_{r+1}\|_2^2}  \leq \sup_{\debug n} 4C_2\sum_{r =0}^{\debug n-1}  \debug\alpha_{r+1}^2 \leq 4C_2 \sum_{r =0}^{\infty}  \debug\alpha_{r+1}^2 < \infty,
		\end{split}
		\end{equation}
		where the last inequality in the second line follows from the martingale property as follows:
		\begin{equation}
		\begin{split}
		&\exof{\langle  \debug\alpha_{i+1} U_{i+1},  \debug\alpha_{j+1} U_{j+1}\rangle} =
		\debug\alpha_{i+1} \debug\alpha_{j+1} \exof{\langle U_{i+1},  U_{j+1}\rangle} \\
		&= \debug\alpha_{i+1} \debug\alpha_{j+1} \exof{\exof{\langle U_{i+1},  U_{j+1}\rangle \mid Y_0, Y_1, \dots, Y_{i+1}}} \\
		& = \debug\alpha_{i+1} \debug\alpha_{j+1} \exof{\langle U_{i+1}, \exof{U_{j+1} \mid Y_0, Y_1, \dots, Y_{i+1}} \rangle} =
		\debug\alpha_{i+1} \debug\alpha_{j+1} \exof{\langle U_{i+1}, 0\rangle} = 0,
		\end{split}
		\end{equation}
		where the  second equality follows from the tower property (and without loss of generality, we have assumed $i < j$, the third equality follows from $U_{i+1}$ is adapted to $  Y_0, Y_1, \dots, Y_{i+1}$ and the second-to-last equality follows
		from $U_{\debug n +1}$ is a martingale difference.
		Consequently, all the cross terms in the second line of Equation~\eqref{eq:reason} are $0$.
		Therefore, by \cref{thm:mg_convergence2}, by taking $p=2$
		$\lim_{\debug n \to \infty} M_{\debug n} = M_{\infty}, \text{ a.s.} $, where $M_{\infty}$ has finite second-moment.
		Further, since in all three cases $s(\debug n) \to \infty$ as $\debug n \to \infty$ (because there is at most a polynomial lag between $s(\debug n)$ and $\debug n$), we have $\lim_{\debug n \to \infty} M_{s(\debug n)} = M_{\infty}, \text{ a.s.} $. Therefore 
		$$\lim_{\debug n \to \infty} \Big\{\sum_{r = 0}^{\debug n-1} \debug\alpha_{r+1} U_{r+1} -\sum_{r = 0}^{s(\debug n)-1} \debug\alpha_{r+1} U_{r+1}\Big\}= \lim_{\debug n \to \infty} \Big\{M_{\debug n} - M_{s(\debug n)}\Big\}= 0, \text{ a.s.},$$
		thereby implying:
		\begin{equation}\label{eq:limit2}
		\lim_{\debug n \to \infty} C_3 \Big\|\sum_{r = 0}^{\debug n-1} \debug\alpha_{r+1} U_{r+1} -\sum_{r = 0}^{s(\debug n)-1} \debug\alpha_{r+1} U_{r+1}\Big\|_2 = 0.
		\end{equation}
		Combining Equation~\eqref{eq:limit1} and Equation~\eqref{eq:limit2} yields
		$\lim_{\debug n \to \infty} \|B_{\debug n}\|_2 = 0, { a.s.}.$
		

		\item The full \ac{DASGD} update is then:
		\begin{eqnarray}
		& \debug X_{\debug n} = \mirror(\debug Y_{\debug n}) \\
		& \debug Y_{\debug n+1} = \debug Y_{\debug n} - \debug\alpha_{\debug n+1}\{\v(\debug X_\debug n) + B_{\debug n} + U_{\debug n+1}\}. 
		\end{eqnarray}
		We now bound the one-step change of the energy function 
		$\fench(\feas^{\ast}, \debug Y_{\debug n+1}) - \fench(\feas^{\ast}, \debug Y_{\debug n})$ (which is now a random quantity) and then telescope the differences.
		
		Pick an arbitrary $\debug x^{\ast} \in \feas^{\ast}$ and apply Lemma~\ref{lem:energy_prop}, we have:
		\begin{equation}\label{eq:energy_changes}
		\begin{split}
		&\fench_{\debug x^*}(\debug Y_{\debug n+1}) - \fench_{\debug x^*}(\debug Y_{\debug n})  \leq -2\alpha_{\debug n+1} \langle  \v(\debug X_\debug n) + B_{\debug n} + U_{\debug n+1}, \debug X_{\debug n} - \debug x^*  \rangle +
		\|\debug Y_{\debug n}-\debug Y_{\debug n+1}\|_2^2 \\
		& = -2 \alpha_{\debug n+1} \langle  \v(\debug X_\debug n) + B_{\debug n} + U_{\debug n+1}, \debug X_{\debug n} - \debug x^*  \rangle +
		\alpha_{\debug n+1}^2 \|  \v(\debug X_\debug n) + B_{\debug n} + U_{\debug n+1}  \|_2^2 \\
		& \leq -2 \alpha_{\debug n+1} \langle  \v(\debug X_\debug n) + B_{\debug n} + U_{\debug n+1}, \debug X_{\debug n} - \debug x^*  \rangle +
		3\alpha_{\debug n+1}^2 \Big\{\|  \v(\debug X_\debug n)\|_2^2 + \|B_{\debug n} \|_2^2+ \|U_{\debug n+1}  \|_2^2 \Big\} \\
		& \leq -2 \alpha_{\debug n+1} \langle  \v(\debug X_\debug n) + B_{\debug n} + U_{\debug n+1}, \debug X_{\debug n} - \debug x^*  \rangle +
		3\alpha_{\debug n+1}^2 \Big\{C_2 + \|B_{\debug n}\|_2^2+ \|U_{\debug n+1}  \|_2^2 \Big\}.
		\end{split}
		\end{equation}
		
		For contradiction purposes assume $\debug X_{\debug n}$ enters $\mathcal{N}(\feas^{\ast}, \epsilon)$ only a finite number of times with positive probability. By starting the sequence at a later index if necessary, we can without loss of generality
		$\debug X_{\debug n}$ never enters $\mathcal{N}(\debug x^*, \epsilon)$ with positive probability.  Then on this event (of $\debug X_{\debug n}$ never entering $\mathcal{N}(\feas^{\ast}, \epsilon)$), we have $ \langle  \v(\debug X_\debug n), \debug X_{\debug n} - \debug x^*  \rangle \ge a > 0$ as before. Telescoping Equation~\eqref{eq:energy_changes} then yields:
		\begin{equation}\label{eq:telescopings}
		\begin{split}
		& -\infty <  - \fench_{\debug x^*}(\debug Y_0)  \leq \fench_{\debug x^*}(\debug Y_{\debug n+1}) - \fench_{\debug x^*}(\debug Y_0)  = \sum_{r=0}^{\debug n} \{\fench_{\debug x^*}(\debug Y_{r+1}) - \fench_{\debug x^*}(\debug Y_{r})\}  \\
		& \leq -2\sum_{r=0}^{\debug n} \alpha_{\debug n+1} \langle  \v(\debug X_\debug n) + B_{\debug n} + U_{\debug n+1}, \debug X_{\debug n} - \debug x^*  \rangle + 3\sum_{r=0}^{\debug n} \alpha_{\debug n+1}^2 \Big\{C_2 + \|B_{\debug n}\|_2^2+ \|U_{\debug n+1}  \|_2^2 \Big\}\\
		& \leq -2\sum_{r=0}^{\debug n} \alpha_{\debug n+1} \Big\{a +\langle B_{\debug n} + U_{\debug n+1}, \debug X_{\debug n} - \debug x^*  \rangle \Big\}+
		3\sum_{r=0}^{\debug n} \alpha_{\debug n+1}^2 \Big\{C_2 + \|B_{\debug n}\|_2^2+ \|U_{\debug n+1}  \|_2^2 \Big\}\\
		& \to -\infty \text{ a.s. } \text{ as } \debug n \to \infty.
		\end{split}
		\end{equation}
		
		We justify the last-line limit of Equation~\eqref{eq:telescopings} by looking at each of its components in turn:
		\begin{enumerate}
			\item Since $\sum_{r=0}^{\debug n} \alpha_{\debug n+1}^2 < \infty$, and per the previous step, $\lim_{\debug n \to \infty} \|B_{\debug n}\|_2^2 = 0, \text{ a.s.}$, we have $3\sum_{r=0}^{\infty} \alpha_{\debug n+1}^2 \Big\{C_2 + \|B_{\debug n}\|_2^2\Big\} = C, \text{ a.s.}$,
			for some constant $C < \infty$.
			
			\item $\sum_{r=0}^{\debug n} \alpha_{\debug n+1}^2 \|U_{\debug n+1}  \|_2^2$ is submartingale that is $L_1$ bounded since:
			\begin{equation}
			\begin{split}
			&\sup_{\debug n}\exof{\sum_{r=0}^{\debug n} \alpha_{\debug n+1}^2 \|U_{\debug n+1}  \|_2^2} \leq \sup_{\debug n} \sum_{r=0}^{\debug n} \alpha_{\debug n+1}^2 \exof{\|U_{\debug n+1}  \|_2^2} \leq \sup_{\debug n} \sum_{r=0}^{\debug n} \alpha_{\debug n+1}^2 \exof{\|U_{\debug n+1}  \|_2^2} \\
			& = 
			\sup_{\debug n} \sum_{r=0}^{\debug n} \alpha_{\debug n+1}^2 \exof{\| \nabla \debug F(\debug X_{s(\debug n)},\debug\omega_{\debug n+1}) -  \v(\debug X_{s(\debug n)})       \|_2^2} 
			\notag\\
			&\leq 2\sup_{\debug n} \sum_{r=0}^{\debug n} \alpha_{\debug n+1}^2\Big\{ \exof{\|\grad \debug F(\debug X_{s(\debug n)}, \debug\omega_{\debug n+1}) \|_2^2} + \exof{\|\grad \debug f(\debug X_{s(\debug n)}) \|_2^2} \Big\} \\
			& \leq 2\sup_{\debug n} \sum_{r=0}^{\debug n} \alpha_{\debug n+1}^2\Big\{ \sup_{x \in \mathcal{\debug X}}\exof{\|\grad \debug F(x, \debug\omega) \|_2^2} + \sup_{x \in \mathcal{\debug X}}\|\grad \debug f(x) \|_2^2 \Big\} \\
			& \leq  2\sup_{\debug n} \sum_{r=0}^{\debug n} 2C_2 \alpha_{\debug n+1}^2 < \infty.
			\end{split}
			\end{equation}
			
			Consequently, by martingale convergence theorem (Lemma~\ref{thm:mg_convergence2} by taking $p=1$), $3\sum_{r=0}^{\debug n} \alpha_{\debug n+1}^2 \|U_{\debug n+1}  \|_2^2 \to R, \text {a.s.}$, for some
			random variable $R$ that is almost surely finite (in fact $\exof{|R|} < \infty$).
			
			\item Since $\|B_{\debug n}\|_2$ converges to $0$ almost surely, its average also converges to $0$ almost surely: $$\sum_{\debug n= 0}^{\infty} \frac{\debug\alpha_{\debug n+1}\|B_{\debug n}\|_2}{\sum_{r=1}^{\debug n} \debug\alpha_{r+1}} = 0, \text{ a.s.},$$
			there by implying that $$\sum_{\debug n= 0}^{\infty} \frac{\debug\alpha_{\debug n+1}\langle B_{\debug n}, \debug X_{\debug n} - \debug x^*\rangle}{\sum_{r=1}^{\debug n} \debug\alpha_{r+1}} = 0, \text{ a.s.},$$
			since $|\langle B_{\debug n}, \debug X_{\debug n} - \debug x^*\rangle| \leq \|B_{\debug n}\|_2\|\debug X_{\debug n} - \debug x^*\|_2 \leq C_4 \|B_{\debug n}\|_2$.
			
			In addition,  $\debug\alpha_{\debug n+1}\langle U_{\debug n+1}, \debug X_{\debug n} - \debug x^*  \rangle$ is a martingale difference that is $L_2$ bounded because $\debug\alpha_{\debug n+1}$ is square summable and
			$$\exof{\|\langle U_{\debug n+1}, \debug X_{\debug n} - \debug x^*  \rangle\|_2^2} \leq \exof{\|U_{\debug n+1}\|_2^2 \|\debug X_{\debug n} - \debug x^* \|_2^2} \leq C_4 \exof{\|U_{\debug n+1}\|_2^2} \leq 4C_4C_2 < \infty.$$ 
			
			Consequently, by applying Lemma~\ref{thm:mg_convergence1} with $p = 2$ and $u_n = \sum_{r=1}^{\debug n} \debug\alpha_{r+1}$ (which is not summable), law of large number for martingales therefore implies:
			$$\sum_{\debug n= 0}^{\infty} \frac{\debug\alpha_{\debug n+1}\langle U_{\debug n+1}, \debug X_{\debug n} - \debug x^*\rangle}{\sum_{r=1}^{\debug n} \debug\alpha_{r+1}} = 0, \text{ a.s.}$$
			Combining the above two limits, we have
			$$\lim_{\debug n \to \infty} 
			\frac{\sum_{r=0}^{\debug n} \alpha_{\debug n+1}\langle B_{\debug n} + U_{\debug n+1}, \debug X_{\debug n} - \debug x^*  \rangle}{\sum_{r=0}^{\debug n} \debug\alpha_{r+1}} = 0, \text { a.s.}$$
			
			Consequently, $-\sum_{r=0}^{\debug n} \alpha_{\debug n+1} \Big\{a +\langle B_{\debug n} + U_{\debug n+1}, \debug X_{\debug n} - \debug x^*  \rangle \Big\} =-\{\sum_{r=0}^{\debug n} \alpha_{\debug n+1}\} \Big\{a +  \frac{\sum_{r=0}^{\debug n} \alpha_{\debug n+1}\langle B_{\debug n} + U_{\debug n+1}, \debug X_{\debug n} - \debug x^*  \rangle}{\sum_{r=0}^{\debug n} \debug\alpha_{r+1}}\Big\} \to -\infty$,
			as $\debug n \to -\infty$.
		\end{enumerate}
	\end{enumerate}


\subsubsection{Proof of Lemma~\ref{lem:odeprop}}


	The first claim follows by computing the derivative of the energy function with respect
	to time (for notational simplicity, here we just use $\debug y(t)$ to denote $\debug P(t, y))$:
	\begin{equation}
	\label{eq:dFench}
	\begin{split}
	&\frac{d}{dt} \debug E_{\debug x^{\ast}}(y(t)) =\frac{d}{dt} \Big\{ \|\debug x^{\ast}\|^2_2 -  \|\mirror(\debug y(t))\|_2^2 + 
	2\langle\debug y(t), \mirror(\debug y(t)) - \debug x^{\ast}\rangle \Big\} \\
	& = \frac{d}{dt} \Big\{  - \|\mirror(\debug y(t)) - \debug y(t)\|_2^2 + \|\debug y(t)\|_2^2+
	2\langle\debug y(t),  - \debug x^{\ast}\rangle \Big\} \\
	& =    2\langle \dot{\debug y}(t), \mirror(\debug y(t)) - \debug y(t)\rangle + 2\langle \debug y(t), \dot{\debug y}(t)\rangle+
	2\langle\dot{\debug y}(t),  - \debug x^{\ast}\rangle \Big\} \\
	& =    - 2\langle \grad\debug f(\debug x(t)), \debug x(t) - \debug y(t)\rangle - 2\langle \grad\debug f(\debug x(t)),  \debug y(t)\rangle-
	2\langle \grad\debug f(\debug x(t)),  - \debug x^{\ast}\rangle \Big\} \\
	&= -\braket{\grad\debug f(\debug x(t))}{x(t) - \debug x^{\ast}}
	\leq 0,
	\end{split}
	\end{equation}
	where the last inequality is strict unless $\mirror(y(t)) = x(t) = \debug x^{\ast}$. Take the infimum over $\debug x^{\ast}$ then yields the result: in particular, if $\mirror(y(t)) = x(t) \notin \feas^{\ast}$, then $\frac{d}{dt} \debug E_{\debug x^{\ast}}(y(t)) \leq -\epsilon < 0, \forall \debug x^{\ast} \in \feas^{\ast}$, hence yielding the strict part of the inequality.
	Note also that even though $\mirror(\debug y(t)) - \debug y(t)$ is not differentiable, 
	$\|\mirror(\debug y(t)) - \debug y(t)\|_2^2$ is; and in computing its derivative, we applied the envelope theorem as given in Lemma~\ref{lem:env}.
	
	
	For the second claim, consider all $\debug y$ that satisfy $\fench (\debug P(t, \debug y)) >  \frac{\delta}{2}$. Fix any $\debug x^{\ast} \in \feas^{\ast}$.
	By the monotonicity property in the first part of the lemma, it follows that
	$\fench_{\debug x^*}(\debug P(s, \debug y)) >  \frac{\delta}{2}, \forall 0\leq s \leq t$.
	Consequently, $P(s, y)$ must be outside some $\epsilon$ neighborhood of $\feas^{\ast}$ for $0 \leq s \leq t$, for otherwise, $\fench_{\debug x^*} (\debug P(t, \debug y))$ would be $0$ for at least some $\debug x^{\ast} \in \feas^{\ast}$, which is a contradiction to
	$\fench (\debug P(t, \debug y)) >  \frac{\delta}{2}$.
	
	This means that there exists some positive constant $a(\delta)$ such that $\forall 0\leq s \leq t, \forall \debug x^{\ast} \in \feas^{\ast}$:
	\begin{equation}
	\label{eq:dFench-far}
	\frac{d}{ds} \debug E_{\debug x^{\ast}}(\debug P(s, y))
	= -\braket{\grad\debug f(x(s))}{x(s) - \debug x^{\ast}}
	\leq -a(\delta).
	\quad
	\end{equation}
	
	Consequently, pick $T(\delta) = \frac{\delta}{2a(\delta)}$,
	Equation~\eqref{eq:dFench-far} implies that for any $t > T(\delta)$:
	\begin{equation}
	\label{eq:SA6}
	\debug E_{\debug x^{\ast}}(\debug P(t, y)) \leq \debug E_{\debug x^{\ast}}(\debug P(T(\delta), y)) 
	\leq \debug E_{\debug x^{\ast}}(y) - T(\delta) a(\delta)
	\leq \debug E_{\debug x^{\ast}}(y) - \frac{\delta}{2}.
	\end{equation}
	Taking the supremum over $\debug x^{\ast} \in \feas^{\ast}$ then yields:
	\begin{equation}\label{eq:SA7}
	\debug E(\debug P(t, y))
	\leq \debug E(y) - \frac{\delta}{2}.
	\end{equation}
	Since Equation~\eqref{eq:SA7} is true for any $\debug y$, taking $\sup$ over $\debug y$ establishes the claim.
	\hfill$\blacksquare$

\subsection{Proof of Lemma~\ref{lem:apt}}
First, from previous analysis, we already know 
that $U_{\debug n+1}$ is a martingale difference sequence with $\sup_{\debug n}\exof{\|U_{\debug n+1}\|_2^2} < \infty$ and that $\debug\alpha_{\debug n}$ is a square-summable sequence. Recall also $\lim_{\debug n \to \infty} B_{\debug n} = 0$ \textup(a.s.\textup)\xspace.
Now, by combining Proposition 4.1 and Proposition 4.2 in \citet{Ben99},  
we obtain the following sufficient condition for APT:

The affine interpolation curve of the iterates generated by the difference equation
	$\debug Y_{\debug n+1} = \debug Y_{\debug n} - \debug\alpha_{\debug n+1}\{G(\debug X_\debug n) + U_{\debug n+1}\}$ is an \ac{APT} for the solution to the ODE $\dot{\debug y} = -G(\debug y)$ if the following three conditions \textbf{all} hold:
	\vspace{2ex}
	\begin{enumerate}

		\item $G$ is Lipschitz continuous and bounded.
			\footnote{This condition is also sufficient for the existence and uniqueness of the ODE solution}.
		\item  $U_{\debug n+1}$ is a martingale difference sequence with $\sup_{\debug n}\exof{\|U_{\debug n+1}\|_2^p} < \infty$
		and $\sum_{\debug n = 0}^{\infty} \debug\alpha_{\debug n+1}^{1 +\frac{p}{2}} < \infty$ for some $p >  2$.
	\end{enumerate}

\vspace{2ex}
By using a similar analysis as in~\citet{Ben99} (which we omit here due to space limitation), we can show that
if the following conditions all hold:
	\begin{enumerate}

	\item $G$ is Lipschitz continuous and bounded.
	\item $\lim_{\debug n \to \infty} B_{\debug n} = 0$ \textup(a.s.\textup)\xspace.
	\item  $U_{\debug n+1}$ is a martingale difference sequence with $\sup_{\debug n}\exof{\|U_{\debug n+1}\|_2^p} < \infty$
	and $\sum_{\debug n = 0}^{\infty} \debug\alpha_{\debug n+1}^{1 +\frac{p}{2}} < \infty$ for some $p >  2$,
\end{enumerate}
then, the affine interpolation curve of the iterates generated by the difference equation
$\debug Y_{\debug n+1} = \debug Y_{\debug n} - \debug\alpha_{\debug n+1}\{G(\debug X_\debug n) + B_{\debug n} + U_{\debug n+1}\}$ is an \ac{APT} for the solution to the ODE $\dot{\debug y} = -G(\debug y)$.

Take $p=2$ and recall $\v(\mirror(\cdot))$ is Lipschitz continuous and bounded.
The above list of three conditions are thus all verified, thereby yielding the result.


\subsubsection{Proof of Theorem~\ref{thm:main}}
	By Proposition~\ref{prop:recur_s}, $\debug Y_{\debug n}$ will get arbitrarily close to $\feas^{\ast}$ infinitely often.	
	It then suffices to show that, after long enough iterations, if $\debug Y_{\debug n}$ ever gets $\epsilon$-close to $\feas^{\ast}$, all the ensuing iterates will be $\epsilon$-close to $\feas^{\ast}$ almost surely.
	The way we show this ``trapping" property is to use the energy function. Specifically, we consider $\fench(\debug x^{\ast},A(t))$ and show that
	no matter how small $\epsilon$ is, for all sufficiently large $t$, if $\fench(\debug x^{\ast},A(t_0))$ is less than $\epsilon$ for some $t_0$, then $\fench(\debug x^{\ast},A(t)) < \epsilon, \forall t> t_0$.
	This would then complete the proof because $A(t)$ actually contains all the \ac{DASGD} iterates, and hence if 
	$\fench(\debug x^{\ast},A(t)) < \epsilon, \forall t> t_0$, then $\fench(\debug x^{\ast},\debug Y_{\debug n}) < \epsilon$ for all sufficiently large $\debug n$. Furthermore, since $A(t)$ contains all the iterates, the hypothesis  that `` if $\fench(\debug x^{\ast},A(t_0))$ is less than $\epsilon$ for some $t_0$" will be satisfied due to Prop~\ref{prop:recur_s}.

	
We now flesh out more details of the proof. Fix any $\epsilon > 0$.
	Since $A(t)$ is an asymptotic pseudotrajectory for $\debug P$,  we have:
	\begin{equation}
	\lim_{t \to \infty} \sup_{0 \leq h \leq T} \dnorm{\debug Y(t+h) - \debug P(h,\debug Y(t))}
	= 0.
	\end{equation}
	Consequently, for any $\delta > 0$,	there exists some $\tau(\delta,T)$ such that $\dnorm{\debug Y(t+h) - \debug P(h,\debug Y(t))} < \delta$ for all $t\geq\tau$ and all $h \in [0, T]$. We therefore have the following 
	chain of inequalities:
	\begin{flalign}
	&\debug E_{\debug x^{\ast}}(A(t+h))
	= \debug E_{\debug x^{\ast}}(\debug P(h,A(t)) + A(t+h) -\debug P(h,A(t)))\\
	&\leq \debug E_{\debug x^{\ast}}(\debug P(h,A(t)))
	+ \braket{A(t+h) - \debug P(h,A(t))}{\mirror(\debug P(h,A(t))) - \debug x^{\ast}}
	+ \frac{1}{2} \dnorm{A(t+h) -\debug P(h,A(t))}^{2} 
	\notag\\
	&\leq \debug E_{\debug x^{\ast}}(\debug P(h,A(t)))
	+ C_4\delta
	+ \frac{1}{2}\delta^2
	= \debug E_{\debug x^{\ast}}(\debug P(h,A(t))) + \frac{\varepsilon}{2},
	\end{flalign}
	where in the last step we have choosen $\delta$ small enough such that $ C_4\delta
	+ \frac{1}{2}\delta^2 =\frac{\varepsilon}{2}$.
	
	Now by Proposition~\ref{prop:recur_s}, there exists some $\tau_{0}$ such that
	$\debug E_{\debug x^{\ast}}(A(\tau_{0}))
	< \frac{\varepsilon}{2}.$
	Our goal is to establish that $\debug E_{\debug x^{\ast}}(A(\tau_{0} + h)) < \varepsilon$ for all $h\in[0, \infty)$.
	To that end, partition the $[0, \infty)$ into disjoint time intervals of the form $[(n-1) T_{\varepsilon},nT_{\varepsilon})$ for some appropriate $T_{\varepsilon}$.
	By \cref{lem:odeprop}, we have:
	\begin{equation}
	\label{eq:bound1}
	\debug E_{\debug x^{\ast}}(\debug P(h,A(\tau_{0})))
	\leq \debug E_{\debug x^{\ast}}(\debug P(0, A(\tau_{0})))
	= \debug E_{\debug x^{\ast}}(A(\tau_{0})) <\frac{\varepsilon}{2}
	\quad
	\text{for all $h \geq 0$}.
	\end{equation}
	Consequently:
	\begin{equation}
	\label{eq:bound2}
	\debug E_{\debug x^{\ast}}(A(\tau_{0}+h))
	< \debug E_{\debug x^{\ast}}(\debug P(h,A(\tau_{0}))) + \frac{\varepsilon}{2}
	< \frac{\varepsilon}{2} + \frac{\varepsilon}{2}
	= \varepsilon,
	\end{equation}
	where the last inequality is a consequence of \eqref{eq:bound1}.
	
	Now, assume inductively that \cref{eq:bound2} holds for all $h \in [(n-1)T_{\varepsilon}, nT_{\varepsilon})$ for some $n\geq1$.
	Then, for all $h\in [(n-1)T_{\varepsilon}, nT_{\varepsilon})$, we have:
	\begin{flalign}
	\label{eq:bound3}
	&\debug E_{\debug x^{\ast}}(A(\tau_{0} + T_{\varepsilon} + h))
	< \debug E_{\debug x^{\ast}}(\debug P(T_{\varepsilon},A(\tau_{0}+h))) + \frac{\varepsilon}{2}
	\leq \max\braces*{\frac{\varepsilon}{2}, \debug E_{\debug x^{\ast}}(A(\tau_{0}+h)) -\frac{\varepsilon}{2}} + \frac{\varepsilon}{2}
	\notag\\
	&\leq \frac{\varepsilon}{2} + \frac{\varepsilon}{2}
	= \varepsilon.
	\end{flalign}
	Consequently, \cref{eq:bound2} holds for all $h \in [nT_{\varepsilon}, (n+1)T_{\varepsilon})$.
This completes the induction. Taking the infimum over $\debug x^{\ast}$ then completes our proof.
	





\subsection{Our Contributions and Related Work}
Our aim in this paper is to establish the convergence of \ac{DASGD} in the presence of large, superlinear delays, in as wide a class of objectives as possible and in the presence of constraints where efficient projection can be performed.
To that end, we focus on the following classes of problems, where different convergence results can be obtained:

\paragraph{General non-convex objectives\afterhead}
We first consider the class of general smooth non-convex functions, with no structural assumption on the objectives. 
In this (difficult) case, \citet{TBA86} showed that, under \emph{sublinear} delays, \ac{DASGD} converges to a level set of the objective which contains a critical point with probability $1$;
in particular, if every such point is a global minimizer (e.g., if the problem is pseudo-convex) and the method is run with an $\Omega(1/n)$ step-size schedule, \ac{DASGD} converges to the problem's solution set.
More recently, \citet{lian2015asynchronous} derived an estimate for the rate of convergence of the surrogate length $\debug n^{-1} \sum_{\debug k=1}^{\debug n} \exof{\norm{\grad \debug f(\debug X_{\debug k})}_{2}^{2}}$ as $\debug n\to\infty$ under the assumption that the delays affecting the algorithm are \emph{bounded}.
Our first contribution is to show that these assumptions on the delays are not needed:
specifically, as we show in \cref{thm:nonconvex}, by tuning the algorithm's step-size appropriately, it is possible to retain this convergence guarantee,
even if delays grow as polynomials of arbitrary degree.


\paragraph{Variationally coherent objectives\afterhead}
Albeit directly applicable to general non-convex stochastic optimization, \cref{thm:nonconvex} only guarantees convergence to stationarity in the mean square sense;
to ensure global optimality, stronger structural assumptions on the objectives must be imposed.
The ``gold standard'' of such assumptions (and by far the most widely studied one) is convexity:
in the context of distributed stochastic convex optimization, recent works by \citet{NIPS2011_4247} and \citet{recht2011hogwild}, have established convergence for \ac{DASGD} under \emph{bounded} delays for each of the two distributed computing architectures, while \citet{chaturapruek2015asynchronous} and \citet{feyzmahdavian2016asynchronous} extended the bounded delays assumption to a setting with finite-mean \acs{iid} delays. To go beyond this framework, we focus the class of \emph{mean variationally coherent} optimization problems~\cite{zhou2017stochastic, zhou2020}, which includes pseudo-, quasi- and/ star-convex problems, as well as many other classes of functions with non-convex profiles.
Our main result here is that, in such problems, the global state parameter $\debug X_{\debug n}$ of \ac{DASGD} converges to a global minimum with probability $1$, even when the delays between gradient updates and requests grow at a polynomial rate (and this, without any distributional assumption on how the underlying delays are generated).

To go beyond this framework, we focus on a class of unimodal problems, which we call \emph{variationally coherent}, and which properly includes all pseudo-, quasi- and/ star-convex problems,
as well as many other classes of functions with highly non-convex profiles.
Our main result here may be stated as follows:
in stochastic variationally coherent problems, provided a lazy descent scheme is used (akin to dual averaging) to mesh with the constraint set in the distributed procedure, the global state parameter $\debug X_{\debug n}$ of \ac{DASGD} converges to a global minimum with probability $1$, even when the delays between gradient updates and requests grow at a polynomial rate (and this, without any distributional assumption on how the underlying delays are generated).

This result extends the works mentioned above in several directions:
specifically, it shows that
\begin{enumerate}
[\indent 1.]
\item
Convexity is not required to obtain global convergence results.
\item 
Constraints do not hinder almost sure convergence under a suitable lazy projection scheme.
\item
The robustness of \ac{DASGD} is guaranteed even under large, superlinear delays.
\end{enumerate}
We find these outcomes particularly appealing because, coupled with the existing rich literature on the topic, they help explain and reaffirm the prolific empirical success of \ac{DASGD} in large-scale machine learning problems, and offer concrete design insights for fortifying the algorithm's distributed implementation against delays and asynchronicities.



\paragraph{Techniques\afterhead}
Our analysis relies on techniques and ideas from stochastic approximation and \mbox{(sub-)}mar\-tin\-gale convergence theory.
A key feature of our approach is that, instead of focusing on the discrete-time algorithm, we first establish the convergence of an underlying, deterministic dynamical system by means of a particular energy (Lyapunov) function which is decreasing along continuous-time trajectories and ``quasi-decreasing'' along the iterates of \ac{DASGD}.
To control this gap, we connect the continuous- and discrete-time frameworks via the theory of \acdefp{APT}, as pioneered by \citet{BH96}.
By itself, the \ac{APT} method does not suffice to establish convergence under delays.
However, if the step-size of the method is chosen appropriately (following a quasi-linear decay rate for polynomially growing delays), it is possible to leverage $L^{p}$ martingale tail convergence results to show that the problem's solution set is recurrent under \ac{DASGD}.
This, combined with the above, allows us to prove our core convergence results.

Even though the \ac{ODE} approximation of discrete-time Robbins\textendash Monro algorithms has been widely studied in control and optimization theory \citep{kushner2013stochastic, ghadimi2013stochastic}, transferring the convergence guarantees of an \ac{ODE} solution trajectory to a discrete-time algorithm is a fairly subtle affair that must be done on a case-by-case basis. 
Further, even if this transfer is complete, the results typically have the nature of convergence-in-probability:
almost-sure convergence is usually much harder to obtain \cite{Bor08}.
Specifically, exisiting stochastic approximation results
cannot be applied to our setting because
\begin{inparaenum}
[\itshape a\upshape)]
\item
the non-invertibility of the projection map makes the underlying dynamical system on the problem's feasible region non-autonomous (so \ac{APT} results do not apply);
\item
unbounded delays only serve to aggravate this issue, as they introduce a further disconnect between the \ac{DASGD} algorithm and its continuous-time version.
\end{inparaenum}
To control the discrepancy between discrete and continuous time requires a more fine-grained analysis, for which we resort to a sharper law of large numbers for $L^{p}$-bounded martingales.
Finally, we also mention that the recent work of \citet{lian2016comprehensive} has also considered distributed zeroth-order methods (where only the function value, rather than the gradient is available)
and used techniques that are different from gradient-based analyses. 


\section{Introduction}
\label{sec:introduction}
\input{Introduction}


\section{Problem Setup}
\label{sec:model}
\input{Model}


\section{General Nonconvex Objectives}
\label{sec:nonconvex}
\input{Nonconvex}


\section{Variationally Coherent Problems}
\label{sec:vc}
\input{VC}


\section{Discussion}
\label{sec:discussion}
\input{Discussion}


\numberwithin{lemma}{section}	
\numberwithin{proposition}{section}	
\numberwithin{corollary}{section}	
\numberwithin{equation}{section}	

\subsection{SGD on Master-Worker Systems}\label{subsec:master_slave}
Here we consider the first distributed computing architecture: the master-worker system. 
The standard way of deploying stochastic gradient descent in such systems \textendash\ and that which we adopt here \textendash\ is for the workers
to asychronously compute stochastic gradients and then send them to the master,%
\footnote{As alluded to before, in machine learning applications, this is done by sampling a subset of the training data, computing the gradient for each datapoint and averaging over all datapoints in the sample.}
while the master updates the global state of the system and pushes the updated  state back to the workers (\cite{NIPS2011_4247,lian2015asynchronous}).
This process is presented in \cref{alg:MSSGD}.


\begin{algorithm}[htbp]
	\caption{Running SGD on a Master-Worker Architecture} 
	\label{alg:MSSGD}
	\begin{algorithmic}[1]
		\REQUIRE $1$ Master and $K$ workers, $k=1,\dotsc,K$
		\STATE Each worker is seeded with the same inital iterate
		\REPEAT 
		\STATE \textbf{Master}:
		
		(a) Receive a stochastic gradient from worker $k$
		
		(b) Update current iterate.
		
		(c) Send updated iterate to worker $k$
		
		\STATE \textbf{Workers}:
		
		(a) Receive iterate
		
		(b) Compute an i.i.d. stochastic gradient (at the received iterate)
		
		(c) Send the computed gradient to master
		
		
		\UNTIL{end}
	\end{algorithmic}
\end{algorithm}

Due to the distributed nature of the master-worker system, 
a gradient received by the master on any given iteration can be stale. As a simple example, consider a fully coordinated update scheme where each worker sends the computed gradient to and receives the updated iterate from the master following a round-robin schedule. In this case, each worker's gradient is received with a delay exactly equal to $K-1$ ($K$ is the number of workers in the system), because by the time the master receives worker $K$'s computed gradient, the master has already applied $K-1$ gradient updates from workers $1$ to $K-1$ (and since the schedule is round-robin, this delay of $K-1$ is true for any one of the $K$ workers).
 
However, delays can be much worse since we allow 
full asynchrony: workers can compute and send (stochastic) gradients to the master without any coordinated schedule.  
In the asynchronous setting, fast workers (workers that are fast in computing gradients) will cause disproprotionately large delays to gradients produced by slow workers (workers that are slow in computing gradients): when a slower worker has finished computing a gradient, a fast worker may have already computed and communicated many gradients to the master. 
Since the master updates the global state of the system (the current iterate), one can gain a clearer representation of this scheme by looking at the master's update. This is given in \cref{subsec:unify}.

\subsection{SGD on Multi-Processor Systems with Shared Memory}\label{subsec:multi_core}

Here we consider the second distributed computing architecture: multi-processor system with shared memory. 
In this architecture, all processors can access a global, shared memory, which holds all the data needed for computing a (stochastic) gradient, as well as the current iterate (the global state of the system).
The standard way of deploying stochastic gradient descent in such systems (\cite{chaturapruek2015asynchronous,lian2015asynchronous}) is for each processor to independently and asychronously read the current global iterate, compute a stochastic gradient~\footnote{This is again done by sampling a subset of the training data in the global memory and computing the gradient at the iterate for each datapoint and averaging over all the comptued gradients in the sample.}, and then update the global iterate in the shared memory.
This process is given \cref{alg:MSSGD1}:


\begin{algorithm}[H]
	\caption{Running SGD on a Multi-Processor System with Shared Memory} 
	\label{alg:MSSGD1}
	\begin{algorithmic}[1]
		\REQUIRE $K$ processors and global (shared) memory.
		\STATE The initial iterate in the global memory. 
		\REPEAT 
		\STATE
		(a) Each proceesor reads the current global iterate.
		
		(b) Each processor reads data from memory and computes a stochastic gradient.
		
		(c) Each processor updates the global iterate.
		
		\UNTIL{end}
	\end{algorithmic}
\end{algorithm}

The key difference from \cref{alg:MSSGD} is that
there is no central entity that updates the global state; instead, each processor 
can both read the global state and update it. Since each processor is performing the operations
asynchronously, different processors may be reading the same global iterate at the same time. Further, the delays
in this case is again caused by the heterogeneity across different processors: if a processor is slow in computing gradients,
then by the time it finishes computing its gradient, the global iterate has been updated by other, faster processors many times over, thereby causing its own gradient stale. Here, we also adopt a common assumption that updating the global iterate is an atomic operation (and hence no two processors will be updating the global iterate at the same time). This is justified\footnote{

	On a related note, we also note that our analysis can be further extended to cases where only one variable or a small block of variables are being updated at a time. We omit this discussion because the resulting notation is quite onerous, and will obscure the main ideas behind the already complex theoretical framework developed here.} because performing gradient update is a simple arithemtic operation, and hence takes negligible time compared to reading data and computing a stochastic gradient, which is the main computational bottleneck in practice. 
However, despite a cheap computation, performing the whole gradient update (typically achieved via locking) does have overhead. In particular, a less stringent model would be only updating one coordinate at a time: this is known as an inconsistent write/read model since different processors are updating different components of the global parameters and hence can read in elements of different ages. For simplicity, we do not consider this case, as our focus in this paper  is on delays. See~\cite{lian2015asynchronous} and \cite{liu2014asynchronous}
for lucid discussions and analyses on this model.
Finally, one can gain a clearer picture of this update scheme by tracking the update to the global iterate in the shared memory. This is given in \cref{subsec:unify}.


\subsection{DASGD: A Unifying Algorithmic Representation}
\label{subsec:unify}
In this subsection, we present a unified algorithmic description,  aptly called \acdef{DASGD}, that formally captures both \cref{alg:MSSGD} and \cref{alg:MSSGD1}, where their differences are reflected in the assumptions of the meta algorithm's parameters. We start with the unconstrained case, see \cref{alg:DASGD_un}.

\begin{algorithm}[htbp]
	\caption{Distributed asynchronous stochastic gradient descent} 
	\label{alg:DASGD_un}
	\begin{algorithmic}[1]
		\REQUIRE Initial state $\debug X_{0}\in \mathbb{R}^\debug d$, step-size sequence $\debug\alpha_{\debug n}$
		\STATE $\debug n\leftarrow0$;
		\REPEAT
		\STATE $\debug X_{\debug n+1} = \debug X_{\debug n} - \debug\alpha_{\debug n+1} \nabla \debug F(\debug X_{s(\debug n)},\debug\omega_{\debug n+1})$;
		\STATE $\debug n\leftarrow \debug n+1$;
		\UNTIL{end}
		\STATE \textbf{return} solution candidate $\debug X_{\debug n}$
	\end{algorithmic}
\end{algorithm}

\begin{algorithm}[htbp]
	\caption{Distributed asynchronous stochastic gradient descent with projection} 
	\label{alg:DASGD}
	\begin{algorithmic}[1]
		\REQUIRE Initial state $\debug Y_{0}\in \mathbb{R}^\debug d$, step-size sequence $\debug\alpha_{\debug n}$
		\STATE $\debug n\leftarrow0$;
		\REPEAT
		\STATE $\debug X_{\debug n} = \mirror(\debug Y_{\debug n})$;
		\STATE $\debug Y_{\debug n+1} = \debug Y_{\debug n} - \debug\alpha_{\debug n+1} \nabla \debug F(\debug X_{s(\debug n)},\debug\omega_{\debug n+1})$;
		\STATE $\debug n\leftarrow \debug n+1$;
		\UNTIL{end}
		\STATE \textbf{return} solution candidate $\debug X_{\debug n}$
	\end{algorithmic}
\end{algorithm}


In more detail, $\debug n$ is a global counter and is incremented every time an update occurs to the current solution candidate $\debug X_{\debug n}$ (the global iterate): in the master-worker systems, the master updates it; in the multi-processor systems, each processor updates it.
Since there are delays in both systems, the gradient applied to the current iterate $\debug X_{\debug n}$ can be a gradient associated with a previous time step. This fact is abstractly captured by Line 3 in \cref{alg:DASGD}.
In full generality, we will write $s(\debug n)$ for the iteration from which the gradient received at time $\debug n$ originated.
In other words, the delay associated with iteration $s(\debug n)$ is $\debug n - s(\debug n)$, since it took $\debug n - s(\debug n)$ iterations for the gradient computed on iteration $s(\debug n)$ to be received at stage $\debug n$. Note that $s(\debug n)$ is always no larger than $\debug n$; and if $\debug n = s(\debug n)$, then there is no delay in iteration $\debug n$.

Now, the difference between the two distributed computing archecitures is reflected in the assumption of $s(\debug n)$.
Specifically, in the master-worker systems, each $s(\cdot)$ is a one-to-one function\footnote{Except initially if all the workers have the same initial point.}, because no two workers will ever receive same iterates from the master per \cref{alg:MSSGD}.
On the other hand, in multi-processor systems, $s(\debug n)$ can be the same for different $\debug n$'s (since different processors may read the current iterate at the same time); however, it is easy to observe that the same $s$ will appear at most $K$ times for different $\debug n$'s, since there are $K$ processors in total.
As an important note, our analysis is \textbf{agnostic} to whether $s(\debug n)$ is one-to-one or not. Consequently, 
in establishing theoretical guarantees for the meta algorithm \ac{DASGD}, we obtain the same guarantees for both architectures simultaneously.

Notation-wise, we will write $d_{\debug n}$ for the delay required to compute a gradient requested at iteration $\debug n$.
This gradient is received at iteration $\debug n+d_{\debug n}$. Following this notation, the delay for a gradient received at $\debug n$ is $d_{s(\debug n)} = \debug n - s(\debug n)$. Note also we have chosen the subscript associated with $\debug\omega$ to be $\debug n+1$: we can do so because $\debug\omega_{\debug n}$'s are \textbf{iid} (and hence the indexing is irrelevant). 
Finally, in constrained optimization case (where $\mathcal{\debug X}$ is a strict subset of $\mathbb{R}^{\debug d}$), projection must be performed. This results in $\ac{DASGD}$ with projection\footnote{This tpe of projection is technically known as lazy projection.}, which is formally given in \cref{alg:DASGD}.
\subsection{Delay Assumption}
Our goal here is to establish convergence guarantees (in mean square) of \ac{DASGD} for a general non-convex objective in the presence of delays. 
In fact, a large family of unbounded delay processes can be tolerated. We state our main assumption regarding the delays and step-sizes:
\begin{assumption}
	\label{asm:delays}
	\quad
	The gradient delay process $d_{\debug n}$ and the step-size sequence $\debug\alpha_{\debug n}$ of \ac{DASGD} (\cref{alg:DASGD_un,alg:DASGD}) satisfy one of the following conditions:
	\vspace{2ex}
	\begin{enumerate}
		\addtolength{\itemsep}{1ex}
		\item
		\emph{Bounded delays:}
		$\sup_{\debug n} d_{\debug n} \le D$ for some positive number $D$
		and
		$\sum_{\debug n=1}^{\infty} \debug\alpha_{\debug n}^{2} < \infty$,
		$\sum_{\debug n=1}^{\infty} \debug\alpha_{\debug n} = \infty$.
		\item
		\emph{Sublinearly growing delays:}
		$d_{\debug n} = O(\debug n^p)$ for some $0 < p < 1$
		and
		$\debug\alpha_{\debug n} \propto 1/\debug n$.
		\item
		\emph{Linearly growing delays:}
		$d_{\debug n} = \bigoh(\debug n)$
		and
		$\debug\alpha_{\debug n} \propto 1/(\debug n \log\debug n)$.
		\item
		\emph{Polynomially growing delays:}
		$d_{\debug n} = \bigoh(\debug n^{q})$ for some $q\geq1$
		and
		$\debug\alpha_{\debug n} \propto 1/(\debug n \log\debug n \log\log\debug n)$.
	\end{enumerate}
\end{assumption}

Note that as delays get larger, we need to use less aggressive step-sizes. This is to be expected, 
because the larger the delays, the more ``averaging" one needs perform in order to remove the staleness that is caused by
the delays; and smaller step-sizes correspond to averaging over a longer horizon. This is a one of the important insights that occur throughout the paper. Another thing to note that is the Assumption~\ref{asm:delays} also highlights the quantitative relationship between the class of delays and the class of step-sizes. For instance, when the delays increase from a linear rate to a polynomial rate, only a factor of $\frac{1}{\log\log \debug n}$ needs to be added (which is effectively a constant).
From a practical standpoint, this means that a step-size on the order of $1/(\debug n \log\debug n)$ will be a good model-agnostic choice and more-or-less sufficient for almost all delay processes.

\subsection{Controlling the Tail Behavior of Second Moments}
We now turn to establish the theoretical convergence guarantees of \ac{DASGD} for general non-convex objectives.
Our first step lies in controlling the tail behavior of the second moments of the gradients that are generated from \ac{DASGD}.
By leveraging the Lipschitz continuity of the gradient, its telescoping sum, appropriate conditionings and a careful analysis of the interplay between delays and step-sizes, 
we show that (next lemma) a particularly weighted tail sum of the second moments are vanishingly small in the limit (see appendix for a detailed proof).


\begin{lemma}\label{lem:tail}\quad
	Under \cref{asm:basic,asm:delays}, if $\inf_{\debug x \in \mathcal{\debug X}}\debug f(\debug x) > - \infty$, then
	\begin{equation}\label{eq:tail}
	\sum_{\debug n = 0}^{\infty} \debug\alpha_{\debug n+1}\exof{\|\grad f(X_{\debug n})\|_2^2} < \infty.
	\end{equation}
\end{lemma}

\begin{remark}\quad
	Since $\sum_{\debug n = 0}^{\infty} \debug\alpha_{\debug n+1} = \infty$, \cref{lem:tail} implies that
	$\lim\inf_{\debug n \rightarrow \infty} \exof{\|\grad f(X_{\debug n})\|_2^2} = 0$ (see Appendix A of~\cite{bertsekas1996neuro}).
Note that  the converse is not true: when a subsequence of $\exof{\|\grad f(X_{\debug n})\|_2^2}$ converges to $0$,
	the sum in Equation~\eqref{eq:tail} need not be finite. As a simple example, consider
	$\debug\alpha_{\debug n+1} = \frac{1}{\debug n}$, and
	\begin{equation}
	\exof{\|\grad f(X_{\debug n})\|_2^2} =
	\begin{cases}
	\frac{1}{\debug n}, \text{ if } \debug n = 2^k\\
	1, \text{ otherwise}.
	\end{cases}
	\end{equation}
	Then the subsequence on indicies $2^{\debug n}$ converges to $0$, but the sum still diverges.
	Consequently, \cref{lem:tail} is stronger than subsequence convergence.
\end{remark}


\subsection{Bounding the Successive Differences}

However, \cref{lem:tail} is still not strong enough to guarantee that 
$\lim_{\debug n \rightarrow \infty} \exof{\|\grad f(X_{\debug n})\|_2^2} = 0$.
This is because the convergent sum given in Equation~\eqref{eq:tail} only limits the tail growth somewhat, but 
not completely. To demonstrate this point, let $c_t$ be the following boolean variable: 
\begin{equation}
c_{\debug n} =
\begin{cases}
1, \text{ if $\debug n$ contains the digit $9$ in its decimal expansion,} \\
0, \text{ otherwise}.
\end{cases}
\end{equation}
For instance, $c_9 = 1$, $c_{11} = 0$.
Now define 
$\debug\alpha_{\debug n+1} = \frac{1}{\debug n}$, and
$	\exof{\|\grad f(X_{\debug n})\|_2^2}=
\begin{cases}
\frac{1}{\debug n}, \text{ if } c_{\debug n} = 1\\
1, \text{ if } c_{\debug n} = 0.
\end{cases}$
As first-year Berkeley Math PhDs delightfully found out during-- or painfully found out after -- their qualifying exam,  $\sum_{\debug n = 1}^{\infty} \debug\alpha_{\debug n+1} 	\exof{\|\grad f(X_{\debug n})\|_2^2}< \infty$ (see Problem 1.3.24 in~\cite{de2012berkeley}), even though the limit $\exof{\|\grad f(X_{\debug n})\|_2^2}$ does not exist. This indicates that to obtain convergence of $\exof{\|\grad f(X_{\debug n})\|_2^2}$, we need to impose more stringent conditions to ensure its sufficient tail decay. Note that one issue that is revealed by the above counter-example is that the distance between two successive terms is always bounded away from $0$, no matter how large $\debug n$ is. This obviously makes it possible
for convergence to occur: a necessary condition for convergence is that the difference converges to $0$.
Note that intuitively, this pathological case cannot occur for gradient descent, because the step-size is shrinking to $0$, hence making the successive difference converge to $0$ (at least in expectation).
Consequently, 
we can bound the difference between every two successive terms in terms of a decreasing sequence that is converging to $0$.
This ensures that $\exof{\|\grad f(X_{\debug n})\|_2^2}$ cannot change two much from iteration to iteration. Further, the change between two successive terms will be vanishing.
This result is formalized in the following lemma (the proof is given in the appendix):

\begin{lemma}\label{lem:diff}\quad
	Under \cref{asm:basic,asm:delays},
	there exists a constant $C > 0$ such that for every $\debug n$, 
	$$\Big|\exof{\|\grad f(X_{\debug n+1})\|_2^2} - 	\exof{\|\grad f(X_{\debug n})\|_2^2}\Big| \le C \alpha_{\debug n+1}.$$
\end{lemma}



\subsection{Main Convergence Result}

Putting the above two characterizations together, we obtain:
\begin{theorem}\label{thm:nonconvex}\quad
	Let $\debug X_{\debug n}$ be the \ac{DASGD} iterates generated from \cref{alg:DASGD_un}.
	Under \cref{asm:basic,asm:delays}, if $\inf_{\debug x \in \mathcal{\debug X}}\debug f(\debug x) > - \infty$, then
	\begin{equation}\label{eq:tail}
	\lim_{\debug n \to \infty} \exof{\|\grad f(X_{\debug n})\|_2^2} = 0.
	\end{equation}
\end{theorem}

\begin{remark}\quad
Three remarks are in order here. 
First, note that the condition $\inf_{\debug x \in \mathcal{\debug X}}\debug f(\debug x) > - \infty$ means that the optimization problem~\eqref{eq:opt}
has a solution. This is necessary, because otherwise, a stationary point may not exist in the first place, and \ac{DASGD} (or even simple gradient descent) may continue decrease the objective value \textit{ad infinitum}.
Note that since $\debug f$ is smooth, a minimum point is necessarily a stationary point.

Second, the above convergence is a fairly strong characterization of the fact that the gradient vanishes. In particular,
it means that the gradient of the \ac{DASGD} iterates converges to $0$ in mean square. Consequently, this implies that the norm of the gradient vanishes in expectation, and that the gradient converges to $0$ with high probability. Note further that 
if we strengthen Assumption~\ref{asm:basic} to require that all stochastic gradients are bounded almost surely (as opoosed to just bounded in second moments), then a similar analysis ensures almost sure convergence of
$\|\grad f(X_{\debug n})\|_2$. We omit the details due to space limitation.
Finally, that \cref{thm:nonconvex} is a direct consequence of \cref{lem:tail,lem:diff}
is a simple excercise in elementary series theory (in particular, \cref{lem:A3}).
\end{remark}

\subsection{Mean Variational Coherence}\label{subsec:VC}

\begin{assumption}
	\label{asm:VC}\quad
	The optimization problem \eqref{eq:opt} is called \emph{variationally coherent in the mean} if
	\begin{equation}
	\label{eq:VC}
	\tag{VC}
	\exof{\braket{\grad\debug F(x;\debug\omega)}{x-\debug x^{\ast}}}
	\ge 0,


	\end{equation}
	for all $\debug x^{\ast}\in\feas^{\ast}$ and all $x \in \mathcal{\debug X}$ with equality only if $x\in\feas^{\ast}$.
\end{assumption} 

Note that we do not impose the ``if" condition for equality: if $x \in \feas^{\ast}$, then the equality may or may not hold.
By \cref{asm:basic}, we can interchange expectation and differentiation in \eqref{eq:VC} to obtain
\begin{equation}
\braket{\nabla\debug f(x)}{x-\debug x^{\ast}}
\ge 0,
\end{equation}
for all $x\in\mathcal{\debug X}$, $\debug x^{\ast}\in\feas^{\ast}$.
As a result, mean variational coherence can be interpreted as an averaged coherence condition for the deterministic optimization problem with objective $\debug f(x)$. Mean variationally coherent optimization problems include convex programs, pseudo-convex programs, non-degenerate quasi-convex programs and star-convex programs as special cases.  For instance, if $g$ is star-convex, then $\braket{\nabla \debug f(x)}{x - \debug x^{\ast}} \geq \debug f(x) - \debug f(\debug x^{\ast})$ for all $x\in\mathcal{\debug X}$, $\debug x^{\ast}\in\feas^{\ast}$.
This is easily seen to be a special case of variational coherence because
$\braket{\nabla \debug f(x)}{x - x^*} \ge \debug f(x) - \debug f(x^*) \ge 0$, with the last inequality strict unless $x \in \mathcal{X}^*$. Note that star-convex functions contain convex functions as a subclass (but not necessarily pseudo/quasi-convex functions). See~\cite{zhou2017stochastic, zhou2020} for a more detailed discussion on why these various classes of optimization problems are special cases.
\cref{fig:plot-global} also provides two more elaborate examples of a variationally coherent optimization problem that are not quasi-convex.
Put together, these examples indicate that variationally coherent objectives can have highly non-convex profiles. 
Nevertheless, for the class of variationally coherent functions, it is possible to establish almost sure convergence guarantees (which we do so in the subsequent sections).



\subsection{Deterministic Analysis: Convergence to Global Optima}
\label{sec:deterministic}

To streamline our presentation and build intuition along the way, we will begin with the deterministic case in this subsection, where there is no randomness in the calculation of a gradient update. 
In this case, \ac{DASGD} boils down to \acdef{DAGD}, as illustrated in \cref{alg:DAGD}:

\begin{algorithm}[H]
	\caption{Distributed Asynchronous Gradient Descent (\ac{DAGD})} 
	\label{alg:DAGD}
	\begin{algorithmic}[1]
		\REQUIRE Initial state $\debug y_{0}\in \mathbb{R}^{\debug d}$, step-size sequence $\debug\alpha_{\debug n}$
		\STATE $\debug n\leftarrow0$
		\REPEAT
		\STATE $\debug x_{\debug n} = \mirror(\debug y_{\debug n})$;
		\STATE $\debug y_{\debug n+1} = \debug y_{\debug n} - \debug\alpha_{\debug n+1} \nabla \debug f(\debug x_{s(n)})$;
		\STATE $\debug n\leftarrow \debug n+1$;
		\UNTIL{end}
		\STATE \textbf{return} solution candidate $\debug x_{\debug n}$
	\end{algorithmic}
\end{algorithm}

\subsubsection{Energy Function}
\label{subsec:energy}

We start with an energy function that measures on how ``optimal" the dual variable $\debug y$ is: the smaller the energy (i.e. the closer it is to $0$), the better the dual variable. 

\begin{definition}
	\label{def:energy}\quad
	Let $\debug x^{\ast} \in \feas^{\ast}$. Define the energy function $\fench\colon\mathbb{R}^{\debug d} \to \mathbb{R}$ as follows:
	\begin{equation}
	\fench(\debug y)
	= \inf_{\debug x^{\ast} \in \feas^{\ast}} \fench_{\debug x^{\ast}}(\debug y), \text{ where }
	\fench_{\debug x^{\ast}}(\debug y)
	= \norm{\debug x^{\ast}}_{2}^{2}
	- \norm{\mirror(\debug y)}_{2}^{2}
	+  2\braket{\debug y}{\mirror(\debug y) - \debug x^{\ast}}.
	\end{equation}




\end{definition}

Note that one can think of $\fench_{\debug x^{\ast}}(\debug y)$ as the energy of $\debug y$ with respect to a fixed optimal solution $\debug x^{\ast}$, while $\fench(\debug y)$ is the best (smallest) energy for a given $\debug y$. We next characterize a few of its useful properties.
\begin{lemma}\quad
	\label{lem:energy}
	For all $\debug y\in\mathbb{R}^{\debug d}$, we have:
	\vspace{2ex}
	\begin{enumerate}
		\addtolength{\itemsep}{1ex}
		\item $\fench(\debug y) \ge 0$ with equality if and only if $\mirror(\debug y) \in \feas^{\ast}$.
		\item Let $\{\debug y_{\debug n}\}_{\debug n=1}^{\infty}$  be a given sequence.
		If $\lim_{\debug n\to\infty} \fench(\debug y_{\debug n}) = 0$, then
		$\mirror(\debug y_{\debug n}) \rightarrow \feas^{\ast}$ as\footnote{Following the convention in point-set topology, a sequence $s_{\debug n}$ converges to a set $\mathcal{S}$ if $\text{dist}(s_{\debug n}, \mathcal{S}) \rightarrow 0$, with $\text{dist}(\cdot, \cdot)$ being the standard point-to-set distance function: $\text{dist}(s_{\debug n}, \mathcal{S}) \triangleq \inf_{s \in \mathcal{S}} \text{dist}(s_{\debug n}, s)$, where $\text{dist}(s_{\debug n}, s)$ is the Euclidean distance between $s_{\debug n}$ and $s$.} $\debug n\to\infty$.
	\end{enumerate}
\end{lemma}

\begin{remark}\label{rem:prop1}\quad
	The proof is given in the appendix, but it is helpful to make a few quick remarks.
	The first statement justifies the terminology of ``energy", as $\fench(\debug y)$ is always non-negative.
	This energy function will also be the tool we use to establish an important component of the global convergence result.
	Further, given that $\fench(\debug y) \ge 0$, it should also be clear that
	when $\mirror(\debug y) \in \feas^{\ast} \subset \mathcal{\debug X}$, we can choose a particular $\debug x^{\ast} = \mirror(\debug y)$ such that $\fench(\debug y) = 0$ (but that $\fench(\debug y) = 0$ implies $\mirror(\debug y) \in \feas^{\ast}$ is less obvious). 
    Statement 2 of the lemma provides us with a way to establish convergence to optimal solutions.
	If we can show that $\fench(\debug y_{\debug n}) \to 0$, then $\debug x_{\debug n} = \mirror(\debug y_{\debug n}) \to \feas^{\ast}$.
	Nevertheless, as we shall see later, unlike many other Lyapunov functions in optimization, $\fench(\debug y_{\debug n})$ does not decrease monotonically; consequently, it is difficult to directly establish $\fench(\debug y_{\debug n}) \to 0$. In fact, a finer-grained analysis is required to characterize the convergence behavior of $\fench(\debug y_{\debug n})$ (see \cref{subsubsec:main_deter}).
\end{remark}

We also assume that the converse of Statement 2 of the above lemma holds:

\begin{assumption}\label{assump:reci}\quad
	If $\lim_{\debug n\to\infty} \fench(\debug y_{\debug n}) = 0$, then
	$\mirror(\debug y_{\debug n}) \rightarrow \feas^{\ast}$ as $n \rightarrow \infty$.
\end{assumption} 

\Cref{assump:reci} can be seen as a primal-dual analogue
of the reciprocity conditions for the Bregman divergence. This assumption usually holds (e.g. when the feasible set $\mathcal{\debug X}$ is a polytope), unless  $\mathcal{\debug X}$ is pathological. 

\begin{lemma}\label{lem:energy_prop}
	\quad Fix any $\debug x^{\ast} \in \feas^{\ast}$.
	\vspace{2ex}
	\begin{enumerate}
		\addtolength{\itemsep}{1ex}
		\item
		$\|\mirror(\debug y) - \hat{\debug y}\|_2^2 - \|\mirror(\hat{\debug y}) -\hat{\debug y}\|_2^2 \le
		\|\debug y-\hat{\debug y}\|_2^2$, for any $\debug y, \hat{\debug y} \in \mathbb{R}^d$.
		\item
		$\fench_{\debug x^{\ast}}(\debug y + \Delta \debug y) - \fench_{\debug x^{\ast}}(\debug y) \le 2\langle   \Delta \debug y, \mirror(\debug y) - \debug x^*  \rangle +
		\|\Delta \debug y\|_2^2$, for any $\debug y,\Delta \debug y \in \mathbb{R}^d$.
	\end{enumerate}
\end{lemma}

\begin{remark}\label{rem:prop2}
	\quad
	The first statement of the above lemma serves as an important intermediate step in proving the second statement, 
	and is established by leveraging the envelop theorem and several important properties of Euclidean projection.
	To see that this is not trivial, consider the quantity $\|\mirror(\debug y) - \hat{\debug y}\|_2 - \|\mirror(\hat{\debug y}) -\hat{\debug y}\|_2$,
	which we know by triangle's inequality satisfies:
	\begin{equation}\label{eq:easy}
	\begin{split}
	&\|\mirror(\debug y) - \hat{\debug y}\|_2 - \|\mirror(\hat{\debug y}) -\hat{\debug y}\|_2 \le \| \mirror(\debug y) -\mirror(\hat{\debug y}) \|_2 \le \|\debug y - \hat{\debug y}\|_2,
	\end{split}
	\end{equation}
	where the last inequality follows from the fact that projection is a non-expansive map.
	However, this inequality is not sufficient for our purposes because in quantifying the perturbation $\fench(\debug y + \Delta \debug y) - \fench(\debug y)$, we also need the squared term  $\|\Delta \debug y\|_2^2$, which is not easily obtainable from Equation~\eqref{eq:easy}.
	In fact, a tighter analysis is needed to establish that $\|\debug y -\hat{\debug y}\|_2^2$ is an upper bound on $\|\mirror(\debug y) - \hat{\debug y}\|_2^2 - \|\mirror(\hat{\debug y}) -\hat{\debug y}\|_2^2$. 
\end{remark}

\subsubsection{Main Convergence Result}\label{subsubsec:main_deter}


An intermediate result, interesting on its own and useful also as a heavy-lifting tool for our convergence analysis is then provided by the following technical result:

\begin{proposition}\quad
	\label{prop:recur_d}
	Under \cref{asm:basic,asm:VC,asm:delays,assump:reci}, \ac{DAGD} admits a subsequence $\debug x_{\debug n_{k}}$ that converges to $\feas^{\ast}$ as $k\to\infty$.
\end{proposition}

We highlight the main steps below and refer the reader to the appendix for the details:
\begin{enumerate}
	\item 
	Letting $b_{\debug n} = \v(x_{s(n)}) - \v(x_\debug n)$,
	we can rewrite the gradient update in \ac{DAGD} as:
	\begin{flalign}
	y_{\debug n+1}
	&= y_{\debug n}
	- \debug\alpha_{\debug n+1} \v(x_{s(n)})
	\notag\\
	&= y_{\debug n}
	- \debug\alpha_{\debug n+1} \v(x_\debug n)
	- \debug\alpha_{\debug n+1} \{\v(x_{s(n)}) - \v(x_\debug n)\}
	\notag\\
	&= y_{\debug n}
	- \debug\alpha_{\debug n+1} (\v(x_\debug n) + b_{\debug n}).
	\end{flalign}
	Recall here that  $s(\debug n)$ denotes the previous iteration count whose gradient becomes available only at the current iteration $\debug n$.
	By bounding the magnitude of $b_{\debug n}$ using the delay sequence through a careful analysis, we establish that under any one of the conditions in \cref{asm:delays}, $\lim_{\debug n \to \infty} \|b_{\debug n}\|_2 = 0$.
	The analysis here, particularly the one for the last three conditions, reveals the following pattern: as the magnitude of the delays gets larger and larger in the order of growth, one needs to use a more conservative step-size sequence in order to mitigate the damage done by the stale gradient information. Intuitively, smaller step-sizes are more helpful in larger delays because they carry a better ``amortization" effect that makes \ac{DAGD} more tolerant to delays.
	
	
	\item 
	With the defintion of $b_{\debug n}$, \ac{DAGD} can be written as:
	\begin{equation}
	\begin{aligned}
	\debug x_{\debug n}
	&= \mirror(\debug y_{\debug n}),\\
	\debug y_{\debug n+1}
	&= y_{\debug n} - \debug\alpha_{\debug n+1} (\v(x_\debug n) + b_{\debug n}).
	\end{aligned}
	\end{equation}
	We then use the energy function to study the behavior of $\debug y_{\debug n}$ and $\debug x_{\debug n}$. More specifically, we look at the quantity $\fench(\debug y_{\debug n+1}) - \fench(\debug y_{\debug n})$ and, using \cref{lem:energy_prop}, we bound this one-step change using the step size $\debug\alpha_{\debug n}$, the $b_{\debug n}$ sequence and 
	the defining quantity $ \langle \v(x_\debug n), \debug x_{\debug n} - \debug x^*  \rangle$ of a variationally coherent function (as well as another term that will prove inconsequential).  
	We then telescope on $\fench(\debug y_{\debug n+1}) - \fench(\debug y_{\debug n})$ to
	obtain an upper bound for $\fench(\debug y_{\debug n+1}) - \fench(\debug y_0)$.
	Since the energy function is always non-negative, 
	$\fench(\debug y_{\debug n+1}) - \fench(\debug y_0)$ is at least $- \fench(\debug y_0)$ 
	for every $\debug n$. Then, utilizing the fact that $b_{\debug n}$ converges to $0$ and that 
	$ \langle \v(x_\debug n), \debug x_{\debug n} - \debug x^*  \rangle$ is always positive (unless the iterate is exactly an optimal solution, in which case it is $0$), we show that the upper bound will approach $-\infty$ if $\debug X_{\debug n}$ only enters $\mathcal{N}(\feas^{\ast}, \epsilon)$, an open $\epsilon$-neighborhood of $\feas^{\ast}$, a finite number of times (for an arbitrary $\epsilon > 0$). This generates an immediate contradiction, and thereby establishes that $\debug X_{\debug n}$ will get arbitrarily close to $\feas^{\ast}$ for an infinite number of times.


	This then implies that there exists a subsequence of \ac{DAGD} iterates that converges to the solution set of \eqref{eq:opt}, i.e, $\debug x_{\debug n_k} \to \feas^{\ast}$ as $k \to \infty$.
\end{enumerate}

\begin{theorem}
	\label{thm:convergence}\quad
	Under \cref{asm:basic,asm:VC,asm:delays,assump:reci}, the global state variable $\debug x_{\debug n}$ of \ac{DAGD} \textpar{\cref{alg:DAGD}} converges to the solution set $\feas^{\ast}$ of \eqref{eq:opt}.
\end{theorem}

We give an outline of the proof below, referring to the appendix for the details.

Fix a $\delta > 0$.
Since $\debug x_{\debug n_k} \to \debug x^*$, as $k \to \infty$ (per \cref{prop:recur_d}),
we have $\fench(\debug y_{\debug n_k}) \to 0$ as $k \to \infty$ per \cref{lem:energy}.
So we can pick an $\debug n$ that is sufficiently large and $\fench(\debug y_{\debug n}) < \delta$.
Our goal here is to show that for $n$ large enough, once $\fench(\debug y_{\debug n}) < \delta$, it will stay so forever: $\fench(\debug y_{m}) < \delta, \forall m \ge n$.
Note in particular that the above statement is not true for all $\debug n$, but only for $\debug n$ large enough.

However, the behavior of $\fench(\debug y_{\debug n})$ is not very regular: it can certainly increase from iteration to iteration for any $\debug n$. Nevertheless, we can precisely quantify how large this increment (if any) can be. This leads us to
break it down to two cases:
\vspace{2ex}
\begin{enumerate}
	\addtolength{\itemsep}{1ex}
	\item 
	Case 1: $\fench(\debug y_{\debug n}) < \delta/2$.
	
	\item
	Case 2: $\delta/2 \leq \fench(\debug y_{\debug n}) < \delta$.
\end{enumerate}

For Case 1, we show in the appendix that 
\begin{equation}
\fench{(\debug y_{\debug n+1})} - \fench(\debug y_{\debug n}) 
\leq 2BC_4\debug\alpha_{\debug n+1}  +  2\debug\alpha_{\debug n+1}^2(C_2 + B^2),
\end{equation}
for suitable constants $B$ and $C$‘s.
Now for $\debug n$ sufficiently large, we can make the right-hand arbitrarily small, and in particular, smaller
than $\delta/2$.
This means $\fench{(\debug y_{\debug n+1})} \le \fench{(\debug y_{\debug n})} + \frac{\delta}{2} < \delta$. Consequently, in this case, the energy stays within the $\delta$ bound in the next iteration.

For Case 2, we show in the appendix that
\begin{equation}
\fench{(\debug y_{\debug n+1})} - \fench(\debug y_{\debug n})
\leq -2\debug\alpha_{\debug n+1} \bracks*{\frac{a}{2}  - \debug\alpha_{\debug n+1}(C_2 + B^2)},
\end{equation}
where $a$ is a positive constant that depends only on $\delta$.
Again, since $\debug n$ is sufficiently large, we can make $\frac{a}{2}  - \debug\alpha_{\debug n+1}(C_2 + B^2)$ positive, thereby making the right-hand side negative.
Consequently, $\fench{(\debug y_{\debug n+1})} < \fench(\debug y_{\debug n}) < \delta$.
Hence, again, the energy stays within the $\delta$ bound in the next iteration.

The key conclusion from the above is that, for large enough $\debug n$, once 
$\fench{(\debug y_{\debug n})}$ is less than $\delta$,
$\fench(\debug y_{\debug n+1})$ is less than $\delta$ as well and so are all the iterates afterwards. Since $\delta$ is arbitrary, it follows $\fench{(\debug y_{\debug n})} \to 0$, and therefore $\debug x_\debug n \to \feas^{\ast}$ by \cref{lem:energy}.

\subsection{Stochastic Analysis: Almost Sure Convergence to Global Optima}
\label{sec:stochastic}


Having established deterministic global convergence of \ac{DAGD}, we now proceed to study stochastic global convergence of \ac{DASGD}. Compared to the deterministic analysis,
the stochastic case is much more involved because randomness can lead to very volatile behavior in the presence of delays; in particular, the simple approach employed in \cref{thm:convergence} (to establish that once  $\fench{(\debug y_{\debug n})}$ is less than $\delta$, it will always remain so) no longer works. To deal with both delays and noise, a much more sophisticated analysis framework needs to be developed, which requires several news ideas.
To streamline our presentation, we break the theoretical development into four subsections, each comprising an important component and step of the overall analysis.

\subsubsection{Recurrence of \ac{DASGD}}

Our first step lies in generalizing \cref{prop:recur_d} to the stochastic case. Specifically, in the presence of noise, we show that the iterates of \ac{DASGD} 
visit any neighorhood of $\feas^{\ast}$ infinitely often almost surely.

\begin{proposition}
	\label{prop:recur_s}\quad
	Under \cref{asm:basic,asm:VC,asm:delays,assump:reci}, \ac{DASGD} admits a subsequence $\debug X_{\debug n_{k}}$ that converges to $\feas^{\ast}$ almost surely: $\debug X_{\debug n_k} \to \feas^{\ast}$ with probability $1$ as $k \to \infty$.
\end{proposition}


	We outline the two main steps of the proof below, referring the reader to the appendix for the details.
	
	\begin{enumerate}
		\item 
		We begin by rewriting the gradient update step in \ac{DASGD} as:
		\begin{flalign}
		\frac{\debug Y_{\debug n+1} - \debug Y_{\debug n}}{\debug\alpha_{\debug n+1}}
		&= -\nabla \debug F(\debug X_{s(\debug n)},\debug\omega_{\debug n+1})
		\notag\\
		&= -\v(\debug X_\debug n)
		\notag\\
		&-\bracks{\v(\debug X_{s(\debug n)})- \v(\debug X_\debug n)}
		\notag\\
		&-\bracks{\nabla \debug F(\debug X_{s(\debug n)},\debug\omega_{\debug n+1}) -  \v(\debug X_{s(\debug n)})}.
		\end{flalign}
		Letting $B_{\debug n} = \v(\debug X_{s(\debug n)})- \v(\debug X_\debug n)$ and $U_{\debug n+1} = \nabla \debug F(\debug X_{s(\debug n)},\debug\omega_{\debug n+1}) -  \v(\debug X_{s(\debug n)})$, we can rewrite the \ac{DASGD} update as
		\begin{equation}
		\label{eq:DASGD_abbre}
		\debug Y_{\debug n+1}
		= \debug Y_{\debug n} - \debug\alpha_{\debug n+1}\{\v(\debug X_\debug n) + B_{\debug n} + U_{\debug n+1}\}.
		\end{equation} 
		
		We then establish the following two facts in this step. 
		First, we verify that 
		$\sum_{r=0}^{\debug n} U_{\debug n+1}$ is a martingale adapted to $\debug Y_1, \debug Y_2 \dots, \debug Y_{\debug n+1}$, where $\{U_{\debug n+1}\}_{\debug n=0}^{\infty}$ is a $L_2$-bounded martingale difference sequence. Second, we show that $\lim_{\debug n \to \infty} \|B_{\debug n}\|_2 = 0, { a.s.}.$
		
		The second claim is done by first giving an upper bound on $\|B_{\debug n}\|_2$ by writing $\v(\debug X_{s(\debug n)})- \v(\debug X_\debug n)$ as a sum of one-step changes ($\v(\debug X_{s(\debug n)})- \v(\debug X_{s(\debug n)+1}) + \v(\debug X_{s(\debug n)+1}) - \dots + \v(\debug X_{\debug n-1})- \v(\debug X_\debug n)$) and analyzing each such successive change. 
		We then break that upper bound into two parts, one deterministic and one stochastic. For the deterministic part, the same analysis in the proof of \cref{prop:recur_d} yields convergence to $0$.
		
		The stochastic part turns out to be the tail of a martingale.
		By leveraging the property of the step-size and a crucial property of martingale differences (two martingale differences at different time steps are uncorrelated), we establish that said martingale is $L_2$-bounded.
		Then, by applying a version of Doob's martingale convergence theorem, it follows that said martingale converges almost surely to a limit random variable with finite second moment (and hence almost surely finite).
		Consequently, writing the tail as a difference between two terms (each of which converges to the same limit variablewith probability $1$), we conclude that the tail converges to $0$ \textup(a.s.\textup)\xspace.
		
		
	
	
		
		\item 
		The full \ac{DASGD} update may then be written as
		\begin{flalign}
		\debug X_{\debug n}
		&= \mirror(\debug Y_{\debug n})
		\notag\\
		\debug Y_{\debug n+1}
		&= \debug Y_{\debug n} - \debug\alpha_{\debug n+1} \bracks{\v(\debug X_\debug n) + B_{\debug n} + U_{\debug n+1}}.
		\end{flalign}
		As in Step 2 of the proof of \cref{prop:recur_d}, we again bound the one-step change of the energy function 
		$\fench(\debug Y_{\debug n+1}) - \fench(\debug Y_{\debug n})$ and then telescope the differences.
		The two distinctions from the determinstic case are: 1) Everything is now a random variable. 2) We have three terms: in addition to the random gradient $\v(\debug X_\debug n)$ and the random drift
		$B_{\debug n}$, we also have a martingale term 
		$U_{\debug n+1}$. Since $B_{\debug n}$ converges to $0$ almost surely (as shown in the previous step),
		its effect can be shown to be negligible.
		Futher, the analysis utilizes law of large numbers for martingale as well as Doob's martingale convergence theorem to bound the effect of the various martingale terms and to establish that the final dominating term converges to $-\infty$ almost surely (which generates a contradiction since the energy function is always positive) unless a subsequence $\debug X_{\debug n_k}$ converges almost surely to $\feas^{\ast}$.
		\hfill$\blacksquare$
	\end{enumerate}	 

Even though recurrence, which can be seen as the counterpart of \cref{prop:recur_d}, holds as per the above proposition, the random iterates $\debug X_{\debug n}$ are much more irregular than their deterministic counterpart $\debug x_{\debug n}$ in \ac{DAGD}. To deal with this complexity, we work with and characterize the sample trajectories ``generated" (to be made precise later) by $\debug X_{\debug n}$ (rather than individual iterates $\debug X_{\debug n}$).
To work towards this general direction, we first push the \ac{DASGD} update into a determinstic ordineary differential equation (ODE), as explained in the next subsection.

\subsubsection{Mean-Field Approximation of \ac{DASGD}}

We can rewrite the \ac{DASGD} update as:
\begin{flalign}
\label{eq:help}
\debug X_{\debug n}
&= \mirror(\debug Y_{\debug n})
\notag\\
\debug Y_{\debug n+1}
&= \debug Y_{\debug n} - \debug\alpha_{\debug n+1}\{\v(\debug X_\debug n) + B_{\debug n} + U_{\debug n+1}\}.
\end{flalign}
Written in this way, \ac{DASGD} can be viewed as a discretization of the ``mean-field'' ODE
\begin{flalign}
\label{eq:ode}
\debug x
&= \mirror(\debug y),
\notag\\
\dot\debug y
&= -\v(\debug x).
\end{flalign}
The intuition is that this ODE provides a ``mean" approximation of the \ac{DASGD} update, because in \eqref{eq:help},
the noise term $U_{\debug n+1}$ has $0$ mean, and the term $B_{\debug n}$ converges to $0$ (and therefore has
negilible effect in the long run).
Thes leaves only the term $\debug Y_{\debug n+1} = \debug Y_{\debug n} - \debug\alpha_{\debug n+1}\v(\debug X_\debug n)$, which can be seen as a Euler discretization of the ODE.
(Of course, that Equation~\eqref{eq:ode} is a good-enough approximation of Equation~\eqref{eq:help} for global almost sure convergence purposes here will be rigorously justified later.) 

Next, writing Equation~\eqref{eq:ode} solely in terms of $\debug y$ yields $\dot{\debug y}= -\v(\mirror(\debug y))$.
Since $\v$ and $\mirror$ are both Lipschitz continuous and $\mathcal{\debug X}$ is a compact set, the composition $\v\circ\mirror$ is itself Lipschitz continuous and bounded.
Standard results from the theory of dynamical systems then show that~\eqref{eq:ode} admits a unique global solution $\debug y(t)$ for any initial condition $\debug y(0)$.
On the other hand, since $\mirror$ is not a one-to-one map, it is not invertible;
consequently, there need not exist a unique solution trajectory for $\debug x(t)$.
By this token, the rest of our analysis will focus on the trajetory of $\debug y(t)$.

With the guarantee of the existence and uniqueness of the $\debug y$ trajectory, let $P\colon\mathbb{R}_{+}\times\mathbb{R}^{\debug d}\to\mathbb{R}^{\debug d}$ be the semiflow\footnote{See \cref{app:first} for a more rigorous definition. Furthermore, $\mathbb{R}_+$ is the set of non-negative reals.} of \eqref{eq:ode}, i.e., $P(t, y_0)$ denotes the state of \eqref{eq:ode} at time $t$ when the initial condition is $y_0$.
In other words, when viewed as a function of time, $P(\cdot, y_0)$ is the solution trajectory to $\dot{\debug y}= -\v(\mirror(\debug y))$.
It is worth pointing out that writing it in this double-argument form also allows us to interpret $P$ as a function of the initial condition: for a fixed $t$, $P(t, \cdot)$ gives different states at $t$ when the ODE starts from different initial conditions (in particular, $P(0, \debug y) = \debug y$). Both views will be useful later.

It turns out that with the energy function introduced here, 
$\fench(\debug P(t, \debug y))$ is always non-increasing. Furthermore, it is also decreasing at a meaingful rate.
We end this subsection with a ``sufficient decrease'' property of the mean dynamics \eqref{eq:ode} (the proof given in the appendix due to space limitation):
\begin{lemma}
	\label{lem:odeprop}\quad
	With notation as above, we have:
	\vspace{2ex}
	\begin{enumerate}
		\addtolength{\itemsep}{1ex}
		\item
		If $\mirror(\debug P(t, \debug y)) \notin \feas^{\ast}$,  then $\fench(\debug P(t, \debug y))$ is strictly decreasing in $t$ for all $y\in\mathbb{R}^{\debug d}$.
		\item
		For all $\delta > 0$, there exists some $T\equiv T(\delta) > 0$ such that, for all $t\geq T$, we have
		\begin{equation}
		\hspace{-2em}
		\sup\nolimits_{\debug y}
		\setdef{\fench(\debug P(t, \debug y)) - \fench(\debug y)}{\fench(\debug P(t,\debug y)) \ge  \delta/2}
		\le -\delta/2.
		\end{equation}
	\end{enumerate}
	\end{lemma}	

\cref{lem:odeprop} essentially says $\fench(\debug P(t, \debug y))$ is strictly and uniformly (across all $\debug y$) decreasing at a non-vanishing rate. 
More specifically, to give some intuition of the second part of \cref{lem:odeprop}, note that $\fench(\debug P(t, \debug y)) - \fench(\debug y)$
is the energy change at time $t$ when starting at $\debug y$.
The first part says this quantity is always non-negative. While the second part
says provided $\fench(\debug P(t,\debug y)) \ge  \delta/2$,
the decrease in energy will be at least $\frac{\delta}{2}$, no matter what the initial point $\debug y$ is. If $\fench(\debug P(t,\debug y)) \ge  \delta/2$ does not hold,
that means the energy at time $t$ is already really small (i.e. $\fench(\debug P(t,\debug y)) <  \delta/2$). Consequently, the mantra of the above lemma can be stated succinctly as follows: either the energy is already close to $0$, or the energy will decrease towards $0$. 

In fact, by some additional analysis, one can further show\footnote{Although this is an interesting conclusion, we do not prove it here because we are mainly concerned with establishing convergence of the \ac{DASGD} iterates, rather than the ODE solution trajectory.} that \cref{lem:odeprop} implies $  \debug P(t, \debug y) \to \feas^{\ast}, \forall y$ as $t \to \infty$. 
Now, despite being an interesting result on its own (which establishes that the continuous dynamics of \ac{DASGD} converges to $\feas^{\ast}$), it is still 
some distance away from our final desideratum: our goal is to establish (almost sure) convergence of the discrete-time process in \ac{DASGD}.
So unless we can somehow relate the discrete-time iterates to the continuous-time trajectories, the convergence of \ac{DASGD} is still uncertain.
We fulfill this taks in the next subsection.

\subsubsection{Relating \ac{DASGD} Iterates to ODE Trajectories}

Our goal here is to establish a quantitative connection between the 
\ac{DASGD} iterates and the ODE trajectory studied in the previous subsection.
Our general idea is that if we show the trajectory generated by the discrete-time iterates of \ac{DASGD} is \textbf{path-by-path} ``close" to the continuous-time trajectory $\debug P(t, \debug y)$, then likely almost-sure convergence of the \ac{DASGD} iterates can be guaranteed as well.

To be more specific, there are two things that need to be more precisely defined from the preceding high-level discussion.
First, what does it mean to be a trajectory generated by the discrete-time iterates of \ac{DASGD}?
Second, what does it mean for this trajectory to be ``close" to the ODE trajectory?
In general, the answers to these questions can vary depend on the specific goal at hand.
In the current context, our goal is to establish global almost sure convergence (a very strong result).
Consequently, we need to choose the answers rather judiciously: on the one hand, the answers must be stringent enough to ensure
global almost sure convergence in the end (for instance, for the second question, a fairly strong notion of ``closeness" is needed); on the other hand, they must also be flexible enough to fit in the current context.

As it turns out, the answer to the first question is rather intuitive: (perhaps) the simplest way to generate a continuous trajectory from a sequence of discrete points is the \textit{affine interpolation}: connect the iterates $\debug Y_0, \debug Y_1,\dots, \debug Y_{\debug n}$ at times $0, \alpha_1, \dots, \sum_{r=1}^{\debug n=1}\alpha_r$.
We call this curve the affine interpolation curve of \ac{DASGD} and denote it by $A(t)$. 
Note that $A(t)$ is a random curve because the \ac{DASGD} iterates $\debug Y_0, \debug Y_1,\dots, \debug Y_{\debug n}$ are random. To avoid confusion, we summarize the three different objects discussed so far:
\vspace{2ex}
\begin{enumerate}
	\addtolength{\itemsep}{1ex}
	\item The \ac{DASGD} iterates $\debug Y_0, \debug Y_1,\dots, \debug Y_{\debug n}$.
	\item The affine interpolation curve $A(t)$ of $\debug Y_{\debug n}$.
	\item The  flow $P(t, \debug y)$ of the ODE \eqref{eq:ode}.
\end{enumerate}

The answer to the second question lies in the notion of an \acdef{APT} , a concept introduced by \citet{BH96} and \citet{BS00}.
Specifically, in our current context, a continuous curve $s(t)$ is considered close
to ODE solution $P(t, \debug y)$ if the following holds:

\begin{definition}\quad
	A continuous function $s: \mathbb{R}_+ \to \mathbb{R}^\debug d$ is an \ac{APT} for $\debug P$ if for every $T > 0$, 
	\begin{equation}
	\lim_{t \to \infty} \sup_{0 \le h \le T} d(s(t+h), \debug P(h, s(t))) = 0,
	\end{equation} 	
	where $d(\cdot, \cdot)$ is the Euclidean metric\footnote{As should be obvious from the definition, APTs can be defined more generally in metric spaces in exactly the same way.} in $\mathbb{R}^\debug d$.
\end{definition}

Intuitively, the definition matches exactly the naming:
$s$ is an \ac{APT} for $\debug P$ if, for sufficiently large $t$, the flow lines of $P$ remain arbitrarily close to $s(t)$ over a time window of any (fixed) length.
More precisely, for each fixed $T > 0$, one can find a large enough $t_0$, such that for all $t > t_0$,
the curve $s(t+h)$ approximates the trajectory $\debug P(h, s(t))$ on the interval $h \in [0, T]$ with any predetermined degree of accuracy. 

With this definition in place, to push through the agenda, we need to establish that 
$A(t)$, the affine interpolation curve of the \ac{DASGD} iterates, is an APT for the flow $P(t, \debug y)$ induced by the ODE~\eqref{eq:ode}. 
More precisely, we establish that $A(t)$ is an APT for the flow $P(t, \debug y)$ almost surely, because as mentioned before, $A(t)$ is a random curve.


\begin{lemma}\label{lem:apt}\quad
Let	$A(t)$ be the random affine interpolation curve generated from the \ac{DASGD}
iterates.
Then $A(t)$ is an \ac{APT} of $P(t,\debug y)$ almost surely.
\end{lemma}

Note that this result means any resulting affine interpolation curve of \ac{DASGD} is close to the ODE trajectory. This also forms the basis for reasoning convergence on a path-by-path scale.
However, more work still remains to be done 
because, unfortunately, the condition that $A(t)$ is an \ac{APT} for $P$ almost surely does not guarantee that the discrete-time iterates of \ac{DASGD} converge to $\feas^{\ast}$ (for many counterexamples in general dynamical systems, see \citet{Ben99}). 
In other words, the notion of \ac{APT} is not sharp enough to ensure direct convergence result.
We fulfill this final gap in the next subsection.



\subsection{Main Convergence Result}

Even though $A(t)$ being an \ac{APT} for $P$ almost surely does not itself guarantee that the discrete-time iterates of \ac{DASGD} converge to $\feas^{\ast}$,
we can use the energy function to further sharpen this result. Specifically, we use $E(A(t))$ to further control the behavior 
of the affine interpolation curve.
In fact, the advantage of working with the affine interpolation curve $A(t)$ is that 
once we show $\fench(A(t))$ is bounded by some $\delta$ almost surely from some point onwards, then we know
$\fench(\debug Y_{\debug n})$ is bounded by $\delta$ almost surely (also from some point onwards): this is because the 
discrete-time iterates and the affine curve coincide at discrete time points.
Consequently, we focus on boudning $\fench(A(t))$, which will enable us to
 obtain our main convergence result:

\begin{theorem}
	\label{thm:main}\quad
	Under \cref{asm:basic,asm:VC,asm:delays,assump:reci}, the global state variable $\debug X_{\debug n}$ of \ac{DASGD} \textpar{\cref{alg:DASGD}} converges \textup(a.s.\textup)\xspace to the solution set $\feas^{\ast}$ of \eqref{eq:opt}.
\end{theorem}

Again, we only give an outline of the proof below.
	By \cref{prop:recur_s}, $\debug X_{\debug n}$ gets arbitrarily close to $\feas^{\ast}$ infinitely often.	
	Thus, it suffices to show that, if $\debug X_{\debug n}$ ever gets $\epsilon$-close to $\feas^{\ast}$, all the ensuing iterates are $\epsilon$-close to $\feas^{\ast}$ \textup(a.s.\textup)\xspace.
The way we show this ``trapping" property is to use the energy function. Specifically, we consider $\fench(A(t))$ and show that
	no matter how small $\epsilon$ is, for all sufficiently large $t_0$, if $\fench(A(t_0))$ is less than $\epsilon$ for some $t_0$, then $\fench(A(t)) < \epsilon, \forall t> t_0$.
	This would then complete the proof because $A(t)$ actually contains all the \ac{DASGD} iterates, and hence if 
	$\fench(A(t)) < \epsilon, \forall t> t_0$, then $\fench(\debug Y_{\debug n}) < \epsilon$ for all sufficiently large $\debug n$. Furthermore, since $A(t)$ contains all the iterates, the hypothesis  that `` if $\fench(A(t_0))$ is less than $\epsilon$ for some $t_0$" will be satisfied due to \cref{prop:recur_s}.
	

	We expand on one more layer of detail and defer the rest into appendix.
	To obtain control $\fench(A(t))$, we control two things: the energy on the ODE path $\fench(\debug P(t, y))$ and the discrepancy between $\fench(\debug P(t, y))$ and $\fench(A(t))$.
	The former can be made arbitrarily small as a result of \cref{lem:odeprop} (we have a direct handle on how the ODE path would behave). The latter can also be made arbitrarily small as a result of \cref{lem:apt}:
	since $A(t)$ is an \ac{APT} for $\debug P$, the two paths are close. Therefore, the discrepancy between $\fench(\debug P)$ and $\fench(A)$ should also be vanishingly small.
	Consequently, since $\fench(A(t)) = \fench(\debug P(t, y)) +  \{\fench(A(t))- \fench(\debug P(t, y))\}$,
	and both terms on the right can be made arbitrarily small, so can $\fench(A(t))$ be made arbitrarily small.



