Generalization of Neural Networks on Second-Order Hypercomplex Numbers

Abstract

The vast majority of existing neural networks operate by rules set within the algebra of real numbers. However, as theoretical understanding of the fundamentals of neural networks and their practical applications grow stronger, new problems arise, which require going beyond such algebra. Various tasks come to light when the original data naturally have complex-valued formats. This situation is encouraging researchers to explore whether neural networks based on complex numbers can provide benefits over the ones limited to real numbers. Multiple recent works have been dedicated to developing the architecture and building blocks of complex-valued neural networks. In this paper, we generalize models by considering other types of hypercomplex numbers of the second order: dual and double numbers. We developed basic operators for these algebras, such as convolution, activation functions, and batch normalization, and rebuilt several real-valued networks to use them with these new algebras. We developed a general methodology for dual and double-valued gradient calculations based on Wirtinger derivatives for complex-valued functions. For classical computer vision (CIFAR-10, CIFAR-100, SVHN) and signal processing (G2Net, MusicNet) classification problems, our benchmarks show that the transition to the hypercomplex domain can be helpful in reaching higher values of metrics, compared to the original real-valued models.

1. Introduction

The exploration of neural networks in complex algebra is a modern area of research. Many tasks in classification, signal pattern recognition, and the generation of new signals inherently present data in a complex format. Recent studies have shown that complex-valued neural networks outclass real-valued ones in terms of representational capacity [1] and learning speed [2]. Other publications aim to develop complex-valued layers and activation functions for neural networks (see the survey in [3]). The progress allows researchers to build effective models with better convergence and accuracy [4,5].

These promising results have inspired us to further advance the generalization of neural networks, extending to the dual $\mathbb{D}$ and double or complex-split $\mathbb{S}$ algebras. The basic mathematical operations for these numbers are as follows:

$$\begin{array}{c}(x_1+\tau y_1)\pm (x_2+\tau y_2)=(x_1\pm x_2)+\tau (y_1\pm y_2)\\ (x_1+\tau y_1)·(x_2+\tau y_2)=(x_1{x}_2+\sigma y_1{y}_2)+\tau (x_1{y}_2+y_1{x}_2)\\ \frac{x_1+\tau y_1}{x_2+\tau y_2}=\frac{x_1{x}_2-\sigma y_1{y}_2}{x_2^2-\sigma y_2^2}+\tau \frac{x_2{y}_1-x_1{y}_2}{x_2^2-\sigma y_2^2}\\ {(x+\tau y)}^{*}=x-\tau y\end{array}$$

In addition to complex numbers, dual and double numbers have found applications in physics, e.g., screw theory [6] or relativistic cosmology [7]. Dual numbers also make it possible to automatically compute derivatives of functions [8,9,10]. Double-valued networks have been used for solving CV tasks [11]. Therefore, using hypercomplex numbers for deep learning is well justified. There have been studies on neural networks based on Clifford algebra [12], specifically on dual numbers [13], and room for research still remains. The goal of this paper is to generalize the approach of using second-order hypercomplex algebras for neural networks.

There is another reason for our interest in dual and double numbers. As real-valued layers expand to the hypercomplex algebras, the number of operations increases. The magnitude of this growth depends on the computational complexity of basic operations on the numbers. For example, one multiplication of two complex numbers is more expensive than four multiplications of real numbers. This influences the computational cost of the operators discussed in Section 2. Thus, the convolution of a complex-valued input $z=x+iy$ with a complex-valued filter matrix $W=A+iB$ can be rewritten as $Wz=(Ax-By)+i(Ay+Bx)$, which implies four real-valued matrix multiplications (see Figure 1). The dual-valued convolution of a dual filter $W=A+\epsilon B$ and a dual input $z=x+\epsilon y$ evaluates $Wz=Ax+\epsilon (Ay+Bx)$ because ${\epsilon }^2=0$. This type of convolution only has three real-valued matrix multiplications (see Figure 1), which gains a $+25\%$ theoretical performance speed-up, compared to the complex-valued convolution, and interconnects x and y in the process. As for double numbers, the cost of the naive implementation of the convolution is the same as for complex numbers because $(A+jB)·(x+jy)=(Ax+By)+j(Bx+Ay)$. However, the features of double numbers enable variable substitution that allows changing to a so-called diagonal basis. Under that new setup, a double-valued convolution is effectively only two real-valued convolutions instead of four. Section 3 of this paper is dedicated to the optimizations of hypercomplex networks, including the diagonal representation of double numbers (Section 3.3). In this paper, we will mainly use the original representation of double numbers, for the sake of following the universal convention.

It is worth mentioning that double numbers in the diagonal representation have a significant flaw: the inability to connect x and y. However, the influence of this problem can be mitigated by a special batch normalization operator and/or an activation function for double numbers, which would make them connect. This adjustment is very important. As the experiments with real-valued networks have shown [14], an appropriate batch normalization helps accelerate the training process and reach better accuracy. We will further discuss the components of hypercomplex-valued neural networks in detail.

In Section 4 of this paper, we show that the adaptation of neural networks to the algebra of hypercomplex numbers allows us to improve precision. The important conclusion here is that the proposed types of neural networks, being barely explored so far, hold great potential for further research.

2. Methodology

This section is dedicated to the definitions of key elements of a hypercomplex neural network.

2.1. MindSpore Hypercomplex Data Representation

