Quaternion Convolutional Neural Networks: Current Advances and Future
  Directions
Since their first applications, Convolutional Neural Networks (CNNs) have
solved problems that have advanced the state-of-the-art in several domains.
CNNs represent information using real numbers. Despite encouraging results,
theoretical analysis shows that representations such as hyper-complex numbers
can achieve richer representational capacities than real numbers, and that
Hamilton products can capture intrinsic interchannel relationships. Moreover,
in the last few years, experimental research has shown that Quaternion-Valued
CNNs (QCNNs) can achieve similar performance with fewer parameters than their
real-valued counterparts. This paper condenses research in the development of
QCNNs from its very beginnings. We propose a conceptual organization of current
trends and analyze the main building blocks used in the design of QCNN models.
Based on this conceptual organization, we propose future directions of
research.
Gerardo Altamirano-Gomez, Carlos Gershenson
# Quaternion Convolutional Neural Networks: Current Advances and Future Directions

###### Abstract

Since their first applications, Convolutional Neural Networks (CNNs) have solved problems that have advanced the state-of-the-art in several domains. CNNs represent information using real numbers. Despite encouraging results, theoretical analysis shows that representations such as hyper-complex numbers can achieve richer representational capacities than real numbers, and that Hamilton products can capture intrinsic interchannel relationships. Moreover, in the last few years, experimental research has shown that Quaternion-Valued CNNs (QCNNs) can achieve similar performance with fewer parameters than their real-valued counterparts. This paper condenses research in the development of QCNNs from its very beginnings. We propose a conceptual organization of current trends and analyze the main building blocks used in the design of QCNN models. Based on this conceptual organization, we propose future directions of research.

 deep learning, quaternion algebra, computer vision, natural language processing

## I Introduction

In the last decade, the use of deep learning models has become ubiquitous for solving difficult and open problems in science and engineering. Convolutional Neural Networks (CNN) were one of the first deep learning models [1, 2, 3], and their success in tackling the large scale object recognition and classification problem (Imagenet challenge) [4], led to its application in other domains.

The core components of a CNN architecture are the convolution and pooling layers. A convolution layer is as a variation of a fully connected layer (Perceptron), as shown in Figure 1. In the former case, a weight-sharing mechanism over locally connected inputs is applied [5]. This technique is inspired by the local receptive fields discovered by Hubel and Wiesel in their experiments with macaques [6].

Formally speaking, the convolution layer applies the mathematical definition of convolution between discrete signals; thus, for a bi-dimensional input:

\[\mathbf{D}=[d(x,y)]\in\mathbb{R}^{N_{1}\times N_{2}}, \tag{1}\]

the convolution with a kernel:

\[\mathbf{W}=[w(x,y)]\in\mathbb{R}^{M_{1}\times M_{2}}, \tag{2}\]

is defined as follows:

\[\mathbf{F} = \mathbf{W}*\mathbf{D} \tag{3}\] \[= \sum_{r=-\infty}^{\infty}\sum_{s=-\infty}^{\infty}\left[w(r,s)d(x- r,y-s)\right];\]

thus, \(\mathbf{F}\in\mathbb{R}^{N_{1}-M_{1}+1\times N_{2}-M_{2}+1}\) is called _feature map_.

The convolution layer is typically followed by a pooling layer; this provides a sort of local invariance to small rotations and translations of the input features [7, 8]. Moreover, T. Poggio proved that the combination of convolution and pooling layers produce an invariant signature to a group of geometric transformations [9, 10, 11]. However, in the design of very deep architectures, researchers have encountered some difficulties, e.g. reducing the number of parameters without losing generalization, and finding fast optimization methods for adjusting millions of parameters avoiding the vanishing and exploding gradient problems [12, 13].

Fundamental theoretical, as well as experimental analysis, have shown that some algebraic systems, different from the real numbers, have the potential to solve these problems. For example, using a complex numbers representation avoids local minima caused by the hierarchical structure [14], exhibits better generalization [15], and faster learning [16]. Because of the Cayley-Dickson construction [17], it could be inferred that these properties would hold on quaternion-valued neural networks. Recent experimental work has favored this conjecture, where quaternion-valued neural networks show a reduction in the number of parameters, and improved classification accuracies compared to its real-valued counterparts [18, 19, 20, 21, 22, 23, 24]. In addition, a quaternion representation can deal with four-dimensional signals as a single entity [25, 26, 27, 28, 29], models efficiently 3D transformations [30, 31], captures internal latent relationships via the Hamilton product [32], among other properties.

Because of the diversity of deep learning models, this paper focuses in those using quaternion convolution as the main component. Consequently, we have identified three dominant conceptual trends: the classic, the geometric, and the equivariant approaches. These differ in the definition and interpretation of the quaternion convolution layer. The main contributions of this paper can be summarized as follows:

1. This paper presents a classification of QCNNs models based on the definition of quaternion convolution.
2. This paper provides a description of all atomic components needed for implementing QCNNs models, the motivation behind each component, the challenges when they are applied in the quaternion domain, and future directions of research for each component.
3. This paper presents an organized overview of the models that have been found in the literature. They are organized by application domain, classified in: classic, geometric, or equivariant approach; and presented by the type of model: recurrent, residual, convNet, generative or CAE.

The organization of the paper is as follows: In Section II,
we introduce fundamental concepts of quaternion algebra, then in Section III we explain the seminal works on each of the three conceptual trends (classic, geometric, and equivariant), followed by a presentation of the key atomic components to construct QCNNs architectures in Section IV. Thereafter, in Section V, we show a classification of current works by application, and present the diverse types of architectures. Finally, Section VI presents open issues and guidelines for future work, followed by the conclusions in Section VII.

## II Quaternion algebra

This mathematical system was developed by W.R. Hamilton (1805-1865) at the middle of the XIX century [33, 34, 35]. His work on this subject started by exploring ratios between geometric elements, consequently he called _quaternion_ to the quotient of two vectors. He also noticed that for two similar triangles lying on a common plane, \(AOB\) and \(COD\), which are similarly turned, see Figure 2a), the following equality holds [34, pp. 112]:

\[\vec{OB}:\vec{OA}=\vec{OD}:\vec{OC} \tag{4}\]

or expressed as _geometric fractions_:

\[\frac{\vec{OB}}{\vec{OA}}=\frac{\vec{OD}}{\vec{OC}}, \tag{5}\]

where \(\vec{OA}\), \(\vec{OB}\), \(\vec{OC}\), and \(\vec{OD}\), are vectors from \(O\) to point \(A\), \(B\), \(C\), and \(D\), respectively.

Then, by making the triangle \(COD\) into \(BOA^{\prime}\), see Figure 2b); where \(\vec{OA}\) and \(\vec{OB}\) represents any two equally long vectors, we obtain the following relationship [34, pp. 130]:

\[\frac{\vec{OB}}{\vec{OA}}=\frac{\vec{OA^{\prime}}}{\vec{OB}} \tag{6}\]

multiplying by \(\vec{OB}/\vec{OA}\):

\[\left(\frac{\vec{OB}}{\vec{OA}}\right)^{2} = \left(\frac{\vec{OA^{\prime}}}{\vec{OB}}\right)\left(\frac{\vec{OB }}{\vec{OA}}\right)\] \[\left(\frac{\vec{OB}}{\vec{OA}}\right)^{2} = \frac{\vec{OA^{\prime}}}{\vec{OA}} \tag{7}\]

Since \(\vec{OA}\) and \(\vec{OA^{\prime}}\) have the same magnitude, but opposite direction, \(\vec{OA}=-\vec{OA^{\prime}}\), we obtain:

\[\left(\frac{\vec{OB}}{\vec{OA}}\right)^{2}=-1. \tag{8}\]

The quotient of two perpendicular equally long vectors, like in this case, is called _right radial quaternion_. Since not particular assumption was made about quotients in Equation (8), e.g. the plane of that quotient is arbitrary, Hamilton concluded that every right radial quaternion, was one of the square roots of negative unity [34, pp. 131]. In addition, for a quaternion: \(\textbf{q}=\vec{v_{1}}:\vec{v_{2}}\), where \(\vec{v_{1}}\) and \(\vec{v_{2}}\) are vectors, we can rewrite it as \(\textbf{q}\vec{v_{1}}=\vec{v_{2}}\). In this case, **q** is called a _versor_, i.e. an element that transforms \(\vec{v_{1}}\) into \(\vec{v_{2}}\) by rotating it [34, pp. 133]. In this way, these results connected quaternion representation with the root \(\sqrt{-1}\), and its geometric meaning.

Now, considering Figure 3, let \(\textbf{q}=\vec{OB}:\vec{OA}\) be a quaternion, and let \(\vec{OB^{\prime}}\) and \(\vec{OB^{\prime\prime}}\) be parallel and perpendicular vectors to \(\vec{OA}\), respectively; such that \(\vec{OB}=\vec{OB^{\prime}}+\vec{OB^{\prime\prime}}\). Then, we can decompose the quaternion, **q**, into:

\[\textbf{q}=\vec{OB^{\prime}}:\vec{OA}+\vec{OB^{\prime\prime}}:\vec{OA} \tag{9}\]

Since \(\vec{OB^{\prime}}\) and \(\vec{OA}\) are parallel vectors, their quotient is just a scale factor representing the projection of \(\vec{OA}\) into \(\vec{OB^{\prime}}\), so first term turns into a scalar. The second term represents the projection of \(\vec{OA}\) on the plane through \(O\), which is perpendicular to \(\vec{OA}\); in addition, this means \(\vec{OB^{\prime\prime}}\) can be obtained from \(\vec{OA}\) by applying a versor transformation. Since \(\vec{OB^{\prime\prime}}\) and \(\vec{OA}\) are perpendicular to each other, its quotient is called _right quaternion_ (note that it differs from right radial

Fig. 1: Fully connected layer (a) vs. convolutional layer (b). For an input array of \(2\times 3\) elements, the Perceptron uses \(12\) weights and \(24\) connections (not locally connected,) while the convolution layer uses \(4\) weights and \(8\) connections. In the convolutional layer, we have a reduction in the number of parameters and connections, because a single weight is connected to several inputs, and the weights are applied over locally connected inputs.
quaternion in the length of the vectors), and can be expressed as a linear combination of right radial quaternions (versors). This leads to the well know result that every quaternion can be reduced to the quadrinomial form [34, pp. 160]:

\[\textbf{q}=q_{R}+q_{I}\hat{i}+q_{J}\hat{j}+q_{K}\hat{k}, \tag{10}\]

where \(q_{R},q_{I},q_{J},q_{K}\) are scalars, and \(\hat{i},\hat{j},\hat{k}\) are three right versors.

In terms of modern mathematics, the quaternion algebra, \(\mathbb{H}\), is: the 4-dimensional vector space over the field of the real numbers, generated by the basis \(\{1,\hat{i},\hat{j},\hat{k}\}\), and endowed with the following multiplication rules (Hamilton product):

\[(1)(1)= 1\] \[(1)(\hat{i})= \hat{j}\hat{k}= -\hat{k}\hat{j}= \hat{i}\] \[(1)(\hat{j})= \hat{k}\hat{i}=-\hat{i}\hat{k}= \hat{j}\] \[(1)(\hat{k})= \hat{i}\hat{j}=-\hat{j}\hat{i}= \hat{k}\] \[\hat{i}^{2}= \hat{j}^{2}= \hat{k}^{2}= -1 \tag{11}\]

