# Convolutional Filters and Neural Networks with
Non Commutative Algebras

Alejandro Parada-Mayorga, Landon Butler and Alejandro Ribeiro Dept of
Electrical and Systems Engineering, University of Pennsylvania. Email:
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Supported by NSF-Simons MoDL Award #2031985

###### Abstract

In this paper we introduce and study the algebraic generalization of non
commutative convolutional neural networks. We leverage the theory of algebraic
signal processing to model convolutional non commutative architectures, and we
derive concrete stability bounds that extend those obtained in the literature
for commutative convolutional neural networks. We show that non commutative
convolutional architectures can be stable to deformations on the space of
operators. We develop the spectral representation of non commutative signal
models to show that non commutative filters process Fourier components
independently of each other. In particular we prove that although the spectral
decompositions of signals in non commutative models are associated to
eigenspaces of dimension larger than one, there exists a trade-off between
stability and selectivity, which is controlled by matrix polynomial functions
in spaces of matrices of low dimension. This tradeoff shows how when the
filters in the algebra are restricted to be stable, there is a loss in
discriminability that is compensated in the network by the pointwise
nonlinearities. The results derived in this paper have direct applications and
implications in non commutative convolutional architectures such as group
neural networks, multigraph neural networks, and quaternion neural networks,
for which we provide a set of numerical experiments showing their behavior
when perturbations are present.

###### Index Terms:

Non commutative convolutional architectures, Algebraic Neural Networks
(AlgNNs), algebraic signal processing (ASP), representation theory of
algebras, non commutative algebras, non commutative operators, non commutative
neural networks, Fréchet differentiability.

##  I Introduction

Deep learning relies on parameterizations given by the composition of layers
which are themselves compositions of linear operators with pointwise
nonlinearities. In problems that involve high dimensional inputs it becomes
necessary to exploit their structure to reduce the complexity of the learning
parametrization. This is often accomplished with the use of particular
instantiations of _convolutional_ filter banks. The most notable examples of
this approach are the use of standard (Euclidean) _convolutional_ filters for
learning with time signals and images, graph _convolutional_ filters for
learning with graphs and graph signals, and group _convolutions_ for signals
with group symmetries. The success of convolutional neural networks of
different types – Euclidean, graph, and group – is in part due to how
convolutions leverage symmetries of the domain and the signals. However, both
filter banks and convolutional neural networks are equally good at leveraging
these symmetries [1, 2, 3]. Yet, it is the latter that are more successful at
learning. This empirical observation prompts a search for properties that
explain the better performance of convolutional neural networks relative to
the corresponding convolutional filter banks.

In the case of Euclidean convolutions and graph convolutions part of the
insight into the relative performance of convolutional filter banks and neural
networks follows from their respective responses to deformations [4, 2] –
deformations in the domain of the signals in [4] and deformations on the space
of graph matrix representations in [2]. These works prove that stability to
deformations requires filters that do not distinguish different high frequency
components. It follows that convolutional neural networks can better trade-off
discriminability and stability as the pointwise nonlinearity mixes frequency
components across layers. Recent analysis conducted on the operator space has
demonstrated that the shared stability properties of Euclidean and graph
convolutions are a consequence of their shared algebraic structure [1]. It
then follows that whenever we use convolutional parameterizations in machine
learning we can expect convolutional neural networks with multiple layers to
outperform convolutional filter banks. These results build on the theory of
algebraic signal processing (ASP) which provides a common language for
describing convolutional filters of different types [5]. They apply to
Euclidean, graph, and group convolutions as well as to a large class of less
ubiquitous convolutional filters and neural networks [1, 6].

Although these results provide valuable insights, they are limited to
commutative operators. Namely, to signal processing architectures in which
filters commute. This is not true in general as non commutative signal
processing arises naturally in multiple scenarios [7, 8, 9, 10, 11, 12, 13,
14, 15]. For instance, multigraphs arise when nodes are related by several
different types of edges [7, 16]. In this case it is natural to define
multigraph filters by combining diffusions across matrix representations of
the graphs defined by each individual type (Example 3). The resulting filters
do not commute except in the rare eventuality that the different matrix
representations themselves commute. Filters that do not commute also arise
with signal models on groups like ${\mathbb{E}}(n)$, or $SO(n)$ [8, 9, 10,
17], or any nontrivial Lie group [11, 12, 17, 18]. Non commutative filters
appear here because filters are associated with group symmetries and it is
common for groups to have symmetries that do not commute – such as rotations
in $SO(3)$.

Existing stability results do not apply to any of these non commutative
settings but empirical results show the usual advantage of layered neural
networks with respect to filter banks [16, 17, 11, 12, 15]. Such an empirical
observation poses the question of whether similar stability results hold for
non commutative convolutional architectures. That is, whether convolutional
versions of architectures such as multigraph neural networks [16, 19], neural
networks on groups [17, 15], quaternion (graph) neural networks [14, 20, 21,
22, 23, 24, 25], hypercomplex algebras neural networks [26], hyperbolic neural
networks [27, 28], octonion neural networks [29], and Clifford algebras neural
networks [30, 31, 32] can be stable to deformations and, more to the point,
whether stability requires spectral restrictions analogous to those of
Euclidean, graph, and commutative group neural networks. The goal of this
paper is to answer this question in the affirmative. To do so our first two
contributions are the following:

  * (C1)

We define convolutional filters (Section II) and neural networks (Section IV)
with non commutative algebras.

  * (C2)

We develop the spectral representation of non commutative signal models
(Section III).

Non commutative filters are just a particular case of algebraic filters [5]
and the construction of algebraic neural networks with non commutative filters
is a close to verbatim extension of algebraic neural networks [1, 6]. The
development of spectral representations, however, has substantial differences.
Commutative signal models define Fourier decompositions as projections on
single dimensional subspaces and, consequently, frequencies are defined as
scalars. In non commutative signal models Fourier decompositions are
projections in multidimensional subspaces and frequencies are, consequently,
matrices of corresponding dimensions (Definition 5). Despite this significant
difference we can still prove the equivalent of an spectral representation
theorem of non commutative algebraic filters:

  * (C3)

Non commutative algebraic filters process Fourier components independently of
each other (Theorem 1).

Contribution (C3) is the equivalent of the claim that commutative algebraic
filters process frequency components independently of each other. As in the
case of commutative filters, Contribution (C3) implies that non commutative
algebraic filters are completely characterized by their frequency responses.
The difference is that in non commutative filters the frequency response is a
_matrix_ polynomial (Definition 6). This is in contrast to the scalar
polynomials that define frequency representations in commutative signal models
and a consequence of the fact that frequencies are matrices, not scalars.

It is this difference between frequency representations of commutative and non
commutative filters – matrix versus scalar polynomials – that prevents the
results derived in [1] to be applied to non commutative convolutional
architectures. The main technical contribution of this paper is to generalize
the analysis of [1] to show that analogous results hold for non commutative
algebraic signal models. In particular, we prove that:

  * (C4)

Non commutative algebraic convolutional filters _can_ be stable to additive
and multiplicative perturbations of the algebraic signal model (Section VI).

  * (C5)

Non commutative algebraic neural networks inherit the stability properties of
algebraic convolutional filters (Section VI-A).

The proof of Contribution (C4) and (C5) requires that filters have two
spectral properties (Definition 9). The first property is a generalized
Lipschitz condition in which changes in the filter’s frequency response are
upper bounded by a linear function of the frequency’s norm. The second
condition is that the Fréchet derivative of the filter’s response acting on a
frequency matrix has bounded norm. These two conditions limit the
discriminability of non commutative filters. They respectively imply that: (i)
The variability of a stable filter must be bounded by the norm of the
difference between frequency matrices. (ii) The variability of stable filters
must decrease with the norm of the frequency matrix. Stability of commutative
filters requires analogous conditions on scalar frequencies. Since imposing
conditions on frequency matrices is more stringent, our stability results
suggest that non commutative filters require more stringent filter
restrictions to attain the same level of stability.

To illustrate our theoretical results we provide a set of numerical
experiments for non commutative convolutional multigraph and quaternion neural
networks (Section VII).

##  II Non Commutative Filters

${\mathcal{M}}$${\mathbf{x}}$$\rho(th){\mathbf{x}}$$\rho(ht){\mathbf{x}}$$\text{End}({\mathcal{M}})$$\rho(th)$$\rho(ht)$${\mathcal{A}}$$th$$ht$$\rho\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ $$\rho\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ $ Figure 1: Non commutative
algebraic signal model. The algebraic filters $th$ and $ht$ are realized
physically in $\text{End}(\mathcal{M})$ to process the signals ${\mathbf{x}}$
which are modeled as elements of ${\mathcal{M}}$.

In this section we introduce the generalization of non commutative signal
models on arbitrary domains under the lens of algebraic signal processing
(ASP). An algebraic signal model (ASM) is defined as the triplet

$({\mathcal{A}},{\mathcal{M}},\rho),$ (1)

where ${\mathcal{A}}$ is a unital associative algebra, ${\mathcal{M}}$ is a
vector space, and $\rho:{\mathcal{A}}\rightarrow\text{End}({\mathcal{M}})$ is
a homomorphism between the algebra ${\mathcal{A}}$ and the set of
endomorphisms in the vector space ${\mathcal{M}}$ [5]. We recall that an
algebra is a vector space with a closed operation of product, the endomorphims
of ${\mathcal{M}}$ are the set of linear maps from ${\mathcal{M}}$ onto itself
and a homomorphism is a linear map between algebras that preserves the product
operation. We refer the reader to Appendix A for a review of these concepts
and the discussion of some examples as well.

For our discussion we consider that ${\mathcal{A}}$ and ${\mathcal{M}}$ as
vector spaces are defined on an algebraically closed field ${\mathbb{F}}$
which for the sake of simplicity will be considered
${\mathbb{F}}={\mathbb{C}}$. However, we remark that the material discussed is
valid for any algebraically closed Field. In the triplet
$({\mathcal{A}},{\mathcal{M}},\rho)$ the pair $({\mathcal{M}},\rho)$
constitutes a representation of ${\mathcal{A}}$ in the context of
representation theory of algebras.

The processing of information in the context of (1) is known as algebraic
signal processing (ASP). The signals are modeled as elements of
${\mathcal{M}}$, the filters are defined as elements of ${\mathcal{A}}$ and
their realization is given by $\rho$. The rules or laws governing the
operations are given by the structure of ${\mathcal{A}}$ and $\rho$ translates
those operations into actions on the elements of ${\mathcal{M}}$. Then, the
filtering of a signal ${\mathbf{x}}\in{\mathcal{M}}$ by an algebraic filter
$a\in{\mathcal{A}}$ is given by

${\mathbf{y}}=\rho(a){\mathbf{x}}.$ (2)

In (2) we have a generalized representation of the convolution operation
between a filter and a signal. Particular instantiations of (2) lead to the
traditional signal processing models of time signals and images and more
sophisticated models such as graph signal processing (GSP), graphon signal
processing (WSP) among others [1, 6, 33].

Since it is ${\mathcal{A}}$ the algebraic object defining the structural rules
of information processing, the non commutativity in an algebraic model is
associated to ${\mathcal{A}}$. This is, any non commutative convolutional
signal model is given by $({\mathcal{A}},{\mathcal{M}},\rho)$ where
${\mathcal{A}}$ is a non commutative algebra. For the reader unfamiliar with
algebraic concepts, we remark that one representative example of this type of
algebra is $M_{n\times n}({\mathbb{C}})$, which is the set of matrices of
dimension $n\times n$ with entries in ${\mathbb{C}}$. It is an algebra, since
it is a vector space over ${\mathbb{C}}$ and the product of two matrices of
size $n\times n$ is again a matrix of the same size.

Now we recall the notion of generators in algebraic signal models.

######  Definition 1 (Generators).

For an associative algebra with unity ${\mathcal{A}}$ we say the set
${\mathcal{G}}\subseteq{\mathcal{A}}$ generates ${\mathcal{A}}$ if all
$a\in{\mathcal{A}}$ can be represented as polynomial functions of the elements
of ${\mathcal{G}}$. We say elements $g\in{\mathcal{G}}$ are generators of
${\mathcal{A}}$ and we denote as $a=p({\mathcal{G}})$ the polynomial that
generates $a$.

Representing the elements of ${\mathcal{A}}$ as polynomial functions of a set
of generators highlights the fact that such generators characterize the
algebra, and can be conceived as a measure of the degrees of freedom
associated to the algebra. Consequently, the realization of the generators by
means of the homomorphism $\rho$, called shift operators, play a central role
when analyzing algebraic signal models. We introduce its formal definition
next.

######  Definition 2 (Shift Operators).

Let $({\mathcal{A}},{\mathcal{M}},\rho)$ be an algebraic signal model. Then,
if ${\mathcal{G}}\subset{\mathcal{A}}$ is a generator set of ${\mathcal{A}}$,
the operators ${\mathbf{S}}=\rho(g)$ with $g\in{\mathcal{G}}$ are called shift
operators. The set of all shift operators is denoted by ${\mathcal{S}}$.

Then, we can express any element $\rho(a)\in\text{End}({\mathcal{M}})$ as a
polynomial function of the shift operators. In particular, if
$a=p({\mathcal{G}})$ we have that

$\rho(a)=p_{\mathcal{M}}\big{(}\rho({\mathcal{G}})\big{)}=p_{\mathcal{M}}\big{(}{\mathcal{S}}\big{)}=p\big{(}{\mathcal{S}}\big{)},$
(3)

where $p_{\mathcal{M}}$ indicates the polynomial of form $p$ but whose
independent variables are the shift operators. For the sake of simplicity we
drop the subindex in $p$ in our subsequent discussion since it is clear where
$p$ is defined from the independent variable.

To facilitate the understanding of these basic concepts we present an example.

$1$$\mathbf{x}(1)$$2$$\mathbf{x}(2)$$3$$\mathbf{x}(3)$$4$$\mathbf{x}(4)$$5$$\mathbf{x}(5)$$6$$\mathbf{x}(6)$$\mathbf{S}_{1}\rightarrow\mathcal{E}_{1}$
$\mathbf{S}_{2}\rightarrow\mathcal{E}_{2}$ 1726354453627171 Figure 2: Example
of a multigraph signal processing model. The signal
${\mathbf{x}}\in{\mathbb{R}}^{6}$ defined on the multigraph
$G=({\mathcal{V}},\\{{\mathcal{E}}_{1},{\mathcal{E}}_{2}\\})$ is associated to
the nodes ${\mathcal{V}}$, while the shift operator ${\mathbf{S}}_{i}$ is
associated to the set of edges $\mathcal{E}_{i}$.

######  Example 1 (Multigraph signal processing).

A multigraph consists of a common set of vertices and multiple separate sets
of edges. Consider the specific case
$G=\left({\mathcal{V}},\\{{\mathcal{E}}_{1},{\mathcal{E}}_{2}\\}\right)$
consisting of a set of $N$ vertices ${\mathcal{V}}$ and two separate edge sets
${\mathcal{E}}_{1}$ and ${\mathcal{E}}_{2}$ – see Fig. 2. Associated with each
edge set we consider matrix representations ${\mathbf{S}}_{1}$ and
${\mathbf{S}}_{2}$. This is an algebraic signal processing model in which the
vector space of signals ${\mathcal{M}}={\mathbb{C}}^{N}$ contains vectors in
${\mathbb{C}}^{N}$ with entries associated with each node of the multigraph
and the algebra ${\mathcal{A}}={\mathbb{C}}[t_{1},t_{2}]$ is the set of non
commutative polynomials of two variables. The algebra
${\mathbb{C}}[t_{1},t_{2}]$ is generated by the monomials $t_{1}$ and $t_{2}$.
The homomorphism $\rho$ is defined by mapping the generators $t_{1}$ and
$t_{2}$ to $\rho(t_{1})={\mathbf{S}}_{1}$ and $\rho(t_{2})={\mathbf{S}}_{2}$,
where the shift operators ${\mathbf{S}}_{i}$ are the matrix representations of
the corresponding set of edges ${\mathcal{E}}_{i}$. Then, if
$p(t_{1},t_{2})=t_{1}^{2}+t_{1}t_{2}+2t_{2}t_{1}+t_{2}^{2}+1$ the filtering in
(2) takes the form

$\rho\left(t_{1}^{2}+t_{1}t_{2}+2t_{2}t_{1}+t_{2}^{2}+1\right){\mathbf{x}}\\\
=\left({\mathbf{S}}_{1}^{2}+{\mathbf{S}}_{1}{\mathbf{S}}_{2}+2{\mathbf{S}}_{2}{\mathbf{S}}_{1}+{\mathbf{S}}_{2}^{2}+1\right){\mathbf{x}}.$
(4)

We refer the reader to Appendix F in the supplementary material for a more
sophisticated example on multigraph signal processing in the light of ASP. In
Appendix B and Appendix C of the paper it is shown how classical group signal
processing and convolutional quaternion signal processing can be seen as
particular cases of algebraic signal models.

######  Remark 1.

It is important to highlight that the definition of an algebraic signal model
is independent of any norm that could be associated to ${\mathcal{M}}$, if
any. As we will discuss in subsequent sections (see Section VI) we will endow
${\mathcal{M}}$ with a norm in order to analyze the size of perturbations, and
doing this will not affect the structure of the triplet
$({\mathcal{A}},{\mathcal{M}},\rho)$. We also remark that the algebra
${\mathcal{A}}$ is not associated to a norm, but it happens that
${\mathcal{A}}$ can be isomorphic to algebras that are naturally endowed with
a norm. This is something we exploit in Section III to characterize subsets of
the algebra and spectral representations.

##  III Spectral Representations of Non Commutative Filters

