Quaternion Recurrent Neural Network
with Real-Time Recurrent Learning
and Maximum Correntropy Criterion

Denote a quaternion, q, as

where $q_{a},q_{b},q_{c},q_{d}\in\mathbb{R}$ and the imaginary units i, j, and k satisfy $i^{2}=j^{2}=k^{2}=ijk=-1$, $ij=-ji=k$, $jk=-kj=i$, $ki=-ik=j$. The set of quaternions is defined as $\mathbb{H}\triangleq\{q=q_{a}+iq_{b}+jq_{c}+kq_{d}\;|\;q_{a},q_{b},q_{c},q_{d}\in\mathbb{R}\}$. The multiplication of two quaternions in $\mathbb{H}$ is noncommutative due to the properties of the imaginary units. The real part of $q$ is defined as $\operatorname{Re}\{q\}=q_{a}$, whereas the imaginary part is $\operatorname{Im}\{q\}=iq_{b}+jq_{c}+kq_{d}$. The conjugate of $q$ is $q^{*}=\operatorname{Re}\{q\}-\operatorname{Im}\{q\}=q_{a}-iq_{b}-jq_{c}-kq_{d}$. The modulus of a quaternion is denoted as $|q|=\sqrt{qq^{*}}$. The inverse of a quaternion $q\neq 0$ is $q^{-1}=q^{*}/|q|^{2}$. For any quaternion, q, and a nonzero quaternion, $\mu$, the transformation [17]

describes a rotation of q.

A quaternion function $f:\mathbb{H}\rightarrow\mathbb{H}$, defined as $f(q)=f_{a}(q_{a},q_{b},$ $q_{c},q_{d})+if_{b}(q_{a},q_{b},q_{c},q_{d})+jf_{c}(q_{a},q_{b},q_{c},q_{d})+kf_{d}(q_{a},q_{b},q_{c},$ $q_{d})$ is called real-differentiable, if $f_{a},\,f_{b},\,f_{c},\,f_{d}$ are differentiable as functions of the real variables $q_{a},\,q_{b},\,q_{c},\,q_{d}$ [18]. If $f:\mathbb{H}\rightarrow\mathbb{H}$ is real-differentiable, then the left GHR derivative of the function $f$ with respect to $q^{\,\mu}$ ($\mu\neq 0,\mu\in\mathbb{H}$) is

where $q=q_{a}+iq_{b}+jq_{c}+kq_{d}$, $q_{a},q_{b},q_{c},q_{d}\in\mathbb{R}$, and $\frac{\partial f}{\partial q_{a}},\frac{\partial f}{\partial q_{b}},$ $\frac{\partial f}{\partial q_{c}},\frac{\partial f}{\partial q_{d}}\in\mathbb{H}$ are the partial derivatives of $f$ with respect to $q_{a},q_{b},q_{c},q_{d}$. If $f:\mathbb{H}\rightarrow\mathbb{H}\,$, $g:\mathbb{H}\rightarrow\mathbb{H}\,$, the product rule is [10]

the chain rule is

and the rotation rule is

Denote the quaternion gradient of a function $f:\mathbb{H}^{n}\rightarrow\mathbb{R}$ as

where $\Bigl{(}\frac{\partial f}{\partial\textbf{{q}}}\Bigr{)}^{T}$ is the transpose of $\frac{\partial f}{\partial\textbf{{q}}}$ [9]. The quaternion Jacobian matrix of $\textbf{f}:\mathbb{H}^{N\times 1}\rightarrow\mathbb{H}^{M\times 1}$ is then defined as

Note the convention that for two vectors $\textbf{f}\in\mathbb{H}^{M\times 1}$ and $\textbf{q}\in\mathbb{H}^{N\times 1}$, $\frac{\partial\textbf{f}}{\partial\textbf{q}}$ is a matrix for which the $(m,n)$th element is $(\partial f_{m}/\partial q_{n})$, thus, the dimension of $\frac{\partial\textbf{f}}{\partial\textbf{q}}$ is $M\times N$ [9].

The aim is to form a quaternion function $z$, based on a sequence $(\textbf{s}_{1},d_{1}),(\textbf{s}_{2},d_{2}),\dots,(\textbf{s}_{N},d_{N})$, where $\textbf{s}_{p}$ is the system input, and $d_{p}$ is the expected response, both quaternion valued, for an instant time $p$. We denote the expected response as $d_{p}=z_{opt}(\textbf{s}_{p})+\xi_{p}$, where $z_{opt}$ is a nonlinear quaternion function to estimate, and $\xi_{p}$ is noise with an arbitrary probability density function that does not need to be Gaussian.