Therefore, this definition makes the quaternion algebra associative and non-commutative. Thus, for two arbitrary quaternions: \(\textbf{p}=p_{R}+p_{I}\hat{i}+p_{J}\hat{j}+p_{K}\hat{k}\) and \(\textbf{q}=q_{R}+q_{I}\hat{i}+q_{J}\hat{j}+q_{K}\hat{k}\), their multiplication is calculated as follows:

\[\textbf{p}\textbf{q} = p_{R}q_{R}-p_{I}q_{I}-p_{J}q_{J}-p_{K}q_{K} \tag{12}\] \[+(p_{R}q_{I}+p_{I}q_{R}+p_{J}q_{K}-p_{K}q_{J})\hat{i}\] \[+(p_{R}q_{J}-p_{I}q_{K}+p_{J}q_{R}+p_{K}q_{I})\hat{j}\] \[+(p_{R}q_{K}+p_{I}q_{J}-p_{J}q_{I}+p_{K}q_{R})\hat{k}\]

Notice that each coefficient of the resulting quaternion, is composed of real and imaginary parts of the factors \(p\) and \(q\). In this way, the Hamilton product capture interchannel relationships between both factors.

Next, there are introduced some useful operations with quaternions.

Let, \(\textbf{q}=q_{R}+q_{I}\hat{i}+q_{J}\hat{j}+q_{K}\hat{k}\), be a quaternion, its _conjugate_ is defined as:

\[\bar{\textbf{q}}=q_{R}-q_{I}\hat{i}-q_{J}\hat{j}-q_{K}\hat{k} \tag{13}\]

And its magnitude is computed as follows:

\[\|\textbf{q}\|=\sqrt{\textbf{q}\bar{\textbf{q}}} \tag{14}\]

As well as complex numbers, quaternions can be represented in polar form [36, 37]:

\[\textbf{q}=\|\textbf{q}\|\left[\cos(\theta)+\sin(\theta)\frac{q_{I}\hat{i}+q_ {J}\hat{j}+q_{K}\hat{k}}{\|q_{I}\hat{i}+q_{J}\hat{j}+q_{K}\hat{k}\|}\right], \tag{15}\]

where:

\[\theta=atan\left(\frac{\sqrt{q_{I}^{2}+q_{J}^{2}+q_{K}^{2}}}{q_{R}}\right). \tag{16}\]

Even though other polar parametrizations has been proposed, see for example [38, 39, 25, 40], current QCNNs apply the former equation.

As was mention before, a quaternion can represent a geometric transformation. Let \(\textbf{w}_{\theta}\) be a unitary versor expressed in polar form:

\[\textbf{w}_{\theta}=\cos(\theta)+\sin(\theta)(w_{I}\hat{i}+w_{J}\hat{j}+w_{K} \hat{k}), \tag{17}\]

and let \(\textbf{q}\) be any quaternion; then, their multiplication:

\[\textbf{p}=\textbf{w}_{\theta}\textbf{q} \tag{18}\]

Fig. 3: Given two vectors \(O\bar{A}\) and \(O\bar{B}\), we construct the following geometric configuration: plane \(\Pi\) is orthogonal to \(O\bar{A}\), vector \(O\bar{B}^{\prime\prime}\) lies on the plane \(OAB\) and is the projection of \(O\bar{B}\) into \(\Pi\), while vector \(O\bar{B}^{\prime}\) is the projection of \(O\bar{B}\) into \(O\bar{A}\). Adapted from [34].

Fig. 2: Two similar triangles, turned and in a common plane. a) General configuration b) Right triangles with equal length cathetus configuration. Adapted from [34].
applies a rotation, with angle \(\theta\), along an axis \(\mathbf{w}=w_{I}\hat{i}+w_{J}\hat{j}+w_{K}\hat{k}\).

Since we have not made any particular assumption on \(\mathbf{q}\), we are applying a rotation in the four-dimensional space of quaternions. Denoting by \(\Pi_{1}\) to the plane defined by the scalar axis and the vector \(w_{I}\hat{i}+w_{J}\hat{j}+w_{K}\hat{k}\), and \(\Pi_{2}\) the plane perpendicular to \(w_{I}\hat{i}+w_{J}\hat{j}+w_{K}\hat{k}\); it can be proved that this 4-dimensional rotation is composed of a simultaneous rotation of the elements on plane \(\Pi_{1}\), and of the elements on plane \(\Pi_{2}\)[37].

Alternatively, we can split the transformation as a sandwiching product:

\[\mathbf{p}=\mathbf{w}_{\frac{q}{2}}\mathbf{q}\mathbf{\bar{w}}_{\frac{q}{2}}; \tag{19}\]

in this case, the angle of each versor, \(\mathbf{w}_{\theta/2}\), is divided to half.

From the group theory perspective, the set of unitary versors lies on a 3-Sphere, \(\mathbb{S}^{3}\), embedded in a 4D Euclidean space [30]; and together with the Hamilton product form a group, which is isomorphic to the 4D rotation group SO(4) [37], and to the special unitary group SU(2) [41, 30]. Moreover, there exist a two to one homomorphism with the rotation group SO(3) [41, 30].

Another operation of interest is quaternion convolution; because the non-commutative property of quaternion multiplication, we have 3 different definitions of discrete quaternion convolution. First, left-side quaternion convolution, defined as follows [42, 43, 44]:

\[(\mathbf{w}\ast\mathbf{q})(x,y)=\sum_{r=-\frac{L}{2}}^{\frac{L}{2}}\sum_{s=- \frac{L}{2}}^{\frac{L}{2}}[\mathbf{w}(r,s)\mathbf{q}(x-r,y-s)] \tag{20}\]

Second, right-side quaternion convolution, defined as follows [43]:

\[(\mathbf{q}\ast\mathbf{w})(x,y)=\sum_{r=-\frac{L}{2}}^{\frac{L}{2}}\sum_{s=- \frac{L}{2}}^{\frac{L}{2}}[\mathbf{q}(x-r,y-s)\mathbf{w}(r,s)] \tag{21}\]

In third place, we have two-sided quaternion convolution [43, 44, 45]:

\[(\mathbf{w_{left}}*\mathbf{q}*\mathbf{w_{right}})(x,y)=\] \[\sum_{r=-\frac{L}{2}}^{\frac{L}{2}}\sum_{s=-\frac{L}{2}}^{\frac{ L}{2}}[\mathbf{w_{left}}(r,s)\mathbf{q}(x-r,y-s)\mathbf{w_{right}}(r,s)] \tag{22}\]

where \(q,w,w_{left},w_{right}\in\mathbb{H}\), \(\ast\) represents the convolution operator, and is applied quaternion product between \(\mathbf{q}\)'s and \(\mathbf{w}\)'s.

Finally, Table I summarizes the notation that will be used in the rest of this paper.

## III Development of the classic, geometric, and equivariant approaches.

In the previous section, it was presented the different definitions of quaternion convolution. From this definition, different conceptual approaches can be obtained by setting some restrictions on the quaternions, e.g. using unitary quaternions, or using just the real part of one of the quaternions. In the following paragraphs, we describe the seminal works that led to the development of the three main conceptual treats found in most of the current works on QCNNs. Since these approaches were not developed in an incremental manner, we focus on tracking the seminal ideas of each approach, describing what components were introduced, and specifying the domain of application in which they were tested. Thereafter, in Section IV-A, it is presented the formal definition of each approach.

