# The $\mathbb{HR}$-Calculus: Enabling Information Processing with Quaternion
Algebra

Danilo P. Mandic, Sayed Pouria Talebi, Clive Cheong Took,
Yili Xia, Dongpo Xu, Min Xiang, and Pauline Bourigault

## Abstract

From their inception, quaternions and their division algebra have proven to be
advantageous in modelling rotation/orientation in three-dimensional spaces and
have seen use from the initial formulation of electromagnetic filed theory
through to forming the basis of quantum filed theory. Despite their impressive
versatility in modelling real-world phenomena, adaptive information processing
techniques specifically designed for quaternion-valued signals have only
recently come to the attention of the machine learning, signal processing, and
control communities. The most important development in this direction is
introduction of the $\mathbb{HR}$-calculus, which provides the required
mathematical foundation for deriving adaptive information processing
techniques directly in the quaternion domain. In this article, the foundations
of the $\mathbb{HR}$-calculus are revised and the required tools for deriving
adaptive learning techniques suitable for dealing with quaternion-valued
signals, such as the gradient operator, chain and product derivative rules,
and Taylor series expansion are presented. This serves to establish the most
important applications of adaptive information processing in the quaternion
domain for both single-node and multi-node formulations. The article is
supported by Supplementary Material, which will be referred to as SM.

##  I Introduction

For the best part of ten year Sir William Rowan Hamilton concerned himself
with constructing a framework for modelling three-dimensional spaces in a
manner analogous to that achieved through complex-valued algebra for two-
dimensional spaces. This came to fruition on the $16^{\text{th}}$ of October
$1843$, when he derived the fundamental elements of what is known today as the
quaternion algebra and carved the formula onto Brougham Bridge [1]. This focus
on deriving a model for three-dimensional spaces, resulted in an algebra that
is most effective in modelling real-world phenomena in our three-dimensional
surroundings. Moreover, the resulting quaternion algebra itself is a four-
dimensional one. Hamilton himself surmised of the use of quaternions for
simultaneous modelling of space and time in a letter to Sir John Herschel [2],
writing; “… how the one of time, of space the three, might in the chain of
symbols girdled be. It is not so much to be wondered at, that they should have
let me to strike out some new lines of research, which former methods have
failed to suggest.” Indeed, quaternions have seen use in quantum mechanics
where, as one example, the quaternion four dimensional algebra allows to
simultaneously formulate rotational and compression gradients [3, 4]. The more
recent use of quaternions for representing quantum fields is somewhat akin to
that used for modelling electromagnetic fields by James Clark Maxwell [5], one
of the main proponents for use of quaternion algebra over vector algebras [1].

Modern applications of the high dimensional division algebra that quaternions
provide include multi-phase power systems [6, 7], body motion tracking [8],
sound processing [9, 10], colour image processing [11], communication
techniques that adopt space-time-polarisation multiplexing [12, 13], and
quantum computing [14]. Before the discussion on quaternions, their associate
division algebra, and the $\mathbb{HR}$-calculus, it is important to introduce
the nomenclature.

Nomenclature: Scalars, column vectors, and matrices are denoted respectively
by lowercase, bold lowercase, and bold uppercase letters, while uppercase bold
italic letters denote multivariate random processes, with calligraphic letters
used to denote sets. The reminder of nomenclature is summarised as follows:

$\mathbb{N}$

natural numbers

$\mathbb{R}$

real-valued numbers

$\mathbb{R}^{+}$

positive real-valued numbers

$\mathbb{C}$

complex numbers

$\mathbb{H}$

quaternion numbers

$\wedge$

logical conjunction

$\Re\left(\cdot\right)$

operator returning the real component

$\Im\left(\cdot\right)$

operator returning the imaginary component

$\\{\imath,\jmath,\kappa\\}$

orthonormal imaginary units spanning the imaginary subspace of $\mathbb{H}$

$\Im_{\chi}\left(\cdot\right)$

operator returning the imaginary component alongside
$\forall\chi\in\mathbb{H}\wedge\chi^{2}=-1$

$\langle\cdot,\cdot\rangle$

inner product operator

$\times$

outer product operator

$\left(\cdot\right)^{*}$

conjugate operator

$\|\cdot\|$

second-order norm

$\left(\cdot\right)^{\mathsf{T}}$

transpose operator

$\left(\cdot\right)^{\mathsf{H}}$

transpose Hermitian operator

$e$

Euler’s number

$\mathsf{E}\left\\{\cdot\right\\}$

statistical expectation operator

$\boldsymbol{\Phi}_{\chi}(\cdot)$

characteristic function of a random vector $\chi$

$\partial$

partial differential operator

$\mathsf{d}$

total differential operator

$\nabla_{\chi}$

gradient operator with respect to $\chi$

$\mathcal{O}(\cdot)$

order of computational complexity

$\mathsf{tr}\left\\{\cdot\right\\}$

trace operator

$\mathbf{I}$

identity matrix of appropriate size

$\mathsf{diag}\left\\{\cdot\right\\}$

constructs a block-diagonal matrix from its entries

$\mathsf{p}\left\\{\cdot\right\\}$

spectral radius operator

###  I-A Quaternion Algebra

A quaternion variable $q\in\mathbb{H}$, consists of a real part, $\Re(q)$, and
a three-dimensional imaginary part, $\Im(q)$. In order to express $\Im\\{q\\}$
in the imaginary subspace of $\mathbb{H}$, three orthonormal imaginary units
that span this space need to be selected. These units are generally denoted by
$\\{\imath,\jmath,\kappa\\}$ and admit the product rules

$\imath\jmath=\kappa,\hskip 3.41418pt\jmath\kappa=\imath,\hskip
3.41418pt\kappa\imath=j,\hskip
3.41418pt\imath^{2}=\jmath^{2}=\kappa^{2}=\imath\jmath\kappa=-1.$ (1)

In this setting, a quaternion $q\in\mathbb{H}$ can be expressed as

$q=\Re\\{q)\\}+\Im\\{q\\}=\Re\\{q\\}+\Im_{\imath}\\{q\\}+\Im_{\jmath}\\{q\\}+\Im_{\kappa}\\{q\\}=q_{r}+\imath
q_{\imath}+\jmath q_{\jmath}+\kappa q_{k}$

where $\\{q_{r},q_{\imath},q_{\jmath},q_{\kappa}\\}\subset\mathbb{R}$.

The product of $q_{1},q_{2}\in\mathbb{H}$ can be found through component-wise
multiplication of their real-valued constituents using the product rules in
(1) and can be alternatively formulated using the inner product,
$\langle\cdot,\cdot\rangle$, and the outer product, $\times$, as

$\displaystyle
q_{1}q_{2}=\Re(q_{1})\Re(q_{2})+\Re(q_{1})\Im(q_{2})+\Re(q_{2})\Im(q_{1})+\Im(q_{1})\times\Im(q_{2})-\big{\langle}\Im(q_{1}),\Im(q_{2})\big{\rangle}.$
(2)

######  Remark 1.

Selection of the orthonormal units that span the imaginary subspace of
$\mathbb{H}$ is not unique. Indeed, a convenient change of these units can be
used to simplify derivation/analysis of quaternion-valued adaptive processing
techniques, see [15, 6, 7] and SM-9A.

######  Remark 2.

As a direct consequence of (1), product operations are not commutative in the
quaternion domain, unless under an exceptional case. From (2), note that if
$q_{1}$ and $q_{2}$ have parallel imaginary parts, by definition
$\Im(q_{1})\times\Im(q_{2})=0$ and thus the product becomes commutative. In
essence, under this circumstance, $q_{1}$ and $q_{2}$ lie in a subspace of
$\mathbb{H}$ spanned by the real axis and the imaginary segment of $q_{1}$ (or
equivalently $q_{2}$). This subspace admits a complex-valued algebra, making
the multiplication commutative. In general, $\forall
q_{1},q_{2}\in\mathbb{H}:q_{1}q_{2}=\left(q^{*}_{2}q^{*}_{1}\right)^{*}$.

The quaternion conjugate is defined as $q^{*}\triangleq\Re(q)-\Im(q)$, which
allows for an inverse to be defined as

$\forall
q\in\mathbb{H}\setminus\\{0\\}:q^{-1}=\frac{q^{*}}{\left\|q\right\|^{2}}$

with $\left\|q\right\|$ denoting the second-order norm of $q$, that is,
$\left\|q\right\|=\sqrt{qq^{*}}=\sqrt{q^{2}_{r}+q^{2}_{\imath}+q^{2}_{\jmath}+q^{2}_{\kappa}}$.

###  I-B Involution and the Augmented Representation

Strictly speaking, an involution is defined as a mapping that is its own
inverse. However, a richer construct has been shown to be useful in the
context of quaternions [16]. For the purposes of this article, an involution
is defined as [16]

$\forall q,\zeta\in\mathbb{H}\hskip 1.0pt\wedge\hskip
1.0pt\\{\Re\\{\zeta\\}=0,\|\zeta\|\neq 0\\}:q^{\zeta}\triangleq\zeta
q\zeta^{-1}$ (3)

and is referred to as the involution of $q$ around $\zeta$. The definition in
(3) is selected as it will allow the real-valued components of quaternions to
be expressed using the linear and invertible mapping111The quaternion
conjugate can also be formulated using involutions as
$q^{*}=\frac{1}{2}\left(q^{\imath}+q^{\jmath}+q^{\kappa}-q\right)$. [17, 18]

$\underbrace{\begin{bmatrix}\mathbf{q}\\\ \mathbf{q}^{\imath}\\\
\mathbf{q}^{\jmath}\\\
\mathbf{q}^{\kappa}\end{bmatrix}}_{\mathbf{q}^{a}}=\underbrace{\begin{bmatrix}\mathbf{I}&\imath\mathbf{I}&\jmath\mathbf{I}&\kappa\mathbf{I}\\\
\mathbf{I}&\imath\mathbf{I}&-\jmath\mathbf{I}&-\kappa\mathbf{I}\\\
\mathbf{I}&-\imath\mathbf{I}&\jmath\mathbf{I}&-\kappa\mathbf{I}\\\
\mathbf{I}&-\imath\mathbf{I}&-\jmath\mathbf{I}&\kappa\mathbf{I}\end{bmatrix}}_{\mathbf{A}}\begin{bmatrix}\mathbf{q}_{r}\\\
\mathbf{q}_{\imath}\\\ \mathbf{q}_{\jmath}\\\
\mathbf{q}_{\kappa}\end{bmatrix}$ (4)

where $\mathbf{q}^{i}$, $\mathbf{q}^{j}$, and $\mathbf{q}^{k}$, denote the
element-wise implementation of the involution,
$\mathbf{A}^{{}^{-1}}=\frac{1}{4}\mathbf{A}^{\mathsf{H}}$, and
$\mathbf{q}^{a}$ is referred to as the augmented quaternion vector.

Finally, in a similar manner to complex-valued numbers, a quaternion
$q\in\mathbb{H}$ can be expressed by its polar presentation as [15]

$q=\|q\|e^{\xi\theta}=\|q\|\big{(}\cos\left(\theta\right)+\xi\sin\left(\theta\right)\big{)}\hskip
7.39755pt\text{where}\hskip 7.39755pt\xi=\frac{\Im(q)}{\|\Im(q)\|}\hskip
7.39755pt\text{and}\hskip
7.39755pt\theta=\text{atan}\left(\frac{\|\Im(q)\|}{\Re(q)}\right).$

Thus, it follows that $\sin(\cdot)$ and $\cos(\cdot)$ functions can be
expressed as

$\sin\left(\theta\right)=\frac{1}{2\xi}\left(e^{\xi\theta}-e^{-\xi\theta}\right)\hskip
14.22636pt\text{and}\hskip 14.22636pt\hskip
5.0pt\cos\left(\theta\right)=\frac{1}{2}\left(e^{\xi\theta}+e^{-\xi\theta}\right)$
(5)

with $\xi^{2}=-1$.

###  I-C Three-Dimensional Rotations

The attitude of a body after having foregone a rotation can be presented using
a single right-hand rotation of that body by an angle of $\theta$ degrees
about an axis $\eta$ parallel to the direction that is unchanged by the
rotation [15]. As shown in Fig. 1, the three-dimensional Cartesian coordinate
system can be considered as representing the imaginary sub-space of a
quaternion skew filed, so that $\eta=\left(\eta_{x},\eta_{y},\eta_{z}\right)$
can be presented as $\eta=\imath\eta_{x}+\jmath\eta_{y}+\kappa\eta_{z}$. In
this setting, the vector presetting the pre- and post-rotation orientation of
the object, $q_{\text{pre}}$ and $q_{\text{post}}$, are related as

$q_{\text{post}}=\mu q_{\text{pre}}\mu^{-1}\hskip 10.243pt\text{with}\hskip
10.243pt\mu=e^{\eta\frac{\theta}{2}}=\cos\left(\frac{\theta}{2}\right)+\eta\sin\left(\frac{\theta}{2}\right).$
(6)

Figure 1: Schematic of a rotation around $\eta$ by an angle of $\theta$, with
$q_{\text{pre}}$ and $q_{\text{post}}$ pointing to the pre- and post-rotation
orientation of the object in question.

######  Remark 3.

The quaternion representation of rotations, as compared to Euler angles and
real-valued rotation matrix algebras, holds a number of advantages. First,
real-valued rotation matrices lose a degree of freedom as one of the rotation
angles reaches $\frac{\pi}{2}$ [15]. However, this is not the case for
quaternions, where only the angle and the axis of rotation are required.
Second, real-valued rotation matrices must have a unit determinant. This is
hard to guarantee after many rotations have been computed due to the finite
accuracy of computer operations. This problem does not arise in the quaternion
representation [19, 15].

######  Remark 4.

Importantly, there seems to be an inherent relation between the formulation in
(6) and the algebraic Lyapunov and Riccati recursions, two of the most
fundamental equations in adaptive signal processing and control [20]. This is
explored in SM-1. Moreover, in Example 1, we show how to leverage the division
algebra of quaternions to model/track status of qubits.

[h]

######  Example 1.

Quaternion-Valued Model of Qubits In quantum computing, the fundamental unit
of information is referred to as a quantum bit or qubit. A qubit, generally
speaking, is a superposition of two bases denoted with $|0\rangle$ and
$|1\rangle$ that are orthogonal to one another and represent state of quantum
particle, e.g., up or down spin of an electron. The qubit can be in any
superposition that satisfies
$\forall\\{\alpha,\beta\\}\subset\mathbb{C}\wedge\|\alpha\|^{2}+\|\beta\|^{2}=1:\hskip
7.39755pt|\psi\rangle=\alpha|0\rangle+\beta|1\rangle$ where $|\psi\rangle$
represents state of the qubit while $\|\alpha\|^{2}$ and $\|\beta\|^{2}$
correspond to the probability of finding the quantum particle in state
$|0\rangle$ or $|1\rangle$. Although use of two complex-vlaued variables
$\alpha$ and $\beta$ might suggest four degrees of freedom for $|\psi\rangle$,
the condition, $\|\alpha\|^{2}+\|\beta\|^{2}=1$, limits the degrees of freedom
to three. Thus, the state of $|\psi\rangle$ can be described via a pure
imaginary quaternion as demonstrated in Fig. 2 where
$|\psi\rangle\leftrightarrow
q=\imath\sin(\theta)\cos(\phi)+\jmath\sin(\theta)\sin(\phi)+\kappa\cos(\theta)$
(7) with $\theta$ and $\phi$ referring to the angles of $|\psi\rangle$ with
the state $|\psi\rangle=|0\rangle$, corresponding to
$\Im_{\kappa}\left\\{q\right\\}$, and angle between the state
$\frac{1}{\sqrt{2}}\left(|0\rangle+|1\rangle\right)$, corresponding to
$\Im_{\imath}\left\\{q\right\\}$, with projection of $|\psi\rangle$ onto the
plane denoted
$\left\\{\frac{1}{\sqrt{2}}\left(|0\rangle+|1\rangle\right),\frac{1}{\sqrt{2}}\left(|0\rangle+\imath|1\rangle\right)\right\\}$,
corresponding to the plane denoted by
$\left\\{\Im_{\imath}\left\\{q\right\\},\Im_{\jmath}\left\\{q\right\\}\right\\}$.
In general, an operation of a quantum computer correspond to manipulations of
the quantum state (groups of quantum states). These operations are effectively
modelled in the quaternion domain using state-space and $\mathbb{HR}$-calculus
tools. For more on this issue see SM-9E.

Figure 2: A qubit representation with Bloch sphere, where
$|\psi\rangle\leftrightarrow
q=\imath\sin(\theta)+\jmath\sin(\theta)\cos(\theta)+\kappa\sin(\theta)\cos(\theta)$.

###  I-D Introduction to Quaternion Statistics

Consider the quaternion-valued random process Q, the statistical description
of which is discernible from the joint characteristic function of its real-
valued components given as

$\mathbf{\Phi}_{\textbf{{Q}}}(\mathbf{s}_{r},\mathbf{s}_{\imath},\mathbf{s}_{\jmath},\mathbf{s}_{\kappa})=\mathsf{E}\left\\{e^{\xi(\mathbf{s}^{\mathsf{T}}_{r}\mathbf{q}_{r}+\mathbf{s}^{\mathsf{T}}_{\imath}\mathbf{q}_{\imath}+\mathbf{s}^{\mathsf{T}}_{\jmath}\mathbf{q}_{\jmath}+\mathbf{s}^{\mathsf{T}}_{\kappa}\mathbf{q}_{\kappa})}\right\\}\cdot$
(8)

where $\xi\in\mathbb{H}\wedge\xi^{2}=-1$. By exploiting the relation in (4),
the characteristic function in (8) is expressed directly in the quaternion
domain as

$\mathbf{\Phi}_{\textbf{{Q}}^{a}}(\mathbf{s}^{a})=\mathsf{E}\left\\{e^{\xi\Re\\{\mathbf{s}^{\mathsf{T}}\mathbf{q}\\}}\right\\}=\mathsf{E}\left\\{e^{(\frac{\xi}{4}\mathbf{s}^{a\mathsf{H}}\mathbf{q}^{a})}\right\\}$
(9)