Until recently, only complex operators have been supported in popular frameworks, such as TensorFlow and PyTorch. The list of supported operations only contains basic manipulations with complex-valued tensors, as well as convolutional and linear layers. As of now, the situation has changed dramatically. MindSpore has become the first platform to support all the core layers needed to build second-order hypercomplex neural networks.

MindSpore allows researchers to create complex, dual, and double tensors, and perform all basic operations, such as slicing, concatenation, conjugation, and accessing both components of tensors. Convolution and linear layers are supported by MindSpore for all types of hypercomplex algebras of the second order. Other essential operators for deep learning, such as activation functions and algebra-specific batch normalization, are also integrated into MindSpore.

Within this framework, the hypercomplex tensor is represented as a two-channel real-valued tensor. Tensors of this shape are used as arguments and all operators return values in the same format. The peculiarities of every type of algebra are considered by a specific realization of network operators inside MindSpore, which interlinks the two channels of the tensors.

2.2. Hypercomplex Operators

The implementation of the operators mentioned above is based on our proposed method of generalizing the process of constructing building blocks for hypercomplex neural networks.

2.2.1. Convolution

Convolution is one of the most important operations in neural networks. In order to explain convolution in the hypercomplex algebra, we exploit the matrix representation of second-order hypercomplex numbers that use real numbers. The algebras of complex, dual, and double numbers $u=x+\tau y$ are all known to be isomorphic to the corresponding algebras of the second-order real matrices of the following form: $\left(\begin{array}{cc}x& y\\ \sigma y& x\end{array}\right).$ Recall that ${\tau }^2=\sigma$. Therefore, the convolution of a hypercomplex filter $W=W_x+\tau W_y$ and a hypercomplex input $u=x+\tau y$ can be expressed as follows:

$$W\ast u\sim \left(\begin{array}{cc}{W}_x& W_y\\ \sigma W_y& W_x\end{array}\right)\ast \left(\begin{array}{cc}x& y\\ \sigma y& x\end{array}\right)=\left(\begin{array}{cc}{W}_x\ast x+\sigma W_y\ast y& W_x\ast y+W_y\ast x\\ \sigma (W_x\ast y+W_y\ast x)& W_x\ast x+\sigma W_y\ast y\end{array}\right).$$

As mentioned earlier, we need to make three real-valued multiplications in the dual algebra ($W_x\ast x$, $W_y\ast x$, and $W_x\ast y$), but four in the cases of complex and double numbers ($W_x\ast x$, $W_y\ast x$, $W_x\ast y$, and $W_y\ast x$). For double numbers, however, changing variables, as shown in Section 3.3, enables moving to the diagonal representation, where the algebra becomes isomorphic to the algebra of real matrices of the following form: $\frac{1}{2}\left(\begin{array}{cc}x+y& x-y\\ x-y& x+y\end{array}\right)$, and the convolution only requires two convolutions in real numbers.

2.2.2. Linear Layer

Linear (fully connected, dense) layers are commonly used in neural networks as the head blocks in the CV (computer vision) domain and as the main operators in NLP (natural language processing) models. To generalize linear layers for the two-dimensional algebras, we use the matrix representation of second-order hypercomplex numbers. In the end, a linear layer with hypercomplex inputs and weights is equivalent to the superposition of real-valued linear layers:

$$\mathbb{H}L(W,B,u)=\mathbb{R}L(W_x,B_x,u_x)+{\tau }^2\mathbb{R}L(W_y,0,u_y)+\tau (\mathbb{R}L(W_x,0,u_y)+\mathbb{R}L(W_y,B_y,u_x)),$$

where $\mathbb{R}L(w,b,x)$ stands for a real-valued linear layer, with w, b, and x denoting weights, bias, and the input, respectively.

2.2.3. Average Pooling

The average pooling operation implies calculating the arithmetic mean for each patch of the feature map. This means collapsing every $n\times n$ square of the feature map into its average value. This is equivalent to convolution, with the stride being equal to the kernel size, where kernel weights are real numbers that are equal to $\frac{1}{n^2}$. Because the kernel $W=W_x+\tau W_y$ in this particular case is purely real, $W_y=0$, and the convolution formula can be simplified as follows:

$$W\ast u=W\ast (x+\tau y)=W_x\ast (x+\tau y)=W_x\ast x+\tau (W_x\ast y).$$

From this expression, average pooling is equivalent to two real-valued average pooling operations, with each independently applied to every component of the input data $u=x+\tau y$:

$$\mathbb{H}AvgPool\left(u\right)=\mathbb{R}AvgPool\left(x\right)+\tau \mathbb{R}AvgPool\left(y\right).$$

2.2.4. ReLU

Activation functions are used to introduce non-linearity into neural networks. There are multiple activation functions based on real numbers, and greater varieties are based on hypercomplex numbers. Among the real-valued activations, there is a family of $ReLU$-type functions, which do not suffer from the vanishing gradient problem. This property makes them quite stable. Functions of this type are used in the complex algebra as well. For example, the authors of [15] applied $ReLU$ to the real and imaginary parts separately. In this paper, we extend this definition and apply it to other algebras:

$$\mathbb{H}ReLU\left(u\right)=\mathbb{R}ReLU\left(x\right)+\tau \mathbb{R}ReLU\left(y\right).$$