In the _classic approach_, the definition of quaternion convolution is a natural extension of the real and complex convolution; its role is to compute the correlation between input data and kernel patterns, in the quaternion domain. First works refer to Altamirano [46], whom based on the work on Quaternion-Valued Multilayer Perceptrons (QMLP's) by Arena _et al._[47, 48, 18], defines the main components of a QCNN: quaternion convolution layers, quaternion pooling layers that use the magnitude of the quaternion, quaternion split-RELU activation function, quaternion fully connected layers, and the quaternion back-propagation training method. As a proof of concept, these components were applied to a simple pattern classification problem. In an independent manner, Gaudet and Maida [19] arrived to a similar definition of convolution layers, but their work is based on Travelsi's research on Complex-Valued CNNs [49]; in addition, they proposed quaternion bath-normalization and weight initialization algorithms. Their models were tested by classifying images of the CIFAR-10 and CIFAR-100 databases [50], as well as the KITTI Road Estimation Benchmark [51]. At the same time, and also following the work of Travelsi [49], T. Parcollet implemented quaternion-valued CNNs [22], and LSTM models [52]; these were applied on a voice recognition task using the TIMIT Dataset [53]. In addition, they proposed a decoding-encoding model that converts grayscale images from the KODAK PhotoCD Dataset [54] to RGB images.

Thereafter, Yin _et al._[55] proposed a quaternion attention mechanism and evaluates it on the image classification task using the CIFAR-10 dataset [50]. Moreover, he proposed a model for detection of double compression JPEG images and tested it on the Uncompressed Color Image Database (UCID) dataset [56]. Their models implement a pooling method that uses the magnitude of the quaternion, quaternion split activation functions, and an alternative bath normalization method, which reduces the computational cost of other methods by using a single variance value instead of a multichannel covariance matrix.

Alternatively, we have the _geometric approach_ which was
constructed based on the previous work of Matsui _et al._[57] and Isokawa _et al._[58]. In this case, the quaternion product applies affine transformations over the input features, consequently the quaternion convolution inherits this property. Thus, Zhu _et al._[24] defines quaternion convolution based on a geometric transformation that applies fixed-axis rotation and scaling (2DoF); moreover, they apply the same concept for constructing quaternion fully connected layers. Their model was tested on the image classification problem using the CIFAR-10, CIFAR-100 [50], and Oxford flowers [59] datasets. In addition, they tested their models for the noise elimination task using the Oxford flowers dataset [59] and a subset of Microsoft COCO dataset [60]. Under the same approach, Hongo _et al._[20] define a quaternion convolution layer that applies affine transformations; in addition, they applied quaternion pooling layers using the magnitude of the quaternion, quaternion split-RELU activation function, and the batch-normalization method of [55]. From these concepts, they implemented QCNNs and residual QCNNs for the image classification problem on the CIFAR-10 dataset [50]. In a recent work, Matsumoto _et al._[61] note the limited expression ability of working with fixed axes, and propose a model that learns general rotation axis (4DoF). They applied this model in a pixel classification task with PolSAR images obtained from the Japan Aerospace Exploration Agency (JAXA).

Finally, the third approach is based on the concept of rotation _equivariance_, i.e. if an input feature is rotated a specific angle, then the output feature produced by the convolution layer is equivalent to: take the non-rotated input feature, apply the convolution layer, and then apply the rotation to the output feature. In order to satisfy the rotation equivariant property, Shen _et al._[23] define a convolution layer involving products between quaternion inputs and real-valued kernels. In addition, they define a rotation equivariant RELU layer and a batch normalization algorithm. Their models were applied to 3D point cloud classification on the ModelNet40 [62] and 3D MNIST [63] datasets. They show that these type of networks are robust to rotated input features, and the feature maps produced by inner layers are invariant to permutations of the input data points.

Another work lying in this category is the one by Jing _et al._[64], who uses the middle element of a convolutional window as a pivot, and define a convolution-like operation that produces rotation equivariant features. By taking the magnitude of the quaternion output, the convolutional blocks can be used to construct rotation invariant classifiers.

Most work on Quaternion-valued CNN is based on these, or lies in one of the preceding categories; in the following section, we will give the formal definition of each approach, and will explain the key atomic components used to construct QCNNs.

## IV QCNN components

In this section, we present the main building blocks for implementing quaternion-valued convolutional deep learning architectures. Future directions for each component are specified at the end of each subsection.

### _Quaternion convolution layers_

Lets assume a dataset, where each sample has dimension \(N\times M\times 4\), i.e. an input can be represented as a 4-channels matrix of real numbers. Then, each sample, \(\mathbf{Q}\), is represented as a \(N\times M\) matrix where each element is a quaternion:

\[\mathbf{Q}=[\mathbf{q}(x,y)]\in\mathbb{H}^{N\times M} \tag{23}\]

then, \(\mathbf{Q}\) can be decomposed in its real and imaginary components:

\[\mathbf{Q}=Q_{R}+Q_{I}\hat{i}+Q_{J}\hat{j}+Q_{K}\hat{k} \tag{24}\]

where \(Q_{R},Q_{I},Q_{J},Q_{K}\in\mathbb{R}^{N\times M}\), and \(\hat{i},\hat{j},\hat{k}\) represent the complex basis of the quaternion algebra.

In the same way, a convolution kernel of size \(L\times L\) is represented by a quaternion matrix, as follows:

\[\mathbf{W}=[\mathbf{w}(x,y)]\in\mathbb{H}^{L\times L} \tag{25}\]

which can be decomposed as:

\[\mathbf{W}=W_{R}+W_{I}\hat{i}+W_{J}\hat{j}+W_{K}\hat{k}, \tag{26}\]

where \(W_{R},W_{I},W_{J},W_{K}\in\mathbb{R}^{L\times L}\), and \(\hat{i},\hat{j},\hat{k}\) represent the basis of the quaternion algebra.

Then, in the _classic approach_, Altamirano [46] and Gaudet and Maida [19] define the convolution layer using left-sided convolution:

\[\mathbf{F}=\mathbf{W}*\mathbf{Q}. \tag{27}\]

Thus, \(\mathbf{F}\in\mathbb{H}^{(N-L+1)\times(M-L+1)}\) represents the output of the layer, i.e. a quaternion feature map, and each element of the tensor is computed as follows:

\[\mathbf{f}(x,y)=(\mathbf{w}*\mathbf{q})(x,y). \tag{28}\]

This approach does not make any particular assumption about quaternions \(\mathbf{w}\) and \(\mathbf{q}\). Thus, the convolution represents the integral transformation of a quaternion function on a quaternion input signal.

On contrast, in the _geometric approach_, Zhu _et al._[24] apply the two-sided convolution definition:

\[\mathbf{F}=\mathbf{W}*\mathbf{Q}*\bar{\mathbf{W}} \tag{29}\]

where \(\mathbf{F}\in\mathbb{H}^{(N-L+1)\times(M-L+1)}\), and each element of the output is computed as follows:

\[\mathbf{f}(x,y)=\sum_{r=-\frac{L}{s}}^{\frac{L}{2}}\sum_{s=-\frac{L}{2}}^{\frac {L}{2}}\frac{\mathbf{w}(r,s)\mathbf{q}(x-r,y-s)\bar{\mathbf{w}}(r,s)}{a_{r,s}}. \tag{30}\]

In addition, each quaternion \(\mathbf{w}\), is represented in its polar form:

\[\mathbf{w}(r,s)=a_{r,s}\left(\cos\frac{\theta_{r,s}}{2}+\sin\frac{\theta_{r,s }}{2}\bar{u}\right) \tag{31}\]

where \(\theta_{r,s}\in[-\pi,\pi]\), \(a_{r,s}\in\mathbb{R}\), and \(\bar{u}\) represents the unitary rotation axis. Thus, quaternion convolution applies rotation and scaling transformations on the quaternion \(\mathbf{q}\).

Similarly, Hongo _et al._[20] use the two-sided quaternion convolution definition, but adds a threshold value for each component of the quaternion:

\[\mathbf{F}=\mathbf{W}*\mathbf{Q}*\bar{\mathbf{W}}+\mathbf{B}. \tag{32}\]
where \(\mathbf{B}\in\mathbb{H}\).

In both cases, quaternion convolution applies geometric transformations over the input data.

Finally, in the _equivariant approach_ Shen _et al._[23] propose a simplified version; they convolutes the quaternion input, \(\mathbf{Q}\in\mathbb{H}\), with a kernel of real numbers, \(\mathbf{W}\in\mathbb{R}\):

\[\mathbf{F} = \mathbf{W}*\mathbf{Q} \tag{33}\] \[= \mathbf{W}*\mathbf{Q_{R}}+\mathbf{W}*\mathbf{Q_{I}}\hat{i}+ \mathbf{W}*\mathbf{Q_{J}}\hat{j}\] \[+\mathbf{W}*\mathbf{Q_{K}}\hat{k}.\]

This version, allows to extract equivariant features, e.i. if an input sample, \(\mathbf{Q}\), produce a feature map, \(g(\mathbf{Q})\), then the rotated input sample, \(R(\mathbf{Q})\), will produce a rotated feature map, \(g(\mathbf{Q})\), thus:

\[g(R(\mathbf{Q}))=R(g(\mathbf{Q})). \tag{34}\]

Independently of the approach that we follows, a related problem when we implement QCNNs, is how to deal with multidimensional inputs. In real-valued CNNs, the common way of dealing with them is: defining 2D convolution kernels with the same number of channels as the input data, apply 2D convolution separately for each channel; then, the resulting 2D outputs are summed over all channels to produce a single-channel feature map. A different approach is to apply N-dimensional convolution; in this case, the multidimensional kernel is convoluted over all the channels of the input data.

For quaternion-valued CNNs, current implementations use variations of the first approach; i.e. the input data is divided into 4-channel sub-inputs, thereafter quaternion convolution is computed for each sub-input. Before explaining the details, some notation is introduced.

Let \(\mathbf{X}\in\mathbb{R}^{N\times M\times C}\) be an input data, \(N\) is the number of rows, \(M\) the number of columns, and \(C\) is the number of channels, where \(C\%4=0\), then \(X\) is partitioned as follows:

\[\mathbf{X}=[\mathbf{Q_{0}},\mathbf{Q_{1}},\dots,\mathbf{Q_{(C/4)-1}}] \tag{35}\]

where each \(\mathbf{Q_{s}}\in\mathbb{H}^{N\times M}\), \(0<s<(C/4)-1\) is a _quaternion channel_.

Let \(\mathbf{V}\in\mathbb{R}^{L\times L\times K}\), with \(K\%4=0\), be the convolution kernel, then:

\[\mathbf{V}=[\mathbf{W_{0}},\mathbf{W_{1}},\dots,\mathbf{W_{(K/4)-1}}] \tag{36}\]

where each \(\mathbf{W_{s}}\in\mathbb{H}^{L\times L}\), \(0<s<(K/4)-1\) is a quaternion channel.

Thus, there are three different ways of dealing with multi-dimensional quaternion inputs, see Figure 4:

1. Autoencoder convolution: Kernel and input have the same number of channels (\(K=C\)). In this case, each quaternion channel input is assigned to one quaternion channel kernel, and convolution between them is computed using Equation (27), (29), (32) or (33); Thus, if we use left-sided convolution, each individual output, \(\mathbf{F_{s}}\in\mathbb{H}^{N-L+1\times M-L+1}\), is computed as follows: \[\mathbf{F_{s}}=\mathbf{W_{s}}*\mathbf{Q_{s}},\] (37) where \(0<i<(C/4)-1\), and the final quaternion feature map is obtained by concatenating all outputs: \[\mathbf{F}=[\mathbf{F_{0}},\mathbf{F_{2}},\dots,\mathbf{F_{(C/4)-1}}]\] (38) This method produces an output with the same number of channels as the input data; and could be used in Convolutional Auto-Encoders (CAE).
2. Pyramidal convolution: Kernel and input have different number of channels (\(K\neq C\)), but they are multiples of \(4\). In [46], it is proposed the computing of feature maps using a pyramidal approach: each kernel \(\mathbf{W_{t}}\), where \(0<t<(K/4)-1\) is convolved with each sub-input \(\mathbf{Q_{s}}\), where \(0<s<(C/4)-1\), hence, it produces \(C/4\) quaternion outputs. Since each quaternion input channel is convolved with each quaternion kernel channel, see Figure 4, a convolution kernel \(\mathbf{V}\in\mathbb{R}^{L\times L\times K}\) will produce \(C*K/16\) quaternion outputs: \[\mathbf{F}=[\mathbf{F_{0}},\mathbf{F_{2}},\dots,\mathbf{F_{C-K/16}}].\] (39) Thus, if left-sided convolution is applied, each quaternion output channel, \(\mathbf{F_{k}}\), is computed as follows: \[\mathbf{F_{(t+C/4)+s}}=\mathbf{W_{t}}*\mathbf{Q_{s}}.\] (40) Even though [46] used left-sided convolution, this approach is valid with any of the other convolution definitions. The intuition behind this approach is to detect a quaternionic pattern in any sub-input; however, its application is impractical Beacuse of the exponential growth of the number of channels in deep architectures. Calculation of summed outputs can alleviate the computational cost.
3. Summed convolution: Similarly to the former method, but in this approach, each quaternion input channel is convolved with a different quaternion kernel channel: \[\mathbf{F_{s}}=\mathbf{W_{s}}*\mathbf{Q_{s}},\] (41) where \(0<s<(C/4)-1\). Then, the final quaternion feature map, \(\mathbf{F}\in\mathbb{H}^{N\times M}\), is obtained by summing all outputs: \[\mathbf{F}=\sum_{s=0}^{(C/4)-1}\mathbf{F_{s}}.\] (42)

_Future directions:_ As was mentioned before, N-dimensional convolution has been applied on Real-Valued CNNs for processing multichannel inputs. In a similar way, the equations presented can be extended to 8-channel inputs using an octonion or hypercomplex algebra, see for example [65, 66, 67, 68]. For larger number of inputs, a geometric algebra [69] representation can be applied. To the best of our knowledge, this type of architecture has not been published to date, but a first approach in this direction can be found in [70, 71]. In these type of deep learning architectures: quaternion, hyper-complex, or geometric, a major concern is the selection of the signature of the algebra, which will embed data into different geometric spaces, and the processing will take distinct meanings accordingly. Thus, a sensible selection of the dimension and signature should be made according to the nature of the problem and the meaning of input data as well as intermediate layers.
### _Quaternion Fully Connected Layers_

Let \(\mathbf{Q}\) be a \(N_{1}\times N_{2}\times N_{3}\) tensor, representing the input to a fully connected layer; then each element of \(\mathbf{Q}\) is a quaternion:

\[\mathbf{Q}=[\mathbf{q}(x,y,z)]\in\mathbb{H}^{N_{1}\times N_{2}\times N_{3}}, \tag{43}\]

where \(N_{1},N_{2},N_{3}\) are the height, width, and number of channels of the input.

Now, for the fully connected layer, it is defined a quaternion kernel, \(\mathbf{W}\), of size \(N_{1}\times N_{2}\times N_{3}\), where each element is a quaternion:

\[\mathbf{W}=[\mathbf{w}(x,y,z)]\in\mathbb{H}^{N_{1}\times N_{2}\times N_{3}}. \tag{44}\]

where \(N_{1},N_{2},N_{3}\) are the height, width, and number of channels of the input.

Note that elements of the input and weight tensors are denoted as \(q(x,y,z)\), and \(w(x,y,z)\), respectively; and the output of the layer will be a quaternion, \(\mathbf{f}\in\mathbb{H}\).

Thus, for classic fully connected layers, the output, \(\mathbf{f}\), is computed as follows [46]:

\[\mathbf{f}=\sum_{r,s,t}^{N_{1},N_{2},N_{3}}\left[\mathbf{w}(r,s,t)\mathbf{q}( r,s,t)\right]. \tag{45}\]

Similarly, for geometric fully connected layers, the output is computed as follows [24]:

\[\mathbf{f}=\sum_{r,s,t}^{N_{1},N_{2},N_{3}}\left[\frac{1}{\|\mathbf{w}(r,s,t) \|}\mathbf{w}(r,s,t)\mathbf{q}(r,s,t)\bar{\mathbf{w}}(r,s,t)\right]. \tag{46}\]

The difference between quaternion-valued and real-valued fully connected layers relies on the application of the Hamilton product, which captures interchannel relationships, and in the former case, the output is a quaternion. Moreover, it should be noticed that for real-based networks, fully connected layers are equivalent to inner product layers, but for quaternion-valued networks, the output of quaternion fully connected layers and quaternion inner product layers are different.

_Future directions:_ Current implementations of fully quaternion layers follows a classic or geometric approach; the equivariance approach should be implemented and tested.

### _Quaternion Pooling_

Most of the current methods for applying quaternion pooling rely on channel-wise pooling, in a similar way as is applied in Real-Valued CNNs. For example, [19, 61] use channel-wise global average pooling layers, and [72, 73, 24] applied channel-wise average as well as channel-wise max pooling layers, while [20, 74, 22] applied just channel-wise max pooling layers, defined as follows:

\[splitPool(\mathbf{Q}) = \underset{(x,y)}{\max}(q_{R}(x,y))+\underset{(x,y)}{\max}(q_{I}(x, y))\hat{i} \tag{47}\] \[+\underset{(x,y)}{\max}(q_{J}(x,y))\hat{j}\] \[+\underset{(x,y)}{\max}(q_{K}(x,y))\hat{k},\]

where \(\mathbf{Q}\) is a quaternion submatrix.

In contrast, [23, 55, 46] use a max-pooling approach, but instead of using the channel-wise maximum value, they select the quaternion with maximum amplitude within a region:

\[FullyPool(\mathbf{Q}) = \mathbf{q}(\bar{x},\bar{y})\ s.t.\ (\bar{x},\bar{y}) \tag{48}\] \[= \underset{(x,y)}{\arg\max}(\|\mathbf{q}(x,y)\|).\]

In addition, since this method can obtain multiple maximum amplitude quaternions, [55] applies the angle cosine theorem to discriminate between them.

_Future directions:_ Future works should introduce novel fully quaternion pooling methods, emphasizing the use of interchannel relationships. For example: using polar representations, introducing quaternion measures [75], or taking inspiration from information theory.

### _Quaternion batch normalization_

Internal covariance shift [76], is a statistical phenomenon that occurs during training: the change of the network parameters causes changes in the statistical distribution of the inputs in hidden layers. Whitening procedures [77, 78], i.e. apply linear transformations to obtain uncorrelated inputs with zero means and unit variances, alleviate this phenomenon. Since whitening the input layers is computationally expensive, because of the computing of covariance matrices, its inverse square root, and the derivatives of these transformations, Ioffe and Szegedy [76] introduced the batch normalization algorithm, which normalizes each dimension independently. It uses mini-batchs to estimate the mean and variance of each channel, and transform the channel to have zero mean and unit variance using the following equation:

\[\tilde{x}=\frac{x-E(x)}{\sqrt{\sigma^{2}+\epsilon}} \tag{49}\]

where \(\epsilon\) is a constant added for numerical stability. In order to maintain the representational ability of the network, an affine

Fig. 4: Visualization of different approaches for computing quaternion convolution on a multichannel input. This example shows a 2D input with 2 quaternion channels; the output of an autoencoder convolution layer is shaded in pink, the output of the pyramidal convolution layer is framed in cyan, and the summed convolution layer is shaded in orange.
transformation with two learnable parameters, \(\gamma\) and \(\beta\), is applied:

\[BN(\tilde{x})=\gamma\tilde{x}+\beta \tag{50}\]

Even though this method does not produce uncorrelated inputs, it improves convergence time by enabling higher learning rates, and has allowed the training of deeper neural network models. On the other hand, uncorrelated inputs reduce overfitting and improve generalization [79].

Since channel-wise normalization does not assure equal variance in the real and imaginary components, Gaudet and Maida [19] proposed a quaternion batch-normalization algorithm using the whitening approach [80], and treat each component of the quaternion as an element of a four dimensional vector. We call this approach _Whitening Quaternion Batch-Normalization (WQBN)_. Let \(\mathbf{x}\) be a quaternion input variable, \(\mathbf{x}=[x_{R},x_{I},x_{J},x_{K}]^{T}\), \(E(\mathbf{x})\) its expected value, and \(V(\mathbf{x})\) its quaternion covariance matrix, both computed over a mini-batch, then:

\[V(\mathbf{x})=\begin{bmatrix}v_{rr}&v_{ri}&v_{rj}&v_{rk}\\ v_{ir}&v_{ii}&v_{ij}&v_{ik}\\ v_{jr}&v_{ji}&v_{jj}&v_{jk}\\ v_{kr}&v_{ki}&v_{kj}&v_{kk}\end{bmatrix}, \tag{51}\]

where subscripts represent the covariance between real or imaginary components of \(\mathbf{x}\), e.g. \(v_{ij}=cov(x_{I},x_{J})\). Thus, the Cholesky decomposition of \(V^{-1}\) is computed, and one of the resulting matrices, \(W\), is selected. Thereafter, the whitened quaternion variable, \(\tilde{\mathbf{x}}\), is calculated with the following matrix multiplication [19]:

\[\tilde{\mathbf{x}}=W(\mathbf{x}-E(\mathbf{x})), \tag{52}\]

and finally:

\[WQBN(\tilde{\mathbf{x}})=\mathbf{\Gamma}\tilde{\mathbf{x}}+\beta \tag{53}\]

where \(\beta\in\mathbb{H}\) is a trainable parameter, and:

\[\mathbf{\Gamma}=\begin{bmatrix}\Gamma_{rr}&\Gamma_{ri}&\Gamma_{rj}&\Gamma_{rk} \\ \Gamma_{ri}&\Gamma_{ii}&\Gamma_{ij}&\Gamma_{ik}\\ \Gamma_{rj}&\Gamma_{ij}&\Gamma_{jj}&\Gamma_{jk}\\ \Gamma_{rk}&\Gamma_{ik}&\Gamma_{jk}&\Gamma_{kk}\end{bmatrix}, \tag{54}\]

is a symmetric matrix with trainable parameters.

In constrast, Yin _et al._[55] applied the quaternion variance definition proposed by [81]:

\[V(\mathbf{x})=\frac{1}{T}\sum_{i=1}^{T}\mathbf{v}\bar{\mathbf{v}}, \tag{55}\]

where \(\mathbf{v}=\mathbf{x}-E(\mathbf{x})\). Note that in this case, the variance is a single real value. We call this approach the _Variance Quaternion Batch Normalization (VQBN)_. Thus, the batch normalization is computed as follows:

\[VQBN(\mathbf{x})=\gamma\frac{\mathbf{x}-E(\mathbf{x})}{\sqrt{V(\mathbf{x})^{2 }+\epsilon}}+\beta, \tag{56}\]

where \(\gamma,\epsilon\in\mathbb{R}\), \(\beta\in\mathbb{H}\) are trainable parameters.

Recently, Grassucci _et al._[82] noted that for proper quaternion random variables [83, 84], the covariance matrix in Equation (51) becomes the diagonal matrix:

\[V(\mathbf{x})=4\sigma^{2}I, \tag{57}\]

and the batch-normalization procedure is simplified to applying Equations (55) and (56) with \(V(\mathbf{x})=2\sigma\).

A third approach is the one proposed by [23], who define the batch-normalization operation to be rotation-equivariance, thus:

\[RQBN(\mathbf{x})=\frac{\mathbf{x}}{\sqrt{E(\|\mathbf{x}\|^{2})+\epsilon}}. \tag{58}\]

Generative Adversarial Networks apply another type of normalization, called Spectral Normalization [85]. In these networks, having a Lipschitz-bounded discriminative function is crucial to mitigate the gradient explosion problem [86, 87, 88]; thus, based on their real-valued counterpart, Grassucci _et al._[82] proposed a Quaternion Spectral Normalization algorithm that constrain the spectral norm of each layer. To explain this procedure, we introduce some definitions.

Let \(f\) be a generic function, it is K-Lipschitz continuous if, for any two points, \(x_{1},x_{2}\), it satisfies:

\[\frac{\|f(x_{1})-f(x_{2})\|}{|x_{1}-x_{2}|}\leq K \tag{59}\]

Let \(\sigma(\cdot)\) be the spectral norm of a matrix, i.e. the largest singular value of a matrix, the Lipschitz norm of a function \(f\), denoted by \(\|f\|_{Lip}\), is defined as follows:

\[\|f\|_{Lip}=\sup_{x}\sigma(\nabla f(x)). \tag{60}\]

Thus, for a generic linear layer, \(f(h)=Wx+b\), whose gradient is \(W\), their Lipschitz norm is:

\[\|f\|_{Lip}=\sigma(W). \tag{61}\]

Now, let \(\mathbf{W}\) be a quaternion matrix; their Lipschitz norm, \(\sigma(\mathbf{W})\), is computed by estimating the largest singular value of \(\mathbf{W}\) via the power iteration method [82]. Then, Quaternion Spectral Normalization is applied, in a split way, using the following equations:

\[\bar{W}_{R} = \frac{W_{R}}{\sigma(W)},\bar{W}_{I}=\frac{W_{I}}{\sigma(W)},\bar{ W}_{J}=\frac{W_{J}}{\sigma(W)},\text{ and}\] \[\bar{W}_{K} = \frac{W_{K}}{\sigma(W)}. \tag{62}\]

For applying the Lipschitz bound to the whole network, it is used the following relationship:

\[\|f_{1}\circ f_{2}\|_{Lip}\leq\|f_{1}\|_{Lip}\cdot\|f_{2}\|_{Lip}. \tag{63}\]

provided the Lipschitz norm of each layer is constrained to \(1\)[82].

_Future directions:_ The WQBN algorithm produces uncorrelated, zero mean, and unit variance inputs, but is computationally expensive. In addition, Kessy _et al._[80] states that there are infinitely many possible matrices satisfying \(V^{-1}=\bar{W}W\). In contrast, VQBN algorithm produces zero mean and unit variance input (according to the Wang _et al._ definition [81], which averages the variance of the four channels), but since we use a single value variance, scaling is isotropic, and the input channels still correlated. However, this approach greatly reduces the computational time, by avoiding the decomposition of the covariance matrix. Further theoretical and experimental analysis is required to grasp its advantages versus independent channel batch-normalization.
### _Activation functions_

Biological neurons produce an output signal if a set of input stimuli surpasses a threshold value within a lapse of time; for artificial neurons, the role of the activation function is to simulate this triggering action. Mathematically, this behavior is modeled by a mapping: in the case of real-valued neurons, the domain and image is the field of real numbers, while complex and quaternion neurons map complex or quaternion inputs to complex or quaternion outputs, respectively.

Since back-propagation has become the standard method for training artificial neural networks, it is required for the activation function to be analytic [89], i.e. the derivative exists at any point. In the complex domain, by the Liouville theorem, it is known that a bounded function, which is analytical at any point, is a constant function; reciprocally, complex non-constant functions have non-bounded images [90]. Accordingly, some common real-valued functions, like the hyperbolic tangent and sigmoid, will diverge to infinity when they are extended to the complex domain [89]; this makes them unsuitable for representing the behavior of a biological neuron. A similar problem arises in the quaternion domain: "the only quaternion function regular with bounded norm in \(E^{4}\) is a constant" [91], where \(E^{4}\) stands for a 4-dimensional Euclidean space.

The relationship between regular functions and the existence of their quaternionic derivatives is stated by the Cauchy-Riemann-Fueter equation [92]; leading to the result that the only globally analytic quaternion functions are some linear and constant functions. Moreover, non-linear activation functions are required for constructing a neural network architecture that works as a universal interpolator of a continuous quaternion valued function [18].

In this manner, a typical approach to circumvent this problem, has been to relax the constraints, and use non-linear non-analytic quaternion functions satisfying input and output properties, where the learning dynamics is built using partial derivatives on the quaternion domain. An example of this approach are split quaternion functions, defined as follows:

\[f(\mathbf{q})=f_{R}(\mathbf{q})+f_{I}(\mathbf{q})\hat{i}+f_{J}(\mathbf{q}) \hat{j}+f_{K}(\mathbf{q})\hat{k}, \tag{64}\]

where \(\mathbf{q}\in\mathbb{H}\), \(f_{R}\), \(f_{I}\), \(f_{J}\), and \(f_{K}\) are mappings over the real numbers: \(f_{*}:\mathbb{R}\rightarrow\mathbb{R}\).

Nowadays, split quaternion functions remain the only type of activation function that has been implemented on QCNNs. Even though any type of quaternion split activation function used in QMLP can be applied, split quaternion ReLU is currently the most common activation function applied on QCNNs. For the sake of completeness, we present the activation functions found in current works.

Let \(\mathbf{q}=q_{R}+q_{I}\hat{i}+q_{J}\hat{j}+q_{K}\hat{k}\in\mathbb{H}\), we have the following activation functions:

1. Split Quaternion Sigmoid. It was introduced in [93, 46], and is defined as follows: \[\mathbb{Q}S(\mathbf{q})=S(q_{R})+S(q_{I})\hat{i}+S(q_{J})\hat{j}+S(q_{K})\hat {k},\] (65) where \(S:\mathbb{R}\rightarrow\mathbb{R}\) is the real-valued sigmoid function: \[S(x)=\frac{1}{1+e^{-x}}.\] (66)
2. Split Quaternion Hyperbolic Tangent. It was introduced in [52, 24], and is defined as follows: \[\mathbb{Q}tanh(\mathbf{q}) = \tanh(q_{R})+\tanh(q_{I})\hat{i}\] (67) \[+\tanh(q_{J})\hat{j}\] \[+\tanh(q_{K})\hat{k},\] where \(\tanh:\mathbb{R}\rightarrow\mathbb{R}\) is the real-valued hyperbolic tangent function: \[\tanh(x) = \frac{\sinh(x)}{\cosh(x)}\] (68) \[= \frac{\exp(2x)-1}{\exp(2x)+1}.\]
3. Split Quaternion Hard Hyperbolic Tangent. It was introduced in [32], and is defined as follows: \[\mathbb{Q}H^{2}T(\mathbf{q}) = H^{2}T(q_{R})+H^{2}T(q_{I})\hat{i}\] (69) \[+H^{2}T(q_{J})\hat{j}\] \[+H^{2}T(q_{K})\hat{k},\] where \(H^{2}T:\mathbb{R}\rightarrow\mathbb{R}\)[94]: \[H^{2}T(x)=\begin{cases}-1&\text{si }x<-1\\ x&\text{si }-1\leq x\leq 1\\ 1&\text{si }x>1.\end{cases}\] For this type of activation functions, Collobert [94] proved that hidden layers of MLP works as local SVM's when the real-value hard hyperbolic tangent is used.
4. Split Quaternion ReLU. It was introduced in [19, 20, 55], and is defined as follows: \[\mathbb{Q}ReLU(\mathbf{q}) = ReLU(q_{R})+ReLU(q_{I})\hat{i}\] (71) \[+ReLU(q_{J})\hat{j}\] \[+ReLU(q_{K})\hat{k},\] where \(ReLU:\mathbb{R}\rightarrow\mathbb{R}\) is the real-valued ReLU function [95, 1]: \[ReLU(x)=max(0,x).\] (72)
5. Split Quaternion Parametric ReLU. It was introduced in [22], and is defined as follows: \[\mathbb{Q}PReLU(\mathbf{q}) = PReLU(q_{R})+PReLU(q_{I})\hat{i}\] (73) \[+PReLU(q_{J})\hat{j}\] \[+PReLU(q_{K})\hat{k},\] where \(PReLU:\mathbb{R}\rightarrow\mathbb{R}\)[96]: \[PRELU(x)=\begin{cases}x&\text{si }x>0\\ \alpha x&\text{si }x\leq 0\end{cases}\] (74) and \(\alpha\) is a parameter learned during the training stage, which controls the slope of the negative side of the ReLU function. Equivalently, we have: \[PReLU(x)=max(0,x)+\alpha min(0,x).\] (75)
6. Split Quaternion Leaky ReLU. It was introduced in [97] as a particular case of the \(\mathbb{Q}\)PReLU function; in this case, \(\alpha\) is a small constant value, e.g. \(0.01\)[98].
An advantage of using split activation functions is that by processing each channel separately, we can adopt existing frameworks without additional modifications of the source code; however, this separate processing does not adequately capture the cross-dynamics of the data channels.

A different approach in the design of quaternion activation functions comes from quaternion analysis [92]. So far, it is clear that the key problem is designing suitable activation functions, and computing its quaternion derivatives; thus, some mathematicians have been working in redefining quaternion calculus.

For example, De Leo and Rotelli [99, 100], as well as Schwartz [101], have introduced the concept of local quaternionic derivative. The trick was to extend "the concept of a derivative operator from those with constant (even quaternionic) coefficients, to one with variable coefficients depending upon the point of application... [thus] the derivative operator passes from a global form to local form" [100]. Using this novel approach, they define _local analyticity_, which, in contrast to global analyticity, does not reduce functions to a trivial class (constant or linear functions). Following this idea, Ujang _et al._[102] define the concept of fully quaternion functions, and the properties they should fulfill to be locally analytic, non-linear, and suitable for gradient-based learning. They displayed better performance over split-quaternion functions when they were applied for designing adaptive filters. In the same train of thought, T. Isokawa proposed a QMLP and a backpropagation algorithm [103, 104]. Opposed to split functions, fully quaternion functions capture interchannel relationships, making them suitable for quaternion-based learning; however, experimental comparison over a standard benchmark remain an open issue.

To the best of our knowledge, currently, the only fully quaternion activation function that has been applied in QCNNs is the _Rotation-Equivariant ReLU_ function. It was proposed by Shen _et al._[23] as part of a model that extracts rotation equivariant features. Let \(\{\mathbf{q_{1}},\mathbf{q_{2}},\ldots,\mathbf{q_{N}}\}\) be a set of quaternions; then, for a quaternion \(\mathbf{q_{k}}\), with \(1\leq k\leq N\), the activation function is defined as follows:

\[\mathbb{Q}REReLU(\mathbf{q_{k}})=\frac{\|\mathbf{q_{k}}\|}{max(\|\mathbf{q_{ k}}\|,c)}\mathbf{q_{k}}; \tag{76}\]

where \(c\) is a positive constant, computed as follows:

\[c=\frac{1}{N}\sum_{j=1}^{N}\|\mathbf{q_{j}}\|. \tag{77}\]

In the same trend of developing novel quaternion calculus tools, but from a different perspective, Mandic _et al._[105] start from the observation that for gradient-based optimization, a common objective is to minimize a positive real function of quaternion variables, e.g. \(J(\mathbf{e},\mathbf{\bar{e}})=\mathbf{e\bar{e}}\), and for that purpose it is used the pseudo-gradient, i.e. the sum of component-wise gradients. Formalization of these ideas is achieved by "establishing the duality between derivatives of quaternion valued functions in \(\mathbb{H}\) and the corresponding quadrivariate real functions in \(\mathbb{R}^{4}\)" [105], leading to what is called Hamiltonian-Real Calculus (HR). In addition, Mandic _et al._[105] proved that for a real function of quaternion vector variable, the maximum change is in the direction of the conjugate gradient, establishing a general framework for quaternion gradient-based optimization. Going further, Xu and Mandic [106] proposed the product and chain rules for computing derivatives, as well as the quaternion counterparts of the mean value and Taylor's theorems; this establishes an alternative framework called Generalize HR calculus (GHR). Thereafter, novel quaternion gradient algorithms using GHR calculus were proposed [107, 108]. Although in [109] is shown that a QMLP trained with a GHR-based algorithm obtains better prediction gains, on the 4D Saito's chaotic signal task, than other quaternion-based learning algorithms [18, 110, 57], further experimental analysis on a standard benchmark and proper comparison with real and complex counterparts is required. To this date, there is no published work on the use of QRH calculus on QCNNs.

_Future directions:_ So far, we have identified the fundamental trends of thought for activation functions: split-quaternion functions, whose derivatives for training are computed in a channel wise manner, and fully quaternion functions, whose derivatives are computed locally, or using partial derivatives. Theoretical analysis, as well as some preliminary experimental results, have indicated a better performance of fully quaternion activation functions over others, as well as better performance of quaternion training methods based on GHR calculus. However, these ideas have not been set in practice on QCNNs. Intended works should focus on introducing novel fully connected activation functions that exploit specific properties of the quaternion representation, and establishing proper benchmarks for comparison of existing functions and training methods. In respect to quaternion calculus, the following section will connect the preceding ideas with the training of QCNNs and its future work.

### _Training_

Currently, all of the methods for training that have been tested on QCNNs, rely on adaptations of the QMLP back-propagation algorithm [48], or an extension of the generalized complex chain rule for real-valued loss functions [49]. Both approaches are equivalent and relax the analyticity condition by using partial derivatives with respect to the real and imaginary parts. Alternative approaches, such as GHR calculus or local derivatives have not been tested in current implementations of QCNNs.

Thus, Gaudet and Maida [19] introduced the Generalized Quaternion Chain Rule for a Real-valued function: Let \(L\) be a real-valued loss function and \(\mathbf{q}=q_{R}+q_{I}\hat{i}+q_{J}\hat{j}+q_{K}\hat{k}\), where \(q_{R},q_{I},q_{J},q_{K}\in\mathbb{R}\), then [19, 92]:

\[\frac{\partial L}{\partial\mathbf{q}}=\frac{\partial L}{\partial q_{R}}+\frac {\partial L}{\partial q_{I}}\hat{i}+\frac{\partial L}{\partial q_{J}}\hat{j}+ \frac{\partial L}{\partial q_{K}}\hat{k}. \tag{78}\]

Now, assuming \(\mathbf{q}\) can be expressed in terms of a second quaternion variable, \(\mathbf{g}=g_{R}+g_{I}\hat{i}+g_{J}\hat{j}+g_{K}\hat{k}\), where \(g_{R},g_{I},g_{J},g_{K}\in\mathbb{R}\). Then, the chain rule is calculated as
follows [19]:

\[\frac{\partial L}{\partial\mathbf{g}} = \frac{\partial L}{\partial q_{R}}\left(\frac{\partial q_{R}}{\partial g _{R}}+\frac{\partial q_{R}}{\partial g_{I}}\hat{i}+\frac{\partial q_{R}}{ \partial g_{J}}\hat{j}+\frac{\partial q_{R}}{\partial g_{K}}\hat{k}\right)+ \tag{79}\] \[\frac{\partial L}{\partial q_{I}}\left(\frac{\partial q_{I}}{ \partial g_{R}}+\frac{\partial q_{I}}{\partial g_{I}}\hat{i}+\frac{\partial q_ {I}}{\partial g_{J}}\hat{j}+\frac{\partial q_{I}}{\partial g_{K}}\hat{k}\right)+\] \[\frac{\partial L}{\partial q_{J}}\left(\frac{\partial q_{J}}{ \partial g_{R}}+\frac{\partial q_{J}}{\partial g_{I}}\hat{i}+\frac{\partial q _{J}}{\partial g_{J}}\hat{j}+\frac{\partial q_{J}}{\partial g_{K}}\hat{k}\right)+\] \[\frac{\partial L}{\partial q_{K}}\left(\frac{\partial q_{K}}{ \partial g_{R}}+\frac{\partial q_{K}}{\partial g_{I}}\hat{i}+\frac{\partial q _{K}}{\partial g_{J}}\hat{j}+\frac{\partial q_{K}}{\partial g_{K}}\hat{k}\right)\]

These equations are applied in the implementation of the backpropagation algorithm. Since this is the most used approach, we call it the _Standard Quaternion Backpropagation Algorithm_[46, 93, 111], summarized as follows:

Let \(\mathbf{x}\in\mathbb{H}\) be the input to a layer \(C\), \(\mathbf{w_{nm}}\in\mathbb{H}\) represents the weight connecting input \(n\) to output \(m\), and \(\mathbf{d}^{top},\mathbf{d}^{bottom}\in\mathbb{H}\) are the error propagated from the top and to the bottom layers, respective, and \(\epsilon\in\mathbb{R}\) be the learning rate; then, for the current convolution layer, \(C\):

1. Update its weights using the following equation: \[\mathbf{w_{nm}}=\mathbf{w_{nm}}+\epsilon\mathbf{d_{m}^{top}}\mathbf{\bar{x}_ {n}}\] (80)
2. Update the bias term: \[\mathbf{b_{m}}=\mathbf{b_{m}}+\epsilon\mathbf{d_{m}^{top}}\] (81)
3. Propagate the error to the bottom layer according to the following equation: \[\mathbf{d_{n}^{bottom}}=\sum_{m}(\mathbf{\bar{w}_{nm}}\mathbf{d_{m}^{top}})\] (82)

Note that Equations (80) and (82) use quaternion products. For an activation layer, the error is propagated to the bottom layer according to the following equation:

\[\mathbf{d_{n}^{bottom}}=\mathbf{d_{n}}^{top}\odot f^{\prime}(\mathbf{x_{n}}) \tag{83}\]

where \(\odot\) is the Hadamard or Schur product (component wise).

Adopting the same definition of the chain rule, but a different definition of convolution, we have the work of Zhu _et al._[24], whom apply the two-sided convolution and a polar representation of the quaternion weights. Since their model applies rotation and scaling on the input features, the quaternion gradient for his model is simplified to a rotation transformation over the same axis, but with a reversed angle. A general model of this approach, with learnable arbitrary axes, is presented by [61]. Therefore, the backpropagation algorithm is similar to the one presented before, but Equation (80) changes accordingly [58, 57, 61]:

\[\frac{\epsilon}{\|\mathbf{w_{nm}}\|}\left\{\frac{\mathbf{d_{m}^{top}}\cdot( \mathbf{w_{nm}}\mathbf{x_{n}}\mathbf{\bar{w}_{nm}})}{\|\mathbf{w_{nm}}\|^{2}} \mathbf{w_{nm}}-2\mathbf{d_{m}^{top}}\mathbf{w_{nm}}\mathbf{\bar{x}_{n}}\right\} \tag{84}\]

and Equation (82) changes to [58, 57, 61]:

\[\mathbf{d_{n}^{bottom}}=\sum_{m}\frac{\mathbf{\bar{w}_{nm}}\mathbf{d_{m}^{top}} \mathbf{w_{nm}}}{\|\mathbf{w_{nm}}\|} \tag{85}\]

From the point of view of supervised learning, the problem of training a QCNN can be stated as an optimization one:

\[\underset{\mathbf{w_{1}},\dots,\mathbf{w_{n}}}{\mathrm{argmin}}\ L(\mathbf{w_{1 }},\dots,\mathbf{w_{n}}) \tag{86}\]

where \(L\) is the loss function, and \(\mathbf{w_{1}},\dots,\mathbf{w_{n}}\) are all the quaternion weights of the network.

Besides the problem of computing derivatives of non-analytic quaternion activation functions; another problem is that algorithms such as gradient descendant and back propagation could be trapped in local minima. Therefore, QCNNs could be trained with alternative methods that not rely on the computing of quaternion derivatives, such as evolutionary algorithms, ant colony optimization, particle swarm optimization, etc. [112]. In this case, general-purpose optimization algorithms will have the same performance in average [113, 114].

_Future directions:_ Novel methods for training should be applied, for example: modern quaternion calculus techniques, such as local or GHR calculus, and meta-heuristic algorithms working in the quaternion domain. In addition, fully quaternion loss functions should be introduced.

### _Quaternion weight initialization_

At the beginning of this century, real-valued neural networks were showing the superiority of deep architectures; however, the standard gradient descent algorithm from random weight initialization performed poorly when used for training these models. To understand the reason of this behavior, Glorot and Bengio [115] established an experimental setup and observed that some activation functions can cause saturation in top hidden layers. In addition, by theoretically analyzing the forward and backward propagation variances (expressed with respect to the input, output and weight initialization randomness), they realized that, because of the multiplicative effect through layers, the variance of the back-propagated gradient might vanish or explode in very deep networks. Consequently, they proposed, and validated experimentally, a weight initialization procedure, which makes the variance dependent on each layer, and maintains activation and back-propagated gradient variances as we move up or down the network. Their method, called _normalized initialization_, uses a scaled uniform distribution for initialization. However, one of its assumptions is that the network is in a linear regime at the initialization step; thus, this method works better for softsign units than for sigmoid or hyperbolic tangent ones [115]. Since the linearity assumption is invalid for activation functions such as ReLU, He _et al._[96] developed an initialization method for non-linear ReLU and Parametric ReLU activation functions. Their initialization method uses a zero-mean Gaussian distribution with a specific standard deviation value, and surpasses the performance of the previous method for training extremely deep models, e.g. 30 convolution layers.

In the case of classic QCNNs, Trabelsi _et al._[49] extended these results to deep complex networks. Similar results are presented by Gaudet and Maida [19] for deep quaternion networks, and by Parcollet _et al._[52] for quaternion recurrent networks. These works treat quaternion weights as 4
-dimensional vectors, whose components are normally distributed, centered at zero, and independent. Hence, it can be proved that the weights and their magnitude follow a 4DOF Rayleigh distribution [19, 52], reducing the weight initialization by selecting a single parameter, \(\sigma\), which indicate the mode of the distribution. If \(\sigma=\frac{1}{\sqrt{2(n_{in}+n_{out})}}\) we have a _quaternion normalized initialization_ which ensures that the variances of the quaternion input, output and their gradients are the same; while \(\sigma=\frac{1}{\sqrt{2n_{in}}}\) is used for the _quaternion ReLU initialization_, where \(n_{in}\), and \(n_{out}\) are the number of neurons of the input and output layers, respectively. The method is presented in Algorithmic 1.

On the other hand, for geometric QCNNs, the quaternion weights represent an affine transformation composed by a scale factor, \(s\), and a rotation angle \(\theta\). Thus, to keep the same variance of the gradients during training, Zhu _et al._[24] proposed a simple initialization procedure using the uniform distribution, \(U[\cdot]\):

\[s_{j} \sim U\left[-\frac{\sqrt{6}}{\sqrt{n_{j}+n_{j+1}}},\frac{\sqrt{6}}{ \sqrt{n_{j}+n_{j+1}}}\right]\] \[\theta \sim U\left[-\frac{\pi}{2},\frac{\pi}{2}\right] \tag{87}\]

_Future directions:_ In the case of geometric QCNNs, further theoretical analysis of weight initialization techniques is required, as well as the extension for ReLU units. For equivariant QCNNs, weight initialization procedures have not been introduced or or are not described in the literature. In addition, current works only consider the case of split quaternion activation functions. The propagation of variance using fully quaternion activation functions should be investigated.

## V Architectural design for applications.

The previous section deals with the individual building blocks of QCNNs; the possible ways we can interconnect them, give rise to numerous models. In this manner, there are three leading factors to consider when implementing applications:

1. Domain of application. Current works are primarily focused on 3 areas: vision, language, and forecasting.
2. Mapping the input data from real numbers to the quaternion domain.
3. Topology. Based on current works, we have the following types: * ConvNets. They use convolution layers without additional tricks. * Residual. They use a shortcut connection from input to forward blocks, allowing to learn with reference to the layer inputs instead of learning unreferenced functions. * Convolution Auto-Encoders (CAE). Its aim is to reconstruct the input feature at the output. * Point-based. They focused on unordered sets of vectors such as point-cloud input data. * Recurrent. They exploit connections that can create cycles, allowing the network to learn long-term temporal information. * Generative. They use an adversarial learning paradigm, where one player generates samples resembling the real data, and the other discriminates between real and fake data; the solution of the game is the model achieving Nash equilibrium.

This section is devoted to show how the blocks presented in Section IV can be used to construct different architectures, and we focus our analysis on the three points previously introduced.

The information is organized in three subsections corresponding to the domains of applications: vision, language, and forecasting. Within each subsection, the models are ordered according to the type of convolution they applied: classic, geometric, or equivariant. In addition, different networks are presented according to their topology: ConvNet, Residual, CAE, Point-based, Recurrent, or Generative. Methods for mapping input data to the quaternion domain are presented in each subsection.

A graphic depiction of several models is shown in Figures 5, 6, and 7. Some architectures available in the literature are presented from a high-level abstraction point of view; showing their topologies and the blocks they use. Information about the number of quaternion kernels, size, and stride was included when available. Further implementation details, such as weight initialization, optimization methods, etc. were omitted to avoid a cumbersome presentation. Readers interested in specific models are encouraged to consult the original sources.

In addition, Tables II to V provide a comparison between different QCNNs. They provide extra information, such as datasets on which the models were tested, performance metrics, comparison to real-valued models, etc. For comparison between quaternion and real-valued models, we only report real-valued architectures that have similar topology. The difference in the number of parameters between real-valued and quaternion-valued networks is due to the following reasons: Some authors prefer to compare networks with the same number of real-valued units, and quaternion-valued units, where the quaternion-valued units have \(4\) times more parameters. Others compare networks with the same number of parameters: they reduce the number of quaternion units to \(25\%\) of the real
valued networks, or quadruple the number of real-valued units. In addition, some quaternion-value networks use real-valued components in some parts of the network.

To conclude, Subsection V-D summarizes some insights obtained from this works.

### _Vision_

For mapping data to the quaternion domain, several methods are available in the literature; in the case of color images, the most common approach is to encode the red, green, and blue channels into the imaginary parts of the quaternion:

\[\mathbf{q}=0+R\hat{i}+G\hat{j}+B\hat{k}, \tag{88}\]

see for example [127, 129, 97, 82, 132, 20, 122, 74, 32, 116, 55, 24].

Another method, proposed by Gaudet and Maida [19], is to use a residual block:

\(BN\to ReLU\to Conv\to BN\to ReLU\to Conv\), where a shortcut connection sums the input to the output of the block.

In contrast, for grayscale images, Beijing _et al._[129] propose to map the grayscale values to the real part of the quaternion as follows:

\[\mathbf{q}=grayscale+0\hat{i}+0\hat{j}+0\hat{k} \tag{89}\]

In the case of POLSAR images containing scattering matrices:

\[\begin{bmatrix}S_{HH}&S_{HV}\\ S_{VH}&S_{VV}\end{bmatrix} \tag{90}\]

Matsumoto _et al._[61] propose to compute a Pauli decomposition from the complex scattering matrix:

\[a=\frac{S_{HH}+S_{VV}}{\sqrt{2}},b=\frac{S_{HH}-S_{VV}}{\sqrt{2}},c=\frac{S_{ Hv}+S_{Vh}}{\sqrt{2}}, \tag{91}\]

and assign the square magnitude of the components \(b\), \(c\), and \(a\) to the red, green, and blue channels, respectively, of an RGB image. Thereafter, the image is mapped to the quaternion domain as usual. Alternatively, Matsumoto _et al._[61] propose to transform the scattering matrix into 3D Stokes vectors normalized by their total power. Thereafter, each component is mapped to the imaginary parts of the quaternion.

Finally, for point clouds, Shen _et al._[23] propose to map the \((x,y,z)\) coordinates of each 3D point to a quaternion as follows:

\[\mathbf{q}=0+x\hat{i}+y\hat{j}+z\hat{k}. \tag{92}\]

#### V-A1 Vision under the classic convolution paradigm

Residual Networks were one of the first models to be proposed in the QCNNs literature. Gaudet and Maida [19] tested deep and shallow quaternion-valued residual nets. For both cases, they use three stages; the shallow network contains 2, 1 and 1 residual blocks in each stage, while the deep network uses 10, 9 and 9 residual blocks in each stage, see Figure 5a). These models were applied on image classification tasks using the CIFAR-10 and CIFAR-100 [50] datasets, and the KITTI Road Estimation benchmark [51]. Recently, Sfikas _et al._[119] proposed a standard residual model for keyword spotting in handwritten manuscripts, see Figure 5c). Their model is applied on the PIOP-DAS [120] and the GRPOLY-DB [121] datasets, containing digitized documents written in modern Greek.