with
$\mathbf{s}=\mathbf{s}_{r}+\imath\mathbf{s}_{\imath}+\jmath\mathbf{s}_{\jmath}+\kappa\mathbf{s}_{\kappa}$.
Since (9) uses the quaternion random vector in its augmented form, it is
referred to as the augmented quaternion characteristic function (AQCF) [15].
This terminology highlights the importance of the augmented quaternion vector
in formulating quaternion statistic and the $\mathbb{HR}$-calculus. Another
example of using the augmented quaternion vector is the minimum mean square
error (MMSE) estimator presented in SM-2.

All statistical moments, can be extracted from the derivatives of the AQCF.
These moments are in turn useful in the derivation of statistical adaptive
processing techniques. However, this requires a robust calculus, a subject of
this article. Importantly, the AQCF in (9) establishes a basis for statistical
analysis in the quaternion domain. For example, it can be shown that full
statistical information of a Gaussian random processes Q is extractable from
its mean $\mathsf{E}\left\\{\mathbf{q}^{a}\right\\}$ and the augmented
covariance $\mathsf{E}\left\\{\mathbf{q}^{a}\mathbf{q}^{a\mathsf{H}}\right\\}$
[17]. For an in-depth treatment of quaternion statistics, the reader is
refereed to our sister article [21].

######  Remark 5.

Augmenting a quaternion-valued variable with its involution around the
$\imath$, $\jmath$, and $\kappa$ axes, as formulated in (4), forms the
foundations for augmented quaternion statistics and the
$\mathbb{HR}$-calculus, as it provides a comprehensive method for four
dimensional quaternion-valued variables to be considered from four orthogonal
rotational perspectives. Although counter intuitive at first, it can be shown
that $q$, $q^{\imath}$, $q^{\jmath}$, $q^{\kappa}$ hold different statistical
information and can be uncorrelated [17, 15]. Specifically regarding the
$\mathbb{HR}$-calculus, $q$, $q^{\imath}$, $q^{\jmath}$, and $q^{\kappa}$,
must also be considered as algebraically independent [19, 15].

##  II The $\mathbb{HR}$-calculus

Traditionally, a quaternion-valued function
$\mathrm{f}(\mathbf{q}):\mathbb{H}^{M}\rightarrow\mathbb{H}$ is considered
differentiable if and only if it satisfies the Cauchy-Riemann-Fueter condition
[19]

$\frac{\partial\mathrm{f}}{\partial
q_{r}}+\imath\frac{\partial\mathrm{f}}{\partial
q_{\imath}}+\jmath\frac{\partial\mathrm{f}}{\partial
q_{\jmath}}+\kappa\frac{\partial\mathrm{f}}{\partial q_{\kappa}}=0$ (10)

which imposes a severe restriction on admissible functions, as only linear
functions can meet this condition [19, 22]. The $\mathbb{HR}$-calculus [19,
22] presents an elegant solution to this problem, where a quaternion function
$\mathrm{f}\left(\mathbf{q}\right):\mathbb{H}^{M}\rightarrow\mathbb{H}$ is
considered in its augmented formulation as

$\mathrm{f}^{a}(\mathbf{q}^{a})=\begin{bmatrix}\mathrm{f}(\mathbf{q}^{a}),\mathrm{f}^{\imath}(\mathbf{q}^{a}),\mathrm{f}^{\jmath}(\mathbf{q}^{a}),\mathrm{f}^{\kappa}(\mathbf{q}^{a})\end{bmatrix}^{\mathsf{T}}:\mathbb{H}^{4M}\rightarrow\mathbb{H}^{4M}.$
(11)

Then, the total derivative of
$\mathrm{f}\left(\mathbf{q}\right)=\mathrm{f}\left(\mathbf{q}_{r},\mathbf{q}_{\imath},\mathbf{q}_{\jmath},\mathbf{q}_{\kappa}\right)$
is given by

$\mathsf{d}\mathrm{f}=\mathsf{d}\mathbf{q}_{r}\frac{\partial\mathrm{f}}{\partial\mathbf{q}_{r}}+\mathsf{d}\mathbf{q}_{\imath}\frac{\partial\mathrm{f}}{\partial\mathbf{q}_{\imath}}+\mathsf{d}\mathbf{q}_{\jmath}\frac{\partial\mathrm{f}}{\partial\mathbf{q}_{\jmath}}+\mathsf{d}\mathbf{q}_{\kappa}\frac{\partial\mathrm{f}}{\partial\mathbf{q}_{\kappa}}$
(12)

while, on the other hand, treating
$\\{\mathbf{q},\mathbf{q}^{\imath},\mathbf{q}^{\jmath},\mathbf{q}^{\kappa}\\}$
as a set of “algebraically” independent variables and considering the
augmented formulation in (11) yields

$\mathsf{d}\mathrm{f}=\mathsf{d}\mathbf{q}\frac{\partial\mathrm{f}}{\partial\mathbf{q}}+\mathsf{d}\mathbf{q}^{\imath}\frac{\partial\mathrm{f}}{\partial\mathbf{q}^{\imath}}+\mathsf{d}\mathbf{q}^{\jmath}\frac{\partial\mathrm{f}}{\partial\mathbf{q}^{\jmath}}+\mathsf{d}\mathbf{q}^{\kappa}\frac{\partial\mathrm{f}}{\partial\mathbf{q}^{\kappa}}$
(13)

for the total derivative of $\mathrm{f}\left(\mathbf{q}^{a}\right)$. Now,
using (4) in conjunction with (12) and (13) gives the
$\mathbb{HR}$-derivatives in the form of

$\begin{bmatrix}\frac{\partial\mathrm{f}(\mathbf{q},\mathbf{q}^{\imath},\mathbf{q}^{\jmath},\mathbf{q}^{\kappa})}{\partial\mathbf{q}}\\\
\frac{\partial\mathrm{f}(\mathbf{q},\mathbf{q}^{\imath},\mathbf{q}^{\jmath},\mathbf{q}^{\kappa})}{\partial\mathbf{q}^{\imath}}\\\
\frac{\partial\mathrm{f}(\mathbf{q},\mathbf{q}^{\imath},\mathbf{q}^{\jmath},\mathbf{q}^{\kappa})}{\partial\mathbf{q}^{\jmath}}\\\
\frac{\partial\mathrm{f}(\mathbf{q},\mathbf{q}^{\imath},\mathbf{q}^{\jmath},\mathbf{q}^{\kappa})}{\partial\mathbf{q}^{\kappa}}\end{bmatrix}=\frac{1}{4}\mathbf{A}^{*}\begin{bmatrix}\frac{\partial\mathrm{f}(\mathbf{q}_{r},\mathbf{q}_{\imath},\mathbf{q}_{\jmath},\mathbf{q}_{\kappa})}{\partial\mathbf{q}_{r}}\\\
\frac{\partial\mathrm{f}(\mathbf{q}_{r},\mathbf{q}_{\imath},\mathbf{q}_{\jmath},\mathbf{q}_{\kappa})}{\partial\mathbf{q}_{\imath}}\\\
\frac{\partial\mathrm{f}(\mathbf{q}_{r},\mathbf{q}_{\imath},\mathbf{q}_{\jmath},\mathbf{q}_{\kappa})}{\partial\mathbf{q}_{\jmath}}\\\
\frac{\partial\mathrm{f}(\mathbf{q}_{r},\mathbf{q}_{\imath},\mathbf{q}_{\jmath},\mathbf{q}_{\kappa})}{\partial\mathbf{q}_{\kappa}}\end{bmatrix}.$
(14)

The $\mathbb{HR}$-calculus in (14) accommodates the derivation of learning
algorithms that reside naturally in the quaternion domain completely bypassing
the need to deal with quaternion-valued functions as real-valued vector
functions.

From the derivatives of $\mathrm{f}(\cdot)$ in (14), of particular interest to
information processing applications is the conjugate derivative given by

$\frac{\partial\mathrm{f}}{\partial\mathbf{q}^{*}}=\frac{1}{4}\left(\frac{\partial\mathrm{f}}{\partial\mathbf{q}_{r}}+\imath\frac{\partial\mathrm{f}}{\partial\mathbf{q}_{\imath}}+\jmath\frac{\partial\mathrm{f}}{\partial\mathbf{q}_{\jmath}}+\kappa\frac{\partial\mathrm{f}}{\partial\mathbf{q}_{\kappa}}\right)$
(15)

which indicates the direction of maximum rate of change in
$\mathrm{f}(\cdot)$. Thus, presenting the $\mathbb{HR}^{*}$-derivatives as the
gradient operator so that

$\begin{bmatrix}\frac{\partial\mathrm{f}(\mathbf{q},\mathbf{q}^{i},\mathbf{q}^{j},\mathbf{q}^{k})}{\partial\mathbf{q}^{*}}\\\
\frac{\partial\mathrm{f}(\mathbf{q},\mathbf{q}^{i},\mathbf{q}^{j},\mathbf{q}^{k})}{\partial\mathbf{q}^{i*}}\\\
\frac{\partial\mathrm{f}(\mathbf{q},\mathbf{q}^{i},\mathbf{q}^{j},\mathbf{q}^{k})}{\partial\mathbf{q}^{j*}}\\\
\frac{\partial\mathrm{f}(\mathbf{q},\mathbf{q}^{i},\mathbf{q}^{j},\mathbf{q}^{k})}{\partial\mathbf{q}^{k*}}\end{bmatrix}=\frac{1}{4}\mathbf{A}\begin{bmatrix}\frac{\partial\mathrm{f}(\mathbf{q}_{r},\mathbf{q}_{i},\mathbf{q}_{j},\mathbf{q}_{k})}{\partial\mathbf{q}_{r}}\\\
\frac{\partial\mathrm{f}(\mathbf{q}_{r},\mathbf{q}_{i},\mathbf{q}_{j},\mathbf{q}_{k})}{\partial\mathbf{q}_{i}}\\\
\frac{\partial\mathrm{f}(\mathbf{q}_{r},\mathbf{q}_{i},\mathbf{q}_{j},\mathbf{q}_{k})}{\partial\mathbf{q}_{j}}\\\
\frac{\partial\mathrm{f}(\mathbf{q}_{r},\mathbf{q}_{i},\mathbf{q}_{j},\mathbf{q}_{k})}{\partial\mathbf{q}_{k}}\end{bmatrix}\hskip
1.70709pt\text{that is}\hskip
1.70709pt\nabla_{\mathbf{q}^{a*}}\mathrm{f}=\frac{1}{4}\mathbf{A}\nabla\mathrm{f}.$
(16)

At this point, the formulation of a gradient operator and quaternion
derivatives allows us to establish a relation between the AQCF and statistical
moments of quaternion random processes, given in SM-3, formulate statistical
gradient descent, given in Section III-A, and introduce a quaternion-valued
Taylor series expansion, given in Example 2. In what follows, a number rules
essential in deriving adaptive processing techniques are presented.

###  II-A Multiplication, Conjugate, and Rotation Rules in the
$\mathbb{HR}$-Calculus

From (15), it is straightforward to show that [22]

$\forall\nu\in\mathbb{H}\hskip 3.41418pt\&\hskip
3.41418pt\xi\in\\{\imath,\jmath,\kappa\\}:\hskip
3.41418pt\frac{\partial\nu\mathrm{f}\left(\mathbf{q}^{a}\right)}{\partial
q^{\xi}}=\nu\frac{\partial\mathrm{f}\left(\mathbf{q}^{a}\right)}{\partial
q^{\nu^{-1}\xi}}\hskip 7.39755pt\text{and}\hskip
7.39755pt\frac{\partial\mathrm{f}\left(\mathbf{q}^{a}\right)\nu}{\partial
q^{\xi}}=\frac{\partial\mathrm{f}\left(\mathbf{q}^{a}\right)}{\partial
q^{\xi}}\nu\cdot$ (17)

In the same vein, we can formulate the effect of multiplication from both
sides, yielding [22]

$\forall\nu\in\mathbb{H}\wedge\nu^{2}=-1\hskip 3.41418pt\&\hskip
3.41418pt\xi\in\\{\imath,\jmath,\kappa\\}:\hskip
3.41418pt\left(\frac{\partial\mathrm{f}\left(\mathbf{q}^{a}\right)}{\partial
q^{\xi}}\right)^{\nu}=\nu\frac{\partial\mathrm{f}\left(\mathbf{q}^{a}\right)}{\partial
q^{\xi}}\nu^{-1}=\frac{\partial\mathrm{f}^{\nu}\left(\mathbf{q}^{a}\right)}{\partial
q^{\nu\xi}}\cdot$ (18)

Once more, considering (15), we have

$\left(\frac{\partial\mathrm{f}\left(\mathbf{q}^{a}\right)}{\partial\mathbf{q}^{*}}\right)^{*}=\frac{1}{4}\left(\frac{\partial\mathrm{f}^{*}\left(\mathbf{q}^{a}\right)}{\partial\mathbf{q}_{r}}-\frac{\partial\mathrm{f}^{*}\left(\mathbf{q}^{a}\right)}{\partial\mathbf{q}_{\imath}}\imath-\frac{\partial\mathrm{f}^{*}\left(\mathbf{q}^{a}\right)}{\partial\mathbf{q}_{\jmath}}\jmath-\frac{\partial\mathrm{f}^{*}\left(\mathbf{q}^{a}\right)}{\partial\mathbf{q}_{\kappa}}\kappa\right)=\frac{\partial_{\text{left}}\mathrm{f}^{*}\left(\mathbf{q}^{a}\right)}{\partial\mathbf{q}}$
(19)

where $\partial_{\text{left}}$ denotes the placement of partial differential
on the left side of the imaginary units $\\{\imath,\jmath,\kappa\\}$. This can
also be seen from (16) which yields
$\left(\nabla_{\mathbf{q}^{a*}}\mathrm{f}\right)^{\mathsf{H}}=\frac{1}{4}\left(\nabla\mathrm{f}\right)^{\mathsf{T}}\mathbf{A}^{\mathsf{H}}$.
These results can now be used to establish a chain derivative rule in the
$\mathbb{HR}$-calculus. This is addressed in Section II-B.

Now, consider
$\\{\mathrm{f}\left(\cdot\right),\mathrm{g}\left(\cdot\right)\\}:\mathbb{H}^{M}\rightarrow\mathbb{H}$;
then,
$\mathrm{h}\left(\mathbf{q}\right)=\mathrm{f}\left(\mathbf{q}\right)\mathrm{g}\left(\mathbf{q}\right)$
can be expressed in terms of the real-valued components of
$\mathrm{f}\left(\cdot\right)$ and $\mathrm{g}\left(\cdot\right)$ as

$\mathrm{h}\left(\mathbf{q}\right)=\left(\Re\left\\{\mathrm{f}\left(\mathbf{q}\right)\right\\}+\dots+\kappa\Im_{\kappa}\left\\{\mathrm{f}\left(\mathbf{q}\right)\right\\}\right)\left(\Re\left\\{\mathrm{g}\left(\mathbf{q}\right)\right\\}+\dots+\kappa\Im_{\kappa}\left\\{\mathrm{g}\left(\mathbf{q}\right)\right\\}\right)$
(20)

with the exact expression given in SM-4. Now,
$\forall\xi\in\\{1,\imath,\jmath,\kappa\\}$,
$\frac{\partial\mathrm{h}\left(\mathbf{q}\right)}{\partial\mathbf{q}^{\xi}}$
can be expressed in terms of its real-valued components, allowing us to use
the framework of (10)-(16) in order to formulate the derivative [22]. In this
setting, after some tedious mathematical manipulations we arrive at