This activation approach is not perfect. First, in contrast to $\mathbb{R}ReLU$ on real numbers, $\mathbb{H}ReLU\left(u\right)$ may not be equal to 0 or u. Second, the real and imaginary parts are processed independently and, therefore, are not fully interlinked. Despite that, the method shows good experimental results, and its obvious advantage is simplicity.

2.3. Norm of Hypercomplex Numbers

In deep learning, the ability to properly update model parameters is important for training a model with better accuracy on a test dataset. The weights may not be efficiently updated if elements of intermediate tensors are not normalized due to the gradient explosion or vanishing problem. In real-valued models, the solution is the batch normalization procedure. The idea is to balance tensor elements based on the variance of the batch. To develop a counterpart operator for neural networks based on second-order algebras, we must find an equivalent of the absolute value function. Here, we define a concept of the norm (or modulus) of hypercomplex numbers. In mathematical or physical literature, the modulus of a complex number is traditionally defined as

$$\mathbb{C}\left|z\right|=\sqrt{z{z}^{\ast }}=\sqrt{(x+iy)(x-iy)}=\sqrt{x^2+y^2}.$$

(1)

However, the generalization of this method to other types of hypercomplex numbers does not work quite well:

$$\mathbb{H}\left|u\right|=\sqrt{u{u}^{\ast }}=\sqrt{(x+\tau y)(x-\tau y)}=\sqrt{x^2-\sigma y^2}=\left\{\begin{array}{c}\left|x\right|,u\in \mathbb{D}\hfill \\ \sqrt{x^2-y^2},u\in \mathbb{S}\hfill \end{array}\right.$$

This result is hardly applicable to our purpose. The dual norm does not depend on the dual part. In the case of double numbers, the function is not well-defined for half of the elements. In [16], the authors derive the norms of dual numbers that are free from the mentioned downsides. Here, we extend the idea to develop the formula for all hypercomplex numbers of the second order, and make sure that it replicates the expression of the standard norm of complex numbers (1). The approach is based on the matrix representation $A=\left(\begin{array}{cc}x& y\\ \sigma y& x\end{array}\right)$ of a hypercomplex number $u=x+\tau y$. We derive the norm of this matrix and associate it with the norm of a hypercomplex number. For the matrix norm, we use the following definition:

$${\parallel A\parallel }^2=sup\left\{{\parallel At\parallel }_2^2:t\in {\mathbb{R}}^{2\times 1},{\parallel t\parallel }_2^2=1\right\}$$

After some straightforward mathematical manipulations (see Appendix A), we obtain the following formula for the hypercomplex number norm:

$$\mathbb{H}\parallel u\parallel =\left\{\begin{array}{c}\sqrt{x^2+y^2},u\in \mathbb{C}\hfill \\ \left|\frac{y}{2}\right|+\sqrt{x^2+{\left(\frac{y}{2}\right)}^2},u\in \mathbb{D}.\hfill \\ \left|x\right|+\left|y\right|,u\in \mathbb{S}\hfill \end{array}\right.$$

(2)

To show the differences among the proposed norms for the hypercomplex algebras, we depict $\parallel u\parallel$ equal to 2, 4, or 6 with contour lines in Figure 2. It is shown that the contour lines of dual and double norms are closer to the center than the corresponding contour lines for complex numbers. Different norm expressions have different effects on batch normalization for the dual, complex, and double algebras.

2.4. Batch Normalization

As mentioned in the previous section, batch normalization is very important for a successful training process and convergence to high-metric values. In its original form [14], it can only be applied to the real-valued models. This approach includes the normalization of each dimension of the k-dimensional input:

$${\check{x}}^{\left(k\right)}=\frac{x^{\left(k\right)}-E\left[x^{\left(k\right)}\right]}{\sqrt{Var\left[x^{\left(k\right)}\right]}}.$$

After this transformation, the mean of ${\check{x}}^{\left(k\right)}$ is equal to zero, and the standard variance is equal to one. However, it is better to make the mean and variance trainable parameters. In practice, it is achieved by performing additional scaling and shifting:

$$\mathbb{R}BN\left[{\check{x}}^{\left(k\right)}\right]={\check{\gamma }}^{\left(k\right)}{\check{x}}^{\left(k\right)}+{\check{\beta }}^{\left(k\right)}.$$

The technique of batch normalization has been extended to complex inputs [15]. The authors suggested transforming the inputs using a covariance matrix so that the real and imaginary parts do not correlate. However, we cannot generalize this method explicitly; this is because—for dual numbers—the pseudo-covariance C cannot be equal to zero:

$$\begin{array}{c}\hfill \mathbb{D}C=E\left[{(z-\mu )}^2\right]=E\left[{((z_x-{\mu }_x)+\epsilon (z_y-{\mu }_y))}^2\right]=\\ \hfill =E\left[{(z_x-{\mu }_x)}^2\right]+\epsilon (2·E\left[(z_x-{\mu }_x)(z_y-{\mu }_y)\right]),\end{array}$$

with at least the real parts always being greater than zero. Therefore, we need to find an alternative way to generalize batch normalization to all types of hypercomplex algebras of the second order. In this paper, we present a method based on the norm concept of hypercomplex numbers, as derived in the previous section, and compare the proposed method for complex numbers with the idea presented in [15].

First, we define the mean value in the traditional channel-wise way, which is the same for all algebras:

$${\mu }^{\left(k\right)}=E\left[u^{\left(k\right)}\right]={\mu }_x^{\left(k\right)}+\tau {\mu }_y^{\left(k\right)}.$$

Then we define the variance that is specific to the type of hypercomplex algebra:

$$Var\left[u^{\left(k\right)}\right]=\frac{1}{m-1}\sum _{i=1}^m{\parallel u_i^{\left(k\right)}-{\mu }^{\left(k\right)}\parallel }^2.$$

where $\parallel ·\parallel$ is determined by (2). Then we transform the input as follows:

$${\widehat{u}}^{\left(k\right)}=\frac{u^{\left(k\right)}-{\mu }^{\left(k\right)}}{\sqrt{Var\left[u^{\left(k\right)}\right]+\delta }},$$

where $\delta$ is a small number needed to avoid division by zero. The last step is hypercomplex channel-wise scaling and shifting:

$$\mathbb{H}BN\left[{\widehat{u}}^{\left(k\right)}\right]={\widehat{\gamma }}^{\left(k\right)}{\widehat{u}}^{\left(k\right)}+{\widehat{\beta }}^{\left(k\right)},$$

(3)

where ${\widehat{\gamma }}^{\left(k\right)}$ and ${\widehat{\beta }}^{\left(k\right)}$ are the hypercomplex weights and bias, accordingly.

In the specific case of complex numbers:

$$\mathbb{C}Var\left[z\right]=E\left[(z-\mu ){(z-\mu )}^{\ast }\right]=E[{(x-{\mu }_x)}^2+{(y-{\mu }_y)}^2]=\mathsf{\Gamma }\left(z\right).$$

So, dividing the centered input by $\sqrt{\mathbb{C}Var\left[z\right]}$ provides the covariance $\mathsf{\Gamma }\left(\widehat{z}\right)=1$. Note that, in general, pseudo-covariance is not equal to zero:

$$C\left(\widehat{z}\right)=E\left[{\widehat{z}}^2\right]=E\left[{(\widehat{x}+i\widehat{y})}^2\right]=E[{\widehat{x}}^2-{\widehat{y}}^2]+i(2·E\left[\widehat{x}\widehat{y}\right])\ne 0.$$

Since $\widehat{x}$ and $\widehat{y}$ may correlate, $E\left[\widehat{x}\widehat{y}\right]\ne 0$. Despite that, as shown in Section 4, using this method results in good metric values. We emphasize that MindSpore is the only framework that has implemented this mathematically justified batch normalization technique.

We conduct an ablation study of the proposed batch normalization technique. For this purpose, we construct a hypercomplex LeNet-5 model and vary the normalization. The following approaches are considered: real-valued channel-wise normalization, covariance-based normalization from [15], and algebra-aware normalization (3) based on the proposed hypercomplex norm (2).

We show that all normalization types lead to about the same results on the MNIST dataset after 20 epochs (see

Table 1. Ablation study of normalization techniques.

Type of BN	Algebra
	Complex	Dual	Double
Covariance	98.12	96.88	98.02
Channel-wise	98.53	98.25	98.40
Complex	98.50	98.38	98.38
Dual	98.35	98.23	98.28
Double	98.26	98.26	98.20

). Nevertheless, the proposed hypercomplex normalization is more beneficial because it can be further fused with a previous convolutional layer (Section 3.4).

2.5. Backward Propagation of Loss Function Gradient

An important part of neural network training is gradient calculation. Computing the loss function gradient relies on the backpropagation algorithm, which exploits the chain rule. In this paper, we explore the gradient propagation problem in the algebra of hypercomplex numbers. The classic definition of the derivative of a function $f\left(u\right)=f(x+\tau y)=v(x,y)+\tau w(x,y)$ of a hypercomplex argument $u=x+\tau y$, where v and w are real-valued functions, is

$$f^{\prime }\left(u\right)=\underset{\Delta u\to 0}{lim}\frac{f(u+\Delta u)-f\left(u\right)}{\Delta u}.$$

This limit can only exist if it is well-defined when $\Delta u$ approaches zero along the real axis ($\Delta u=\Delta x$) or the imaginary axis ($\Delta u=\tau \Delta y$). In either case, it should provide the same result. Equating these special cases, we obtain the general equivalents of Cauchy–Riemann equations:

$$\frac{\partial v}{\partial x}=\frac{\partial w}{\partial y},\phantom{AA}\frac{\partial v}{\partial y}={\tau }^2\frac{\partial w}{\partial x}=\sigma \frac{\partial w}{\partial x}.$$

Functions that satisfy these equations are called holomorphic. In practice, this is a strong restriction, and most of the existing operators do not satisfy the Cauchy–Riemann criteria. To overcome this, we use the approach invented by Wirtinger for complex numbers [17]. It applies variable substitution in order to rewrite a complex-variable function $f\left(z\right)$ as a holomorphic function of two arguments $f(z,z^{\ast })$. This idea was applied in [18] to derive the gradient of double numbers. Here, we extend this approach to all second-order algebras:

$$x=\frac{u+u^{\ast }}{2},y=\frac{u-u^{\ast }}{2\tau }.$$

For complex and double numbers, we can easily eliminate $\frac{1}{\tau }$ by multiplying both the numerator and denominator by $\tau$ and we end up with ${\tau }^2=\sigma \in \mathbb{R}$ in the denominator. However, we cannot perform the same trick for dual numbers because—in that case—$\sigma$ is equal to zero. Thereby, we must keep $\epsilon$ in the denominator. We must also keep the expressions with ${\epsilon }^2$ in the numerator, because they will contribute to the component with $\epsilon$ if they cancel one of the $\epsilon$ with $\frac{1}{\epsilon }$.

Here, we clarify what exactly must be calculated in order to update hypercomplex weights. With complex-valued networks, researchers usually treat the real $x_n$ and imaginary $y_n$ parts of weights $z_n$ as separate real-valued channels and update them using the real-valued derivative of a loss function:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{x}_{n+1}=x_n-\alpha \frac{\partial L}{\partial x}\hfill \\ y_{n+1}=y_n-\alpha \frac{\partial L}{\partial y}\hfill \end{array}\Rightarrow z_{n+1}=z_n-\alpha \left(\frac{\partial L}{\partial x}+i\frac{\partial L}{\partial y}\right)\right.\end{array}$$

We generalize these calculations for all hypercomplex numbers of the second order, considering:

$$u_{n+1}=u_n-\alpha ·\left(\frac{\partial L}{\partial x}+\tau \frac{\partial L}{\partial y}\right)$$

We define the expression inside the parenthesis as a hypercomplex gradient and calculate it via the gradient of a hypercomplex operator f (see Appendix B) as follows:

$$\frac{\partial L}{\partial x}+\tau \frac{\partial L}{\partial y}=\frac{\partial L}{\partial f}\left(\frac{\partial f}{\partial x}+\tau \frac{\partial f}{\partial y}\right)+{\left(\frac{\partial L}{\partial f}\right)}^{\ast }{\left(\frac{\partial f}{\partial x}-\tau \frac{\partial f}{\partial y}\right)}^{\ast }$$

It is noteworthy that $\frac{\partial L}{\partial x}+\tau \frac{\partial L}{\partial y}$ can be expressed in terms of u and $u^{\ast }$, i.e.:

$$\frac{\partial f}{\partial x}+\tau \frac{\partial f}{\partial y}=\left\{\begin{array}{c}2·\frac{\partial f}{\partial u^{\ast }},u\in \mathbb{C}\hfill \\ 2·\frac{\partial f}{\partial u},u\in \mathbb{S}\hfill \\ \frac{\partial f}{\partial u}+\frac{\partial f}{\partial u^{\ast }}+{\epsilon }^2\left(\frac{\partial f}{\partial u}-\frac{\partial f}{\partial u^{\ast }}\right),u\in \mathbb{D}\hfill \end{array}\right.$$

With dual numbers, we keep the part with ${\epsilon }^2$, despite the fact that it must be equal to zero, as stated above.

We implemented this formula for dual algebra. This method yields results consistent with the calculations from two real-valued derivatives. Our experiments show that the training time is approximately the same as well. However, the performance can be improved if operations on hypercomplex numbers are implemented as instructions on the hardware level.

3. Optimizations

In this section, we discuss technical methods to speed up the execution of testing and/or training models based on hypercomplex numbers. While this is an important engineering challenge, it is not as crucial as achieving the best possible accuracy. We focus on searching for equivalent transformations that speed up hypercomplex operators and do not affect the network output or metric values. These methods are not universal, so the final acceleration depends on the model architecture.

3.1. Hypercomplex Convolution Grouping

We express a hypercomplex convolution through a certain number of real-valued convolutions (Section 2.2). For complex- and dual-valued convolutions, we have to conduct extra addition and/or subtraction (see Figure 1). It is possible to “insert” these arithmetic operations into the real-valued convolutions. In order to do that, we concatenate the components of the input and perform a similar operation for the weights (see Figure 3). This approach increases memory usage and necessitates tensor duplication, using the best-optimized solution. Note that, in the case of the complex-valued convolution, we need to invert the imaginary part of the weights when we calculate the real part of the output. Thus, we implement subtraction. Although the negation operation must be cheap in terms of execution time, it is possible, if needed, to store the inverted weights in advance to speed up the inference at the cost of a larger model size.

3.2. Karatsuba’s Algorithm

Karatsuba’s algorithm is a well-known technique used for the optimization of complex number multiplication. It can also be applied to convolutions and linear complex-valued operators. For example, a linear operator can be rewritten as three real-valued linear operators:

$$\mathbb{C}L(w_x+i{w}_y,x+iy,b_x+i{b}_y)=L_1-L_2+i(L_3-L_1-L_2),$$

where $L_1=\mathbb{R}L(w_x,x,b_x)$, $L_2=\mathbb{R}L(w_y,y,0)$, and $L_3=\mathbb{R}L(w_x+w_y,x+y,b_x+b_y)$. Fewer matrix multiplications come at the expense of a greater number of additions and subtractions. Therefore, Karatsuba’s algorithm performs better for convolutions with large batch sizes and weight widths, where multiplications are more expensive.

3.3. Diagonal Representation of Double Numbers

In Section 1, we briefly mentioned the representation of double numbers in the diagonal basis. Here is its mathematical definition:

$$\begin{array}{c}\hfill x^{\prime }=x+y,y^{\prime }=x-y\\ \hfill z=x+jy=x^{\prime }\frac{1+j}{2}+y^{\prime }\frac{1-j}{2}\end{array}$$

The computation complexity of multiplication is halved in this form:

$$(A\frac{1+j}{2}+B\frac{1-j}{2})·(x\frac{1+j}{2}+y\frac{1-j}{2})=Ax\frac{1+j}{2}+By\frac{1-j}{2}$$

Section 2.2 shows the isomorphism between hypercomplex numbers and matrices of the second order. Then, a convolution of double numbers in the diagonal form can be expressed as follows:

$$\begin{array}{c}\hfill W\ast u\sim \frac{1}{2}\left(\begin{array}{cc}{W}_x+W_y& W_x-W_y\\ W_x-W_y& W_x+W_y\end{array}\right)\ast \frac{1}{2}\left(\begin{array}{cc}x+y& x-y\\ x-y& x+y\end{array}\right)=\\ \hfill =\frac{1}{2}\left(\begin{array}{cc}{W}_x\ast x+W_y\ast y& W_x\ast x-W_y\ast y\\ W_x\ast x-W_y\ast y& W_x\ast x+W_y\ast y\end{array}\right),\end{array}$$

which implies the convolution calculation as two component-wise real-valued convolutions $W_x\ast x$ and $W_y\ast y$. Thus, we define new variables as $j_1=\frac{1+j}{2}$ and $j_2=\frac{1-j}{2}$, so that $z=j_{1}x+j_{2}y$. Batch normalization can be updated accordingly. Moving from the regular representation to the diagonal one is identical to the affine transformation of the plane of double numbers. The norm of a double number $z=j_{1}x+j_{2}y$ is expressed as $max\left(\right|x|,|y\left|\right)$.

As discussed before, the downside of using the diagonal representation is that x and y are poorly interlinked and are almost equivalent to two independent real numbers. The batch normalization based on the norm of double numbers partially solves this issue. Average pooling and linear functions are calculated for both components independently, because they are special cases of convolution. The relationship between two components is introduced through non-linear activation functions. The ReLU implementation

$$\mathbb{S}ReLU(x+jy)=\mathbb{R}ReLU\left(x\right)+j\mathbb{R}ReLU\left(y\right)$$

works well with the classic representation of double numbers. It makes sense to reuse this by transforming the intermediate tensor to the classic representation before applying the component-wise activation. In the end, the formula looks as follows:

$$\begin{array}{c}\hfill \mathbb{S}ReLU(j_{1}x+j_{2}y)=j_1\left[\mathbb{R}ReLU\left(\frac{x+y}{2}\right)+\mathbb{R}ReLU\left(\frac{x-y}{2}\right)\right]+\\ \hfill +j_2\left[\mathbb{R}ReLU\left(\frac{x+y}{2}\right)-\mathbb{R}ReLU\left(\frac{x-y}{2}\right)\right].\end{array}$$

This activation function requires more operations compared to the base ReLU implementation. However, in combination with the enhanced performance of the convolution, it allows the model to save precision with significant acceleration.

3.4. Equivalent Transformations

Real-valued convolutional neural networks can be accelerated by means of a fusing convolution with a subsequent batch normalization layer into a convolution with new parameters, as described by the following equations:

$$W_{fused}=\frac{\alpha }{\sqrt{Var}}\times W,$$

(4)

$$b_{fused}=\frac{\alpha }{\sqrt{Var}}\times (b-\mu )+\beta ,$$

(5)

where $W_{fused}$ and $b_{fused}$ are weights and biases of the fused convolution, W and b are weights and biases of the original convolution, $\mu$, $Var$, $\alpha$, $\beta$ are the mean value, variance, scale, and shift of the batch normalization layer, respectively. The × operator is an element-wise multiplication.

This transformation leads to improved inference time and memory conservation. In the case of real models, some frameworks have already implemented this optimization. For example, the ONNX simplifier does that automatically for models in the ONNX graph format.

Hypercomplex operators require more resources, which makes this optimization even more critical. However, so far, there is no automated tool that could apply it to hypercomplex models. For example, an attempt to export a model to ONNX and use the ONNX simplifier would fail. Currently, the ONNX simplifier cannot fuse a sequence consisting of hypercomplex convolutional and batch normalization layers. So, we introduce a new layer called CBN (Conv+BN), which behaves differently depending on the mode. During training, it behaves as the original sequence of two layers, and in the evaluation mode, it works as a convolution with fused parameters (Figure 4). These parameters are pre-calculated according to Formulas (4) and (5), assuming that all multiplications follow the rules of arithmetic of hypercomplex numbers.

Our CBN layer is implemented in MindSpore to improve the inference times of hypercomplex models.

4. Results

To explore the generalization of neural networks on hypercomplex numbers, we develop a demonstrative network to predict values of noised hypercomplex functions. In addition, we conduct a series of experiments with hypercomplex models for computer vision and signal processing tasks. We compare the results of these networks with the results of real-valued models of the same architecture.

4.1. Hypercomplex Toy Net

We start with a demonstrative neural network designed for predicting the values of noisy functions of a hypercomplex argument. Our toy network consists of two linear layers and a sigmoid-type activation function (see Appendix C) in between. This architecture remains the same for all the algebras, but for every algebra, we use its implementation of the operators.

In order to show the advantage of hypercomplex models, we compare their results with the result of a real-valued model of the same architecture. The mean square error (MSE) is used as a loss function. We train these four models for 1000 epochs with the same number of parameters to predict values of two functions: $Ai\left(u\right)$—the Airy function of the first kind and $J_3\left(u\right)$—the Bessel function of the first kind of the third order. We also add noise with a normal distribution to the values of functions in a training set (and a testing set). The results are shown in

Table 2. Toy net—value of the loss function.

Algebra	$\mathit{Ai}\left(\mathit{u}\right)$			${\mathit{J}}_3\left(\mathit{u}\right)$
	$\mathit{u}\in \mathbb{C}$	$\mathit{u}\in \mathbb{D}$	$\mathit{u}\in \mathbb{S}$	$\mathit{u}\in \mathbb{C}$	$\mathit{u}\in \mathbb{D}$	$\mathit{u}\in \mathbb{S}$
Real	0.026	0.045	0.013	0.018	0.017	0.015
Complex	0.008	0.015	0.015	0.009	0.012	0.011
Dual	0.017	0.009	0.017	0.013	0.009	0.012
Double	0.023	0.050	0.010	0.017	0.018	0.010

.

It is clear from Figure 5 and Table 2 that the lowest mean square error is achieved by a neural network belonging to the same type of hypercomplex number as the function argument. For example, a complex-valued function is best predicted by a complex-valued model, and so on. Thus, we conclude that networks based on hypercomplex numbers can learn dependencies or patterns that lie in the corresponding algebra.

4.2. Classical CV Problems

Before we move on to solving CV classification problems, it is necessary to clarify the process of data preparation. To convert a real-valued image to the hypercomplex format in this research, we use a more general method compared to the one proposed in [5]. In [5], the authors convert a real-valued image to a complex form by $[R,G,B]\Rightarrow [R+iG,G+iB]$, claiming that this type of encoding captures the channel correlations and hue shift. In this paper, we generalize this idea and use the following transformation $[R,G,B]\Rightarrow [R+\tau G,G+\tau B,B+\tau R]$, which we call the color combination. However, there is another possibility to generate an imaginary component by processing the real input data using several layers [15].

We take ResNet18 architecture [19] as the basis of our model that is generalized to hypercomplex algebras of the second order by replacing real-valued operators with their counterparts in hypercomplex algebras. We use stochastic gradient descent with 0.9 momentum to optimize real-valued loss functions by treating the real and hypercomplex parts as separate real-valued channels. To transform hypercomplex features into the real ones, we concatenate both components of the output into a double-sized real-valued tensor before feeding it to the terminating fully connected layer. The cross-entropy loss between the input and the target is used as a criterion for these problems. The total number of epochs is 200 for every experiment. The learning rate is scheduled by the cosine annealing rule, starting from $0.1$. The size of a mini-batch is 8. Every model is trained 5 times with different seed numbers, and the median value of the accuracy is taken to negate the effects of fluctuation. Our image classification results on CIFAR-10, CIFAR-100, and SVHN are presented in

Table 3. The accuracy (%) of our models for classical CV problems.

Algebra	Dataset
	CIFAR-100	CIFAR-10	SVHN
Real	78.63	95.33	96.47
Dual	79.09	95.51	96.60
Complex	79.23	95.66	96.60
Double	78.88	95.51	96.45

.

One can observe in Table 3 that all models based on the second-order numbers achieve higher accuracy compared to the real one. In addition, the complex-valued neural network shows better metric values compared to other hypercomplex models.

Table 4. Effect of batch normalization fusion on the ResNet18 inference time (ms).

Algebra	Device
	CPU		GPU
	No Fusion	With CBN	No Fusion	With CBN
Real	8.56	6.54	1.72	1.42
Complex	34.97	26.49	6.73	5.13
Dual	28.00	20.78	5.45	4.17
Double	34.83	26.62	6.67	5.08

Software as follows: Python 3.8.10, MindSpore 1.10.0; Hardware as follows: CPU: Intel(R) Core (TM) i9-11900KF @ 3.50 GHz, GPU: NVIDIA GeForce RTX 3090.

demonstrates how remarkably inference time is reduced by using the CBN equivalent transformation described in Section 3.4. The effect is so dramatic that, from now on, it is worth considering models with CBN only. Figures from the CBN column of Table 4 are used as the baseline for the assessment of other optimizations, which are discussed below. These figures also indicate that transferring from a real model to a hypercomplex one significantly increases computational complexity.

Table 5. Average inference time (ms) of ResNet18 models.

Algebra	Optimization	GPU Inference Time, ms
Real	Baseline	1.42
Dual	Baseline	4.17
	Conv Grouping	4.04
Complex	Baseline	5.13
	Karatsuba	4.79
	Conv Grouping	4.91
Double	Baseline	5.08
	Diagonal	3.29

Software as follows: Python 3.8.10, MindSpore 1.10.0; Hardware as follows: CPU: Intel(R) Core (TM) i9-11900KF @ 3.50 GHz, GPU: NVIDIA GeForce RTX 3090.

shows the effect of optimization on the second-order hypercomplex models.

The worst performance is shown by the naive double-valued network, but the transition to the diagonal representation accelerates it dramatically. Still, the double-valued network is not an optimal model because it shows the worst accuracy among all hypercomplex models (see Table 3). As for the complex- and dual-valued networks, grouping convolutions helps to improve the inference time, and for the complex-valued network, the implementation of Karatsuba’s algorithm helps even more.

Our main idea is that dual-valued neural networks, as well as complex-valued ones, allow different channels to interact with each other due to the internal connection between real and dual parts. This inner complexity should be exploited; possibly, more sophisticated input preprocessing can lead to even better accuracy for these more complicated networks.

4.3. Gravitational Waves Detection

This part is dedicated to the solution of a signal detection task using hypercomplex networks for the G2Net dataset. This dataset consists of simulated noised signals that are similar to gravitational waves recorded by a system of three ground-based laser interferometers: LIGO Hanford, LIGO Livingston, and Virgo [20,21,22]. Generally, gravitational waves are emitted during cosmic events, such as the merging of black holes [20]. The G2Net dataset contains records of emulated events of the same nature.

To classify the original signal, we preprocess data to an image and run it through our neural networks. The aim of preprocessing is to build a representative frequency map of the original signal. The CQT algorithm [23] is considered efficient for analyzing gravitational waves.

This preprocessing transfers the time series into the frequency portrait (see Figure 6), which is treated as an image in further steps. Thus, the variations of frequency characteristics of the signal at certain moments are reflected as visual features, such as specific shapes and colors, in the resulting image. To classify the image obtained after CQT, we again use the ResNet18 model, whose operators are changed to the corresponding ones in ad hoc algebras.

The model is optimized using the stochastic gradient descent algorithm with a momentum of 0.9 and weight decay of $5\times {10}^{-5}$ for $L_2$ regularization. The regularization is enhanced by adding dropout with a value of 0.2. The scheduler for the learning rate uses the cosine annealing policy with $\mathit{Tmax}=6$, and the initial value of lr is set to $5\times {10}^{-2}$. Training lasts 30 epochs.

From

Table 6. Average values of metrics (%) of our models for gravitational wave detection.

Algebra	Accuracy, %	AUC ROC
Real	76.45	0.82
Complex	78.73	0.84
Dual	79.24	0.85
Double	77.41	0.84

, one can see that all models on the hypercomplex algebra outclass the real-valued model in both accuracy and AUC ROC (area under the receiver operating characteristic curve) metrics. The best result is achieved by a dual-valued neural network.

4.4. Music Transcription Task

In this section, we present the results of training hypercomplex networks on an automatic music transcription task. The experiments are performed on the MusicNet dataset [24]. In order to convert the input signal into a complex-valued tensor, we use the Fourier transformation. We reuse the DeepConvNet architecture developed in [15], with its building blocks replaced with the corresponding layers designed for specific algebras.

Table 7. DeepConvNet—average precision.

Algebra	Avg Precision, %
Real	70.0
Complex	73.2
Dual	73.4
Double	73.2

shows that hypercomplex networks significantly outclass the same model on real numbers. The best figures are shown by a dual-valued network.

Next, we present the results of inference time measurement on GPU (see

Table 8. DeepConvNet—average GPU inference time (ms).

Algebra	Not Optimized, ms	Optimized, ms
Real	0.66	0.53
Complex	2.26	1.39
Dual	1.88	1.36
Double	2.26	1.53

Software as follows: Python 3.8.10, MindSpore 1.10.0; Hardware as follows: CPU: Intel(R) Core (TM) i9-11900KF @ 3.50 GHz, GPU: NVIDIA GeForce RTX 3090.

). The first column presents the measured inference time of the baseline model, and the second one shows the inference time of the same models after optimizations.

The optimizations of the complex-valued network include convolution grouping and Karatsuba’s algorithm for the most expensive linear operator. For dual numbers, we only use convolution grouping, which provides a small performance gain. The double-valued operators are redesigned to take advantage of the diagonal representation. In all cases, batch normalization operators are fused with the preceding convolutions. The experiments are performed with a batch size of 8. The dual-valued model outperforms all other second-order hypercomplex models. The double-valued network, despite its diagonal representation, yields the poorest results due to the significant overhead associated with activation functions.

5. Discussion

5.1. Summary

Improving model accuracy is the cornerstone of many modern research studies in the field of deep learning. This target can be achieved by different methodologies, such as dataset expansion, changing the training procedure, variation of neural network architecture, and so on. In this study, we offer a take-and-go recipe on how to enhance the model’s generalization ability. Namely, one needs to replace real-valued layers with the corresponding complex, dual, or double counterparts, without any modifications to the overall architecture. Generally speaking, it is impossible to predict which kind of algebra will be better for a particular task. To overcome this obstacle, we define a universal approach to developing operators based on any kind of hypercomplex algebra of the second order. This allows researchers to design an architecture of neural networks without heavy dependence on a specific kind of algebra.

In our experiments on the CV, G2Net, and MusicNet tasks, complex and dual-valued models have shown the best results. Nevertheless, neural network construction is a multi-criteria problem. If inference and training times are important, then dual-valued models are often preferable, whereas the complex-valued ones can provide extra accuracy when the timing is not so critical.

5.2. Limitations and Outlook

A limitation of the present study is our sole method of converting real-valued data to hypercomplex numbers for CV problems. It is possible that suboptimal data preprocessing resulted in small gains on these tasks. We plan to explore new methods of real input transformation.

Another unsolved issue is how to identify the type of algebra—for a particular task—that would have the best performance among the hypercomplex algebras of the second order. We will be looking for some heuristic methods that predict the best algebra without training the models with all types of algebras.

In conclusion, the presented methodology exhibits significant promise in enhancing the efficiency of transformer-based architectures. Also, we plan to develop a hypercomplex network compression methodology because we did not find any compression algorithms for non-real models.