In this section we discuss spectral representations of filters and signals. We
will show that spectral representations in non commutative signal models are
determined by matrix polynomials and that frequencies themselves are described
by matrices instead of scalars. This introduces technical challenges, and yet
we will also show that despite the significant differences with respect to the
commutative scenario, spectral filtering is an analog of spectral filtering in
commutative signal models. Additionally, this analog behavior extends to the
norms of the operators, in particular the norm of non commutative filters is
determined by the norm of the spectral responses which can be characterized in
analogous ways to those in commutative signal models when considering induced
norm on direct sums of spaces.

${\mathcal{M}}$${\mathcal{U}}_{1}$$\bigoplus$${\mathcal{U}}_{2}$$\bigoplus$${\mathcal{U}}_{3}$$\phi_{1}(a)$$\text{End}({\mathcal{U}}_{1})$$\phi_{2}(a)$$\text{End}({\mathcal{U}}_{2})$$\phi_{3}(a)$$\text{End}({\mathcal{U}}_{3})$$\text{End}({\mathcal{M}})$$\rho(a)$$\cong$${\mathcal{A}}$$a$$\phi_{3}$$\phi_{2}$$\phi_{1}$$\rho\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ $ Figure 3: Spectral decomposition
in an ASM: The maps $\phi_{i}$ indicate the frequencies while each
${\mathcal{U}}_{i}$ is an eigenspace. Notice that if $a\in{\mathcal{A}}$ is a
polynomial in terms of the generators of ${\mathcal{A}}$, then $\rho(a)$ and
$\phi_{i}(a)$ are polynomial operators in terms of the shift operators in
$\text{End}(\mathcal{M})$ and $\text{End}({\mathcal{U}}_{i})$, respectively.
When $\dim({\mathcal{U}}_{i})=1$, $\phi_{i}(a)$ are polynomials in terms of
the $i$th eigenvalue of the operator $\rho(a)$.

The notion of spectral decomposition descends from the concepts of irreducible
and completely reducible subrepresentations of an algebra [5, 34, 35, 36, 37].
Widely known notions of spectral decompositions used in graph signal
processing (GSP), discrete signal processing (DSP), and discrete time signal
processing (DTSP) among others, are obtained as particular cases of a
decomposition of representations of an algebra as a sum of irreducible
subrepresentations. Then, in order to present a natural generalization of
frequency decompositions in non commutative ASM, we discuss the concepts of
subrepresentation, irreducibility and decomposability. We restrict our
attention to algebras that have a finite number of generators.

######  Definition 3.

Let $(\mathcal{M},\rho)$ be a representation of $\mathcal{A}$. Then, a
representation $(\mathcal{U},\rho)$ of $\mathcal{A}$ is a subrepresentation of
$(\mathcal{M},\rho)$ if $\mathcal{U}\subseteq\mathcal{M}$ and $\mathcal{U}$ is
invariant under all operators $\rho(a)\leavevmode\nobreak\ \text{for
all}\leavevmode\nobreak\ a\in\mathcal{A}$, i.e. $\rho(a)u\in\mathcal{U}$ for
all $u\in\mathcal{U}$ and $a\in\mathcal{A}$.

The notion of subrepresentation is tied to the property of invariance since
subrepresentations are invariant under the action of the instantiations of
elements of ${\mathcal{A}}$ in $\text{End}(\mathcal{M})$. Notice that when
considering representations of an algebra with a single generator, the spaces
generated by subsets of eigenvectors of $\rho(a)$ determine the
subrepresentations of $({\mathcal{M}},\rho)$.

Now, we introduce formally irreducible subrepresentations which provide a
minimal structural unit of invariance.

######  Definition 4.

A representation $(\mathcal{M},\rho)$ (with ${\mathcal{M}}\neq 0$) is
irreducible or simple if the only subrepresentations of $(\mathcal{M},\rho)$
are $(0,\rho)$ and $(\mathcal{M},\rho)$.

The irreducibility property of a subrepresentation implies that there is not a
subspace that has its own invariance under the action of the elements of the
algebra. Additionally, notice that one dimensional subrepresentations are
always irreducible. It is important to point out that the calculation of
irreducible subrepresentations entails a significant computational cost – see
Appendix K in the supplementary material.

As we will show in the next subsection, irreducibility can be used to write a
general representation as a decomposition in terms of irreducible
subrepresentations. To show this we introduce the notion of direct sum of
representations. Given two representations $({\mathcal{M}}_{1},\rho_{1})$ and
$({\mathcal{M}}_{1},\rho_{2})$ of an algebra ${\mathcal{A}}$ we can obtain a
new representation, called the direct sum representation, as
$({\mathcal{M}}_{1}\oplus{\mathcal{M}}_{1},\rho)$ where $\rho(a)(u_{1}\oplus
u_{2})=\rho_{1}(a)u_{1}\oplus\rho_{2}(a)u_{2}$.

###  III-A Fourier Decompositions and Spectral Representation of Filters

Using the ideas and concepts discussed above we introduce the notion of
Fourier decomposition in algebraic signal processing.

######  Definition 5 (Fourier Decomposition).

For an algebraic signal model $(\mathcal{A},\mathcal{M},\rho)$ we say that
there is a spectral or Fourier decomposition if

$(\mathcal{M},\rho)\cong\bigoplus_{(\mathcal{U}_{i},\phi_{i})\in\text{Irr}\\{\mathcal{A}\\}}(\mathcal{U}_{i},\phi_{i}),$
(5)

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

$\Delta:\mathcal{M}\to\bigoplus_{(\mathcal{U}_{i},\phi_{i})\in\text{Irr}\\{\mathcal{A}\\}}\mathcal{U}_{i}$
(6) $\mathbf{x}\mapsto\hat{\mathbf{x}},$

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)$.

In Definition 5 it is assumed that each individual subrepresentation
$({\mathcal{U}}_{i},\phi_{i})$ cannot be expressed as a direct sum of
irreducible subrepresentations isomorphic to $({\mathcal{U}}_{i},\phi_{i})$.
This assumption is analogous to the assumption of not having repeated
eigenvalues when considering spectral decompositions where the irreducible
subrepresentations have dimension equal to 1. This choice is done for the sake
of simplicity and we refer the reader to Appendix K of the supplementary
material for a more sophisticated formulation of the Fourier decompositions.
We remark that the decomposition in (5) and 6 is unique up to isomorphism –
see Appendix K in the supplementary material.

It is worth pointing out that the homomorphisms $\phi_{i}$ associated to each
non isomorphic irreducible representation define the frequency associated to
the vector space ${\mathcal{U}}_{i}$ [5] – see Fig. 3. In commutative
scenarios, like for instance in GSP, we have $\dim({\mathcal{U}}_{i})=1$. In
this case the term $\phi_{i}(a)$ is a scalar value which corresponds to the
eigenvalues of $\rho(a)$ while ${\mathcal{U}}_{i}$ is the space spanned by the
$i$-th eigenvector.111We refer the reader to Appendix G in the supplementary
material where we show with more details how the Fourier transform in GSP can
be obtained as a particular case of Definition 5.. If
$\dim({\mathcal{U}}_{i})>1$ then $\phi_{i}(a)$ is a matrix. This last scenario
is typical when non commutative algebras are involved.

$\text{End}(\mathcal{M})$$\tilde{\mathcal{A}}$$\tilde{\mathcal{C}}$$\mathcal{A}$$\mathcal{C}$$\rho$$\iota$
Figure 4: Schematic representation of the spectral representation
$\tilde{{\mathcal{C}}}$ of a subset of filters
${\mathcal{C}}\subset{\mathcal{A}}$ algebraic signal model. The properties of
$\rho({\mathcal{C}})$ are determined by ${\mathcal{C}}$ which is characterized
by the elements in $\tilde{{\mathcal{C}}}$.

We synthesize the role of Fourier decompositions and the filtering operation
in the following theorem.

######  Theorem 1.

(Filtering Spectral Theorem) Let $({\mathcal{A}},{\mathcal{M}},\rho)$ be an
algebraic signal model where ${\mathcal{A}}$ has generators
$\\{g_{k}\\}_{k=1}^{m}$, and $\rho(a)$ is the realization of
$a\in{\mathcal{A}}$ in $\text{End}(\mathcal{M})$ by means of $\rho$. If
$(\mathcal{M},\rho)\cong\bigoplus_{(\mathcal{U}_{i},\phi_{i})\in\text{Irr}\\{\mathcal{A}\\}}(\mathcal{U}_{i},\phi_{i})$
then

$\rho\left(p(g_{1},\ldots,g_{m})\right)=p\left({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m}\right)=\\\
\sum_{i}p\left(\boldsymbol{\Lambda}_{1,i},\ldots,\boldsymbol{\Lambda}_{m,i}\right)\mathrsfso{P}_{i},$
(7)

where $\mathrsfso{P}_{i}$ is the projection operator on ${\mathcal{U}}_{i}$
and $\boldsymbol{\Lambda}_{k,i}=\phi_{i}(g_{k})$. Additionally, if
${\mathbf{y}}=\rho\left(p(g_{1},\ldots,g_{m})\right){\mathbf{x}}$, then

$\hat{{\mathbf{y}}}_{i}=p\left(\boldsymbol{\Lambda}_{1,i},\ldots,\boldsymbol{\Lambda}_{m,i}\right)\hat{{\mathbf{x}}}_{i},$
(8)

where the sub index $i$ indicates a projection on ${\mathcal{U}}_{i}$.

###### Proof.

See Appendix G in the supplementary material. ∎

Theorem 1 exhibits similarties in form with the classical spectral theorem for
commutative signal models but also substantial differences. Indeed, Theorem 1
is a generalization of the classical spectral theorem where decompositions of
operators are expressed in terms of projections on spaces spanned by more than
just one vector. As indicated in Appendix G, Theorem 1 follows from the fact
that the restriction of the homomorphism $\rho$ to an invariant subspace
${\mathcal{U}}_{i}$ given by $\phi_{i}$ is again a homomorphism. Therefore, if
$p(g_{1},\ldots,g_{m})$ is a polynomial so it is $\rho(p(g_{1},\ldots,g_{m}))$
and $\phi_{i}(p(g_{1},\ldots,g_{m}))$ – and they have the same coefficients.
Additionally, the dimensions of ${\mathcal{M}}$ and ${\mathcal{U}}_{i}$ will
determine the dimensions of the independent variables in
$\rho(p(g_{1},\ldots,g_{m}))$ and $\phi_{i}(p(g_{1},\ldots,g_{m}))$.

The algebraic filter $a\in{\mathcal{A}}$ determines the form and properties of
the spectral response indicated by the homomorphisms $\phi_{i}$. Therefore,
classes of filters in ${\mathcal{A}}$ lead to specific classes of filters
defined in the spectral domain. We use this fact and an auxiliary matrix
algebra, $\tilde{{\mathcal{A}}}$, to provide a concrete characterization of
subsets of filters in ${\mathcal{A}}$. $\tilde{{\mathcal{A}}}$ is endowed with
a norm and is isomorphic to ${\mathcal{A}}$ – see Fig. 4.

In $\tilde{{\mathcal{A}}}$ the generators can be considered as variables
taking specific values on a vector space. In particular, given any
$({\mathcal{A}},{\mathcal{M}},\rho)$ with ${\mathcal{A}}$ generated by
$\\{g_{i}\\}_{i=1}^{m}$ we can associate
$\tilde{{\mathcal{A}}}=\left\\{\left.p:M_{r,r}^{m}(\mathbb{C})\mapsto
M_{r,r}(\mathbb{C})\right|p:\text{polynomial}\right\\},$ where
$M_{r,r}(\mathbb{C})$ is the space of matrices of size $r\times r$ whose
entries belong to $\mathbb{C}$, while $M_{r,r}^{m}(\mathbb{C})$ is the
$m$-times cartesian product of $M_{r,r}(\mathbb{C})$. Additionally, we endow
$\tilde{{\mathcal{A}}}$ with a norm. It is possible to see that there is a
natural isomorphism $\iota$ between ${\mathcal{A}}$ and
$\tilde{{\mathcal{A}}}$ given by $\iota(g_{i})=\boldsymbol{\Lambda}_{i}\in
M_{r,r}(\mathbb{C})$. We attribute the properties of a set of filters realized
in the set $\tilde{{\mathcal{C}}}\subset\tilde{{\mathcal{A}}}$ to the set
${\mathcal{C}}=\left\\{a\in{\mathcal{A}}|\quad\iota(a)\in\tilde{{\mathcal{C}}}\right\\}.$
The value of $r$ is selected as $r=\max\left\\{d_{i}\right\\},$ where $d_{i}$
are the dimensions of the irreducible subrepresentations of
$({\mathcal{M}},\rho)$. It is worth pointing out that the value of the $d_{i}$
depends on $\rho$, and in general $d_{i}\ll\dim({\mathcal{M}})$ [38, 35, 1,
5]. Taking into account this, we formally introduce the spectral
representation of a filter.

######  Definition 6.

Let $({\mathcal{A}},{\mathcal{M}},\rho)$ be an ASM where ${\mathcal{A}}$ is an
algebra with generators $\\{g_{i}\\}_{i=1}^{m}$. Let $\tilde{{\mathcal{A}}}$
the matrix algebra isomorphic to ${\mathcal{A}}$ given by

$\tilde{{\mathcal{A}}}=\left\\{\left.p(\boldsymbol{\Lambda}_{1},\ldots,\boldsymbol{\Lambda}_{m}):M_{r,r}^{m}(\mathbb{C})\mapsto
M_{r,r}(\mathbb{C})\right|p:\text{polynomial}\right\\},$ (9)

where $r=\max\left\\{d_{i}\right\\},$ and $d_{i}$ are the dimensions of the
irreducible subrepresentations of $({\mathcal{M}},\rho)$. Then, we say that
$p(\boldsymbol{\Lambda}_{1},\ldots,\boldsymbol{\Lambda}_{m})\in\tilde{{\mathcal{A}}}$
is the spectral representation of $p(g_{1},\ldots,g_{m})\in{\mathcal{A}}$.

######  Remark 2.

It is important to remark that a substantial part of the technical challenges
that have made the analysis of non commutative signal models and architectures
elusive, is precisely the difficulty associated to the use of the Fourier
representations. As we just discussed, some frequencies are not described by
scalars but instead by matrices. Part of our contribution in this paper is
precisely the fact that we have overcome these technical challenges and as the
reader can corroborate in Appendix D the formal proof of our stability results
make use of tools different from those used in [1] for commutative signal
models where frequencies are associated to scalars.

As a final comment it is worth highlighting that the notion of aliasing on the
spectral domain for non commutative signal models is analog to that one of
commutative models. The difference is that the representation of the aliased
information is expressed in terms the basis of the ${\mathcal{U}}_{i}$ spaces
where aliasing takes place.

##  IV Algebraic Neural Networks with Non Commutative Algebras

$({\mathcal{A}}_{1},{\mathcal{M}}_{1},\rho_{1})$$({\mathcal{A}}_{2},{\mathcal{M}}_{2},\rho_{2})$$({\mathcal{A}}_{3},{\mathcal{M}}_{3},\rho_{3})$$\mathbf{x}$$\displaystyle{\mathbf{y}_{1}=\rho_{1}({\color[rgb]{0.265625,0.00390625,0.328125}\definecolor[named]{pgfstrokecolor}{rgb}{0.265625,0.00390625,0.328125}a_{1}})\,\mathbf{x}}$$\displaystyle{\mathbf{z}_{1}={\eta_{1}}\Big{[}\,\mathbf{y}_{1}\,\Big{]}}$$\displaystyle{\mathbf{x}_{1}={P_{1}}\Big{[}\,\mathbf{z}_{1}\,\Big{]}}$$\mathbf{y}_{1}$$\mathbf{z}_{1}$$\displaystyle{\mathbf{y}_{2}=\rho_{2}({\color[rgb]{0.265625,0.00390625,0.328125}\definecolor[named]{pgfstrokecolor}{rgb}{0.265625,0.00390625,0.328125}a_{2}})\,\mathbf{x}_{1}}$$\displaystyle{\mathbf{z}_{2}={\eta_{2}}\Big{[}\,\mathbf{y}_{2}\,\Big{]}}$$\displaystyle{\mathbf{x}_{2}={P_{2}}\Big{[}\,\mathbf{z}_{2}\,\Big{]}}$$\mathbf{y}_{2}$$\mathbf{z}_{2}$$\displaystyle{\mathbf{y}_{3}=\rho_{3}({\color[rgb]{0.265625,0.00390625,0.328125}\definecolor[named]{pgfstrokecolor}{rgb}{0.265625,0.00390625,0.328125}a_{3}})\,\mathbf{x}_{2}}$$\displaystyle{\mathbf{z}_{3}={\eta_{3}}\Big{[}\,\mathbf{y}_{3}\,\Big{]}}$$\displaystyle{\mathbf{x}_{3}={P_{3}}\Big{[}\,\mathbf{z}_{3}\,\Big{]}}$$\mathbf{y}_{3}$$\mathbf{z}_{3}$$\mathbf{x}_{1}$$\mathbf{x}_{1}$$\mathbf{x}_{2}$$\mathbf{x}_{2}$$\mathbf{x}_{3}$
Figure 5: Algebraic Neural Network
$\\{({\mathcal{A}}_{\ell},{\mathcal{M}}_{\ell},\rho_{\ell};\sigma_{\ell})\\}_{\ell=1}^{3}$
with three layers indicating how the input signal $\mathbf{x}$ is processed
and mapped into $\mathbf{x}_{3}$. In each layer the information is transformed
by a convolutional filter $\rho_{\ell}(a_{\ell})$, followed by a pointwise
nonlinearity $\eta_{\ell}$ and a pooling operator $P_{\ell}$.