$\forall\xi\in\\{1,\imath,\jmath,\kappa\\}:\left\\{\begin{aligned}
\frac{\partial\mathrm{h}\left(\mathbf{q}\right)}{\partial\mathbf{q}^{\xi}}=&\mathrm{f}\left(\mathbf{q}\right)\frac{\partial\mathrm{g}\left(\mathbf{q}\right)}{\partial\mathbf{q}^{\mathrm{f}^{*}\left(\mathbf{q}\right)\xi}}+\frac{\partial\mathrm{f}\left(\mathbf{q}\right)}{\partial\mathbf{q}^{\xi}}\mathrm{g}\left(\mathbf{q}\right)\\\
\frac{\partial\mathrm{h}\left(\mathbf{q}\right)}{\partial\mathbf{q}^{\xi*}}=&\mathrm{f}\left(\mathbf{q}\right)\frac{\partial\mathrm{g}\left(\mathbf{q}\right)}{\partial\left(\mathbf{q}^{\mathrm{f}^{*}\left(\mathbf{q}\right)\xi}\right)^{*}}+\frac{\partial\mathrm{f}\left(\mathbf{q}\right)}{\partial\mathbf{q}^{\xi*}}\mathrm{g}\left(\mathbf{q}\right).\end{aligned}\right.$
(21)

######  Remark 6.

From (19), note that an equivalent calculus can be derived via replacing
$\partial$ with $\partial_{\text{left}}$, see [22]. Although both approaches
yield equivalent results, a skilful use of both approaches can simplify
derived algorithms, see SM-5.

[h]

######  Example 2.

Quaternion-Valued Taylor Series Expansion Consider the first-order Taylor
series expansion in its real-valued formulation given by
$\mathrm{f}\left(\mathbf{x}+\Delta\mathbf{x}\right)=\mathrm{f}\left(\mathbf{x}\right)+\left(\Delta\mathbf{x}\right)^{\mathsf{T}}\left(\nabla_{\mathbf{x}}\mathrm{f}\left(\mathbf{x}\right)\right)+\mathcal{O}\left(\mathrm{f}\left(\mathbf{x}\right)\right)$
(22) where
$\mathbf{x}^{\mathsf{T}}=\left[\mathbf{x}^{\mathsf{T}}_{r},\mathbf{x}^{\mathsf{T}}_{\imath},\mathbf{x}^{\mathsf{T}}_{\jmath},\mathbf{x}^{\mathsf{T}}_{\kappa}\right]$
and $\Delta\mathbf{x}$ represents the change in $\mathbf{x}$. Now, mapping
$\mathbf{x}$, $\Delta\mathbf{x}$, and
$\nabla_{\mathbf{x}}\mathrm{f}\left(\mathbf{x}\right)$ on to the quaternion
domain using (29) and (16) allows (22) to be formulated in the quaternion
domain as [19]
$\mathrm{f}\left(\mathbf{x}^{a}+\Delta\mathbf{x}^{a}\right)=\mathrm{f}\left(\mathbf{x}^{a}\right)+\left(\Delta\mathbf{x}^{a}\right)^{\mathsf{H}}\left(\nabla_{\mathbf{x}^{a*}}\mathrm{f}\right)+\mathcal{O}\left(\mathrm{f}\left(\mathbf{x}^{a}\right)\right).$
One can formulate higher-order expansion of the series in a similar manner
[22].

###  II-B Chain Derivative Rule in the $\mathbb{HR}$-Calculus

Consider the quaternion-valued compound function
$\mathrm{f}\left(\mathrm{g}\left(\cdot\right)\right):\mathbb{H}^{M}\rightarrow\mathbb{H}$,
where $\mathrm{f}\left(\cdot\right):\mathbb{H}^{M}\rightarrow\mathbb{H}$ and
$\mathrm{g}\left(\cdot\right):\mathbb{H}^{M}\rightarrow\mathbb{H}^{M}$. One
can consider the relation between $\mathrm{g}^{a}\left(\mathbf{q}\right)$ and
the real-valued components of $\mathrm{g}\left(\mathbf{q}\right)$, akin to
that given in (4) and (14), yielding

$\begin{bmatrix}\frac{\partial\mathrm{f}\left(\mathrm{g}^{a}\left(\mathbf{q}\right)\right)}{\partial\mathrm{g}\left(\mathbf{q}\right)}\\\
\frac{\partial\mathrm{f}\left(\mathrm{g}^{a}\left(\mathbf{q}\right)\right)}{\partial\mathrm{g}^{\imath}\left(\mathbf{q}\right)}\\\
\frac{\partial\mathrm{f}\left(\mathrm{g}^{a}\left(\mathbf{q}\right)\right)}{\partial\mathrm{g}^{\jmath}\left(\mathbf{q}\right)}\\\
\frac{\partial\mathrm{f}\left(\mathrm{g}^{a}\left(\mathbf{q}\right)\right)}{\partial\mathrm{g}^{\kappa}\left(\mathbf{q}\right)}\end{bmatrix}=\frac{1}{4}\mathbf{A}^{*}\begin{bmatrix}\frac{\partial\mathrm{f}\left(\mathrm{g}^{a}\left(\mathbf{q}\right)\right)}{\partial\mathrm{g}_{r}\left(\mathbf{q}\right)}\\\
\frac{\partial\mathrm{f}\left(\mathrm{g}^{a}\left(\mathbf{q}\right)\right)}{\partial\mathrm{g}_{\imath}\left(\mathbf{q}\right)}\\\
\frac{\partial\mathrm{f}\left(\mathrm{g}^{a}\left(\mathbf{q}\right)\right)}{\partial\mathrm{g}_{\jmath}\left(\mathbf{q}\right)}\\\
\frac{\partial\mathrm{f}\left(\mathrm{g}^{a}\left(\mathbf{q}\right)\right)}{\partial\mathrm{g}_{\kappa}\left(\mathbf{q}\right)}\end{bmatrix}\hskip
7.39755pt\text{with}\hskip
7.39755pt\begin{bmatrix}\frac{\partial\mathrm{g}_{r}\left(\mathbf{q}\right)}{\partial\mathbf{q}}\\\
\frac{\partial\mathrm{g}_{\imath}\left(\mathbf{q}\right)}{\partial\mathbf{q}}\\\
\frac{\partial\mathrm{g}_{\jmath}\left(\mathbf{q}\right)}{\partial\mathbf{q}}\\\
\frac{\partial\mathrm{g}_{\kappa}\left(\mathbf{q}\right)}{\partial\mathbf{q}}\end{bmatrix}=\frac{1}{4}\mathbf{A}^{\mathsf{H}}\begin{bmatrix}\frac{\partial\mathrm{g}\left(\mathbf{q}\right)}{\partial\mathbf{q}}\\\
\frac{\partial\mathrm{g}^{\imath}\left(\mathbf{q}\right)}{\partial\mathbf{q}}\\\
\frac{\partial\mathrm{g}^{\jmath}\left(\mathbf{q}\right)}{\partial\mathbf{q}}\\\
\frac{\partial\mathrm{g}^{\kappa}\left(\mathbf{q}\right)}{\partial\mathbf{q}}\end{bmatrix}\cdot$
(23)

Following the framework set in (10)-(16) gives

$\displaystyle\forall\xi\in\\{\imath,\jmath,\kappa\\}:\hskip
3.41418pt\frac{\partial\mathrm{f}\left(\mathrm{g}\left(\mathbf{q}\right)\right)}{\partial\mathbf{q}^{\xi}}=$
$\displaystyle\frac{\partial\mathrm{g}_{r}\left(\mathbf{q}\right)}{\partial\mathbf{q}^{\xi}}\frac{\partial\mathrm{f}\left(\mathrm{g}\left(\mathbf{q}\right)\right)}{\partial\mathrm{g}_{r}\left(\mathbf{q}\right)}+\frac{\partial\mathrm{g}_{\imath}\left(\mathbf{q}\right)}{\partial\mathbf{q}^{\xi}}\frac{\partial\mathrm{f}\left(\mathrm{g}\left(\mathbf{q}\right)\right)}{\partial\mathrm{g}_{\imath}\left(\mathbf{q}\right)}$
(24)
$\displaystyle+\frac{\partial\mathrm{g}_{\jmath}\left(\mathbf{q}\right)}{\partial\mathbf{q}^{\xi}}\frac{\partial\mathrm{f}\left(\mathrm{g}\left(\mathbf{q}\right)\right)}{\partial\mathrm{g}_{\jmath}\left(\mathbf{q}\right)}+\frac{\partial\mathrm{g}_{\kappa}\left(\mathbf{q}\right)}{\partial\mathbf{q}^{\xi}}\frac{\partial\mathrm{f}\left(\mathrm{g}\left(\mathbf{q}\right)\right)}{\partial\mathrm{g}_{\kappa}\left(\mathbf{q}\right)}\cdot$

Upon substituting (23) into (24) and using the rotation rule in (18) we arrive
at

$\displaystyle\forall\xi\in\\{\imath,\jmath,\kappa\\}:\hskip
3.41418pt\frac{\partial\mathrm{f}\left(\mathrm{g}\left(\mathbf{q}\right)\right)}{\partial\mathbf{q}^{\xi}}=$
$\displaystyle\frac{\partial\mathrm{g}\left(\mathbf{q}\right)}{\partial\mathbf{q}^{\xi}}\frac{\partial\mathrm{f}\left(\mathrm{g}\left(\mathbf{q}\right)\right)}{\partial\mathrm{g}\left(\mathbf{q}\right)}+\frac{\partial\mathrm{g}^{\imath}\left(\mathbf{q}\right)}{\partial\mathbf{q}^{\xi}}\frac{\partial\mathrm{f}\left(\mathrm{g}\left(\mathbf{q}\right)\right)}{\partial\mathrm{g}^{\imath}\left(\mathbf{q}\right)}$
(25)
$\displaystyle+\frac{\partial\mathrm{g}^{\jmath}\left(\mathbf{q}\right)}{\partial\mathbf{q}^{\xi}}\frac{\partial\mathrm{f}\left(\mathrm{g}\left(\mathbf{q}\right)\right)}{\partial\mathrm{g}^{\jmath}\left(\mathbf{q}^{\xi}\right)}+\frac{\partial\mathrm{g}^{\kappa}\left(\mathbf{q}\right)}{\partial\mathbf{q}^{\xi}}\frac{\partial\mathrm{f}\left(\mathrm{g}\left(\mathbf{q}\right)\right)}{\partial\mathrm{g}^{\kappa}\left(\mathbf{q}\right)}$

with the direction of steepest change given by

$\displaystyle\frac{\partial\mathrm{f}\left(\mathrm{g}\left(\mathbf{q}\right)\right)}{\partial\mathbf{q}^{*}}=$
$\displaystyle\frac{\partial\mathrm{g}^{*}\left(\mathbf{q}\right)}{\partial\mathbf{q}^{*}}\frac{\partial\mathrm{f}\left(\mathrm{g}\left(\mathbf{q}\right)\right)}{\partial\mathrm{g}^{*}\left(\mathbf{q}\right)}+\frac{\partial\mathrm{g}^{\imath*}\left(\mathbf{q}\right)}{\partial\mathbf{q}^{*}}\frac{\partial\mathrm{f}\left(\mathrm{g}^{\imath}\left(\mathbf{q}\right)\right)}{\partial\mathrm{g}^{\imath*}\left(\mathbf{q}\right)}$
(26)
$\displaystyle+\frac{\partial\mathrm{g}^{\jmath*}\left(\mathbf{q}\right)}{\partial\mathbf{q}^{*}}\frac{\partial\mathrm{f}\left(\mathrm{g}^{\jmath}\left(\mathbf{q}\right)\right)}{\partial\mathrm{g}^{\jmath*}\left(\mathbf{q}\right)}+\frac{\partial\mathrm{g}^{\kappa*}\left(\mathbf{q}\right)}{\partial\mathbf{q}^{*}}\frac{\partial\mathrm{f}\left(\mathrm{g}^{\kappa}\left(\mathbf{q}\right)\right)}{\partial\mathrm{g}^{\kappa*}\left(\mathbf{q}\right)}\cdot$

Most useful in adaptive information processing applications is the case where
$\mathrm{g}\left(\cdot\right):\mathbb{H}^{M}\rightarrow\mathbb{R}$ and
$\mathrm{f}\left(\cdot\right):\mathbb{R}\rightarrow\mathbb{R}$. In this case,
from (24)-(26), we have

$\forall\xi\in\\{\imath,\jmath,\kappa\\}:\hskip
3.41418pt\frac{\partial\mathrm{f}\left(\mathrm{g}\left(\mathbf{q}\right)\right)}{\partial\mathbf{q}^{\xi}}=\frac{\partial\mathrm{g}\left(\mathbf{q}\right)}{\partial\mathbf{q}^{\xi}}\frac{\partial\mathrm{f}\left(\mathrm{g}\left(\mathbf{q}\right)\right)}{\partial\mathrm{g}\left(\mathbf{q}\right)}\hskip
3.41418pt\text{and}\hskip
3.41418pt\frac{\partial\mathrm{f}\left(\mathrm{g}\left(\mathbf{q}\right)\right)}{\partial\mathbf{q}^{\xi*}}=\frac{\partial\mathrm{g}\left(\mathbf{q}\right)}{\partial\mathbf{q}^{\xi*}}\frac{\partial\mathrm{f}\left(\mathrm{g}\left(\mathbf{q}\right)\right)}{\partial\mathrm{g}\left(\mathbf{q}\right)}\cdot$
(27)

The chain derivative rule allows for the derivation of quaternion-valued
adaptive learning techniques. The case of linear (cf. nonlinear) adaptive
filters is presented in SM-5 (cf. SM-6), while its application in quaternion-
vlaued backpropagation demonstrated in Section III-C.

##  III Learning and Decision Making in the Quaternion Domain

###  III-A Statistical Gradient Descent

Consider the problem of tracking the state of an evolving system, whether a
physical system (aircraft) or state of adaptive processing machinery (weights
of a quaternion-valued neural network). The state of this system is
represented via the state vector sequence
$\\{\mathbf{x}_{\left[n\right]}:n\in\mathbb{N}\\}$ so that

$\mathbf{x}_{\left[n+1\right]}=\mathrm{f}\left(\mathbf{x}_{\left[n\right]},\mathbf{v}_{\left[n\right]},\mathbf{u}_{\left[n\right]}\right)$
(28)

where $\mathrm{f}(\cdot)$, $\mathbf{u}_{\left[n\right]}$, and
$\mathbf{v}_{\left[n\right]}$, are the evolution function, a known control
input, and a noise vector at time instant $n$. Informative observations
regarding the state vector are assumed available through observations

$\mathbf{y}_{\left[n\right]}=\mathrm{h}\left(\mathbf{x}_{\left[n\right]},\mathbf{w}_{\left[n\right]}\right)$
(29)

where $\mathrm{h}(\cdot)$ denotes the observation function, while
$\mathbf{y}_{\left[n\right]}$ and $\mathbf{w}_{\left[n\right]}$ are the
observation and observation noise at time instant $n$. The aim is to find the
estimate of $\mathbf{x}_{\left[n\right]}$ given the available information set
$\mathbfcal{I}_{\left[n\right]}=\\{\mathbf{u}_{\left[1:n\right]},\mathbf{y}_{\left[1:n\right]}\\}$,
formulated as
$\hat{\mathbf{x}}_{\left[n\right]}=\mathsf{E}\left\\{\mathbf{x}_{\left[n\right]}|\mathbfcal{I}_{\left[n\right]}\right\\}$.

From (9) and using (15), we have

$\hat{\mathbf{x}}^{a}_{\left[n\right]}=\frac{4}{\xi}\frac{\partial\boldsymbol{\Phi}_{\textbf{{x}}^{a}_{\left[n\right]}|\mathbfcal{I}_{\left[n\right]}}(\mathbf{s}^{a})}{\partial\mathbf{s}^{a*}}{{{|}}}_{\mathbf{s=0}}=\mathsf{E}\left\\{\mathbf{x}^{a}_{\left[n\right]}|\mathbfcal{I}_{\left[n\right]}\right\\}\cdot$

In the same vain, the a priori estimate is defined as

$\boldsymbol{\psi}^{a}_{\left[n+1\right]}=\frac{4}{\xi}\frac{\partial\boldsymbol{\Phi}_{\textbf{{x}}^{a}_{\left[n+1\right]}|\mathbfcal{I}_{\left[n\right]}}(\mathbf{s}^{a})}{\partial\mathbf{s}^{a*}}{{{|}}}_{\mathbf{s=0}}=\mathsf{E}\left\\{\mathbf{x}^{a}_{\left[n+1\right]}|\mathbfcal{I}_{\left[n\right]}\right\\}\cdot$
(30)

where using (28) yields

$\boldsymbol{\Phi}_{\textbf{{x}}^{a}_{\left[n+1\right]}|\mathbfcal{I}_{\left[n\right]}}(\mathbf{s}^{a})=\mathsf{E}\left\\{e^{\frac{\xi}{4}\mathbf{s}^{a\mathsf{H}}\mathbf{x}^{a}_{\left[n+1\right]}}|\mathbfcal{I}_{\left[n\right]}\right\\}=\mathsf{E}\left\\{e^{\frac{\xi}{4}\mathbf{s}^{a\mathsf{H}}\mathrm{f}^{a}(\mathbf{x}_{\left[n\right]},\mathbf{v}_{\left[n\right]},\mathbf{u}_{\left[n\right]})}|\mathbf{y}_{\left[1:n\right]}\right\\}\cdot$
(31)

A substitution of (30) into (31) gives

$\boldsymbol{\psi}^{a}_{\left[n+1\right]}=\mathsf{E}\left\\{\mathrm{f}^{a}(\mathbf{x}_{\left[n\right]},\mathbf{v}_{\left[n\right]},\mathbf{u}_{\left[n\right]})|\mathbf{y}_{\left[1:n\right]}\right\\}\approx\mathsf{E}\left\\{\mathrm{f}^{a}(\mathbf{x}_{\left[n\right]},\mathbf{u}_{\left[n\right]})|\mathbf{y}_{\left[1:n\right]}\right\\}$
(32)

allowing for an approximation of $\boldsymbol{\psi}^{a}_{\left[n+1\right]}$ to
be formulated as

$\boldsymbol{\psi}^{a}_{\left[n+1\right]}\approx\mathsf{E}\left\\{\mathbf{A}_{\left[n\right]}\hat{\mathbf{x}}^{a}_{\left[n\right]}+\mathcal{O}_{\mathrm{f}(\mathbf{x}_{\left[n\right]},\mathbf{u}_{\left[n\right]})}|\mathbf{y}_{\left[1:n\right]}\right\\}\approx\mathbf{A}_{\left[n\right]}\hat{\mathbf{x}}^{a}_{\left[n\right]}+\mathsf{E}\left\\{\mathcal{O}_{\mathrm{f}(\mathbf{x}_{\left[n\right]},\mathbf{u}_{\left[n\right]})}|\mathbf{y}_{\left[1:n\right]}\right\\}$
(33)

where $\mathbf{A}_{\left[n\right]}$ represents the linearisation of
$\mathrm{f}(\cdot)$ at
$\\{\hat{\mathbf{x}}^{a}_{\left[n\right]},\mathbf{u}^{a}_{\left[n\right]}\\}$
obtainable through the gradient operator (16). In addition, we have

$\hat{\mathbf{x}}^{a}_{\left[n+1\right]}=\frac{4}{\xi}\frac{\partial\boldsymbol{\Phi}_{\mathbf{x}^{a}_{\left[n+1\right]}|\mathbfcal{I}_{\left[n+1\right]}}(\mathbf{s})}{\partial\mathbf{s}^{a*}}$

with
$\boldsymbol{\Phi}_{\mathbf{x}^{a}_{\left[n+1\right]}|\mathbfcal{I}_{\left[n+1\right]}}(\mathbf{s})=\mathsf{E}\left\\{e^{\frac{\xi}{4}\mathbf{s}^{a\mathsf{H}}\mathrm{f}^{a}(\mathbf{x}_{\left[n\right]},\mathbf{u}_{\left[n\right]})}e^{\frac{\xi}{4}\mathbf{s}^{a\mathsf{H}}\mathcal{O}_{(\mathrm{f}(\mathbf{x}_{\left[n\right]},\mathbf{v}_{\left[n\right]},\mathbf{u}_{\left[n\right]}))}}|\mathbf{y}_{\left[1:n+1\right]}\right\\}$.

Assuming that dependency of
$\mathcal{O}_{(\mathrm{f}(\mathbf{x}_{\left[n\right]},\mathbf{v}_{\left[n\right]},\mathbf{u}_{\left[n\right]}))}$
on $\mathrm{f}(\mathbf{x}_{\left[n\right]},\mathbf{u}_{\left[n\right]})$ is
negligible, (33) can be reformulated to give

$\boldsymbol{\Phi}_{\mathbf{x}^{a}_{\left[n+1\right]}|\mathbfcal{I}_{\left[n+1\right]}}(\mathbf{s}^{a})=\boldsymbol{\Phi}_{\mathbf{x}^{a}_{\left[n+1\right]}|\mathbfcal{I}_{\left[n\right]}}(\mathbf{s}^{a})\mathsf{E}\left\\{e^{\frac{\xi}{4}\mathbf{s}^{a\mathsf{H}}\mathcal{O}_{(\mathrm{f}(\mathbf{x}_{\left[n\right]},\mathbf{v}_{\left[n\right]},\mathbf{u}_{\left[n\right]})}}|\mathbf{y}_{\left[n+1\right]}\right\\}$

where
$\boldsymbol{\Phi}_{\mathbf{x}^{a}_{\left[n+1\right]}|\mathbfcal{I}_{\left[n+1\right]}}(\mathbf{s}^{a})$
consists of two parts, one representing the a priori,
$\boldsymbol{\Phi}_{\mathbf{x}^{a}_{\left[n+1\right]}|\mathbfcal{I}_{\left[n\right]}}(\mathbf{s}^{a})$,
and the other, representing the a posteriori,
$\mathsf{E}\left\\{e^{j\mathbf{s}^{\mathsf{T}}\mathcal{O}_{(\mathrm{f}_{\left[n\right]}(\mathbf{x}_{\left[n\right]},\mathbf{v}_{\left[n\right]},\mathbf{u}_{\left[n\right]}))}}|\mathbf{y}_{\left[n+1\right]}\right\\}$,
correction. Moreover, from (32), we have

$\hat{\mathbf{x}}^{a}_{\left[n+1\right]}=\boldsymbol{\psi}^{a}_{\left[n+1\right]}+\mathsf{E}\left\\{\mathcal{O}_{(\mathrm{f}_{\left[n\right]}(\mathbf{x}_{\left[n\right]},\mathbf{v}_{\left[n\right]},\mathbf{u}_{\left[n\right]}))}|\mathbf{y}_{\left[n+1\right]}\right\\}\cdot$
(34)

Although the estimate, $\hat{\mathbf{x}}^{a}_{\left[n+1\right]}$, in (34) can
be approximated using particle based inference techniques, the
$\mathbb{HR}$-calculus can provide a more computationally friendly and
mathematically tractable solutions. From (28) and (29), the following error
terms are defined

$\boldsymbol{\epsilon}_{\left[n\right]}=\mathbf{x}_{\left[n\right]}-\hat{\mathbf{x}}_{\left[n\right]}\hskip
6.82881pt\text{and}\hskip
6.82881pt\tilde{\boldsymbol{\epsilon}}_{\left[n\right]}=\mathbf{y}_{\left[n\right]}-\hat{\mathbf{y}}_{\left[n\right]}$
(35)

where
$\hat{\mathbf{y}}_{\left[n+1\right]}=\mathrm{h}(\boldsymbol{\psi}_{\left[n+1\right]})$
is projection of $\boldsymbol{\psi}_{\left[n+1\right]}$ onto the observation
space. Now, let $\mathrm{d}(\cdot)$ be a quaternion-valued metric function.
Then, $\\{\boldsymbol{\epsilon}_{\left[n\right]}:n\in\mathbb{N}\\}$ being
exponentially stable sequence implies that
$\\{\tilde{\boldsymbol{\epsilon}}_{\left[n\right]}:n\in\mathbb{N}\\}$ should
also be exponentially stable.222Note that; i) converse of the statement might
not hold true, and ii) the statement implicitly implies that
$\mathrm{h}(\cdot)$ and $\mathrm{f}(\cdot)$ are bounded. In this setting,
$\nabla_{\boldsymbol{\psi}^{a*}_{\left[n+1\right]}}\mathrm{d}(\tilde{\boldsymbol{\epsilon}}_{\left[n+1\right]})$
indicates a direction of ascent for
$\mathrm{d}(\tilde{\boldsymbol{\epsilon}}_{\left[n+1\right]})$, and hence, an
adaptive filtering structure is formulated as