A different type of model is ConvNets. The first proposal came from Yin _et al._[55], whom implemented a basic QCNN model, see Figure 6f). Since quaternion models extract \(4\) times more features, they propose to use an attention module, whose purpose is to filter out redundant features. These models were applied to the image classification problem on the CIFAR-10 dataset [50]. In addition, they propose a similar model for Double JPEG compression detection on the UCID dataset [56]. Thereafter, Jin _et al._[122] propose to include a deformable layer, leading to a Deformable Quaternion Gabor CNN (DQG-CNN), see Figure 6h). They apply it on a facial recognition task using the Oulu-CASIA [123], MMI [124], and SFEW [125] datasets. Beijing _et al._[127] propose an architecture based on quaternion versions of Fully Convolutional Neural Networks (FCN) [139]. Their model is applied for the color image splicing localization problem, and tested on CASIA v1.0, CASIA v2.0 [128], and Columbia color DVMM [130] datasets. The basis of their architecture is the Quaternion-valued Fully Convolutional Network (QFCN), see Figure 6d). Then, the final model is composed of three QFCNs working in parallel, each one has different up-sampling layers: The first network does not have extra connections, the second network has a shortcut connection fusing the results of the fifth pooling layer with the output of the last layer, while the third network combines results of the third, and fourth pooling layers with the output of the last layer. In addition, each network uses a Super-pixel-enhanced pairwise Condition Random Field module to improve the results from the QCNN. This fusion of the outputs allows to work with different scales of image contents [140]. Thereafter, Beijing _et al._[129] propose the quaternion-valued two-stream Region-CNN. They extend the real-valued RGB-N [141] to the quaternion domain, improve it for pixel-level processing, and implement two extra-modules: an Attention Region Proposal Network, based on CBAM [142], for enhancing the features of important regions; and a Feature Pyramid Network, based on Quaternion-valued ResNet, to extract multi-scale features. Their results improved the ones previously published on [127]. Because of the complexity of these models, it is recommended to read it directly from [127, 129].