Algebraic neural networks (AlgNNs) are stacked layered structures (see Fig. 5)
where the processing of information in each layer is carried out by means of
filters of an algebraic signal model and pointwise nonlinear operators. For a
training set ${\mathcal{T}}=\\{{\mathbf{x}},{\mathbf{y}}\\}$ with inputs
${\mathbf{x}}$ and outputs ${\mathbf{y}}$, it is possible to learn the
algebraic filters $a_{\ell}\in{\mathcal{A}}_{\ell}$ in each layer of the AlgNN
to produce a mapping representation. This allows an estimation of the output
to an unseen input $\tilde{{\mathbf{x}}}\in{\mathcal{T}}$. The data from the
training set ${\mathcal{T}}$ is used to find subsets of filters
${\mathcal{P}}_{\ell}\subset{\mathcal{A}}_{\ell}$ that minimize a cost
function of the form
$\sum_{({\mathbf{x}},{\mathbf{y}})\in{\mathcal{T}}}f_{{\mathcal{P}}_{\ell}}({\mathbf{x}},{\mathbf{y}}),$
where $f_{{\mathcal{P}}_{\ell}}({\mathbf{x}},{\mathbf{y}})$ is a fitting
metric that penalizes the difference between ${\mathbf{y}}$ and the output of
the AlgNN produced when the input is ${\mathbf{x}}$. In the $\ell$-th layer of
the AlgNN, an incoming signal ${\mathbf{x}}_{\ell-1}$ from the layer $\ell-1$
is filtered by means of the convolution
$\rho_{\ell}(a_{\ell}){\mathbf{x}}_{\ell-1}$. Then, a pointwise nonlinear
operator $\eta_{\ell}:{\mathcal{M}}_{\ell}\rightarrow{\mathcal{M}}_{\ell}$ is
applied and finally a pooling operator
$P_{\ell}:{\mathcal{M}}_{\ell}\rightarrow{\mathcal{M}}_{\ell+1}$ matches
information between ${\mathcal{M}}_{\ell}$ and ${\mathcal{M}}_{\ell+1}$. The
output signal of the layer $\ell$ can be written as

${\mathbf{x}}_{\ell}=\sigma_{\ell}\left(\rho_{\ell}(a_{\ell}){\mathbf{x}}_{\ell-1}\right)=\Phi(\mathbf{x}_{\ell-1},\mathcal{P}_{\ell},\mathcal{S}_{\ell}).$
(10)

We use the symbol
$\Phi(\mathbf{x}_{\ell-1},\mathcal{P}_{\ell},\mathcal{S}_{\ell})$ to make
emphasis in the fact that the filters used in each layer belong to specific
subsets of the algebra and that a specific family of shift operators is being
used. If several features per layer are used, we use the notation

${\mathbf{x}}_{\ell}^{f}=\sigma_{\ell}\left(\sum_{g=1}^{F_{\ell}}\rho_{\ell}\left(a_{\ell}^{gf}\right){\mathbf{x}}_{\ell-1}^{g}\right),$
(11)

where the super index $f$ indicates the $f$th feature and $\ell$ the layer
were information is being processed. To denote an AlgNN with $L$ layers we use
$\left\\{\left({\mathcal{A}}_{\ell},{\mathcal{M}}_{\ell},\rho_{\ell};\sigma_{\ell}\right)\right\\}_{1}^{L}$,
where $\sigma_{\ell}=P_{\ell}\eta_{\ell}$. As pointed out in [1], traditional
neural networks (CNNs), graph neural networks (GNNs), and graphon neural
networks (WNNs) among others, can be obtained as particular instantiations of
a general AlgNN with a commutative algebra. For our discussion we consider
that $\sigma_{\ell}=P_{\ell}\eta_{\ell}$ is $C_{\ell}$-Lipschitz and
$\sigma_{\ell}(0)=0$. The role of $\eta_{\ell}$ is crucial for the performance
of any AlgNN, since it is the nonlinearity of $\eta_{\ell}$ what allows the
AlgNN to redistribute spectral information and compensate restrictions imposed
on the filters. There are many possible choices for $\eta_{\ell}$ that could
improve the generalization capacity of the AlgNN. Typically a low
computational cost $\eta_{\ell}$ is selected. The operator
$P_{\ell}:{\mathcal{M}}_{\ell}\rightarrow{\mathcal{M}}_{\ell+1}$ performs an
operation of dimensionality reduction. In particular instantiations of AlgNNs
like GNNs, such operator can be associated to graph coarsening techniques or
optimal sampling strategies for bandlimited signals [39]. In what follows we
present some examples of concrete non commutative convolutional architectures
instantiated as particular cases of non commutative algebraic neural networks
(AlgNNs).

######  Remark 3.

Note that the action of the pointwise nonlinearity, $\eta_{\ell}$, and the
pooling operator, $P_{\ell}$, is defined in terms of a basis of
${\mathcal{M}}_{\ell}$. This is, given a basis
$\displaystyle\left\\{b_{i}^{(\ell)}\right\\}_{i}$ of ${\mathcal{M}}_{\ell}$
we have
$\eta_{\ell}\left(\sum_{i}\alpha_{i}b_{i}^{(\ell)}\right)=\sum_{i}\eta_{\ell}(\alpha_{i})b_{i}^{(\ell)}$,
and
$P_{\ell}:\text{span}(\\{b_{i}^{(\ell)}\\}_{i})\to\text{span}(\\{b_{i}^{(\ell+1)}\\}_{i})$.
In our analysis, the specifics of such basis are not relevant as long as the
operator $\sigma_{\ell}=P_{\ell}\eta_{\ell}$ is Lipschitz and
$\sigma_{\ell}(0)=0$. Notice also that since $P_{\ell}$ performs
dimensionality reduction, there are multiple choices for linear and nonlinear
versions of $P_{\ell}$ depending on the properties of ${\mathcal{M}}_{\ell}$
[40, 41, 42] and we consider those choices for which
$\sigma_{\ell}=P_{\ell}\eta_{\ell}$ is Lipschitz and $\sigma_{\ell}(0)=0$. We
point out however, that the optimal selection of such basis and the operators
$\eta_{\ell}$ and $P_{\ell}$ opens up an interesting future research
direction.

######  Example 2 (Multigraph neural networks).

Let us consider a multigraph
$G=({\mathcal{V}},\\{{\mathcal{E}}_{r}\\}_{r=1}^{m})$ with set of nodes
${\mathcal{V}}$, $|{\mathcal{V}}|=N$ and multiset of edges
$\\{{\mathcal{E}}_{r}\\}_{r=1}^{m}$. Let ${\mathbf{S}}_{i}$ be a matrix
representation of the multigraph on the set of edges ${\mathcal{E}}_{i}$,
which could be the adjacency matrix or a Laplacian matrix of the graph
$({\mathcal{V}},{\mathcal{E}}_{i})$. Then, the $\ell$-th layer of a multigraph
neural network is composed by convolutional operators followed by a pointwise
nonlinearity and a pooling operator. The convolutional filters are given by
the multivariable polynomials (see Example 3):
$p({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m}).$ Then, given the input
${\mathbf{x}}_{\ell}$ to the $\ell$-th layer of the multigraph neural network,
we leverage symmetries on the multigraph to obtain
${\mathbf{y}}_{\ell}=p({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m}){\mathbf{x}}_{\ell}$.
After this, we apply a pointwise nonlinear operator
$\eta_{\ell}:{\mathbb{R}}^{N}\rightarrow{\mathbb{R}}^{N}$ to obtain
${\mathbf{z}}_{\ell}=\eta_{\ell}\left({\mathbf{y}}_{\ell}\right)$, with
${\mathbf{z}}_{\ell}(u)=\max\left\\{{\mathbf{y}}_{\ell}(u),0\right\\}$, and
where ${\mathbf{z}}_{\ell}(u)$ is the $u$-th component of
${\mathbf{z}}_{\ell}$. Besides the convolutional filter and the pointwise
nonlinear operator, an operation of pooling may be considered to reduce the
computational cost of processing with multiple filters. For multigraph neural
networks, the pooling operator may be defined using the zeroing approach used
for GNNs in [39], where the information is forced to be zero in a subset of
nodes. This is done under the hypothesis that the information in those nodes
is less relevant or redundant. Then, if zeroing is used the pooling operator
is given by $P_{\ell}:{\mathbb{R}}^{N}\rightarrow{\mathbb{R}}^{N}$ with
${\mathbf{x}}_{\ell+1}=P_{\ell}({\mathbf{z}}_{\ell})$ and
${\mathbf{x}}_{\ell+1}(u)=0$ for any $u\in{\mathcal{U}}\subset{\mathcal{V}}$,
and where the subset ${\mathcal{U}}$ could be chosen according to a specific
heuristic or the minimization of a cost function. Notice that for the sake of
simplicity, in this example we considered one filter per layer, but we can
indeed use several filters.

We refer the reader to Appendix B and Appendix C for a discussion on how
convolutional group neural networks and quaternion neural networks can be seen
as particular cases of a generic AlgNN.

######  Remark 4.

The cost of the convolution in the $\ell$-th layer is
$\mathcal{O}(N_{\ell}^{2}m^{K}F_{\ell}G_{\ell})$, where
$N_{\ell}=\text{dim}({\mathcal{M}}_{\ell})$ is the dimension of the vector
space, $m$ is the number of generators in the algebra, $F_{\ell}$ is the
number of input features, $G_{\ell}$ is the number of output features, and $K$
is the order of the polynomial filters. The number of learnable parameters in
each layer is $\mathcal{O}(m^{K}F_{\ell}G_{\ell})$, which does not depend on
the dimension of the vector space. It is important to remark that the
complexity associated with the learnable parameters can be reduced when
adapting pruning algorithms like the one proposed in [16, 19], which reduces
efficiently the number of monomials in an polynomial operator.

##  V Algebraic Perturbation Models

For our discussion we consider perturbations of the generic algebraic signal
model $({\mathcal{A}},{\mathcal{M}},\rho)$ determined by a perturbation of
$\rho$. As discussed in [1], if the realization of algebraic filters in the
algebra is achieved by $\rho$, it is natural to consider that mismatches in
the model occur on $\rho$. In the following definition we state formally the
notion of perturbation in the context of ASP.

######  Definition 7.

(ASP Model Perturbation [1]) Let $({\mathcal{A}},{\mathcal{M}},\rho)$ be an
ASP model with algebra elements generated by $g\in{\mathcal{G}}$ (Definition
1) and recall the definition of the shift operators ${\mathbf{S}}=\rho(g)$
(Definition 2). We say that $({\mathcal{A}},{\mathcal{M}},\tilde{\rho})$ is a
perturbed ASP model if for all $a=p({\mathcal{G}})$ we have that

$\tilde{\rho}(a)=p\big{(}\tilde{\rho}(g)\big{)}=p\big{(}\tilde{\mathcal{S}}\big{)},$
(12)

where $\tilde{\mathcal{S}}$ is a set of perturbed shift operators of the form

${\tilde{\mathbf{S}}}={\mathbf{S}}+{\mathbf{T}}({\mathbf{S}}),$ (13)

for all shift operators ${\mathbf{S}}\in{\mathcal{S}}$.

From Definition 7 we can see that the effect of a perturbation on the
homomorphism $\rho$ is expressed in terms of perturbations of the shift
operators, which at the same time produce a perturbed version of any algebraic
filter. Notice that $\tilde{\rho}$ is a general map that is not necessarily a
homomorphism, but it could be. As shown in [1] this notion of perturbation is
associated to practical scenarios in GSP, WSP and group signal processing.

The notion of stability is tied to the concept of size of a deformation. To
measure the size of those deformations on the space of operators we use norms.
In particular, we use the operator norms induced by a norm associated to
${\mathcal{M}}$. If $\text{dim}({\mathcal{M}})<\infty$ the specifics of the
norm in ${\mathcal{M}}$ are not relevant since all norms are equivalent [43].
If ${\mathcal{M}}$ is infinite dimensional we select a norm that guarantees
that every bounded operator in $\text{End}({\mathcal{M}})$ is Hilbert-Schmidt
and that the induced norm in a direct sum of replicas of ${\mathcal{M}}$
satisfies the maximum property – see Appendix H in the supplementary material.
We remark that this does not affect the structure of any algebraic signal
model $({\mathcal{A}},{\mathcal{M}},\rho)$. We now state the formal definition
of stability considered in our analysis.

######  Definition 8 (Operator Stability [1]).

Given operators $p(\mathbf{S})$ and $p(\tilde{\mathbf{S}})$ defined on the
models $({\mathcal{A}},{\mathcal{M}},\rho)$ and
$({\mathcal{A}},{\mathcal{M}},\tilde{\rho})$ (cf. Definition 7) we say the
operator $p(\mathbf{S})$ is Lipschitz stable if there exist constants
$C_{0},C_{1}>0$ such that

$\left\|p(\mathbf{S})\mathbf{x}-p(\tilde{\mathbf{S}})\mathbf{x}\right\|\leq\\\
\left[C_{0}\sup_{{\mathbf{S}}\in{\mathcal{S}}}\|\mathbf{T}(\mathbf{S})\|+C_{1}\sup_{{\mathbf{S}}\in{\mathcal{S}}}\big{\|}D_{{\mathbf{T}}}({\mathbf{S}})\big{\|}+\mathcal{O}\left(\|\mathbf{T}(\mathbf{S})\|^{2}\right)\right]\big{\|}{\mathbf{x}}\big{\|},$
(14)

for all ${\mathbf{x}}\in{\mathcal{M}}$. In (14)
$D_{{\mathbf{T}}}({\mathbf{S}})$ is the Fréchet derivative of the perturbation
operator ${\mathbf{T}}$.

The right hand side of (14) provides a measure of the deformation produced by
${\mathbf{T}}({\mathbf{S}})$. Then, Definition 8 states that a given operator
is stable to a perturbation ${\mathbf{T}}({\mathbf{S}})$ if the deformation
induced in the operator is proportional to the size of the deformation. For
our discussion we consider the perturbation model given by

${\mathbf{T}}({\mathbf{S}}_{i})={\mathbf{T}}_{0,i}+{\mathbf{T}}_{1,i}{\mathbf{S}}_{i},$
(15)

which is composed of an absolute or additive perturbation ${\mathbf{T}}_{i,0}$
and a relative or multiplicative perturbation
${\mathbf{T}}_{i,1}{\mathbf{S}}_{i}$. The family of perturbations is ruled by
the condition

$\left\|{\mathbf{T}}_{i,r}\right\|_{F}\leq\delta\left\|{\mathbf{T}}_{i,r}\right\|,$
(16)

where $\delta>0$. This is, the Frobenius norm of the perturbation operators is
bounded by a scalar factor of the operator norm.

######  Remark 5.

The notion of stability discussed here is stability to deformations on the
_operator_ space. This is the same kind of deformation that is studied in [1]
for commutative algebraic filters and in [2] for graph filters. It is not the
same as the _domain_ deformations studied in [4]. Although different in
principle, both notions can be related [1]. For instance, a small deformation
of the time axis for time signals implies that signals will not have
translation symmetry but instead quasi translation symmetry. This can be seen
as a perturbation of the time delay operator. It is important to point out
that in all scenarios – domain and operator deformations – a stable operators
must have changes that are proportional in size to the given deformations.
This is done by using metrics used to measure the size of the diffeomorphisms
involved in each case [44, 45, 46].

##  VI Stability Theorems

In the following definitions we state properties of algebraic filters in
${\mathcal{A}}$ in terms of their spectral representations in
$\tilde{{\mathcal{A}}}$ (see Def. 6). This will be used in the derivation of
the stability results. We start introducing the notions of Lipschitz and
integral Lipschitz filters.

######  Definition 9.

Let $({\mathcal{A}},{\mathcal{M}},\rho)$ be an ASM where ${\mathcal{A}}$ has
generators $\\{g_{i}\\}_{i=1}^{m}$. Let $\tilde{{\mathcal{A}}}$ the algebra of
matrices containing the spectral representations of the elements in
${\mathcal{A}}$ (see Def. 6). We say that
$p(g_{1},\ldots,g_{m})\in{\mathcal{A}}$ is $L_{0}$-Lipschitz if there exists
$L_{0}>0$ such that

$\left\|p(x_{1},\ldots,x_{m})-p(\tilde{x}_{1},\ldots,\tilde{x}_{m})\right\|\\\
\leq
L_{0}\left\|(x_{1},\ldots,x_{m})-(\tilde{x}_{1},\ldots,\tilde{x}_{m})\right\|,$
(17)

for all $x_{i},\tilde{x}_{i}$ and where
$p(x_{1},\ldots,x_{m})\in\tilde{{\mathcal{A}}}$. Additionally, it is said that
$p(g_{1},\ldots,g_{m})\in{\mathcal{A}}$ is $L_{1}$-integral Lipschitz if there
exists $L_{1}>0$ such that

$\left\|{\mathbf{D}}_{p|x_{i}}(x_{1},\ldots,x_{m})\left\\{\left(\cdot\right)x_{i}\right\\}\right\|\leq
L_{1},$ (18)

$\text{for all}\leavevmode\nobreak\ x_{i}$, where
${\mathbf{D}}_{p|x_{i}}(x_{1},\ldots,x_{m})$ is the partial Fréchet derivative
of $p(x_{1},\ldots,x_{m})\in\tilde{{\mathcal{A}}}$, and $\|\cdot\|$ is the
operator norm.

From now on we denote the set of algebraic Lipschitz filters by
${\mathcal{A}}_{L_{0}}$ and the set of algebraic integral Lipschitz filters by
${\mathcal{A}}_{L_{1}}$.

Now we introduce the first stability theorem for algebraic signal models with
non commutative algebras with multiple generators.

######  Theorem 2.