$\hat{\mathbf{x}}^{a}_{\left[n+1\right]}\approx\boldsymbol{\psi}^{a}_{\left[n+1\right]}-\mathbf{G}_{\left[n\right]}\mathsf{E}\left\\{\nabla_{\boldsymbol{\psi}^{a*}_{\left[n+1\right]}}\mathrm{d}(\tilde{\boldsymbol{\epsilon}}_{\left[n+1\right]})\right\\}\approx\boldsymbol{\psi}^{a}_{\left[n+1\right]}-\mathbf{G}_{\left[n\right]}\nabla_{\boldsymbol{\psi}^{a*}_{\left[n+1\right]}}\mathrm{d}(\tilde{\boldsymbol{\epsilon}}_{\left[n+1\right]})$
(36)

where $\\{\mathbf{G}_{\left[n\right]}:n\in\mathbb{N}\\}$ is a gain sequence
while
$\mathsf{E}\left\\{\nabla_{\boldsymbol{\psi}^{a*}_{\left[n+1\right]}}\mathrm{d}(\tilde{\boldsymbol{\epsilon}}_{\left[n+1\right]})\right\\}$
is replaced with its instantaneous value,
$\nabla_{\boldsymbol{\psi}^{a*}_{\left[n+1\right]}}\mathrm{d}(\tilde{\boldsymbol{\epsilon}}_{\left[n+1\right]})$,
as in most cases these expectations are not available.

######  Remark 7.

Most adaptive filtering techniques can be formalised within the framework
derived in this section. As an example, the stochastic gradient descent
technique is used to derive the quaternion-valued Kalman filter in Example 3.

[h]

######  Example 3.

Quaternion-Valued Kalman Filtering From Section III-A, let the evolution,
observation, and metric functions be
$\displaystyle\mathrm{f}\left(\mathbf{x}_{\left[n\right]},\mathbf{v}_{\left[n\right]},\mathbf{u}_{\left[n\right]}\right)=$
$\displaystyle\mathbf{A}\mathbf{x}^{a}_{\left[n\right]}+\mathbf{B}\mathbf{u}^{a}_{\left[n\right]}+\mathbf{v}^{a}_{\left[n\right]}$
(37)
$\displaystyle\mathrm{h}\left(\mathbf{x}_{\left[n\right]},\mathbf{w}_{\left[n\right]}\right)=$
$\displaystyle\mathbf{H}\mathbf{x}^{a}_{\left[n\right]}+\mathbf{w}^{a}_{\left[n\right]}$
$\displaystyle\mathrm{d}\left(\tilde{\boldsymbol{\epsilon}}_{\left[n\right]}\right)=$
$\displaystyle
4\left\|\tilde{\boldsymbol{\epsilon}}_{\left[n\right]}\right\|^{2}=\mathsf{tr}\left\\{\tilde{\boldsymbol{\epsilon}}^{a}_{\left[n\right]}\tilde{\boldsymbol{\epsilon}}^{a\mathsf{H}}_{\left[n\right]}\right\\}.$
In this setting, from the framework of the $\mathbb{HR}$-calculus and (36), we
have
$\hat{\mathbf{x}}^{a}_{\left[n+1\right]}=\underbrace{\mathbf{A}\hat{\mathbf{x}}^{a}_{\left[n\right]}}_{\boldsymbol{\psi}^{a}_{\left[n+1\right]}}+\mathbf{G}_{\left[n\right]}\left(\mathbf{y}_{\left[n\right]}-\underbrace{\mathbf{H}\hat{\mathbf{x}}^{a}_{\left[n\right]}}_{\hat{\mathbf{y}}_{\left[n\right]}}\right).$
(38) Using (37) and (38), the evolution of
$\boldsymbol{\epsilon}^{a}_{\left[n\right]}$ can be formulate as
$\boldsymbol{\epsilon}^{a}_{\left[n+1\right]}=\left(\mathbf{I}-\mathbf{G}_{\left[n\right]}\mathbf{H}\right)\mathbf{A}\boldsymbol{\epsilon}^{a}_{\left[n\right]}+\left(\mathbf{I}-\mathbf{G}_{\left[n\right]}\mathbf{H}\right)\mathbf{v}^{a}_{\left[n\right]}+\mathbf{G}_{n}\mathbf{w}^{a}_{\left[n\right]}$
which in turn yields
$\displaystyle\boldsymbol{\epsilon}^{a}_{\left[n+1\right]}\boldsymbol{\epsilon}^{a\mathsf{H}}_{\left[n+1\right]}=$
$\displaystyle\left(\mathbf{I}-\mathbf{G}_{\left[n\right]}\mathbf{H}\right)\mathbf{A}\boldsymbol{\epsilon}^{a}_{\left[n\right]}\boldsymbol{\epsilon}^{a\mathsf{H}}_{\left[n\right]}\mathbf{A}^{\mathsf{H}}\left(\mathbf{I}-\mathbf{G}_{\left[n\right]}\mathbf{H}\right)^{\mathsf{H}}$
(39)
$\displaystyle+\left(\mathbf{I}-\mathbf{G}_{\left[n\right]}\mathbf{H}\right)\mathbf{v}^{a}_{\left[n\right]}\mathbf{v}^{a\mathsf{H}}_{\left[n\right]}\left(\mathbf{I}-\mathbf{G}_{\left[n\right]}\mathbf{H}\right)^{\mathsf{H}}+\mathbf{G}_{n}\mathbf{w}^{a}_{\left[n\right]}\mathbf{w}^{a\mathsf{H}}_{\left[n\right]}\mathbf{G}^{\mathsf{H}}_{\left[n\right]}.$
Taking the statistical expectation of (39) and solving for
$\partial\mathsf{E}\left\\{\mathrm{d}\left(\boldsymbol{\epsilon}_{\left[n\right]}\right)\right\\}/\partial\mathbf{G}_{\left[n\right]}=0$
yields the optimal gain as
$\mathbf{G}_{\left[n\right]}=\mathbf{M}_{\left[n\right]}\mathbf{H}^{\mathsf{H}}\boldsymbol{\Sigma}_{\mathbf{w}^{a}_{\left[n\right]}}^{-1}$
(40) where
$\boldsymbol{\Sigma}_{\mathbf{w}^{a}_{\left[n\right]}}=\mathsf{E}\left\\{\mathbf{w}^{a}_{\left[n\right]}\mathbf{w}^{a\mathsf{H}}_{\left[n\right]}\right\\}$
and
$\mathbf{M}_{n}=\mathsf{E}\left\\{\boldsymbol{\epsilon}^{a}_{\left[n\right]}\boldsymbol{\epsilon}^{a\mathsf{H}}_{\left[n\right]}\right\\}$
are positive semi-definite matrices. Furthermore, replacing (40) into (39)
gives the Riccati recursions
$\mathbf{M}_{\left[n+1\right]}=\left(\left(\mathbf{A}\mathbf{M}_{\left[n\right]}\mathbf{A}^{\mathsf{H}}+\boldsymbol{\Sigma}_{\mathbf{v}^{a}_{\left[n\right]}}\right)^{-1}+\mathbf{H}^{\mathsf{H}}\boldsymbol{\Sigma}^{-1}_{\mathbf{w}^{a}_{\left[n\right]}}\mathbf{H}\right)^{-1}.$
(41) Note that (39) and (41) are Lyapunov and Riccati recursions that express
how the elements of the error vector are expected to behave, that is, rotate
and scale, that is on par with what is presented in SM-1. Finally, if
$\\{\mathbf{A},\boldsymbol{\Sigma}^{\frac{1}{2}}_{\mathbf{v}^{a}_{\left[n\right]}}\\}$
(cf. $\\{\mathbf{A},\mathbf{H}\\}$) are stabilisable (cf. detectable) in the
isomorphic matrix algebra of quaternions (see SM-1 and [25]), for any positive
definite initial matrix $\mathbf{M}_{\left[1\right]}$, (41) has a unique
stabilising solution $\mathbf{M}$ [23], which can be used for formulating a
time-independent gain as,
$\mathbf{G}=\mathbf{M}\mathbf{H}^{\mathsf{H}}\boldsymbol{\Sigma}_{\mathbf{w}^{a}_{\left[n\right]}}$
that guarantees
$\mathsf{p}\left\\{\left(\mathbf{I}-\mathbf{G}_{\left[n\right]}\mathbf{H}\right)\mathbf{A}\right\\}<1$.
For a detailed treatment of quaternion-valued Kalman filter the reader is
referred to [23, 15].

###  III-B Sensor Fusion

Consider the set of sensors $\\{\mathsf{S}_{l}:l=1,\ldots,M\\}$, where each
sensor, $\mathsf{S}_{l}$, has estimate
$\hat{\mathbf{x}}_{l,\left[n\right]}=\mathsf{E}\left\\{\mathbf{x}_{\left[n\right]}|\mathbf{y}_{l,\left[1:n\right]}\right\\}$
with $\mathbf{y}_{l,\left[1:n\right]}$ repressing the observation of sensor
$l$ up to time instant $n$. The aim is to find
$\hat{\mathbf{x}}_{\left[n\right]}=\mathsf{E}\left\\{\mathbf{x}_{\left[n\right]}|\hat{\mathbf{x}}_{1:M,\left[1:n\right]}\right\\}$.
This is on par with the problem considered in Section III-A so that; i)
$\\{\mathrm{f}\left(\cdot\right),\mathrm{h}\left(\cdot\right)\\}$ become
identity functions, ii) $\boldsymbol{\psi}_{\left[n\right]}$ becomes a vector
incorporating any a priori information, and iii)
$\begin{bmatrix}\hat{\mathbf{x}}^{a\mathsf{T}}_{1,\left[n\right]}&\ldots&\hat{\mathbf{x}}^{a\mathsf{T}}_{m,\left[n\right]}\end{bmatrix}^{\mathsf{T}}$
becomes the available observation at time instant $n$, i.e.,
$\mathbf{y}_{\left[n\right]}$. Using the the same machinery as in Section
III-A, from (36) it follows that

$\hat{\mathbf{x}}^{a}_{\left[n+1\right]}=\psi^{a}_{\left[n+1\right]}-\mathbf{G}_{\left[n\right]}\nabla_{\psi^{a}_{\left[n+1\right]}}\mathrm{d}\left(\psi^{a}_{\left[n+1\right]}-\mathbf{y}^{a}_{\left[n\right]}\right).$
(42)

Furthermore, it is generally assumed that; i) sensor estimates are independent
with Gaussian distributions, ii) $\mathrm{d}\left(\cdot\right)=\|\cdot\|^{2}$.
Thus, through the framework presented in Example 3, (42) simplifies into

$\hat{\mathbf{x}}^{a}_{\left[n+1\right]}=\psi^{a}_{\left[n+1\right]}-\mathbf{G}_{\left[n\right]}\left(\sum^{M}_{l=1}\boldsymbol{\Sigma}^{-1}_{\left(\mathbf{x}^{a}_{\left[n\right]}-\hat{\mathbf{x}}^{a}_{l,\left[n\right]}\right)}\hat{\mathbf{x}}^{a}_{\left[n\right]}\right)$
(43)

where
$\boldsymbol{\Sigma}_{\left(\mathbf{x}^{a}_{\left[n\right]}-\hat{\mathbf{x}}^{a}_{l,\left[n\right]}\right)}=\mathsf{E}\left\\{\left(\mathbf{x}^{a}_{\left[n\right]}-\hat{\mathbf{x}}^{a}_{l,\left[n\right]}\right)\left(\mathbf{x}^{a}_{\left[n\right]}-\hat{\mathbf{x}}^{a}_{l,\left[n\right]}\right)^{\mathsf{H}}\right\\}$
and $\mathbf{G}_{\left[n\right]}$ is selected to minimise
$\boldsymbol{\Sigma}_{\left(\mathbf{x}^{a}_{\left[n\right]}-\hat{\mathbf{x}}^{a}_{\left[n\right]}\right)}$.
Although the required augmented covariance matrices,
$\\{\boldsymbol{\Sigma}_{\left(\mathbf{x}^{a}_{\left[n\right]}-\hat{\mathbf{x}}^{a}_{l,\left[n\right]}\right)}:l=1,\ldots,M\\}$,
are tractable is some cases (see Example 3), it is generally assumed that
$\psi^{a}_{\left[n+1\right]}\rightarrow 0$ and $\forall
l:\boldsymbol{\Sigma}_{\left(\mathbf{x}^{a}_{\left[n\right]}-\hat{\mathbf{x}}^{a}_{l,\left[n\right]}\right)}\rightarrow
w_{l}\mathbf{I}$, further simplifying (43) into

$\mathbf{x}^{a}_{\left[n+1\right]}=\sum^{M}_{n=1}w_{l}\mathbf{x}^{a}_{l,\left[n\right]}$
(44)

where we have assumed $\sum^{M}_{l=1}\mathbf{w}_{l}=1$ to enforce333Note that
since the combination weight $w_{l}$ is taking the role of a positive definite
covariance; then, $w_{l}\in\mathbb{R}^{+}$.
$\mathbf{G}_{\left[n\right]}\rightarrow\mathbf{I}$. Application of sensor
fusion in distributed estimation and federated learning is elaborated upon in
SM-7.

###  III-C Backpropagation in the Quaternion Domain

Consider a neural network of $L$ layers with $N_{l}$ neurons in the $l$th
layer. Output of the $n$th neuron on the $l$th layer can be formulated as

$y_{l,n}=\sum_{m=1}^{N_{l}}w_{l,n,m}x_{l-1,m}+b_{l,n}$

where $w_{l,n,m}$ are quaternion-valued weights and $b_{l,n}$ is the bias
term, while

$x_{l,n}=\mathrm{f}\left(y_{l,n}\right)$

