\section{Introduction}\label{sec:introduction}
Operator theory has been shown to offer a valuable formalism for the analysis  of optimization algorithms, both in traditional static settings \cite{ryu_primer_2016,combettes_monotone_2018,bauschke_convex_2017,davis_convergence_2016} and in online optimization and learning \cite{shalev-shwartz_online_2011,simonetto_time-varying_2017,simonetto_time-varying_2020,Jadbabaie2015}. The key insight is that  algorithmic steps can be interpreted as operators (or maps when working in an Euclidean space), with fixed points of the operators coinciding with (or being utilized to compute) optimal solutions of the optimization problem \cite{bauschke_convex_2017,combettes_monotone_2018,davis_convergence_2016}. This link between optimization and operator theory gives access to powerful tools to characterize the convergence of existing algorithms \cite{bauschke_convex_2017,ryu_primer_2016,davis_convergence_2016}, and it further  inspires the design of new ones, \emph{e.g.} \cite{beck_fast_2009,davis_three-operator_2017,themelis_supermann_2019}.

In many applications, ranging from machine learning \cite{bottou2018optimization,liu_primer_2020} distributed optimization \cite{bastianello_asynchronous_2020}, to data-driven optimization, algorithms may feature sources of stochasticity, thus leading to \emph{stochastic} recursions. Examples of stochastic algorithms include first-order methods with inexact gradient information (including zero-order methods, stochastic gradient methods, feedback optimization, \emph{etc.}), algorithms with inexact projections, and randomized algorithms where only a subset of coordinates are updated at each iteration \cite{combettes_stochastic_2015,combettes_stochastic_2019,berinde_iterative_2007,bianchi_coordinate_2016,salzo_parallel_2021}. It is typical to handle stochastic errors in the algorithmic map or operator by either assuming that their norm vanishes asymptotically, or by multiplying them by a vanishing parameter \cite{berinde_iterative_2007,borkar_concentration_2021,combettes_stochastic_2015,combettes_stochastic_2019}. While this choice allows one to prove almost sure convergence to a fixed point, in many practical applications the additive error is not (or cannot be made) vanishing; moreover, in online optimization methods algorithmic parameters (\emph{e.g.}, the step size) are usually constant. Convergence is typically established in mean and in an almost sure sense, depending on the underlying assumptions; a high probability bound was recently developed in \cite{borkar_concentration_2021} under suitable assumptions on the additive errors and by employing a decaying parameter sequence. 

In the following, we briefly discuss two motivating examples of applications where models based on stochastic operators are necessary, although we emphasize that the use of stochastic operator frameworks is widely applicable to many optimization problems in engineering and science.

\begin{example}[Asynchronous distributed optimization]\label{ex:asynchronous-distributed}
The line of research into stochastic operators that perform updates randomly has found fruitful application in distributed optimization, where a set of agents cooperate towards the solution of a problem, but the updates may be performed asynchronously and with different precision \cite{bianchi_coordinate_2016,peng_arock_2016,alaviani_distributed_2020,bastianello_asynchronous_2020}. Consider a distributed system with $N$ agents, where each  agent has a local cost function $f_i : \mathbb{R} \to \mathbb{R}$. Assuming that agents are connected through a communication network, one may be interested in developing a distributed algorithm for solving the following optimization problem: 
\begin{equation}\label{eq:distributed-problem}
	\mathbold{x}^* = \argmin_{x_i \in \mathbb{R}} \sum_{i = 1}^N f_i(x_i) \qquad \text{s.t.} \ x_i = x_j \ \text{if $i$, $j$ connected}
\end{equation}
with the prerequisite that local functions are not shared and 
only node-to-node exchanges of information can occur.  Many algorithms have been proposed to this end, ranging from distributed gradient methods \cite{xin_general_2020,alghunaim_decentralized_2020} to the alternating direction method of multipliers (ADMM) \cite{boyd_distributed_2010,shi_linear_2014,peng_arock_2016} -- see \emph{e.g.} the surveys \cite{yang_survey_2019,notarstefano_distributed_2019}. In this setting, algorithms acquire a stochastic flavor because of  \emph{asynchrony} in the local computations  (\emph{e.g.} due to the agents performing computations at different rates),   \emph{communication failures} (\emph{e.g.} due to unreliable  
communication channels connecting the agents)~\cite{peng_arock_2016,bianchi_coordinate_2016,chang_asynchronous_2016,hannah_unbounded_2018,bastianello_asynchronous_2020}, and because they may employ an approximate first-order information (\emph{e.g.} in the case of stochastic gradient descent or zeroth order methods)~\cite{tang2020zeroth,bottou2018optimization}. \hfill\proofbox
\end{example}

\begin{example}[Online inexact gradient]
\label{ex:gradient}
Consider the online optimization problem characterized by the sequence $\mathbold{x}^{*,\ell} = \argmin_{\mathbold{x} \in \mathbb{R}^n} f^\ell(\mathbold{x})$, $\ell \in \mathbb{N}$, where $f$ is strongly convex and smooth uniformly in time. The problem can be solved using the online gradient method \cite{shalev-shwartz_online_2011}, provided that the gradient of $f$ can be computed.

\noindent However, in many applications access to the exact gradient $\nabla f^\ell$ is not available \cite{dixit_online_2019,liu_primer_2020,bottou2018optimization}, and the algorithmic steps must rely on an approximation $\hat{\nabla} f^\ell$. The online, inexact gradient method is then characterized by
$$
    \mathbold{x}^{\ell+1} = \mathbold{x}^\ell - \alpha \hat{\nabla} f^\ell(\mathbold{x}^\ell), \quad \ell \in \mathbb{N},
$$
where \emph{e.g.} $\hat{\nabla} f^\ell$ could be built using only functional evaluations \cite{duchi_optimal_2015,nesterov_random_2017,berahas_theoretical_2019}. \hfill\proofbox
\end{example}


Abstracting now from the particular application, let $\mathcal{T} \mathbold{x} = (\mathcal{T}_1 \mathbold{x}, \ldots, \mathcal{T}_n \mathbold{x})$ be a separable operator, defined over the direct sum of possibly infinite-dimensional Hilbert spaces (technical assumptions will be provided in \cref{sec:preliminaries}). In this paper, we are interested in studying the following \emph{stochastic Banach-Picard} iteration
\begin{equation}\label{eq:bp-prototype}
    x_i^{\ell+1} = \begin{cases}
        \hat{\mathcal{T}}_i \mathbold{x}^\ell & \text{w.p.} \ p_i \\
        x_i^\ell & \text{w.p.} \ 1 - p_i
    \end{cases}, \quad \ell \in \mathbb{N}