Let $(\mathcal{A},{\mathcal{M}},\rho)$ be an ASM where ${\mathcal{A}}$ is
generated by $\\{g_{i}\\}_{i=1}^{m}$ and let
$\rho(g_{i})={\mathbf{S}}_{i}\in\text{End}(\mathcal{M})$ for all $i$. Let
$\tilde{\rho}(g_{i})=\tilde{{\mathbf{S}}}_{i}\in\text{End}(\mathcal{M})$,
where $({\mathcal{A}},\mathcal{M},\tilde{\rho})$ is a perturbed version of
$({\mathcal{A}},\mathcal{M},\rho)$ and $\tilde{{\mathbf{S}}}_{i}$ is related
with ${\mathbf{S}}_{i}$ by the perturbation model in (13). Then, for any
$p\in\mathcal{A}$ we have

$\left\|p({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\mathbf{x}-p(\tilde{{\mathbf{S}}}_{1},\ldots,\tilde{{\mathbf{S}}}_{m})\mathbf{x}\right\|\leq\\\
\|\mathbf{x}\|\sum_{i=1}^{m}\left(\left\|D_{p|\mathbf{S}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\mathbf{T}(\mathbf{S}_{i})\right\|+\mathcal{O}\left(\|\mathbf{T}(\mathbf{S}_{i})\|^{2}\right)\right),$
(19)

where $D_{p|\mathbf{S}_{i}}(\mathbf{S})$ is the partial Fréchet derivative of
$p$ on $\mathbf{S}_{i}$.

###### Proof.

See Appendix D-A. ∎

It is worth pointing out that in (19) the upper bound adds the contributions
of the deformation in each direction ${\mathbf{S}}_{i}$ given by the shift
operators. Each individual contribution associated to the perturbation of
${\mathbf{S}}_{i}$ is determined by the Fréchet derivative of the filters
acting on the perturbation, and this is true no matter what function
${\mathbf{T}}({\mathbf{S}}_{i})$ is being considered.

In the following theorems we provide the basic stability result for algebraic
filters showing how with a restriction of the filters shaped by the functional
form of the right hand side of (19) leads to stability.

######  Theorem 3.

Let $(\mathcal{A},{\mathcal{M}},\rho)$ be an ASM where ${\mathcal{A}}$ is a
non commutative algebra with $m$ generators $\\{g_{i}\\}_{i=1}^{m}$. Let
$({\mathcal{A}},\mathcal{M},\tilde{\rho})$ a perturbed version of
$({\mathcal{A}},\mathcal{M},\rho)$ by means of the perturbation model in (15).
Then, if $p\in\mathcal{A}_{L_{0}}\cap\mathcal{A}_{L_{1}}$ it holds that

$\left\|D_{p|\mathbf{S}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m}){\mathbf{T}}({\mathbf{S}}_{i})\right\|\leq\delta
L_{0}\sup_{{\mathbf{S}}_{i}\in\mathcal{S}}\|{\mathbf{T}}({\mathbf{S}}_{i})\|\\\
+\delta
L_{1}\sup_{{\mathbf{S}}_{i}\in\mathcal{S}}\|D_{{\mathbf{T}}}({\mathbf{S}}_{i})\|.$
(20)

###### Proof.

See Appendix D-B. ∎

From Theorems 2 and 3 we can state the stability results for filters in
algebraic models with multiple generators in the following corollary.

######  Corollary 1.

Let $({\mathcal{A}},{\mathcal{M}},\rho)$ be a non commutative ASM where
${\mathcal{A}}$ has generators $\\{g_{i}\\}_{i=1}^{m}$. Let
$({\mathcal{A}},{\mathcal{M}},\tilde{\rho})$ be a perturbed version of
$({\mathcal{A}},{\mathcal{M}},\rho)$ associated to the perturbation model in
(15). If $p\in\mathcal{A}_{L_{0}}\cap\mathcal{A}_{L_{1}}\subset{\mathcal{A}}$,
the operator $p({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})$ is stable in the
sense of Definition 8 with $C_{0}=m\delta L_{0}$ and $C_{1}=m\delta L_{1}$.

###### Proof.

Replacing (20) from Theorem 3 into (19) from Theorem 2 and organazing the
terms.∎

Notice that the stability property of the algebraic filters comes at the
expense of their selectivity. A restriction on the subsets of the algebra
naturally limits the type of frequency representations. Additionally, it is
worth pointing out that it is the functional form of
${\mathbf{T}}({\mathbf{S}}_{i})$ the one feature that dictates the types of
restrictions imposed on the algebra, necessary to guarantee stability.

######  Remark 6.

We remark that the non commutativity imposes multidimensional frequency
representations, but even in this scenario we have a trade-off between
stability and selectivity. We could therefore say that although the non
commutativity of the algebra does not change the fact that filters can be
stable and selective, the multidimensional nature of the frequency
representations leads to restrictions that can be captured only by means of
matrix algebras. This in itself implies that two filters with identical
functional form – symbolic expression – do not have the same stability
properties if one of them is associated to a commutative model and the other
to a non commutative one.

###  VI-A Stability of non commutative Algebraic Neural Networks

In this section we extend the stability results obtained for algebraic filters
to operators representing AlgNNs. Let
$\left\\{\left(\mathcal{A}_{\ell},\mathcal{M}_{\ell},\rho_{\ell};\sigma_{\ell}\right)\right\\}_{\ell=1}^{L}$
be an AlgNN with $L$ layers, and whose perturbed version is represented as
$\left\\{\left(\mathcal{A}_{\ell},\mathcal{M}_{\ell},\tilde{\rho}_{\ell};\sigma_{\ell}\right)\right\\}_{\ell=1}^{L}$.
We start stating formally how the algebraic operators are affected by the
functions that map information between layers.

######  Theorem 4.

Let
$\left\\{\left(\mathcal{A}_{\ell},\mathcal{M}_{\ell},\rho_{\ell};\sigma_{\ell}\right)\right\\}_{\ell=1}^{L}$
be an algebraic neural network and
$\left\\{\left(\mathcal{A}_{\ell},\mathcal{M}_{\ell},\tilde{\rho}_{\ell};\sigma_{\ell}\right)\right\\}_{\ell=1}^{L}$
its perturbed version by means of the perturbation model in (15). We consider
one feature per layer and non commutative algebras $\mathcal{A}_{\ell}$ with
$m$ generators. 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 $\ell$-th mapping operators of the AlgNN and its perturbed
version, it follows that

$\left\|\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\|\leq\\\
C_{\ell}\delta\|\mathbf{x}_{\ell-1}\|m\left(L_{0}^{(\ell)}\sup_{{\mathbf{S}}_{i,\ell}}\|{\mathbf{T}}^{(\ell)}({\mathbf{S}}_{i,\ell})\|\right.\\\
\left.+L_{1}^{(\ell)}\sup_{{\mathbf{S}}_{i,\ell}}\|D_{\mathbf{T^{(\ell)}}}({\mathbf{S}}_{i,\ell})\|\right)$
(21)

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$.

###### Proof.

See Appendix J-A1 in the supplementary material. ∎

Theorem 4 highlights the role of the maps $\sigma_{\ell}$ in the stability of
the algebraic operators. In particular, we can see that the effect of such
functions is only to scale the stability bound of the algebraic operators. If
$C_{\ell}=1$, the stability bounds for the operators in the layers of the
AlgNN are identical to the stability bounds for the algebraic filters.
However, it is important to remark that the discriminability of the layer
operators in the AlgNN is enriched by the pointwise nonlinearities.

Notice that the resemblance between the results provided in Theorem 4 and
previous results for commutative architectures is rooted in the fact that the
same notion of stability and the same type of perturbations are considered in
both scenarios. However, there is a substantial difference between the
stability bounds. In the non commutative scenario we have that the
restrictions on the subsets of the algebra are dictated by matrix functions
and not scalar maps. Additionally, as can be seen in the proof of the theorems
above, the estimate of the upper bounds on the change of the operators
requires a more general representation of spectral responses not given by the
traditional spectral theorem.

Now we are ready to state the stability theorem for a general AlgNN.

######  Theorem 5.

Let
$\left\\{\left(\mathcal{A}_{\ell},\mathcal{M}_{\ell},\rho_{\ell};\sigma_{\ell}\right)\right\\}_{\ell=1}^{L}$
be an algebraic neural network and
$\left\\{\left(\mathcal{A}_{\ell},\mathcal{M}_{\ell},\tilde{\rho}_{\ell};\sigma_{\ell}\right)\right\\}_{\ell=1}^{L}$
its perturbed version by means of the perturbation model in (15). We consider
one feature per layer and non commutative algebras $\mathcal{A}_{\ell}$ with
$m$ generators. 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 operator and its perturbed version, it follows that

$\left\|\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\|\\\
\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\|\mathbf{x}\right\|,$
(22)

where $C_{\ell}$ is the Lipschitz constant of $\sigma_{\ell}$,
$\|\rho_{\ell}(a)\|\leq B_{\ell}\leavevmode\nobreak\ \text{for
all}\leavevmode\nobreak\ a\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

$\boldsymbol{\Delta}_{\ell}=\delta
m\left(L_{0}^{(\ell)}\sup_{{\mathbf{S}}_{i,\ell}}\|{\mathbf{T}}^{(\ell)}({\mathbf{S}}_{i,\ell})\|\right.\left.+L_{1}^{(\ell)}\sup_{{\mathbf{S}}_{i,\ell}}\|D_{\mathbf{T^{(\ell)}}}({\mathbf{S}}_{i,\ell})\|\right),$
(23)

with the index $(\ell)$ indicating quantities and constants associated to the
layer $\ell$.

###### Proof.

See Section J-A2 in the supplementary material. ∎

This final result highlights that the stability of an AlgNN is inherited from
the stability properties associated to the operators in each layer. Each layer
of the AlgNN contributes to increase the size of the stability constants, and
we can observe that with the appropriate normalization of $C_{\ell}$ and
$B_{\ell}$, we obtain a stability bound that is in essence the same derived
for filters and for the mapping operators in each layer. However, the
discriminability power associated to the AlgNN is by far larger than the one
related with the operators in the layers and the filters. This is a
consequence of the pointwise nonlinearities mapping information between
layers. It is important to highlight that the $\sigma_{\ell}$ that map
information between layers are identical for commutative and non commutative
architectures. However, they way they redistribute spectral information is
tied to the nature of the spectral representations. While in commutative
architectures such redistribution is always done between spaces of the same
dimension – one dimensional representations –, in non commutative AlgNNs the
redistribution of spectral information occurs in general between spaces of
different dimension.

##  VII Numerical Experiments

In order to provide numerical evidence of the stability results derived for
AlgNN with non commutative algebras, we consider two architectures: Multigraph
Neural Networks (MultiGNN) and Quaternion Convolutional Networks (QCN). For
the first, we propose an architecture for learning on multigraphs and show
that if the learned filters are penalized to be Integral Lipschitz, the model
is much more resilient to perturbations enforced on the graph. This
architecture is used for a rating prediction task on the MovieLens-100K
dataset [47]. We then apply the quaternion architecture introduced in [14] to
an image classification task using the MNIST dataset [48] and observe that
model performance quickly diminishes when additive and relative perturbations
are applied to the quaternion filters. As we will elaborate in Subsection
VII-B, due to the cyclic nature of the generators in the quaternion algebra,
no filter in (26) can be integral Lipschitz. As a consequence, no filter can
be selected to mitigate the perturbations that severely affect the performance
of the QCN.

(a)

(b)

Figure 6: Stability to perturbations caused by estimation error for the Movie
Recommendation problem. We demonstrate in (a) the difference in evaluation
measure as the ratio of training samples is increased, where the penalized
MultiGNN maintains the smallest difference for all ratios. This trend is also
held in (b) where we showcase the difference in the output of the
convolutional layer. Then, MultiGNN with IL filters provide more stable and
consistent rating predictions with respect to estimation error perturbations.

###  VII-A Multigraph Neural Networks

We first consider the application of the Multigraph Neural Network as a movie
recommendation system [49]. Using the MovieLens-100K dataset containing
100,000 ratings from 943 users, we model the system as a multigraph. Nodes
represent movies which are connected by two classes of edges: one measuring
the rating similarity and the other measuring the genre similarity. The
signals supported by this multigraph are user ratings, with the task of
inferring unseen ratings based on ratings given to other movies. On this task,
we carry out stability experiments which measure robustness to estimation
error from the training set. We find that Multigraph Neural Networks that are
penalized to employ integral Lipschitz filters are most stable to
perturbations.

Multigraph Formation: From the 1582 movies available, we use the 200 with the
highest number of ratings. The movie _Toy Story_ is used as the target movie
for predicted ratings. Of users who rated this movie, we define a node signal
as the rating 1-5 they gave to each movie or 0 if they did not rate it. These
signals are divided with a 90%/10% train/test split. Using the training set,
we add rating edge weights between movies using the Pearson correlation
coefficient and genre edge weights through the Jaccard similarity of
overlapping genres. To appeal to sparsity, the 20 highest weighted edges of
both edge classes are kept for each node, and the rest of the edges are
removed.

Stability: In contrast with the scenario of the quaternion algebras,
multigraph signal models are associated with a very rich algebra: the regular
algebra of polynomials with multiple independent variables. This guarantees
the existence of subclasses of filters that are Lipschitz and integral
Lipschitz. Therefore, it is possible to learn filters – Lipschitz and integral
Lipschitz – that can mitigate the effect of deformations on the shift
operators. We illustrate this in our experiments analyzing the magnitude of
change in the filters and operators networks when subjected to deformations.

Architectures and Training: In our stability experiments, we consider three
architectures. The first, MultiFilter, is a learned linear multigraph filter.
We also train a Multigraph Neural Network (MultiGNN) and another which is
regularized by an estimate of the filter’s integral Lipschitz constant
(MultiGNN IL), seeking to train filters that are integral Lipschitz. All
architectures employ a convolutional layer with 3 filter taps and 64 output
features, a ReLU nonlinearity function, and a local linear readout layer
mapping the output rating to a single scalar rating estimate. We minimize the
smooth L1 loss using an ADAM optimizer for 40 epochs, and evaluate using the
root mean squared error (RMSE). Results are reported as the average and
standard deviation of performance across 5 random splits of the dataset.

Estimation Error Experiment: To simulate perturbations as a result of
estimation error, we first train each architecture on 90% of the training set.
At evaluation, we replace the underlying multigraph support with an estimate
generated by a training set ranging from 10% to 90% of the size of the overall
dataset. First, we find the difference in RMSE between the trained and
evaluated models. The regularized MultiGNN maintains the smallest difference
for each training set size as demonstrated in Figure 6(a). We also compare the
norm of the difference in output of the convolutional layer. As can be seen in
Figure 6(b), the MultiGNN with integral Lipschitz filters is approximately an
order of magnitude more stable than the unregularized MultiGNN.

###  VII-B Quaternion Convolutional Networks

We now consider a non commutative AlgNN proposed in literature: Quaternion
Convolutional Networks. Quaternions have shown to be useful tools for modeling
spatial transformations with diverse applications to computer graphics,
quantum mechanics, and signal processing. Gaudet and Maida [14] propose a deep
convolutional architecture for learning through use of quaternion
convolutions, quaternion batch-normalization, and a quaternion weight
initialization scheme. We use this architecture as the basis for our second
stability experiment.

Architecture: Quaternion convolutions occur through convolving a quaternion
filter matrix $\mathbf{W}=\mathbf{A}+i\mathbf{B}+j\mathbf{C}+k\mathbf{D}$ with
a quaternion vector
$\mathbf{h}=\mathbf{w}+i\mathbf{x}+j\mathbf{y}+k\mathbf{z}$ according to (29)
– see Appendix C. In forming a Quaternion Convolutional Network, we compose
layers of quaternion convolutions, applying a nonlinear function to the output
of each layer. We then append quaternion dense layers, which are traditional
fully-connected neural networks applied component-wise. The resultant
embedding is shaped into the desired output shape using a readout layer. Our
particular architecture is composed of three quaternion convolutional layers
(with 32 channels each), two quaternion dense layers (with 82 and 48), and a
readout layer, using the ReLU function as the nonlinearity.

Stability: We recall that one of the requirements to guarantee stability to
the perturbations in (15) is the existence of subclasses of filters in (26)
that are Lipschitz and integral Lipschitz. Although in principle any filter in
(26) can be written by considering arbitrary powers on the generators of the
algebra, due to the cyclic nature of
$\\{\boldsymbol{1},\boldsymbol{i},\boldsymbol{j},\boldsymbol{k}\\}$, such
representations reduce to a linear combination of the generators. Therefore,
the subclass of integral Lipschitz filters in (26) is empty. Additionally,
while the filters in (26) can still be Lipschitz, the Lipschitz constants will
be large. Thus, the performance of the QCN architecture – as defined in [14] –
is expected to be severely affected by perturbations, without the potential to
design stable filters. We demonstrate this vulnerability in our numerical
experiments.

Synthetic Experiment: To corroborate these theoretical results, we first train
an unperturbed architecture. Then, during evaluation, we inject noise in the
first quaternion convolutional layer by perturbing $\mathbf{A}$, a layer of
the quaternion filter matrix $\mathbf{W}_{1}$. More precisely, we first
consider additive perturbations by replacing $\mathbf{A}$ with
$\tilde{\mathbf{A}}_{1}=\mathbf{A}+\mathbf{T}_{1}$, where $\mathbf{T}_{1}$ has
its components sampled from the uniform distribution
$[-\varepsilon_{1},\varepsilon_{1}]$. We also observe relative perturbations
by substituting $\mathbf{A}$ with
$\tilde{\mathbf{A}}_{2}=\mathbf{A}+\mathbf{T}_{2}\mathbf{A}$, where
$\mathbf{T}_{2}$ has its values drawn from the uniform distribution
$[\varepsilon_{2},\varepsilon_{2}]$.

| Additive Perturbation |  | Relative Perturbation
---|---|---|---
$\varepsilon_{1}$ | First Con. | Last Con. | Acc. | $\varepsilon_{2}$ | First Con. | Last Con. | Acc.
0.005 | 27,813 | 4,415 | 0.991 | 0.05 | 24,986 | 6,316 | 0.991
0.01 | 55,467 | 9,542 | 0.991 | 0.1 | 49,961 | 13,100 | 0.990
0.05 | 289,610 | 78,092 | 0.968 | 0.5 | 256,754 | 79,819 | 0.956
0.1 | 604,553 | 213,046 | 0.772 | 1 | 533,462 | 158,565 | 0.825
0.5 | 3,303,611 | 234,771 | 0.222 | 5 | 2,975,365 | 1,361,751 | 0.321
Table 1. Stability of Quaternion Convolutional Networks to additive and
relative perturbation.

We apply this architecture to the image classification task using the MNIST
dataset [48]. Following [50], we form the quaternion input for each pixel by
repeating the grayscale value along each of the three imaginary components and
zeroing out the real component. We train the unperturbed architecture for 10
epochs using a batch size of 500. After training, we first evaluate the
performance of the unperturbed architecture on the test samples. We then apply
additive and multiplicative perturbation to the first quaternion convolutional
layer for various values of $\varepsilon_{1}$ and $\varepsilon_{2}$, and
observe the new accuracy (Acc.) and norm of the difference between the
perturbed and unperturbed architectures of the output from the first
quaternion convolutional layer (First Con.) and the last quaternion
convolutional layer (Last Con.). Average results across 10 noise samplings are
reported in Table 1.

Our results find that Quaternion Convolutional Networks are particularly
sensitive to small additive and relative perturbations applied in the first
layer, as performance quickly diminishes as $\varepsilon_{1}$ and
$\varepsilon_{2}$ increases. These results indicate that an architecture
trained to strictly minimize loss not result in stable filters.

##  VIII Discussion and Conclusions

We analyzed non commutative algebraic neural networks (AlgNNs) providing
general stability results applicable to convolutional architectures such as
multigraph-CNNs, quaternion-CNNs, quiver-CNNs, and Group-NNs (non commutative)
among others. We have shown that non commutative convolutional neural networks
can be stable. This property has an analogous form to the stability bounds of
commutative signal models but with substantial differences. As we showed in
previous sections, the restrictions necessary to guarantee stability are
defined by algebras of matrices, which highlights the fact that the stability
of filters in non commutative signal models is not a trivial extension of the
stability results for commutative models. This is a consequence of the fact
that spectral representations in non commutative signal models are given by
matrix and not scalar functions.

As shown in Theorem 2, the size of the deformation of a filter operator
$p({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})$ is bounded by the norm of the
Fréchet derivative of $p$ acting on ${\mathbf{T}}({\mathbf{S}}_{i})$. This
holds for any arbitrary perturbation ${\mathbf{T}}({\mathbf{S}}_{i})$
indicated in (13) and it highlights the versatility of the deformation model
for the representation of a large variety of perturbations.

Notice that non commutative models naturally embed properties of those signal
frameworks where there are not shift operator invariances, i.e.
${\mathbf{S}}_{i}p({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\neq
p({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m}){\mathbf{S}}_{i}$. The stability
results derived in this paper show that despite the lack of this invariance,
stability is possible and it comes at a price on selectivity and with stronger
restrictions on the spectral responses of the filters.

We point out that our results can be extended to more elaborated versions of
${\mathbf{T}}({\mathbf{S}}_{i})$. This opens up interesting research
applications where ${\mathbf{T}}({\mathbf{S}}_{i})$ can be tailored to more
specific scenarios where stable filters could be designed for
${\mathbf{T}}({\mathbf{S}}_{i})=\sum_{r}{\mathbf{T}}_{r}{\mathbf{S}}_{i}^{r}$,
and the restrictions imposed on the algebra will be determined by the
properties of the operators ${\mathbf{T}}_{r}$.

##  Appendix A Basic Algebra Concepts

Algebra: We say that the vector space ${\mathcal{A}}$ is an algebra if
${\mathcal{A}}$ is equipped with a product. This is, there is an operation of
product between two elements of ${\mathcal{A}}$ that results in another
element of ${\mathcal{A}}$ [35]. Let us denote by $ab$ the product between $a$
and $b$ with $a,b\in{\mathcal{A}}$. Then, we say that the algebra
${\mathcal{A}}$ is unital if there exists a multiplicative unit element
$\mathbf{1}\in{\mathcal{A}}$ such that $a\mathbf{1}=\mathbf{1}a=a$ for all
$a\in{\mathcal{A}}$. Typical examples of algebras are:

  * •

$\mathbb{C}$: The set of complex numbers is a vector space over $\mathbb{C}$
itself. Using as algebra product in $\mathbb{C}$ the ordinary product between
complex numbers we can see $\mathbb{C}$ as an algebra.

  * •

$\mathbb{C}[t]$: The set of polynomials with coefficients in $\mathbb{C}$ and
independent variable $t$ is a vector space over $\mathbb{C}$. Using as algebra
product the ordinary product between polynomials we can see $\mathbb{C}[t]$ as
an algebra.

  * •

$M_{r\times r}(\mathbb{C})$: The space of matrices of size $r\times r$ with
entries in $\mathbb{C}$ is a vector space over $\mathbb{C}$. Using as algebra
product the ordinary product between matrices we can see $M_{r\times
r}(\mathbb{C})$ as an algebra.

It is important to remark that the formal definition of a product is given by
a bilinear map. We say that the map
$m_{{\mathcal{A}}}:{\mathcal{A}}\times{\mathcal{A}}\to{\mathcal{A}}$ is
bilinear if $m_{{\mathcal{A}}}$ is linear on each argument. Then, talking
about an operation that behaves like a product we talk about a bilinear map
that assigns the value of the given product [35].

Endomorphisms: Let ${\mathcal{M}}$ be a vector space and let
$\text{End}({\mathcal{M}})$ be the space of linear maps from ${\mathcal{M}}$
onto itself. If $\text{dim}({\mathcal{M}})=n$, then
$\text{End}({\mathcal{M}})$ is isomorphic to the space of matrices of size
$n\times n$. This is, whenever the dimension of ${\mathcal{M}}$ is finite we
can think about $\text{End}({\mathcal{M}})$ as a space of matrices.

##  Appendix B Group Signal Processing and Group Neural Networks

Classical convolutional information processing on locally compact groups is
given by the algebraic model $({\mathcal{A}}_{G},{\mathcal{M}},\rho)$, where
${\mathcal{M}}=\left\\{\sum_{g\in
G}{\mathbf{x}}(g)g,\quad{\mathbf{x}}(g)\in\mathbb{C}\right\\}$ is the set of
complex valued functions defined on the group $G$. The algebra
${\mathcal{A}}_{G}$ is the group algebra, which is given by
${\mathcal{A}}_{G}={\mathcal{M}}$. The homomorphism $\rho$ is given by
$\rho({\mathbf{a}})=L_{{\mathbf{a}}}$ with
$L_{{\mathbf{a}}}{\mathbf{x}}={\mathbf{a}}{\mathbf{x}}$. Then, the convolution
between a group filter and a group signal takes the form

$\rho\left(\sum_{g\in G}{\mathbf{a}}(g)g\right){\mathbf{x}}=\sum_{g\in
G}{\mathbf{a}}(g)g{\mathbf{x}}=\sum_{g\in G}\sum_{h\in
G}{\mathbf{a}}(g){\mathbf{x}}(h)gh.$ (24)

Notice that the representation of the convolution can also be re written
making $u=gh$, which leads to [36]

$\sum_{g,h\in G}{\mathbf{a}}(g){\mathbf{x}}(h)gh=\sum_{u,h\in
G}{\mathbf{a}}(uh^{-1}){\mathbf{x}}(h)u.$ (25)

We emphasize that the role of this example is to show that classical group
convolutions – as considered in [17] – are one particular case of algebraic
convolutions. What is more, it is one particular case of the many algebraic
convolutions that can be associated to the same group algebra. This is, given
a non commutative group $G$ and its associated group algebra
${\mathcal{A}}_{G}$, it is possible to define different convolutional models
$({\mathcal{A}}_{G},{\mathcal{M}},\rho)$ when different choices of
${\mathcal{M}}$ and $\rho$ are made [15]. This implies that it is possible to
leverage the symmetries of the group under a convolutional signal model in
different ways and in vector spaces of arbitrary dimension.

Notice that when convolutions are computed as in [17], a lifting process is
required to map the signals – usually not defined on the group – on the group
$G$. If $G$ is continuous – a Lie group – computing the group convolution
requires the approximation of a Haar integral that generalizes (25) to locally
compact groups [34]. In [15] the computing of the convolutions – including
convolutions with Lie groups – is performed considering an ASM
$({\mathcal{A}}_{G},{\mathcal{M}},\rho)$, where ${\mathcal{A}}_{G}$ is the
group algebra associated to the group $G$, ${\mathcal{M}}$ is a vector space
of arbitrary dimension and $\rho$ is the homomorphism that will translate
approximations of elements in ${\mathcal{A}}_{G}$ into concrete operators. In
[15] lifting is not necessary since the implementation of the filters in
${\mathcal{A}}_{G}$ is carried out by means of an efficient approximation of
$\rho$.

A group convolutional neural network is a stacked layered structure where
information in each layer is processed using group convolutions. With
convolutions performed as in [15], in each layer of the architecture
information is processed according to the triplet
$({\mathcal{A}}_{G},{\mathcal{M}}_{\ell},\rho_{\ell})$. The information is
mapped between layers by means of $\sigma_{\ell}=P_{\ell}\eta_{\ell}$ where
$\eta_{\ell}$ is a pointwise nonlinearity and $P_{\ell}$ is a dimensionality
reduction operator. When dealing with discrete scenarios we have
$\text{dim}({\mathcal{M}})<\infty$ and the action of $\eta_{\ell}$ is
expressed in terms of a canonical basis in
$\mathbb{R}^{\text{dim}({\mathcal{M}})}$.

Notice that the processing of information with [17] can be seen as a
particular case of the convolutional architectures proposed in [15] with the
advantage that in [15] there is no need for the discretization of the Haar
integral since the implementation of the filters relies on building an
accurate approximation of $\rho_{\ell}$ in
$({\mathcal{A}}_{G},{\mathcal{M}}_{\ell},\rho_{\ell})$.

##  Appendix C Quaternion Signal Processing and Quaternion Neural Networks

Let ${\mathcal{Q}}$ be the algebra on ${\mathbb{C}}$ generated by the symbols
$\boldsymbol{i},\boldsymbol{j},\boldsymbol{k}$ and the unit element
$\mathbf{1}$, where
$\boldsymbol{i}^{2}=\boldsymbol{j}^{2}=\boldsymbol{k}^{2}=-\mathbf{1},\quad\boldsymbol{ij}=-\boldsymbol{ji}=\boldsymbol{k},\quad\boldsymbol{jk}=-\boldsymbol{kj}=\boldsymbol{i},\quad\boldsymbol{ki}=-\boldsymbol{ik}=\boldsymbol{j}.$
${\mathcal{Q}}$ is known as the algebra of quaternions [38]. Using
${\mathcal{Q}}$ is possible to build quaternion-based algebras that can be
used to define convolutional models. In [14] such models are build upon the
algebra

${\mathcal{A}}=\left\\{\left.p_{1}(t)\boldsymbol{1}+p_{2}(t)\boldsymbol{i}+p_{3}(t)\boldsymbol{j}+p_{4}(t)\boldsymbol{k}\right|\boldsymbol{1},\boldsymbol{i},\boldsymbol{j},\boldsymbol{k}\in{\mathcal{Q}}\right\\},$
(26)

and where $p_{i}(t)=\sum_{k=0}^{k}h_{i,k}t^{k}$ are elements of the algebra of
polynomials of a single variable. We can see that ${\mathcal{A}}$ is a vector
space when considering the sum component wise on the quaternion elements, and
the multiplication by scalars. The product algebra in ${\mathcal{A}}$ follows
the distribution product of the quaternions. Then, we can define the
quaternion convolutional model $({\mathcal{A}},{\mathcal{M}},\rho)$, where

${\mathcal{M}}=\left\\{\left.{\mathbf{x}}_{1}\boldsymbol{1}+{\mathbf{x}}_{2}\boldsymbol{i}+{\mathbf{x}}_{3}\boldsymbol{j}+{\mathbf{x}}_{4}\boldsymbol{k}\right|\boldsymbol{1},\boldsymbol{i},\boldsymbol{j},\boldsymbol{k}\in{\mathcal{Q}};\quad{\mathbf{x}}_{i}\in{\mathbb{R}}^{N}\right\\},$
(27)

while the homomorphism $\rho$ is given according to

$\rho\left(p_{1}(t)\boldsymbol{1}+p_{2}(t)\boldsymbol{i}+p_{3}(t)\boldsymbol{j}+p_{4}(t)\boldsymbol{k}\right)\\\
=p_{1}({\mathbf{C}})\boldsymbol{1}+p_{2}({\mathbf{C}})\boldsymbol{i}+p_{3}({\mathbf{C}})\boldsymbol{j}+p_{4}({\mathbf{C}})\boldsymbol{k},$
(28)

where $\rho(t)={\mathbf{C}}\in{\mathbb{R}}^{N\times N}$ is the discrete time
delay operator. Then, the convolutional action of a filter
$\rho(f)={\mathbf{F}}$ on $\boldsymbol{u}\in{\mathcal{M}}$ is given by

${\mathbf{F}}\boldsymbol{u}=\left(p_{1}({\mathbf{C}}){\mathbf{x}}_{1}-p_{2}({\mathbf{C}}){\mathbf{x}}_{2}-p_{3}({\mathbf{C}}){\mathbf{x}}_{3}-p_{4}({\mathbf{C}}){\mathbf{x}}_{4}\right)\boldsymbol{1}+\\\
\left(p_{1}({\mathbf{C}}){\mathbf{x}}_{2}-p_{2}({\mathbf{C}}){\mathbf{x}}_{1}-p_{3}({\mathbf{C}}){\mathbf{x}}_{4}-p_{4}({\mathbf{C}}){\mathbf{x}}_{3}\right)\boldsymbol{i}+\\\
\left(p_{1}({\mathbf{C}}){\mathbf{x}}_{3}-p_{2}({\mathbf{C}}){\mathbf{x}}_{4}-p_{3}({\mathbf{C}}){\mathbf{x}}_{1}-p_{4}({\mathbf{C}}){\mathbf{x}}_{2}\right)\boldsymbol{j}+\\\
\left(p_{1}({\mathbf{C}}){\mathbf{x}}_{4}-p_{2}({\mathbf{C}}){\mathbf{x}}_{3}-p_{3}({\mathbf{C}}){\mathbf{x}}_{2}-p_{4}({\mathbf{C}}){\mathbf{x}}_{1}\right)\boldsymbol{k}.$
(29)

Notice that $p_{i}({\mathbf{C}}){\mathbf{x}}_{i}$ is the ordinary Euclidean
convolution in ${\mathbb{R}}^{N}$ between a filter and the signal
${\mathbf{x}}_{i}$. We emphasize that the particular form of (29) is a
consequence of the non commutativity of the filter operators from
${\mathcal{A}}$ in (26), which is inherited from the quaternion algebra.

A quaternion neural network as introduced in [14] is a stacked layered
structure where information is processed in each layer according to (29). The
information processed in the layer $\ell$ has the form
$\boldsymbol{u}_{\ell}={\mathbf{x}}_{1}\boldsymbol{1}+\boldsymbol{i}{\mathbf{x}}_{2}+\boldsymbol{j}{\mathbf{x}}_{3}+\boldsymbol{k}{\mathbf{x}}_{4}$,
while the convolutional filters are of the form
${\mathbf{F}}={\mathbf{A}}\boldsymbol{1}+\boldsymbol{i}{\mathbf{B}}+\boldsymbol{j}{\mathbf{C}}+\boldsymbol{k}{\mathbf{D}}$.
The action of the filter ${\mathbf{F}}$ on the signal $\boldsymbol{u}$ follows
the ordinary product between quaternions considering that
${\mathbf{A}},{\mathbf{B}},{\mathbf{C}},{\mathbf{D}}$ act as Euclidean
convolution operators on ${\mathbf{x}}_{i}$ as indicated in (29). Then, given
an input $\boldsymbol{u}_{\ell}$ to the $\ell$-th layer we leverage the
symmetries given by the quaternion algebra to obtain
${\mathbf{y}}_{\ell}={\mathbf{F}}\boldsymbol{u}_{\ell}$. After this a point-
wise nonlinear operator $\eta_{\ell}$ is applied to obtain
${\mathbf{z}}_{\ell}=\eta_{\ell}({\mathbf{y}}_{\ell})$ with
${\mathbf{z}}_{\ell}(u)=\max\left\\{{\mathbf{y}}_{\ell}(u),0\right\\}$, and
where ${\mathbf{z}}_{\ell}(u)$ is the $u$-th component of
${\mathbf{z}}_{\ell}$. A conventional pooling operator
$P_{\ell}:{\mathbb{R}}^{N}\rightarrow{\mathbb{R}}^{M}$ performs traditional
sampling on a quaternion signal according to
$P_{\ell}(\boldsymbol{u}_{\ell})=\boldsymbol{1}P_{\ell}({\mathbf{x}}_{1})+\boldsymbol{i}P_{\ell}({\mathbf{x}}_{2})+\boldsymbol{j}P_{\ell}({\mathbf{x}}_{3})+\boldsymbol{k}P_{\ell}({\mathbf{x}}_{4})$.

##  Appendix D Proof of Theorems

###  D-A Proof of Theorem 2

###### Proof.

To simplify notation we use
${\mathbf{S}}=({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})$, and we start taking
into account the definition of the Fréchet derivative

$p(\mathbf{S}+\boldsymbol{\xi})=p(\mathbf{S})+D_{p}(\mathbf{S})\left\\{\boldsymbol{\xi}\right\\}+o(\|\boldsymbol{\xi}\|).$
(30)

Taking the norm on both sides of (30) and applying the triangle inequality it
follows that

$\left\|p(\mathbf{S}+\boldsymbol{\xi})-p(\mathbf{S})\right\|\leq\left\|D_{p}(\mathbf{S})\left\\{\boldsymbol{\xi}\right\\}\right\|+\mathcal{O}\left(\|\boldsymbol{\xi}\|^{2}\right)$
(31)

for all
$\boldsymbol{\xi}=(\boldsymbol{\xi}_{1},\ldots,\boldsymbol{\xi}_{m})\in\text{End}(\mathcal{M})^{m}$.
Now, we use the properties of the total derivative (see [51] pages 69-70) to
obtain
$\|D_{p}(\mathbf{S})\left\\{\boldsymbol{\xi}\right\\}\|\leq\sum_{i=1}^{m}\left\|D_{p|\mathbf{S}_{i}}(\mathbf{S})\left\\{\boldsymbol{\xi}_{i}\right\\}\right\|,$
which leads to

$\left\|p(\mathbf{S}+\boldsymbol{\xi})-p(\mathbf{S})\right\|\leq\sum_{i=1}^{m}\left\|D_{p|\mathbf{S}_{i}}(\mathbf{S})\left\\{\boldsymbol{\xi}_{i}\right\\}\right\|+\mathcal{O}\left(\|\boldsymbol{\xi}\|^{2}\right).$
(32)

Then, taking into account that

$\left\|p(\mathbf{S}+\boldsymbol{\xi})\mathbf{x}-p(\mathbf{S})\mathbf{x}\right\|\leq\|\mathbf{x}\|\left\|p(\mathbf{S}+\boldsymbol{\xi})-p(\mathbf{S})\right\|$
(33)

we have

$\left\|p(\mathbf{S}+\boldsymbol{\xi}){\mathbf{x}}-p(\mathbf{S}){\mathbf{x}}\right\|\\\
\leq\left\|{\mathbf{x}}\right\|\left(\sum_{i=1}^{m}\left\|D_{p|\mathbf{S}_{i}}(\mathbf{S})\left\\{\boldsymbol{\xi}_{i}\right\\}\right\|+\mathcal{O}\left(\|\boldsymbol{\xi}\|^{2}\right)\right),$
(34)

and selecting $\boldsymbol{\xi}_{i}=\mathbf{T}(\mathbf{S}_{i})$ we complete
the proof. ∎

###  D-B Proof of Theorem 3

###### Proof.

We start evaluating the expression ${\mathbf{D}}_{p|{\mathbf{S}}_{i}}$ in
terms of the specific form of the perturbation to obtain

${\mathbf{D}}_{p|{\mathbf{S}}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\left\\{{\mathbf{T}}({\mathbf{S}}_{i})\right\\}={\mathbf{D}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\left\\{{\mathbf{T}}_{i,0}\right\\}\\\
+{\mathbf{D}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\left\\{{\mathbf{T}}_{i,1}{\mathbf{S}}_{i}\right\\}.$
(35)

Taking the norm in (35) and using the triangular inequality we have

$\left\|{\mathbf{D}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\left\\{{\mathbf{T}}({\mathbf{S}}_{i})\right\\}\right\|\leq\\\
\left\|{\mathbf{D}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\left\\{{\mathbf{T}}_{i,0}\right\\}\right\|\\\
+\left\|{\mathbf{D}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\left\\{{\mathbf{T}}_{i,1}{\mathbf{S}}_{i}\right\\}\right\|.$
(36)

Now we analyze each term in (36). First, we start that the Fréchet derivative
of $p({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})$ with respect to
${\mathbf{S}}_{i}$ acting on ${\mathbf{T}}_{0,i}$ is a linear operator
${\mathbf{D}}_{p|{\mathbf{S}}_{i}}$ from $\text{End}(\mathcal{M})$ into
$\text{End}(\mathcal{M})$ that can be expressed as

${\mathbf{D}}_{p|{\mathbf{S}}_{i}}\\{{\mathbf{T}}_{0,i}\\}=\sum_{r=1}^{\infty}f_{r}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m}){\mathbf{T}}_{0,i}h_{r}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m}),$
(37)

where $f_{r}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})$ and
$h_{r}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})$ are polynomial functions of
${\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m}$. As stated in [43] (p. 267), when
the operator is assumed to be Hilbert-Schmidt (finite Frobenius norm) there is
a unique linear operator
$\overline{{\mathbf{D}}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\left\\{\cdot\right\\}\in\text{End}({\mathcal{M}})^{\ast}\otimes\text{End}({\mathcal{M}})$
with
$\|\overline{{\mathbf{D}}}_{p|\mathbf{S}_{1}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\left\\{\cdot\right\\}\|=\|{\mathbf{D}}_{p|\mathbf{S}_{1}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\left\\{\cdot\right\\}\|$,
acting on a vectorized version of ${\mathbf{T}}_{0,i}$. In particular we have

$\overline{{\mathbf{D}}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})=\sum_{r=1}^{\infty}h_{r}^{\ast}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\otimes
f_{r}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m}).$ (38)

Now, we recall that since the representation $({\mathcal{M}},\rho)$ can be
expressed as a direct sum of the irreducible subrepresentations. Then, any
polynomial from ${\mathcal{A}}$ mapped into $\text{End}(\mathcal{M})$ by
$\rho$ can be expressed in terms of such direct sum by means of the spectral
theorem representation (Theorem 1), this is

$f_{r}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})=\sum_{k=1}^{q}f_{r}\left(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k}\right)\mathrsfso{P}_{k},$
(39)
$h_{r}^{\ast}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})=\sum_{k=1}^{q}h_{r}^{\ast}\left(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k}\right)\mathrsfso{P}_{k}.$
(40)