When finding $z$, we ideally would like to solve the empirical risk minimization problem, that is, to minimize

where $Q$ is a quaternion vector space. Notably, the use of MSE can result in large variations for the weights, or shift in the output, when the noise distribution contains outlying points, is non-symmetric, or has a nonzero mean.

The correntropy is defined as $V_{\sigma}(X,Y)=E\left[\kappa_{\sigma}(X-Y)\right]$, where $X$ and $Y$ are quaternion random variables, and $\kappa_{\sigma}(\cdot)$ is a symmetric positive-definite quaternion kernel, with kernel size $\sigma$ [19]. For real-valued inputs, the Gaussian kernel can be expressed as

where $X$ and $Y$ are real valued [20]. An analogous Gaussian-based kernel for quaternion data may be expressed as

where $X$ and $Y$ are the quaternions $X=X_{r}+iX_{i}+jX_{j}+kX_{k}$ and $Y=Y_{r}+iY_{i}+jY_{j}+kY_{k}$ [20] (see Fig. 1). The factor of 4 is introduced to normalize the effect of the four components in the quaternion as opposed to a single real component. In the complex case, the factor of 2 would be included for normalization.

Let $\textbf{h}_{a}^{(l)}(n)$ denote the hidden state of the $a$th neuron in the $l$th layer at time $n$, defined as

where

Here, $\textbf{u}_{ab}^{(l)}$ represents the vector of recurrent quaternion weights, $\Phi$ is the quaternion split activation function (see Appendix VI-B), $\textbf{v}_{b}^{(l-1)}(n)$ represents the output of the $b$-th neuron in layer $l-1$ at time $n$, $\textbf{h}_{b}^{(l-1)}(n-1)$ is the hidden state of the $b$-th neuron in layer $l-1$ at time $n-1$, and $\boldsymbol{\otimes}$ is the Hamilton product.

Denote the error vector by $\mathbf{e}(n)=\mathbf{d}(n)-\mathbf{h}(n)$. Given the quaternion Gaussian kernel in (11), the correntropy between $\textbf{{h}}(n)$ and the desired output $\mathbf{\textbf{d}}(n)$ becomes

where

Usually, the joint probability density function is unknown and only a finite set of samples is available. To estimate the correntropy, the expectation can be approximated based on

The correntropy should account for the even moments between $\textbf{d}(n)$ and $\textbf{h}(n)$, especially the magnitude $|\textbf{d}(n)-\textbf{h}(n)|$. As $\sigma$ increases, higher-order terms should diminish, making the second-order term dominant. This facilitates using a gradient based method to maximize correntropy.

As we aim to use gradient descent, instead of maximising $\hat{V}_{N,\sigma}(\textbf{d}(n),\textbf{h}(n))$, the objective of the network training is to minimise the loss function $-\hat{V}_{N,\sigma}(\textbf{d}(n),\textbf{h}(n))$, expressed as

where $n$ is the current time index and $N$ is the total number of time steps. Using the correntropy with a Gaussian kernel function as a cost function gives our problem kernel Hilbert space attributes, which differs from the linear MSE that only considers the Euclidean distance of errors. Essentially, the correntropy represents a weighted sum of error distances, while the kernel function modifies the influence of each error by mapping it from the input space to the reproducing kernel Hilbert space.

We now derive the quaternion real-time recurrent learning algorithm for the QRNN. The quaternion gradient of the correntropy cost function can be expressed as

To find $\frac{\partial J(n)}{\partial\textbf{q}}$, we use the GHR chain rule [10] to obtain

The first term $\frac{\partial\kappa_{\sigma}(\textbf{e}(n))}{\partial\textbf{e}^{\mu}(n)}$ can be evaluated as

Upon substituting (20) into the formula for $\frac{\partial J(n)}{\partial\textbf{q}}$, we obtain

Then, upon employing the GHR chain rule [10], the derivative of $|\textbf{e}(n)|^{2}$ can be calculated as

whereby the GHR product rule gives

Given that $\textbf{e}(n)=\textbf{d}(n)-\textbf{h}(n)$, we obtain

