\section{Perturbation Model}
\label{appendix_upperboundPr}
\begin{proposition}\label{prop_commutation_factor}
Let $\mathbf{T}_{r}, \mathbf{S}\in\text{End}(\mathcal{M})$ two compact normal operators acting on $\mathcal{M}$. Then, $\mathbf{S}\mathbf{T}_{r}=\mathbf{T}_{cr}\mathbf{S}+\mathbf{S}\mathbf{P}_{r}$ with $\mathbf{T}_{cr}=\sum_{i}\mu_{i}\mathbf{u}_{i}\langle\mathbf{u}_{i},\rangle$ and $\mathbf{P}_{r}\in\text{End}(\mathcal{M})$. Additionally, if $\delta = \max_{r} \frac{\| \bbP_{r}\|}{\| \bbT_r\|}$ then
\begin{equation}
 \delta \leq \max_r \sum_{i,j}\left\vert\frac{\lambda_{i}}{\lambda_{max}}\right\vert\left\Vert\bbT_{\mathbf{u}_{i}}\bbT_{\mathbf{v}_{j}}-\bbT_{\mathbf{u}_{i}}\bbT_{\mathbf{u}_{j}}\right\Vert,
\label{eq:deltaigenvectors0}
\end{equation}
where $\bbT_{\mathbf{u}_{i}}=\mathbf{u}_{i}\left\langle\mathbf{u}_{i},\cdot\right\rangle$,  $\bbT_{\mathbf{v}_{j}}=\mathbf{v}_{j}\langle\mathbf{v}_{j},\cdot\rangle$, $\mathbf{u}_{i}$ is the $i$th eigenvector of $\bbS$, $\mathbf{v}_{j}$ is the $j$th eigenvector of $\bbP_{r}$, $\langle\cdot,\cdot\rangle$ represents the inner product, $\lambda_{i}$ is the $i$th eigenvalue of $\bbS$, $\mu_{i}$ is the $i$th eigenvalue of $\mathbf{T}_{r}$ and $\lambda_{max}=\sup_{i}\{\vert\lambda_{i}\vert\}$.
\end{proposition}
\begin{proof}
Let $\mu_{i}$ the $i$th eigenvalue of $\mathbf{T}_{r}$. Taking into account that $\Vert\mathbf{P}_{r}\Vert$ can be computed as $\Vert\mathbf{P}_{r}\Vert=\sup_{\mathbf{S}}\frac{\Vert\mathbf{S}\mathbf{P}_{r}\Vert}{\Vert\mathbf{S}\Vert}$, we have that
\begin{multline}
\frac{\Vert\mathbf{S}\mathbf{P}_{r}\Vert}{\Vert\mathbf{S}\Vert}=
\frac{\left\Vert\mathbf{S}\mathbf{T}_{r}-\mathbf{T}_{cr}\mathbf{S}\right\Vert}{\lambda_{max}}
\\
=\frac{1}{\lambda_{max}}\left\Vert\sum_{i,j}\lambda_{i}\mu_{j}\mathbf{T}_{\mathbf{u}_{i}}\mathbf{T}_{\mathbf{v}_{j}}-\sum_{i,j}\lambda_{i}\mu_{j}\mathbf{T}_{\mathbf{u}_{i}}\mathbf{T}_{\mathbf{u}_{j}}\right\Vert
\end{multline}
and using triangle inequality it follows that:
\begin{multline}
\frac{\Vert\mathbf{S}\mathbf{P}_{r}\Vert}{\Vert\mathbf{S}\Vert}
\leq
\sum_{i,j}\left\vert\frac{\lambda_{i}}{\lambda_{max}}\right\vert\vert\mu_{j}\vert
\left\Vert\mathbf{T}_{\mathbf{u}_{i}}\mathbf{T}_{\mathbf{v}_{j}}-\mathbf{T}_{\mathbf{u}_{i}}\mathbf{T}_{\mathbf{u}_{j}}\right\Vert
\\
\leq
\Vert\mathbf{T}_{r}\Vert\sum_{i,j}\left\vert\frac{\lambda_{i}}{\lambda_{max}}\right\vert
\left\Vert\mathbf{T}_{\mathbf{u}_{i}}\mathbf{T}_{\mathbf{v}_{j}}-\mathbf{T}_{\mathbf{u}_{i}}\mathbf{T}_{\mathbf{u}_{j}}\right\Vert
\end{multline}
\end{proof}
Notice that the right hand side of eqn.~(\ref{eq:deltaigenvectors0}) provides a measure of the differences between the eigenvectors of $\mathbf{S}$ and $\mathbf{T}_{r}$ that is weighted by the relative size of the eigenvalues $\lambda_{i}$. As expected, when the changes in the eigenvectors are associated to eigenvectors whose eigenvalues are small, value of $\delta$ is small.
\section{Frechet Derivative $D_{p\vert\mathbf{S}_{i}}(\mathbf{S})$}
\label{appendix:howtofindDH}
First, notice that $p(\mathbf{S})=\sum_{k_{1},\ldots,k_{m}=0}^{\infty}h_{k_{1}\ldots k_{m}}\mathbf{S}_{1}^{k_{1}}\ldots\mathbf{S}_{m}^{k_{m}}
=
\sum_{k_{i}=0}^{\infty}\mathbf{S}_{i}^{k_{i}}\mathbf{A}_{k_{i}}$,
where $\mathbf{A}_{k_{i}}=\sum_{\substack{\{k_{j}\}=0\\ j\neq i}}^{\infty}h_{k_{1},\ldots,k_{m}}\prod_{\substack{j=1\\j\neq i}}^{m}\mathbf{S}_{j}^{k_{j}}.$
Then, it follows that
\begin{multline}
p(\mathbf{S}+\boldsymbol{\xi})-p(\mathbf{S})=
\sum_{k_{i}=0}^{\infty}\left(\mathbf{S}_{i}+\boldsymbol{\xi}_{i}\right)^{k_{i}}\mathbf{A}_{k_{i}}-
\sum_{k_{i}=0}^{\infty}\mathbf{S}_{i}^{k_{i}}\mathbf{A}_{k_{i}}
\label{eq:howtofingDH1}
\end{multline}
for $\boldsymbol{\xi}=(\mathbf{0},\ldots,\boldsymbol{\xi}_{i},\ldots,\mathbf{0})$. Considering the expansion $(\mathbf{S}_{i}+\boldsymbol{\xi}_{i})^{k_{i}}=\mathbf{S}^{k_{i}}_{i}+\boldsymbol{\xi}_{i}^{k_{i}}+\sum_{r=1}^{k-1}\boldsymbol{\pi}(r\mathbf{S}_{i},(k_{i}-r)\boldsymbol{\xi}_{i})$ for $k_{i}\geq 2$, eqn.~(\ref{eq:howtofingDH1}) takes the form
\begin{multline}
p(\mathbf{S}+\boldsymbol{\xi})-p(\mathbf{S})=
\\
\sum_{k_{i}=1}^{\infty}\sum_{r=1}^{k_{i}-1}\boldsymbol{\pi}\left(r\boldsymbol{\xi}_{i},(k_{i}-r)\mathbf{S}_{i}\right)\mathbf{A}_{k_{i}}+\sum_{k_{i}=1}^{\infty}\boldsymbol{\xi}_{i}^{k_{i}}\mathbf{A}_{k_{i}}.
\label{eq:howtofingDH2}
\end{multline}
Separating the linear terms on $\boldsymbol{\xi}_{i}$ eqn.~(\ref{eq:howtofingDH2}) leads to
\begin{multline}
p(\mathbf{S}+\boldsymbol{\xi})-p(\mathbf{S})=
\sum_{k_{i}=1}^{\infty}\boldsymbol{\pi}\left(\boldsymbol{\xi}_{i},(k_{i}-1)\mathbf{S}_{i}\right)\mathbf{A}_{k_{i}}
\\
+\sum_{k_{i}=2}^{\infty}\sum_{r=2}^{k_{i}-1}\boldsymbol{\pi}\left(r\boldsymbol{\xi}_{i},(k_{i}-r)\mathbf{S}_{i}\right)\mathbf{A}_{k_{i}}
+\sum_{k_{i}=2}^{\infty}\boldsymbol{\xi}^{k_{i}}\mathbf{A}_{k_{i}}.
\end{multline}
Therefore, taking into account the definition of Fr\'echet derivative (see Section~\ref{theorem:HvsFrechet}) it follows that
\begin{equation}
D_{p\vert\mathbf{S}_{i}}(\mathbf{S})\left\lbrace\boldsymbol{\xi}_{i}\right\rbrace=\sum_{k_{i}=1}^{\infty}\boldsymbol{\pi}\left(\boldsymbol{\xi}_{i},(k_{i}-1)\mathbf{S}_{i}\right)\mathbf{A}_{k_{i}}
\end{equation}