Then, substituting (39) and (40) in (38) we have

$\overline{{\mathbf{D}}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})=\sum_{r=1}^{\infty}\\\
\left(\sum_{k=1}^{q}h_{r}^{\ast}\left(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k}\right)\mathrsfso{P}_{k}\otimes\sum_{k=1}^{q}f_{r}\left(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k}\right)\mathrsfso{P}_{k}\right).$
(41)

Taking into account the properties of the tensor product it follows that:

$\overline{{\mathbf{D}}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})=\sum_{k,\ell=1}^{q}\\\
\left(\sum_{r=1}^{\infty}\left(h_{r}^{\ast}\left(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k}\right)\otimes
f_{r}\left(\boldsymbol{\Lambda}_{1,\ell},\ldots,\boldsymbol{\Lambda}_{m,\ell}\right)\right)\left(\mathrsfso{P}_{k}\otimes\mathrsfso{P}_{\ell}\right)\right).$
(42)

where we have that

$\sum_{r=1}^{\infty}\left(h_{r}^{\ast}\left(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k}\right)\otimes
f_{r}\left(\boldsymbol{\Lambda}_{1,\ell},\ldots,\boldsymbol{\Lambda}_{m,\ell}\right)\right)\\\
=\left\\{\begin{matrix}\overline{{\mathbf{D}}}_{p|\boldsymbol{\Lambda}_{i,k}}(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k})&\text{if}&k=\ell\\\
\boldsymbol{\Gamma}_{k,\ell}&\text{if}&k\neq\ell\end{matrix}\right.,$ (43)