\end{equation}
where $\mathbold{x} = (x_1, \ldots, x_n)$ and $\hat{\mathcal{T}}_i$ is an inexact version of the operator $\mathcal{T}_i$. The stochasticity in~\cref{eq:bp-prototype} has two sources: \emph{(i)} each ``coordinate'' of the operator is updated with probability $p_i$ (capturing, \emph{e.g.}, asynchronous updates as in \cref{ex:asynchronous-distributed}), and \emph{(ii)} the update is performed using an inexact version of the operator (as in \cref{ex:gradient}). Additionally, we will be interested in the \emph{(iii)} \emph{online} scenario in which the operator changes (possibly, randomly) over time to model applications in online optimization and learning; that is, at time $\ell+1$ we apply $\hat{\mathcal{T}}_i^{\ell+1}$ (if an update is performed), as in \cref{ex:gradient}. The iteration~\cref{eq:bp-prototype} will be studied under the assumption that the operator is either contractive or averaged (in which case~\cref{eq:bp-prototype} can be seen as a \emph{stochastic Krasnosel'ski\u{\i}-Mann\xspace}).

As in \emph{e.g.}, \cite{combettes_stochastic_2015,combettes_stochastic_2019}, we will model inexact operators as operators subject to an additive source of errors; \emph{i.e.},
$$
    \hat{\mathcal{T}}_i \mathbold{x} := \mathcal{T}_i \mathbold{x} + e_i
$$
with $e_i$ a random vector. In particular, the framework proposed in this paper relies on a model where additive errors are random vectors with a norm that follows a  sub-Weibull distribution.  That is, the norm of $e_i$ satisfies
\begin{equation}\label{eq:tail-bound-prototype}
    \pr{\norm{e_i} \geq \epsilon} \leq 2 \exp\left( - \left( \frac{\epsilon}{\nu} \right)^{1 / \theta} \right), \quad \forall \,\, \epsilon > 0
\end{equation}
for some parameters $\theta \geq 0$, $\nu > 0$. The class of sub-Weibull r.v.s allows  us to consider distributions for the additive errors that may have heavy tails \cite{vladimirova_subweibull_2020,zhang_concentration_2020,kuchibhotla_moving_2020,wong_lasso_2020}; this class is also general and it includes the  sub-Gaussian and sub-exponential classes as sub-cases~\cite{boucheron_concentration_2013,vershynin_high-dimensional_2018}.  As explained shortly, the sub-Weibull model will allow us to provide high probability bounds for the convergence of the considered stochastic Banach-Picard method by deriving pertinent concentration results. As defined in~\cref{eq:tail-bound-prototype}, sub-Weibull random variables are defined by a tail parameter $\theta$, that defines the decay rate of the tails. Besides the convenient properties of the sub-Weibull class (\emph{e.g.} closure under scaling, sum, and product) that allow for a unified theoretical analysis, there is growing evidence that high-tailed distributions arise in machine learning applications, see \cite{vladimirova_subweibull_2020,vladimirova_understanding_2019} and the pre-prints \cite{simsekli_heavy-tailed_2019,gurbuzbalaban_heavy-tail_2021}.

\vspace{.2cm}

\noindent \textbf{Contributions.} Overall, this paper offers the following contributions:
\begin{enumerate}
    \item We study the convergence of the stochastic Banach-Picard iteration, and provide convergence results  in mean and in high probability. In particular,  bounds are offered for the distance from the fixed point (if the operator is contractive) or to the \emph{cumulative fixed point residual} (if the operator is averaged).  These bounds hold for any iteration $\ell \in \mathbb{N}$ (and not only asymptotically, as is usual in the literature \cite{combettes_stochastic_2015,combettes_stochastic_2019,bianchi_coordinate_2016}) with a given, arbitrary high probability. We also show that the high-probability bounds -- in the form $ \pr{\norm{x} \leq \epsilon(\delta)} \geq 1-\delta$ for the random variable $x$ -- scale with a factor $\sim \log(1/\delta)$, as opposed to a scaling $\sim 1/\delta$  that one would obtain via Markov's inequality.    
    
    \item The framework we propose leverages a sub-Weibull model \cite{vladimirova_subweibull_2020,zhang_concentration_2020,kuchibhotla_moving_2020,wong_lasso_2020} for the norm of the additive errors. To the best of our knowledge, this is the first time that sub-Weibull models are employed in combination with operator-theoretic tools. By using sub-Weibulls in the proposed framework, we are able to model a very broad class of random variables, fitting different practical applications. 
    
    \item As mentioned above, the convergence results proposed in this paper hold with high probability, as opposed to the almost sure convergence results of \emph{e.g.} \cite{combettes_stochastic_2015}. We discuss the difference between the two viewpoints, and also further characterize the convergence of the stochastic Banach-Picard in almost sure terms (under some additional assumptions).
    
    \item We extend the analysis to the \emph{online} stochastic Banach-Picard, which is characterized by an operator that varies over time. To the best of our knowledge, the paper presents the first convergence analysis for online algorithms with sub-Weibull distributed additive noise. 
\end{enumerate}

We also mention that, in the case of sequences of errors that are vanishing, existing results can be recovered. 

\vspace{.2cm}

\noindent \textbf{Organization.} In \cref{sec:preliminaries} we review some preliminaries in operator theory, and introduce the sub-Weibull formalism, in which we derive some lemmas. In \cref{sec:framework} we formalize the stochastic Banach-Picard framework with precise assumptions, and discuss them. In \cref{sec:static} and \cref{sec:dynamic} we present and discuss convergence results in the static and online scenarios, respectively.


\section{Preliminaries}\label{sec:preliminaries}

\subsection{Operators}
Let $\mathbb{H}_i$, $i = 1, \ldots, n$ be (possibly infinite-dimensional) Hilbert spaces with inner product $\langle \cdot, \cdot \rangle_i$, induced norm $\norm{\cdot}_i$ and identity $\mathcal{I}_i : \mathbb{H}_i \to \mathbb{H}_i$. We consider the direct sum space $\mathbb{H} = \mathbb{H}_1 \oplus \cdots \oplus \mathbb{H}_n$, whose elements are $\mathbold{x} = (x_1, \ldots, x_n)$ with $x_i \in \mathbb{H}_i$ for all $i = 1, \ldots, n$. For $\mathbold{x}, \mathbold{y} \in \mathbb{H}$, we define the inner product  as $\langle \mathbold{x}, \mathbold{y} \rangle = \sum_{i = 1}^n \langle x_i, y_i \rangle_i$, and denote the induced norm and identity of $\mathbb{H}$ as $\norm{\cdot}$ and $\mathcal{I} : \mathbb{H} \to \mathbb{H}$, respectively. We consider operators $\mathcal{T} : \mathbb{H} \to \mathbb{H}$ defined as \begin{equation}
    \mathcal{T} \mathbold{x} = (\mathcal{T}_1 \mathbold{x}, \ldots, \mathcal{T}_n \mathbold{x})
\end{equation}
where $\mathcal{T}_i : \mathbb{H} \to \mathbb{H}_i$ for all $i = 1, \ldots, n$. A remark is in order. We note that each of the operators $\mathcal{T}_i$ in the definition of $\mathcal{T}$ is evaluated at the same point $\mathbold{x}$; hence, the components of $\mathcal{T} \mathbold{x}$ are computed interdependently. A particular case is when $\mathcal{T} \mathbold{x} = (\mathcal{T}_1 x_1, \ldots, \mathcal{T}_n x_n)$; in this case, we can study the convergence  of each individual component $\{ \mathcal{T}_i \}_{i = 1}^n$ in parallel. We are, however, interested in the more general class of operators $\mathcal{T} \mathbold{x} = (\mathcal{T}_1 \mathbold{x}, \ldots, \mathcal{T}_n \mathbold{x})$ as motivated in the previous section.


A central theme of the paper is to  compute fixed points of a given operator via iterative algorithms. To this end, pertinent definitions are provided in the following.

\begin{definition}[Fixed points]
Let $\mathcal{T} : \mathbb{H} \to \mathbb{H}$ be an operator. The point $\mathbold{x}^* \in \mathbb{H}$ is a \emph{fixed point} of $\mathcal{T}$ if $\mathbold{x}^* = \mathcal{T} \mathbold{x}^*$. We denote the \emph{fixed set} of $\mathcal{T}$ as $\fix(\mathcal{T}) = \{ \mathbold{x} \in \mathbb{H} \ | \mathbold{x} = \mathcal{T} \mathbold{x} \}$.
\end{definition}

\noindent Notice that a fixed point $\mathbold{x}^* = (x_1^*, \ldots, x_n^*)$ by definition is such that $x_i^* = \mathcal{T}_i \mathbold{x}^*$ for all $i = 1, \ldots, n$. In the following, we extend well-known definitions and properties of operators to the separable framework of interest (see also, \emph{e.g.}, \cite{bauschke_convex_2017}). 

\begin{definition}[Quasi-non-expansive, quasi-contractive operator]\label{def:lipschitz-operators}
The operator  $\mathcal{T}$ is \emph{$\zeta$-quasi-Lipschitz continuous}, $\zeta > 0$, if
\begin{equation}
	\norm{\mathcal{T}_i \mathbold{x} - y_i}_i \leq \zeta \norm{x_i - y_i}_i, \quad \forall \mathbold{x} \in \mathbb{H}, \mathbold{y} \in \fix(\mathcal{T}).
\end{equation}
The operator is \emph{quasi-non-expansive} if $\zeta = 1$ and \emph{quasi-contractive} if $\zeta \in (0, 1)$.
\end{definition}

\begin{definition}[Quasi-averaged operator]\label{def:averaged-operators}
Let $\alpha \in (0, 1)$, the operator $\mathcal{T}$ is \emph{$\alpha$-quasi-averaged} if
\begin{equation}
	\norm{\mathcal{T}_i \mathbold{x} - y_i}_i^2 \leq \norm{x_i - y_i}_i^2 - \frac{1 - \alpha}{\alpha} \norm{(\mathcal{I}_i - \mathcal{T}_i) \mathbold{x}}_i^2, \quad \forall \mathbold{x} \in \mathbb{H}, \mathbold{y} \in \fix(\mathcal{T}).
\end{equation}
Alternatively, $\mathcal{T}$ is $\alpha$-quasi-averaged if there exist $\alpha \in (0, 1)$ and $\mathcal{R} : \mathbb{H} \to \mathbb{H}$ quasi-non-expansive such that $\mathcal{T} = (1 - \alpha) \mathcal{I} + \alpha \mathcal{R}$ \cite[Proposition~4.35]{bauschke_convex_2017}.
\end{definition}

\begin{remark}[Terminology]
The definitions follow the usual terminology in operator theory, according to which the modifier ``quasi'' denotes a definition enunciated w.r.t. a fixed point \cite{bauschke_convex_2017}. Of course, the ``non-quasi'' version of the definitions imply the ``quasi'' versions, but not vice-versa.

\noindent We further remark that operators with Lipschitz constant strictly smaller than one are sometimes called ``contractive'' (as will be the convention in this paper) \cite{ryu_primer_2016,giselsson_linear_2017}, and sometimes ``strictly non-expansive'' \cite{bauschke_convex_2017,combettes_stochastic_2019}.
\end{remark}

\begin{remark}[Coordinate-wise properties]
The \emph{component-wise} \cref{def:lipschitz-operators} of a quasi-contractive operator differs from the standard definition, which requires
$$
	\norm{\mathcal{T} \mathbold{x} - y_i}^2 = \sum_{i = 1}^n \norm{\mathcal{T}_i \mathbold{x} - y_i}_i^2 \leq \zeta^2 \norm{\mathbold{x} - \mathbold{y}}^2 = \zeta^2 \sum_{i = 1}^n \norm{x_i - y_i}_i^2.
$$
In \cref{def:lipschitz-operators} each component must be a contraction; on the other hand, in the standard definition above the operator $\mathcal{T}$ as a whole must be a contraction.
\end{remark}

We discuss now the existence of fixed points of quasi-non-expansive operators. The following lemma can be found in \cite[Corollary~4.24]{bauschke_convex_2017}.

\begin{lemma}[Fixed set properties]\label{lem:fixed-set-properties}
Let $\mathcal{T} : \mathbb{H} \to \mathbb{H}$ be a quasi-non-expansive operator, then $\fix(\mathcal{T})$ is closed and convex. 
\end{lemma}

\cref{lem:fixed-set-properties} does not guarantee that $\fix(\mathcal{T})$ is non-empty (indeed, it can be empty; for example, in the case of a translation). The following result from \cite[Theorem~4.29]{bauschke_convex_2017} shows that restricting the domain ensures existence of at least a fixed point.

\begin{lemma}[Browder's theorem]\label{lem:browder-theorem}
Let $\mathbb{D} \subset \mathbb{H}$ be a non-empty bounded closed convex subset, and let $\mathcal{T} : \mathbb{D} \to \mathbb{D}$ be quasi-non-expansive, then $\fix(\mathcal{T}) \neq \emptyset$ .
\end{lemma}

\begin{lemma}[Banach-Picard]\label{lem:banach-picard}
Let $\mathcal{T} : \mathbb{H} \to \mathbb{H}$ be a quasi-contractive operator, then $\fix(\mathcal{T})$ is a singleton. The unique fixed point is the limit of the sequence generated by \cite[Theorem~1.50]{bauschke_convex_2017}:
$$
	\mathbold{x}^{\ell+1} = \mathcal{T} \mathbold{x}^\ell, \quad \ell \in \mathbb{N}.
$$
\end{lemma}

If we want to compute a fixed point of the quasi-non-expansive operator $\mathcal{T}$, the following results can be used.

\begin{lemma}
Let $\mathcal{T} : \mathbb{H} \to \mathbb{H}$ be a quasi-non-expansive operator with $\fix(\mathcal{T}) \neq \emptyset$. Then for any $\alpha \in (0, 1)$ we have $\fix\left( (1-\alpha) \mathcal{I} + \alpha \mathcal{T} \right) = \fix(\mathcal{T})$.
\end{lemma}
\begin{proof}
Let $\mathbold{x}^* \in \fix(\mathcal{T})$, by definition $\mathbold{x}^* = \mathcal{T} \mathbold{x}^*$. Therefore given $\alpha \in (0, 1)$ we have
$
	(1 - \alpha) \mathbold{x}^* + \alpha \mathcal{T} \mathbold{x}^* = \mathbold{x}^* - \alpha \mathbold{x}^* + \alpha \mathbold{x}^* = \mathbold{x}^*.
$
\end{proof}

\begin{lemma}[Krasnosel'ski\u{\i}-Mann\xspace]\label{lem:krasnoselskii-mann}
Let $\mathbb{D} \subset \mathbb{H}$ be a non-empty closed convex subset, and let $\mathcal{T} : \mathbb{D} \to \mathbb{D}$ be a quasi-non-expansive operator with $\fix(\mathcal{T}) \neq \emptyset$. Let $\alpha \in (0, 1)$, then the Krasnosel'ski\u{\i}-Mann\xspace iteration
$$
	\mathbold{x}^{\ell+1} = \mathcal{T}_\alpha \mathbold{x}^\ell, \quad \ell \in \mathbb{N},
$$
where
$$
	\mathcal{T}_\alpha = (1 - \alpha) \mathcal{I} + \alpha \mathcal{T},
$$
guarantees that $\norm{(\mathcal{I} - \mathcal{T}) \mathbold{x}^\ell} \to 0$ as $\ell \to \infty$ \cite[Theorem~5.15]{bauschke_convex_2017}.
\end{lemma}

Notice that the Krasnosel'ski\u{\i}-Mann\xspace iteration corresponds to the Banach-Picard iteration of \cref{lem:banach-picard} applied to the $\alpha$-averaged operator $(1 - \alpha) \mathcal{I} + \alpha \mathcal{T}$. In the remainder of the paper we will thus focus on studying the convergence of the Banach-Picard in two cases: i) when the operator is quasi-contractive, and ii) 
 when $\mathcal{T}$ is quasi-averaged.


\subsection{Probability}\label{subsec:background-probability}
In this section, we provide some definitions and results in probability theory that will be used in this paper to derive  convergence results in high-probability. Throughout the paper, the underlying probability space will be $(\Omega, \mathcal{F}, \mathbb{P})$.

In the following, we introduce the definition of sub-Weibull random variables, which are characterized in terms of their tail probabilities. We refer to \cite{vladimirova_subweibull_2020} and the pre-print \cite{kuchibhotla_moving_2020} for a discussion (with the latter presenting sub-Weibull r.v.s in a slightly different manner) and to the pre-print \cite{zhang_concentration_2020}.

\begin{definition}[Sub-Weibull random variable]\label{def:sub-weibull}
A random variable $x$ is said to be sub-Weibull if  $\exists \,\, \theta \geq 0, \nu > 0$ such that
\begin{equation}\label{eq:tail-probability}
	\pr{|x| \geq \epsilon} \leq 2 \exp\left( - \left( \frac{\epsilon}{\nu} \right)^{1 / \theta} \right), \quad \forall \,\, \epsilon > 0.
\end{equation}
\end{definition}

The tails of a sub-Weibull r.v. become heavier as the parameter $\theta$ grows larger. Moreover, setting $\theta = 1/2$ and $\theta = 1$ yields the class of sub-Gaussians and sub-exponential random variables, respectively; see, \emph{e.g.}, \cite{boucheron_concentration_2013,vershynin_high-dimensional_2018}.



A sub-Weibull r.v. can be equivalently characterized using the following properties \cite[Theorem~1]{vladimirova_subweibull_2020}. We recall that $\norm{x}_k = \ev{|x|^k}^{1 / k}$, for any $k \geq 1$.

\begin{proposition}[Equivalent sub-Weibull r.v. definition]\label{pr:equivalent-sub-weibull}
Given $\theta \geq 0$, the following properties are equivalent:
\begin{enumerate}
	\item[(i)] $\exists \,\, \nu_1 > 0$ s.t. $\pr{|x| \geq \epsilon} \leq 2 \exp\left( - \left( \epsilon / \nu_1 \right)^{1 / \theta} \right)$, $\forall \epsilon > 0$;
	
	\item[(ii)] $\exists \,\, \nu_2 > 0$ s.t. $\norm{x}_k \leq \nu_2 k^\theta$, $\forall k \geq 1$;
	
	\item[(iii)] $\exists \,\, \nu_3 > 0$ s.t. $\ev{\exp\left( (\lambda |x|)^{1 / \theta} \right))} \leq \exp\left( (\lambda \nu_3)^{1 / \theta} \right)$, $\forall \lambda \in (0, 1/\nu_3]$;
	
	\item[(iv)] $\exists \,\, \nu_4 > 0$ s.t. $\evexp{\left( |x| / \nu_4 \right)^{1 / \theta}} \leq 2$.
\end{enumerate}

\noindent The parameters $\nu_1, \nu_2, \nu_3, \nu_4$ differ each by a constant that only depends on $\theta$.
\end{proposition}

\cref{pr:equivalent-sub-weibull} indicates that we can define the class of sub-Weibull random variables using any of the four equivalent properties. In the following we leverage the property \emph{(ii)} to characterize sub-Weibull random variables, and we write $x \sim \mathrm{subW}(\theta, \nu)$ if $\norm{x}_k \leq \nu k^\theta$, $\forall k \geq 1$. This choice allows us to readily characterize the properties of the class of sub-Weibull random variables (\emph{e.g.} that it is closed to sums), while using the other properties may lead to a less clean characterization.

As mentioned in \cref{pr:equivalent-sub-weibull}, the $\nu_i$, $i = 1, 2, 3, 4$, parameters differ by a constant; we remark in particular that if property \emph{(ii)} holds with parameter $\nu_2$, then property \emph{(i)} holds with $\nu_1 = \left( 2 e / \theta \right)^\theta \nu_2$ \cite[Lemma~5]{wong_lasso_2020}, and in the following we denote $c(\theta) := \left( 2 e / \theta \right)^\theta$. The following lemma shows how the tail probability equation~\cref{eq:tail-probability} can be used to give high probability bounds.

\begin{lemma}[High probability bound]\label{lem:high-probability-bound}
Let $x \sim \mathrm{subW}(\theta, \nu)$ according to \cref{pr:equivalent-sub-weibull}~(ii), then for any $\delta \in (0, 1)$, w.p. $1 - \delta$ we have the bound:
$$
    |x| \leq \nu \log^\theta(2 / \delta) c(\theta).
$$
\end{lemma}
\begin{proof}
By \cref{pr:equivalent-sub-weibull}~\emph{(ii)} and its equivalence with definition~\emph{(i)} we have, for any $\epsilon > 0$:
$$
    \pr{|x| \geq \epsilon} \leq 2 \exp\left( - \left( \frac{\epsilon}{c(\theta) \nu} \right)^{1 / \theta} \right).
$$
Setting the RHS equal to $\delta$ and solving for $\epsilon$ we get $\epsilon = c(\theta) \nu \log^\theta(2 / \delta)$ which implies that, w.p. $1 - \delta$ we have $|x| \leq \epsilon$. The thesis follows using the expression for $\epsilon$.
\end{proof}


With this premise, we can now characterize the properties of the class of sub-Weibull random variables.

\begin{proposition}[Inclusion]\label{pr:inclusion-sub-weibull}
Let $x \sim \mathrm{subW}(\theta, \nu)$ and let $\theta'$ and $\nu'$ such that $\theta' \geq \theta$, $\nu' \geq \nu$. Then, $x \sim \mathrm{subW}(\theta', \nu')$.
\end{proposition}
\begin{proof}
By assumption we have $\norm{x}_k \leq \nu k^\theta$. Using the fact that $\nu k^\theta \leq \nu' k^{\theta'}$ (which holds since $k \geq 1$) yields the thesis; cf. \cite[Proposition~1]{vladimirova_subweibull_2020}.
\end{proof}

\begin{proposition}[Closure of sub-Weibull class]\label{pr:closure-sub-weibull-class}
The class of sub-Weibull random variables is closed w.r.t. product by a scalar, sum, product, and exponentiation, according to the following rules.
\begin{enumerate}	
	\item \emph{Product by scalar:} let $x \sim \mathrm{subW}(\theta, \nu)$ and $a \in \mathbb{R}$, then $a x \sim \mathrm{subW}(\theta, |a| \nu)$;
	
	\item \emph{Sum:} let $x_i \sim \mathrm{subW}(\theta_i, \nu_i)$, $i = 1, 2$, possibly dependent, then $x_1 + x_2 \sim \mathrm{subW}(\max\{ \theta_1, \theta_2 \}, \nu_1 + \nu_2)$;
	
	\item \emph{Product (independence):} let $x_i \sim \mathrm{subW}(\theta_i, \nu_i)$, $i = 1, 2$, independent, then $x_1 x_2 \sim \mathrm{subW}(\theta_1 + \theta_2, \nu_1 \nu_2)$.
	
	\item \emph{Product (dependence):} let $x_i \sim \mathrm{subW}(\theta_i, \nu_i)$, $i = 1, 2$, possibly dependent, then $x_1 x_2 \sim \mathrm{subW}(\theta_1 + \theta_2, c(\theta_1, \theta_2) \nu_1 \nu_2)$, $c(\theta_1, \theta_2) = (\theta_1 + \theta_2)^{\theta_1 + \theta_2} / (\theta_1^{\theta_1} \theta_2^{\theta_2})$.
	
	\item \emph{Power:} let $x \sim \mathrm{subW}(\theta, \nu)$ and $a > 0$, then $x^a \sim \mathrm{subW}(a \theta, \nu^a \max\{ 1, a^{a \theta} \})$.
\end{enumerate}
\end{proposition}
\begin{proof} The proofs of 1-5 are presented in the following. 

\emph{Proof of 1.} The result follows by $\norm{a x}_k = |a| \norm{x}_k \leq |a| \nu k^\theta$.

\emph{Proof of 2.} For completeness we report the proof provided in \cite[Proposition~3]{vladimirova_subweibull_2020}. Using the triangle inequality we write
$$
	\norm{x_1 + x_2}_k \leq \norm{x_1}_k + \norm{x_2}_k \overset{\text{(i)}}{\leq} \nu_1 k^{\theta_1} + \nu_2 k^{\theta_2} \overset{\text{(ii)}}{\leq} (\nu_1 + \nu_2) k^{\max\{ \theta_1, \theta_2 \}},
$$
where (i) holds by the assumption that $x_i$ are sub-Weibull, and (ii) holds since $k \geq 1$.

\emph{Proof of 3.} By definition of $\norm{\cdot}_k$ we can write
$$
	\norm{x y}_k \leq \ev{|x|^k |y|^k}^{1/k} \overset{\text{(i)}}{=} \ev{|x|^k}^{1/k} \ev{|y|^k}^{1/k} \overset{\text{(ii)}}{\leq} \nu_1 k^{\theta_1} \nu_2 k^{\theta_2} = \nu_1 \nu_2 k^{\theta_1 + \theta_2}
$$
where (i) holds by independence and (ii) by sub-Weibull assumption.

\emph{Proof of 4.} By definition of $\norm{\cdot}_k$ we can write
	$$
		\ev{|x y|^k} \leq \ev{|x|^k |y|^k} \overset{\text{(i)}}{\leq} \ev{|x|^{k p}}^{1/p} \ev{|y|^{k q}}^{1/q}
	$$
	where (i) used H\"older's inequality with $p, q > 1$, $1/p + 1/q = 1$. Taking the $k$-th root we get
	$$
		\norm{x y}_k \leq \norm{x}_{p k} \norm{y}_{q k} \overset{\text{(ii)}}{\leq} \nu_1 \nu_2 k^{\theta_1 + \theta_2} p^{\theta_1} q^{\theta_2}
	$$
	where (ii) holds by sub-Weibull assumption. Finally, the thesis follows minimizing $p^{\theta_1} q^{\theta_2}$ subject to $q = p / (p - 1)$.


\emph{Proof of 5.} By definition of $\norm{\cdot}_k$ we have $\norm{x^a}_k = \ev{|x^a|^k}^{1/k} = \ev{|x|^{a k}}^{1/k}$. Now, we distinguish two cases: if $0 < a < 1$ then by Jensen's inequality we have
$$
    \ev{|x|^{a k}}^{1/k} \leq \left( \ev{|x|^k}^{1/k} \right)^a \leq (\nu k^\theta)^a;
$$
instead if $a \geq 1$ then we can write
$$
    \ev{|x|^{a k}}^{1/k} = \left( \ev{|x|^{a k}}^{1/a k} \right)^a \overset{\text{(i)}}{\leq} \left( \nu (a k)^\theta \right)^a = \nu^a a^{a \theta} k^\theta
$$
where (i) holds by the fact that $a k \geq 1$ and that $x$ is sub-Weibull.
\end{proof}

\begin{remark}[Square of sub-Weibull]\label{rem:square-sub-weibull}
A consequence of \cref{pr:closure-sub-weibull-class} (by property 5., or equivalently by 4.) is that the square of $x \sim \mathrm{subW}(\theta, \nu)$ is itself a sub-Weibull, characterized by $x^2 \sim \mathrm{subW}(2 \theta, 4^\theta \nu^2)$. A particular case is the well known fact that the square of a sub-Gaussian ($\theta = 1/2$) is sub-exponential ($\theta = 1$), see \emph{e.g.} \cite[Lemma~2.7.6]{vershynin_high-dimensional_2018}.
\end{remark}


Notice that the definition of sub-Weibulls and their properties does not require that their mean be zero. Moreover, we can easily characterize centered r.v.s in the sub-Weibull framework with the following lemma, cf. \cite[Lemma~2.6.8]{vershynin_high-dimensional_2018}.

\begin{lemma}[Centered random variables]\label{lem:centered-rv}
Let $x \sim \mathrm{subW}(\theta, \nu)$; then its centered version is sub-Weibull with $x - \ev{x} \sim \mathrm{subW}(\theta, 2 \nu)$.
\end{lemma}
\begin{proof}
By triangle inequality we have $\norm{x - \ev{x}}_k \leq \norm{x}_k + \norm{\ev{x}}_k$ for any $k \geq 1$, and we want to bound $\norm{\ev{x}}_k$. We have
$$
	\norm{\ev{x}}_k = \ev{|\ev{x}|^k}^{1/k} \overset{\text{(i)}}{\leq} \ev{\ev{|x|^k}}^{1/k} = \ev{|x|^k}^{1/k} = \norm{x}_k
$$
where (i) holds by Jensen's inequality. Therefore we have $\norm{x - \ev{x}}_k \leq 2 \norm{x}_k \leq 2 \nu k^\theta$ which implies the thesis.
\end{proof}


We can see that bounded r.v.s are sub-Weibull with $\theta = 0$; indeed, let $x$ be a r.v. such that a.s. $|x| \leq b$, then $\norm{x}_k \leq b = b k^0$, $k \geq 1$. This characterization is ``optimal'' in terms of $\theta$, since $\theta = 0$ corresponds to the lightest possible tail. However, it is sub-optimal in the other parameter, $\nu$, which does not reflect the overall distribution of $x$, only its maximum absolute value. The following Lemma provides a sub-Gaussian ($\theta = 1/2$) characterization of bounded r.v.s based on Hoeffding's inequality (see \emph{e.g.} \cite[Theorem~2.2.6]{vershynin_high-dimensional_2018}).

\begin{lemma}[Bounded random variables]\label{lem:bounded-rv}
Let $x$ be a random variable with mean $\mu := \ev{x}$, such that $x \in [a, b]$ a.s.. Then $x - \mu \sim \mathrm{subW}(1/2, (b - a) / \sqrt{2})$.
\end{lemma}
\begin{proof}
The proof is adapted from \cite[Proposition~D.1]{kuchibhotla_moving_2020}, which uses Bernstein's inequality rather than Hoeffding's. By Hoeffding's inequality \cite[Theorem~2.2.6]{vershynin_high-dimensional_2018} we have
$$
	\pr{|x - \mu| \geq \epsilon} \leq 2 \exp\left( - \left( \frac{\epsilon}{\nu} \right)^2 \right)
$$
with $\nu = (b - a) / \sqrt{2}$. Considering then $k \geq 2$, we can write
\begin{align*}
	\ev{|x - \mu|^k} &= \int_0^\infty k \epsilon^{k - 1} \pr{|x - \mu| \geq \epsilon} d\epsilon 
	\overset{\text{(i)}}{\leq} 2 \int_0^\infty k \epsilon^{k-1} \exp\left( - \left( \frac{\epsilon}{\nu} \right)^2 \right) d\epsilon \\
	&\overset{\text{(ii)}}{=} 2 \left( \frac{\nu}{\sqrt{1}} \right)^k \int_0^\infty k z^{k-1} \exp(- z^2 / 2) dz
	\overset{\text{(iii)}}{=} 2 \nu^k \Gamma(1 + k/2)
\end{align*}
where (i) holds by Hoeffding's inequality, (ii) by the change of variable $z = \sqrt{2} \epsilon / \nu$, and (iii) by \cite[Ex.~3.3.4(a)]{gine_mathematical_2016}. Consequently, we have $\norm{x - \mu}_k \leq 2^{1/k} \nu \Gamma(1 + k/2)^{1/k}$. Finally, using Stirling's inequality (since $k \geq 2$) we have the bound
$$
	2^{1/k} \Gamma(1 + k/2)^{1/k} \leq 2^{1/k} \left( \frac{e}{2} \right)^{1/k} k^{1/k} \left( \frac{k}{2 e} \right)^{1/2} \overset{\text{(iv)}}{\leq} k^{1/2}
$$
where (iv) follows by the fact that $(e k)^{1/k} \leq (2 e)^{1/2}$, and the thesis is proved.
\end{proof}


\section{Framework}\label{sec:framework}
In this section, we formalize the framework of the stochastic Banach-Picard iteration to compute the fixed point(s) of inexact operators. We also present the online version of the stochastic Banach-Picard, characterized by a time-varying sequence of operators (and, hence, a time-varying sequence of sets of fixed points).

\subsection{Stochastic Banach-Picard}

The focus of this paper is the study of the Banach-Picard iteration within a stochastic framework, and in \cref{sec:introduction} we provided some motivations. Mathematically, we are interested in analyzing the following update:
\begin{align}
	x_i^{\ell+1} &= \begin{cases}
		\hat{\mathcal{T}}_i \mathbold{x}^\ell = \mathcal{T}_i \mathbold{x}^\ell + e_i^\ell & \text{w.p.} \ p_i \\
		x_i^\ell & \text{w.p.} \ 1 - p_i
	\end{cases} , \hspace{.5cm} i = 1, \ldots, n \nonumber \\
	&= (1 - u_i^\ell) x_i^\ell + u_i^\ell \mathcal{T}_i \mathbold{x}^\ell + u_i^\ell e_i^\ell , \hspace{1.4cm} i = 1, \ldots, n \label{eq:stochastic-update}
\end{align}
where $\ell \in \mathbb{N}$ is the iteration index, $u_i^\ell$ are Bernoulli random variables that indicate whether a coordinate is updated or not at iteration $\ell$, and $\mathbold{e}^\ell = (e_1^\ell, \ldots, e_n^\ell)$ is a random vector of additive noise. Throughout the paper, the initial condition $\mathbold{x}^0 \in \mathbb{H}$ will be treated as a deterministic quantity. The following assumptions are used to formalize this stochastic framework.

\begin{assumption}[Stochastic framework]\label{as:stochastic-framework} The following is assumed regarding the random vectors $\mathbold{u}^\ell$ and the points $\mathbold{e}^\ell$. 

\begin{itemize}
    \item[(i)] The vector $\mathbold{u}^\ell = (u_1^\ell, \ldots, u_n^\ell)$ is a realization of the $\{ 0, 1 \}^n$-valued random vector $\mathbold{u}$ such that $p_i := \pr{u_i = 1} > 0$ for any $i = 1, \ldots, n$. The vectors $\{ \mathbold{u}^\ell \}_{\ell \in \mathbb{N}}$ are i.i.d..
    
    \item[(ii)] The components of the additive error $\mathbold{e}^\ell \in \mathbb{H}$ are realizations of $\mathbb{H}_i$-valued random variables such that their norm $\norm{e_i^\ell}_i$ is sub-Weibull with parameters $\theta$ and $\nu_i^\ell$, $i = 1, \ldots, n$, and with mean $\mu_i^\ell := \ev{\norm{e_i^\ell}_i} < +\infty$.
    
    \item[(iii)] The random processes $\{ \mathbold{u}^\ell \}_{\ell \in \mathbb{N}}$ and $\{ \mathbold{e}^\ell \}_{\ell \in \mathbb{N}}$ are independent of each other.
\end{itemize}
\end{assumption}

Remarks regarding these assumptions on the additive errors are provided below. In the next sections we study the stochastic Banach-Picard under one of the following assumptions.

\begin{assumption}[Quasi-contractive]\label{as:contractive}
The operator $\mathcal{T} : \mathbb{H} \to \mathbb{H}$ in update~\cref{eq:stochastic-update} is $\zeta$-quasi-contractive.
\end{assumption}

\begin{assumption}[Quasi-averaged]\label{as:averaged}
Let $\mathbb{D} = \mathbb{D}_1 \oplus \cdots \oplus \mathbb{D}_n$ be a non-empty bounded closed convex subset of $\mathbb{H}$, the operator $\mathcal{T} : \mathbb{D} \to \mathbb{D}$ in update~\cref{eq:stochastic-update} is $\alpha$-quasi-averaged.
\end{assumption}

By \cref{as:contractive} the operator has a unique fixed point; this case will be considered in~\cref{pr:contractive}. On the other hand, under \cref{as:averaged} Browder's theorem guarantees that $\fix(\mathcal{T}) \neq \emptyset$; this case will be considered in~\cref{pr:averaged}. 
Hereafter we will denote by $\diam(\mathbb{D}) < \infty$ the diameter of $\mathbb{D}$.

\smallskip

In the following sections we discuss the assumptions that characterize the stochastic framework of~\cref{eq:stochastic-update}.


\subsubsection{Update model}
We note that \cref{as:stochastic-framework}(i) does not required independence among the components of $\mathbold{u}^\ell$. Independence is only assumed between $\mathbold{u}^\ell$ and $\mathbold{u}^h$, for any pair $\ell, h \in \mathbb{N}$, $\ell \neq h$. This is the case, for example, for distributed optimization algorithms described in \cref{ex:asynchronous-distributed}, where $\mathbold{u}^\ell$ models communication drops; these drops may be independent across time. 


We also note that under \cref{as:stochastic-framework}(i) there may exist times during which \emph{none} of the coordinates are updated, because $u_i^\ell = 0$ for all $i = 1, \ldots, n$. This is different from the framework used in \emph{e.g.} \cite{combettes_stochastic_2015,combettes_stochastic_2019}, where at least one $u_i^\ell$ must always be different from zero. One may argue that in \cite{combettes_stochastic_2015} iterations that do not see an update are not counted, thus the event $u_i^\ell = 0$ for all $i = 1, \ldots, n$ can be excluded. However, even an iteration without updates requires some time (in terms of computational time elapsed), and choosing not to count them would skew the convergence rate evaluation; this considerations will be important when we will study the convergence of online algorithms later in this paper. 


\subsubsection{Additive error model}
The additive error model defined in \cref{as:stochastic-framework}(ii) is very general, since it imposes only a tail condition on the error norm. For example, by the discussion leading to \cref{lem:bounded-rv} we know that if $\norm{e_i^\ell}_i < b < + \infty$ a.s., then $\norm{e_i^\ell}_i$ is sub-Weibull. Furthermore, in finite dimensional spaces we have the following result, which can be seen as an extension of \cite[Theorem~3.1.1]{vershynin_high-dimensional_2018} to the sub-Weibull case.

\begin{lemma}[Norm of sub-Weibull vectors]
Let $\mathbold{e} = [e_1, \ldots, e_d]^\top$ be a random vector in $\mathbb{R}^d$, $d < +\infty$, such that $e_i \sim \mathrm{subW}(\theta, \nu)$, $i = 1, \ldots, d$. Then the Euclidean norm of $\mathbold{e}$ is sub-Weibull with
$$
    \norm{\mathbold{e}} \sim \mathrm{subW}(\theta, 2^\theta \sqrt{d} \nu).
$$
\end{lemma}
\begin{proof}
We want to characterize $\norm{\mathbold{e}} = \left( \sum_{i = 1}^d e_i^2 \right)^{1/2}$ as a sub-Weibull. By \cref{pr:closure-sub-weibull-class} we know that $e_i^2 \sim \mathrm{subW}(2 \theta, 2^{2\theta} \nu^2)$, and it follows that $\sum_{i = 1}^d e_i^2 \sim \mathrm{subW}(2 \theta, d 2^{2\theta} \nu^2)$. Finally, taking the square root implies $\norm{\mathbold{e}} \sim \mathrm{subW}(\theta, \sqrt{d} 2^\theta \nu)$.
\end{proof}

We also remark that for simplicity, we take the error norms to have equally heavy tails, that is, they have the same parameter $\theta$. This choice can be relaxed, allowing $\norm{e_i^\ell}_i \sim \mathrm{subW}(\theta_i, \nu_i^\ell)$, and the convergence results of the following section hold by using $\theta = \max_i \theta_i$.

Moreover, the assumption does not require independence of the error components at a fixed time $\ell$, or between the errors drawn at different times $\ell$ and $h$. Finally, we remark that the errors could be \emph{biased}, that is, they could have mean different from zero (either positive or negative). The sub-Weibull model that we use indeed allows for biased r.v.s, see \emph{e.g.} the discussion in \cite{vladimirova_subweibull_2020} and \cref{lem:centered-rv}.


\subsection{Extension to online Banach-Picard}\label{subsec:dynamic-banach-picard-framework}
The stochastic Banach-Picard in \cref{eq:stochastic-update} is designed to (approximately) compute one fixed point of the time-invariant operator $\mathcal{T}$. However, in many applications we may be interested in computing the fixed points of a \emph{sequence of operators} \cite{simonetto_time-varying_2017}, \emph{e.g.} in online (or time-varying) optimization. In this section we formalize the online version of~\cref{eq:stochastic-update}, whose convergence will be analyzed in \cref{sec:dynamic}.

Let $\{ \mathcal{T}^\ell \}_{\ell \in \mathbb{N}}$ be a sequence of quasi-non-expansive operators, we define the \emph{online Banach-Picard} (also referred to as ``dynamic'' or ``running'') as
\begin{equation}\label{eq:dynamic-stochastic-update}
	x_i^{\ell+1} = (1 - u_i^\ell) x_i^\ell + u_i^\ell \mathcal{T}_i^{\ell+1} \mathbold{x}^\ell + u_i^\ell e_i^\ell
\end{equation}
which differs from the static update since the (inexact) operator may change at each step. As the convergence results will highlight, the dynamic nature of the operator introduces an additional source of inexactness. This is due to the fact that a fixed point of $\mathcal{T}^\ell$ in general may not approximate well a fixed point of the subsequent operator $\mathcal{T}^{\ell+1}$.

As in the previous section, we distinguish two cases: each operator in the sequence $\{ \mathcal{T}^\ell \}_{\ell \in \mathbb{N}}$ is (i) quasi-contractive or (ii) quasi-averaged. In case (i), the operators have each a unique fixed point, denoted by $\mathbold{x}^{*,\ell}$, and the analysis of~\cref{eq:dynamic-stochastic-update} is in terms of tracking the fixed point sequence $\{ \mathbold{x}^{*,\ell} \}_{\ell \in \mathbb{N}}$. In case (ii) each operator has a set -- non-empty, closed and convex -- of fixed points, $\fix(\mathcal{T}^\ell)$, and the Banach-Picard is tracking any one of the possible sequences $\{ \mathbold{x}^{*,\ell} \in \fix(\mathcal{T}^\ell) \}_{\ell \in \mathbb{N}}$.

We can now characterize the time variability of the operator in terms of the distance between consecutive fixed points. The following assumption formalizes this for the case (i).

\begin{assumption}[Online quasi-contractive]\label{as:dynamic-contractive}
Let $\{ \mathcal{T}^\ell \}_{\ell \in \mathbb{N}}$ be the sequence of quasi-contractive operators in~\cref{eq:dynamic-stochastic-update}, we define the distance between (the $i$-th components of) consecutive fixed points as
$$
	\sigma_i^\ell := \norm{x_i^{*,\ell+1} - x_i^{*,\ell}}_i \, .
$$
Further, we assume that $\sigma_i^\ell \sim \mathrm{subW}(\varphi, \gamma_i^\ell)$.
\end{assumption}

Before stating the assumption on the variability of the fixed-points for the quasi-on-expansive operator in case (ii), we introduce the following definitions.

\begin{definition}[Set-to-set distance]\label{def:set-set-distance}
Let $\mathbb{C}, \mathbb{D} \subset \mathbb{H}$ be non-empty convex sets. We define the \emph{minimal distance} of sets $\mathbb{C}, \mathbb{D}$ as
$$
	\ldist(\mathbb{C}, \mathbb{D}) = \inf_{\mathbold{x} \in \mathbb{C}, \mathbold{y} \in \mathbb{D}} \norm{\mathbold{x} - \mathbold{y}}.
$$

\noindent We define the \emph{maximal distance} -- commonly called \emph{Pompeiu-Hausdorff distance}, see \emph{e.g.} \cite[Section~3.1]{doncev_implicit_2014} -- of sets $\mathbb{C}, \mathbb{D}$ as
$$
	\udist(\mathbb{C}, \mathbb{D}) = \max\left\{ \sup_{\mathbold{x} \in \mathbb{C}} \inf_{\mathbold{y} \in \mathbb{D}} \norm{\mathbold{x} - \mathbold{y}}, \sup_{\mathbold{y} \in \mathbb{D}} \inf_{\mathbold{x} \in \mathbb{C}} \norm{\mathbold{x} - \mathbold{y}} \right\}.
$$
\end{definition}

\begin{assumption}[Online quasi-averaged]\label{as:dynamic-averaged}
Let $\{ \mathcal{T}^\ell \}_{\ell \in \mathbb{N}}$ be the sequence of quasi-averaged operators in~\cref{eq:dynamic-stochastic-update}, we define the minimal and maximal distance between (the $i$-th components of) consecutive fixed point sets as
$$
	0 \leq \underline{\sigma}_i^{\ell} := \ldist\left( \fix(\mathcal{T}_i^{\ell+1}), \fix(\mathcal{T}_i^\ell) \right) \ \leq \ \udist\left( \fix(\mathcal{T}_i^{\ell+1}), \fix(\mathcal{T}_i^\ell) \right) =: \overline{\sigma}_i^\ell \, .
$$
We assume $\underline{\sigma}_i^\ell \sim \mathrm{subW}(\underline{\varphi}, \underline{\gamma}_i^\ell)$ and $\overline{\sigma}_i^\ell \sim \mathrm{subW}(\overline{\varphi}, \overline{\gamma}_i^\ell)$.
\end{assumption}


\begin{remark}
We remark that minimal and maximal distances in \cref{as:dynamic-averaged} are finite, since the operators have the bounded $\mathbb{D}$ as domain and range. By definition we also can see that, for any pair $x_i^{*,\ell+1} \in \fix(\mathcal{T}_i^{\ell+1})$ and $x_i^{*,\ell} \in \fix(\mathcal{T}_i^\ell)$, it holds
$$
	\underline{\sigma}_i^\ell \leq \norm{x_i^{*,\ell+1} - x_i^{*,\ell}} \leq \overline{\sigma}_i^\ell.
$$
Moreover, $\underline{\sigma}_i^\ell = 0$ if $\fix(\mathcal{T}_i^{\ell+1}) \cap \fix(\mathcal{T}_i^\ell) \neq \emptyset$, since there is at least a common fixed point; but, unless $\fix(\mathcal{T}_i^{\ell+1}) \equiv \fix(\mathcal{T}_i^\ell)$, then $\overline{\sigma}_i^\ell > 0$.
\end{remark}


\begin{remark}
In the online Banach-Picard~\cref{eq:dynamic-stochastic-update}, the $(\ell + 1)$-th iterate $x_i^{\ell+1}$ is computed applying the $(\ell + 1)$-th operator. This implies that first we observe the operator $\mathcal{T}_i^{\ell+1}$, \emph{e.g.} obtaining data from a streaming source \cite{dallanese_optimization_2020,simonetto_time-varying_2020}, and then apply the Banach-Picard. An alternative scheme, closer in spirit to online learning frameworks, would be to compute $x_i^{\ell+1}$ using the $\ell$-th operator; that is, using $\mathcal{T}_i^\ell$ as a \emph{prediction} of $\mathcal{T}_i^{\ell+1}$, before the latter is made available. After $x_i^{\ell+1}$ is computed, the following operator $\mathcal{T}_i^{\ell+1}$ is observed, and we incur a \emph{regret} due to the discrepancy of the fixed points of $\mathcal{T}_i^{\ell}$ and $\mathcal{T}_i^{\ell+1}$; see \emph{e.g.} the discussion in \cite[section~4.1]{bastianello_primal_2020}. The convergence results presented in \cref{sec:dynamic} can be adapted with minimal changes to this alternative framework.
\end{remark}


\section{Convergence Analysis: Static Operators}\label{sec:static}

\subsection{Main results}

The following proposition establishes convergence in mean and high-probability of the stochastic Banach-Picard iteration~\cref{eq:stochastic-update}, when the operator $\mathcal{T}$ is quasi-contractive (and we recall that $\ell \in \mathbb{N}$ is the iteration index, and $c(\theta')$ is defined in \cref{subsec:background-probability}). The convergence results are then discussed in \cref{subsec:discussion-convergence}.

\begin{proposition}[Static and quasi-contractive operator]\label{pr:contractive}
Let \cref{as:stochastic-framework,as:contractive} hold, and let $\{ \mathbold{x}^\ell \}_{\ell \in \mathbb{N}}$ be the sequence generated by~\cref{eq:stochastic-update}. The following holds. 

\noindent  (i) For each coordinate, $i = 1, \ldots, n$, the mean of $\norm{x_i^\ell - x_i^*}_i$ can be upper bounded by:
$$
	\ev{\norm{x_i^\ell - x_i^*}_i} \leq \chi_i^\ell \norm{x_i^0 - x_i^*}_i + p_i \frac{1 - \chi_i^\ell}{1 - \chi_i} \sup_\ell \mu_i^\ell
$$
where $\chi_i := 1 - p_i + p_i \zeta \geq \zeta$.

\noindent (ii) With probability $1 - \delta$, $\delta \in (0, 1)$, we have that:
$$
	\norm{x_i^\ell - x_i^*}_i \leq \log^{\theta'}(2/\delta) \ c(\theta') \left( \eta_i(\ell) \norm{x_i^0 - x_i^*}_i + \frac{1 - \zeta^\ell}{1 - \zeta} \sup_\ell \nu_i^\ell \right)
$$
where $\theta' = \max\{ 1/2, \theta \}$, and $\eta_i(\ell)$ is a monotonically decreasing function of $\ell$, see \cref{ssubsec:zeta-beta-subweibull}.
\end{proposition}
\begin{proof}
By the triangle inequality, we have:
\begin{align}
	\norm{x_i^{\ell+1} - x_i^*}_i &\leq \norm{(1 - u_i^\ell) x_i^\ell + u_i^\ell \mathcal{T}_i \mathbold{x}^\ell - x_i^*}_i + u_i^\ell \norm{e_i^\ell}_i \nonumber \\
	&\leq \zeta^{u_i^\ell} \norm{x_i^\ell - x_i^*}_i + u_i^\ell \norm{e_i^\ell}_i \label{eq:component-wise}
\end{align}
where we used the quasi-contractiveness of $\mathcal{T}$. Starting from inequality~\cref{eq:component-wise}, we now bound the mean of the error and provide a high probability bound.

\paragraph{Part (i)}
Taking the expectation of~\cref{eq:component-wise} yields
\begin{equation}\label{eq:component-wise-mean}
	\ev{\norm{x_i^{\ell+1} - x_i^*}_i} \leq \chi_i \ev{\norm{x_i^\ell - x_i^*}_i} + p_i \mu_i^\ell
\end{equation}
where we defined $\chi_i := 1 - p_i + p_i \zeta$, and notice that $\chi_i \geq \zeta$. Recursively applying~\cref{eq:component-wise-mean}, using the fact that $\mu_i^\ell \leq \sup_\ell \mu_i^\ell$, and using the geometric sum yields the thesis.


\paragraph{Part (ii)}
Iterating~\cref{eq:component-wise} we have
\begin{align}
	\norm{x_i^\ell - x_i^*}_i &\leq \prod_{h = 0}^{\ell-1} \zeta^{u_i^\ell} \norm{x_i^0 - x_i^*}_i + \sum_{h = 0}^{\ell-1} \prod_{j = h+1}^{\ell-1} \zeta^{u_i^j} u_i^h \norm{e_i^h}_i \nonumber \\
	&\overset{\text{(i)}}{=} \zeta^{\beta_i(\ell)} \norm{x_i^0 - x_i^*}_i + \sum_{h = 0}^{\ell} \zeta^{(\sum_{j = h+1}^{\ell} u_i^j)} u_i^h \norm{e_i^h}_i \label{eq:component-wise-iterated}
\end{align}
where in (i) we defined $\beta_i(\ell) = \sum_{j = 0}^{\ell-1} u_i^j$, which is a binomial with $\beta_i(\ell) \sim \mathcal{B}(\ell, p_i)$.

Let us now define the sub-sequence $\{ h_k \}_{k = 0}^{\beta_i(\ell)-1}$ with $u_i^{h_k} = 1$ for $k = 0, \ldots, \beta_i(\ell)-1$, that is, the indices of the iterations that see an update of the $i$-th coordinate. We can then rewrite the RHS of~\cref{eq:component-wise-iterated} as follows:
$$
	\sum_{h = 0}^{\ell-1} \zeta^{(\sum_{j = h+1}^{\ell-1} u_i^j)} u_i^h \norm{e_i^h}_i = \sum_{k = 0}^{\beta_i(\ell)-1} \zeta^{(\sum_{j = h_k+1}^{\ell-1} u_i^j)} \norm{e_i^{h_k}}_i.
$$
By definition of $\{ h_k \}_{k = 0}^{\beta_i(\ell)-1}$ we know that between the times $h_k + 1$ and $\ell$ the total number of updates that are performed is $\beta_i(\ell) - (k+1)$. As a consequence we rewrite~\cref{eq:component-wise-iterated} as
\begin{align}
	\norm{x_i^\ell - x_i^*}_i &\leq \zeta^{\beta_i(\ell)} \norm{x_i^0 - x_i^*}_i + \sum_{k = 0}^{\beta_i(\ell)-1} \zeta^{\beta_i(\ell) - k - 1} \norm{e_i^{h_k}}_i \nonumber \\
	&\leq \zeta^{\beta_i(\ell)} \norm{x_i^0 - x_i^*}_i + \sum_{k = 0}^{\ell-1} \zeta^k \norm{e_i^{\beta_i(\ell) - 1 -h_k}}_i \label{eq:component-wise-iterated-2}
\end{align}
where we added $\ell - \beta(\ell) - 1$ terms to the sum to remove the dependence on $\beta_i(\ell)$, with the additional additive errors having the same distribution as $\norm{e_i^{\beta_i(\ell)}}_i$.

We observe now that $\zeta^{\beta_i(\ell)}$ is $\mathrm{subW}(1/2, \eta_i(\ell))$, where $\eta_i(\ell)$ is a monotonically decreasing function of $\ell$, see the discussion in \cref{ssubsec:zeta-beta-subweibull}. This fact, combined with \cref{as:stochastic-framework}(ii) for which $\norm{e_i}_i \sim \mathrm{subW}(\theta, \nu_i^\ell)$ and using \cref{pr:closure-sub-weibull-class}, shows that the RHS in~\cref{eq:component-wise-iterated-2} is sub-Weibull with parameters
$$
	\theta' = \max\{ 1/2, \theta \} \quad \text{and} \quad \nu_i'(\ell) = \eta_i(\ell) \norm{x_i^0 - x_i^*}_i + \frac{1 - \zeta^\ell}{1 - \zeta} \sup_\ell \nu_i^\ell,
$$
where we used the fact that by \cref{pr:inclusion-sub-weibull} $\norm{e_i^\ell}_i \sim \mathrm{subW}(\theta, \nu_i^\ell)$ implies $\norm{e_i^\ell}_i \sim \mathrm{subW}(\theta, \sup_\ell \nu_i^\ell)$. Finally, using \cref{lem:high-probability-bound} the thesis follows.
\end{proof}


The following result studies the convergence in mean and high probability under the assumption that the operator is quasi-averaged.

\begin{proposition}[Averaged case]\label{pr:averaged}
Let \cref{as:stochastic-framework,as:averaged} hold, and let $\{ \mathbold{x}^\ell \}_{\ell \in \mathbb{N}}$ be the trajectory generated by the stochastic Banach-Picard.

For each coordinate, $i = 1, \ldots, n$, the mean cumulative fixed point residual is upper bounded as follows:
\begin{align*}
	&\ev{\frac{1}{\ell+1} \sum_{h = 0}^\ell u_i^h \norm{(\mathcal{I}_i - \mathcal{T}_i) \mathbold{x}^h}_i^2} \leq \\ &\qquad \frac{\alpha}{1 - \alpha} \Bigg( \frac{1}{\ell+1} \norm{x_i^0 - x_i^*}_i^2 + p_i \sup_\ell \left(\ev{\norm{e_i^\ell}_i^2} + 2 \diam(\mathbb{D}_i) \mu_i^\ell \right) \Bigg).
\end{align*}

Moreover, with probability $1 - \delta$, $\delta \in (0, 1)$, we have, for any $i = 1, \ldots, n$:
\begin{align*}
	\frac{1}{\ell+1} \sum_{h = 0}^\ell u_i^h \norm{(\mathcal{I}_i - \mathcal{T}_i) \mathbold{x}^h}_i^2 &\leq \frac{\alpha}{1 - \alpha} \Bigg( \frac{1}{\ell+1} \norm{x_i^0 - x_i^*}_i^2 + \\ &+ \log^{2\theta}(2 / \delta) c(2\theta) \sup_\ell \left( 2^{2\theta} (\nu_i^\ell)^2 + 2 \diam(\mathbb{D}_i) \nu_i^\ell \right) \Bigg)
\end{align*}
where $c(2\theta)$ is defined in \cref{subsec:background-probability}.
\end{proposition}
\begin{proof}
Let $\mathbold{x}^* \in \fix(\mathcal{T})$; we have for each coordinate $i = 1, \ldots, n$:
\begin{align}
	\norm{x_i^{\ell+1} - x_i^*}_i^2 &\overset{\text{(i)}}{\leq} \norm{(1 - u_i^\ell) x_i^\ell + u_i^\ell \mathcal{T}_i \mathbold{x}^\ell - x_i^*}_i^2 + \nonumber \\ &\qquad + u_i^\ell \left( \norm{e_i^\ell}_i^2 + 2 \langle (1 - u_i^\ell) x_i^\ell + u_i^\ell \mathcal{T}_i \mathbold{x}^\ell - x_i^*, e_i^\ell \rangle_i \right)\nonumber \\
	&\overset{\text{(ii)}}{\leq} \norm{(1 - u_i^\ell) x_i^\ell + u_i^\ell \mathcal{T}_i \mathbold{x}^\ell - x_i^*}_i^2 + u_i^\ell \left( \norm{e_i^\ell}_i^2 + 2 \diam(\mathbb{D}_i) \norm{e_i^\ell}_i \right) \nonumber \\
	&\overset{\text{(iii)}}{\leq} \norm{x_i^\ell - x_i^*}_i^2 - u_i^\ell \frac{1 - \alpha}{\alpha} \norm{(\mathcal{I}_i - \mathcal{T}_i) \mathbold{x}^\ell}_i^2 + u_i^\ell m_i^\ell \label{eq:error-averaged}
\end{align}
where (i) follows by definition of norm and using the fact that $u_i^\ell = (u_i^\ell)^2$, (ii) holds by Cauchy-Schwarz inequality and boundedness of $\mathbb{D}_i$, and (iii) by the quasi-averagedness of $\mathcal{T}$ and defining $m_i^\ell := \norm{e_i^\ell}_i^2 + 2 \diam(\mathbb{D}_i) \norm{e_i^\ell}_i$.

Rearranging~\cref{eq:error-averaged} with $u_i^\ell \frac{1 - \alpha}{\alpha} \norm{(\mathcal{I}_i - \mathcal{T}_i) \mathbold{x}^\ell}_i^2$ on the LHS, and averaging over time we get:
\begin{equation}\label{eq:cumulative-fpr}
	\frac{1}{\ell+1} \sum_{h = 0}^\ell u_i^h \norm{(\mathcal{I}_i - \mathcal{T}_i) \mathbold{x}^h}_i^2 \leq \frac{1}{\ell+1} \frac{\alpha}{1 - \alpha} \left( \norm{x_i^0 - x_i^*}_i^2 + \sum_{h = 0}^\ell u_i^h m_i^h \right)
\end{equation}
where all error terms on the RHS cancel out, except for the initial one $\norm{x_i^0 - x_i^*}_i^2$, and  for $- \norm{x_i^{\ell+1} - x_i^*}_i^2$, which we removed for simplicity.


\paragraph{Convergence in mean}
Taking the expected value of~\cref{eq:cumulative-fpr} we get
$$
	\ev{\frac{1}{\ell+1} \sum_{h = 0}^\ell u_i^h \norm{(\mathcal{I}_i - \mathcal{T}_i) \mathbold{x}^h}_i^2} \leq \frac{1}{\ell+1} \frac{\alpha}{1 - \alpha} \left( \norm{x_i^0 - x_i^*}_i^2 + p_i \sum_{h = 0}^\ell \ev{m_i^h} \right)
$$
and upper bounding $\ev{m_i^h} \leq \sup_\ell \ev{m_i^\ell}$ the thesis follows.


\paragraph{Convergence in high probability}
In~\cref{eq:cumulative-fpr} we see that the first term on the RHS is deterministic, so we need only bound the second term. By \cref{pr:closure-sub-weibull-class} we know that $m_i^h \sim \mathrm{subW}(2\theta, 2^{2\theta} (\nu_i^h)^2 + 2 \diam(\mathbb{D}_i) \nu_i^h)$, which, by \cref{pr:inclusion-sub-weibull} implies $m_i^h \sim \mathrm{subW}(2\theta, \sup_\ell (2^{2\theta} (\nu_i^\ell)^2 + 2 \diam(\mathbb{D}_i) \nu_i^\ell))$. Finally, using \cref{pr:closure-sub-weibull-class} we easily see that
$$
	\frac{1}{\ell+1} \sum_{h = 0}^\ell u_i^h m_i^h \leq \frac{1}{\ell+1} \sum_{h = 0}^\ell m_i^h \sim \mathrm{subW}\left( 2\theta, \sup_\ell (2^{2\theta} (\nu_i^\ell)^2 + 2 \diam(\mathbb{D}_i) \nu_i^\ell) \right)
$$
and by \cref{lem:high-probability-bound} the thesis follows.
\end{proof}


\subsubsection{Sub-Weibull characterization of $\zeta^{\beta_i(\ell)}$}\label{ssubsec:zeta-beta-subweibull}
The goal of this section is to characterize the r.v. $\zeta^{\beta_i(\ell)}$ as a sub-Weibull, where we define $\beta_i(\ell) := \sum_{h = 0}^{\ell - 1} u_i^h$. This characterization is instrumental in proving \cref{pr:contractive} and, later, \cref{pr:dynamic-contractive}. Our goal then is to find $\theta'$ and $\eta_i(\ell)$ such that $\zeta^{\beta_i(\ell)} \sim \mathrm{subW}(\theta', \eta_i(\ell))$.

First of all, we can see that $\beta(\ell) \sim \mathcal{B}(\ell, p_i)$, since it is the sum of $\ell$ Bernoulli trials with probability $p_i$. Moreover, by the fact that $\zeta \in (0, 1)$, we know that $\zeta^{\beta_i(\ell)} \in [\zeta^\ell, 1]$, which means that $\zeta^{\beta_i(\ell)}$ is a bounded r.v.. As a consequence, we can model $\zeta^{\beta_i(\ell)}$ as a sub-Gaussian r.v., that is, a sub-Weibull with $\theta' = 1/2$, as proved in \cref{lem:bounded-rv}. On the other hand, by using the particular distribution of $\beta_i(\ell)$ we can give the following result.

\begin{lemma}
Let $\beta_i(\ell) \sim \mathcal{B}(p_i, \ell)$ and $\zeta \in (0, 1)$, then for any $k \geq 1$:
\begin{equation}\label{eq:k-norm-zeta_beta}
	\norm{\zeta^{\beta_i(\ell)}}_k = \left( 1 - p_i + p_i \zeta^k \right)^{\ell / k}.
\end{equation}
\end{lemma}
\begin{proof}
By the binomial distribution of $\beta_i(\ell)$ we have
\begin{align*}
	\norm{\zeta^{\beta_i(\ell)}}_k^k &= \ev{(\zeta^k)^{\beta_i(\ell)}} = \sum_{j = 0}^\ell (\zeta^k)^j \binom{\ell}{j} p_i^j (1 - p_i)^{\ell - j} \\
	&= \sum_{j = 0}^\ell \binom{\ell}{j} \left( \zeta^k p_i \right)^j (1 - p_i)^{\ell - j} = (1 - p_i + p_i \zeta^k)^\ell
\end{align*}
where the last equality holds by the binomial theorem. Taking the $k$-th root yields the thesis.
\end{proof}

We can clearly see from~\cref{eq:k-norm-zeta_beta} that the $k$-norm of $\zeta^{\beta_i(\ell)}$, for a fixed $k$, decays to zero as $\ell \to \infty$. Therefore we expect the sub-Weibull parameter $\eta_i(\ell)$ to be a decreasing function of $\ell$ as well. On the other hand, if we fix a finite $\ell \in \mathbb{N}$, we can see that the maximum value of $\norm{\zeta^{\beta_i(\ell)}}_k$ is $1$, attained when $k \to \infty$

However, a precise theoretical expression of $\eta_i(\ell)$ is, to the best of our knowledge, not possible. On the other hand, it is straightforward to characterize its value numerically. By the definition~(ii) in \cref{pr:equivalent-sub-weibull} of sub-Weibulls, we must have
$$
	\eta_i(\ell) \geq \frac{\left( 1 - p_i + p_i \zeta^k \right)^{\ell / k}}{\sqrt{k}}, \quad \forall k \geq 1.
$$
As a consequence of the observations above, $\eta_i(\ell)$ is a decreasing function of $\ell$, and it has a value in $[0, 1)$, where $0$ is attained when $k \to \infty$. Therefore, for any given $\ell$ we can choose $\eta_i(\ell) = \max_k \left\{ \left( 1 - p_i + p_i \zeta^k \right)^{\ell / k} / \sqrt{k} \right\} \in [0, 1)$. And, although this value does not have a closed form expression, numerically it is always possible to compute it.


\paragraph{An alternative approach}
The characterization above of $\zeta^{\beta_i(\ell)}$ as a sub-Weibull r.v. does not have a closed form, which may be unpalatable to some readers. Therefore, in the following we discuss an alternative approach that gives an exact value for $\eta_i(\ell)$, as well as an alternative to \cref{pr:contractive}.

As we observed above, $\zeta^{\beta_i(\ell)}$ is a bounded r.v. in $[\zeta^\ell, 1]$, and so according to \cref{lem:bounded-rv} it is such that $\zeta^{\beta_i(\ell)} - \chi_i^\ell \sim \mathrm{subW}(1/2, (1 - \zeta^\ell) / \sqrt{2})$. Using this fact, the following result holds.

\begin{proposition}[Contractive case -- alternative]\label{pr:alternative-contractive}
Let \cref{as:stochastic-framework,as:contractive} hold, and let $\{ \mathbold{x}^\ell \}_{\ell \in \mathbb{N}}$ be the trajectory generated by the stochastic Banach-Picard.

With probability $1 - \delta$, $\delta \in (0, 1)$, we have:
\begin{align*}
	\norm{x_i^\ell - x_i^*}_i &\leq \chi_i^\ell \norm{x_i^0 - x_i^*}_i + \frac{1 - \chi_i^\ell}{1 - \chi_i} \sup_\ell \mu_i^\ell + \\ &+ \log^{\theta'}(2/\delta) \ c(\theta') \left( \frac{1 - \zeta^\ell}{\sqrt{2}} \norm{x_i^0 - x_i^*}_i + 2 \frac{1 - \zeta^\ell}{1 - \zeta} \sup_\ell \nu_i^\ell \right)
\end{align*}
where $\theta' = \max\{ 1/2, \theta \}$ and $c(\theta')$ is defined in \cref{subsec:background-probability}.
\end{proposition}
\begin{proof}
By the inequality~\cref{eq:component-wise-iterated-2} from the proof of \cref{pr:contractive} we have
\begin{equation}\label{eq:component-wise-iterated-3}
	\norm{x_i^\ell - x_i^*}_i \leq \zeta^{\beta_i(\ell)} \norm{x_i^0 - x_i^*}_i + \sum_{k = 0}^{\ell-1} \zeta^k \norm{e_i^{\beta_i(\ell) - 1 -h_k}}_i.
\end{equation}
Now, by \cref{lem:centered-rv} we know that $\norm{e_i^\ell}_i \sim \mathrm{subW}(\theta, \nu_i^\ell)$ implies $\norm{e_i}_i^\ell - \mu_i^\ell \sim \mathrm{subW}(\theta, 2 \nu_i^\ell)$. Moreover, by \cref{lem:bounded-rv} we know that $\zeta^{\beta_i(\ell)} - \chi_i^\ell \sim \mathrm{subW}(1/2, (1 - \zeta^\ell) / \sqrt{2})$. As a consequence
$$
	\cref{eq:component-wise-iterated-3} - \ev{\cref{eq:component-wise-iterated-3}} \sim \mathrm{subW}\left( \max\{ 1/2, \theta \}, \frac{1 - \zeta^\ell}{\sqrt{2}} \norm{x_i^0 - x_i^*}_i + 2  \frac{1 - \chi_i^\ell}{1 - \chi_i} \sup_\ell \nu_i^\ell \right)
$$
and using \cref{lem:high-probability-bound} the thesis follows.
\end{proof}

\cref{pr:alternative-contractive} states that the distance from the fixed error $\norm{x_i^\ell - x_i^*}_i$ is upper bounded by the evolution of the mean distance plus a term due to the fact that the distance is sub-Weibull. We remark first of all that this bound involves the mean because of \cref{lem:bounded-rv} applied to $\zeta^{\beta_i(\ell)}$, which is stated for centered r.v.s. Secondly, we observe that the second part of the bound is similar to the high probability bound of \cref{pr:contractive}, with the difference that \emph{the dependence on the initial condition does not decay to zero as $\ell \to \infty$}. This is the reason why the sub-Weibull characterization of $\zeta^{\beta_i(\ell)}$ discussed above is preferable, since it correctly reflects the fact that the initial condition does not affect the asymptotic behavior.


\subsection{Discussion of convergence results}\label{subsec:discussion-convergence}

\subsubsection{Interpretation of high probability convergence results}
For simplicity we discuss the result of \cref{pr:contractive}, but the same reasoning holds for \cref{pr:averaged}.

The high probability bound of \cref{pr:contractive} states that, w.p. $1 - \delta$, the error satisfies:
$$
	\norm{x_i^\ell - x_i^*}_i \leq \log^{\theta'}(2/\delta) \ c(\theta') \left( \eta_i(\ell) \norm{x_i^0 - x_i^*}_i + \frac{1 - \zeta^\ell}{1 - \zeta} \sup_\ell \nu_i^\ell \right).
$$
We can see that the RHS is multiplied by $\log^{\theta'}(2/\delta)$, which grows as $\delta$ decreases. This implies that the more confidence we want in the bound (which requires $\delta$ to be smaller), the looser the bound becomes. Intuitively, smaller $\delta$ requires that we enlarge the bound to include more trajectory realizations $\{ x_i^\ell \}_{\ell \in \mathbb{N}}$.

\noindent Similarly, if $\theta' = \max\{ 1/2, \theta \}$ is larger, then the bound becomes looser. This is a consequence of the fact that the larger $\theta$ is, the heavier the tails of the additive noise are.

We can further observe that the high probability bound (multiplicative factors notwithstanding) has a similar structure to the mean bound:
$$
	\ev{\norm{x_i^\ell - x_i^*}_i} \leq \chi_i^\ell \norm{x_i^0 - x_i^*}_i + p_i \frac{1 - \chi_i^\ell}{1 - \chi_i} \sup_\ell \mu_i^\ell.
$$
Indeed, both bounds have a first term that depends on the initial condition and decays to zero as $\ell \to \infty$, and a second term that bounds the asymptotic distance from the fixed point. This second term depends on the additive noise, either via the bound on the mean $\sup_\ell \mu_i^\ell$ or the bound on its sub-Weibull parameter $\sup_\ell \nu_i^\ell$.

\noindent Interestingly, the convergence rate of the mean error is $\chi_i > \zeta$, while the convergence rate appearing in the high probability bound is $\zeta$ itself.


\paragraph{Sub-Weibull \emph{v.} Markov's inequality}
The fact that the dependence on $1 / \delta$ appears through its logarithm in the bound is an important motivation for the use of the proposed sub-Weibull framework. This feature of the bounds ensures that the RHS grows relatively slowly when we ask for increasing confidence (that is, $\delta \to 0$).

\noindent Consider the following alternative high probability bound, which is based on Markov's inequality: with this approach, the dependence on the RHS is with $1 / \delta$ and not its logarithm. Although Markov's inequality holds for a more general class of random variables than sub-Weibull, this result makes it clear that using the sub-Weibull framework allows to derive sharper bounds, while still considering a number of distributions as sub-cases. 

\begin{lemma}
Let \cref{as:stochastic-framework,as:contractive} hold, and let $\{ \mathbold{x}^\ell \}_{\ell \in \mathbb{N}}$ be the sequence generated by~\cref{eq:stochastic-update}. The following holds. 

For each coordinate, $i = 1, \ldots, n$, with probability $1 - \delta$, $\delta \in (0, 1)$, we have that
\begin{equation}\label{eq:bound-using-markov}
    \norm{x_i^\ell - x_i^*}_i \leq \frac{1}{\delta} \left( \chi_i^\ell \norm{x_i^0 - x_i^*}_i + p_i \frac{1 - \chi_i^\ell}{1 - \chi_i} \sup_\ell \mu_i^\ell \right)
\end{equation}
where $\chi_i := 1 - p_i + p_i \zeta \geq \zeta$.
\end{lemma}
\begin{proof}
By \cref{pr:contractive} we have
$$
    \ev{\norm{x_i^\ell - x_i^*}_i} \leq \chi_i^\ell \norm{x_i^0 - x_i^*}_i + p_i \frac{1 - \chi_i^\ell}{1 - \chi_i} \sup_\ell \mu_i^\ell.
$$
Then, by Markov's inequality we have that $\pr{\norm{x_i^\ell - x_i^*}_i \geq \ev{\norm{x_i^\ell - x_i^*}_i} / \delta} \leq \delta$, which yields~\cref{eq:bound-using-markov}.
\end{proof}


\subsubsection{Difference from deterministic convergence}

\paragraph{Contractive case}
In the contractive case, we know by \cite[Theorem~1.50]{bauschke_convex_2017} that the deterministic Banach-Picard converges linearly to the fixed point. In particular, we have that
$$
	\norm{x_i^\ell - x_i^*} \leq \zeta^\ell \norm{x_i^0 - x_i^*}.
$$
Comparing this fact with \cref{pr:contractive} we notice that the introduction of random coordinate updates degrades the convergence rate from $\zeta$ to $\chi_i \geq \zeta$, in accordance with the results of \cite{combettes_stochastic_2019}. Moreover, the introduction of additive errors implies that we do not reach exact convergence to the fixed point, but rather to a neighborhood thereof. We will return to this point in \cref{subsec:hp-vs-as}.

A further observation is that the additive error may result in an \emph{expansive} stochastic update, that is, one for which at some $\ell \in \mathbb{N}$:
$$
	\norm{x_i^{\ell+1} - x_i^*}_i \geq \norm{x_i^\ell - x_i^*}_i.
$$
Nonetheless, the underlying contractiveness of the (exact) operator ensures that expansive updates yield an inexact -- but bounded -- convergence, instead of an unbounded one.


\paragraph{Averaged case}
The convergence of~\cref{eq:stochastic-update} when $\mathcal{T}$ is averaged can be proved only in terms of the \emph{cumulative fixed point residual}, due to this property being weaker than contractiveness. The cumulative FPR, defined, recall, as
$$
	\frac{1}{\ell+1} \sum_{h = 0}^\ell u_i^h \norm{(\mathcal{I}_i - \mathcal{T}_i) \mathbold{x}^h}_i^2
$$
averages the FPR achieved by each term in the trajectory $\{ x_i^\ell \}_{\ell \in \mathbb{N}}$, where of course a smaller FPR implies that a point $x_i^\ell$ is closer to being a fixed point of $\mathcal{T}_i$.

We recall that for a (non-stochastic) Banach-Picard applied to $\mathcal{T}$ averaged, we can prove that $\{ \norm{(\mathcal{I} - \mathcal{T}) \mathbold{x}^\ell} \}_{\ell \in \mathbb{N}}$ \cite[Theorem~5.15]{bauschke_convex_2017} is a monotonically decreasing sequence that converges to zero. This is no longer the case in the stochastic framework of this paper, due to the presence of additive errors. In this scenario, the FPR is no longer monotonically decreasing, and to account for this fact we analyze the overall evolution of the FPR.

We remark that the concept of cumulative FPR is similar to that of \emph{regret} in convex optimization, and in particular to that of dynamic regret in online optimization \cite{mokhtari_online_2016}, see also \cref{sec:dynamic}. Moreover, it includes as a particular case the different concept of regret based on the residual of the proximal gradient method proposed in \cite{hallak_regret_2020}.


\subsection{Convergence without additive error}
The results (both in mean and high probability) of \cref{pr:contractive,pr:averaged} hold even when the additive error is absent. However, in this particular scenario it is actually possible to give more refined bounds, as discussed in this section.

In particular, if the only source of stochasticity are random coordinate updates (that is, we choose $e_i = 0$ and consequently $\nu_i = 0$) we have the following result.

\begin{proposition}[Without additive noise]\label{pr:no-additive-noise}
Let \cref{as:stochastic-framework,as:contractive} hold, with $e_i = 0$ a.s. for all $i = 1, \ldots, n$; let $\{ \mathbold{x}^\ell \}_{\ell \in \mathbb{N}}$ be the trajectory generated by the stochastic Banach-Picard.

Let $\epsilon \in (0, p_i]$, then for each coordinate, $i = 1, \ldots, n$, with probability $1 - \delta(\epsilon, \ell)$ we have:
$$
	\norm{x_i^\ell - x_i^*}_i \leq \zeta^{\ell (p_i - \epsilon)} \norm{x_i^0 - x_i^*}_i
$$
where
\begin{equation}\label{eq:sanov-theorem-parameters}
\begin{split}
	\delta(\epsilon, \ell) &= \exp\left( - \ell D(p_i - \epsilon || p_i) \right), \\
	D(p_i - \epsilon || p_i) &= (p_i - \epsilon) \log\left( 1 - \frac{\epsilon}{p_i} \right) + (1 - p_i + \epsilon) \log\left( 1 + \frac{\epsilon}{1 - p_i} \right).
\end{split}
\end{equation}
\end{proposition}
\begin{proof}
Recalling the bound~\cref{eq:component-wise-iterated-2} from the proof of \cref{pr:contractive}, we have that
$$
	\norm{x_i^\ell - x_i^*}_i \leq \zeta^{\beta_i(\ell)} \norm{x_i^0 - x_i^*}_i
$$
where $\beta_i(\ell) \sim \mathcal{B}(\ell, p_i)$.

Using Sanov's theorem \cite[Theorem~D.3]{mohri_foundations_2018} (in particular its symmetric form in \cite[eq.~(D.7)]{mohri_foundations_2018}) we know that
$$
	\pr{\beta(\ell) \leq \ell (p_i - \epsilon)} \leq \exp\left( - \ell D(p_i - \epsilon || p_i) \right)
$$
and, since $\pr{\zeta^{\beta_i(\ell)} \geq \zeta^{\ell (p_i - \epsilon)}} = \pr{\beta(\ell) \leq \ell (p_i - \epsilon)}$, the thesis follows.
\end{proof}

We observe in \cref{pr:no-additive-noise} that the convergence rate is characterized by $\zeta^{p_i - \epsilon}$, which is larger than $\zeta$, that is, the convergence rate achieved by the non-stochastic Banach-Picard \cite[Theorem~1.50]{bauschke_convex_2017}. This is similar to the result of \cref{pr:contractive} where the convergence rate is $\chi_i \geq \zeta$. Notice also that since $\epsilon \in (0, p_i]$ then $\zeta^{p_i-\epsilon} < \chi_i$.

\noindent We further notice that $\delta(\epsilon, \ell)$ is a decreasing function of $\ell$, with $\delta(\epsilon, \ell) \to 0$ as $\ell \to \infty$. Thus this implies that asymptotically, the bound holds with probability $\lim_{\ell \to \infty} 1 - \delta(\epsilon, \ell) = 1$, thus \emph{yielding the known almost sure convergence results}, \emph{e.g.} of \cite{combettes_stochastic_2015,bianchi_coordinate_2016}. But \cref{pr:no-additive-noise} provides additional information w.r.t. the asymptotic results \cite{combettes_stochastic_2015,bianchi_coordinate_2016}, since the bound holds also for any finite $\ell$, albeit within a confidence level tuned by $\epsilon$ and not almost surely.

A similar result can be derived for the averaged case.

\begin{proposition}[Without additive noise -- averaged]\label{pr:no-additive-noise-averaged}
Let \cref{as:stochastic-framework,as:averaged} hold, with $e_i = 0$ a.s. for all $i = 1, \ldots, n$; let $\{ \mathbold{x}^\ell \}_{\ell \in \mathbb{N}}$ be the trajectory generated by the stochastic Banach-Picard.

Let $\epsilon \in (0, p_i]$, then for each coordinate, $i = 1, \ldots, n$, with probability $1 - \delta(\epsilon, \ell+1)$ we have:
$$
	\norm{(\mathcal{I}_i - \mathcal{T}_i) \mathbold{x}^\ell}_i^2 \leq \frac{1}{(\ell + 1) (p_i - \epsilon)} \frac{\alpha}{1 - \alpha} \norm{x_i^0 - x_i^*}_i^2
$$
where $\delta(\epsilon, \ell+1)$ and $D(p_i - \epsilon || p_i)$ are defined as in~\cref{eq:sanov-theorem-parameters}.
\end{proposition}
\begin{proof}
First of all, from~\cref{eq:error-averaged} in the proof of \cref{pr:averaged} we can see that, if $u_i^\ell = 1$, then
$$
	\norm{x_i^{\ell+1} - x_i^*}_i^2 \leq \norm{x_i^\ell - x_i^*}_i^2 - \frac{1 - \alpha}{\alpha} \norm{(\mathcal{I}_i - \mathcal{T}_i) \mathbold{x}^\ell}_i^2,
$$
while if $u_i^\ell = 0$ we have $\norm{x_i^{\ell+1} - x_i^*}_i^2 = \norm{x_i^\ell - x_i^*}_i^2$. This implies that the operator is \emph{stochastic Fej\'er monotone}, see \cite{combettes_stochastic_2015,combettes_stochastic_2019}, as well as \cite{bauschke_convex_2017} for the notion of deterministic Fej\'er monotonicity; as a consequence, the sequence of fixed point residuals $\{ \norm{(\mathcal{I}_i - \mathcal{T}_i) \mathbold{x}^\ell}_i \}_{\ell \in \mathbb{N} | u_i^\ell = 1}$ is monotonically decreasing.

Therefore, taking~\cref{eq:cumulative-fpr} and using the Fej\'er monotonicity we can write
$$
	\beta_i(\ell+1) \norm{(\mathcal{I}_i - \mathcal{T}_i) \mathbold{x}^\ell}_i^2 \leq \sum_{h = 0}^\ell u_i^h \norm{(\mathcal{I}_i - \mathcal{T}_i) \mathbold{x}^h}_i^2 \leq \frac{\alpha}{1 - \alpha} \norm{x_i^0 - x_i^*}_i^2
$$
where we used the fact that the sequence $\{ \norm{x_i^{\ell+1} - x_i^\ell} \}_{\ell \in \mathbb{N} | u_i^\ell = 1}$ has $\beta_i(\ell+1)$ non-zero terms. Dividing by $\beta_i(\ell+1)$ on both sides we get
$$
	\norm{x_i^{\ell+1} - x_i^\ell}_i^2 \leq \frac{1}{\beta_i(\ell+1)} \frac{\alpha}{1 - \alpha} \norm{x_i^0 - x_i^*}_i^2.
$$
We have now the following fact
\begin{align*}
	&\pr{\frac{1}{\beta_i(\ell+1)} \norm{x_i^0 - x_i^*}_i^2 \geq \frac{1}{(p_i - \epsilon) (\ell + 1)} \norm{x_i^0 - x_i^*}_i^2} = \\ &\hspace{3cm} = \pr{\beta_i(\ell+1) \leq (\ell + 1) (p_i - \epsilon)} \overset{\text{(i)}}{\leq} \exp\left( - (\ell + 1) D(p_i - \epsilon || p_i)\right)
\end{align*}
where (i) holds by Sanov's theorem \cite[Theorem~D.3]{mohri_foundations_2018} (cf. \cite[eq.~(D.7)]{mohri_foundations_2018}), and the thesis follows.
\end{proof}

Recalling that the fixed point residual in a deterministic scenario convergences as \cite{davis_convergence_2016}
$$
	\norm{(\mathcal{I}_i - \mathcal{T}_i) \mathbold{x}^\ell}_i^2 \leq \frac{1}{\ell+1} \frac{\alpha}{1 - \alpha} \norm{x_i^0 - x_i^*}_i^2
$$
we see once again that introducing random coordinate updates -- which implies that the RHS is divided by $p_i - \epsilon < 1$ -- degrades the speed of convergence.

Similarly to the contractive case, a consequence of this result is the almost sure convergence, with the difference that \cref{pr:no-additive-noise-averaged} provides a non-asymptotic, high probability bound as well, differently from \cite{combettes_stochastic_2015}.


\subsection{High probability \emph{v.} almost sure convergence}\label{subsec:hp-vs-as}
The previous sections presented convergence results in terms of mean and high probability bounds, and only in the absence of additive error could we derive almost sure convergence results. In this section then we discuss the question: \emph{can we achieve almost sure convergence in the framework defined by~\cref{eq:stochastic-update}?}

First of all, consider the quasi-contractive scenario of \cref{as:contractive}. By the proof of \cref{pr:contractive} we know that $\norm{x_i^{\ell+1} - x_i^*}_i \leq \zeta^{u_i^\ell} \norm{x_i^\ell - x_i^*}_i + u_i^\ell \norm{e_i^\ell}_i$, which implies that the always present source of randomness $\norm{e_i^\ell}_i$ would not allow $\norm{x_i^{\ell} - x_i^*}_i$ to settle on a single value -- and hence $x_i^\ell$ cannot reach a fixed point.

The following result, on the other hand, shows that we do achieve almost sure convergence to a neighborhood of the fixed point.

\begin{corollary}[Convergence to neighborhood]
Let \cref{as:stochastic-framework,as:contractive} hold, and let $\{ \mathbold{x}^\ell \}_{\ell \in \mathbb{N}}$ be the sequence generated by~\cref{eq:stochastic-update}.

\noindent Then for each coordinate, $i = 1, \ldots, n$ it holds that:
$$
	\limsup_{\ell \to \infty} \norm{x_i^\ell - x_i^*}_i \leq \frac{\sup_\ell \mu_i^\ell}{1 - \zeta} \quad \text{a.s.}.
$$
\end{corollary}
\begin{proof}
By \cref{pr:contractive} we know that
$$
	\ev{\norm{x_i^\ell - x_i^*}_i} \leq \chi_i^\ell \norm{x_i^0 - x_i^*}_i + p_i \frac{1 - \chi_i^\ell}{1 - \chi_i} \sup_\ell \mu_i^\ell,
$$
and, defining the random variable
$
    y_i^\ell := \max\left\{ 0, \norm{x_i^\ell - x_i^*}_i - p_i \frac{1 - \chi_i^\ell}{1 - \chi_i} \sup_\ell \mu_i^\ell \right\},
$
this fact implies $\ev{y_i^\ell} \leq \chi_i^\ell \norm{x_i^0 - x_i^*}_i$.

By Markov's inequality then we have that, for any $\epsilon > 0$:
$$
    \pr{y_i^\ell \geq \epsilon} \leq \frac{\ev{y_i^\ell}}{\epsilon} \leq \frac{\norm{x_i^0 - x_i^*}_i}{\epsilon} \chi_i^\ell
$$
and we further have:
$$
    \sum_{\ell = 0}^\infty \pr{y_i^\ell \geq \epsilon} \leq \frac{\norm{x_i^0 - x_i^*}_i}{\epsilon} \frac{1}{1 - \chi_i} < \infty.
$$
By the Borel-Cantelli lemma, this fact implies that almost surely
$
    \limsup_{\ell \to \infty} y_i^\ell \leq \epsilon,
$
and, since the inequality holds for any $\epsilon > 0$, the thesis is proved.
\end{proof}


We conclude this section remarking that it is possible to prove almost sure convergence of the stochastic Banach-Picard in another particular case: when the additive noise is diminishing -- roughly, when the sequence $\{ \norm{e_i^h}_i \}_{h \in \mathbb{N}}$ is summable, as proved in \emph{e.g.} \cite{combettes_stochastic_2015,combettes_stochastic_2019}. The following example shows how the proposed framework can be applied to study this scenario.

\begin{corollary}[Vanishing additive error]\label{cor:vanishing-error}
Consider the stochastic Banach-Picard~\cref{eq:stochastic-update} under \cref{as:contractive}, and suppose that: (i) $p_i = 1$, $i = 1, \ldots, n$, and (ii) the additive errors are such that $\norm{e_i^\ell}_i \sim \mathrm{subW}(\theta, \nu_i^\ell)$, with $\nu_i^\ell \to 0$ as $\ell \to +\infty$; then:
$$
    \limsup_{\ell \to \infty} \norm{x_i^\ell - x_i^*}_i = 0 \quad \text{a.s.}.
$$
\end{corollary}
\begin{proof}
By \cref{pr:contractive} we know that $\norm{x_i^\ell - x_i^*}_i$ is upper bounded by a sub-Weibull with parameters $\theta$ and $N_i^\ell = \zeta^\ell \norm{x_i^0 - x_i^*}_i + \sum_{h = 0}^{\ell-1} \zeta^{\ell - h - 1} \nu_i^h$.

\noindent We can see that $\lim_{\ell \to \infty} N_i^\ell = 0$, where the second term converges to zero by using \cite[Lemma~3.1(a)]{sundharram_distributed_2010}. As a consequence $\lim_{\ell \to \infty} \norm{x_i^\ell - x_i^*} \sim \mathrm{subW}(\theta, 0)$, and, since $\norm{x_i^\ell - x_i^*} \geq 0$ for all $\ell \in \mathbb{N}$, this implies that almost surely $\lim_{\ell \to \infty} \norm{x_i^\ell - x_i^*} = 0$.
\end{proof}

From the discussion above, we can clearly see that almost sure convergence can be achieved -- either to a set, in the general case, or to the fixed point, in particular scenarios. However, this type of convergence only provides \emph{asymptotic} results, contrary to the high probability convergence analysis, which provides bounds that hold for all $\ell \in \mathbb{N}$.



\section{Convergence Analysis: Online Banach-Picard}\label{sec:dynamic}

In this section we analyze the convergence of the dynamic Banach-Picard introduced in \cref{subsec:dynamic-banach-picard-framework}, under the assumptions of quasi-contractiveness and quasi-averagedness.


\begin{proposition}[Dynamic contractive case]\label{pr:dynamic-contractive}
Let \cref{as:stochastic-framework,as:contractive,as:dynamic-contractive} hold, and let $\{ \mathbold{x}^\ell \}_{\ell \in \mathbb{N}}$ be the trajectory generated by the stochastic Banach-Picard.

For each coordinate, $i = 1, \ldots, n$, the mean tracking error is upper bounded as follows:
$$
	\ev{\norm{x_i^\ell - x_i^{*,\ell}}_i} \leq \chi_i^\ell \ev{\norm{x_i^0 - x_i^{*,0}}_i} + \frac{1 - \chi_i^\ell}{1 - \chi_i} \sup_\ell \left( p_i \mu_i^\ell + \ev{\sigma_i^\ell} \right)
$$
where recall that $\chi_i := 1 - p_i + p_i \zeta \geq \zeta$.

\noindent Moreover, with probability $1 - \delta$, $\delta \in (0, 1)$, we have:
$$
	\norm{x_i^\ell - x_i^{*,\ell}}_i \leq \log^{\theta'}(2/\delta) \ c(\theta') \left(
	\eta_i(\ell) \norm{x_i^0 - x_i^{*,0}}_i + \frac{1 - \zeta^\ell}{1 - \zeta} \sup_\ell \left( \nu_i^\ell + \gamma_i^\ell / p_i \right) \right)
$$
where $\theta' = \max\{ 1/2, \theta, \varphi \}$, $c(\theta')$ is defined in \cref{subsec:background-probability}, and $\eta_i(\ell)$ is a monotonically decreasing function of $\ell$, see \cref{ssubsec:zeta-beta-subweibull}.
\end{proposition}
\begin{proof}
Using the triangle inequality, coordinate-wise we have:
\begin{align}
	\norm{x_i^{\ell+1} - x_i^{*,\ell+1}}_i &\leq \norm{(1 - u_i^\ell) x_i^\ell + u_i^\ell \mathcal{T}_i^{\ell+1} \mathbold{x}^\ell - x_i^{*,\ell+1}}_i + u_i^\ell \norm{e_i^\ell}_i \nonumber \\
	&\overset{\text{(i)}}{\leq} \zeta^{u_i^\ell} \norm{x_i^\ell - x_i^{*,\ell+1}}_i + u_i^\ell \norm{e_i^\ell}_i \nonumber \\
	&\overset{\text{(ii)}}{\leq} \zeta^{u_i^\ell} \norm{x_i^\ell - x_i^{*,\ell}}_i + u_i^\ell \norm{e_i^\ell}_i + \sigma_i^\ell \label{eq:dynamic-component-wise}
\end{align}
where (i) follows by quasi-contractiveness of $\mathcal{T}_i^{\ell+1}$, and (ii) by the triangle inequality and the definition of distance between consecutive fixed points, $\sigma_i^\ell$.

\paragraph{Convergence in mean}
Taking the expectation of~\cref{eq:dynamic-component-wise} yields
\begin{equation}\label{eq:dynamic-component-wise-mean}
	\ev{\norm{x_i^{\ell+1} - x_i^{*,\ell+1}}_i} \leq \chi_i \ev{\norm{x_i^\ell - x_i^{*,\ell}}_i} + p_i \mu_i^\ell + \ev{\sigma_i^\ell}
\end{equation}
Recursively applying~\cref{eq:dynamic-component-wise-mean} and using the geometric sum yields the thesis.


\paragraph{Convergence in high probability}
Iterating~\cref{eq:dynamic-component-wise} we have
\begin{align}
	\norm{x_i^\ell - x_i^{*,\ell}}_i &\overset{\text{(i)}}{\leq} \zeta^{\beta_i(\ell)} \norm{x_i^0 - x_i^{*,0}}_i + \sum_{k = 0}^{\beta_i(\ell)-1} \zeta^k \norm{e_i^{\beta_i(\ell) - h_k - 1}}_i + \sum_{h = 0}^{\ell - 1} \zeta^{\sum_{j = h+1}^{\ell-1} u_i^j} \sigma_i^h \nonumber \\
	&\overset{\text{(ii)}}{\leq} \zeta^{\beta_i(\ell)} \norm{x_i^0 - x_i^{*,0}}_i + \sum_{k = 0}^{\ell-1} \zeta^k \norm{e_i^{\ell - h_k - 1}}_i + \sum_{h = 0}^{\ell - 1} \zeta^{\sum_{j = h+1}^{\ell-1} u_i^j} \sigma_i^h \label{eq:dynamic-component-wise-iterated}
\end{align}
where we used the fact that $u_i^\ell$ multiplies the error to simplify the first sum in (i), and in (ii) we upper bounded it by $\sum_{k = 0}^{\ell-1} \zeta^k \norm{e_i^{\ell - h_k - 1}}_i$.

Similarly to the static case, we can characterize the first two terms of~\cref{eq:dynamic-component-wise-iterated} as a sub-Weibull, but we need to do the same for the third term. First of all, we notice that $\sum_{h = 0}^{\ell - 1} \zeta^{\sum_{j = h+1}^{\ell-1} u_i^j} \sigma_i^h$ is \emph{not} a geometric series, since $u_i^\ell$ does not multiply $\sigma_i^\ell$. Rather, the sum can contain several $\sigma_i^h$'s that are multiplied by the same power of $\zeta$. To see this, we can write
$$
	\sum_{h = 0}^{\ell - 1} \zeta^{\sum_{j = h+1}^{\ell-1} u_i^j} \sigma_i^h = \sum_{k = 0}^{\beta_i(\ell) - 1} \zeta^k \sum_{h = h_k}^{h_{k+1} - 1} \sigma_i^h \leq \sum_{k = 0}^{\ell - 1}  \zeta^k \sum_{h = h_k}^{h_{k+1} - 1} \sigma_i^h
$$
which becomes a geometric series with the r.v.s $\sum_{h = h_k}^{h_{k+1} - 1} \sigma_i^h$ having different distributions, since they sum a variable number of $\sigma_i^h$'s.

Our goal then is to characterize $\sum_{h = h_k}^{h_{k+1} - 1} \sigma_i^h$ as a sub-Weibull. Using the law of total expectation we can write
\begin{align*}
	\norm{\sum_{h = h_k}^{h_{k+1} - 1} \sigma_i^h}_j &= \ev{\norm{\sum_{h = h_k}^{h_{k+1} - 1} \sigma_i^h \ \Bigg| \ h_k, h_{k+1} }_j} \overset{\text{(i)}}{\leq} \ev{\sum_{h = h_k}^{h_{k+1} - 1} \norm{\sigma_i^h \ \big| \ h_k, h_{k+1}}_j} \\
	&\overset{\text{(ii)}}{\leq} \ev{(h_{k+1} - h_k) \sup_\ell \gamma_i^\ell j^\varphi} \overset{\text{(iii)}}{=} \frac{\sup_\ell \gamma_i^\ell}{p_i} j^\varphi
\end{align*}
where (i) holds by triangle inequality, (ii) by sub-Weibull distribution of the $\sigma_i^h$'s, and by using the upper bound $\gamma_i^\ell \leq \sup_\ell \gamma_i^\ell$, and (iii) because $h_{k+1} - h_k$ is a geometric random variable with mean $1 / p_i$, since it models the consecutive number of missed updates, which are Bernoulli trials.

Overall we have that~\cref{eq:dynamic-component-wise-iterated} is sub-Weibull with
$$
	\cref{eq:dynamic-component-wise-iterated} \sim \mathrm{subW}\left( \max\{1/2, \theta, \varphi\}, \eta_i(\ell) \norm{x_i^0 - x_i^{*,0}}_i + \frac{1 - \zeta^\ell}{1 - \zeta} \sup_\ell \left( \nu_i + \gamma_i / p_i \right). \right)
$$
and using \cref{lem:high-probability-bound} the thesis follows.
\end{proof}

\begin{remark}[Path-length]
In online optimization, a widely used concept is that of \emph{path-length}, that is, the cumulative distance between consecutive optimizers. This concept can be straightforwardly extended to the operator theoretical set-up, by defining it as the cumulative distance between consecutive fixed points:
$$
	c^\ell := \sum_{h = 0}^{\ell-1} \ev{\norm{x_i^{*,h+1} - x_i^{*,h}}_i} \leq \ell \sup_\ell \ev{\sigma_i^\ell}.
$$
The path-length often appears on the RHS of regret bounds, see \emph{e.g.} \cite{mokhtari_online_2016}, but notice that it \emph{grows linearly with $\ell$}. What instead appears in our bound is a \emph{weighted} path-length
$$
	c_\mathrm{w}^\ell := \sum_{h = 0}^{\ell-1} \chi_i^{\ell-h-1} \ev{\norm{x_i^{*,h+1} - x_i^{*,h}}_i} \leq \frac{1 - \chi_i^\ell}{1 - \chi_i} \sup_\ell \ev{\sigma_i^\ell}
$$
which asymptotically reaches a fixed value, independent of $\ell$. This is due to the fact that we use the quasi-contractiveness of the operator, allowing to reach a tighter bound.
\end{remark}


\begin{proposition}[Dynamic averaged case]\label{pr:dynamic-averaged}
Let \cref{as:stochastic-framework,as:averaged} hold, and let $\{ \mathbold{x}^\ell \}_{\ell \in \mathbb{N}}$ be the trajectory generated by the stochastic Banach-Picard.

For each coordinate, $i = 1, \ldots, n$, the mean cumulative fixed point residual is upper bounded as follows:
\begin{align*}
	&\ev{\frac{1}{\ell+1} \sum_{h = 0}^\ell u_i^h \norm{(\mathcal{I}_i - \mathcal{T}_i^{h+1}) \mathbold{x}^h}_i^2} \leq \frac{\alpha}{1 - \alpha} \Bigg( \frac{1}{\ell+1} \norm{x_i^0 - x_i^*}_i^2 + \\ &\qquad + \sup_\ell \left( p_i \left(\ev{\norm{e_i^\ell}_i^2} + 2 \diam(\mathbb{D}_i) \mu_i^\ell \right) + \ev{(\sigma_i^\ell)^2} + 2 \diam(\mathbb{D}_i) \ev{\sigma_i^\ell} \right) \Bigg)
\end{align*}
where $\sigma_i^h$ can be either $\underline{\sigma}_i^h$ or $\overline{\sigma}_i^h$.

Moreover, with probability $1 - \delta$, $\delta \in (0, 1)$, we have, for any $i = 1, \ldots, n$:
\begin{align*}
	\frac{1}{\ell+1} \sum_{h = 0}^\ell u_i^h & \norm{(\mathcal{I}_i - \mathcal{T}_i^{h+1}) \mathbold{x}^h}_i^2 \leq \frac{\alpha}{1 - \alpha} \Bigg( \frac{1}{\ell+1} \norm{x_i^0 - x_i^{*,0}}_i^2 + \\ +& \log^{\theta'}(2 / \delta) c(\theta') \sup_\ell \left( 2^{2\theta'} \left( (\nu_i^\ell)^2 + (\gamma_i^\ell)^2 \right) + 2\diam(\mathbb{D}_i) \left( \nu_i^\ell + \gamma_i^\ell \right) \right) \Bigg)
\end{align*}
where $\theta' = 2 \max\{ \theta, \varphi \}$ and $c(\theta')$ is defined in \cref{subsec:background-probability}, and where we choose the sub-Weibull distribution of either $\underline{\sigma}_i$ or $\overline{\sigma}_i$.
\end{proposition}
\begin{proof}
Let $\mathbold{x}^* \in \fix(\mathcal{T})$; similarly to the static case, we have for each coordinate $i = 1, \ldots, n$:
\begin{equation}\label{eq:dynamic-error-averaged}
\begin{split}
	\norm{x_i^{\ell+1} - x_i^{*,\ell+1}}_i^2 &\leq \norm{x_i^\ell - x_i^{*,\ell+1}}_i^2 - u_i^\ell \frac{1 - \alpha}{\alpha} \norm{(\mathcal{I}_i - \mathcal{T}_i^{\ell+1}) \mathbold{x}^\ell}_i^2 + \\ &\qquad + u_i^\ell \left( \norm{e_i^\ell}_i^2 + 2 \diam(\mathbb{D}_i) \norm{e_i^\ell}_i \right).
\end{split}
\end{equation}
Denoting $\sigma_i^\ell = \norm{x_i^{*,\ell+1} - x_i^{*,\ell}}_i$, and using Cauchy-Schwarz we have
\begin{equation}\label{eq:dynamic-error-bound}
	\norm{x_i^\ell - x_i^{*,\ell+1}}_i^2 \leq \norm{x_i^\ell - x_i^{*,\ell}}_i^2 + (\sigma_i^\ell)^2 + 2 \diam(\mathbb{D}_i) \sigma_i^\ell;
\end{equation}
$\sigma_i^\ell$ can be substituted by either $\underline{\sigma}_i^\ell$ or $\overline{\sigma}_i^\ell$ to give a minimal and maximal bound to the cumulative FPR. We substitute~\cref{eq:dynamic-error-bound} into~\cref{eq:dynamic-error-averaged}, and rearranging and averaging over time we get:
\begin{equation}\label{eq:dynamic-cumulative-fpr}
	\frac{1}{\ell+1} \sum_{h = 0}^\ell u_i^h \norm{(\mathcal{I}_i - \mathcal{T}_i^{h+1}) \mathbold{x}^h}_i^2 \leq \frac{1}{\ell+1} \frac{\alpha}{1 - \alpha} \left( \norm{x_i^0 - x_i^{*,0}}_i^2 + \sum_{h = 0}^\ell m_i^h \right)
\end{equation}
where we denote $m_i^\ell = u_i^\ell \left( \norm{e_i^\ell}_i^2 + 2 \diam(\mathbb{D}_i) \norm{e_i^\ell}_i \right) + (\sigma_i^\ell)^2 + 2 \diam(\mathbb{D}_i) \sigma_i^\ell$.


\paragraph{Convergence in mean}
Taking the expected value of~\cref{eq:dynamic-cumulative-fpr} we get
$$
	\ev{\frac{1}{\ell+1} \sum_{h = 0}^\ell u_i^h \norm{(\mathcal{I}_i - \mathcal{T}_i^{h+1}) \mathbold{x}^h}_i^2} \leq \frac{1}{\ell+1} \frac{\alpha}{1 - \alpha} \left( \norm{x_i^0 - x_i^*}_i^2 + \sum_{h = 0}^\ell \ev{m_i^h} \right)
$$
and hence the thesis.


\paragraph{Convergence in high probability}
By the properties of sub-Weibulls, we can see that $\sum_{h = 0}^\ell m_i^h / (\ell + 1)$ is sub-Weibull with parameters
$$
	\theta' = 2 \max\{ \theta, \varphi \} \quad \text{and} \quad \nu' = 2^{2\theta'} \sup_\ell \left( \left( (\nu_i^\ell)^2 + (\gamma_i^\ell)^2 \right) + 2\diam(\mathbb{D}_i) \left( \nu_i^\ell + \gamma_i^\ell \right) \right)
$$
and using \cref{lem:high-probability-bound} the thesis follows.
\end{proof}


\begin{remark}
We remark that in this framework the time variability of the operators is a source of error that is always present, even if the update is not performed. This models the fact that even when we do not perform an update, the operators are changing, and the tracking error will be worse the longer we go without updating. This is the case in time-varying optimization, where \emph{e.g.} data are streaming whether we observe them or not.

Alternatively, we could choose a framework in which the operator changes only when we apply an update; this would imply that in~\cref{eq:dynamic-component-wise} also $\sigma_i^\ell$ is multiplied by $u_i^\ell$. In this scenario, the error due to operator variability can be treated like the additive noise, and the results of the static case hold straightforwardly.
\end{remark}


\begin{remark}[Deterministic $\sigma_i$]
In \cref{as:dynamic-contractive,as:dynamic-averaged} we defined a stochastic model for the distance(s) between consecutive fixed points. This model is very general, since it includes many applications of interest in online optimization, for example when the variability comes from a changing distribution of the data.

In the literature, a widely used alternative model for $\sigma_i$ is to assume that it is a deterministic (albeit usually unknown) quantity that admits an upper bound. We remark that the proposed stochastic approach subsumes this deterministic model, since any (bounded) deterministic value is sub-Weibull with parameters $\theta = 0$ and $\nu = \sigma_i$. Therefore the results above specialize to the deterministic case in a straightforward fashion.
\end{remark}



\bibliographystyle{siamplain}