\subsubsection{#1}\vspace{-3\baselineskip}\color{black}\medskip{\noindent \bf \thesubsubsection. #1.}}

\newcommand{\myparagraph}[1]{\needspace{1\baselineskip}\medskip\noindent {\bf #1}}

\newcommand{\myindentedparagraph}[1]{\needspace{1\baselineskip}\medskip \hangindent=11pt \hangafter=0 \noindent{\it #1.}}

\newcommand{\myparagraphtc}[1]{\needspace{1\baselineskip}\medskip\noindent {\it #1.}\addcontentsline{toc}{subsubsection}{\qquad\qquad\quad#1}}

\section{Perturbations}\label{sec:perturbandstability}
In an ASP triple $(\ccalA,\ccalM,\rho)$, signals $\bbx\in\ccalM$ are observations of interest and the algebra $\ccalA$ defines the operations that are to be performed on signals. The homomorphism $\rho$ ties these two objects and, as such, is one we can consider as subject to model mismatch. In this paper we consider perturbations adhering to the following model.

\begin{definition}\label{def:perturbmodel}(ASP Model Perturbation)
Let $(\ccalA,\ccalM,\rho)$ be an ASP model with algebra elements generated by $g\in\ccalG$ (Definition \ref{def_generators}) and recall the definition of the shift operators $\bbS = \rho(g)$ (Definition \ref{def_shift_operators}). We say that $(\ccalA,\ccalM,\tdrho)$ is a perturbed ASP model if for all $a=p_{\ccalA}(\ccalG)$ we have that
\begin{equation}\label{eqn_def_perturbation_model_10}
   \tdrho(a) = p_\ccalM\big(\tdrho(g)\big)
             = p_\ccalM\big(\tilde\ccalS\big)
             =p\big(\tilde\ccalS\big),
\end{equation}
where $\tilde\ccalS$ is a set of perturbed shift operators of the form
\begin{equation}\label{eqn_def_perturbation_model_20}
   \tbS = \bbS + \bbT(\bbS),
\end{equation}
for all shift operators $\bbS\in\ccalS$ \end{definition}

As per Definition \ref{def:perturbmodel}, an ASP perturbation model, is a perturbation of the homomorphism $\rho$ defined by a perturbation of the shift operators $\bbS$. Each shift operator $\bbS$ is perturbed to the shift operator $\tbS$ according to \eqref{eqn_def_perturbation_model_20} and this perturbation propagates to the filter $\rho(a)$ according to \eqref{eqn_def_perturbation_model_10}. An important technical remark is that the resulting mapping $\tdrho$ is not required to be a homomorphism -- although it can be, indeed, often is.

We point out that Definition \ref{def:perturbmodel} limits the perturbation of the homomorphism $\rho$ to perturbations of the shift operators. This is a somewhat arbitrary choice but one that is motivated by practice as we illustrate for graph signals (Example \ref{ex_GSP_perturbation}), discrete time signals (Example \ref{ex_DTSP_perturbation}), signals on groups (Example \ref{ex_groupSP_perturbation}) and graphon signals (Example \ref{ex_WSP_perturbation}). The model is also arbitrary in that it perturbs the homomorphism $\rho$ but not the Algebra $\ccalA$ or the signal $\bbx$. This is also justifiable. Notice, first that the algebra $\ccalA$ is a choice of operations and therefore not naturally subject to perturbation. Perturbing signals $\bbx$ is of interest but a transformation of $\bbx$ can be reinterpreted as a transformation of $\bbS$. Indeed, if we are given signals $\bbx$ and $T\bbx$ we can always define a shift operator $\tilde{\bbS}=\bbS T$ and write
\begin{equation}\label{eqn_def_perturbation_from_x_to_S}
   \bbS (T\bbx) = (\bbS T) \bbx =  \tilde{\bbS}\bbx.
\end{equation}
Thus, we understand the effect of the perturbation $T\bbx$ by studying the processing of $\bbx$ with the perturbed shift operator $\tilde{\mathbf{S}}=\bbS T$. We discuss this further for
discrete time signals in Example \ref{ex_DTSP_perturbation} and for signals on groups in Example \ref{ex_groupSP_perturbation} in the following. We also discuss the model in \eqref{def:perturbmodel} particularized to graphs (Example \ref{ex_GSP_perturbation}) and graphons (Example \ref{ex_WSP_perturbation}).

\subsection{Perturbation Models}
In our subsequent analysis we consider perturbation models of the form
\begin{equation}\label{eqn_perturbation_model_absolute_plus_relative}
   \bbT(\bbS)=\bbT_{0} + \bbT_{1}\bbS,
\end{equation}
which is a generic model of small perturbations of a shift operator that involve an absolute perturbation $\bbT_{0}$ and a relative perturbation $\bbT_{1}\bbS$;  see \cite{fern2019stability}.  The $\bbT_{r}$ are compact normal operators with operator norm $\Vert\bbT_{r}\Vert<1$. Requiring $\Vert\bbT_{r}\Vert<1$ is a minor restriction as we are interested in small perturbations with $\Vert\bbT_{r}\Vert\ll1$.

For the model in \eqref{eqn_perturbation_model_absolute_plus_relative} it is important to describe the commutativity of the shift operator $\bbS$ and the perturbation model matrices $\bbT_r$. To that end, we write
\begin{equation}\label{eq:PrvsTr}
\mathbf{S}\mathbf{T}_{r}=\mathbf{T}_{cr}\mathbf{S}+\mathbf{S}\mathbf{P}_{r},
\end{equation}
where $\mathbf{T}_{cr}=\sum_{i}\mu_{i}\mathbf{u}_{i}\langle\mathbf{u}_{i},\cdot\rangle$, $\mu_{i}$ is the $i$th eigenvalue of $\mathbf{T}_{r}$, $\mathbf{u}_{i}$ is the $i$th eigenvector of $\mathbf{S}$, and $\langle,\rangle$ represents the inner product operation. As a consequence, we have that $\mathbf{S}\mathbf{T}_{rc}=\mathbf{T}_{rc}\mathbf{S}$ and $\Vert\mathbf{T}_{cr}\Vert=\Vert\mathbf{T}_{r}\Vert$. We define the commutation factor
\begin{equation}\label{eqn_commutation_factor}
   \delta = \max_{r} \frac{\| \bbP_{r}\|}{\| \bbT_r\|},
\end{equation}
which is a measure of how far the operators $\bbS$ and $\bbT_r$ are from commuting with each other. Notice that $\delta=0$ implies $\mathbf{P}_{r}=\mathbf{0}$ and $\mathbf{T}_{r}=\mathbf{T}_{cr}$. The commutation factor $\delta$ in \eqref{eqn_commutation_factor} can be bounded in terms of the spectra of $\bbT_r$ and $\bbS$ as we show in Proposition \ref{prop_commutation_factor}. The specifics of this bound are not central to the results of Section \ref{sec:stabilitytheorems}.

\begin{remark}\normalfont
Notice that $\delta$ is a property of the perturbation and not from the filters (see remark~\ref{remark:interpSpecFilt}). As we will show in the next section a large value of $\delta$ does not modify the property of being stable, although it increases the values of the constants associated to stability. Studying strategies to reduce the effects of $\delta$ in the stability constants is an interesting future research direction.
\end{remark}

Notice that when representations of an algebra $\mathcal{A}$ with multiple generators $\{g_{i}\}_{i=1}^{m}$ are considered, we have that for $a\in\mathcal{A}$ the operator
$p(\rho(a))\in\text{End}(\mathcal{M})$ is a function of  $\rho(g_{i})=\mathbf{S}_{i}\in\text{End}(\mathcal{M})$ and therefore can be seen as the function $p: \text{End}(\mathcal{M})^{m}\rightarrow\text{End}(\mathcal{M})$, where $\text{End}(\mathcal{M})^{m}$ is the $m$-times cartesian product of $\text{End}(\mathcal{M})$.  In this scenario we use the notation $p(\mathbf{S})=p(\mathbf{S}_{1},\ldots,\mathbf{S}_{m})$  and  when considering the perturbation model in eqn.~(\ref{eqn_perturbation_model_absolute_plus_relative}) acting on $\mathbf{S}=(\mathbf{S}_{1},\ldots,\mathbf{S}_{m})$ we use the following notation $\mathbf{T}(\mathbf{S})=\left(\mathbf{T}(\mathbf{S}_{1}),\ldots,\mathbf{T}(\mathbf{S}_{m})\right)$ where $\mathbf{T}(\mathbf{S}_{i})=\mathbf{T}_{0,i}+\mathbf{T}_{1,i}\mathbf{S}_{i}$.
\section{Stability Theorems}\label{sec:stabilitytheorems}
The filters in Section \ref{sec_alg_filters} and the algebraic neural networks in Section \ref{sec_Algebraic_NNs} are operators acting on the space $\ccalM$. These operators are of the form $p(\mathbf{S})$, and their outputs depend on a filter set $\ccalH\subset\mathcal{A}$ which is denoted as $p_{\mathcal{A}}\in\mathcal{H}\subset\mathcal{A}$, and the set of shift operators $\ccalS$, where $\mathbf{S}\in\mathcal{S}$. When we perturb the processing model according to Definition \ref{def:perturbmodel}, these operators are perturbed as well. The goal of this paper is to analyze these perturbations. In particular, our goal is to identify conditions for filters and algebraic neural networks to be stable in the sense of the following definition.

\begin{definition}[Operator Stability]\label{def:stabilityoperators1} Given operators $p(\mathbf{S})$ and $p(\tilde{\mathbf{S}})$ defined on the processing models $(\ccalA,\ccalM,\rho)$ and $(\ccalA,\ccalM,\tdrho)$ (cf. Definition \ref{def:perturbmodel}) we say the operator $p(\mathbf{S})$ is Lipschitz stable if there exist constants $C_{0}, C_{1}>0$ such that
\begin{multline}\label{eq:stabilityoperators1}
\left\Vert p(\mathbf{S})\mathbf{x}  - p(\tilde{\mathbf{S}})\mathbf{x}\right\Vert
\leq
\\
\left[
C_{0} \sup_{\bbS\in\ccalS}\Vert\mathbf{T}(\mathbf{S})\Vert + C_{1}\sup_{\bbS\in\ccalS}\big\|D_{\bbT}(\bbS)\big\|
+\mathcal{O}\left(\Vert\mathbf{T}(\mathbf{S})\Vert^{2}\right)
\right] \big\| \bbx \big\|,
\end{multline}
for all $\bbx\in\ccalM$. In \eqref{eq:stabilityoperators1} $D_{\bbT}(\bbS)$ is the Fr\'echet derivative of the perturbation operator $\bbT$.
\end{definition}

When the perturbation value $\bbT(\bbS)$ and its derivative $D_{\bbT}(\bbS)$ are small, the inequality in \eqref{eq:stabilityoperators1} states that the operators $p(\mathbf{S})$ and $ p(\tilde{\mathbf{S}})$ are close uniformly across all inputs $\bbx$. Our stability theorems are presented in the next section, but at this point it is important to remark that algebraic filters are not always stable in the sense of \eqref{eq:stabilityoperators1}. We know that this is true because unstable counterexamples are known in the case of graph signal processing \cite{gamagnns} and the processing of time signals \cite{mallat_ginvscatt}.

\begin{example}[Graph Signal Processing]\label{ex_GSP_perturbation}\normalfont In graph signal processing the shift operator $\bbS$ is a matrix representation of a graph. The perturbation model in \eqref{def:perturbmodel} simply states that $\tbS$ is a matrix representation of a different graph. Definition \ref{def:stabilityoperators1} defines a stable operator $p(\mathbf{S})$ as one that doesn't change much when run on graphs that are close and related by perturbations that are sufficiently smooth in the space of matrix representations of graphs. The notion of absolute graph perturbation considered for the stability analysis of GNNs in~\cite{fern2019stability} is a particular case of eqn.~(\ref{eqn_perturbation_model_absolute_plus_relative}) with $\mathbf{T}_{1}=\mathbf{0}$. Additionally, the notion of relative graph perturbations can be obtained from eqn.~(\ref{eqn_perturbation_model_absolute_plus_relative}) considering $\mathbf{T}_{0}=\mathbf{0}$.
\end{example}

\begin{example}[Discrete Time Signal Processing]\label{ex_DTSP_perturbation} \normalfont
In discrete time signal processing the shift operator is the discrete time shift $S$. In this case the perturbation model is more meaningful if interpreted as a perturbation of the signal $\bbx$. The processing induced by \eqref{eqn_DTSP_filter_outputs} is invariant to shifts and therefore adequate to processing signals that are shift invariant. In general, signals are close to shift invariant but not exactly so. That is, a given signal $\bbx$ is invariant with respect to a shift operator $\tdS$ that is close to the time shift $S$. If the stability property in \eqref{eq:stabilityoperators1} holds we can guarantee that processing the signal $\bbx$ with the operator $\tdS$ is not far from processing the signal with the operator $S$. The latter represents the operations we perform -- since we choose to use $S$ in the processing of time signals. The former represents the processing we should undertake to respect the actual invariance properties of the signal $\bbx$ -- which are characterized by $\tdS$.
\end{example}

\begin{example}[Graphon Signal Processing]\label{ex_WSP_perturbation} \normalfont
Similar to Example \ref{ex_GSP_perturbation} the graphon $W(u,v)$ is a limit object that represents a family of random graphs. The perturbed graphon $\tdW(u,v)$ represents a different family of random graphs. The perturbation of the graphon generates a corresponding perturbation of the shift operator defined in \eqref{eqn_graphon_shift}. If the condition stated in \eqref{eq:stabilityoperators1} is satisfied, then filtering graphon signals using the perturbed shift operator associated to $\tdW(u,v)$  will lead to similar results to the ones obtained with the unperturbed operator and the differences are proportional to the size of the perturbation acting on $W(u,v)$.
\end{example}

\begin{example}[Group Signal Processing]\label{ex_groupSP_perturbation} \normalfont
Similar to Example \ref{ex_DTSP_perturbation}, the filters in \eqref{eqn_group_filters} are invariant to the action of the group. Actual signals $\bbx$ are invariant to actions of operators that are close to actions of the group -- e.g., a signal is close to invariant to rotations and symmetries. If \eqref{eq:stabilityoperators1} is true, processing the signal with operators $g$ -- as we choose to do -- is not far from processing the signal with operators $\tdg$ -- as we should do to leverage the actual invariance of the signal $\bbx$. We remark that when we perturb $g$ to $\tdg$ the resulting shift operators will not, in general, be representations of a homomorphism.
\end{example}

\subsection{Stability of Algebraic Filters}
Taking into account that the notion of stability is meant to be satisfied by subsets of filters of the algebra and not necessarily the whole algebra, it is important to have a characterization of these subsets in simple terms. In the case of commutative algebras with $m$ generators this characterization can be stated using $m$-multivariate functions. To do so, we use an algebra $\mathcal{A}_{\mathbb{C}}$ isomorphic to $\mathcal{A}$, where the symbols are variables taking values in $\mathbb{C}$ and then define conditions in this new algebra.  For instance,  if $\mathcal{A}$ is the algebra generated by $g$ with coefficients in $\mathbb{R}$, we have the isomorphism $\iota:\mathcal{A}\rightarrow\mathcal{A}_{\mathbb{C}}$ given by $\iota(g)=\lambda$ where $\lambda\in\mathbb{C}$. Then, $\iota\left(p_{\mathcal{A}}(g)\right)=p(\lambda)$,  and stating conditions on $p(\lambda)$ we can define subsets of $\mathcal{A}$ via $\iota^{-1}$. In the following paragraphs we introduce some definitions used to characterize subsets of algebras with one generator. Results associated to multiple generators will be discussed in the next subsection.

\begin{definition}
Let $p:\mathbb{C}\rightarrow\mathbb{C}$ be  single variable function. Then, it is said that $p$ is Lipschitz if there exists $L_{0}>0$ such that
\begin{equation}
\vert p(\lambda)-p(\mu)\vert\leq L_{0}\vert\lambda-\mu\vert
\end{equation}
for all $\lambda, \mu\in\mathbb{C}$. Additionally, it is said that $p(\lambda)$ is Lipschitz integral if there exists $L_{1}>0$ such that
\begin{equation}
\left\vert
\lambda\frac{dp(\lambda)}{d\lambda}
\right\vert
\leq L_{1}
~\forall~\lambda.
\end{equation}
\end{definition}

In what follows, when considering subsets of a commutative algebra $\mathcal{A}$, we denote by $\mathcal{A}_{L_{0}}$ the subset of elements in $\mathcal{A}$ that are Lipschitz with constant $L_{0}$ and by $\mathcal{A}_{L_{1}}$ the subset of element of $\mathcal{A}$ that are Lipschitz integral with constant $L_{1}$. Additionally, for the sake of simplicity we will not make reference to $\iota$.

We start with a result for operators in algebraic models with a single generator, highlighting the role of the Fr\'echet derivative of the map that relates the operator and its perturbed version.
\begin{theorem}\label{theorem:HvsFrechet}
Let $\mathcal{A}$ be an algebra generated by $g$ and let $(\mathcal{M},\rho)$ be a representation of $\mathcal{A}$ with $\rho(g)=\bbS\in\text{End}(\mathcal{M})$. Let $\tilde{\rho}(g)=\tilde{\bbS}\in\text{End}(\mathcal{M})$ where the pair
$(\mathcal{M},\tilde{\rho})$ is a perturbed version of $(\mathcal{M},\rho)$ and $\tilde{\bbS}$ is related to $\bbS$ by the perturbation model in eqn.~(\ref{eqn_def_perturbation_model_20}). Then, for any $p_{\mathcal{A}}\in\mathcal{A}$ we have
\begin{equation}
\left\Vert p(\bbS)\mathbf{x}-p(\tilde{\bbS})\mathbf{x}\right\Vert
\leq
\Vert\mathbf{x}\Vert
\left(
\left\Vert D_{p}(\mathbf{S})\left\lbrace\mathbf{T}(\mathbf{S})\right\rbrace\right\Vert + \mathcal{O}\left(\Vert\mathbf{T}(\mathbf{S})\Vert^{2}
\right)\right)
\label{eq:HSoptbound}
\end{equation}
where $D_{p}(\mathbf{S})$ is the Fr\'echet derivative of $p$ on $\mathbf{S}$.
\end{theorem}
\begin{proof}
See Section~\ref{prooftheoremHvsFrechet}
\end{proof}
Theorem~\ref{theorem:HvsFrechet} highlights an important point, the difference between two operators obtained from the same elements in the algebra is bounded by the Fr\'echet derivative of $p(\bbS)$ which depends of the properties of the elements in $\mathcal{A}$. In particular, we can see that an upper bound in the term $\left\Vert D_{p}(\mathbf{S})\bbT(\bbS)\right\Vert$ depends on how the the operator $D_{p}(\mathbf{S})$ acts on the perturbation $\bbT(\bbS)$. Then, $D_{p}(\mathbf{S})$ will determine whether $p(\bbS)$ is stable under the effect of $\bbT(\bbS)$, or in other words the properties of $p_{\mathcal{A}}$ act on the perturbation via the operator $D_{p}(\mathbf{S})$. Additionally, notice that eqn.~(\ref{eq:HSoptbound}) is satisfied for any $\bbT(\bbS)$ if $D_{p}(\mathbf{S})$ exists.

In the following theorems we show how these terms are related to $\bbT(\bbS)$ and its Fr\'echet derivative $D_{\bbT}$.

\begin{theorem}\label{theorem:uppboundDH}
Let $\mathcal{A}$ be an algebra with one generator element $g$ and let $(\mathcal{M},\rho)$ be a finite or countable infinite dimensional representation of $\mathcal{A}$. Let  $(\mathcal{M},\tilde{\rho})$ be a perturbed version of $(\mathcal{M},\rho)$ associated to the perturbation model in eqn.~(\ref{eqn_perturbation_model_absolute_plus_relative}).  If $p_{\mathcal{A}}\in\mathcal{A}_{L_{0}}\cap\mathcal{A}_{L_{1}}$, then
\begin{equation}
\left\Vert D_{p}\bbT(\bbS)\right\Vert
\leq(1+\delta)
\left(L_{0}\sup_{\bbS}\Vert\bbT(\bbS)\Vert +L_{1}\sup_{\bbS}\Vert D_{\bbT}(\bbS)\Vert\right)
\label{eq:DHTS}
\end{equation}
 \end{theorem}
\begin{proof}See Section~\ref{prooftheorem:uppboundDH}\end{proof}

It is worth pointing out that the constants involved in the upper bound of eqn.~(\ref{eq:DHTS}) depend on the properties of the filters and the difference between the eigenvectors of $\bbS$ and $\bbT_{r}$. Therefore, the difference between the eigenvectors of these operators do not determine if $p(\bbS)$ is stable or not, although the absolute value of the stability constants increase proportionally to $\delta$.

From theorems~\ref{theorem:HvsFrechet} and~\ref{theorem:uppboundDH} we can state the notion of stability for algebraic filters in the following corollary.

\begin{corollary}\label{corollary:DHboundsprcases}
Let $\mathcal{A}$ be an algebra with one generator element $g$ and let $(\mathcal{M},\rho)$ be a finite or countable infinite dimensional representation of $\mathcal{A}$. Let  $(\mathcal{M},\tilde{\rho})$ be a perturbed version of $(\mathcal{M},\rho)$ related by the perturbation model in eqn.~(\ref{eqn_perturbation_model_absolute_plus_relative}). Then, if $p_{\mathcal{A}}\in\mathcal{A}_{L_{0}}\cap\mathcal{A}_{L_{1}}$ the operator $p(\mathbf{S})$ is stable in the sense of definition~\ref{def:stabilityoperators1} with $C_{0}=(1+\delta)L_{0}$ and $C_{1}=(1+\delta)L_{1}$.
\end{corollary}
\begin{proof}
Replace (\ref{eq:DHTS}) from Theorem \ref{theorem:uppboundDH} into (\ref{eq:HSoptbound}) from Theorem \ref{theorem:HvsFrechet} and reorder terms
\end{proof}

\subsection{Stability of Algebraic Filters in Algebras with multiple geneators}
The stability results presented in previous subsection can be extended naturally to operators associated to representations of algebras with multiple generators. Before doing so, we extend the notions used to characterize subsets of algebras. If $\mathcal{A}$ is the algebra generated by $\{g_{i}\}_{i=1}^{m}$ and $g_{i}g_{j}=g_{j}g_{i}$ for all $i,j$ with coefficients in $\mathbb{R}$, we use the isomorphism $\iota:\mathcal{A}\rightarrow\mathcal{A}_{\mathbb{C}}$ given by $\iota(g_{i})=\lambda_{i}\in\mathbb{C}$. Then, we can use $\iota\left(p_{\mathcal{A}}(g_{1},\ldots,g_{m})\right)=p(\boldsymbol{\lambda})$ where $\boldsymbol{\lambda}=(\lambda_{1},\ldots,\lambda_{m})\in\mathbb{C}^{m}$ to characterize elements in $\mathcal{A}$. Therefore, stating conditions on $p(\boldsymbol{\lambda})$ we can define subsets of $\mathcal{A}$ via $\iota^{-1}$. In the following definition we extend the concepts presented in previous subsection.
\begin{definition}
Let $p:\mathbb{C}^{m}\rightarrow\mathbb{C}$ be  a multivariate function. Then, it is said that $p$ is Lipschitz if there exists $L_{0}>0$ such that
\begin{equation}
\vert p(\boldsymbol{\lam})-p(\boldsymbol{\mu})\vert\leq L_{0}\Vert\boldsymbol{\lam}-\boldsymbol{\mu}\Vert
\end{equation}
for all $\boldsymbol{\lambda},\boldsymbol{\mu}\in\mathbb{C}^{m}$. Additionally, it is said that $p(\boldsymbol{\lambda})$ is Lipschitz integral if there exists $L_{1}>0$ such that
\begin{equation}
\left\vert \lambda_{i}\frac{\partial p(\boldsymbol{\lambda})}{\partial \lambda_{i}}\right\vert\leq L_{1}\quad\forall~i\in\{1,\ldots m \},
\end{equation}
where $\boldsymbol{\lambda}=(\lambda_{1},\ldots,\lambda_{m})$.
\end{definition}
With this notion at hand, we are ready to  extend the stability theorems.
\begin{theorem}\label{theorem:HvsFrechetmultgen}
Let $\mathcal{A}$ be an algebra generated by $\{g_{i}\}_{i=1}^{m}$ and let $(\mathcal{M},\rho)$ be a representation of $\mathcal{A}$ with $\rho(g_{i})=\bbS_{i}\in\text{End}(\mathcal{M})$ for all $i$. Let $\tilde{\rho}(g_{i})=\tilde{\bbS}_{i}\in\text{End}(\mathcal{M})$ where the pair
$(\mathcal{M},\tilde{\rho})$ is a perturbed version of $(\mathcal{M},\rho)$ and $\tilde{\bbS}_{i}$ is related to $\bbS_{i}$ by the perturbation model in eqn.~(\ref{eqn_def_perturbation_model_20}). Then, for any $p_{\mathcal{A}}\in\mathcal{A}$ we have
\begin{equation}
\left\Vert p(\bbS)\mathbf{x}-p(\tilde{\bbS})\mathbf{x}\right\Vert\leq
\Vert\mathbf{x}\Vert\sum_{i=1}^{m}\left(\left\Vert D_{p\vert\mathbf{S}_{i}}(\mathbf{S})\mathbf{T}(\mathbf{S}_{i})\right\Vert+\mathcal{O}\left(\Vert\mathbf{T}(\mathbf{S}_{i})\Vert^{2}\right)\right)
\label{eq:HSoptboundmultgen}
\end{equation}
where $D_{p\vert\mathbf{S}_{i}}(\mathbf{S})$ is the partial Fr\'echet derivative of $p$ on $\mathbf{S}_{i}$.
\end{theorem}
\begin{proof}
See Section~\ref{prooftheoremHvsFrechet}
\end{proof}
Notice that in eqn.~(\ref{eq:HSoptboundmultgen}) we naturally add the contribution associated to each generator. Therefore, to guarantee stability we must have stability in each generator. Now, we show how the Fr\'echet derivative of $\mathbf{T}(\mathbf{S})$ is involved in the stability properties when considering multiple generators.

\begin{theorem}\label{theorem:uppboundDHmultg}
Let $\mathcal{A}$ be an algebra with $m$ generators $\{g_{i}\}_{i=1}^{m}$ and $g_{i}g_{j}=g_{j}g_{i}$ for all $i,j\in\{1,\ldots m\}$. Let $(\mathcal{M},\rho)$ be a finite or countable infinite dimensional representation of $\mathcal{A}$ and $(\mathcal{M},\tilde{\rho})$ a perturbed version of $(\mathcal{M},\rho)$ related by the perturbation model in eqn.~(\ref{eqn_perturbation_model_absolute_plus_relative}).  Then, if $p_{\mathcal{A}}\in\mathcal{A}_{L_{0}}\cap\mathcal{A}_{L_{1}}$ it holds that
\begin{multline}
\left\Vert D_{p\vert\mathbf{S}_{i}}(\mathbf{S})\bbT(\bbS_{i})\right\Vert
\\
\leq (1+\delta)
\left(
L_{0}\sup_{\bbS_{i}\in\mathcal{S}}\Vert\bbT(\bbS_{i})\Vert +L_{1}\sup_{\bbS_{i}\in\mathcal{S}}\Vert D_{\bbT}(\bbS_{i})\Vert\right)
\label{eq:DHTSmultg}
\end{multline}
 \end{theorem}
 \begin{proof}See Section~\ref{prooftheorem:uppboundDH}\end{proof}
 %
 It is important to remark that the upper bound in eqn.~(\ref{eq:DHTSmultg}) is defined by the largest perturbation in a given generator although the constants associated are determined completely by the properties of the filters.

From theorems~\ref{theorem:HvsFrechetmultgen} and~\ref{theorem:uppboundDHmultg} we can state the stability results for filters in algebras with multiple generators in the following corollary.
\begin{corollary}\label{corollary:stabmultgenfilt}
Let $\mathcal{A}$ be an algebra with generators $\{ g_{i}\}_{i=1}^{m}$ and $g_{i}g_{j}=g_{j}g_{i}$ for all $i,j$. Let $(\mathcal{M},\rho)$ be a finite or countable infinite dimensional representation of $\mathcal{A}$ and $(\mathcal{M},\tilde{\rho})$ be a perturbed version of $(\mathcal{M},\rho)$ related by the perturbation model in eqn.~(\ref{eqn_perturbation_model_absolute_plus_relative}). Then, if $p_{\mathcal{A}}\in\mathcal{A}_{L_{0}}\cup\mathcal{A}_{L_{1}}$ the operator $p(\mathbf{S})$ is stable in the sense of definition~\ref{def:stabilityoperators1} with $C_{0}=m(1+\delta)L_{0}$ and $C_{1}=m(1+\delta)L_{1}$.
\end{corollary}
\begin{proof} Replacing eqn.~(\ref{eq:DHTSmultg}) from theorem~\ref{theorem:uppboundDHmultg} into eqn.~(\ref{eq:HSoptboundmultgen}) from theorem~\ref{theorem:HvsFrechetmultgen} and organazing the terms.\end{proof}

\subsection{Stability of Algebraic Neural Networks}
The results in Theorems~\ref{theorem:HvsFrechet} to~\ref{theorem:uppboundDHmultg} and corollaries~\ref{corollary:DHboundsprcases} and~\ref{corollary:stabmultgenfilt} can be extended to operators representing AlgNNs. We say that for a given AlgNN,  $\Xi=\left\lbrace (\mathcal{A}_{\ell},\mathcal{M}_{\ell},\rho_{\ell})\right\rbrace_{\ell=1}^{L}$, a perturbed version of $\Xi$ is given by  $\tilde{\Xi}=\left\lbrace (\mathcal{A}_{\ell},\mathcal{M}_{\ell},\tilde{\rho}_{\ell})\right\rbrace_{\ell=1}^{L}$ where $(\mathcal{A}_{\ell},\mathcal{M}_{\ell},\tilde{\rho}_{\ell})$ is a perturbed version of $(\mathcal{A}_{\ell},\mathcal{M}_{\ell},\rho_{\ell})$. For the sake of simplicity we present a theorem for AlgNNs with algebras with a single generator, but notice that these results can be easily stated for
AlgNNs with multiple generators directly from theorems~\ref{theorem:HvsFrechetmultgen} and~\ref{theorem:uppboundDHmultg}. To do so, we start highlighting in the following theorem the stability properties of the operators in the layer $\ell$ of an AlgNN.

\begin{theorem}\label{theorem:stabilityAlgNN0}
Let $\Xi=\left\lbrace (\mathcal{A}_{\ell},\mathcal{M}_{\ell},\rho_{\ell})\right\rbrace_{\ell=1}^{L}$ be an algebraic neural network  with $L$ layers, one feature per layer and algebras $\mathcal{A}_{\ell}$ with a single generator.  Let  $\tilde{\Xi}=\left\lbrace (\mathcal{A}_{\ell},\mathcal{M}_{\ell},\tilde{\rho}_{\ell})\right\rbrace_{\ell=1}^{L}$ be the perturbed version of $\Xi$ by means of the perturbation model in eqn.~(\ref{eqn_perturbation_model_absolute_plus_relative}). Then, if  $\Phi\left(\mathbf{x}_{\ell-1}, \mathcal{P}_{\ell},\mathcal{S}_{\ell}\right)$ and
$\Phi\left(\mathbf{x}_{\ell-1},\mathcal{P}_{\ell},\tilde{\mathcal{S}}_{\ell}\right)$ represent the mapping operators associated to $\Xi$ and $\tilde{\Xi}$ in the layer $\ell$ respectively, we have
\begin{multline}
\left\Vert
\Phi\left(\mathbf{x}_{\ell-1},\mathcal{P}_{\ell},\mathcal{S}_{\ell}\right)-
\Phi\left(\mathbf{x}_{\ell-1},\mathcal{P}_{\ell},\tilde{\mathcal{S}}_{\ell}\right)
\right\Vert
\leq
\\
C_{\ell}(1+\delta_{\ell})\left(L_{0}^{(\ell)} \sup_{\bbS_{\ell}}\Vert\bbT^{(\ell)}(\bbS_{\ell})\Vert
+L_{1}^{(\ell)}\sup_{\bbS_{\ell}}\Vert D_{\mathbf{T^{(\ell)}}}(\bbS_{\ell})\Vert\right)\Vert\mathbf{x}_{\ell-1}\Vert
\label{eq:theoremstabilityAlgNN0}
\end{multline}
where $C_{\ell}$ is the Lipschitz constant of $\sigma_{\ell}$, and $\mathcal{P}_{\ell}=\mathcal{A}_{L_{0}}\cap\mathcal{A}_{L_{1}}$ represents the domain of $\rho_{\ell}$.
The index $(\ell)$ makes reference to quantities and constants associated to the layer $\ell$.
\end{theorem}
\begin{proof}See Section~\ref{prooftheorem:stabilityAlgNN0}\end{proof}
This result, although simple, highlights the role of the maps $\sigma_{\ell}$ when perturbations are considered in each layer. In particular, we see that the effect of $\sigma_{\ell}$  is to scale $\boldsymbol{\Delta}_{\ell}$ by a constant but it does not change the nature or mathematical form of the perturbation. Notice also that $\sigma_{\ell}$ plays the role of a mixer that allows an AlgNN to provide selectivity without affecting the stability (see Section~\ref{sec:discussion}).

Now we present in the following theorem the stability result for a general AlgNN with commutative algebras.

\begin{theorem}\label{theorem:stabilityAlgNN1}
Let $\Xi=\left\lbrace (\mathcal{A}_{\ell},\mathcal{M}_{\ell},\rho_{\ell})\right\rbrace_{\ell=1}^{L}$ be an algebraic neural network  with $L$ layers, one feature per layer and algebras $\mathcal{A}_{\ell}$ with a single generator.  Let  $\tilde{\Xi}=\left\lbrace (\mathcal{A}_{\ell},\mathcal{M}_{\ell},\tilde{\rho}_{\ell})\right\rbrace_{\ell=1}^{L}$ be the perturbed version of $\Xi$ by means of the perturbation model in eqn.~(\ref{eqn_perturbation_model_absolute_plus_relative}). Then, if  $\Phi\left(\mathbf{x},\{ \mathcal{P}_{\ell} \}_{1}^{L},\{ \mathcal{S}_{\ell}\}_{1}^{L}\right)$ and
$\Phi\left(\mathbf{x},\{ \mathcal{P}_{\ell} \}_{1}^{L},\{ \tilde{\mathcal{S}}_{\ell}\}_{1}^{L}\right)$ represent the mapping operators associated to $\Xi$ and $\tilde{\Xi}$ respectively, we have
\begin{multline}
\left\Vert
\Phi\left(\mathbf{x},\{ \mathcal{P}_{\ell} \}_{1}^{L},\{ \mathcal{S}_{\ell}\}_{1}^{L}\right)-
\Phi\left(\mathbf{x},\{ \mathcal{P}_{\ell} \}_{1}^{L},\{ \tilde{\mathcal{S}}_{\ell}\}_{1}^{L}\right)
\right\Vert
\\
\leq
\sum_{\ell=1}^{L}\boldsymbol{\Delta}_{\ell}\left(\prod_{r=\ell}^{L}C_{r}\right)\left(\prod_{r=\ell+1}^{L}B_{r}\right)
\left(\prod_{r=1}^{\ell-1}C_{r}B_{r}\right)\left\Vert\mathbf{x}\right\Vert
\label{eq:theoremstabilityAlgNN1}
\end{multline}
where $C_{\ell}$ is the Lipschitz constant of $\sigma_{\ell}$, $\Vert\rho_{\ell}(\xi)\Vert\leq B_{\ell}~\forall~\xi\in\mathcal{P}_{\ell}$, and $\mathcal{P}_{\ell}=\mathcal{A}_{L_{0}}\cap\mathcal{A}_{L_{1}}$ represents the domain of $\rho_{\ell}$. The functions $\boldsymbol{\Delta}_{\ell}$ are given by
\begin{equation}\label{eq:varepsilonl}
\boldsymbol{\Delta}_{\ell}=(1+\delta_{\ell})\left(L_{0}^{(\ell)} \sup_{\bbS_{\ell}}\Vert\bbT^{(\ell)}(\bbS_{\ell})\Vert
+L_{1}^{(\ell)}\sup_{\bbS_{\ell}}\Vert D_{\mathbf{T^{(\ell)}}}(\bbS_{\ell})\Vert\right)
\end{equation}
with the index $(\ell)$ indicating quantities and constants associated to the layer $\ell$.
\end{theorem}
\begin{proof}See Section~\ref{prooftheorem:stabilityAlgNN1}\end{proof}
Theorem~\ref{theorem:stabilityAlgNN1} states how an AlgNN can be made stable by the selection of an appropriate subset of filters in the algebra, for a given perturbation model. It is worth pointing out that conditions like the ones obtained in~\cite{fern2019stability} for GNNs can be considered particular instantiations of the conditions in Theorem~\ref{theorem:stabilityAlgNN1}. Additionally, notice that Theorem~\ref{theorem:stabilityAlgNN1} can be easily extended to consider several features per layer, the reader can check the details of the proof of the theorem in Section~\ref{prooftheorem:stabilityAlgNN1} where the analysis is performed considering multiple features.
\begin{remark}\normalfont
It is important to highlight the fact that the perturbation model in eqn.~(\ref{eqn_perturbation_model_absolute_plus_relative}) is \textit{smooth} in the space of admissible $\bbS$ and this smoothness allows a consistent calculation of the Fr\'echet derivative of $\bbT(\bbS)$. This can be considered as a consequence of the fact that deformations between arbitrary spaces can be measured according to the topology of the space. In particular, if a diffeomorphism is used to produce deformation in the domain of the signals of interest, it is possible to find an equivalent associated map that produces deformation of the set of operators acting on the signal. If notions of differentiability are used to measure the size of the original diffeomorphism, it is natural to find similar notions involved on the map acting on the operators, but with the difference that the differentiability is measured according to the topology of the new space.
\end{remark}

\section{Algebraic Neural Networks}
\subsection{Convolution Feaures}
\subsection{Pooling}
\section{Conclusions}\label{sec:conclusions}

We considered algebraic neural networks (AlgNN) with commutative algebras as a tool to unify convolutional architectures like CNNs and GNNs, synthesizing the algebraic structure by exploiting results from the representation theory of algebras and algebraic signal processing. Within this framework, we showed that AlgNNs can, in general,  be stable to different types of perturbations, and the conditions under which the AlgNN operators are stable are determined by subsets of the algebra. We pointed out that the perturbations of the domain of the signals can be equivalently modeled as a perturbation of the representation or the signal model, and the degree of this perturbation can be measured by means of the Fr\'echet derivative of two functions, the image of the homomorhisms in $\text{End}(\mathcal{M})$ and the perturbation model $\mathbf{T}(\mathbf{S})$.  The perturbation model considered provides enough expressive power to represent a wide variety of perturbations affecting the domain of the signals or the operator themselves directly. In particular, when considering the algebraic model for GNNs, the absolute and the relative perturbation models can be considered particular cases of the perturbation model considered in this work.

\section{Discussion}\label{sec:discussion}
\begin{figure}[t]
    \centering
    \input{./source_figs/plots_stability/graph_dilation_unstable.tex}\medskip
    \input{./source_figs/plots_stability/graph_dilation_stable.tex}
   \caption{Filter properties and stability for algebraic operators considering algebras with a single generator. (Top) We depict a Lipschitz filter where it is possible to see that an arbitrary degree of selectivity can be achieved in any part of the spectrum. (bottom) We depict a Lipschitz integral filter where we can see how the magnitude of the filters tends to a constant value as the size of $\vert\lambda_{i}\vert$ grows. As a consequence there is no discriminability in one portion of the spectrum.}
 \label{fig_graph_dilation}
\end{figure}
The mathematical form of the notion of stability introduced in definition~(\ref{def:stabilityoperators1}), eqn.~(\ref{eq:stabilityoperators1}) is uncannily similar to the expressions associated to the stability conditions stated in~\cite{mallat_ginvscatt,bruna_iscn} when the perturbation operator $\tau$  considered was affecting directly the domain of the signals.  This is consistent with the fact that  the size of the perturbation on the operators is the size of an induced diffeomorphism $\mathbf{T}$ acting on $\text{End}(\mathcal{M})$. Measuring the size of perturbations in this way, although less intuitive, provides an alternative way to handle and interpret perturbations on irregular domains.

The nature and severity of the perturbations, imposes restrictions on the behavior of the filters needed to guarantee stability. The more complex and severe the perturbation is the more conditions  on the filters are necessary to guarantee stability. This in particular has implications regarding to the selectivity of the filters in some specific frequency bands. The trade-off between stability and selectivity in the filters of the AlgNN can be measured by the size of the upper bounds in Theorems~\ref{theorem:HvsFrechet} up to~\ref{theorem:stabilityAlgNN1}, in particular, the size of $L_{0}$ and $L_{1}$ associated to the boundedness of the derivatives of the elements in $\mathcal{A}_{L0}$ and $\mathcal{A}_{L1}$. The smaller the value of $L_{0}, L_{1}$ the more stable the operators but the less selectivity we have. In Fig.~(\ref{fig_graph_dilation}) the properties in frequency of Lipschtiz and Lipschitz integral filters are depicted, where it is possible to see how the selectivity on portions of the spectrum is affected by properties that at the same time provide stability conditions for the perturbation models considered.

It is important to remark that the function $\sigma_{\ell}=P_{\ell}\circ\eta_{\ell}$ composed by the projection operator $P_{\ell}$ and the nonlinearity function $\eta_{\ell}$ relocates information from one layer to the other performing a mapping between different portions of the spectrum associated to each of the spaces $\mathcal{M}_{\ell}$.  As $\eta_{\ell}$ maps elements of $\mathcal{M}_{\ell}$ onto itself,  we can see in light of the decomposition of $\mathcal{M}_{\ell}$ in terms of irreducible representations that
$\eta_{\ell}$ is nothing but a relocator of information from one portion of the spectrum to the other. Additionally, the simplicity of $\eta_{\ell}$ provides  a rich variety of choices that can be explored in future research.

The notion of differentiability between metric spaces or Banach spaces can be considered also using the notion of \textit{G\^ateaux derivative} which is considered a \text{weak} notion of differentiability. Although Gateaux differentiability is in general different from Fr\'echet  differentiability, it is possible to show that when $\text{dim}\left(\text{End}(\mathcal{M})\right)<\infty$ both notions are equivalent for Lipschtiz functions, but substantial differences may exist if $\text{dim}\left(\text{End}(\mathcal{M})\right)=\infty$ even if the functions are Lipschitz~\cite{benyamini2000geometric,lindenstrauss2012frechet}.

\section{Introduction}
The numerical evidence that shows the benefits and goodness of using convolutional neural networks (CNNs) and graph neural networks (GNNs) for machine learning applications grows year after year, rising interest among researchers on finding solid and consistent explanations for the performance and the properties of these structures. In this context, the stability analysis of the operators representing the networks plays a central role, and some insights exist in the literature for CNNs~\cite{mallat_ginvscatt,bruna_iscn,bruna_groupinvrepconvnn} and GNNs~\cite{zou_stability,gamabruna_diffscattongraphs,fern2019stability}, obtained independently but similar in form and nature. These results opened questions about whether these notions descend from a common notion of stability, and finding a connection or a general framework where these concepts can be unified have become of interest.

The notion of stability in CNNs is rooted in the notion of Lipschitz-continuity to the action of diffeomorphisms introduced by Mallat in~\cite{mallat_ginvscatt} for the analysis of translation-invariant operators acting on $L^{2}(\mathbb{R}^{n})$. This concept is adapted to invariant scattering convolutional networks~\cite{bruna_iscn}, where the stability condition also implies global translation invariance. Mild variants of these notions have been considered in the literature for the analysis of convolutional networks. In~\cite{bruna_groupinvrepconvnn} a version of the notion introduced in~\cite{mallat_ginvscatt} is used considering the action of groups, formulating the Lipschitz stability in terms of the distance between a given group $G$ and the invariance group $G_{0}$, while in~\cite{bietti2017_stabilitycnn} the stability condition is a variant of~\cite{mallat_ginvscatt} considering signals whose domain is restricted to a compact subset of $\mathbb{R}^{n}$. For GNNs the problem of formulating stability conditions has been considered in~\cite{zou_stability}, analyzing two scenarios, first the stability is measured with respect to changes in the amplitude of the signal, and second the perturbation is measured with respect to changes in the graphs. In~\cite{gamabruna_diffscattongraphs} the notion of stability on graphs is considered in depth pointing out that the generalization of the conditions stated in~\cite{mallat_ginvscatt, bruna_iscn} is not straightforward for non smooth, non Euclidean domains, and as a way to quantify stability in GNNs the notion of \textit{metric stability} is considered using a diffusion operator to measure the perturbations or changes in the graphs. In~\cite{fern2019stability} a related notion of stability is used to provide concrete results about the stability on GNNs.

Although apparently different in nature, the stability results for CNNs and GNNs are uncannily similar not only in their final formulation but also in the analysis that lead to such results. Therefore, one may ask: is there any common framework where these proofs and analysis can be derived as particular cases?, can these stability results be extended to other convolutional architectures on different domains and with different notions of convolutions?. This paper addresses these questions, analyzing  \textit{algebraic neural networks (AlgNN)} to provide an affirmative answer. In particular, we analyze the stability properties of AlgNNs deriving general results whose particular instances lead to those results known for GNNs and CNNs. An AlgNN is a layered structure where the processing of the data in each layer is done in the context of an \textit{algebraic signal model} in the sense of~\cite{algSP0}. In particular, each layer is defined by an algebra $\mathcal{A}$, a vector space $\mathcal{M}$ where the signals, as elements of $\mathcal{M}$, are processed by elements of the set of endomorphisms $\text{End}(\mathcal{M})$, which are selected according to an homomorphism $\rho:\mathcal{A}\rightarrow\text{End}(\mathcal{M})$. Using the map that relates an operator and its perturbed version, we introduce a notion of stability that uses the Fr\'echet derivative of this map and we provide stability bounds which depend on the algebraic structure of the signal models and whose particular instantiations lead to results like the one provided in~\cite{fern2019stability}. For this analysis we consider a perturbation model that contains the already studied perturbation models analyzed in~\cite{zou_stability,gamabruna_diffscattongraphs,fern2019stability}. We remark as the most important contributions of this paper:
\begin{itemize}
\item The introduction of AlgNNs to unify stability analysis of neural networks where formal convolutions are used. In particular, we introduce the notion of perturbations and stability to algebraic signal models~(definition~\ref{def:stabilityoperators1}), where the size of the perturbation is measured taking into account the Fr\'echet derivative of an operator that relates a filter with its perturbed version. This notion is compatible with the notion of stability considered in~\cite{mallat_ginvscatt,bruna_iscn,bruna_groupinvrepconvnn,zou_stability,gamabruna_diffscattongraphs,fern2019stability}.
\item Stability results for algebraic filters considering commutative algebras~(theorems~\ref{theorem:HvsFrechet} to~\ref{theorem:uppboundDHmultg} and corollaries~\ref{corollary:DHboundsprcases} and~\ref{corollary:stabmultgenfilt}). As a consequence of these results we obtain also stability results for AlgNNs
(theorems~\ref{theorem:stabilityAlgNN0} and~\ref{theorem:stabilityAlgNN1}), whose particular instantiations lead to known results for CNNs and GNNs.
\end{itemize}

The results presented in theorems~\ref{theorem:HvsFrechet} to~\ref{theorem:uppboundDHmultg} and corollaries~\ref{corollary:DHboundsprcases} and~\ref{corollary:stabmultgenfilt} have several important implications among which are the following. First, they show that the stability properties in convolutional architectures can be expressed in terms of the algebraic laws (structure of $\mathcal{A}$) that govern the signal model in each layer. This fact highlights that the property of stability is universal always that the algebra remains the same, and it also explains the remarkable similarities between the results obtained for GNNs and CNNs. Second, the changes in the filters acting on a given signal are determined by the Fr\'echet derivative of the filters as functions of the images of the generators in the algebra, and this is true no matter which algebra is used and what particular representation is being considered. At the same time, this derivative as an operator acting on the perturbation determines whether the filters are stable or not. Third, although the perturbation is considered on the images of the algebra on the space of endomorphisms, the stability to a given perturbation is determined by the existence of certain subsets of the algebra itself, i.e. the stability relies on restrictions that are imposed on the filters and as a consequence the AlgNN can resolve the frequency selectivity without loosing stability. We remark that these ideas contribute to a deeper understanding of convolutional architectures and explain behaviors, results and analysis that exhibit similarities but that were conceived by different means and tools due to an apparent disconnection between them, for instance like with CNNs and GNNs.

This paper is organized as follows. In Section~\ref{sec_alg_filters} we start discussing the basics about algebraic signal processing and algebraic filtering introducing the basics about representation theory of algebras and its connection to the general algebraic signal model, we provide a set of examples to make basic concepts and ideas clear to the reader. In Section~\ref{sec_Algebraic_NNs} we introduce algebraic neural networks discussing how this notion generalizes CNNs, Group-CNNs, Graphon-CNNs, and GNNs which are considered as particular instantiations. The notions of perturbations and stability are discussed in Section~\ref{sec:perturbandstability} where the formal definition of stability is introduced with the perturbation model considered for our discussion. Additionally,  several stability theorems are presented in Section~\ref{sec:stabilitytheorems} considering filters in the algebraic signal model and operators representing AlgNNs as well. In Section~\ref{sec:spectopt} we discuss the basics about the spectral representation of algebraic filters, relevant to Section~\ref{sec:proofofTheorems} where the proofs of the theorems in Section~\ref{sec:perturbandstability} are presented. Section~\ref{sec:discussion} offers a discussion about the results pointing out the role of subsets of the algebra and the stability to perturbations given a perturbation model. Finally, some conclusions are provided in Section~\ref{sec:conclusions}.

\section{Multiple Generators}

\begin{theorem}\label{theorem:uppboundDHmultg}
Let $\mathcal{A}$ be an algebra with $m$ generators $\{g_{i}\}_{i=1}^{m}$ and $g_{i}g_{j}=g_{j}g_{i}$ for all $i,j\in\{1,\ldots m\}$. Let $(\mathcal{M},\rho)$ be a finite or countable infinite dimensional representation of $\mathcal{A}$ and $(\mathcal{M},\tilde{\rho})$ a perturbed version of $(\mathcal{M},\rho)$ related by the perturbation model in eqn.~(\ref{eq:perturbmod1}).  Then, if $p_{\mathcal{A}}\in\mathcal{A}_{L_{0}}\cap\mathcal{A}_{L_{1}}$ it holds that
\begin{multline}
\left\Vert D_{p_{\mathcal{M}}\vert\mathbf{S}_{i}}(\mathbf{S})\bbT(\bbS_{i})\right\Vert\leq
\\
(1+\delta_{0})L_{0}\sup_{\bbS_{i}\in\mathcal{S}}\Vert\bbT(\bbS_{i})\Vert +\left(1+q\delta\right)E_{q}L_{1}\sup_{\bbS_{i}\in\mathcal{S}}\Vert D_{\bbT}(\bbS_{i})\Vert
\label{eq:DHTSmultg}
\end{multline}
 \end{theorem}

--- CALCULATING the Frechet Derivative-------
\begin{proof}
$p(\mathbf{S})=p(\mathbf{S}_{1},\mathbf{S}_{2},\ldots,\mathbf{S}_{m})$, then

\begin{multline}
p(\mathbf{S}+\boldsymbol{\xi})-p(\mathbf{S})
\\
=\sum_{k_{1},\ldots,k_{m}=0}^{\infty}h_{k_{1},\ldots k_{m}}(\mathbf{S}_{1}+\boldsymbol{\xi}_{1})^{k_{1}}(\mathbf{S}_{2}+\boldsymbol{\xi}_{2})^{k_{2}}\cdots (\mathbf{S}_{m}+\boldsymbol{\xi}_{m})^{k_{m}}
\\
-\sum_{k_{1},\ldots,k_{m}=0}^{\infty}h_{k_{1},\ldots k_{m}}\mathbf{S}_{1}^{k_{1}}\mathbf{S}_{2}^{k_{2}}\ldots\mathbf{S}_{m}^{k_{m}}
\\
=D_{p\vert\mathbf{S}_{1}}(\mathbf{S})+\ldots +D_{p\vert\mathbf{S}_{m}}(\mathbf{S})+N(\boldsymbol{\xi})
\end{multline}
\end{proof}
}
\section{Proof of Theorems}\label{sec:proofofTheorems}
Let us start defining some notation. Let $\boldsymbol{\pi}(a_{1}\mathbf{A}_{1},\ldots,a_{r}\mathbf{A}_{r})$ the operator that represents the sum of all the products of the operators $\mathbf{A}_{1},\ldots,\mathbf{A}_{r}$ that appear $a_{1}, a_{2}, \ldots, a_{r}$ respectively. For instance, $\boldsymbol{\pi}(2\mathbf{\mathbf{A}},\mathbf{B})=\mathbf{AAB}+\mathbf{ABA}+\mathbf{BAA}$. Additionally, when considering all summation and product symbols the following convention is used $\sum_{i=a}^{b}F(i)=0$ if $b<a$, and $\prod_{i=a}^{b}F(i)=0$ if $b<a$.
\subsection{Proof of Theorems~\ref{theorem:HvsFrechet} and~\ref{theorem:HvsFrechetmultgen}}
\label{prooftheoremHvsFrechet}
\begin{proof}
We say that $p(\mathbf{S})$ as a function of $\mathbf{S}$ is Fr\'echet differentiable at $\mathbf{S}$ if there exists a bounded linear operator $D_{p}:\text{End}(\mathcal{M})^{m}\rightarrow\text{End}(\mathcal{M})$  such that~\cite{benyamini2000geometric,lindenstrauss2012frechet}
\begin{equation}
\lim_{\Vert\boldsymbol{\xi}\Vert\to 0}\frac{\left\Vert p(\mathbf{S}+\boldsymbol{\xi})-p(\mathbf{S})-D_{p}(\mathbf{S})\left\lbrace\boldsymbol{\xi}\right\rbrace\right\Vert}{\Vert\boldsymbol{\xi}\Vert}=0
\end{equation}
which in Landau notation can be written as
\begin{equation}
p(\mathbf{S}+\boldsymbol{\xi})=p(\mathbf{S})+D_{p}(\mathbf{S})\left\lbrace\boldsymbol{\xi}\right\rbrace+o(\Vert\boldsymbol{\xi}\Vert).
\label{eq:frechetDdef}
\end{equation}
Calculating the norm in eqn.~(\ref{eq:frechetDdef}) and applying the triangle inequality we have:
 $\left\Vert p(\mathbf{S}+\boldsymbol{\xi})-
p(\mathbf{S})\right\Vert\leq
\left\Vert D_{p}(\mathbf{S})\left\lbrace\boldsymbol{\xi}\right\rbrace\right\Vert+\mathcal{O}\left(\Vert\boldsymbol{\xi}\Vert^{2}\right)$
for all $\boldsymbol{\xi}=(\boldsymbol{\xi}_{1},\ldots,\boldsymbol{\xi}_{m})\in\text{End}(\mathcal{M})^{m}$. Now, taking into account that (see~\cite{berger1977nonlinearity} pages 69-70)
\begin{equation}
\Vert D_{p}(\mathbf{S})\left\lbrace\boldsymbol{\xi}\right\rbrace\Vert\leq
\sum_{i=1}^{m}\left\Vert D_{p\vert\mathbf{S}_{i}}(\mathbf{S})\left\lbrace\boldsymbol{\xi}_{i}\right\rbrace\right\Vert
\end{equation}
we have
\begin{equation*}
\left\Vert p(\mathbf{S}+\boldsymbol{\xi})-
p(\mathbf{S})\right\Vert\leq
\sum_{i=1}^{m}\left\Vert D_{p\vert\mathbf{S}_{i}}(\mathbf{S})\left\lbrace\boldsymbol{\xi}_{i}\right\rbrace\right\Vert+\mathcal{O}\left(\Vert\boldsymbol{\xi}\Vert^{2}\right),
\end{equation*}
where $D_{p\vert\mathbf{S}_{i}}(\mathbf{S})$ is the partial Frechet derivative of $p(\mathbf{S})$ on $\mathbf{S}_{i}$. Then, taking into account that
$\left\Vert p(\mathbf{S}+\boldsymbol{\xi})\mathbf{x}-
p(\mathbf{S})\mathbf{x}\right\Vert\leq\Vert\mathbf{x}\Vert\left\Vert p(\mathbf{S}+\boldsymbol{\xi})-
p(\mathbf{S})\right\Vert$ and selecting $\boldsymbol{\xi}_{i}=\mathbf{T}(\mathbf{S}_{i})$ we complete the proof.
\end{proof}
\subsection{Proof of Theorem~\ref{theorem:uppboundDH}}
\label{prooftheorem:uppboundDH}
\begin{proof}
Taking into account the definition of the Fr\'echet derivative of $p$ on $\mathbf{S}_{i}$ (see Appendix~\ref{appendix:howtofindDH}) we have
\begin{equation*}
\left\Vert D_{p\vert\mathbf{S}_{i}}(\mathbf{S})\left\lbrace\mathbf{T}(\mathbf{S}_{i})\right\rbrace\right\Vert=\left\Vert\sum_{k_{i}=1}^{\infty}\mathbf{A}_{k_{i}}\boldsymbol{\pi}\left(\mathbf{T}(\mathbf{S}_{i}),(k_{i}-1)\mathbf{S}_{i}\right)\right\Vert,
\end{equation*}
and re-organizating terms we have
\begin{equation*}
\left\Vert D_{p\vert\mathbf{S}_{i}}(\mathbf{S})\left\lbrace\mathbf{T}(\mathbf{S}_{i})\right\rbrace\right\Vert=\left\Vert\sum_{\ell=1}^{\infty}\mathbf{S}_{i}^{\ell-1}\mathbf{T}(\mathbf{S}_{i})\sum_{k_{i}=\ell}^{\infty}\mathbf{A}_{k_{i}}\mathbf{S}_{i}^{k_{i}-\ell}\right\Vert.
\end{equation*}
Taking into account eqn.~(\ref{eq:PrvsTr}), it follows that
\begin{multline}
\left\Vert D_{p\vert\mathbf{S}_{i}}(\mathbf{S})\left\lbrace\mathbf{T}(\mathbf{S}_{i})\right\rbrace\right\Vert=
\\
\left\Vert\sum_{\ell=1}^{\infty}\left(\mathbf{T}_{0c,i}\mathbf{S}_{i}^{\ell-1}+\mathbf{S}^{\ell-1}_{i}\mathbf{P}_{0,i}\right)\sum_{k_{i}=\ell}^{\infty}\mathbf{A}_{k_{i}}\mathbf{S}_{i}^{k_{i}-\ell}\right.
\\
\left.+\sum_{\ell=1}^{\infty}\left(\mathbf{T}_{1c,i}\mathbf{S}_{i}^{\ell}
+\mathbf{S}_{i}^{\ell-1}\mathbf{P}_{1,i}\mathbf{S}_{i}\right)
\sum_{k=\ell}^{\infty}\mathbf{A}_{k_{i}}\mathbf{S}^{k_{i}-\ell}\right\Vert.
\label{eq:auxDHT1}
\end{multline}
Applying the triangle inequality and distribuiting the sum we have
\begin{multline}\label{eq:DpDTboundeqn1}
\left\Vert D_{p\vert\mathbf{S}_{i}}(\mathbf{S})\left\lbrace\mathbf{T}(\mathbf{S}_{i})\right\rbrace\right\Vert
\leq
\left\Vert\mathbf{T}_{0c,i}\sum_{\ell=1}^{\infty}\sum_{k_{i}=\ell}^{\infty}\mathbf{S}_{i}^{k_{i}-1}\mathbf{A}_{k_{i}}\right\Vert
\\
+\left\Vert D_{p\vert\mathbf{S}_{i}}(\mathbf{S})\left\lbrace\mathbf{P}_{0,i}\right\rbrace\right\Vert
+\left\Vert\mathbf{T}_{1c,i}\sum_{\ell=1}^{\infty}\sum_{k_{i}=\ell}^{\infty}\mathbf{S}^{k_{i}}\mathbf{A}_{k_{i}}\right\Vert
\\
+\left\Vert D_{p\vert\mathbf{S}_{i}}(\mathbf{S})\left\lbrace\mathbf{P}_{0,i}\mathbf{S}_{i}\right\rbrace\right\Vert
\end{multline}
Now, we analyze term by term in eqn.~(\ref{eq:DpDTboundeqn1}). For the first term we take into account that $\sum_{\ell=1}^{\infty}\sum_{k_{i}=\ell}^{\infty}\mathbf{S}_{i}^{k_{i}-1}\mathbf{A}_{k_{i}}=\sum_{k_{i}=1}^{\infty}k_{i}\mathbf{A}_{k_{i}}\mathbf{S}_{i}^{k_{i}-1}$ and we apply the product norm property taking into account that the filters belong to $\mathcal{A}_{L_{0}}$, which leads to
\begin{multline}
\left\Vert\mathbf{T}_{0c,i}\sum_{\ell=1}^{\infty}\sum_{k_{i}=\ell}^{\infty}\mathbf{S}_{i}^{k_{i}-1}\mathbf{A}_{k_{i}}\right\Vert
\\
\leq
\left\Vert\mathbf{T}_{0c,i}\right\Vert\left\Vert\sum_{k_{i}=1}^{\infty}k_{i}\mathbf{A}_{k_{i}}\mathbf{S}_{i}^{k_{i}-1}\right\Vert\leq L_{0}\Vert\mathbf{T}_{0,i}\Vert.
\end{multline}
For the second term in eqn.~(\ref{eq:DpDTboundeqn1}) we take into account that (see \cite{benyamini2000geometric} page 84, \cite{lindenstrauss2012frechet} page 158 and~\cite{Lindenstrauss2010FrchetDO} page 386): $\left\Vert D_{p\vert\mathbf{S}_{i}}(\mathbf{S})\left\lbrace\mathbf{P}_{0,i}\right\rbrace\right\Vert
\leq L_{0}\Vert\mathbf{P}_{0,i}\Vert$ if $p$ is G\^ateaux differentiable, which is always true because $p$ is Fr\'echet differentiable and with the fact that $\Vert\mathbf{P}_{0,i}\Vert\leq\delta\Vert\mathbf{T}_{0,i}\Vert$, we have $\left\Vert D_{p\vert\mathbf{S}_{i}}(\mathbf{S})\left\lbrace\mathbf{P}_{0,i}\right\rbrace\right\Vert \leq L_{0}\delta\Vert\mathbf{T}_{0,i}\Vert$.

For the third term in eqn.~(\ref{eq:DpDTboundeqn1}), we take into account that $\sum_{\ell=1}^{\infty}\sum_{k_{i}=\ell}^{\infty}\mathbf{S}_{i}^{k_{i}}\mathbf{A}_{k_{i}}=\sum_{k_{i}=1}^{\infty}k_{i}\mathbf{A}_{k_{i}}\mathbf{S}_{i}^{k_{i}}$ and we apply the norm product property taking into account that the filters belong to $\mathcal{A}_{L_{1}}$, which leads to
\begin{multline}
\left\Vert\mathbf{T}_{1c,i}\sum_{\ell=1}^{\infty}\sum_{k_{i}=\ell}^{\infty}\mathbf{S}_{i}^{k_{i}}\mathbf{A}_{k_{i}}\right\Vert
\\
\leq
\left\Vert\mathbf{T}_{1c,i}\right\Vert\left\Vert\sum_{k_{i}=1}^{\infty}k_{i}\mathbf{A}_{k_{i}}\mathbf{S}_{i}^{k_{i}}\right\Vert\leq L_{1}\Vert\mathbf{T}_{1,i}\Vert.
\end{multline}
Finally, for the fourth term we use the notation $\tilde{D}(\mathbf{S})\left\lbrace\mathbf{P}_{1,i}
\right\rbrace
=D_{p\vert\mathbf{S}_{i}}(\mathbf{S})\left\lbrace\mathbf{P}_{1,i}\mathbf{S}_{i}\right\rbrace$. We start taking into account that  (see~\cite{higham2008functions} pages 61 and 331) the eigenvalues of the operator $\tilde{D}(\mathbf{S})$ represented as $\zeta_{pq}$ are given by
\begin{equation}
\zeta_{pq}=\left\lbrace
\begin{array}{ccc}
\frac{p(\lambda_{p})-p(\lambda_{q})}{\lambda_{p}-\lambda_{q}}\lambda_{q} & \text{if} & \lambda_{p}\neq\lambda_{q}\\
\lambda_{p}p^{'}(\lambda_{p}) & \text{if} & \lambda_{p}=\lambda_{q}
\end{array}
\right. ,
\end{equation}
then, taking into account that the filters belong to $\mathcal{A}_{L_{1}}$ we have $\Vert\tilde{D}(\mathbf{S})\Vert\leq L_{1}$, therefore $\left\Vert
D_{p\vert\mathbf{S}_{i}}(\mathbf{S})\left\lbrace\mathbf{P}_{1,i}\mathbf{S}_{i}\right\rbrace\right\Vert=\left\Vert\tilde{D}(\mathbf{S})\left\lbrace\mathbf{P}_{1,i}
\right\rbrace\right\Vert\leq L_{1}\Vert\mathbf{P}_{1,i}\Vert$. Additionally, with $\Vert\mathbf{P}_{1,i}\Vert\leq\delta\Vert\mathbf{T}_{1,i}\Vert$ it follows that $\left\Vert
D_{p\vert\mathbf{S}_{i}}(\mathbf{S})\left\lbrace\mathbf{P}_{1,i}\mathbf{S}_{i}\right\rbrace\right\Vert\leq L_{1}\delta\Vert\mathbf{T}_{1,i}\Vert$.

Putting all these results together into eqn.~(\ref{eq:DpDTboundeqn1}) we reach
\begin{multline*}
\left\Vert D_{p\vert\mathbf{S}_{i}}(\mathbf{S})\left\lbrace\mathbf{T}(\mathbf{S}_{i})\right\rbrace\right\Vert
\leq
(1+\delta)L_{0}\Vert\mathbf{T}_{0,i}\Vert+(1+\delta)L_{1}\Vert\mathbf{T}_{1,i}\Vert
\\
\leq
(1+\delta)\left(L_{0}\sup_{\mathbf{S}_{i}\in\mathcal{S}}\Vert\mathbf{T}(\mathbf{S}_{i})\Vert+L_{1}\sup_{\mathbf{S}_{i}\in\mathcal{S}}\Vert D_{\mathbf{T}}(\mathbf{S}_{i})\Vert\right)
\end{multline*}
\end{proof}
\subsection{Proof of Theorems~\ref{theorem:stabilityAlgNN0} and~\ref{theorem:stabilityAlgNN1}}
\subsubsection{Proof of Theorem~\ref{theorem:stabilityAlgNN0}}
\label{prooftheorem:stabilityAlgNN0}
\begin{proof}
Taking into account eqns.~(\ref{eq:xl}), and~(\ref{eq:interlayeropalgnn}) and the fact that the maps $\sigma_{\ell}$ are Lipschitz with constant $C_{\ell}$ we have that
\begin{equation}
\left\Vert
\sigma_{\ell}\left(p(\mathbf{S}_{\ell})\mathbf{x}_{\ell-1}\right)-\sigma_{\ell}\left(p(\tilde{\mathbf{S}}_{\ell})\mathbf{x}_{\ell-1}\right)
\right\Vert
\leq
C_{\ell}\boldsymbol{\Delta}_{\ell}\Vert\mathbf{x}_{\ell-1}\Vert,
\end{equation}
where $\boldsymbol{\Delta}_{\ell}=\Vert p(\mathbf{S}_{\ell})-p(\tilde{\mathbf{S}}_{\ell})\Vert$, and whose value is determined by theorems~\ref{theorem:HvsFrechet} and~\ref{theorem:uppboundDH}.
\end{proof}

\subsubsection{Proof of Theorem~\ref{theorem:stabilityAlgNN1}}
\label{prooftheorem:stabilityAlgNN1}
\begin{proof}
Before starting the calculations let us introduce some notation. Let $\varphi_{f}(\ell,g)=\rho_{\ell}\left(\xi^{fg}\right)$ denote the image of the filter $\xi^{fg}\in\mathcal{A}_{\ell}$ that process the $g$th feature coming from the layer $\ell-1$  and that is associated to $f$th feature in layer $\ell$. As indicated before, $\sigma_{\ell}$ indicates the Lipschitz mapping from layer $\ell$ to layer $\ell+1$. The term $\mathbf{x}_{\ell-1}^{g}$ indicates the $g$th feature in the layer $\ell-1$. Then, we have that:
\begin{multline}
\left\Vert\mathbf{x}_{\ell}^{f}-\tilde{\mathbf{x}}_{\ell}^{f}\right\Vert\leq
\left\Vert\sigma_{\ell-1}\sum_{g_{\ell-1}}\varphi_{g_{\ell}}(\ell-1,g_{\ell-1})\sigma_{\ell-2}\right.\\
\left.\sum_{g_{\ell-2}}\varphi_{g_{\ell-1}}(\ell-2,g_{\ell-2})
\cdots\sigma_{1}\sum_{g_{1}}\varphi_{g_{2}}\mathbf{x}\right.-\\
\left. \sigma_{\ell-1}\sum_{g_{\ell-1}}\tilde{\varphi}_{g_{\ell}}(\ell-1,g_{\ell-1})\sigma_{\ell-2}\sum_{g_{\ell-2}}\tilde{\varphi}_{g_{\ell-1}}(\ell-2,g_{\ell-2})\right.\\
\left.\cdots\sigma_{1}\sum_{g_{1}}\tilde{\varphi}_{g_{2}}(1,g_{1})\mathbf{x}
\right\Vert.
\label{eq:longalgNNexp1}
\end{multline}
In order to exapand eqn.~(\ref{eq:longalgNNexp1}) we start pointing out that:
\begin{multline}
A_{\ell+1}\sigma_{\ell}(a)-\tilde{A}_{\ell+1}\sigma_{\ell}(\tilde{a})=
\\
(A_{\ell+1}-\tilde{A}_{\ell+1})\sigma_{\ell}(a)+\tilde{A}_{\ell+1}(\sigma_{\ell}(a)-\sigma_{\ell}(\tilde{a}))
\end{multline}
where $A_{\ell+1}$ and $\tilde{A}_{\ell+1}$ indicate filter operators and their perturbed versions, respectively. Now, noticing that $\Vert\sigma_{\ell}(a)-\sigma_{\ell}(b)\Vert\leq C_{\ell}\Vert a-b\Vert$, $\Vert A_{\ell+1}-\tilde{A}_{\ell+1}\Vert\leq\boldsymbol{\Delta}_{\ell}$ and $\Vert A_{\ell+1}\Vert\leq B_{\ell+1}$ we have the following relations
\begin{equation}
\sum_{g_{k}}\Vert\alpha-\tilde{\alpha}\Vert\leq\sum_{g_{k}}\left(
\boldsymbol{\Delta}_{k}\Vert\sigma_{k-1}(\alpha)\Vert+B_{k}C_{k-1}\Vert\beta-\tilde{\beta}\Vert
\right)
\label{eq:aux_longeq1}
\end{equation}
\begin{equation}
\sum_{g_{k}}\Vert\beta-\tilde{\beta}\Vert\leq\sum_{g_{k}}\sum_{g_{k-1}}\Vert\alpha-\tilde{\alpha}\Vert
\label{eq:aux_longeq2}
\end{equation}
\begin{equation}
\sum_{g_{k}}\Vert\sigma_{k-1}(\alpha)\Vert\leq\left(\prod_{r=1}^{k}F_{r}\right)\left(\prod_{r=1}^{k-1}C_{r}B_{r}\right)\Vert\mathbf{x}\Vert
\label{eq:aux_longeq3}
\end{equation}
where $\alpha$ and $\tilde{\alpha}$ represent sequences of symbols in eqn.~(\ref{eq:longalgNNexp1}) that start with a symbol of the type $\varphi$, while $\beta$ and $\tilde{\beta}$ indicate a sequence of symbols that start with a summation  symbol, and the tilde makes reference to symbols that are associated to the perturbed representations. The term $\boldsymbol{\Delta}_{\ell}$ is associated to the difference between the operators and their perturbed versions (see definition~\ref{def:stabilityoperators1}) in the layer $\ell$ and whose values are given in Theorems~\ref{theorem:HvsFrechet} and~\ref{theorem:uppboundDH}. Combining eqns.~(\ref{eq:aux_longeq1}),~(\ref{eq:aux_longeq2}) and~(\ref{eq:aux_longeq3}) we have:
\begin{multline}
\left\Vert\mathbf{x}_{L}^{f}-\tilde{\mathbf{x}}_{L}^{f}
\right\Vert\leq
\sum_{\ell=1}^{L}\boldsymbol{\Delta}_{\ell}\left(\prod_{r=\ell}^{L}C_{r}\right)\left(\prod_{r=\ell+1}^{L}B_{r}\right)\\
\left(\prod_{r=\ell}^{L-1}F_{r}\right)\left(\prod_{r=1}^{\ell-1}C_{r}F_{r}B_{r}\right)\left\Vert\mathbf{x}\right\Vert,
\end{multline}
where the products $\prod_{r=a}^{b}F(r)=0$ if $b<a$. Now taking into account that
\begin{equation*}
\left\Vert
\Phi\left(\mathbf{x},\{ \mathcal{P}_{\ell} \}_{1}^{L},\{ \mathcal{S}_{\ell}\}_{1}^{L}\right)-
\Phi\left(\mathbf{x},\{ \mathcal{P}_{\ell} \}_{1}^{L},\{ \tilde{\mathcal{S}}_{\ell}\}_{1}^{L}\right)
\right\Vert^{2}
=
\sum_{f=1}^{F_{L}}\left\Vert\mathbf{x}_{L}^{f}-\tilde{\mathbf{x}}_{L}^{f}
\right\Vert^{2}
\end{equation*}
we have
\begin{multline}
\left\Vert
\Phi\left(\mathbf{x},\{ \mathcal{P}_{\ell} \}_{1}^{L},\{ \mathcal{S}_{\ell}\}_{1}^{L}\right)-
\Phi\left(\mathbf{x},\{ \mathcal{P}_{\ell} \}_{1}^{L},\{ \tilde{\mathcal{S}}_{\ell}\}_{1}^{L}\right)
\right\Vert
\\
\leq
\sqrt{F_{L}}
\sum_{\ell=1}^{L}\boldsymbol{\Delta}_{\ell}\left(\prod_{r=\ell}^{L}C_{r}\right)\left(\prod_{r=\ell+1}^{L}B_{r}\right)
\left(\prod_{r=\ell}^{L-1}F_{r}\right)
\\
\left(\prod_{r=1}^{\ell-1}C_{r}F_{r}B_{r}\right)\left\Vert\mathbf{x}\right\Vert
\end{multline}
\end{proof}

\section{Algebraic Filters} \label{sec_alg_filters

Algebraic signal processing (ASP) provides a framework for understanding and generalizing traditional signal processing exploiting the representation theory of algebras~\cite{algSP0, algSP1, algSP2, algSP3}. In ASP, a signal model is defined as the triple
\begin{equation}\label{eqn_ASP_signal_model}
   (\ccalA,\ccalM,\rho),
\end{equation}
in which $\ccalA$ is an associative algebra with unity, $\ccalM$ is a vector space with inner product, and $\rho:\ccalA\rightarrow\text{End}(\ccalM)$ is a homomorphism between the algebra $\ccalA$ and the set of endomorphisms in the vector space $\ccalM$. The elements in \eqref{eqn_ASP_signal_model} are tied together by the notion of a representation which we formally define next.

\begin{definition}[Representation] \label{def_representation} A representation of the associative algebra $\ccalA$ is a vector space $\ccalM$ equipped with a homomorphism $\rho:\ccalA\rightarrow\text{End}(\ccalM)$, i.e., a linear map preserving multiplication and unit.
\end{definition}
The class of representations of  an algebra $\mathcal{A}$ is denoted as $\mathsf{Rep}\{\mathcal{A}\}$.

In an ASP model, signals  are elements of the vector space $\ccalM$, and filters are elements of the algebra $\ccalA$. Thus, the vector space $\ccalM$ deteremines the objects of interest and the algebra $\ccalA$ the rules of the operations that define a filter. The homomorphism $\rho$ translates the abstract operators $a\in\ccalA$ into concrete operators $\rho(a)$ that act on signals $\bbx$ to produce filter outputs
\begin{equation}\label{eqn_ASP_filter_outputs}
   \bby = \rho(a) \bbx.
\end{equation}
The algebraic filters in \eqref{eqn_ASP_filter_outputs} generalize the convolutional processing of time signals -- see Example \ref{ex_DTSP}. Our goal in this paper is to use them to generalize convolutional neural networks (Section \ref{sec_Algebraic_NNs}) and to study their fundamental stability properties (Section \ref{sec:perturbandstability}). Generators, which we formally define next, are important for the latter goal.

\begin{definition}[Generators] \label{def_generators} For an associative algebra with unity $\ccalA$ we say the set $\ccalG\subseteq\ccalA$ generates $\ccalA$ if all $a\in\ccalA$ can be represented as polynomial functions of the elements of $\ccalG$. We say elements $g\in\ccalG$ are generators of $\ccalA$ and we denote as $a=p_{\ccalA}(\ccalG)$  the polynomial that generates $a$.
\end{definition}

Definition \ref{def_generators} states that elements $a\in\ccalA$ can be built from the generating set as polynomials using the operations of the algebra. Given the role of representations in connecting the algebra $\ccalA$ with the signals $\bbx$ as per Definition \ref{def_representation} the representation $\rho(g)$ of a generator $g\in\ccalG$ will be of interest. In the context of ASP, these representations are called shift operators as we formally define next.

\begin{definition}[Shift Operators]\label{def_shift_operators}
Let $(\ccalM,\rho)$ be a representation of the algebra $\ccalA$. Then, if $\ccalG\subset\ccalA$ is a generator set of $\ccalA$, the operators $\bbS = \rho(g)$ with $g\in\ccalG$ are called shift operators. The set fo all shift operators is denoted by $\ccalS$
\end{definition}

Given that elements $a$ of the algebra are generated from elements $g$ of the generating set, it follows that filters $\rho(a)$ are generated from the set of shift operators $\bbS=\rho(g)$. In fact, if we have that $a=p_{\ccalA}(\ccalG)$ is the polynomial that generates $a$, the filter $\rho(a)$ is generated by
\begin{equation}\label{eqn_filter_representation}
   \rho(a) = p_\ccalM \big(\rho(\ccalG)\big)
           = p_\ccalM\big(\ccalS\big),
\end{equation}
where the subindex $\ccalM$ signifies that the operations in \eqref{eqn_filter_representation} are those of the vector space $\ccalM$ -- in contrast to the polynomial $a=p_{\ccalA}(\ccalG)$ whose operations are those of the algebra $\ccalA$.

To clarify ideas it is instructive to consider some examples.

\red{%
\begin{equation}\label{eqn_DTSP_filter_outputs}
   a = \sum_{k=0}^{K-1} h_k g^k
\end{equation}
}

\newcounter{dtsp}
\setcounter{dtsp}{\theexample}
\setcounter{example}{\thedtsp}

\begin{example}[Discrete Time Signal Processing]\label{ex_DTSP}\normalfont Let $\ccalM = \ell^2$ be the space of square summable sequences $\bbx = \{x_n\}_{n\in\mbZ}$ and $\ccalA$ the algebra of polynomials with elements $a = \sum_{k=0}^{K-1} h_k t^k$. We consider the shift operator $S$ such that $S\bbx$ is the sequence with entries $(S\bbx)_n = x_{n-1}$. With the homomorphism $\rho(a) = \sum_{k=0}^{K-1} h_k S^k$ the filters in \eqref{eqn_ASP_filter_outputs} are of the form
\begin{equation}\label{eqn_DTSP_filter_outputs}
   \bby  \,=\, \rho\left(\sum_{k=0}^{K-1} h_k t^k\right) \bbx
         \,=\, \left(\sum_{k=0}^{K-1} h_k S^k\right) \bbx
         \,=\, \sum_{k=0}^{K-1} h_k S^k \bbx.
\end{equation}
This is a representation of time convolutional filters \cite{algSP0}. Observe that the algebra of polynomials is generated by $g=t$ and that since $\rho(t) = S$ the shift operator $S$ generates the set of convolutional filters as it follows from \eqref{eqn_DTSP_filter_outputs}.\end{example}

\newcounter{gsp}
\setcounter{gsp}{\theexample}
\setcounter{example}{\thegsp}

\begin{example}[Graph Signal Processing]\label{ex_GSP}\normalfont We retain the algebra of polynomials as in Example \ref{ex_DTSP} but we change the space of signals to the set of complex vectors with $N$ entries, $\ccalM =\mbC^N$. We interpret components $x_n$ of $\bbx\in\ccalM = \mbC^N$ as being associated with nodes of a graph with matrix representation $\bbS\in\mbC^{N\times N}$. With the homomorphism $\rho(a) = \sum_{k=0}^{K-1} h_k \bbS^k$ the filters in~\eqref{eqn_ASP_filter_outputs} are of the form
\begin{equation}\label{eqn_GSP_filter_outputs}
   \bby \,=\, \rho\left(\sum_{k=0}^{K-1} h_k t^k\right) \bbx
        \,=\, \left(\sum_{k=0}^{K-1} h_k \bbS^k\right) \bbx
        \,=\, \sum_{k=0}^{K-1} h_k \bbS^k \bbx.
\end{equation}
This is a representation of the graph convolutional filters used in graph signal processing (GSP)~\cite{gsp_graphfilters}. Observe that \eqref{eqn_GSP_filter_outputs} and \eqref{eqn_DTSP_filter_outputs} are similar but represent different operations. In \eqref{eqn_GSP_filter_outputs} $\bbx$ is a vector and $\bbS^k$ a matrix power. In \eqref{eqn_DTSP_filter_outputs} $\bbx$ is a sequence and $S^k$ is the composition of the shift operator $S$. Their similarity arises from their common use of the algebra of polynomials. Their differences are because they use different vector spaces $\ccalM$ and different homomorphism functions $\rho$. \end{example}

\newcounter{dsp}
\setcounter{dsp}{\theexample}
\setcounter{example}{\thedsp}
\red{Discuss this with professor Alejandro. I wrote origininally here that $\mathcal{A}$ is the set of polynomials modulo $t^{N}-1$. I think we should review how the example is written.}
\begin{example}[Discrete Signal Processing]\label{ex_DSP}\normalfont Make $\ccalM=\mbC^N$ and $\ccalA$ the algebra of polynomials modulo $N$. Elements of $\ccalA$ are of the form $a = \sum_{k=0}^{K-1} h_k t^k$ but we must have $K<N$ and polynomial products use the rule $t^k = t^{k \mod N}$. Define the cyclic shift operator $C$ as one whose application to $\bbx\in\ccalM=\mbC^N$ yields a vector with components
$(C \bbx)_n = x_{(n-1)\mod N}$. With the homomorphism $\rho(a) = \sum_{k=0}^{K-1} h_k C^k$ the filters in~\eqref{eqn_ASP_filter_outputs} are of the form
\begin{equation}\label{eqn_DSP_filter_outputs}
   \bby \,=\, \rho\left(\sum_{k=0}^{K-1} h_k t^k\right) \bbx
        \,=\, \left(\sum_{k=0}^{K-1} h_k C^k\right) \bbx
        \,=\, \sum_{k=0}^{K-1} h_k C^k \bbx.
\end{equation}
This is the representation of (cyclic) convolutional filters in discrete signal processing~\cite{algSP0}. We emphasize that the homomorphism $\rho(a) = \sum_{k=0}^{K-1} h_k C^k$ is indeed a homomorphism because the cyclic shift operator $C$ satisfies $C^k = C^{k \mod N}$. \end{example}

As is clear from Examples \ref{ex_DTSP}-\ref{ex_DSP}, the effect of the operator $\rho(a)$ on a given signal $\bbx$ is determined by two factors: The filter $a\in\ccalA$ and the homomorphism $\rho$. The filter $a\in\ccalA$ indicates the \textit{laws and rules} to be used to manipulate the signal $\bbx$ and $\rho$ provides a \textit{physical realization} of the filter $a$ on the space $\ccalM$ to which $\bbx$ belongs. For instance, in these three examples the filter $a=1 + 2t$ indicates that the signal is to be added to a transformed version of the signal scaled by coefficient $2$. The homomorphism $\rho$ in Example \ref{ex_DTSP} dictates that the physical implementation of this transformation is a time shift. The homomorphism $\rho$ in Example \ref{ex_GSP} defines a transformation as a multiplication by $\bbS$ and in Example \ref{ex_DSP} the homomorphism entails a cyclic shift. We remark that in order to specify the physical effect of a filter it is always sufficient to specify the physical effect of the generators. In all three examples, the generator of the algebra is $g=t$. The respective effects of an arbitrary filter $a$ are determined once we specify that $\rho(t) = S$ in Example \ref{ex_DTSP}, $\rho(t) = \bbS$ in Example \ref{ex_GSP}, or $\rho(t) = C$ in Example \ref{ex_DSP}.

The flexibility in the choice of algebra and homomorphism allows for a very rich variety of signal processing frameworks. We highlight this richness with three more examples.

\newcounter{groupsp}
\setcounter{groupsp}{\theexample}
\setcounter{example}{\thegroupsp}

\begin{example}[Signal Processing on Groups]\normalfont \label{ex_groupsp} Let $\ccalM=\{\bbx: G\rightarrow\mathbb{C} \}=\{ \sum_{g\in G}\bbx(g)g\}$ be the set of functions defined on the group $G$ with values in $\mathbb{C}$ and $\ccalA=\ccalM$ the \textit{group algebra} with elements $\sum_{g\in G}\bba(g)g$. With a homomorphism given by $\rho(\bba)=\bba$ the filtering in~\eqref{eqn_ASP_filter_outputs} takes the form
\begin{equation}\label{eqn_group_filters}
\rho\left(\sum_{g\in G}\bba(g)g\right)\bbx=\sum_{g\in G}\bba(g)g\bbx=\sum_{g\in G}\sum_{h\in G}\bba(g)\bbx(h)gh,
\end{equation}
and making $u=gh$ we have
\begin{equation}\label{eqn_group_convolution}
   \sum_{g,h\in G} \bba(g)\bbx(h)gh
             =\sum_{u,h\in G}\bba(uh^{-1})\bbx(h)u=\bba\ast\bbx .
\end{equation}
This is the standard representation of convolution of signals on groups~\cite{steinbergrepg,terrasFG}.
\end{example}

\newcounter{wsp}
\setcounter{wsp}{\theexample}
\setcounter{example}{\thewsp}

\begin{example}[Graphon Signal Processing ]\normalfont \label{ex_wsp} Let $\ccalM=L^2([0,1])$ be the space of signals and $\ccalA$ the algebra of polynomials with elements $a = \sum_{k=0}^{K-1} h_k t^k$ with a homomorphism given by
$\rho(a)=\sum_{k=0}^{K-1} h_k S^k$ and
\begin{equation}\label{eqn_graphon_shift}
   \left(S\bbx\right)(u)=\int_{0}^{1}W(u,v)\bbx(v)dv,
\end{equation}
where $W(u,v):[0,1]^{2}\rightarrow [0,1]$ is a graphon, i.e. bounded symmetric measurable function. Then, the filtering in~\eqref{eqn_ASP_filter_outputs} is of the form~\eqref{eqn_DTSP_filter_outputs}.
This is a representation of convolutional filters in the context of graphon signal processing~\cite{ruiz2020graphon}.
\end{example}

\begin{example}[Lie Algebra Signal Processing ]\normalfont
Let $\ccalM=\mathbb{R}^{N}$ be the space of signals and
\begin{equation}
\ccalA=\frac{\mathbb{C}\langle x,y,z\rangle}{(xy-yx-2y,xz-zx-2z,yz-zy-x)}
\end{equation}
be the \textit{universal enveloping algebra}\footnote{The symbol $\mathbb{C}\langle x,y,z\rangle$ is the free algebra, i.e. the set of all words obtained by sequences of the symbols $x,y,z$ and their sums} of the Lie algebra $\mathfrak{sl}_{2}\mathbb{C}=\{A\in M_{2,2}(\mathbb{C})\vert \text{Tr}(A)=0\}$ with Lie bracket $[A,B]=AB-BA$. The elements in $\ccalA$ can be written as $a=\sum_{m,i,j,k}\alpha_{m,i,j,k}\boldsymbol{\pi}(x^{i},y^{j},z^{k})$ where $\alpha_{m,i,j,k}\in\mathbb{C}$ and $\boldsymbol{\pi}(x^{i},y^{j},z^{k})$ represents a permutation of the terms in $x^{i}y^{j}z^{k}$. The homomorphism $\rho$ is defined according to $\rho(x)=X$, $\rho(y)=Y$, $\rho(z)=Z$ where $X,Y,Z\in M_{N,N}(\mathbb{C})$, $XY-YX=2Y$, $XZ-ZX=2Z$ and $YZ-ZY=X$. As indicated in~\cite{hall_liealg,erdmann2006introduction} these sets of equations lead to a unique solution. Then, the filtering in~\eqref{eqn_ASP_filter_outputs} is of the form
\begin{multline}\label{eqn_liesp_filter_outputs}
\rho\left(\sum_{m,i,j,k}\alpha_{m,i,j,k}\boldsymbol{\pi}\left(x^{i},y^{j},z^{k}\right)\right)\bbx\\
=\sum_{m,i,j,k}\alpha_{m,i,j,k}\boldsymbol{\pi}\left(X^{i},Y^{j},Z^{k}\right)\bbx .
\end{multline}
These filters define convolutional signal processing in Lie Algebras.
\end{example}

The choice of $\ccalA$ and $\rho$ provides means to leverage our knowledge of the signal's domain in its processing. The convolutional filters in \eqref{eqn_DTSP_filter_outputs} leverage the shift invariance of time signals and the filters in \eqref{eqn_DSP_filter_outputs} the cyclic invariance of periodic signals. The group convolutional filters in \eqref{eqn_group_convolution} generalize shift invariance to invariance with respect to the an arbitrary group. The graph convolutional filters in \eqref{eqn_GSP_filter_outputs} engender signal processing that is independent of node labeling \cite{gamagnns} and the graphon filters in Example \ref{ex_wsp} a generalization of this notion to dense domains~\cite{ruiz2020graphon}. Leveraging this structure is instrumental in achieving scalable information processing. In the following section we explain how neural network architectures combine algebraic filters as defined in \eqref{eqn_ASP_filter_outputs} with pointwise nonlinearities to attain signal processing that inherits the invariance properties of the respective algebraic filters.

\begin{figure}
\centering\input{../figures/gnn_block_diagram_single_feature_uncolored.tex}
  \caption{Algebraic Neural Network $\Xi=\{(\ccalA_{\ell},\ccalM_{\ell},\rho_{\ell},\sigma_{\ell})\}_{\ell=1}^{3}$ with three layers indicating how the mapping $\mathbf{\Phi}_{\Xi}$ acts over an input signal $\bbx$.}
  \label{fig_6}
\end{figure}

\section{Algebraic Neural Networks}\label{sec_Algebraic_NNs}

With the concept of algebraic filtering at hand we define an algebraic neural network (AlgNN) as a stacked layered structure (see Fig.~\ref{fig_6})
 in which each layer is composed by the quadruple $(\ccalA_{\ell},\ccalM_{\ell},\rho_{\ell},\sigma_{\ell})$, where
\begin{itemize}
\item $(\ccalA_{\ell},\ccalM_{\ell},\rho_{\ell},\cdot)$ represents an algebraic signal model associated to each layer. Notice that $(\ccalM_{\ell},\rho_{\ell})$ is a representation of $\ccalA_{\ell}$.
\item $\sigma_{\ell}$ is and operator that performs that operations of pointwise nonlinearity and pooling. $\sigma_{\ell}$ maps signals from
$\ccalM_{\ell}$ into $\ccalM_{\ell+1}$.
\end{itemize}
Then, the ouput from the layer $\ell$ in the AlgNN is given by
\begin{equation}
\bbx_{\ell}=\sigma_{\ell}\left(\rho_{\ell}(a_{\ell})\bbx_{\ell-1}\right)
\end{equation}
where $a_{\ell}\in\ccalA_{\ell}$. The AlgNNs will be typically represented by the symbol $\Xi=\{(\ccalA_{\ell},\ccalM_{\ell},\rho_{\ell},\sigma_{\ell})\}_{\ell=1}^{L}$ and the mapping associated to $\Xi$ by $\boldsymbol{\Phi}_{\Xi}$.
\subsubsection{Convolutional Features}
The processing in each layer can be performed by means of several families of filters, which will lead to several \textit{features}. In particular the feature $f$ obtained in the layer $\ell$ is given by
\begin{equation}\label{eqn:alg_feat}
\bbx_{\ell}^{f}=\sigma_{\ell}\left(\sum_{g=1}^{F_{\ell}}\rho_{\ell}\left(a_{\ell}^{gf}\right)\bbx_{\ell-1}^{g}\right),
\end{equation}
where $a_{\ell}^{fg}$ is the filter in $\ccalA_{\ell}$ used to process the $g$-th feature $\bbx_{\ell-1}^{g}$ obtained from layer $\ell-1$ and $F_{\ell}$ is the number of features.
\subsubsection{Pooling}
As stated in~\cite{deeplearning_book} the pooling operation in CNNs helps to keep representations approximately invariant to small translations of an input signal, and also helps to to improve the computational efficiency. In this work this operation is attributed  to the operator $\sigma_{\ell}$. In particular, we consider $\sigma_{\ell}=P_{\ell}\circ\eta_{\ell}$ where $P_{\ell}$ is a pooling operator and $\eta_{\ell}$ is a pointwise nonlinearity. The only property assumed from $\sigma_{\ell}$ is to be Lipschitz and to have zero as a fixed point, i.e. $\sigma_{\ell}(0)=0$. It is important to point out that $P_{\ell}$ is a \textit{projection operator} that maps elements from a given vector space into another.

\begin{example}[Traditional CNNs]\normalfont
Traditional CNNs rely on the use of typical signal processing models and can be considered a particular case of an AlgNN  where the algebraic signal model is the same as in example~\ref{ex_DSP}. Consequently,
the $f$th feature in layer $\ell$ is given by
\begin{equation}
\bbx_{\ell}^{f}=\sigma_{\ell}\left(\sum_{g=1}^{F_{\ell}}\sum_{i=1}^{K}h_{i}^{gf}\bbS_{\ell}^{i}\bbx_{\ell-1}^{g}\right),
\label{eq:featureCNN}
\end{equation}
where $\rho_{\ell}(t)=\bbS_{\ell}$. In this case $P_{\ell}$ is a sampling operator while typically $\eta_{\ell}(u)=\max\{0,u\}$.
\end{example}

\begin{example}[Graphon Neural Networks]\normalfont
In graphon neural networks the algebraic signal model in each layer corresponds to the one discussed in example~\ref{ex_wsp}. Consequently, the $f$th feature in layer $\ell$ has the form given in~\eqref{eq:featureCNN}. It is worth noticing that in this case the signals considered are meant to be represented using a finite number of eigenvectors of the shift operator. Then, if we denote by $\boldsymbol{\varphi}_{i}^{(\ell)}$ the eigenvectors of the graphon $W_{\ell}(u,v)$ associated to the $\ell$th layer $P_{\ell}$ is given by
\begin{equation}
P_{\ell}:\text{span}\left\lbrace \boldsymbol{\varphi}_{i}^{(\ell)}\right\rbrace_{i\in I}\rightarrow\text{span}\left\lbrace\boldsymbol{\varphi}_{j}^{(\ell+1)}\right\rbrace_{j\in J},
\end{equation}
where $I$ and $J$ are finite subsets of $\mathbb{N}$. Then, the signals in each $\ccalM_{\ell}$ are always functions of a continuous independent variable but with \textit{lower degrees of freedom} as $\ell$ increases.
\end{example}

\begin{example}[Group Neural Networks]\normalfont
In group neural networks the algebraic model is the same as specified in example~\ref{ex_groupsp}. Therefore, the $f$th feature in layer $\ell$ is given by
\begin{equation}
\bbx_{\ell}^{f}=\sigma_{\ell}\left(\sum_{n=1}^{F_{\ell}}\sum_{u,h\in G_{\ell}}\bba^{nf}_{(\ell)}(uh^{-1})\bbx_{\ell-1}(h)u\right).
\label{eq:featureGroupNN}
\end{equation}
where $G_{\ell}$ is the group associated to the $\ell$th layer and $\bba^{nf}_{(\ell)}$ are the coefficients of the filter associated to the feature $f$ in layer $\ell$. In this case $P_{\ell}:L^{2}(G_{\ell})\rightarrow L^{2}(G_{\ell+1})$,  where $L^{2}(G)$ is the set of signals of finite energy defined on the group $G$. If the groups $G_{\ell}$ are finite $P_{\ell}$ can be conceived as a typical projection mapping between $\mathbb{R}^{\vert G_{\ell}\vert}\rightarrow\mathbb{R}^{\vert G_{\ell+1}\vert}$.
\end{example}

\section{Algebraic Filters} \label{sec_alg_filters

Algebraic signal processing (ASP) provides a framework for understanding and generalizing traditional signal processing exploiting the representation theory of algebras~\cite{algSP0, algSP1, algSP2, algSP3}. In ASP, a signal model is defined as the triple
\begin{equation}\label{eqn_ASP_signal_model}
   (\ccalA,\ccalM,\rho),
\end{equation}
in which $\ccalA$ is an associative algebra with unity, $\ccalM$ is a vector space with inner product, and $\rho:\ccalA\rightarrow\text{End}(\ccalM)$ is a homomorphism between the algebra $\ccalA$ and the set of endomorphisms in the vector space $\ccalM$. The elements in \eqref{eqn_ASP_signal_model} are tied together by the notion of a representation which we formally define next. In this paper we further restrict attention to commutative algebras $\ccalA$.

\begin{definition}[Representation] \label{def_representation} A representation of the associative algebra $\ccalA$ is a vector space $\ccalM$ equipped with a homomorphism $\rho:\ccalA\rightarrow\text{End}(\ccalM)$, i.e., a linear map preserving multiplication and unit.
\end{definition}
When denoting a representation of a given algebra we will use a pair consisting of the vector space and the homomorphism, for instance with the notation in definition~\ref{def_representation} the representation of the algebra $\mathcal{A}$ is denoted by $(\mathcal{M},\rho)$. The class of representations of  an algebra $\mathcal{A}$ is denoted as $\mathsf{Rep}\{\mathcal{A}\}$.

In an ASP model, signals  are elements of the vector space $\ccalM$, and filters are elements of the algebra $\ccalA$. Thus, the vector space $\ccalM$ deteremines the objects of interest and the algebra $\ccalA$ the rules of the operations that define a filter. The homomorphism $\rho$ translates the abstract operators $a\in\ccalA$ into concrete operators $\rho(a)$ that act on signals $\bbx$ to produce filter outputs
\begin{equation}\label{eqn_ASP_filter_outputs}
   \bby = \rho(a) \bbx.
\end{equation}
The algebraic filters in \eqref{eqn_ASP_filter_outputs} generalize the convolutional processing of time signals -- see Example \ref{ex_DTSP}. Our goal in this paper is to use them to generalize convolutional neural networks (Section \ref{sec_Algebraic_NNs}) and to study their fundamental stability properties (Section \ref{sec:perturbandstability}). Generators, which we formally define next, are important for the latter goal.

\begin{definition}[Generators] \label{def_generators} For an associative algebra with unity $\ccalA$ we say the set $\ccalG\subseteq\ccalA$ generates $\ccalA$ if all $a\in\ccalA$ can be represented as polynomial functions of the elements of $\ccalG$. We say elements $g\in\ccalG$ are generators of $\ccalA$ and we denote as $a=p_{\ccalA}(\ccalG)$ the polynomial that generates $a$.
\end{definition}
Definition \ref{def_generators} states that elements $a\in\ccalA$ can be built from the generating set as polynomials using the operations of the algebra. Given the role of representations in connecting the algebra $\ccalA$ with the signals $\bbx$ as per Definition \ref{def_representation} the representation $\rho(g)$ of a generator $g\in\ccalG$ will be of interest. In the context of ASP, these representations are called shift operators as we formally define next.

\begin{definition}[Shift Operators]\label{def_shift_operators}
Let $(\ccalM,\rho)$ be a representation of the algebra $\ccalA$. Then, if $\ccalG\subset\ccalA$ is a generator set of $\ccalA$, the operators $\bbS = \rho(g)$ with $g\in\ccalG$ are called shift operators. The set of all shift operators is denoted by $\ccalS$
\end{definition}

Given that elements $a$ of the algebra are generated from elements $g$ of the generating set, it follows that filters $\rho(a)$ are generated from the set of shift operators $\bbS=\rho(g)$. In fact, if we have that $a=p_{\ccalA}(\ccalG)$ is the polynomial that generates $a$, the filter $\rho(a)$ is generated by
\begin{equation}\label{eqn_filter_representation}
   \rho(a) = p_\ccalM \big(\rho(\ccalG)\big)
           = p_\ccalM\big(\ccalS\big)
           = p\big(\ccalS\big),
\end{equation}
where the subindex $\ccalM$ signifies that the operations in \eqref{eqn_filter_representation} are those of the vector space $\ccalM$ -- in contrast to the polynomial $a=p_{\ccalA}(\ccalG)$ whose operations are those of the algebra $\ccalA$. In the last equality we drop the subindex $\ccalM$ to simplify notation as it is generally understood from context.

To clarify ideas it is instructive to consider some examples.

\newcounter{dtsp}
\setcounter{dtsp}{\theexample}
\setcounter{example}{\thedtsp}

\begin{example}[Discrete Time Signal Processing]\label{ex_DTSP}\normalfont Let $\ccalM = \ell^2$ be the space of square summable sequences $\bbx = \{x_n\}_{n\in\mbZ}$ and $\ccalA$ the algebra of polynomials with elements $a = \sum_{k=0}^{K-1} h_k t^k$. We consider the shift operator $S$ such that $S\bbx$ is the sequence with entries $(S\bbx)_n = x_{n-1}$. With the homomorphism $\rho(a) = \sum_{k=0}^{K-1} h_k S^k$ the filters in \eqref{eqn_ASP_filter_outputs} are of the form
\begin{equation}\label{eqn_DTSP_filter_outputs}
   \bby  \,=\, \rho\left(\sum_{k=0}^{K-1} h_k t^k\right) \bbx
         \,=\, \left(\sum_{k=0}^{K-1} h_k S^k\right) \bbx
         \,=\, \sum_{k=0}^{K-1} h_k S^k \bbx.
\end{equation}
This is a representation of time convolutional filters \cite{algSP0}. Observe that the algebra of polynomials is generated by $g=t$ and that since $\rho(t) = S$ the shift operator $S$ generates the set of convolutional filters as it follows from \eqref{eqn_DTSP_filter_outputs}.\end{example}

\newcounter{gsp}
\setcounter{gsp}{\theexample}
\setcounter{example}{\thegsp}

\begin{example}[Graph Signal Processing]\label{ex_GSP}\normalfont We retain the algebra of polynomials as in Example \ref{ex_DTSP} but we change the space of signals to the set of complex vectors with $N$ entries, $\ccalM =\mbC^N$. We interpret components $x_n$ of $\bbx\in\ccalM = \mbC^N$ as being associated with nodes of a graph with matrix representation $\bbS\in\mbC^{N\times N}$. With the homomorphism $\rho(a) = \sum_{k=0}^{K-1} h_k \bbS^k$ the filters in~\eqref{eqn_ASP_filter_outputs} are of the form
\begin{equation}\label{eqn_GSP_filter_outputs}
   \bby \,=\, \rho\left(\sum_{k=0}^{K-1} h_k t^k\right) \bbx
        \,=\, \left(\sum_{k=0}^{K-1} h_k \bbS^k\right) \bbx
        \,=\, \sum_{k=0}^{K-1} h_k \bbS^k \bbx.
\end{equation}
This is a representation of the graph convolutional filters used in graph signal processing (GSP)~\cite{gsp_graphfilters,ortega_gsp}. Observe that \eqref{eqn_GSP_filter_outputs} and \eqref{eqn_DTSP_filter_outputs} are similar but represent different operations. In \eqref{eqn_GSP_filter_outputs} $\bbx$ is a vector and $\bbS^k$ a matrix power. In \eqref{eqn_DTSP_filter_outputs} $\bbx$ is a sequence and $S^k$ is the composition of the shift operator $S$. Their similarity arises from their common use of the algebra of polynomials. Their differences are because they use different vector spaces $\ccalM$ and different homomorphism functions $\rho$. \end{example}

\newcounter{dsp}
\setcounter{dsp}{\theexample}
\setcounter{example}{\thedsp}
\begin{example}[Discrete Signal Processing]\label{ex_DSP}\normalfont Make $\ccalM=\mbC^N$ and $\ccalA$ the algebra of polynomials in the variable $t$ and modulo $t^{N}-1$. Elements of $\ccalA$ are of the form $a = \sum_{k=0}^{K-1} h_k t^k$ but we must have $K\leq N$ and polynomial products use the rule $t^k = t^{k \mod N}$. Define the cyclic shift operator $C$ as one whose application to $\bbx\in\ccalM=\mbC^N$ yields a vector with components
$(C \bbx)_n = x_{(n-1)\mod N}$. With the homomorphism $\rho(a) = \sum_{k=0}^{K-1} h_k C^k$ the filters in~\eqref{eqn_ASP_filter_outputs} are of the form
\begin{equation}\label{eqn_DSP_filter_outputs}
   \bby \,=\, \rho\left(\sum_{k=0}^{K-1} h_k t^k\right) \bbx
        \,=\, \left(\sum_{k=0}^{K-1} h_k C^k\right) \bbx
        \,=\, \sum_{k=0}^{K-1} h_k C^k \bbx.
\end{equation}
This is the representation of (cyclic) convolutional filters in discrete signal processing~\cite{algSP0}. We emphasize that the homomorphism $\rho(a) = \sum_{k=0}^{K-1} h_k C^k$ is indeed a homomorphism because the cyclic shift operator $C$ satisfies $C^k = C^{k \mod N}$. \end{example}

As is clear from Examples \ref{ex_DTSP}-\ref{ex_DSP}, the effect of the operator $\rho(a)$ on a given signal $\bbx$ is determined by two factors: The filter $a\in\ccalA$ and the homomorphism $\rho$. The filter $a\in\ccalA$ indicates the \textit{laws and rules} to be used to manipulate the signal $\bbx$ and $\rho$ provides a \textit{physical realization} of the filter $a$ on the space $\ccalM$ to which $\bbx$ belongs. For instance, in these three examples the filter $a=1 + 2t$ indicates that the signal is to be added to a transformed version of the signal scaled by coefficient $2$. The homomorphism $\rho$ in Example \ref{ex_DTSP} dictates that the physical implementation of this transformation is a time shift. The homomorphism $\rho$ in Example \ref{ex_GSP} defines a transformation as a multiplication by $\bbS$ and in Example \ref{ex_DSP} the homomorphism entails a cyclic shift. We remark that in order to specify the physical effect of a filter it is always sufficient to specify the physical effect of the generators. In all three examples, the generator of the algebra is $g=t$. The respective effects of an arbitrary filter $a$ are determined once we specify that $\rho(t) = S$ in Example \ref{ex_DTSP}, $\rho(t) = \bbS$ in Example \ref{ex_GSP}, or $\rho(t) = C$ in Example \ref{ex_DSP}.

The flexibility in the choice of algebra and homomorphism allows for a very rich variety of signal processing frameworks. We highlight this richness with three more examples.

\newcounter{groupsp}
\setcounter{groupsp}{\theexample}
\setcounter{example}{\thegroupsp}

\begin{example}[Signal Processing on Groups]\normalfont \label{ex_groupsp} Let $\ccalM=\{\bbx: G\rightarrow\mathbb{C} \}=\{ \sum_{g\in G}\bbx(g)g\}$ be the set of functions defined on the group $G$ with values in $\mathbb{C}$ and $\ccalA=\ccalM$ the \textit{group algebra} with elements $\sum_{g\in G}\bba(g)g$. With a homomorphism given by $\rho(\bba)=\bba$ the filtering in~\eqref{eqn_ASP_filter_outputs} takes the form
\begin{equation}\label{eqn_group_filters}
\rho\left(\sum_{g\in G}\bba(g)g\right)\bbx=\sum_{g\in G}\bba(g)g\bbx=\sum_{g\in G}\sum_{h\in G}\bba(g)\bbx(h)gh,
\end{equation}
and making $u=gh$ we have
\begin{equation}\label{eqn_group_convolution}
   \sum_{g,h\in G} \bba(g)\bbx(h)gh
             =\sum_{u,h\in G}\bba(uh^{-1})\bbx(h)u=\bba\ast\bbx .
\end{equation}
This is the standard representation of convolution of signals on groups~\cite{steinbergrepg,terrasFG,fulton1991representation}. We point out that \eqref{eqn_group_filters} and \eqref{eqn_group_convolution} hold for any group but that not all group algebras are commutative. Results in Section \ref{sec:stabilitytheorems} apply only when the group group algebra is commutative.
\end{example}

\newcounter{wsp}
\setcounter{wsp}{\theexample}
\setcounter{example}{\thewsp}

\begin{example}[Graphon Signal Processing ]\normalfont \label{ex_wsp} Let $\ccalM=L^2([0,1])$ be the space of signals and $\ccalA$ the algebra of polynomials with elements $a = \sum_{k=0}^{K-1} h_k t^k$ with a homomorphism given by
$\rho(a)=\sum_{k=0}^{K-1} h_k S^k$ where
\begin{equation}\label{eqn_graphon_shift}
   \left(S\bbx\right)(u)=\int_{0}^{1}W(u,v)\bbx(v)dv,
\end{equation}
where $W(u,v):[0,1]^{2}\rightarrow [0,1]$ is a graphon, i.e. bounded symmetric measurable function~\cite{lovaz2012large,segarragraphon}. Then, the filtering in~\eqref{eqn_ASP_filter_outputs} is of the form~\eqref{eqn_DTSP_filter_outputs}.
This is a representation of convolutional filters in the context of graphon signal processing~\cite{ruiz2020graphon}. \end{example}

The choice of $\ccalA$ and $\rho$ provides means to leverage our knowledge of the signal's domain in its processing. The convolutional filters in \eqref{eqn_DTSP_filter_outputs} leverage the shift invariance of time signals and the filters in \eqref{eqn_DSP_filter_outputs} the cyclic invariance of periodic signals. The group convolutional filters in \eqref{eqn_group_convolution} generalize shift invariance to invariance with respect to an arbitrary group action. The graph convolutional filters in \eqref{eqn_GSP_filter_outputs} engender signal processing that is independent of node labeling \cite{gamagnns} and the graphon filters in Example \ref{ex_wsp} a generalization of this notion to dense domains~\cite{ruiz2020graphon}. Leveraging this structure is instrumental in achieving scalable information processing. In the following section we explain how neural network architectures combine algebraic filters as defined in \eqref{eqn_ASP_filter_outputs} with pointwise nonlinearities to attain signal processing that inherits the invariance properties of the respective algebraic filters.
\begin{figure}
\centering\input{./figures/gnn_block_diagram_single_feature_uncolored_01.tex}
  \caption{Algebraic Neural Network $\Xi=\{(\ccalA_{\ell},\ccalM_{\ell},\rho_{\ell})\}_{\ell=1}^{3}$ with three layers indicating how the input signal $\mathbf{x}$ is processed by $\Xi$ and mapped into $\mathbf{x}_{3}$.}
  \label{fig_6}
\end{figure}
\section{Algebraic Neural Networks}\label{sec_Algebraic_NNs}
With the concept of algebraic filtering at hand we define an algebraic neural network (AlgNN) as a stacked layered structure (see Fig.~\ref{fig_6})
 in which each layer is composed by the triple $(\ccalA_{\ell},\ccalM_{\ell},\rho_{\ell})$, which is an algebraic signal model associated to each layer. Notice that $(\ccalM_{\ell},\rho_{\ell})$ is a representation of $\ccalA_{\ell}$. The mapping between layers is performed by the maps $\sigma_{\ell}:\mathcal{M}_{\ell}\rightarrow\mathcal{M}_{\ell+1}$ that perform those operations of point-wise nonlinearity and pooling. Then, the ouput from the layer $\ell$ in the AlgNN is given by
\begin{equation}\label{eq:xl}
\bbx_{\ell}=\sigma_{\ell}\left(\rho_{\ell}(a_{\ell})\bbx_{\ell-1}\right)
\end{equation}
where $a_{\ell}\in\ccalA_{\ell}$, which can be represented equivalently as
\begin{equation}\label{eq:interlayeropalgnn}
\mathbf{x}_{\ell}=\Phi (\mathbf{x}_{\ell-1},\mathcal{P}_{\ell-1},\mathcal{S}_{\ell-1}),
\end{equation}
where $\mathcal{P}_{\ell}\subset\mathcal{A}_{\ell}$ highlights the properties of the filters and $\mathcal{S}_{\ell}$ is the set of shifts associated to $(\mathcal{M}_{\ell},\rho_{\ell})$. Additionally, the term $\Phi\left(\mathbf{x},\{ \mathcal{P}_{\ell} \}_{1}^{L},\{ \mathcal{S}_{\ell}\}_{1}^{L}\right)$ represents the total map associated to an AlgNN acting on a signal $\mathbf{x}$.
\subsubsection{Convolutional Features}
The processing in each layer can be performed by means of several families of filters, which will lead to several \textit{features}. In particular the feature $f$ obtained in the layer $\ell$ is given by
\begin{equation}\label{eqn:alg_feat}
\bbx_{\ell}^{f}=\sigma_{\ell}\left(\sum_{g=1}^{F_{\ell}}\rho_{\ell}\left(a_{\ell}^{gf}\right)\bbx_{\ell-1}^{g}\right),
\end{equation}
where $a_{\ell}^{fg}$ is the filter in $\ccalA_{\ell}$ used to process the $g$-th feature $\bbx_{\ell-1}^{g}$ obtained from layer $\ell-1$ and $F_{\ell}$ is the number of features.
\subsubsection{Pooling}
As stated in~\cite{deeplearning_book} the pooling operation in CNNs helps to keep representations approximately invariant to small translations of an input signal, and also helps to improve the computational efficiency. In this work this operation is attributed  to the operator $\sigma_{\ell}$. In particular, we consider $\sigma_{\ell}=P_{\ell}\circ\eta_{\ell}$ where $P_{\ell}$ is a pooling operator and $\eta_{\ell}$ is a pointwise nonlinearity. The only property assumed from $\sigma_{\ell}$ is to be Lipschitz and to have zero as a fixed point, i.e. $\sigma_{\ell}(0)=0$. It is important to point out that $P_{\ell}$ projects elements from a given vector space into another.

\begin{example}[Traditional CNNs]\normalfont
Traditional CNNs rely on the use of typical signal processing models and can be considered a particular case of an AlgNN  where the algebraic signal model is the same as in example~\ref{ex_DSP}. Consequently,
the $f$th feature in layer $\ell$ is given by
\begin{equation}
\bbx_{\ell}^{f}=\sigma_{\ell}\left(\sum_{g=1}^{F_{\ell}}\sum_{i=1}^{K}h_{i}^{gf}\bbS_{\ell}^{i}\bbx_{\ell-1}^{g}\right),
\label{eq:featureCNN}
\end{equation}
where $\rho_{\ell}(t)=\bbS_{\ell}$. In this case $P_{\ell}$ is a sampling operator while typically $\eta_{\ell}(u)=\max\{0,u\}$.
\end{example}

\begin{example}[Graphon Neural Networks]\normalfont
In graphon neural networks the algebraic signal model in each layer corresponds to the one discussed in example~\ref{ex_wsp}. Consequently, the $f$th feature in layer $\ell$ has the form given in eqn.~\eqref{eq:featureCNN} considering a shift $S$ given according to eqn.~(\ref{eqn_graphon_shift}). It is worth noticing that in this case the signals in the inner layers of the network are meant to be represented using a finite number of eigenvectors of the shift operator. Then, if we denote by $\boldsymbol{\varphi}_{i}^{(\ell)}$ the eigenvectors of the graphon $W_{\ell}(u,v)$ associated to the $\ell$th layer $P_{\ell}$ is given by
\begin{equation}
P_{\ell}:\text{span}\left\lbrace \boldsymbol{\varphi}_{i}^{(\ell)}\right\rbrace_{i\in I}\rightarrow\text{span}\left\lbrace\boldsymbol{\varphi}_{j}^{(\ell+1)}\right\rbrace_{j\in J},
\end{equation}
where $I$ and $J$ are finite subsets of $\mathbb{N}$. Then, the signals in each $\ccalM_{\ell}$ are always functions of a continuous independent variable but with \textit{lower degrees of freedom} as $\ell$ increases.
\end{example}

\begin{example}[Group Neural Networks]\normalfont
In group neural networks the algebraic model is the same as specified in example~\ref{ex_groupsp}. Therefore, the $f$th feature in layer $\ell$ is given by
\begin{equation}
\bbx_{\ell}^{f}=\sigma_{\ell}\left(\sum_{n=1}^{F_{\ell}}\sum_{u,h\in G_{\ell}}\bba^{nf}_{(\ell)}(uh^{-1})\bbx_{\ell-1}(h)u\right).
\label{eq:featureGroupNN}
\end{equation}
where $G_{\ell}$ is the group associated to the $\ell$th layer and $\bba^{nf}_{(\ell)}$ are the coefficients of the filter associated to the feature $f$ in layer $\ell$. In this case $P_{\ell}:L^{2}(G_{\ell})\rightarrow L^{2}(G_{\ell+1})$,  where $L^{2}(G)$ is the set of signals of finite energy defined on the group $G$. If the groups $G_{\ell}$ are finite $P_{\ell}$ can be conceived as a typical projection mapping between $\mathbb{R}^{\vert G_{\ell}\vert}\rightarrow\mathbb{R}^{\vert G_{\ell+1}\vert}$.
\end{example}

\section{Stability Properties in Algebraic Neural Networks}

\subsection{Frechet differentiability}
\subsection{Finite Dimensional Representations}
\subsubsection{Absolute Perturbations}
\subsubsection{Relative Perturbations}

\subsection{Infinite Dimensional Representations (countable)}

\subsection{Infinite Dimensional Representations (uncountable)}

\subsection{Tradeoff between Stability and Selectivity}
\section{Spectral Operators}\label{sec:spectopt}
Before we proof the theorems presented in previous section, it is important to discuss the concept of \textit{Fourier decomposition} in the context of representation theory of algebras. We want to highlight specially the role of the filters when a representation is compared to its perturbed version. In particular, we will show that there are essentially two factors that can cause  differences between operators and their perturbed versions, the eigenvalues\footnote{As we will show later, this is indeed a particular case of a general notion of homomorphism between the algebra and an irreducible subrepresentation of $\mathcal{A}$.} and eigenvectors. Additionally, we show how the algebra can only affect one of those sources. This is consistent with the fact that differences in the eigenvectors of the operators only affect the constants that are associated to the stability bounds.

We highlight that the concept of spectral decomposition exploited extensively in traditional signal processing can be considered as a particular case of the concept of direct sums\footnote{When considering representations of non-compact groups this concept is insufficient and it is replaced for the more general notion of \textit{direct integral of Hilbert spaces}~\cite{folland2016course,deitmar2014principles}} of \textit{irreducible} representations. In this section we provide some basic definitions from the representation theory of algebras and its connection with the notion of spectral decomposition of signals in algebraic neural networks.  This will be useful and necessary in Section~\ref{sec:proofofTheorems} as the norm of the operators considered can be calculated in the terms of an spectral representation. We start with the following definition.
\begin{definition}
Let $(\mathcal{M},\rho)$ be a representation of $\mathcal{A}$. Then, a representation $(\mathcal{U},\rho)$ of $\mathcal{A}$ is a \textbf{subrepresentation} of  $(\mathcal{M},\rho)$ if $\mathcal{U}\subseteq\mathcal{M}$ and $\mathcal{U}$ is invariant under all operators $\rho(a)~\forall~a\in\mathcal{A}$, i.e. $\rho(a)u\in\mathcal{U}$ for all $u\in\mathcal{U}$ and $a\in\mathcal{A}$. A representation $(\mathcal{M}\neq 0,\rho)$ is \textbf{irreducible} or simple if the only subrepresentations of $(\mathcal{M}\neq 0,\rho)$ are $(0,\rho)$ and $(\mathcal{M},\rho)$.
\end{definition}
The class of irreducible representations of an algebra $\mathcal{A}$ is denoted by $\text{Irr}\{\mathcal{A}\}$. Notice that the zero vector space and $\mathcal{M}$ induce themselves subrepresentations of $(\mathcal{M},\rho)$. In order to state a comparison between representations the concept of \textit{homomorphism between representations} is introduced in the following definition.
\begin{definition}\label{def:homrep}
Let $(\mathcal{M}_{1},\rho_{1})$ and $(\mathcal{M}_{2},\rho_{2})$ be two representations of an algebra $\mathcal{A}$. A homomorphism or \textbf{interwining operator} $\phi:\mathcal{M}_{1}\rightarrow\mathcal{M}_{2}$ is a linear operator which commutes with the action of $\mathcal{A}$, i.e.
\begin{equation}
\phi(\rho_{1}(a)v)=\rho_{2}(a)\phi(v).
\label{eq:interwinop}
\end{equation}
A homomorphism $\phi$ is said to be an isomorphism of representations if it is an isomorphism of vectors spaces.
\end{definition}
Notice from definition~\ref{def:homrep}  a substantial difference between the concepts of isomorphism of vector spaces and  isomorphism of representations. In the first case we can consider that two arbitray vector spaces of the same dimension (finite) are isomorphic, while for representations that condition is required but still the condition in eqn.~(\ref{eq:interwinop}) must be satisfied. For instance, as pointed out in~\cite{repthysmbook} all the irreducible 1-dimensional representations of the polynomial algebra $\mathbb{C}[t]$ are non isomorphic.

As we have discussed before, the vector space $\mathcal{M}$ associated to $(\mathcal{M},\rho)$  provides the space where the signals are modeled. Therefore, it is of central interest to determine whether it is possible or not to \textit{decompose} $\mathcal{M}$ in terms of simpler or smaller spaces consistent with the action of $\rho$.  We remark that for any two representations $(\mathcal{M}_{1},\rho_{1})$ and $(\mathcal{M}_{2},\rho_{2})$ of an algebra $\mathcal{A}$, their direct sum is given by the representation $(\mathcal{M}_{1}\oplus\mathcal{M}_{2},\rho)$ where $\rho(a)(\mathbf{x}_{1}\oplus\mathbf{x}_{2})=(\rho_{1}(a)\mathbf{x}_{1}\oplus\rho_{2}(a)\mathbf{x}_{2})$. We introduce the concept of indecomposability in the following definition.
\begin{definition}\label{def:indecomprep}
A nonzero representation $(\mathcal{M},\rho)$ of an algebra $\mathcal{A}$ is said to be \textbf{indecomposable} if it is not isomorphic to a direct sum of two nonzero representations.
\end{definition}
Indecomposable representations provide the \textit{minimum units of information} that can be extracted from signals in a given space when the filters have a specific structure (defined by the algebra)~\cite{barot2014introduction}. The following theorem provides the basic building block for the decomposition of finite dimensional representations.
\begin{theorem}[Krull-Schmit, \cite{repthybigbook}]\label{theorem:krullschmitdecomp}
Any finite dimensional representation of an algebra can be decomposed into a finite direct sum of indecomposable subrepresentations and this decomposition is unique up to the order of the summands and up to isomorphism.
\end{theorem}
The uniqueness in this result means that if $(\oplus_{r=1}^{r}V_{i},\rho)\cong(\oplus_{j=1}^{s}W_{j},\gamma)$ for indecomposable representations $(V_{j},\rho_{j}),(W_{j},\gamma_{j})$, then $r=s$ and there is a permutation $\pi$ of the indices such that $(V_{i},\rho_{i})\cong (W_{\pi(j)},\gamma_{\pi(j)})$~\cite{repthybigbook}. Although theorem~\ref{theorem:krullschmitdecomp} provides the guarantees for the decomposition of representation in terms of indecomposable representations, it is not applicable when infinite dimensional representations are considered. However, it is possible to overcome this obstacle taking into account that irreducible representations are indecomposable~\cite{repthysmbook,repthybigbook}, and they can be used then to build representations that are indecomposable. In particular, \textit{irreducibility} plays a central role to decompose the invariance properties of the images of $\rho$ on $\text{End}(\mathcal{M})$~\cite{repthybigbook}. Representations that allow a decomposition in terms of subrepresentations that are irreducible are called \textit{completely reducible} and its formal description is presented in the following definition.
\begin{definition}[~\cite{repthybigbook}]\label{def:complredu}
 A representation $(\mathcal{M},\rho)$ of the algebra $\mathcal{A}$ is said to be \textbf{completely reducible} if $(\mathcal{M},\rho)=\bigoplus_{i\in I}(\mathcal{U}_{i},\rho_{i})$ with irreducible subrepresentations $(\mathcal{U}_{i},\rho_{i})$. The \textbf{length} of $(\mathcal{M},\rho)$ is given by $\vert I\vert$.
\end{definition}

For a given $(\mathcal{U},\rho_{\mathcal{U}})\in\text{Irr}\{\mathcal{A}\}$ the sum of all irreducible subrepresentations of $(V,\rho_{V})$ that are equivalent (isomorphic) to $(\mathcal{U},\rho_{\mathcal{U}})$ is represented by $V(\mathcal{U})$ and it is called the $\mathcal{U}$-homogeneous component of $(V,\rho_{V})$. This sum is a direct sum, therefore it has a length that is well defined and whose value is called the \textit{multiplicity} of $(\mathcal{U},\rho_{\mathcal{U}})$ and is represented by $m(\mathcal{U},V)$~\cite{repthybigbook}. Additionally, the sum of all irreducible subrepresentations of $(V,\rho_{V})$ will be denoted as $\text{soc}\{V\}$. It is possible to see that a given representation $(V,\rho_{V})$ is completely reducible if and only if $(V,\rho_{V})=\text{soc}\{S\}$~\cite{repthybigbook}. The connection between $\text{soc}\{V\}$ and $V(\mathcal{U})$ is given by the following proposition.
\begin{proposition}[Proposition 1.31~\cite{repthybigbook}]\label{prop:repdirectsum1}
Let $(V,\rho_{V})\in\mathsf{Rep}\{A\}$. Then $\text{soc}\{V\}=\bigoplus_{S\in\text{Irr}\{\mathcal{A}\}}V(S)$.
\end{proposition}
Now, taking into account that any homogeneous component $V(\mathcal{U})$ is itself a direct sum we have that
\begin{equation}
\text{soc}\{V\}\cong\bigoplus_{S\in\text{Irr}\{\mathcal{A}\}}S^{\oplus m(\mathcal{U},V)}.
\label{eq:repdirectsumtot}
\end{equation}
Equation~(\ref{eq:repdirectsumtot}) provides the building block for the definition of Fourier decompositions in algebraic signal processing~\cite{algSP1}. With all these concepts at hand we are ready to introduce the following definition.
\begin{definition}[Fourier Decomposition]\label{def:foudecomp}
For an algebraic signal model $(\mathcal{A},\mathcal{M},\rho)$ we say that there is a spectral or Fourier decomposition of $(\mathcal{M},\rho)$ if
\begin{equation}
(\mathcal{M},\rho)\cong\bigoplus_{(\mathcal{U}_{i},\phi_{i})\in\text{Irr}\{\mathcal{A}\}}(\mathcal{U}_{i},\phi_{i})^{\oplus m(\mathcal{U}_{i},\mathcal{M})}
 \label{eq:foudecomp1}
\end{equation}
where the $(\mathcal{U}_{i},\phi_{i})$ are irreducible subrepresentations of $(\mathcal{M},\rho)$. Any signal $\mathbf{x}\in\mathcal{M}$ can be therefore represented by the map $\Delta$ given by
\begin{equation}
\begin{matrix}
\Delta: & \mathcal{M} \to \bigoplus_{(\mathcal{U}_{i},\phi_{i})\in\text{Irr}\{\mathcal{A}\}}(\mathcal{U}_{i},\phi_{i})^{\oplus m(\mathcal{U}_{i},\mathcal{M})} \\ & \mathbf{x}\mapsto \hat{\mathbf{x}}
 \end{matrix}
 \label{eq:foudecomp2}
\end{equation}
known as the Fourier decomposition of $\mathbf{x}$ and the projection of $\hat{\mathbf{x}}$ in each $\mathcal{U}_{i}$ are the Fourier components represented by $\hat{\mathbf{x}}(i)$.
\end{definition}
Notice that in eqn.~(\ref{eq:foudecomp1}) there are two sums, one dedicated to the non isomorphic subrepresentations (external) and another one (internal) dedicated to subrepresentations that are isomorphic. In this context, the sum for non isomorphic representations indicates the sum on the \textit{frequencies} of the representation while the sum for isomorphic representations a sum of components associated to a given frequency. It is also worth pointing out that $\Delta$ is an interwining operator, therefore, we have that $\Delta(\rho(a)\mathbf{x})=\rho(a)\Delta(\mathbf{x})$. As pointed out in~\cite{algSP0} this can be used to define a convolution operator as $\rho(a)\mathbf{x}=\Delta^{-1}(\rho(a)\Delta(\mathbf{x}))$. The projection of a filtered signal $\rho(a)\mathbf{x}$ on each $\mathcal{U}_{i}$ is given by $\phi_{i}(a)\hat{\mathbf{x}}(i)$ and the collection of all this projections is known as the \textit{spectral representation} of the operator $\rho(a)$. Notice that $\phi_{i}(a)\hat{\mathbf{x}}(i)$ translates to different operations depending on the dimension of $\mathcal{U}_{i}$. For instance, if $\text{dim}(\mathcal{U}_{i})=1$, $\hat{\mathbf{x}}(i)$ and $\phi_{i}(a)$ are scalars while if $\text{dim}(\mathcal{U}_{i})>1$ and finite $\phi_{i}(a)\hat{\mathbf{x}}(i)$ is obtained as a matrix product.

\begin{remark}\label{remark:interpSpecFilt}\normalfont
The spectral representation of an operator indicated as $\phi_{i}(a)\hat{\mathbf{x}}(i)$ and eqns.~(\ref{eq:foudecomp1}) and~(\ref{eq:foudecomp2}) highlight one important fact that is essential for the discussion of the results in Section~\ref{sec:proofofTheorems}. For a completely reducible representation $(\mathcal{M},\rho)\in\mathsf{Rep}\{\mathcal{A}\}$ the connection between the algebra $\mathcal{A}$ and the spectral representation is \textit{exclusively} given by $\phi_{i}(a)$ which is acting on $\hat{\mathbf{x}}(i)$, therefore, it is not possible by the selection of elements or subsets of the algebra to do any modification on the spaces $\mathcal{U}_{i}$ associated to the irreducible components in eqn.(\ref{eq:foudecomp1}). As a consequence, when measuring the similarities between two operators $\rho(a)$ and $\tilde{\rho}(a)$ associated to $(\mathcal{M},\rho)$ and $(\mathcal{M},\tilde{\rho})$, respectively, there will be two sources of error. One source of error that can be modified by the selection of $a\in\mathcal{A}$ and another one that will be associated with the differences between spaces $\mathcal{U}_{i}$ and $\tilde{\mathcal{U}}_{i}$, which are associated to the direct sum decomposition of $(\mathcal{M},\rho)$ and $(\mathcal{M},\tilde{\rho})$, respectively. This point was first elucidated in~\cite{fern2019stability} for the particular case of GNNs, but it is part of a much more general statement that becomes more clear in the language of algebraic signal processing.
\end{remark}
\begin{example}[Discrete signal processing]\normalfont
In CNNs the filtering is defined by the polynomial algebra $\mathcal{A}=\mathbb{C}[t]/(t^{N}-1)$, therefore, in a given layer the spectral representation of the filters is given by
\begin{multline*}
\rho(a)\mathbf{x}=\sum_{i=1}^{N}\phi_{i}\left(\sum_{k=0}^{K-1}h_{k}t^{k}\right)\hat{\mathbf{x}}(i)\mathbf{u}_{i}\\
=\sum_{i=1}^{N}\sum_{k=0}^{K-1}h_{k}\phi_{i}(t)^{k}\hat{\mathbf{x}}(i)\mathbf{u}_{i}
=\sum_{i=1}^{N}\sum_{k=0}^{K-1}h_{k}\left(e^{-\frac{2\pi ij}{N}}\right)^{k}\hat{\mathbf{x}}(i)\mathbf{u}_{i},
\end{multline*}
with $a=\sum_{k=0}^{K-1}h_{k}t^{k}$ and where the $\mathbf{u}_{i}(v)=\frac{1}{\sqrt{N}}e^{\frac{2\pi jvi}{N}}$ are the column vectors of the traditional DFT matrix, while $\phi_{i}(t)=e^{-\frac{2\pi ji}{N}}$ is the eigenvalue associated to $\mathbf{u}_{i}$. Here $\hat{\mathbf{x}}$ represents the DFT of $\mathbf{x}$.
\end{example}
\begin{example}[Graph signal processing]\normalfont
Taking into account that the filtering in each layer of a GNN is defined by a polynomial algebra, the spectral representation of the filter is given by
\begin{multline}
\rho(a)\mathbf{x}=\sum_{i=1}^{N}\phi_{i}\left(\sum_{k=0}^{K-1}h_{k}t^{k}\right)\hat{\mathbf{x}}(i)\mathbf{u}_{i}\\
=\sum_{i=1}^{N}\sum_{k=0}^{K-1}h_{k}\phi_{i}(t)^{k}\hat{\mathbf{x}}(i)\mathbf{u}_{i}=
\sum_{i=1}^{N}\sum_{k=0}^{K-1}h_{k}\lambda_{i}^{k}\hat{\mathbf{x}}(i)\mathbf{u}_{i}
\label{eq:specopgnn}
\end{multline}
with $a=\sum_{k=0}^{K-1}h_{k}t^{k}$, and where the $\mathbf{u}_{i}$ are given by the eigenvector decomposition of $\rho(t)=\mathbf{S}$, where $\mathbf{S}$ could be the adjacency matrix or the Laplacian of the graph, while $\phi_{i}(t)=\lambda_{i}$ being $\lambda_{i}$ the eigenvalue associated to $\mathbf{u}_{i}$. The projection of $\mathbf{x}$ in each subspace $\mathcal{U}_{i}$ is given by $\hat{\mathbf{x}}(i)=\langle\mathbf{u}_{i}, \mathbf{x}\rangle$, and if $\mathbf{U}$ is the matrix of eigenvectors of $\mathbf{S}$ we have the widely known representation $\hat{\mathbf{x}}=\mathbf{U}^{T}\mathbf{x}$~\cite{gamagnns}.
\end{example}
\begin{example}[Group signal processing]\label{example:specopGroupCNNs}
\normalfont
Considering the Fourier decomposition on  general groups~\cite{steinbergrepg,terrasFG,fulton1991representation}, we obtain the spectral representation of the algebraic filters as
\begin{equation*}
\boldsymbol{a}\ast\mathbf{x}=\sum_{u,h\in G}\boldsymbol{a}(uh^{-1})\sum_{i,j,k}\frac{d_{k}}{\vert G\vert}\hat{\mathbf{x}}\left(\boldsymbol{\varphi}^{(k)}\right)_{i,j}\boldsymbol{\varphi}_{i,j}^{(k)}(h)hu,
\end{equation*}
where $\hat{\mathbf{x}}(\boldsymbol{\varphi}^{(k)})$ represents the Fourier components associated to the $k$th irreducible representation with dimension $d_{k}$ and $\boldsymbol{\varphi}^{(k)}$ is the associated unitary element. We can see that the $k$th element in this decomposition is $\sum_{i,j}\mathbf{x}(\varphi^{(k)})_{i,j}\sum_{u,h}\frac{d_{k}}{\vert G\vert}\boldsymbol{a}(uh^{-1})\boldsymbol{\varphi}_{i,j}^{(k)}(h)hu$.
\end{example}
\begin{example}[Graphon signal processing]\normalfont
According to the spectral theorem~\cite{conway1994course,aliprantis2002invitation}, it is possible to represent the action of a compact normal operator $S$ as $S\mathbf{x}=\sum_{i}\lambda_{i}\langle\boldsymbol{\varphi}_{i},\mathbf{x}\rangle\boldsymbol{\varphi}_{i}$ where $\lambda_{i}$ and $\boldsymbol{\varphi}_{i}$ are the eigenvalues and eigenvectors of $S$, respectively, and $\langle\cdot\rangle$ indicates an inner product. Then, the spectral representation of the filtering of a signal in the layer $\ell$ is given by
\begin{equation*}
\rho_{\ell}\left(p(t)\right)\mathbf{x}
=
\sum_{i}p(\lambda_{i})\langle\mathbf{x},\boldsymbol{\varphi}_{i}\rangle\boldsymbol{\varphi}_{i}=\sum_{i}\phi_{i}(p(t))\hat{\mathbf{x}}_{i}\boldsymbol{\varphi}_{i},
\end{equation*}
where $\phi_{i}(p(t))=p(\lambda_{i})$.
\end{example}

\section{Stability Theorems}\label{sec:stabilitytheorems}
It is important to remark that in Theorem~\ref{theorem:normbounds2} the value of $\delta$ is meant to have a maximum fixed value for each perturbation, therefore eqns.~(\ref{eq:theoremnormbounds2a}) and~(\ref{eq:theoremnormbounds2b}) indicate proportionality with respect to $\sup_{\mathbf{S}}\left\Vert D_{\mathcal{T}}(\mathbf{S})\right\Vert$.

Theorems~\ref{theorem:normbounds0},~\ref{theorem:normbounds2} and~\ref{theorem:HvsFrechet} show that it is possible to have stable representations with respect to perturbations when the domain of the homomorphisms is restricted to subsets of the algebra with specific properties, which are related with the type of perturbation.

We are finally ready to state the stability results for the operators associated to an AlgNN.
Theorem~\ref{theorem:stabilityAlgNN1} states how an AlgNN can be stable by the selection of an appropriate subset of filters in the algebra, for a given perturbation model. It is worth pointing out that conditions like the ones obtained in~\cite{fern2019stability} for GNNs can be considered particular instantiations of the conditions in Theorem~\ref{theorem:stabilityAlgNN1}. Notice also that the conditions imposed to guarantee the stability are not unique. Additionally, notice that Theorem~\ref{theorem:stabilityAlgNN1} can be easily extended to consider several features per layer, the reader is invited to check the details of the proof of the theorem in appendix~\ref{appendix_theoremstabilityAlgNN1} where the analysis is performed considering multiple features per layer.