where

$\overline{{\mathbf{D}}}_{p|\boldsymbol{\Lambda}_{i,k}}(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k})=\\\
\sum_{r=1}^{\infty}\left(h_{r}^{\ast}\left(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k}\right)\otimes
f_{r}\left(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k}\right)\right),$
(44)

is uniquely associated to
${\mathbf{D}}_{p|\boldsymbol{\Lambda}_{i,k}}(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k})$,
which is the Fréchet derivative of
$p\left(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k}\right)$
with respect o $\boldsymbol{\Lambda}_{i,k}$. Additionally, we point that the
terms in the sum of (42) are unique, and consequently can be expressed as a
direct sum. This is,

$\overline{{\mathbf{D}}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})=\bigoplus_{k,\ell=1}^{q}\sum_{r=1}^{\infty}\\\
\left(h_{r}^{\ast}\left(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k}\right)\otimes
f_{r}\left(\boldsymbol{\Lambda}_{1,\ell},\ldots,\boldsymbol{\Lambda}_{m,\ell}\right)\right)\left(\mathrsfso{P}_{k}\otimes\mathrsfso{P}_{\ell}\right).$
(45)

Then, if we calculate the norm of
$\overline{{\mathbf{D}}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})$
and we take into account the maximum property (Definition 12), it follows that

$\left\|\overline{{\mathbf{D}}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\right\|=\max_{(k,\ell)}\left\\{\left\|\overline{{\mathbf{D}}}_{p|\boldsymbol{\Lambda}_{i,k}}\right\|,\|\boldsymbol{\Gamma}_{k,\ell}\|\right\\}.$
(46)

First, we recall that
$\left\|\overline{{\mathbf{D}}}_{p|\boldsymbol{\Lambda}_{i,k}}\right\|=\left\|{\mathbf{D}}_{p|\boldsymbol{\Lambda}_{i,k}}\right\|$
and since $p$ is $L_{0}$-Lipschitz by Theorem 8 we have that
$\left\|\overline{{\mathbf{D}}}_{p|\boldsymbol{\Lambda}_{i,k}}\right\|\leq
L_{0}$. We focus now on estimating the norm of $\boldsymbol{\Gamma}_{k\ell}$.
First, we recall that

$\boldsymbol{\Gamma}_{k\ell}=\sum_{r=1}^{\infty}\left(h_{r}^{\ast}\left(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k}\right)\otimes
f_{r}\left(\boldsymbol{\Lambda}_{1,\ell},\ldots,\boldsymbol{\Lambda}_{m,\ell}\right)\right),$
(47)

with $k\neq\ell$. Now, we notice that it is possible to define a homomorphism
$\theta$ between the algebras of operators
$\text{End}({\mathcal{U}}_{k})^{\ast}\otimes\text{End}({\mathcal{U}}_{k})$ and
$\text{End}({\mathcal{U}}_{k})^{\ast}\otimes\text{End}({\mathcal{U}}_{\ell})$
with $\ell<k$, which is given by

$\theta\left(\sum_{r=1}^{\infty}\left(h_{r}^{\ast}\left(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k}\right)\otimes
f_{r}\left(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k}\right)\right)\right)\\\
=\sum_{r=1}^{\infty}\left(h_{r}^{\ast}\left(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k}\right)\otimes
f_{r}\left(\boldsymbol{\Lambda}_{1,\ell},\ldots,\boldsymbol{\Lambda}_{m,\ell}\right)\right).$
(48)

By means of Theorem 7 it follows that
$\left\|\boldsymbol{\Gamma}_{k\ell}\right\|\leq\left\|\overline{{\mathbf{D}}}_{p|\boldsymbol{\Lambda}_{i,k}}\right\|,$
and therefore leads to $\left\|\boldsymbol{\Gamma}_{k\ell}\right\|\leq L_{0}.$
Then, with the calculations above we have
$\left\|{\mathbf{D}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\right\|\leq
L_{0}.$ We recall that for bounded operators we have

$\displaystyle\left\|{\mathbf{D}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\left\\{{\mathbf{T}}_{0,i}\right\\}\right\|\leq\left\|{\mathbf{D}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\left\\{{\mathbf{T}}_{0,i}\right\\}\right\|_{F},$

and as stated in [43] (p. 267) we also have that

$\displaystyle\left\|{\mathbf{D}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\dots,{\mathbf{S}}_{m})\left\\{{\mathbf{T}}_{0,i}\right\\}\right\|_{F}\leq\left\|{\mathbf{D}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\right\|\left\|{\mathbf{T}}_{0,i}\right\|_{F}.$

Then, taking into account that
$\left\|{\mathbf{T}}_{0,i}\right\|_{F}\leq\delta\left\|{\mathbf{T}}_{0,i}\right\|,$
we reach

$\left\|{\mathbf{D}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\left\\{{\mathbf{T}}_{0,i}\right\\}\right\|\leq
L_{0}\delta\sup_{{\mathbf{S}}_{i}}\|{\mathbf{T}}({\mathbf{S}}_{i})\|.$ (49)

Now, we turn our attention to the term
$\left\|{\mathbf{D}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\left\\{{\mathbf{T}}_{1,i}{\mathbf{S}}_{i}\right\\}\right\|$
in (36). We start using the notation

${\mathbf{D}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\left\\{{\mathbf{T}}_{1,i}{\mathbf{S}}_{i}\right\\}=\widetilde{{\mathbf{D}}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\left\\{{\mathbf{T}}_{1,i}\right\\}.$
(50)

Now, we point out that

$\widetilde{{\mathbf{D}}}_{p|{\mathbf{S}}_{i}}\\{{\mathbf{T}}_{1,i}\\}=\sum_{r=1}^{\infty}f_{r}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m}){\mathbf{T}}_{1,i}h_{r}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m}),$
(51)

where $f_{r}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})$ and
$h_{r}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})$ are polynomial functions of
${\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m}$. When the operator is assumed to be
Hilbert-Schmidt there is a unique linear operator given by
$\overline{\widetilde{{\mathbf{D}}}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\left\\{\cdot\right\\}\in\text{End}({\mathcal{M}})^{\ast}\otimes\text{End}({\mathcal{M}})$
with
$\|\overline{\widetilde{{\mathbf{D}}}}_{p|\mathbf{S}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\left\\{\cdot\right\\}\|=\|\widetilde{{\mathbf{D}}}_{p|\mathbf{S}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\left\\{\cdot\right\\}\|$,
acting on a vectorized version of ${\mathbf{T}}_{1,i}$ [43] (p. 267). In
particular we have

$\overline{\widetilde{{\mathbf{D}}}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})=\sum_{r=1}^{\infty}h_{r}^{\ast}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\otimes
f_{r}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m}).$ (52)

Since $({\mathcal{M}},\rho)$ can be expressed as a direct sum of the
irreducible subrepresentations. Then, any polynomial from ${\mathcal{A}}$
mapped into $\text{End}(\mathcal{M})$ by $\rho$ can be expressed in terms of
such direct sum by means of the spectral theorem representation (Theorem 1),
this is

$f_{r}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})=\sum_{k=1}^{q}f_{r}\left(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k}\right)\mathrsfso{P}_{k},$
(53)
$h_{r}^{\ast}({\mathbf{S}}_{1},\dots,{\mathbf{S}}_{m})=\sum_{k=1}^{q}h_{r}^{\ast}\left(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k}\right)\mathrsfso{P}_{k}.$
(54)

Then, substituting (53) and (54) in (52) we have

$\overline{\widetilde{{\mathbf{D}}}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})=\sum_{r=1}^{\infty}\\\
\left(\sum_{k=1}^{q}h_{r}^{\ast}\left(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k}\right)\mathrsfso{P}_{k}\otimes\sum_{k=1}^{q}f_{r}\left(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k}\right)\mathrsfso{P}_{k}\right).$

Taking into account the properties of the tensor product it follows that:

$\overline{\widetilde{{\mathbf{D}}}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})=\sum_{k,\ell=1}^{q}\\\
\left(\sum_{r=1}^{\infty}\left(h_{r}^{\ast}\left(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k}\right)\otimes
f_{r}\left(\boldsymbol{\Lambda}_{1,\ell},\ldots,\boldsymbol{\Lambda}_{m,\ell}\right)\right)\left(\mathrsfso{P}_{k}\otimes\mathrsfso{P}_{\ell}\right)\right),$
(55)

where we have that

$\sum_{r=1}^{\infty}\left(h_{r}^{\ast}\left(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k}\right)\otimes
f_{r}\left(\boldsymbol{\Lambda}_{1,\ell},\ldots,\boldsymbol{\Lambda}_{m,\ell}\right)\right)\\\
=\left\\{\begin{matrix}\overline{\widetilde{{\mathbf{D}}}}_{p|\boldsymbol{\Lambda}_{i,k}}(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k})&\text{if}&k=\ell\\\
\boldsymbol{\Gamma}_{k,\ell}&\text{if}&k\neq\ell\end{matrix}\right.,$ (56)

with

$\overline{\widetilde{{\mathbf{D}}}}_{p|\boldsymbol{\Lambda}_{i,k}}(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k})\\\
=\sum_{r=1}^{\infty}\left(h_{r}^{\ast}\left(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k}\right)\otimes
f_{r}\left(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k}\right)\right),$
(57)

is uniquely associated to
${\mathbf{D}}_{p|\boldsymbol{\Lambda}_{i,k}}(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k})\\{(\cdot)\boldsymbol{\Lambda}_{i,k}\\}$.
Like in the scenarios discussed above, we point that the terms in the sum of
(55) are unique, and consequently can be expressed as a direct sum as follows

$\overline{\widetilde{{\mathbf{D}}}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})=\bigoplus_{k,\ell=1}^{q}\sum_{r=1}^{\infty}\\\
\left(h_{r}^{\ast}\left(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k}\right)\otimes
f_{r}\left(\boldsymbol{\Lambda}_{1,\ell},\ldots,\boldsymbol{\Lambda}_{m,\ell}\right)\right)\left(\mathrsfso{P}_{k}\otimes\mathrsfso{P}_{\ell}\right).$
(58)