In the case of CAE's, which aim is to reconstruct the input feature at the output, Parcollet _et al._[32] propose a quaternion convolutional encoder-decoder (QCAE), see Figure 7a), and tested on the KODAK PhotoCD dataset [54].

Another type of CAE's is Variational Autoencoder, which estimates the probabilistic relationship between input and latent spaces. The quaternion-valued version was proposed by Grassucci _et al._[97], see Figure 7b), and applied it to reconstruct and generate faces of the CelebFaces Attributes Dataset (CelebA) [126].

A more sophisticated type of architecture is the generative models. Sfikas _et al._[116] propose a Quaternion Generative Adversarial Network for text detection of inscriptions found on byzantine monuments [118, 117]. The generator is a
U-Net [143] like model, and the discriminator a cascade of convolution layers, see Figures 7e) and 7f), where the activation function of the last layer sums the output of real and imaginary parts of the quaternion to produce a real-valued output.

Grassucci _et al._[82] adapted a Spectral Normalized GAN's (SNGAN) [144, 145] to the quaternion domain, and apply it on an image to image translation task using the CelebA-HQ [131] and 102 Oxford Flowers [59] datasets. The model is presented in Figure 7c) and 7d). Thereafter, Grassucci _et al._[132] proposed the quaternion-valued version of the StarGANv2 model [146]. It is composed of the generator, mapping, encoding, and discriminator networks; this model was evaluated on an image to image translation task using the CelebA-HQ dataset [131]. Because of the complexity of the model, it is recommended to consult them directly from [132].