is the output of neuron $n$ in layer $l$ with $\mathrm{f}\left(\cdot\right)$
represents a nonlinear quaternion-valued activation function. The aim is to
find the gradient of the cost function

$\mathrm{J}=\frac{1}{2}\sum_{m=1}^{N_{L}}\|d_{m}-x_{L,m}\|^{2}$ (45)

with respect to the quaternion-vlaued weights, where $d_{m}$ denoting a
desired output. Given the cost function $\mathrm{J}$ in (45), the gradient is
given by

$\nabla_{\mathbf{W}^{*}_{L}}\mathrm{J}=-\frac{1}{2}\sum_{\forall\zeta\in\\{1,\imath,\jmath,\kappa\\}}\left(\nabla_{\mathbf{W}^{\zeta}_{L}}\mathrm{J}^{\mathsf{H}}\boldsymbol{\epsilon}\right)^{\zeta}$

with

$\mathbf{W}_{L}=\begin{bmatrix}w_{L,1,1}&\dots&w_{L,1,N_{L-1}}\\\
\vdots&\ddots&\vdots\\\
w_{L,N_{L},1}&\dots&w_{L,N_{L},N_{L-1}}\end{bmatrix}\hskip
7.39755pt\text{and}\hskip
7.39755pt\boldsymbol{\epsilon}=\begin{bmatrix}d_{1}\\\ \dots\\\
d_{N_{L}}\end{bmatrix}-\underbrace{\begin{bmatrix}x_{L,1}\\\ \dots\\\
x_{L,N_{L}}\end{bmatrix}}_{\mathbf{x}_{L}}\cdot$

This gradient can be described in a layer-by-layer manner so that