If we calculate the norm of
$\overline{\widetilde{{\mathbf{D}}}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})$
and we use the maximum property (Definition 12), it follows that

$\left\|\overline{\widetilde{{\mathbf{D}}}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\right\|=\max_{(k,\ell)}\left\\{\left\|\overline{\widetilde{{\mathbf{D}}}}_{p|\boldsymbol{\Lambda}_{i,k}}\right\|,\|\boldsymbol{\Gamma}_{k,\ell}\|\right\\}.$
(59)

Since
$\left\|\overline{\widetilde{{\mathbf{D}}}}_{p|\boldsymbol{\Lambda}_{i,k}}\right\|=\left\|\widetilde{{\mathbf{D}}}_{p|\boldsymbol{\Lambda}_{i,k}}\right\|$,
and $p$ is $L_{1}$-integral Lipschitz, we have that
$\left\|\overline{\widetilde{{\mathbf{D}}}}_{p|\boldsymbol{\Lambda}_{i,k}}\right\|\leq
L_{1}$. Now we analyze the norm of $\boldsymbol{\Gamma}_{k\ell}$. First, we
recall that

$\boldsymbol{\Gamma}_{k\ell}=\sum_{r=1}^{\infty}\left(h_{r}^{\ast}\left(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k}\right)\otimes
f_{r}\left(\boldsymbol{\Lambda}_{1,\ell},\ldots,\boldsymbol{\Lambda}_{m,\ell}\right)\right),$
(60)

with $k\neq\ell$. We point out that it is possible to define a homomorphism
$\theta$ between the algebras of operators
$\text{End}({\mathcal{U}}_{k})^{\ast}\otimes\text{End}({\mathcal{U}}_{k})$ and
$\text{End}({\mathcal{U}}_{k})^{\ast}\otimes\text{End}({\mathcal{U}}_{\ell})$
with $k<\ell$, which is given by

$\theta\left(\sum_{r=1}^{\infty}\left(h_{r}^{\ast}\left(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k}\right)\otimes
f_{r}\left(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k}\right)\right)\right)\\\
=\sum_{r=1}^{\infty}\left(h_{r}^{\ast}\left(\boldsymbol{\Lambda}_{1,k},\ldots,\boldsymbol{\Lambda}_{m,k}\right)\otimes
f_{r}\left(\boldsymbol{\Lambda}_{1,\ell},\ldots,\boldsymbol{\Lambda}_{m,\ell}\right)\right).$
(61)

By Theorem 7, it follows that
$\left\|\boldsymbol{\Gamma}_{k\ell}\right\|\leq\left\|\overline{\widetilde{{\mathbf{D}}}}_{p|\boldsymbol{\Lambda}_{i,k}}\right\|,$
which leads to $\left\|\boldsymbol{\Gamma}_{k\ell}\right\|\leq L_{1},$ and
therefore
$\left\|\widetilde{{\mathbf{D}}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\right\|\leq
L_{1}.$ If
$\widetilde{{\mathbf{D}}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})$
is a bounded operator, we have