where $J_{\textbf{q}^{\mu}(n)}\triangleq\frac{\partial\textbf{h}(n)}{\partial\textbf{q}^{\mu}}$ is the Jacobian matrix of $\textbf{h}(n)$. Thus,

which allows us to arrive at

This can be simplified into

Finally, to find the quaternion gradient, we take the Hermitian (conjugate transpose) of (27), that is

Denote the quaternion error terms at time $n$ as $\delta^{(l)}(n)$, that is

with $\bar{\Phi}$ as the derivative of the activation $\Phi$, $\odot$ as the Hadamard product, and $(\textbf{W}^{(l+1)})^{H}$ and $(\textbf{U}^{(l)})^{H}$ as the Hermitians of the weight matrices for the $(l+1)$th and the $l$th layer, respectively. Then, the weight and bias updates for the output layer are

The weight and bias updates for the hidden layers are

Here, $\alpha>0$ is the learning rate, and $\textbf{U}^{(l)}$ represents the matrix of recurrent quaternion weights for layer $l$.

In the context of lung cancer radiotherapy, tracking the position of infrared-reflective markers on the chest is a method used to approximate the location of the tumour. Nonetheless, the precision of radiation delivery in radiotherapy systems is constrained by the inherent latency due to limitations in robotic control. Compensating for delays is crucial to reduce the damage to healthy tissues. Employing RNN with online learning can facilitate predictions that adapt to the dynamic nature of respiratory signals, potentially mitigating this issue [21]. Unlike offline methods, which do not change after initial training, online methods continually update the synaptic weights of the network with each new data input. This continuous updating allows the neural network to adjust to the patient’s evolving breathing patterns, offering improved resilience against complex movements. Online learning is particularly advantageous in medical settings, where acquiring extensive training datasets can be challenging. It provides a way to adapt to data not included in the initial training set. Adaptive or dynamic learning techniques have been frequently utilized in radiotherapy. Studies have shown the effectiveness of these approaches compared to static models [22, 23, 24]. One such dynamic method is RTRL [13], which has been applied in various systems like the Cyberknife Synchrony system [24] and the SyncTraX system [25]. It has also been used to track the position of internal points within the chest [21], demonstrating its broad applicability in medical technology.

While prior studies in the field of respiratory motion prediction primarily concentrated on univariate signal analysis, our simulation focuses on the prediction of 3D marker position, represented as pure quaternions, setting a prediction horizon of two seconds. This is achieved through the use of a QRNN trained with RTRL. We compare its effectiveness, when associated with the MCC loss, against other forecasting algorithms. Furthermore, we conducted predictions for all three markers concurrently, enabling the QRNN to further identify and leverage the interrelationships in their movements.

The data used for this study was collected from nine instances of 3D positional tracking of three external markers placed on the chest and abdomen of subjects resting on a HexaPOD treatment couch. The tracking was accomplished using an NDI Polaris infrared camera. Each recorded sequence varied in length, ranging from 73 to 320 seconds, and was sampled at a frequency of 10 Hz. The movement trajectories of the markers were recorded in three dimensions: superior-inferior (ranging from 6 mm to 40 mm), left-right (2 mm to 10 mm), and antero-posterior (18 mm to 45 mm). Out of these recordings, five depict normal breathing patterns, while the remaining four capture subjects engaging in activities such as talking or laughing. Further specifics about the public dataset are elaborated in [26, 27]. The data from subjects was divided into two groups to assess the robustness of the algorithms.