$\boldsymbol{\delta}_{l}=\left\\{\begin{array}[]{ll}\boldsymbol{\epsilon}&l=L\\\
\left(\mathbf{W}_{l+1}\mathbf{x}_{l+1}\right)^{\mathsf{H}}\boldsymbol{\delta}_{l+1}&l\neq
L.\end{array}\right.$

The weight and bias update rules for the output layer then become

$\mathbf{W}_{l}\leftarrow\mathbf{W}_{l}+\gamma\boldsymbol{\delta}_{l}\mathbf{x}^{\mathsf{H}}_{l-1}\hskip
7.39755pt\text{and}\hskip
7.39755pt\begin{bmatrix}b_{l,1}&\vdots&b_{l,N_{L}}\end{bmatrix}^{\mathsf{T}}\leftarrow\begin{bmatrix}b_{l,1}&\vdots&b_{l,N_{L}}\end{bmatrix}^{\mathsf{T}}+\gamma\boldsymbol{\delta}_{l}$

while the weight and bias update rules for the hidden layers are given by

$\mathbf{W}_{l}\leftarrow\mathbf{W}_{l}+\gamma\sum_{\forall\zeta\in\\{1,\imath,\jmath,\kappa\\}}\left(\boldsymbol{\delta}_{l}\right)^{\zeta}\mathbf{x}^{\mathsf{H}}_{l-1}\hskip
7.39755pt\text{and}\hskip 7.39755pt\begin{bmatrix}b_{l,1}\\\ \vdots\\\
b_{l,N_{l}}\end{bmatrix}\leftarrow\begin{bmatrix}b_{l,1}\\\ \vdots\\\
b_{l,N_{l}}\end{bmatrix}+\gamma\sum_{\forall\zeta\in\\{1,i,j,k\\}}\left(\boldsymbol{\delta}_{l}\right)^{\zeta}$

where $\gamma\in\mathbb{R}^{+}$ is the learning rate.

###  III-D Linear Dynamic Programming and Initial Elements of Reinforcement
Learning

To revisit the fundamental problem at the heart of reinforcement learning [24]
and control [20] applications, consider the optimisation problem formulated
as:

$\displaystyle\text{minimise}:$ $\displaystyle\hskip
4.55254ptJ=\mathbf{x}^{a\mathsf{H}}_{\left[N\right]}\mathbf{T}^{a}\mathbf{x}^{a}_{\left[N\right]}+\sum^{N-1}_{n=1}\left(\mathbf{x}^{a\mathsf{H}}_{\left[n\right]}\mathbf{Q}\mathbf{x}^{a}_{\left[n\right]}+\mathbf{u}^{a\mathsf{H}}_{\left[n\right]}\mathbf{R}\mathbf{u}^{a}_{\left[n\right]}\right)$
(46) so that: $\displaystyle\hskip
4.55254pt\mathbf{u}_{\left[n\right]}=\mathbf{G}_{\left[n\right]}\mathbf{x}^{a}_{n}\hskip
7.39755pt\text{and}\hskip
7.39755pt\mathbf{x}^{a}_{\left[n+1\right]}=\mathbf{A}\mathbf{x}^{a}_{\left[n\right]}+\mathbf{B}\mathbf{u}^{a}_{\left[n\right]}$
(47)

where $N$ represents the time horizon for calculation of inputs
$\left\\{\mathbf{u}_{\left[n\right]}:n\in\mathbb{N}\right\\}$, matrices
$\\{\mathbf{T},\mathbf{R},\mathbf{Q}\\}$ are Hermitian transpose symmetric and
positive semi-definite, that is, $\forall\mathbf{z}\in\mathbb{H}\hskip
3.41418pt\&\hskip
3.41418pt\mathbf{M}\in\\{\mathbf{T},\mathbf{R},\mathbf{Q}\\}:\hskip
4.55254pt\mathbf{z}^{a\mathsf{H}}\mathbf{M}\mathbf{z}^{a}\in\mathbb{R}^{+}\cup\\{0\\}$,
while $\mathbf{G}_{\left[n\right]}$ is the feedback gain matrix (optimal
policy). The matrices $\mathbf{T}$, $\mathbf{Q}$, and $\mathbf{R}$ penalise
the final deviation from the desired state, deviation from the desisted state
during the process, and implementation of an input.

The goal of the problem is to find the optimal feedback gain matrix at each
time instant. To this end, the cost function in (46) is sectioned into time
specific cost functions given by

$J_{\left[n\right]}=\left\\{\begin{array}[]{ll}\mathbf{x}^{a\mathsf{H}}_{\left[N\right]}\mathbf{T}\mathbf{x}^{a}_{\left[N\right]}&\text{if
}n=N,\\\
\mathbf{x}^{a\mathsf{H}}_{\left[n\right]}\mathbf{Q}\mathbf{x}^{a}_{\left[n\right]}+\mathbf{u}^{a\mathsf{H}}_{\left[n\right]}\mathbf{R}\mathbf{u}^{a}_{\left[n\right]}&\text{otherwise.}\end{array}\right.$
(48)

Now, upon inserting (47) into (48) we have

$J_{\left[n\right]}=\mathbf{x}^{a\mathsf{H}}_{\left[n\right]}\mathbf{Q}\mathbf{x}^{a}_{\left[n\right]}+\mathbf{u}^{a\mathsf{H}}_{\left[n\right]}\mathbf{R}\mathbf{u}^{a}_{\left[n\right]}=\mathbf{x}^{a\mathsf{H}}_{\left[n\right]}\boldsymbol{\Psi}_{\left[n\right]}\mathbf{x}^{a}_{\left[n\right]}$
(49)

with
$\boldsymbol{\Psi}_{\left[n\right]}=\mathbf{Q}+\mathbf{G}^{\mathsf{H}}_{\left[n\right]}\mathbf{R}\mathbf{G}_{\left[n\right]}$.
The expressions in (48) and (49) allow the time accumulative cost function to
be defined as

$\mathscr{J}_{\left[n\right]}=\sum^{N}_{m=n}J_{\left[m\right]}.$ (50)

This indicates that there exists a positive definite Hermitian symmetric
matrix $\mathbf{P}_{\left[n\right]}$, the existence of which is proven in
SM-8, such that

$\mathscr{J}_{\left[n\right]}=\mathbf{x}^{a\mathsf{H}}_{\left[n\right]}\mathbf{P}_{\left[n\right]}\mathbf{x}^{a}_{\left[n\right]}$
(51)

which extends the Hamilton-Jacobi-Bellman (see [20]) for solving quaternion-
valued optimisation problems.

Moving backwards in time with the implicit assumption that inputs for future
time instances have been obtained, the task becomes that of calculating the
input at the previous time instant resulting in444All elements in (52) have
been scaled by a factor of $2$ in order to compensate for the scales generated
after differentiation.

$2\mathscr{J}_{\left[n\right]}=\underset{\mathbf{u}_{\left[n\right]}}{\text{min}}\left\\{2\mathbf{x}^{a\mathsf{H}}_{\left[n\right]}\mathbf{Q}\mathbf{x}^{a}_{\left[n\right]}+2\mathbf{u}^{\mathsf{H}}_{\left[n\right]}\mathbf{R}\mathbf{u}_{\left[n\right]}+2\mathscr{J}_{\left[n\right]}\right\\}.$
(52)

The minimum of (52) occurs at

$\mathbf{u}_{\left[n\right]}=-\left(\mathbf{R}+\mathbf{B}^{\mathsf{H}}\mathbf{P}_{\left[n+1\right]}\mathbf{B}\right)^{-1}\mathbf{B}^{\mathsf{H}}\mathbf{P}_{\left[n+1\right]}\mathbf{A}^{a}\mathbf{x}^{a}_{\left[n\right]}.$
(53)

Thus, from (53) and (47) it becomes clear that
$\mathbf{G}_{\left[n\right]}=-\left(\mathbf{R}+\mathbf{B}^{\mathsf{H}}\mathbf{P}_{\left[n+1\right]}\mathbf{B}\right)^{-1}\mathbf{B}^{\mathsf{H}}\mathbf{P}_{\left[n+1\right]}\mathbf{A}$.
Crucially, replacing (53) into (52) and after some mathematical manipulation,
we have

$\mathbf{S}_{\left[n\right]}=\left(\left(\mathbf{A}^{\mathsf{H}}\mathbf{S}_{\left[n+1\right]}\mathbf{A}+\mathbf{Q}\right)^{-1}+\mathbf{B}\mathbf{R}^{-1}\mathbf{B}^{\mathsf{H}}\right)^{-1}$
(54)

where
$\mathbf{P}_{\left[n\right]}=\mathbf{A}^{\mathsf{H}}\mathbf{S}_{\left[n+1\right]}\mathbf{A}+\mathbf{Q}$.
Note that (54) is a quaternion-valued algebraic Riccati recursion, see SM-1,
that converges to a unique solution as $n\rightarrow\infty$ given
$\\{\mathbf{A}^{\mathsf{H}},\mathbf{Q}^{\frac{1}{2}}\\}$ (cf.
$\\{\mathbf{A},\mathbf{B}^{\mathsf{H}}\\}$) are stabilisable (cf. detectable),
see Example 3 and [15, 23]. Finally, the inputs can be formulated as

$\mathbf{u}_{\left[n\right]}=-\mathbf{R}^{-1}\mathbf{B}^{\mathsf{H}}\mathbf{S}_{\left[n+1\right]}\mathbf{A}\mathbf{x}^{a}_{\left[n\right]}.$
(55)

##  IV Conclusion

The article has revisited the $\mathbb{HR}$-calculus and its application in
processing quaternion-vlaued signals. The main tools in this area, such as the
chain derivative rule and the multiplication derivate rule, have been derived
in a straightforward manner. These rules have been in turn used to derive a
stochastic gradient descent algorithm that forms the basis for adaptive
filtering techniques in the quaternion domain. For rigour, application of the
$\mathbb{HR}$-calculus in Kalman filtering and quaternion-vlaued linear
dynamic programming have also been explored. This enables the use of
quaternions in statistical learning and decision making applications, a number
of which are presented in the supplementary martial accompanying the
manuscript.

## References

  * [1] J. Voight,  Quaternion Algebras,  Springer, 2021.
  * [2] T. L. Hankins,  “Triplets and triads: Sir William Rowan Hamilton on the metaphysics of mathematics,”  Isis, vol. 68, no. 2, pp. 175–193, 1977.
  * [3] D. Finkelstein, J. M. Jauch, S. Schiminovich, and D. Speiser,  “Foundations of Quaternion Quantum Mechanics,”  Journal of Mathematical Physics, vol. 3, no. 2, pp. 207–220, December 2004.
  * [4] M. Danielewski and S. Lucjan,  “Foundations of the quaternion quantum mechanics,”  Entropy, vol. 22, no. 12, 2020.
  * [5] J. C. Maxwell,  A Treatise on Electricity and Magnetism,  Dover Publications, New York, 1873,  ISBN 0-486-60636-8 (Vol. 1) and 0-486-60637-6 (Vol. 2).
  * [6] S. P. Talebi and Danilo P. Mandic,  “A quaternion frequency estimator for three-phase power systems,”  In Proceedings of IEEE International Conference on Acoustics, Speech and Signal Processing, pp. 3956–3960, 2015.
  * [7] S. P. Talebi, S. Kanna, and D. P. Mandic,  “A distributed quaternion Kalman filter with applications to smart grid and target tracking,”  IEEE Transactions on Signal and Information Processing over Networks, vol. 2, no. 4, pp. 477–488, December 2016.
  * [8] F. A. Tobar and D. P. Mandic,  “Quaternion reproducing kernel Hilbert spaces: Existence and uniqueness conditions,”  IEEE Transactions on Information Theory, vol. 60, no. 9, pp. 5736–5749, 2014.
  * [9] D. Comminiello, M. Lella, S. Scardapane, and A. Uncini,  “Quaternion convolutional neural networks for detection and localization of 3D sound events,”  In Proceedings of IEEE International Conference on Acoustics, Speech and Signal Processing, pp. 8533–8537, 2019.
  * [10] S. Miron, N. Le Bihan, and J. I. Mars,  “Quaternion-MUSIC for vector-sensor array processing,”  IEEE Transactions on Signal Processing, vol. 54, no. 4, pp. 1218–1229, April 2006.
  * [11] T. A. Ell and S. J. Sangwine,  “Hypercomplex Fourier transforms of color images,”  IEEE Transactions on Image Processing, vol. 16, no. 1, pp. 22–35, January 2007.
  * [12] S. Stern and R. F. H. Fischer,  “Quaternion-valued multi-user MIMO transmission via dual-polarized antennas and QLLL reduction,”  In Proceedings of International Conference on Telecommunications, pp. 63–69, 2018.
  * [13] J. Seberry, K. Finlayson, S. S. Adams, T. A. Wysocki, T. Xia, and B. J. Wysocki,  “The theory of quaternion orthogonal designs,”  IEEE Transactions on Signal Processing, vol. 56, no. 1, pp. 256–265, january 2008.
  * [14] V. Kliuchnikov and J. T. Yard,  “A framework for exact synthesis,”  Tech. Rep., 2015,  arXiv preprint arXiv:1504.04350.
  * [15] S. P. Talebi,  Adaptive filtering algorithms for quaternion-valued signals,  Ph.D. thesis, Imperial College London, 2016.
  * [16] T. A. Ell and S. J. Sangwine,  “Quaternion involutions and anti-involutions,”  Computers & Mathematics with Applications, vol. 53, no. 1, pp. 137–143, 2007.
  * [17] C. C. Took and D. P. Mandic,  “Augmented second-order statistics of quaternion random signals,”  Signal Processing, vol. 91, no. 2, pp. 214–224, 2011.
  * [18] D. P. Mandic, C. Jahanchahi, and C. C. Took,  “A quaternion gradient operator and its applications,”  IEEE Signal Processing Letters, vol. 18, no. 1, pp. 47–50, January 2011.
  * [19] C. Jahanchahi,  Quaternion valued adaptive signal processing,  Ph.D. thesis, Imperial College London, 2013.
  * [20] R. F. Stengel,  Optimal Estimation and Control,  Dover Publications, 1994.
  * [21] C. Cheong Took, S. P. Talebi, R. M. Fernandez Alcala, and D. P. Mandic,  “Augmented statistics of quaternion random variables: A lynchpin of quaternion learning machines,”  IEEE Signal Processing Magazine, Submitted.
  * [22] D. Xu D, C. Jahanchahi, C. C. Took, and D. P. Mandic,  “Enabling quaternion derivatives: The generalized HR calculus,”  Royal Society Open Science, vol. 2, no. 8, August 2015.
  * [23] S. P. Talebi, S. Werner, and D. P. Mandic,  “Quaternion-valued distributed filtering and control,”  IEEE Transactions on Automatic Control, vol. 65, no. 10, pp. 4246–4257, 2020.
  * [24] Richard S. Sutton and Andrew G. Barto,  Reinforcement Learning: An Introduction,  MIT Press, second edition, 2018.
  * [25] D. Alfsmann, H. G. Gockler, S. J. Sangwine, and T. A. Ell,  “Hypercomplex algebras in digital signal processing: Benefits and drawbacks,”  In Proceedings of European Signal Processing Conference, pp. 1322–1326, September 2007.
  * [26] V. P. Brasil, A. L. F. Filho, and J. Y. Ishihara,  “Electrical three phase circuit analysis using quaternions,”  International Conference on Harmonics and Quality of Power, pp. 1–6, 2018.
  * [27] O. V. Nos,  “The quaternion model of doubly-fed induction motor,”  International Forum on Strategic Technology, pp. 32–36, 2016.
  * [28] J. L. Crassidis, F. L. Markley, and Y. Cheng,  “Survey of nonlinear attitude estimation methods,”  Journal of Guidance, Control, and Dynamics, vol. 30, no. 1, pp. 12–28, January 2007.

Danilo Mandic<EMAIL_ADDRESS>received the Ph.D. degree in nonlinear
adaptive signal processing in 1999 from Imperial College, London, where he is
now a Professor. He specialises in Statistical Learning Theory, Machine
Intelligence, and Statistical Signal Processing. Recently, he was invited as a
distinguished lecturer by the IEEE Signal Processing Society to talk on a wide
range of topics ranging from graph theory to innovation in education.

Sayed Pouria Talebi<EMAIL_ADDRESS>received his doctorate
degree from Imperial College London, where is research focused on quaternion-
valued signal processing. He has served as a researcher at a number of
institutions including University of Cambridge. He is currently a lecturer at
University of Roehampton. Some of his contributions to this area include
formulation of the theory of quaternion-valued multi-agent systems,
formulation of the basis of quaternion-valued control theory, and quaternion-
valued nonlinear/non-Gaussian signal processing. In his work, he has
demonstrated the usefulness of quaternions for modelling power distributions
systems and passive target tracking.

Clive Cheong Took<EMAIL_ADDRESS>received his doctorate from
Cardiff University. He is a senior lecturer at Royal Holloway, University of
London, U.K. He has been working on quaternion signal processing since 2008
and delivered a tutorial on this topic in the conference EUSIPCO 2011. Some of
his contributions to this area include augmented statistics, adaptive
learning, and novel matrix factorisations. Over the years, Clive has
demonstrated the usefulness of quaternions to address various problems such as
source separation, prediction or estimation, and classification for various
biomedical applications. Recently, his research in quaternions has nurtured
his interest in three-dimensional rotations for which quaternions are popular
for. He is a Senior Member of IEEE.

Yili Xia<EMAIL_ADDRESS>received the Ph.D. degree in adaptive signal
processing from Imperial College London in 2011. Since 2013, he has been an
Associate Professor in signal processing with the School of Information
Science and Engineering, Southeast University, where he is currently the
Deputy Head of the Department of Information and Signal Processing
Engineering. His research interests include complex and hyper-complex
statistical analysis. He was a recipient of the Best Student Paper Award at
the International Symposium on Neural Networks (ISNN) in 2010 (coauthor) and
the Education Innovation Award at the IEEE International Conference on
Acoustics, Speech, and Signal Processing (ICASSP) in 2019. He is an Associate
Editor for the IEEE Transactions on Signal Processing.

Dongpo Xu<EMAIL_ADDRESS>received the Ph.D. degree in computational
mathematics from the Dalian University of Technology, Dalian, China in 2009.
He was a Lecturer with the College of Science, Harbin Engineering University
and a Visiting Scholar with the Department of Electrical and Electronic
Engineering, Imperial College London. He is currently a Professor with the
School of Mathematics and Statistics, Northeast Normal University, Changchun,
China. His current research interests include neural networks, machine
learning, and signal processing. Dr. Xu has been an Associate Editor of the
Frontiers in Artificial Intelligence and the Frontiers in Big Data.

Min Xiang<EMAIL_ADDRESS>received the Ph.D. degree in electrical
engineering from Imperial College London, London, U.K., in 2018. He was a
Research Associate with Imperial College London from 2018 to 2020. Since 2020,
he has been an Associate Professor with Beihang University, Beijing, China.
His research interests include statistical signal processing, specifically
quaternion-valued signal processing.

Pauline Bourigault<EMAIL_ADDRESS>is currently pursuing her
doctoral degree in AI and Machine Learning enrolled in the UKRI Centre for
Doctoral Training in AI for Healthcare at Imperial College London. Her current
research focuses on 3D human pose estimation.

## Supplementary Material

### SM-1: Quaternion Rotations and Algebraic Lyapunov/Riccati Recursions

Consider the evolution of a quaternion-valued sequence
$\\{q_{\left[n\right]}:n\in\mathbb{N}\\}$ that encapsulates all three
transformations (rotation, scaling, and shift) of a real-world systems
operating in the three-dimensional space, e.g., robot arm, aircraft/vehicle,
or an particle. Based on the quaternion model of rotation given in (6), the
evolutions of $q_{\left[n\right]}$ are formulated as

$q_{\left[n+1\right]}=a_{\left[n\right]}\mu_{\left[n\right]}q_{\left[n\right]}\mu^{-1}_{\left[n\right]}+q_{\text{shift}_{\left[n\right]}}=\sqrt{a_{\left[n\right]}}\mu_{\left[n\right]}q_{\left[n\right]}\sqrt{a_{\left[n\right]}}\mu^{-1}+q_{\text{shift}_{\left[n\right]}}$
(56)

where $a_{\left[n\right]}\in\mathbb{R}$ and
$q_{\text{shift}_{\left[n\right]}}$ represent the scaling and shift factors.
Moreover, as an associative hypercomplex algebra, quaternions have an
isomorphic matrix algebra [25], so that

$\forall
q\in\mathbb{H}\leftrightarrow\mathbf{Q}=\begin{bmatrix}q_{r}&-q_{\imath}&-q_{\jmath}&q_{\kappa}\\\
q_{\imath}&q_{r}&-q_{\kappa}&-q_{\jmath}\\\
q_{\jmath}&q_{\kappa}&q_{r}&q_{\imath}\\\
-q_{\kappa}&q_{\jmath}&-q_{\imath}&q_{r}\end{bmatrix}\in\mathbb{R}^{4}\times\mathbb{R}^{4}$
(57)

where $\mathbf{Q}$ is the representation of $q$ in the isomorphic matrix
algebra [25]. Thus, using the duality in (57) the expression in (56) can be
transformed into the real-valued Lyapunov style recursion

$\mathbf{Q}_{\left[n+1\right]}=\mathbf{M}_{\left[n\right]}\mathbf{Q}_{\left[n+1\right]}\mathbf{M}^{\mathsf{H}}_{\left[n\right]}+\mathbf{Q}_{\text{shift}_{\left[n\right]}}$

with $\mathbf{Q}_{\left[n\right]}$, $\mathbf{M}_{\left[n\right]}$, and
$\mathbf{Q}_{\text{shift}_{\left[n\right]}}$ representing the real-valued
matrix duals of $q_{\left[n\right]}$, $\sqrt{a}\mu$, and
$q_{\text{shift}_{\left[n\right]}}$. In addition, consider the Riccati
recursions in its information matrix formulation given as

$\mathbf{Q}_{\left[n+1\right]}=\underbrace{\left(\underbrace{\left(\mathbf{M}_{\left[n\right]}\mathbf{Q}_{\left[n\right]}\mathbf{M}^{\mathsf{H}}_{\left[n\right]}+\mathbf{Q}_{\text{shift}_{\left[n\right]}}\right)^{-1}}_{\text{inner
Lyapunov
operator}}+\mathbf{H}^{\mathsf{H}}_{\left[n\right]}\mathbf{R}^{-1}_{\left[n\right]}\mathbf{H}_{\left[n\right]}\right)^{-1}}_{\text{outer
Lyapunov operator}}$ (58)

for general matrices $\mathbf{H}_{\left[n\right]}$ and
$\mathbf{R}_{\left[n\right]}$. Then, it becomes clear that (58) consists of an
inner and an outer Lyapunov style operators that interact, which essentially
models two interacting systems of the from in (56). These dualities highlight
an important issue. The Lyapunov and Riccati recursions used to model the
involutions of physical systems are essentially representing the rotation,
scaling, and shift, that these systems experience from one time instant to the
next, which are more effectively dealt with using quaternions. Finally, the
example given for a scalar quaternion-valued random variable is extendable to
higher dimensional systems. For more details on the use and behaviour of
quaternion-valued Lyapunov and Riccati recursions, the reader is referred to
[15, 23].

## SM-2: Quaternion-Valued Minimum Mean Square Error (MMSE) Estimator

Consider the MMSE estimator of a variable, $y$, conditioned on the
observation, $z$, that is, $\hat{y}=\mathsf{E}\left\\{y|z\right\\}$. For the
case where $y$ and $z$ are quaternion-valued zero-mean and jointly Gaussian
random variables, the MMSE estimator has to be expressed according to the
real-valued components of the quaternion-valued random variables so that

$\displaystyle\hat{y}=$
$\displaystyle\mathsf{E}\left\\{y_{r}|z_{r},z_{\imath},z_{\jmath},z_{\kappa}\right\\}+\imath\mathsf{E}\left\\{y_{\imath}|z_{r},z_{\imath},z_{\jmath},z_{\kappa}\right\\}+$
$\displaystyle+\jmath\mathsf{E}\left\\{y_{\jmath}|z_{r},z_{\imath},z_{\jmath},z_{\kappa}\right\\}+\kappa\mathsf{E}\left\\{y_{\kappa}|z_{r},z_{\imath},z_{\jmath},z_{\kappa}\right\\}$

whereby using (4) to replace
$\\{z_{r}$,$z_{\imath}$,$z_{\jmath}$,$z_{\kappa}\\}$ with
$\\{z,z^{\imath},z^{\jmath},z^{\kappa}\\}$ gives

$\displaystyle\hat{y}=$
$\displaystyle\mathsf{E}\left\\{y_{r}|z,z^{\imath},z^{\jmath},z^{\kappa}\right\\}+\imath\mathsf{E}\left\\{y_{\imath}|z,z^{\imath},z^{\jmath},z^{\kappa}\right\\}$
$\displaystyle+\jmath\mathsf{E}\left\\{y_{\jmath}|z,z^{\imath},z^{\jmath},z^{\kappa}\right\\}+\kappa\mathsf{E}\left\\{y_{\kappa}|z,z^{\imath},z^{\jmath},z^{\kappa}\right\\}.$

Therefore, the MMSE solution in this case will take the form of an augmented
estimator given by [17]

$\hat{y}=\mathbf{g}^{\mathsf{T}}\mathbf{z}+\mathbf{h}^{\mathsf{T}}\mathbf{z}^{\imath}+\mathbf{u}^{\mathsf{T}}\mathbf{z}^{\jmath}+\mathbf{v}^{\mathsf{T}}\mathbf{z}^{\kappa}=\begin{bmatrix}\mathbf{g}^{\mathsf{T}},\mathbf{h}^{\mathsf{T}},\mathbf{u}^{\mathsf{T}},\mathbf{v}^{\mathsf{T}}\end{bmatrix}\mathbf{z}^{a}$
(59)

where $\\{\mathbf{g},\mathbf{h},\mathbf{u},\mathbf{v}\\}$ are quaternion-
valued coefficient vectors and $\mathbf{z}$ is the regressor vector. The
estimator in (59) can also be expressed as

$\hat{\mathbf{y}}^{a}=\mathbf{W}\mathbf{z}^{a}$ (60)

where $\mathbf{W}$ is the matrix constructed appropriately from
$\\{\mathbf{g},\mathbf{h},\mathbf{u},\mathbf{v}\\}$ and their involutions
[17].

### SM-3: Relation Between the AQCF and Statistical Moments

Consider the AQCF in the formulation given in (9). Substituting the
exponential with its power series representation yields

$\displaystyle\boldsymbol{\Phi}_{\textbf{{Q}}^{a}}(\mathbf{s}^{a})=\mathsf{E}\left\\{e^{\left(\frac{\xi}{4}\mathbf{s}^{a\mathsf{H}}\mathbf{q}^{a}\right)}\right\\}=$
$\displaystyle\mathsf{E}\left\\{\sum^{\infty}_{m=0}\frac{\xi^{m}\left(\mathbf{s}^{a\mathsf{H}}\mathbf{q}^{a}\right)^{m}}{(4^{m})m!}\right\\}$
(61) $\displaystyle=$ $\displaystyle
1+\mathsf{E}\left\\{\xi\frac{\mathbf{s}^{a\mathsf{H}}\mathbf{q}^{a}}{4}\right\\}-\mathsf{E}\left\\{\frac{\left(\mathbf{s}^{a\mathsf{H}}\mathbf{q}^{a}\right)^{2}}{32}\right\\}\cdots$

Using the $\mathbb{HR}$-calculus and derivatives in (14) one can verify from
(61) that

$\displaystyle\frac{\partial\boldsymbol{\Phi}_{\textbf{{Q}}^{a}}(\mathbf{s}^{a})}{\partial\mathbf{s}^{*\zeta}}=$
$\displaystyle\frac{\partial\boldsymbol{\Phi}_{\textbf{{Q}}^{a}}(\mathbf{s}^{a})}{\partial\mathbf{s}^{a\mathsf{H}}\mathbf{q}^{a}}\frac{\partial\mathbf{s}^{a\mathsf{H}}\mathbf{q}^{a}}{\partial\mathbf{s}^{*\zeta}}=\mathsf{E}\left\\{\sum^{\infty}_{m=0}\frac{\xi^{m+1}\left(\mathbf{s}^{a\mathsf{H}}\mathbf{q}^{a}\right)^{m}}{(4^{m+1})m!}\mathbf{q}^{\zeta}\right\\}$
(62) $\displaystyle=$
$\displaystyle\mathsf{E}\left\\{\xi\frac{\mathbf{q}^{\zeta}}{4}\right\\}-E\left[\frac{\mathbf{s}^{a\mathsf{H}}\mathbf{q}^{a}}{16}\mathbf{q}^{\zeta}\right]-\mathsf{E}\left\\{\xi\frac{\left(\mathbf{s}^{a\mathsf{H}}\mathbf{q}^{a}\right)^{2}}{128}\mathbf{q}^{\zeta}\right\\}\cdots$

with $\zeta\in\\{1,\imath,\jmath,\kappa\\}$. Evaluating (62) at
$\mathbf{s}^{a}=0$ gives

$\left.\frac{\partial\boldsymbol{\Phi}_{\textbf{{Q}}^{a}}(\mathbf{s}^{a})}{\partial\mathbf{s}^{*\zeta}}\right|_{\mathbf{s}^{a}=0}=\mathsf{E}\left\\{\xi\frac{\mathbf{q}^{\zeta}}{4}\right\\}$

and hence we have

$\mathsf{E}\left\\{\mathbf{q}\right\\}=\left(\frac{4}{\xi}\left.\frac{\partial\boldsymbol{\Phi}_{\textbf{{Q}}^{a}}(\mathbf{s}^{a})}{\partial\mathbf{s}^{*\zeta}}\right|_{\mathbf{s}^{a}=0}\right)^{\zeta}.$

The second-order derivatives of $\mathbf{\Phi}_{\textbf{{Q}}^{a}}$ can be
expressed as

$\frac{\partial^{2}\mathbf{\Phi}_{\textbf{{Q}}^{a}}(\mathbf{s}^{a})}{\partial\mathbf{s}^{*\zeta_{1}}\partial\mathbf{s}^{*\zeta_{2}}}=\frac{\partial}{\partial\mathbf{s}^{*\zeta_{1}}}\left(\frac{\partial\mathbf{\Phi}_{\textbf{{Q}}^{a}}(\mathbf{s}^{a})}{\partial\mathbf{s}^{*\zeta_{2}}}\right)=\frac{\partial}{\partial\mathbf{s}^{a\mathsf{H}}\mathbf{q}^{a}}\left(\frac{\partial\mathbf{\Phi}_{\textbf{{Q}}^{a}}(\mathbf{s}^{a})}{\partial\mathbf{s}^{a\mathsf{H}}\mathbf{q}^{a}}\frac{\partial\mathbf{s}^{a\mathsf{H}}\mathbf{q}^{a}}{\partial\mathbf{s}^{*\zeta_{2}}}\right)\frac{\partial\mathbf{s}^{a\mathsf{H}}\mathbf{q}^{a}}{\partial\mathbf{s}^{*\zeta_{1}}}$
(63)

where $\zeta_{1},\zeta_{2}\in\\{1,\imath,\jmath,\kappa\\}$. Upon replacing
(62) into (63) we have

$\frac{\partial^{2}\mathbf{\Phi}_{\textbf{{Q}}^{a}}(\mathbf{s}^{a})}{\partial\mathbf{s}^{*\zeta_{1}}\partial\mathbf{s}^{*\zeta_{2}}}=\mathsf{E}\left\\{\sum^{\infty}_{m=0}\frac{-\xi^{m}\left(\mathbf{s}^{a\mathsf{H}}\mathbf{q}^{a}\right)^{m}}{(4^{m+2})m!}\mathbf{q}^{\zeta_{2}}\left(\mathbf{q}^{\zeta_{1}}\right)^{\mathsf{T}}\right\\}$

that evaluated at $\mathbf{s}^{a}=0$ yields

$\left.\frac{\partial^{2}\mathbf{\Phi}_{\textbf{{Q}}^{a}}(\mathbf{s}^{a})}{\partial\mathbf{s}^{*\zeta_{1}}\partial\mathbf{s}^{*\zeta_{2}}}\right|_{\mathbf{s}^{a}=0}=\frac{-1}{16}\mathsf{E}\left\\{\mathbf{q}^{\zeta_{2}}\left(\mathbf{q}^{\zeta_{1}}\right)^{\mathsf{T}}\right\\}$

giving the cross-correlation between $\mathbf{q}^{\zeta_{1}}$ and
$\mathbf{q}^{\zeta_{2}}$ multiplied by the factor $(\xi/4)^{2}$. Note that all
higher-order statistical moments of Q are obtainable in an analogous manner.

### SM-4: Real-Valued Components of (20)

Given
$\\{\mathrm{f}\left(\cdot\right),\mathrm{g}\left(\cdot\right)\\}:\mathbb{H}^{M}\rightarrow\mathbb{H}$
the real-valued components of
$\mathrm{h}\left(\mathbf{q}\right)=\mathrm{f}\left(\mathbf{q}\right)\mathrm{g}\left(\mathbf{q}\right)$
are

$\displaystyle\Re\left\\{\mathrm{h}\left(\mathbf{q}\right)\right\\}=$
$\displaystyle\left(\Re\left\\{\mathrm{f}\left(\mathbf{q}\right)\right\\}\right)\left(\Re\left\\{\mathrm{g}\left(\mathbf{q}\right)\right\\}\right)+\left(\Im_{\imath}\left\\{\mathrm{f}\left(\mathbf{q}\right)\right\\}\right)\left(\Im_{\jmath}\left\\{\mathrm{g}\left(\mathbf{q}\right)\right\\}\right)$
$\displaystyle+\left(\Im_{\jmath}\left\\{\mathrm{f}\left(\mathbf{q}\right)\right\\}\right)\left(\Im_{\jmath}\left\\{\mathrm{g}\left(\mathbf{q}\right)\right\\}\right)+\left(\Im_{\kappa}\left\\{\mathrm{f}\left(\mathbf{q}\right)\right\\}\right)\left(\Im_{\kappa}\left\\{\mathrm{g}\left(\mathbf{q}\right)\right\\}\right)$
$\displaystyle\Im_{\imath}\left\\{\mathrm{h}\left(\mathbf{q}\right)\right\\}=$
$\displaystyle\left(\Re\left\\{\mathrm{f}\left(\mathbf{q}\right)\right\\}\right)\left(\Im_{\imath}\left\\{\mathrm{g}\left(\mathbf{q}\right)\right\\}\right)+\left(\Im_{\imath}\left\\{\mathrm{f}\left(\mathbf{q}\right)\right\\}\right)\left(\Re\left\\{\mathrm{g}\left(\mathbf{q}\right)\right\\}\right)$
$\displaystyle+\left(\Im_{\jmath}\left\\{\mathrm{f}\left(\mathbf{q}\right)\right\\}\right)\left(\Im_{\kappa}\left\\{\mathrm{g}\left(\mathbf{q}\right)\right\\}\right)+\left(\Im_{\kappa}\left\\{\mathrm{f}\left(\mathbf{q}\right)\right\\}\right)\left(\Im_{\jmath}\left\\{\mathrm{g}\left(\mathbf{q}\right)\right\\}\right)$
$\displaystyle\Im_{\jmath}\left\\{\mathrm{h}\left(\mathbf{q}\right)\right\\}=$
$\displaystyle\left(\Re\left\\{\mathrm{f}\left(\mathbf{q}\right)\right\\}\right)\left(\Im_{\jmath}\left\\{\mathrm{g}\left(\mathbf{q}\right)\right\\}\right)+\left(\Im_{\imath}\left\\{\mathrm{f}\left(\mathbf{q}\right)\right\\}\right)\left(\Im_{\kappa}\left\\{\mathrm{g}\left(\mathbf{q}\right)\right\\}\right)$
$\displaystyle+\left(\Im_{\jmath}\left\\{\mathrm{f}\left(\mathbf{q}\right)\right\\}\right)\left(\Re\left\\{\mathrm{g}\left(\mathbf{q}\right)\right\\}\right)+\left(\Im_{\kappa}\left\\{\mathrm{f}\left(\mathbf{q}\right)\right\\}\right)\left(\Im_{\imath}\left\\{\mathrm{g}\left(\mathbf{q}\right)\right\\}\right)$
$\displaystyle\Im_{\kappa}\left\\{\mathrm{h}\left(\mathbf{q}\right)\right\\}=$
$\displaystyle\left(\Re\left\\{\mathrm{f}\left(\mathbf{q}\right)\right\\}\right)\left(\Im_{\kappa}\left\\{\mathrm{g}\left(\mathbf{q}\right)\right\\}\right)+\left(\Im_{\imath}\left\\{\mathrm{f}\left(\mathbf{q}\right)\right\\}\right)\left(\Im_{\jmath}\left\\{\mathrm{g}\left(\mathbf{q}\right)\right\\}\right)$
$\displaystyle+\left(\Im_{\jmath}\left\\{\mathrm{f}\left(\mathbf{q}\right)\right\\}\right)\left(\Im_{\kappa}\left\\{\mathrm{g}\left(\mathbf{q}\right)\right\\}\right)+\left(\Im_{\kappa}\left\\{\mathrm{f}\left(\mathbf{q}\right)\right\\}\right)\left(\Re\left\\{\mathrm{g}\left(\mathbf{q}\right)\right\\}\right).$

### SM-5:Quaternion-Valued Least Mean Square Algorithm

Based on the regressor in (59), the aim is to find the weight vector
$\mathbf{w}$ so that the estimate
$\hat{y}_{n}=\mathbf{w}^{\mathsf{T}}_{\left[n\right]}\mathbf{z}^{a}_{\left[n\right]}$
tracks the observation sequence $\\{y_{\left[n\right]}:n\in\mathbb{N}\\}$ in a
manner that minimises the cost

$J_{\left[n\right]}=\|\epsilon_{\left[n\right]}\|^{2}\hskip
7.39755pt\text{with}\hskip
7.39755pt\epsilon_{\left[n\right]}=y_{\left[n\right]}-\hat{y}_{\left[n\right]}.$

This is performed in the lest mean square (LMS) fashion via iterative updates
of the weight vector as

$\mathbf{w}_{\left[n+1\right]}=\mathbf{w}_{\left[n\right]}-\gamma\nabla_{\mathbf{w}^{*}}J_{\left[n\right]}$

where $\gamma\in\mathbb{R}^{+}$ is an adaptation gain. Now, considering that
$\|\epsilon_{\left[n\right]}\|\in\mathbb{R}$ and using the chain derivative
rule of the $\mathbb{HR}$-calculus given in (27), we have

$\nabla_{\mathbf{w}^{*}}J_{\left[n\right]}=\frac{\partial\|\epsilon_{\left[n\right]}\|^{2}}{\partial\|\epsilon_{\left[n\right]}\|}\nabla_{\mathbf{w}^{*}}\|\epsilon_{\left[n\right]}\|=-\frac{1}{2}\epsilon_{\left[n\right]}\mathbf{z}^{a*}_{\left[n\right]}.$

Thus, the weight update term becomes

$\mathbf{w}_{\left[n+1\right]}=\mathbf{w}_{\left[n\right]}+\gamma\epsilon_{\left[n\right]}\mathbf{z}^{a*}_{\left[n\right]}$
(64)

with the $\frac{1}{2}$ factor absorbed into the adaptation gain.

######  Remark 8.

Using the multiplication rule in (21), we can also have

$\displaystyle\nabla_{\mathbf{w}^{*}}J_{\left[n\right]}=$
$\displaystyle\nabla_{\mathbf{w}^{*}}\|\epsilon_{\left[n\right]}\|^{2}=\nabla_{\mathbf{w}^{*}}\left(\epsilon^{*}_{\left[n\right]}\epsilon_{\left[n\right]}\right)=\epsilon^{*}_{\left[n\right]}\left(\nabla_{\mathbf{w}^{*}}\epsilon_{\left[n\right]}\right)+\left(\nabla_{\mathbf{w}^{\epsilon_{\left[n\right]*}}}\epsilon^{*}_{\left[n\right]}\right)\epsilon_{\left[n\right]}$
(65) $\displaystyle=$
$\displaystyle\frac{1}{2}\epsilon^{*}_{\left[n\right]}\mathbf{z}^{a*}_{\left[n\right]}-\mathbf{z}^{a*}_{\left[n\right]}\Re\left\\{\epsilon_{\left[n\right]}\right\\}=-\frac{1}{2}\epsilon_{\left[n\right]}\mathbf{z}^{a*}_{\left[n\right]}.$

There are different derivations available, e.g., see [22, 19, 18]. Although
these derivation might seem different from those in (64) and (65), they
operate in the same manner.

Finally, if we assume that $\mathbf{w}^{\text{opt}}$ exists so that
$y_{\left[n\right]}=\mathbf{w}^{\text{opt}^{\mathsf{T}}}\mathbf{z}^{a}_{\left[n\right]}\hskip
7.39755pt\text{and}\hskip
7.39755pt\boldsymbol{\varepsilon}_{\left[n\right]}=\mathbf{w}^{\text{opt}}-\mathbf{w}_{\left[n\right]}$.
Then, from (64) it follows that
$\boldsymbol{\varepsilon}^{\mathsf{H}}_{\left[n+1\right]}=\boldsymbol{\varepsilon}^{\mathsf{H}}_{\left[n\right]}-\gamma\mathbf{z}^{a}_{\left[n\right]}\mathbf{z}^{a\mathsf{H}}_{\left[n\right]}\boldsymbol{\varepsilon}^{\mathsf{H}}_{\left[n\right]}$
which in turn yields the quaternion-valued Lyapunov function

$\boldsymbol{\Sigma}_{\boldsymbol{\varepsilon}^{a}_{\left[n+1\right]}}=\left(\mathbf{I}-\gamma\mathbf{G}_{\left[n\right]}\right)\boldsymbol{\Sigma}_{\boldsymbol{\varepsilon}_{\left[n\right]}}\left(\mathbf{I}-\gamma\mathbf{G}_{\left[n\right]}\right)^{\mathsf{H}}$

where
$\boldsymbol{\Sigma}_{\boldsymbol{\varepsilon}^{a}_{\left[n\right]}}=\mathsf{E}\left\\{\boldsymbol{\varepsilon}^{a}_{\left[n\right]}\boldsymbol{\varepsilon}^{a\mathsf{H}}_{\left[n\right]}\right\\}$,
while

$\boldsymbol{\Sigma}_{\mathbf{z}^{a}_{\left[n\right]}}=\mathsf{E}\left\\{\mathbf{z}^{a}_{\left[n\right]}\mathbf{z}^{a\mathsf{H}}_{\left[n\right]}\right\\}\hskip
7.39755pt\text{and}\hskip
7.39755pt\mathbf{G}_{\left[n\right]}=\text{diag}\left\\{\boldsymbol{\Sigma}_{\mathbf{z}^{a}_{\left[n\right]}},\boldsymbol{\Sigma}^{\imath}_{\mathbf{z}^{a}_{\left[n\right]}},\boldsymbol{\Sigma}^{\jmath}_{\mathbf{z}^{a}_{\left[n\right]}},\boldsymbol{\Sigma}^{\kappa}_{\mathbf{z}^{a}_{\left[n\right]}}\right\\}.$

Thus, if $\gamma<\mathsf{p}\left\\{\mathbf{G}_{\left[n\right]}\right\\}$ and
$\boldsymbol{\Sigma}_{\mathbf{z}^{a}_{\left[n\right]}}$ is full rank; then, as
$n\rightarrow\infty$ one has
$\mathsf{p}\left\\{\boldsymbol{\Sigma}_{\mathbf{z}^{a}_{\left[n\right]}}\right\\}\rightarrow
0$, resulting in a convergent filtering operation.

### SM-6: Quaternion-Valued Nonlinear Adaptive Filtering and Learning

Consider the operations of a computational element (neuron), formulated as

$\hat{y}_{n}=\mathrm{h}\left(\mathbf{w}^{\mathsf{T}}_{\left[n\right]}\mathbf{z}^{a}_{\left[n\right]}\right)$

where $\mathbf{z}_{\left[n\right]}$ is a vector presenting the inputs,
$\mathbf{w}_{\left[n\right]}$ is a vector of the corresponding weights, and
$\mathrm{h}\left(\cdot\right):\mathbb{H}\rightarrow\mathbb{H}$ is a
differentiable function in $\mathbb{H}$. Note that for mathematical
presentation the bias element is simply considered as an element of
$\mathbf{z}_{\left[n\right]}$. The aim is to find optimal wight vector,
$\mathbf{w}^{\text{opt}}$, so that the output, $\hat{y}_{\left[n\right]}$, is
a close fit to desired output, $y_{\left[n\right]}$, in a manner that
minimises a selected metric function $\mathrm{d}\left(\cdot\right)$ so that
$\left\\{\mathbb{H},\mathrm{d}\left(\cdot\right)\right\\}$ constitutes a
convex metric space.

The weight vector is updated iteratively in order to reduce the cost

$\mathrm{g}\left(\mathbf{w}_{\left[n\right]}\right)=\mathrm{d}\left(\hat{y}_{\left[n\right]}-y_{\left[n\right]}\right).$

Specifically, we would like to achieve

$\mathrm{g}\left(\mathbf{w}_{\left[n+1\right]}\right)<\mathrm{g}\left(\mathbf{w}_{\left[n\right]}\right)$
(66)

with
$\mathbf{w}_{\left[n+1\right]}=\mathbf{w}_{\left[n\right]}+\Delta\mathbf{w}_{\left[n\right]}$
where $\Delta\mathbf{w}_{\left[n\right]}$ represents a change in the weight
vector form time instant $n$ to time instant $n+1$. Taking the quaternion-
valued Taylor series expansion of (66) yields

$\mathrm{g}\left(\mathbf{w}_{\left[n\right]}\right)+\left(\Delta\mathbf{w}^{a}_{\left[n\right]}\right)^{\mathsf{H}}\nabla_{\mathbf{w}^{a*}_{\left[n\right]}}\mathrm{g}(\mathbf{w}_{\left[n\right]})+\mathcal{O}\left(\mathrm{g}\left(\mathbf{w}_{n}\right)\right)<\mathrm{g}\left(\mathbf{w}_{n}\right).$
(67)

Neglecting the higher order terms in (67) gives
$\left(\Delta\mathbf{w}^{a}_{\left[n\right]}\right)^{\mathsf{H}}\nabla_{\mathbf{w}^{a*}_{\left[n\right]}}\mathrm{g}(\mathbf{w}_{\left[n\right]})<0$.
Therefore, considering
$\Delta\mathbf{w}^{a}_{\left[n\right]}\propto\mu\nabla_{\mathbf{w}^{a*}_{\left[n\right]}}\mathrm{g}(\mathbf{w}_{\left[n\right]})$
allows the formulation of a weight update procedure as

$\mathbf{w}^{a}_{\left[n+1\right]}=\mathbf{w}^{a}_{\left[n\right]}-\gamma\nabla_{\mathbf{w}^{a*}_{\left[n\right]}}\mathrm{g}\left(\mathbf{w}_{\left[n\right]}\right)$
(68)

where $\gamma\in\mathbb{R}^{{}^{+}}$ is an adaptation gain.

######  Remark 9.

For the special case of $\mathrm{d}\left(\cdot\right)=\|\cdot\|^{2}_{2}$ and
using the Taylor series expansion of
$\mathrm{g}\left(\mathbf{w}_{\left[n+1\right]}\right)$ around
$\mathbf{w}_{\left[n\right]}$ suggests that
$\gamma\in\left(0,\frac{1}{8\|\mathbf{z}^{a}_{\left[n\right]}\|^{2}_{2}\|\nabla_{\mathbf{w}^{a*}_{\left[n\right]}}\mathrm{h}\left(\mathbf{w}^{T}\mathbf{z}^{a}_{\left[n\right]}\right)\|^{2}_{2}}\right)$
ensures convergence.

### SM-7: Distributed and Federated Learning

1-Distributed Estimation

In distributed estimation settings, multiple agents cooperate over an ad-hoc
communication network in order to collaboratively estimate/track a state
vector. In this setting, each agent implements a two stage learning process.
First, a local estimate is formed using an agents own observations based on
(36). Second, the agents fuse their local estimates with those of their
neighbours using (44). This allows agents to learn from observations of other
agents in the network, without the need for these observation to be shared.
This scenario is shown in Fig. 3.

Figure 3: Distributed learning scenario with a typical network shown on the
left hand side and operations of an agent shown on the right hand side. The
agent formulates and shares intermediate estimates
$\boldsymbol{\psi}_{l,\left[n\right]}$, while its fusion process produces
estimates $\hat{\mathbf{x}}_{l,\left[n\right]}$.

2-Federated Learning

One of the most popular multi-agent scenarios for learning from data sets that
are distributed among a multitude of agents is federated learning. In its most
fundamental level, at each time instant, agents (refereed to commonly as edge
devices) that have new data sets available, update their estimate using (36),
and then, send these updated estimates to a fusion centre which fuses these
estimates with its existing information using (43). Then, the fusion centre
pushes these unified estimates of the learnt parameters back to the agents.
This scenario is shown in Fig. 4.

Figure 4: The federated learning scenario is shown with a mix of active and
inactive edge devices.

### SM-8: Quaternion-Valued Hamilton-Jacobi-Bellman Type Recursion

The aim is to recover matrix sequence
$\\{\mathbf{P}_{\left[n\right]}:n\in\mathcal{N}\\}$ that stratifies (50) and
(51). For the case of $n=N$, from (50) it follows that
$\mathbf{P}_{\left[N\right]}=\mathbf{T}$. Now, we assume that a positive
definite Hermitian symmetric matrix $\mathbf{P}_{\left[n\right]}$ exist, so
that
$\mathscr{J}_{\left[n\right]}=\mathbf{x}^{a\mathsf{H}}_{\left[n\right]}\mathbf{P}_{\left[n\right]}\mathbf{x}^{a}_{\left[n\right]}$.
Then, for $\mathscr{J}_{\left[n-1\right]}$, we have

$\displaystyle\mathscr{J}_{\left[n-1\right]}=$
$\displaystyle\sum^{N}_{m=n-1}J_{\left[m\right]}=J_{\left[n-1\right]}+\underbrace{\sum^{N}_{m=n}J_{\left[m\right]}}_{\mathscr{J}_{\left[n\right]}}=\mathbf{x}^{a\mathsf{H}}_{\left[n-1\right]}\boldsymbol{\Psi}_{\left[n-1\right]}\mathbf{x}^{a}_{\left[n-1\right]}+\mathbf{x}^{a\mathsf{H}}_{\left[n\right]}\mathbf{P}^{a}_{\left[n\right]}\mathbf{x}^{a}_{\left[n\right]}.$
(69)

Substituting $\mathbf{x}^{a}_{\left[n\right]}$ from (47), into (69), we arrive
at

$\mathscr{J}_{\left[n-1\right]}=\mathbf{x}^{a\mathsf{H}}_{\left[n-1\right]}\mathbf{P}_{\left[n-1\right]}\mathbf{x}^{a}_{\left[n-1\right]}$
(70)

with

$\mathbf{P}_{\left[n-1\right]}=\boldsymbol{\Psi}_{\left[n-1\right]}+\left(\mathbf{A}^{a\mathsf{H}}+\mathbf{G}^{\mathsf{H}}_{\left[n-1\right]}\mathbf{B}^{\mathsf{H}}\right)\mathbf{P}^{a}_{\left[n\right]}\left(\mathbf{A}^{a}+\mathbf{B}\mathbf{G}_{\left[n-1\right]}\right).$

Therefore, via induction, the expressions (69) and (70) show that there exists
a positive definite matrix $\mathbf{P}_{\left[n-1\right]}$ so that
$\mathscr{J}_{\left[n-1\right]}=\mathbf{x}^{a\mathsf{H}}_{\left[n-1\right]}\mathbf{P}_{\left[n-1\right]}\mathbf{x}^{a}_{\left[n-1\right]}$.

### SM-9: Application Examples

###  IV-A Three-Phase Power Distribution Systems

The power grid is based on three-phase power production and distribution. To
this end, instantaneous voltages can be combined into a quaternion signal

$\displaystyle q_{\left[n\right]}=$ $\displaystyle\imath
V_{a,\left[n\right]}\sin\left(2\pi f\Delta
Tn+\phi_{a,\left[n\right]}\right)+\jmath V_{b,\left[n\right]}\sin\left(2\pi
f\Delta Tn+\phi_{b,\left[n\right]}+\frac{2\pi}{3}\right)$ (71)
$\displaystyle+\kappa V_{c,\left[n\right]}\sin\left(2\pi f\Delta
Tn+\phi_{c,\left[n\right]}+\frac{4\pi}{3}\right)$

where $\\{V_{a,\left[n\right]},V_{b,\left[n\right]},V_{c,\left[n\right]}\\}$
(cf.
$\\{\phi_{a,\left[n\right]},\phi_{b,\left[n\right]},\phi_{c,\left[n\right]}\\}$)
are the instantaneous amplitudes (cf. instantaneous phases) on phase $a$, $b$,
and $c$, while $f$ is the system frequency555Frequency is common among all
phases due to structure of power generation systems. and $\Delta T$ is the
sampling interval.

For a three-phase power system operating under optimal (balanced) conditions,
that is, $V_{a,\left[n\right]}=V_{b,\left[n\right]}=V_{c,\left[n\right]}$ and
$\phi_{a,\left[n\right]}=\phi_{b,\left[n\right]}=\phi_{c,\left[n\right]}$, the
system can be treated in the same manner as single phase operations. However,
for a more general presentation and considering that all elements of (71) have
the same frequency, it was shown in [15] that $q_{\left[n\right]}$ traces an
ellipse in a plane of the three-dimensional imaginary subspace of
$\mathbb{H}$. This is shown in Fig. 5, where the system voltages of a balanced
and an unbalanced three-phase system are presented. Now, consider a new set of
imaginary units, $\\{\zeta,\zeta^{\prime},\zeta^{\prime\prime}\\}$ such that

$\zeta\zeta^{\prime}=\zeta^{\prime\prime},\hskip
5.0pt\zeta^{\prime}\zeta^{\prime\prime}=\zeta,\hskip
5.0pt\zeta^{\prime\prime}\zeta=\zeta^{\prime}$ (72)

and $\zeta$ -$\zeta^{\prime}$ reside in the same plane as $q_{\left[n\right]}$
(resulting in $\zeta^{\prime\prime}$ being normal to the plane). An arbitrary
ellipse in the $\zeta-\zeta^{\prime}$ plane can then be expressed as

$q_{\left[n\right]}=\big{(}A_{\left[n\right]}\sin(2\pi f\Delta
Tn+\phi_{\zeta,\left[n\right]})+\zeta^{\prime\prime}B_{\left[n\right]}\sin(2\pi
f\Delta Tn+\phi_{\zeta^{\prime},\left[n\right]})\big{)}\zeta.$ (73)

where $A_{\left[n\right]},B_{\left[n\right]}\in\mathbb{R}$, are instantaneous
amplitudes and
$\phi_{\zeta,\left[n\right]},\phi_{\zeta^{\prime},\left[n\right]}$ are
instantaneous phases. Replacing the $\text{sin}(\cdot)$ and
$\text{cos}(\cdot)$ functions in (73) with their polar representations from
(5) and after some mathematical manipulation, we have

$\displaystyle q_{\left[n\right]}=$
$\displaystyle\underbrace{\Big{(}\frac{A_{\left[n\right]}e^{\zeta^{\prime\prime}(\phi_{\zeta,\left[n\right]})}}{2\zeta^{\prime\prime}}+\frac{B_{\left[n\right]}e^{\zeta^{\prime\prime}(\phi_{\zeta^{\prime},\left[n\right]})}}{2}\Big{)}e^{\zeta^{\prime\prime}(2\pi
f\Delta Tn)}\zeta}_{q^{+}_{\left[n\right]}}$ (74)
$\displaystyle-\underbrace{\Big{(}\frac{A_{\left[n\right]}e^{-\zeta^{\prime\prime}(\phi_{\zeta,\left[n\right]})}}{2\zeta^{\prime\prime}}+\frac{B_{\left[n\right]}e^{-\zeta^{\prime\prime}(\phi_{\zeta^{\prime},\left[n\right]})}}{2}\Big{)}e^{-\zeta^{\prime\prime}(2\pi
f\Delta Tn)}\zeta}_{q^{-}_{\left[n\right]}}$

where

$q^{+}_{\left[n\right]}=e^{\zeta^{\prime\prime}(2\pi f\Delta
T)}q^{+}_{\left[n-1\right]}\hskip 9.10509pt\text{and}\hskip
9.10509ptq^{-}_{\left[n\right]}=e^{-\zeta^{\prime\prime}(2\pi f\Delta
T)}q^{-}_{\left[n-1\right]}.$ (75)

denote the two counter-rotating signals corresponding to the active and
reactive elements of the power signal. These two counter rotating ellipses are
shown in Fig. 6. Note that ideally, $q^{-}_{\left[n\right]}=0$ with a signal
only travelling in one direction (generation to consumption). From (75), state
of the system is represented via the vector

$\begin{bmatrix}\varphi_{\left[n+1\right]}\\\ q^{+}_{\left[n+1\right]}\\\
q^{-}_{\left[n+1\right]}\end{bmatrix}=\begin{bmatrix}\varphi_{\left[n\right]}\\\
\varphi_{\left[n\right]}q^{+}_{\left[n\right]}\\\
\varphi^{*}_{\left[n\right]}q^{-}_{\left[n\right]}\end{bmatrix}$ (76)

where $\varphi_{\left[n\right]}=e^{\zeta^{\prime\prime}2\pi
f_{\left[n\right]}\Delta T}$ with $f_{\left[n\right]}$ denoting the
instantaneous frequency of the system. The quaternion model has two key
advantages over its complex-valued duals. First, it allows to distinguish
between different classes of faults [15]. Second, it produces estimates of
system frequency with very low bias. This is demonstrated in Fig. 7, where
state of the system is tracked using the derived model in (76) with frequency
of the system is found from

$\hat{f}_{\left[n\right]}=\frac{1}{2\pi\Delta
T}\text{atan}\left(\frac{\|\Im(\hat{\varphi}_{n})\|}{\Re(\hat{\varphi}_{n})}\right)\cdot$
(77)

Figure 5: Geometric view of the system voltages, $q_{\left[n\right]}$, and the
corresponding phasor diagrams of a balanced and an unbalanced three-phase
system. Note that p.u. stands for per unit. Figure 6: System voltage,
$q_{\left[n\right]}$, positive sequence element, $q^{+}_{\left[n\right]}$, and
negative sequence element, $q^{-}_{\left[n\right]}$, of an unbalanced three-
phase system suffering from an 80% drop in the amplitude of one phase and $20$
degree shifts in the other phases. Note that p.u. stands for per unit.

Figure 7: Frequency estimation performance using the derived quaternion-vlaued
model, quaternion frequency estimator (QFE), as compared to complex-valued
linear (CLFE) and widely liner (CWLFE) techniques detailed in [15]; a) mean
square error (MSE) performance for balanced system, b) MSE performance for the
unbalanced system in Figure 6, c) bias performance for balanced system, and d)
bias performance for the unbalanced system in Figure 6.