#### Iv-A2 Vision under the geometric paradigm

The only residual model lying in this paradigm is the one of Hongo _et al._[20], whom proposed a residual model, based on ResNet34 [147], see Figure 5a). This is applied to the image classification task using the CIFAR-10 dataset [50].

On ConvNet models, we have the work of Zhu _et al._[24], whom proposed shallow QCNN and quaternion VGG-S [148] models for image classification problems. These models are shown in Figures 6b) and 6c); the former was tested on the CIFAR-10 dataset [50], and the latter on the 102 Oxford flower dataset [59]. Hongo _et al._[20] proposed a different QCNN model, see Figure 6g), and tested on the
image classification task using the CIFAR-10 dataset [50]. Moreover, Matsumoto _et al._[61] propose QCNN models for classifying pixels of PolSAR images. This type of images contain additional experimental features given in complex scattering matrices. For their experiments, they labeled each pixel of two images in one of 4 classes: water, grass, forest, or town. Then, they converted the complex scattering matrices into PolSAR pseudocolor features, or into normalized Stoke vectors. By testing similar models under these two different representations, they found out that the classification results largely depend on input features. The proposed model is shown in Figure 6j).

In the context of CAE's, Zhu _et al._[24] propose a U-Net-like encoder-decoder network [143] for the color image denoising problem. The model was tested for a denoising task on images of the 102 Oxford flower dataset [59], and on a subset of the COCO dataset [60], see details of the model in [149, 143, 150, 24].