We compared the performance of the QRNN equipped with RTRL and the MCC loss function (Sec. III) against its counterpart that employs the MSE loss function (see Appendix VI-A). The performances of the standard RNN with RTRL and MSE loss as described in [21] or MCC loss as defined in (10), as well as the Quaternion Least Mean Square (QLMS) algorithm [28], and the real-valued Least Mean Square (LMS) algorithm, were also included in the comparison of the prediction methods. Note that the predictions for all three markers were conducted concurrently for the standard RNNs with RTRL. The QRNNs and RNNs were characterized by one hidden layer. The activation function was the hyperbolic tangent function for RNNs, and the split hyperbolic tangent function for QRNNs (see Appendix VI-B). The cross-validation metric was the RMSE. The ranges of hyper-parameters for cross-validation with grid search were as follows. For QRNNs and RNNs with RTRL, the learning rate range was $\eta\in\{0.02,0.05,0.1,0.2\}$. The range in the number of hidden units was $h\in\{10,20,30,40,45\}$ for QRNNs and $h\in\{20,40,60,80,90\}$ for RNNs. For QLMS and LMS, the parameter range was $\eta\in\{0.002,0.005,0.01,0.02,0.05,0.1,0.2\}$. For all prediction methods, the regressor range defined as the number of time steps was $l\in\{10,30,50,70,90\}$. There were 50 runs successively performed for cross-validation, and 300 runs for evaluation. Updating QRNNs and RNNs with the gradient rule may induce numerical instability. The gradient norm was therefore clipped to avoid large weight updates, while the optimization method was stochastic gradient descent. The training and validation of these models were conducted in the first minute of each recorded sequence, both for 30 seconds. The four error types, presented in Table I, used to evaluate the prediction performance are implemented following [27].

The QRNN equipped with RTRL and the MCC loss produced the lowest mean average error, RMSE, and nRMSE averaged over all the sequences as well as in the context of irregular breathing (Table I). The QRNN with RTRL and the MSE loss performed slightly better than its counterpart with the MCC loss for predicting sequences of regular breathing, in terms of the RMSE and nRMSE metrics. The robustness of the QRNN with RTRL and the MCC loss against irregular movements was further confirmed as its nRMSE exhibited an increase of 8.2% between regular and irregular breathing patterns, which is the smallest among the algorithms assessed. More generally, the QRNNs with RTRL displayed consistently better performance for forecasting the breathing sequences than the standard RNNs with RTRL, QLMS, and LMS algorithms. The compact 95% confidence intervals accompanying the performance metrics presented in Table I suggest that choosing 300 test runs was adequate to yield precise results. In addition, the jitter was evaluated. It refers to the fluctuation or oscillation amplitude in the predicted signal, and its extent can significantly affect robotic control during treatment. The QRNN with RTRL and MSE displayed a slightly lower jitter overall than its counterpart with the MCC loss. The QRNN with RTRL and MCC exhibited the lowest jitter in the context of irregular breathing, significantly outperforming the other algorithms, including the RNN with RTRL and MCC (0.77 mm against 0.88 mm) and the QLMS and LMS algorithms (0.77 mm against 1.73 mm).

The simulation has certain limitations, particularly regarding the quantity and length of the sequences utilized, which are relatively modest. Nonetheless, the dataset used encompasses a broad spectrum of respiratory patterns, including various shifts, drifts, slow and abrupt irregular movements, and both active and calm breathing states. A notable aspect of the QRNNs with RTRL and the MCC or MSE loss function we examined is their ability to perform online learning, which does not require extensive pre-existing data to make precise predictions. This is evident from the high accuracy achieved with just one minute of training data. Based on these factors, we believe our results could be extrapolated to larger datasets. Another strength of the simulations is the open-source dataset we used, which promotes reproducibility of our findings. This contrasts with many prior studies in this field, which often rely on private datasets, complicating comparative analyses. We also addressed challenging scenarios like laughing and talking, which are generally controlled in clinical settings. Analyzing performance under such conditions provides insights into other less predictable events that may occur during treatment, such as yawning, hiccupping, or coughing. The current clinical practice is to halt radiation when such anomalies are detected. By differentiating between regular and irregular breathing patterns, we could objectively assess and measure the resilience of the compared algorithms. Given that almost half of our dataset consisted of irregular breathing sequences, the average numerical error metrics across all nine sequences might be higher than what would be encountered in more standard scenarios. Moreover, the QRNN models do take roughly 50% longer to train than the RNN models (Table II). The QRNN with RTRL and MCC demonstrated an average computation time per time step of 184.5 ms. Note that this computation time is less than the approximate marker position sampling interval, which is around 400 ms.

Denote the split activation function as

with $\Phi_{\zeta}(\cdot)$ representing any real-valued activation function. To perform backpropagation using split activation functions, [29] proposed a ”pseudo-gradient” update wherein the gradient is computed component-wise. Accordingly, the ”pseudo-derivative” of $\Phi(\cdot)$ in (43), written as $\bar{\Phi}(\cdot)$, is given by

On the other hand, the compact GHR derivative of any split activation function $\Phi(\cdot)$ is defined as [9]