$\displaystyle\left\|\widetilde{{\mathbf{D}}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\left\\{{\mathbf{T}}_{1,i}\right\\}\right\|\leq\left\|\widetilde{{\mathbf{D}}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\left\\{{\mathbf{T}}_{1,i}\right\\}\right\|_{F},$

and as stated in [43] (p. 267) we also have that

$\displaystyle\left\|\widetilde{{\mathbf{D}}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\left\\{{\mathbf{T}}_{1,i}\right\\}\right\|_{F}\leq\left\|\widetilde{{\mathbf{D}}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\right\|\left\|{\mathbf{T}}_{1,i}\right\|_{F}.$

Taking into account that
$\left\|{\mathbf{T}}_{1,i}\right\|_{F}\leq\delta\left\|{\mathbf{T}}_{1,i}\right\|,$
we reach

$\displaystyle\left\|{\mathbf{D}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\left\\{{\mathbf{T}}_{1,i}{\mathbf{S}}_{i}\right\\}\right\|\leq
L_{1}\delta\left\|{\mathbf{T}}_{1,i}\right\|,$ (62)

and since ${\mathbf{D}}_{{\mathbf{T}}}({\mathbf{S}}_{i})={\mathbf{T}}_{1,i}$,
this leads to

$\displaystyle\left\|{\mathbf{D}}_{p|{\mathbf{S}}_{i}}({\mathbf{S}}_{1},\ldots,{\mathbf{S}}_{m})\left\\{{\mathbf{T}}_{1,i}{\mathbf{S}}_{i}\right\\}\right\|\leq
L_{1}\delta\sup_{{\mathbf{S}}_{i}}\left\|{\mathbf{D}}_{{\mathbf{T}}}({\mathbf{S}}_{i})\right\|.$
(63)

∎

## References

  * [1] A. Parada-Mayorga and A. Ribeiro.  Algebraic neural networks: Stability to deformations.  IEEE Transactions on Signal Processing, 69:3351–3366, 2021.
  * [2] Fernando Gama, Joan Bruna, and Alejandro Ribeiro.  Stability properties of graph neural networks.  IEEE Transactions on Signal Processing, 68:5680–5695, 2020.
  * [3] Fernando Gama, Joan Bruna Estrach, and Alejandro Ribeiro.  Diffusion scattering transforms on graphs.  1 2019.  7th International Conference on Learning Representations, ICLR 2019 ; Conference date: 06-05-2019 Through 09-05-2019.
  * [4] S. Mallat.  Group invariant scattering.  Communications on Pure and Applied Mathematics, 65(10):1331–1398, 2012.
  * [5] M. Püschel and J. M. F. Moura.  Algebraic signal processing theory.  ArXiv, abs/cs/0612077, 2006.
  * [6] A. Parada-Mayorga and A. Ribeiro.  Stability of algebraic neural networks to small perturbations.  In ICASSP 2021 - 2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 5205–5209, 2021.
  * [7] B. Knyazev, X. Lin, M. R. Amer, and G. W. Taylor.  Spectral multigraph networks for discovering and fusing relationships in molecules.  ArXiv, abs/1811.09595, 2018.
  * [8] M. Weiler and G. Cesa.  General e(2)-equivariant steerable cnns.  ArXiv, abs/1911.08251, 2019.
  * [9] M. Weiler, M. Geiger, M. Welling, W. Boomsma, and T. Cohen.  3d steerable cnns: Learning rotationally equivariant features in volumetric data.  In NeurIPS, 2018.
  * [10] D. E. Worrall, S. J. Garbin, D. Turmukhambetov, and G. J. Brostow.  Harmonic networks: Deep translation and rotation equivariance.  2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 7168–7177, 2017.
  * [11] M. Finzi, S. Stanton, P. Izmailov, and A. Gordon Wilson.  Generalizing convolutional neural networks for equivariance to lie groups on arbitrary continuous data.  In ICML, 2020.
  * [12] M. Finzi, M. Welling, and A. Gordon Wilson.  A practical method for constructing equivariant multilayer perceptrons for arbitrary matrix groups.  ArXiv, abs/2104.09459, 2021.
  * [13] J. Hoffmann, S. Schmitt, S. Osindero, K. Simonyan, and E. Elsen.  Algebranets.  ArXiv, abs/2006.07360, 2020.
  * [14] C. J. Gaudet and A. S. Maida.  Deep quaternion networks.  In 2018 International Joint Conference on Neural Networks (IJCNN), pages 1–8. IEEE, 2018.
  * [15] Harshat Kumar, Alejandro Parada-Mayorga, and Alejandro Ribeiro.  Algebraic convolutional filters on lie group algebras.  2022\.
  * [16] Landon Butler, Alejandro Parada-Mayorga, and Alejandro Ribeiro.  Convolutional learning on multigraphs.  IEEE Transactions on Signal Processing, 71:933–946, 2023.
  * [17] Michael M. Bronstein, Joan Bruna, Taco Cohen, and Petar Velivckovi’c.  Geometric deep learning: Grids, groups, graphs, geodesics, and gauges.  ArXiv, abs/2104.13478, 2021.
  * [18] Maurice Weiler, Patrick Forré, Erik Verlinde, and Max Welling.  Coordinate independent convolutional networks–isometry and gauge equivariant convolutions on riemannian manifolds.  arXiv preprint arXiv:2106.06020, 2021.
  * [19] Landon Butler, Alejandro Parada-Mayorga, and Alejandro Ribeiro.  Learning with multigraph convolutional filters.  In ICASSP 2023 - 2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 1–5, 2023.
  * [20] Dai Quoc Nguyen, Tu Dinh Nguyen, and Dinh Q. Phung.  Quaternion graph neural networks.  In ACML, 2020.
  * [21] Titouan Parcollet, Mirco Ravanelli, Mohamed Morchid, Georges Linarès, Chiheb Trabelsi, Renato De Mori, and Yoshua Bengio.  Quaternion recurrent neural networks.  ArXiv, abs/1806.04418, 2018.
  * [22] Xinchi Qiu, Titouan Parcollet, Mirco Ravanelli, Nicholas D. Lane, and Mohamed Morchid.  Quaternion neural networks for multi-channel distant speech recognition.  In Interspeech, 2020.
  * [23] Yili Xia, Min Xiang, Zhe Li, and Danilo P Mandic.  Echo state networks for multidimensional data: Exploiting noncircularity and widely linear models.  In Adaptive Learning Methods for Nonlinear System Modeling, pages 267–288. Elsevier, 2018.
  * [24] Titouan Parcollet, Mohamed Morchid, and Georges Linarès.  A survey of quaternion neural networks.  Artificial Intelligence Review, 53(4):2957–2982, 2020.
  * [25] Eleonora Grassucci, Edoardo Cicero, and Danilo Comminiello.  Quaternion generative adversarial networks.  In Generative Adversarial Learning: Architectures and Applications, pages 57–86. Springer, 2022.
  * [26] Eleonora Grassucci, Aston Zhang, and Danilo Comminiello.  Phnns: Lightweight neural networks via parameterized hypercomplex convolutions.  2021\.
  * [27] S. Buchholz and G. Sommer.  A hyperbolic multilayer perceptron.  In Proceedings of the IEEE-INNS-ENNS International Joint Conference on Neural Networks. IJCNN 2000. Neural Computing: New Challenges and Perspectives for the New Millennium, volume 2, pages 129–133 vol.2, 2000\.
  * [28] Tohru Nitta and Sven Buchholz.  On the decision boundaries of hyperbolic neurons.  In 2008 IEEE International Joint Conference on Neural Networks (IEEE World Congress on Computational Intelligence), pages 2974–2980. IEEE, 2008\.
  * [29] Călin-Adrian Popa.  Octonion-valued neural networks.  In International Conference on Artificial Neural Networks, pages 435–443. Springer, 2016.
  * [30] Eduardo Bayro-Corrochano, Refugio Vallejo, and Nancy Arana-Daniel.  Geometric preprocessing, geometric feedforward neural networks and Clifford support vector machines for visual learning.  Neurocomputing, 67:54–105, 2005.
  * [31] Sven Buchholz and Gerald Sommer.  On Clifford neurons and Clifford multi-layer perceptrons.  Neural Networks, 21(7):925–935, 2008.
  * [32] Sven Buchholz and Gerald Sommer.  Clifford algebra multilayer perceptrons.  In Geometric Computing with Clifford Algebras, pages 315–334. Springer, 2001.
  * [33] A. Parada-Mayorga, H. Riess, A. Ribeiro, and R. Ghrist.  Quiver signal processing (qsp).  ArXiv, abs/2010.11525, 2020.
  * [34] G.B. Folland.  A Course in Abstract Harmonic Analysis.  Textbooks in Mathematics. CRC Press, 2016.
  * [35] M. Lorenz.  A Tour of Representation Theory.  Graduate studies in mathematics. American Mathematical Society, 2018.
  * [36] A. Terras.  Fourier Analysis on Finite Groups and Applications.  London Mathematical Society Student Texts. Cambridge University Press, 1999.
  * [37] P. Diaconis, P.W. Diaconis, Institute of Mathematical Statistics, and JSTOR (Organization).  Group Representations in Probability and Statistics.  IMS Lecture Notes. Institute of Mathematical Statistics, 1988.
  * [38] P.I. Etingof, O. Golberg, S. Hensel, T. Liu, A. Schwendner, D.V. , and E. Yudovina.  Introduction to Representation Theory.  Student mathematical library. American Mathematical Society, 2011.
  * [39] F. Gama, A. G. Marques, G. Leus, and A. Ribeiro.  Convolutional neural network architectures for signals supported on graphs.  IEEE Transactions on Signal Processing, 67(4):1034–1049, 2019.
  * [40] Y-Lan Boureau, Jean Ponce, and Yann LeCun.  A theoretical analysis of feature pooling in visual recognition.  In Proceedings of the 27th International Conference on International Conference on Machine Learning, ICML’10, page 111–118, Madison, WI, USA, 2010. Omnipress.
  * [41] Y-Lan Boureau, Nicolas Le Roux, Francis Bach, Jean Ponce, and Yann LeCun.  Ask the locals: Multi-way local pooling for image recognition.  In 2011 International Conference on Computer Vision, pages 2651–2658, 2011.
  * [42] Alejandro Parada-Mayorga, Zhiyang Wang, and Alejandro Ribeiro.  Graphon pooling for reducing dimensionality of signals and convolutional operators on graphs.  ArXiv, abs/2212.08171, 2022.
  * [43] J.B. Conway.  A Course in Functional Analysis.  Graduate Texts in Mathematics. Springer New York, 1994.
  * [44] M.W. Hirsch.  Differential Topology.  Graduate Texts in Mathematics. Springer New York, 1997.
  * [45] J.A. Leslie.  On a differential structure for the group of diffeomorphisms.  Topology, 6(2):263–271, 1967.
  * [46] A. Banyaga.  The Structure of Classical Diffeomorphism Groups.  Mathematics and Its Applications. Springer US, 2013.
  * [47] F.M. Harper and J. A. Konstan.  The movielens datasets: History and context.  Acm transactions on interactive intelligent systems (tiis), 5(4):1–19, 2015.
  * [48] L. Deng.  The mnist database of handwritten digit images for machine learning research.  IEEE Signal Processing Magazine, 29(6):141–142, 2012.
  * [49] W. Huang, A. G. Marques, and A. R. Ribeiro.  Rating prediction via graph signal processing.  IEEE Transactions on Signal Processing, 66(19):5066–5081, 2018\.
  * [50] Titouan Parcollet, Mohamed Morchid, and Georges Linarès.  Quaternion convolutional neural networks for heterogeneous image processing.  In ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 8514–8518. IEEE, 2019.
  * [51] M.S. Berger.  Nonlinearity and Functional Analysis: Lectures on Nonlinear Problems in Mathematical Analysis.  Pure and Applied Mathematics. Elsevier Science, 1977.
  * [52] P. Lancaster and H. K. Farahat.  Norms on direct sums and tensor products.  Mathematics of Computation, 26(118):401–414, 1972.
  * [53] R.A. Horn and C.R. Johnson.  Matrix Analysis.  Matrix Analysis. Cambridge University Press, 2012.
  * [54] K. Lux and H. Pahlings.  Representations of Groups: A Computational Approach.  Cambridge Studies in Advanced Mathematics. Cambridge University Press, 2010.
  * [55] P. Dräxler, G. Michler, and C.M. Ringel.  Computational Methods for Representations of Groups and Algebras: Euroconference in Essen (Germany), April 1–5, 1977.  Progress in Mathematics. Birkhäuser Basel, 2012.
  * [56] M. Barot.  Introduction to the Representation Theory of Algebras.  Springer International Publishing, 2014.

##  Appendix E Preliminaries Irreducible Subrepresentations

In the following definition we define homomorphisms and isomorphisms of
representations, which will allow us to state comparisons between
representations.

######  Definition 10.

Let $(\mathcal{M}_{1},\rho_{1})$ and $(\mathcal{M}_{2},\rho_{2})$ be two
representations of an algebra $\mathcal{A}$. A homomorphism or interwining
operator $\phi:\mathcal{M}_{1}\rightarrow\mathcal{M}_{2}$ is a linear operator
which commutes with the action of $\mathcal{A}$ under $\rho_{1}$ and
$\rho_{2}$, i.e.

$\phi(\rho_{1}(a)v)=\rho_{2}(a)\phi(v)\quad v\in{\mathcal{M}}_{1}.$ (64)

A homomorphism $\phi$ is said to be an isomorphism of representations if it is
an isomorphism of vectors spaces.

It is worth pointing out that even when two representations have the same
dimension, it might still happens that they are not homomorphic since the
condition in (64) might not be satisfied. This is the typical scenario we have
when considering representations of the polynomial algebra ${\mathbb{F}}[t]$,
where we have irreducible representations that are one dimensional but in
general non isomorphic.

Now, we highlight that irreducible representations are indecomposable – the
opposite is in general false – and therefore, irreducible representations can
be used to build unique decompositions. We focus our attention on the class of
representations that can be decomposed as a direct sum of irreducible
representations.

##  Appendix F Expanded example multigraphs

######  Remark 7 (Non commutative polynomial notation).

Before we continue our presentation of non commutative signal models we make
some remarks about the notation for non commutative polynomials. If a non
commutative algebra ${\mathcal{A}}$ is generated by
$\\{g_{1},\ldots,g_{m}\\}$, we define the set
$\Pi_{j_{1},\ldots,j_{m}}(g_{1},\ldots,g_{m})$ of normalized monomials built
from products of the generators and where the term $g_{i}$ appears $j_{i}$
times. For instance, if the algebra has two generators $g_{1}$ and $g_{2}$ we
have

$\Pi_{1,2}(g_{1},g_{2})=\left\\{g_{1}g_{2}^{2},g_{2}g_{1}g_{2},g_{2}^{2}g_{1}\right\\}.$
(65)

Given an arbitrary ordering of the elements in
$\Pi_{j_{1},\ldots,j_{m}}(g_{1},\ldots,g_{m})$ we will represent its $r$-th
element by $\boldsymbol{\pi}_{j_{1},\ldots,j_{m}}^{r}(g_{1},\ldots,g_{m})$.
Then, any polynomial in the non commutative algebra ${\mathcal{A}}$ with $m$
generators can be written as

$p(g_{1},\ldots,g_{m})=\sum_{j_{1},\ldots,j_{m}=0}\sum_{r}h_{j_{1},\ldots,j_{m};r}\boldsymbol{\pi}_{j_{1},\ldots,j_{m}}^{r}(g_{1},\ldots,g_{m}),$
(66)

where $h_{j_{1},\ldots,j_{m};r}$ indicates the coefficients of the polynomial.

######  Example 3 (Multigraph signal processing).

A multigraph consists of a common set of vertices and multiple separate sets
of edges. Consider the specific case
$G=\left({\mathcal{V}},\\{{\mathcal{E}}_{1},{\mathcal{E}}_{2}\\}\right)$
consisting of a set of $N$ vertices ${\mathcal{V}}$ and two separate edge sets
${\mathcal{E}}_{1}$ and ${\mathcal{E}}_{2}$ – see Fig. 2. Associated with each
edge set we consider matrix representations ${\mathbf{S}}_{1}$ and
${\mathbf{S}}_{2}$. This is an algebraic signal processing model in which the
vector space of signals ${\mathcal{M}}={\mathbb{C}}^{N}$ contains vectors in
${\mathbb{C}}^{N}$ with entries associated with each node of the graph and the
algebra ${\mathcal{A}}={\mathbb{C}}[t_{1},t_{2}]$ is the set of non
commutative polynomials of two variables. Any polynomial in
${\mathbb{C}}[t_{1},t_{2}]$ can be written as

$p(t_{1},t_{2})=\sum_{j_{1},j_{2}=0}^{d}\sum_{r}h_{j_{1},j_{2};r}\boldsymbol{\pi}_{j_{1},j_{2}}^{r}(t_{1},t_{2}).$
(67)

The algebra ${\mathbb{C}}[t_{1},t_{2}]$ is generated by the monomials $t_{1}$
and $t_{2}$. The homomorphism $\rho$ is defined by mapping the generators
$t_{1}$ and $t_{2}$ to $\rho(t_{1})={\mathbf{S}}_{1}$ and
$\rho(t_{2})={\mathbf{S}}_{2}$, where the shift operators ${\mathbf{S}}_{i}$
are the matrix representations of the corresponding set of edges
${\mathcal{E}}_{i}$. Then, the filtering in (2) takes the form

$\rho\left(\sum_{j_{1},j_{2}=0}^{d}\sum_{r}h_{j_{1},j_{2};r}\boldsymbol{\pi}_{j_{1},j_{2}}^{r}(t_{1},t_{2})\right){\mathbf{x}}\\\
=\left(\sum_{j_{1},j_{2}=0}^{d}\sum_{r}h_{j_{1},j_{2};r}\boldsymbol{\pi}_{j_{1},j_{2}}^{r}({\mathbf{S}}_{1},{\mathbf{S}}_{2})\right){\mathbf{x}}.$
(68)

##  Appendix G Proof of Theorem 1

Since $({\mathcal{M}},\rho)$ can be written as a direct sum of irreducible
subrepresentations $({\mathcal{U}},\phi_{i})$, we have that for any
$a\in{\mathcal{A}}$ the action of $\rho(a)$ on ${\mathbf{x}}\in{\mathcal{M}}$
can be written as $\rho(a)=\sum_{i=1}^{q}\phi_{i}(a)\mathrsfso{P}_{i}.$ Then,

$\rho(p(g_{1},\ldots,g_{m}))=p\left(\rho(g_{1}),\ldots,\rho(g_{m})\right)=\\\
p\left(\sum_{i=1}^{q}\phi_{i}(g_{1})\mathrsfso{P}_{i},\ldots,\sum_{i=1}^{q}\phi_{i}(g_{m})\mathrsfso{P}_{i}\right).$
(69)

Since

$\mathrsfso{P}_{j}\phi_{i}(a)\mathrsfso{P}_{i}=\left\\{\begin{matrix}\phi_{i}(a)\mathrsfso{P}_{i}&&i=j\\\
0&&i\neq j\end{matrix}\right.,$ (70)

it follows that

$p\left(\sum_{i=1}^{q}\phi_{i}(g_{1})\mathrsfso{P}_{i},\ldots,\sum_{i=1}^{q}\phi_{i}(g_{m})\mathrsfso{P}_{i}\right)\\\
=\sum_{i=1}^{q}p\left(\boldsymbol{\Lambda}_{1,i},\ldots,\boldsymbol{\Lambda}_{m,i}\right)\mathrsfso{P}_{i},$
(71)

where $\boldsymbol{\Lambda}_{k,i}=\phi_{i}(g_{k})$.

##  Appendix H Norms in Direct Sums of Spaces

In this section we discuss the properties of operator and matrix norms when
working on direct sums of spaces. This plays a central role on how we
calculate the norms of those operators associated to completely reducible
representations. We start introducing the notion of absolute norm.

######  Definition 11 ([52]).

Let ${\mathcal{X}}_{1},\ldots{\mathcal{X}}_{q}$ be normed linear spaces and
let ${\mathcal{X}}$ be given by
${\mathcal{X}}={\mathcal{X}}_{1}\oplus{\mathcal{X}}_{2}\oplus\ldots\oplus{\mathcal{X}}_{q}.$
If $x,y\in{\mathcal{X}}$ and $x=(x_{1},x_{2},\ldots,x_{q})$,
$y=(y_{1},y_{2},\ldots,y_{q})$, we say that a norm on ${\mathcal{X}}$ is
absolute if $\|x_{i}\|=\|y_{i}\|$, $j=1,\ldots,n$, implies $\|x\|=\|y\|$.

Now, we introduce a property of norms of direct sums of spaces that will
become relevant in our proofs.

######  Definition 12 ([52]).

A norm $\|\cdot\|$ on ${\mathcal{X}}=\bigoplus_{i=1}^{q}{\mathcal{X}}_{i}$ has
the maximum property if, for any linear operator
${\mathbf{T}}=\bigoplus_{i=1}^{q}{\mathbf{T}}_{i}$, we have that
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|{\mathbf{T}}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}=\max_{i}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|{\mathbf{T}}_{i}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$,
where
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|{\mathbf{T}}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$
indicates the operator norm induced by $\|\cdot\|$, and
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|{\mathbf{T}}_{i}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$
is the operator norm induced by the norm in ${\mathcal{X}}_{i}$.

Notice that the maximum property is defined on $\|\cdot\|$ but measured or
captured with
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$.
The following results provides a formal connections between absolute norms and
the maximum property.

######  Theorem 6 ([52]).

A norm on ${\mathcal{X}}$ is absolute iff it has the maximum property.

As pointed out in [53] (p. 342) the Frobenius norm $\|\cdot\|_{F}$ is an
absolute norm. This plays a central role when considering operators on
completely reducible representations. To see this, we recall that if the
representation $({\mathcal{M}},\rho)$ – of an algebra ${\mathcal{A}}$ – is
completely reducible, it follows that
$({\mathcal{M}},\rho)\cong\bigoplus_{i=1}^{q}({\mathcal{U}}_{i},\phi_{i})$.
Then, any operator $\rho(a)$ can be expressed as
$\rho(a)=\bigoplus_{i=1}^{q}\phi_{i}(a)$ acting on
$\bigoplus_{i=1}^{q}{\mathcal{U}}_{i}$. Therefore, when computing the
Frobenius norm of $\rho(a)$, we will have that
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\rho(a)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}=\max_{i}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\phi_{i}(a)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$.
Notice that this is consistent with the scenario of commutative algebras where
the subpresentations are 1-dimensional and where $\phi_{i}$ is a scalar
polynomial evaluated at the $i$th eigenvalue.

##  Appendix I Relationship between $\|a\|$ and $\|\rho(a)\|$

######  Theorem 7.

Let ${\mathcal{A}}$ and ${\mathcal{B}}$ be two polynomial algebras of
operators endowed with norms satisfying that $\|aa^{\ast}\|=\|a\|^{2}$ and
where ∗ indicates the adjoint operation. Let $\theta$ be a homomorphism
between ${\mathcal{A}}$ and ${\mathcal{B}}$ with
$\theta(a^{\ast})=\theta(a)^{\ast}$. Then,
$\left\|\theta(a)\right\|\leq\left\|a\right\|.$

###### Proof.

First, notice that the unit element in ${\mathcal{A}}$ given by
$1_{{\mathcal{A}}}$ is mapped into the unit element of ${\mathcal{B}}$ given
by $1_{{\mathcal{B}}}$. Then, if $\mu\in\mathbb{C}$ is not in the spectrum of
$a$, the term $\mu 1_{{\mathcal{A}}}-a$ has an inverse $\left(\mu
1_{{\mathcal{A}}}-a\right)^{-1}\in{\mathcal{A}}$. Taking into account that
$\theta\left(1_{{\mathcal{A}}}\right)=1_{{\mathcal{B}}}$ we have that
$\theta\left(\mu 1_{{\mathcal{A}}}-a\right)\theta\left(\left(\mu
1_{{\mathcal{A}}}-a\right)^{-1}\right)=\theta\left(1_{{\mathcal{A}}}\right)=1_{{\mathcal{B}}}.$
Then $\theta\left(\mu 1_{{\mathcal{A}}}-a\right)=\mu
1_{{\mathcal{B}}}-\theta(a)$ has an inverse in ${\mathcal{A}}$ which implies
that $\mu$ is not in the spectrum of $\theta(a)$. Consequently,
$\text{spec}\left(\theta(a)\right)\subset\text{spec}\left(a\right),$ where
$\text{spec}\left(\theta(a)\right)$ and $\text{spec}\left(a\right)$ are the
spectrum sets of $\theta(a)$ and $a$ respectively. Now, we take into account
that $\left\|a\right\|^{2}=\left\|aa^{\ast}\right\|=r\left(aa^{\ast}\right),$
$\left\|\theta(a)\right\|^{2}=\left\|\theta(a)\theta(a)^{\ast}\right\|=r\left(\theta(a)\theta(a)^{\ast}\right),$
where $r(\cdot)$ represents the spectral radius. Then, we have that
$r\left(\theta(a)\theta(a)^{\ast}\right)\leq r\left(aa^{\ast}\right)$ and
therefore $\|\theta(a)\|\leq\|a\|$. ∎

##  Appendix J Fréchet Derivative of Lipschitz Functions

######  Theorem 8.

Let ${\mathbf{F}}:{\mathcal{V}}\to{\mathcal{U}}$ a Lipschitz map between two
Banach spaces ${\mathcal{V}}$ and ${\mathcal{U}}$ with Lipschitz constant
$L_{0}$. Let ${\mathbf{D}}_{{\mathbf{F}}|{\mathbf{S}}}({\mathbf{S}})$ be the
Fréchet derivative of ${\mathbf{F}}$ at ${\mathbf{S}}\in{\mathcal{V}}$. Then,
if ${\mathbf{D}}_{{\mathbf{F}}|{\mathbf{S}}}({\mathbf{S}})$ exists we have
that
$\left\|{\mathbf{D}}_{{\mathbf{F}}|{\mathbf{S}}}({\mathbf{S}})\right\|\leq
L_{0}.$

###### Proof.

We start taking into account the definition of the Fréchet derivative of
${\mathbf{F}}({\mathbf{S}})$,

${\mathbf{D}}_{{\mathbf{F}}|{\mathbf{S}}}\\{\boldsymbol{\xi}\\}={\mathbf{F}}({\mathbf{S}}+\boldsymbol{\xi})-{\mathbf{F}}({\mathbf{S}})-o(\boldsymbol{\xi}).$
(72)

Taking the norm on both sides of (72) and applying the triangular inequality
this leads to

$\left\|{\mathbf{D}}_{{\mathbf{F}}|{\mathbf{S}}}\\{\boldsymbol{\xi}\\}\right\|\leq\left\|{\mathbf{F}}({\mathbf{S}}+\boldsymbol{\xi})-{\mathbf{F}}({\mathbf{S}})\right\|+\left\|o(\boldsymbol{\xi})\right\|.$
(73)

Since ${\mathbf{F}}({\mathbf{S}})$ is $L_{0}$-Lipschitz, it follows that

$\left\|{\mathbf{D}}_{{\mathbf{F}}|{\mathbf{S}}}\\{\boldsymbol{\xi}\\}\right\|\leq\left\|\boldsymbol{\xi}\right\|\left(L_{0}+\frac{\left\|o(\boldsymbol{\xi})\right\|}{\|\boldsymbol{\xi}\|}\right).$
(74)

Now, we recall one of the definitions of the operator norm [43] (p. 12) we
have $\|{\mathbf{D}}_{{\mathbf{F}}|{\mathbf{S}}}\|\leq
L_{0}+\frac{\left\|o(\boldsymbol{\xi})\right\|}{\|\boldsymbol{\xi}\|}.$
Finally, taking into account that
$\left\|o(\boldsymbol{\xi})\right\|/\|\boldsymbol{\xi}\|\rightarrow 0$ as
$\|\boldsymbol{\xi}\|\rightarrow 0$, it follows that
$\|{\mathbf{D}}_{{\mathbf{F}}|{\mathbf{S}}}\|\leq L_{0}.$ ∎

###  J-A Proof of Theorems 4 and 5

####  J-A1 Proof of Theorem 4

###### Proof.

Considering (10), applying the operator norm, and taking into account that the
maps $\sigma_{\ell}$ are Lipschitz with constant $C_{\ell}$, we have

$\left\|\sigma_{\ell}\left(p({\mathbf{S}}_{1,\ell},\ldots,{\mathbf{S}}_{m,\ell})\mathbf{x}_{\ell-1}\right)-\sigma_{\ell}\left(p(\tilde{\mathbf{S}}_{1,\ell},\ldots,\tilde{{\mathbf{S}}}_{m,\ell})\mathbf{x}_{\ell-1}\right)\right\|\\\
\leq C_{\ell}\boldsymbol{\Delta}_{\ell}\|\mathbf{x}_{\ell-1}\|,$ (75)

where
$\boldsymbol{\Delta}_{\ell}=\|p(\mathbf{S}_{1,\ell},\ldots,{\mathbf{S}}_{m,\ell})-p(\tilde{\mathbf{S}}_{1,\ell},\ldots,\tilde{{\mathbf{S}}}_{m,\ell})\|$,
and whose value is determined by theorems 2 and 3. ∎

####  J-A2 Proof of Theorem 5

We start writing the norm of the difference between the $f$-th feature and its
perturbed version in the $\ell$-th layer.

###### Proof.

$\left\|\mathbf{x}_{\ell}^{f}-\tilde{\mathbf{x}}_{\ell}^{f}\right\|\leq\left\|\sigma_{\ell-1}\sum_{g_{\ell-1}}\rho_{\ell-1}\left(\xi^{g_{\ell}g_{\ell-1}}\right)\sigma_{\ell-2}\right.\\\
\left.\sum_{g_{\ell-2}}\rho_{\ell-2}\left(\xi^{g_{\ell-1}g_{\ell-2}}\right)\cdots\sigma_{1}\sum_{g_{1}}\rho_{1}\left(\xi^{g_{2}g_{1}}\right)\mathbf{x}\right.-\\\
\left.\sigma_{\ell-1}\sum_{g_{\ell-1}}\tilde{\rho}_{\ell-1}\left(\xi^{g_{\ell}g_{\ell-1}}\right)\sigma_{\ell-2}\sum_{g_{\ell-2}}\tilde{\rho}_{\ell-2}\left(\xi^{g_{\ell-1}g_{\ell-2}}\right)\right.\\\
\left.\cdots\sigma_{1}\sum_{g_{1}}\tilde{\rho}_{1}\left(\xi^{g_{2}g_{1}}\right)\mathbf{x}\right\|.$
(76)

To expand and simplify (76), we point that if $A_{\ell+1}$ is the filter in
the layer $\ell+1$, we have

$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})),$
(77)

where the tilde indicates perturbation. Additionally, we take into account
that $\|\sigma_{\ell}(a)-\sigma_{\ell}(b)\|\leq C_{\ell}\|a-b\|$,
$\|A_{\ell+1}-\tilde{A}_{\ell+1}\|\leq\boldsymbol{\Delta}_{\ell}$ and
$\|A_{\ell+1}\|\leq B_{\ell+1}$. We synthesize these facts to write:

$\sum_{g_{k}}\|\alpha-\tilde{\alpha}\|\leq\sum_{g_{k}}\left(\boldsymbol{\Delta}_{k}\|\sigma_{k-1}(\alpha)\|+B_{k}C_{k-1}\|\beta-\tilde{\beta}\|\right),$
(78)
$\sum_{g_{k}}\|\beta-\tilde{\beta}\|\leq\sum_{g_{k}}\sum_{g_{k-1}}\|\alpha-\tilde{\alpha}\|,$
(79)
$\sum_{g_{k}}\|\sigma_{k-1}(\alpha)\|\leq\left(\prod_{r=1}^{k}F_{r}\right)\left(\prod_{r=1}^{k-1}C_{r}B_{r}\right)\|\mathbf{x}\|,$
(80)

where $\alpha$ and $\tilde{\alpha}$ represent sequences of terms in (76) that
start with a symbol of the type $\rho_{\ell}$, and $\beta$ and $\tilde{\beta}$