#### Iv-A3 Vision under the equivariant paradigm

Inspired by the PointNet model [151] for processing point clouds, Shen _et
al._[23] modified it by exchanging all its layers for Rotation Equivariant Quaternion Modules, and remove its Spatial Transformer module since it discards rotation information. The rotation equivariant properties of the modules are evaluated on the ShapeNet dataset [152]. They experimentally proved that point clouds reconstructed using the synthesized quaternion features had the same orientations as point clouds generated by directly rotating the original point cloud. In addition, they modify other models, i.e. PointNet++ [153], DGCNN [154], and PointConv [155], by replacing their components into its equivariance counterparts, and tested on a 3D shape classification task on the ModelNet40 [62] and 3D MNIST [63] datasets. Because the variety of components and interconnections, the reader is directed to [151, 153, 23, 154, 155] to consult the details of these models.

### _Language_

For language applications, all works that have been proposed are ConvNets.

Under the classic paradigm, Parcollet _et al._[22] propose a Connectionist Temporal Classification CTC-QCNN model, see Figure 6a), and tested on a phoneme recognition task, with TIMIT dataset [53]. The mapping of input signals to the quaternion domain is achieved by transforming the raw audio into 40-dimensional log mel-filter-bank coefficients with deltas, delta-deltas, and energy terms, and arranging the resulting vector into the components of an acoustic quaternion signal:

\[\mathbf{q}(\mathbf{f},\mathbf{t})=0+e(f,t)\hat{i}+\frac{\partial e(f,t)}{ \partial t}\hat{j}+\frac{\partial^{2}e(f,t)}{\partial t^{2}}\hat{k}. \tag{93}\]

Another problem related to language is joint 3D Sound Event Localization and Detection (SELD); solving this task can be helpful for activity recognition, assisting hearing impaired people, among other applications. The SELD problem consists in simultaneous solving: the Sound Event Detection (SED) problem, i.e. in a set of overlapping sound events, detecting temporal activities of each sound event and associating a label, and the Sound Localization Problem, which consists of estimating the spatial localization trajectory of each sound. The latter task could be simplified to determining the orientation of a sound source with respect to the microphone, which is called Direction-of-Arrival (DOA). For this problem, also following a classical approach, Comminiello _et al._[133] propose a quaternion-valued recurrent network (QSELD-net), based on the real-valued model of Adavanne _et al._[156]. The model has a first processing stage, where output is processed in parallel by two branches: the first one performs a multi-label classification task (SED), while the second one performs a multi-output regression task (DOA estimation). The model is shown in Figure 6e), where \(N\) is the number of sound event classes to be detected; for DOA estimation, we have three times more outputs, since they represent \((x,y,z)\) coordinates for each sound event class. This model was tested on the Ambisonic, Anechoic and Synthetic Impulse Response (ANSYN), and the Ambisonic, Reverberant and Synthetic Impulse Response (RESYN) datasets [134], consisting of spatially located sound events in an anechoic/reverberant environment synthesized using artificial impulse responses. Each dataset comprises three subsets: no temporally overlapping sources (O1), maximum two temporally overlapping sources (O2) and maximum three temporally overlapping sources (O3). In this application, the input data is a multichannel audio signal; real-valued SELD-net [156] as well as the quaternion-valued counterpart [133] apply the same feature extraction method: the spectrogram is computed using the \(M\)-point discrete Fourier transform
to obtained a feature sequence of \(T\) frames containing the magnitude and phase components for each channel.