######  Remark 10.

The models in (74) and (75) have been used for analysis and fault detection in
three-phase power systems [26, 7, 15, 6, 27].

###  IV-B Fly-by-Wire Flight Controller

In their most general formulation, the equation of motion governing the
rotation of solid bodies can be summarised as

$\frac{\partial\varphi}{\partial t}=\phi\hskip 11.38092pt\text{ and }\hskip
11.38092pt\frac{\partial\phi}{\partial t}=\nu$ (78)

where $\nu$ represents the input to the system, typically torque or a related
variable concerning angular acceleration), while $\varphi$ represents the
angular speed with $\phi$ indicating its rate of change. Transforming the
equations in (78) into its discrete-time formulation yields

$\begin{bmatrix}\varphi_{\left[n+1\right]}\\\
\phi_{\left[n+1\right]}\end{bmatrix}=\mathrm{f}\left(\begin{bmatrix}\varphi_{\left[n+1\right]}\\\
\phi_{\left[n+1\right]}\end{bmatrix},u_{\left[n\right]}\right)=\begin{bmatrix}1&\Delta
T\\\ 0&1\end{bmatrix}\begin{bmatrix}\varphi_{\left[n\right]}\\\
\phi_{\left[n\right]}\end{bmatrix}+\begin{bmatrix}\Delta T^{2}/2\\\ \Delta
T\end{bmatrix}u_{\left[n\right]}$ (79)

where $\Delta T=\text{$0.04$\,s}$ denotes the sampling interval.

In order to present an illustrative example, consider the input, $u_{n}$, to
be linear feedback regulator found via the framework in Section III-D, with

$\mathbf{B}=\begin{bmatrix}\Delta T^{2}/2\\\ \Delta
T\end{bmatrix}\text{and}\hskip 3.41418pt\hskip
3.41418pt\mathbf{A}=\begin{bmatrix}1&\Delta T\\\ 0&1\end{bmatrix}$

while $\mathbf{Q}=\mathbf{I}$, $\mathbf{T}=50\times\mathbf{I}$, and
$\mathbf{R}=10\times\mathbf{I}-1.875\times\left(\mathbf{1}\mathbf{1}^{\mathsf{T}}\right))$.
Inputs were calculated for $1.6$ s long segments with the first half of the
calculated sequences implemented. Then, input sequences were re-calculated
using the new state-vector information. Fig. 8 shows the smooth trajectory of
$\varphi$.

######  Remark 11.

Note that quaternions have been an integral part of flight control systems
[28]. Thus, this is not a novel application. However, quaternions have
previously been decomposed into the real-part and a vector presenting the
imaginary part. The advantage of the framework derived in this article is two-
fold. First, the so-called augmented quaternion framework allows useful tools
from algebra, calculus, and statistics to be recovered in the quaternion
domain. The section on calculus is addressed in this article, while [21] deals
with quaternion statistics. Second, the approach presented in this article
removes the need for transformations, under which the elegance, physical
interpretation, and division algebra of quaternions is lost or obscured.

Figure 8: Trajectory of $\varphi$ under the proposed controller scheme.

###  IV-C Body Motion Tracking

The three Euler angles, as shown in Fig. 9, can be recorded using modern
motion sensors. In order to make this data suitable for processing each angle
was converted into a quaternion using the following transformations

$\begin{array}[]{c}\alpha\rightarrow\big{(}\text{cos}(\alpha),\text{sin}(\alpha)\big{)}\rightarrow
e^{i\frac{\alpha}{2}}\\\
\beta\rightarrow\big{(}\text{cos}(\beta),\text{sin}(\beta)\big{)}\rightarrow
e^{j\frac{\beta}{2}}\\\
\gamma\rightarrow\big{(}\text{cos}(\gamma),\text{sin}(\gamma)\big{)}\rightarrow
e^{k\frac{\gamma}{2}}\end{array}$ (80)

with each term respectively representing the roll, pitch, and yaw rotations in
the three-dimensional space. The three quaternion terms in (80) are now
combined into one quaternion-valued variable given by

$q=\big{(}e^{k\frac{\gamma}{2}}\big{)}\big{(}e^{j\frac{\beta}{2}}\big{)}\big{(}e^{i\frac{\alpha}{2}}\big{)}$
(81)

representing the overall rotation.

######  Remark 12.

Note that since the mapping in (80) is invertible, the resulting quaternion-
vlaued variable in (81) preserves dynamics of the signal.

Figure 9: Setting of the inertial body motion sensor with the fixed coordinate
system (blue), sensor coordinate system (red), and Euler angles (green). The
“N” axis is used as a visual guide to indicate the yaw angle. Figure taken
from [15].

As an illustrative example, we recover the experiment in [15] using the
frameworks in this article. Fig. 10 shows the recorded yaw angle measurements
containing discontinuities and its corresponding continuous transformation in
(80) which is ready for processing. The linear adaptive filter derived in SM-5
was considered in a one-step-ahead prediction setting to track body rotations.
Fig. 11 shows the estimates of the absolute value of the phase for the body
motion signal obtained by a quaternion-valued filter and its real-valued
quadrivariate counterpart, that is, the case where quaternion-valued body
motion signal was considered as a four-element vector with the real-valued
least mean square (LMS) algorithm used to track the signal. Notice that the
quaternion-valued filter had a significantly better performance than the
quadrivariate LMS algorithm.

Figure 10: Recorded yaw angle measurements with discontinuities at $\pm\pi$
(top) and its continuous transform (bottom). Figure 11: Absolute value of the
phase of quaternion-valued body motion signal employing the filter from SM-5
and a quadrivariate real-vlaued LMS approach.

###  IV-D Bearings-only-Tracking

The problem of tracking the position of a manoeuvring target using sensors
that can only measure the bearings of the target is considered. Commonly
referred to as bearings-only-tracking, this problem is often encountered in
passive radar, sonar, or infrared tracking applications. Since sensors can
produce estimates of the target bearing and not of its range, none of the
sensors can track the target location. Traditional solutions are based on
congregating the information of all sensors at a central processing unit and
resolving the location of the target through triangularisation. An elegant
solution to this problem can be found via local fusion of sensor estimates
over and ad-hoc network.

To this end, consider the target location at time instant $n$ to be presented
via the pure imaginary quaternion variable, $q_{\left[n\right]}$, in the
manner given in Section I-C. Then, dynamics of the target can be formulated as

$\displaystyle q_{\left[n+1\right]}=$ $\displaystyle
q_{\left[n\right]}+v_{\left[n\right]}\Delta
T+u_{\left[n\right]}\frac{1}{2}\left(\Delta T\right)^{2}$ (82) $\displaystyle
v_{\left[n+1\right]}=$ $\displaystyle
v_{\left[n\right]}+u_{\left[n\right]}\Delta T$

where $v_{\left[n\right]}$ and $u_{\left[n\right]}$ represent the speed and
acceleration (normally modelled as zero-mean Gaussian noise) of the target,
with $\Delta T$ representing the sampling interval. Moreover, sensor $l$ will
have observations

$y_{l,\left[n\right]}=\frac{q_{\left[n\right]}-L_{\mathsf{S}_{l}}}{\|q_{\left[n\right]}-L_{\mathsf{S}_{l}}\|}+w_{l,\left[n\right]}$
(83)

where $w_{l,\left[n\right]}$ and $L_{\mathsf{S}_{l}}$ represent observation
noise and location of sensor $l$.

As an illustrative example, the network shown in Fig. 12 with its nodes
distributed uniformly inside a $24\times 24\times 24$ cube was used to track a
target moving inside the cube through bearings-only measurements. The sampling
interval was set to $\Delta T=0.04$ s. The acceleration was considered as a
zero-mean Gaussian noise sequence with independent components and variance of
$10$, while observation noise was considered to be zero-mean Gaussian with
covariance of $10^{-4}$ and independent components. Agents used the framework
derived in Section III-A to formulate a local estimates of the target which
were then fused with those of the neighbouring agents via the framework of
Section III-B. As a worst case scenario the tracking solution of the agent
with only one connection is shown in Fig. 13. Note the accurate tracking
solution achieved via the division algebra of quaternions. Although a similar
approach can be implemented in vector algebras, the lack of division operation
means normalisation and sensor fusion operations require more complex
computations.

Figure 12: Network of $20$ nodes and 43 connections used in simulations.
Figure 13: Tracking solution of one agent in the bearings-only-tracking
example.

###  IV-E Quantum Computing and the Compiler Problem

Recall the quaternion-valued model of a qubit given in Example 1. Although the
model of a qubit with quaternions is not unique to that presented in this
article. However, considering that the input and output of quantum computers
are unit pure imaginary quaternions, the setting derived in this article
allows all quantum computing operations to be modelled effectively using the
division algebra of quaternions, e.g., involutions and rotations, operations.
To this end, consider the a general quantum computing operation represented by
the mapping

$\mathbf{q}_{\text{out}}\leftarrow\mathrm{f}\left(\mathbf{q}_{\text{in}}\right)$

where $\mathrm{f}\left(\cdot\right)$ is a quaternion-valued function. On the
most general level, we would like to emulate this operation using a quantum
computer the size and possible operations of which are represented as
$\hat{\mathbf{q}}_{\text{out}}\leftarrow\hat{\mathrm{f}}\left(\mathbf{q}_{\text{in}},\mathbf{W}\right)$
with $\mathbf{W}$ denoting a parameter matrix with fixed values, representing
fixed operations. This problem can be formulated as

$\displaystyle\min_{\mathbf{W}}$
$\displaystyle\left\\{\|\mathrm{f}\left(\mathbf{q}_{\text{in}}\right)-\hat{\mathrm{f}}\left(\mathbf{q}_{\text{in}},\mathbf{W}\right)\|^{2}+\mathrm{g}\left(\mathbf{W}\right)\right\\}$
(84) so that: $\displaystyle\hskip 6.82881pt\mathbf{W}\in\mathcal{W}$ (85)

where $\mathrm{g}\left(\cdot\right):\mathbb{H}^{M\times
M}\rightarrow\mathbb{R}^{+}$ is a function representing the cost of
implementing a certain structure and $\mathcal{W}$ is the set of admissible
values for $\mathbf{W}$. The solution to this problem has been considered in
the form of factorisation and search algorithms [14]. These techniques,
however, can be accelerated using the $\mathbb{HR}$-calculus, to indicate the
steepest direction of reduction in cost. On a simpler level, the minima of
(84) is calculable via gradient descent. This solution can then be projected
to the closest $\mathbf{W}\in\mathcal{W}$ in order to meet the condition in
(85).

######  Remark 13.

The most important application of the setting derived in this article for
modelling qubit is that it allows for most concepts in linear/nonlinear
dynamic programming to be used in the context of quantum computing. Note that
these concepts are expandable to the quaternion domain using the
$\mathbb{HR}$-calculus.