Under the geometric paradigm, Muppidi and Radfar [74] proposed a QCNN for the emotion classification from speech task and evaluate it on the RAVDESS [135], IEMOCAP [136], and EMO-DB [137] datasets. They convert speech waveform inputs from the time domain to the frequency domain using Fourier transform, thereafter compute its Mel-Spectrogram,

Fig. 5: Residual QNN models. Some notation details: For convolution, parameters are presented in the format (number of kernels, size of kernel, stride); if the size of kernel or stride is the same for each dimension, it is shown just a single number, e.g. instead of \(3\times 3\) it is shown \(3\). Symbol \(+\) inside a circle means summing. Dropout and flatten procedures were omitted. All blocks represent operations in the quaternion domain.
and convert it to an RGB image. Finally, RGB images are processed with the model shown in Figure 6i), where the Reset-ReLU activation function resets invalid values to the nearest point in color space.

### _Forecasting_

Neshat _et al._[138] proposed an hybrid forecasting model, composed of QCNNs and Bi-directional LSTM recurrent networks, for prediction of short-term and long-term wind

Fig. 6: Quaternion ConvNet models. Some notation details: For convolution, parameters are presented in the format (number of kernels, size of kernel, stride); if the size of kernel or stride is the same for each dimension, it is shown just a single number, e.g. instead of \(3\times 3\) it is shown \(3\). LRN stands for Local Response Normalization. Dropout and flatten procedures were omitted. All blocks represent operations in the quaternion domain, except those preceded by the word “real”.
speeds. Historical meteorological wind data was collected from Lesvos and Samothraki Greek islands located in the North Aegean Sea, and obtained from the Institute of Environmental Research and Sustainable Development (IERSD), and the National Observatory of Athens (NOA). This work integrates classic standard QCNNs within a complete system, and shows that better results are obtained when using QCNNs than other AI models or handcrafted techniques. Because of the lack of data, the QCNNs are not able to provide a better understanding of the data, and the results are not available.

Fig. 7: Quaternion CAE and Generative models. Some notation details: For convolution, parameters are presented in the format (number of kernels, size of kernel, stride); if the size of kernel or stride is the same for each dimension, it is shown just a single number, e.g. instead of \(3\times 3\) it is shown \(3\). Symbol \(+\) inside a circle means summing while \(\cup\) inside a circle means concatenation. TConv stands for Decoding Transposed Convolution. Dropout and flatten procedures were omitted. All blocks represent operations in the quaternion domain, except those preceded by the word “real”.
the complexity of the model, it is recommended to consult the details directly from [138].

### _Insights_

A fair comparison between the performance of the models and their individual blocks is difficult because of the fact that most of them are applied on different problems or using different datasets; moreover, most of the quaternion-valued models are based on their real-valued counterparts, and there are very few works based on incremental improvement over previous quaternion models. However, from the works reviewed, the following insights can be established:

1. Quaternion-valued models achieve, at least, comparable results to their real-valued counterparts, and with a lower number of parameters. As was stated by [119]: "\(\ldots\)an operation that is written as \(y=Wx\), where \(y\in\mathbb{R}^{4K}\), \(x\in\mathbb{R}^{4L}\), and \(W\in\mathbb{R}^{4K\times 4L}\) is thus mapped to a quaternionic operation \(\mathbf{y}=\mathbf{W}\mathbf{x}\), where \(\mathbf{y}\in\mathbb{H}^{K}\), \(\mathbf{x}\in\mathbb{H}^{L}\) and \(\mathbf{W}\in H^{K\times L}\). Parameter vector \(\mathbf{W}\) only contains \(4\times K\times L\) parameters compared to \(4\times 4\times K\times L\) of \(W\), hence we have a \(4\)x saving." In addition, some works report faster convergence on the training stage.
2. For the image classification task, the CIFAR-10 dataset has become the standard benchmark. Although the classification errors reported in different works are not standardized, for example: some authors use data-augmentation techniques, the best reported result is from Gaudet and Maida [19], whom uses a classic approach and a deep residual model. In addition, residual models present the best results for image classification tasks.
3. Attention modules filter out redundant features, and their use can improve the performance of the network [55].
4. Variance Quaternion Batch Normalization (VQBN) experimentally outperforms Split Batch Normalization [55].
5. The integration of quaternion Gabor filters and convolution layers enhances the abilities of the model to capture features in the aspects of spatial localization, orientation selectivity, and spatial frequency selectivity [122, 157, 158].
6. Quaternion generative models demonstrate better abilities to learn the properties of the color space, and generate RGB images with more defined objects than real-valued models [82, 32].

Finally, Table VI presents information about the available source code. It only includes source code published by the authors of the papers; other implementations can be found in the popular website _paperswithcode_[159].

## VI Further discussion

Section IV presented each individual block, discussion of current issues, knowledge gaps, and future works for improvement. Here, we discuss future directions of research and open questions for any model of QCNNs, independently of its components. Even more, some topics involve any type of quaternion-valued deep learning model.

### _On proper comparisons_

When looking at the different models that have been proposed, one of the questions that arises is: Is it fair to compare a real-valued architecture with a quaternion-valued one using the same topology and number of parameters? Note that in the first case, the optimization of the parameters occurs in an Euclidean \(4\)-dimensional space, while in the second case, the optimization occurs in the quaternion space, which can be connected to geometric spaces different to the Euclidean one. From this observation, we could argue that a more reasonable way to compare is: an optimal real value network vs the optimal quaternion-value network, even though they have different topologies, connections, and number of weights. Now, it is clear that a single real-valued convolution layer cannot capture interchannel relationships, but could a real-valued multilayer or recurrent architecture capture interchannel relationships without the use of Hamilton products? When are Hamilton products _really_ needed?

Something that could shed light on these questions is the reflection of Sfikas _et al._[119], whom believe that the effectiveness of quaternion networks is because of "navigation on a much compact parameter space during learning" [119], as well as an extra parameter sharing trait.

### _Mapping from data to quaternions_

Another question that arises is: how to find an optimal mapping from the input data to the quaternion domain? Common sense would say that finding a suitable mapping, and connecting with the right geometry could seriously improve the performance of quaternion networks; however, theoretical and experimental work is required in this direction. For example, for images, we can adopt different color models with a direct geometric relation to quaternion space, and point clouds as well as deformable models can avoid algebraic singularities when the quaternion representation is used. Moreover, models for language tasks are in their infancy, and novel methods that take advantage of the quaternion space should be proposed for signal processing, speech recognition, and other tasks.

### _Extension to hypercomplex and geometric algebras_

Current QCNNs rely on processing 4-dimensional input data, which is mapped to the quaternion domain. If the input data have less than \(4\) channels, a common solution is to apply zero-padding on some channels, or define a mapping to a \(4\)-dimensional space.

In contrast, for more than \(4\) channels, there are several alternatives:

The first one is to map the input data to an \(n\)-dimensional space, where \(n\) is a multiple of four. Then, we can process the input data by defining a quaternion kernel for each 4-channel input, see Figure 4. Alternatively, it can be applied an extension of the quaternion convolution to the octonions, sedenions or generalized hypercomplex numbers [67]. Moreover, Geometric algebras [69, 160] can be used to generalize hypercomplex convolution to general \(n\)-dimensional spaces, not restricted to multiples of \(4\), and to connect with different geometries.
A partial approach in this direction is to parametrize the hypercomplex convolution [70, 71]. Let \(W\) be a convolution kernel, of size \(k\times k\), \(x\) a \(s\)-dimensional input, and \(y\) a \(d\)-dimensional output, then:

\[y=W*x. \tag{94}\]

Thus, the kernel is decomposed into the sum of its Kronecker products:

\[y=\sum_{i=1}^{n}A_{i}\otimes F_{i}*x. \tag{95}\]

where \(A_{i}\in\mathbb{R}^{n\times n}\), with \(i=1,\ldots,n\) are the matrices containing the algebra multiplication rules, and \(F_{i}\in\mathbb{R}^{\frac{\epsilon}{n}\times\frac{\epsilon}{n}\times k\times k}\) are the filters that compose the final weight matrix. The parameter \(n\) defines the dimension of the algebra; then, for \(n=2\) we are working in the complex domain, while a value \(n=4\) leads us to the quaternion space. Matrices \(A_{i}\), and \(F_{i}\) are obtained during training.

### _Quaternions on the frequency domain_

A natural implementation of QCNNs in the frequency domain would be to compute the Quaternion Fourier transform [161, 25, 162] of the input and kernels, and multiply them in the frequency domain. Moreover, the computing of the convolution could be accelerated using Fast QFT algorithms [163]. To this date, works following this approach have not been published. However, using a bio-inspired approach, Moya-sanchez _et al._[158] proposed a monogenic convolution layer. A monogenic signal [164], \(I_{M}\), is a mapping of the input signal, which simultaneously encodes local space and frequency characteristics:

\[I_{M}=I^{\prime}+I_{1}\hat{i}+I_{2}\hat{j} \tag{96}\]

where \(I_{1}\) and \(I_{2}\) are the Riesz transforms of the input in the \(x\) and \(y\) directions. In [158], the monogenic signal is computed in the Fourier domain, \(I^{\prime}\) is computed as a convolution of the input with quaternion log-Gabor filters with learnable parameters, and local phase and orientation are computed for achieving contrast invariance. In addition, their model provides crucial sensitivity to orientations, resembling properties of V1 cortex layer.

### _Might the classic, geometric, and equivariant models be special cases of a unified general model?_

To answer this question, let us recall that a quaternion is an element of a 4-dimensional space, and unitary versors together with the Hamilton product are isomorphic to the group SO(4). Thus, it is our point of view that the classic model, using the four components, is the most general one. By dividing the product as a sandwiching product and using an equivalent polar representation, a generalized version of the geometric model could be obtained (not Euclidean or affine, but 3D projective model); in addition, for the equivariant model, the kernel is reduced merely to its real components. However, the connection to invariant theory should be investigated for each of these models, so that a stratified organization at the light of geometry might be achieved. For example, if we link a unified general model to a particular geometry and its invariant properties, we could go down on the hierarchy of geometries by setting restrictions on the general model, until we obtain the rotation equivariant model or the others.

## VII Conclusions

We have recollected the substantial majority of the ideas that have been published on QCNNs in recent years, and we presented a comprehensive guide for its application. In particular, being the convolution layer the core component of this type of models, we presented a sounded organization of QCNNs based on the definition of convolution; therefore, we proposed three classifications: classic, geometric, and rotation equivariant. For other components, we presented the purpose of each block, the key problem to solve in its design, knowledge gaps, if any, and ideas for future improvement. In addition, a review of the models by application domain and topology type was presented, including available source code. Further ideas for model design were also discussed.

Finally, most of the ideas that have been implemented are extensions of the work on real-valued CNNs or complex CNNs to the quaternion domain [165, 49], or adaptations from quaternion neural networks models [166]. Further work is required in developing novel ideas, and exploiting the particularities of quaternion representation and its connection to geometry, topology, functional analysis, or invariant theory. This paper presents, in an organized manner, the current advances in the development of QCNNs, in the hope that it serves as a starting point for subsequent research, as well as for those interested in implementing applications.

## Acknowledgments

G.A.G. received a Postdoctoral Fellowship _Estancias Posdoctorales por Mexico_ from CONACYT. C.G. acknowledges support from UNAM-PAPIIT (IN107919, IV100120, IN105122) and from the PASPA program from UNAM-DGAPA.