\chapter{Deep Cyclic-group Networks}
\label{chp:dcgn_method}

The peculiarity of CNNs is that it learns input invariance by using kernel maps and outputs feature maps as detected features. Elements in a feature map and a kernel map represent hidden neurons and weights respectively. It uses the same weights for each neuron. Thus, CNNs can be regarded as a kind of DNN architecture using sharing weights for each of hidden neurons in the same layer.
In previous work, DCN and DQN extend each neuron and weight to two- and four-dimensions respectively, using complex and quaternion algebras. However, they lack the ability to deal with input data of arbitrary dimensions. The idea of our proposed DCGN is to extend the capacity of each neuron from scalars to vectors with cyclic group algebras. Here, the scalar multiplications in the conventional convolutional procedure are changed into vector multiplications along mode-$3$ fibers with the help of circular convolutions, which can also be formulated as vector-matrix multiplications with circulant matrices~\cite{farenick01algebras,davis12circulant}. In other words, our work replaces scalar multiplications with vector multiplications through circulant matrices.


In DCGNs, both the kernel and feature maps are three-way tensors comprising $N$-dimensional mode-$3$ fibers. Each fiber in a feature map stands for a $N$-dimensional vector-valued neuron and its output is also a $N$-dimensional fiber. The circular convolution is operated between a vector-valued neuron and a fiber coming from a weight kernel map. The convolution procedure is illustrated in Figure~\ref{fig:dcgn_explain}. 
 

In what follows, we firstly derive the learning algorithm from the point of view of a vector neuron,  and then describe the overall DCGN model for regression problems and its extension to classification problems. Below, matrices are represented by bold uppercase letters, and vectors by bold lowercase letters. Matrix slices of tensors are represented as matrices and vector slices of matrix are represented as vectors. We provide the simplified version of backpropagation learning algorithm and treat vectors as row vectors. Detailed version will be released along with the code.  

\begin{figure}[t]
\begin{center}
\includegraphics[trim=9cm 6cm 9cm 3cm, clip, scale=0.55]{ch3_dcgn_method/dcgn_explain.pdf}
\caption{Illustration of the convolution procedure when using cyclic group algebra (best viewed in color). 
Each cube represents a scalar. 
Collectively, the $N$ matrices of cubes represent a three-way tensor of mode-$3$ dimension $N$, each matrix is a frontal slice of the tensor, and each mode-$3$ fiber is an $N$-dimensional vector. 
The feature map and the kernel map are composed of blue cubes and green cubes, respectively. 
In the convolution procedure, each mode-$3$ fiber of the kernel map performs circular convolution with a mode-$3$ fiber of the feature map, by rotating itself along the mode-$3$ dimension, as the green textured cubes represent. Therefore, via the kernel map we can learn the association between the elements across the $N$ dimensions
by using a linear kernel and weighted sum. 
The kernel map moves along the (mode-$1$, mode-$2$) plane, each time at a distance of a stride.}
\label{fig:dcgn_explain}
\end{center}
\end{figure}

%---------------------------------------------------------------
\section{Operation of Circular Convolutions}
\label{sec:ch3_operation_of_circular_convolution}
To perform the learning process between a three-order tensor feature map and a three-order tensor kernel map, we leverage the circular convolution to be an basic operator in the network. Given two $N$-dimensional vectors, $\mathbf{p}=\left[p_1\textrm{,}~p_2\textrm{,} \dots\textrm{,} p_N\right]  \in \mathbb{R}^N$ and $\mathbf{q}=[q_1\textrm{,}~q_2\textrm{,} \dots\textrm{,}q_N]  \in \mathbb{R}^N$,
we can compute the circular convolution between them, denoted as $\mathbf{p}\ast\mathbf{q}$, with vector-matrix multiplication by using the circulant matrix $[\mathbf{q}]_{\ast} = \text{circ}(\mathbf{q})$ associated to $\mathbf{q}$ as follows:
\begin{equation}
\renewcommand{\arraystretch}{0.60}
\label{eq:circular_convolution_formulation}
\setlength{\arraycolsep}{3pt}
\mathbf{p}\ast\mathbf{q}
=\mathbf{p}[\mathbf{q}]_{\ast} =
\begin{bmatrix}
p_{1}  \\
p_{2}  \\
p_{3}  \\
\vdots \\
p_{N}
\end{bmatrix}^T
\begin{bmatrix}
q_{1}    & q_{2}   & q_{3}    & \dots     & q_{N}      \\ 
q_{N}    & q_{1}   & q_{2}    & \dots     & q_{N-1}    \\
q_{N-1}  & q_{N}   & q_{1}    & \dots     & q_{N-2}    \\
\vdots   & \vdots  & \vdots   & \ddots    & \vdots     \\
q_{2}    & q_{3}   & q_{4}    & \dots     & q_{1}      \\
\end{bmatrix} \,.
\end{equation}
We note that $[\mathbf{q}]_{\ast} \in \mathbb{R}^{N \times N}$ and $\mathbf{p}\ast\mathbf{q} \in \mathbb{R}^N$.
Algebra of circulant matrices is the vector's isomorphism. The result of $\mathbf{p}\ast\mathbf{q}$ is identical to the definition of circular convolution operation. In our proposed model, the feedforward and backpropagation processes will be described by using above the formulation and the matrix representation.



\section{Feedforward process for DCGN}
\label{sec:ch3_feedforward}
In vector neural learning with cyclic group algebras, the inputs, outputs, weights and biases are all $N$-dimensional vectors. 
Specifically, for a network with $L$ layers, the relation between the output of a neuron $i$ in layer $l$ ($1\leq l\leq L$), $\mathbf{a}_i^{l} \in \mathbb{R}^N$, and the output of a neuron $j$ in the preceding layer $l-1$, $\mathbf{a}_j^{l-1} \in \mathbb{R}^N$, can be written as:
%Suppose there is a $L$-layer network. For the hidden layer $l$ $(1 \leq l \leq L)$, the input vector $\mathbf{z}_i^l$ of the hidden neuron $i$ is defined as:
\begin{equation}
\label{eq:feedforward_equation}
\mathbf{a}_i^l\triangleq\phi(\mathbf{z}_i^l)=\phi\left(\sum_{j=1}^{J^{l-1}}\mathbf{w}_{ij}^{l}\ast\mathbf{a}_j^{l-1}+\mathbf{b}_{i}^{l}\right)~\textrm{,}
\end{equation}
where $\phi(\cdot)$ is a differentiable activation function, $\mathbf{w}_{ij}^l \in \mathbb{R}^N$ the weight vector connecting the two neurons, $J^{l-1}$ the number of neuron in layer $l-1$, and $\mathbf{b}_{i}^{l} \in \mathbb{R}^N$ a bias vector.
%We also call the intermediate result $\mathbf{z}_i^{l} \in \mathbb{R}^N$ an input vector of layer $l$.
For $l=1$, we have $J^{l-1}=1$ and $\mathbf{a}^{l-1}\triangleq\mathbf{x} \in \mathbb{R}^N$, %=\mathbf{a}^{0}
the input data vector of the network.
For $l=L$, we have $J^{l}=1$, and we define $\mathbf{y} \in \mathbb{R}^N$ as $\mathbf{a}^{L}$, %(dropping the subscript) 
the output of the network.




%---------------------------------------------------------------
\section{Backpropagation algorithm for DCGN}
\label{sec:ch3_backpropagation}
During the training phase, we are given the target output $\hat{\mathbf{y}} \in \mathbb{R}^N$ for every input data instance, so that we can compare $\hat{\mathbf{y}}$ and $\mathbf{y}$ via some cost function $\epsilon(\mathbf{y},\hat{\mathbf{y}})$ to know how good the model parameters $\mathbf{\Theta}\triangleq\{\mathbf{w}_{ij}^l, \mathbf{b}_{i}^{l}\}$ are. Here, we use $\mathbf{\Theta}$ to denote the collection of all the trainable parameters (i.e., weights and biases) 
across the $L$ layers. Given the empirical cost $C$ calculated over all the training data instances, we can use backpropagation to update the model parameters $\mathbf{\Theta}$. Finally, given the training data and the cost function $C(\mathbf{\Theta})$, the neural network is used to perform nonlinear mapping between the inputs and the outputs. The objective of the training process is to estimate parameters minimizing the cost function.

From the previous sections, the feedforward process of the proposed architecture operates with the circular convolution, which is formulated as matrix-vector multiplication. In the backpropagation process, the learning algorithm is derived in the same way. To use backpropagation, we need to take the derivative of $C$ with respect to each parameter $\mathbf{w}_{ij}^l$ or $\mathbf{b}_{i}^{l}$. Unlike ordinary CNNs, all the parameters here are vectors. Therefore, in calculating the derivatives we need to compute the Jacobian matrix. For the bias vectors $\mathbf{b}_{i}^{l}$, this amounts to calculating the $N$-dimensional gradient vector: 
%\renewcommand{\arraystretch}{2}
\begin{equation}
\label{eq:derivative_b_fun}
\frac{\partial C}{\partial \mathbf{b}_{i}^l}=\left[\frac{\partial C}{\partial b_{i1}^l}, \frac{\partial C}{\partial b_{i2}^l} \textrm{,}~\dots, \frac{\partial C}{\partial b_{iN}^l}\right]\textrm{,}
\end{equation}
where 
\begin{equation}
\frac{\partial C}{\partial b_{ik}^l}= \sum_{n=1}^N \frac{\partial C}{\partial z_{in}^l} \frac{\partial z_{in}^l}{\partial b_{ik}^l}\textrm{.}
\end{equation}
It can be shown that $\frac{\partial z_{in}^l}{\partial b_{ik}^l}=1$ only when $n=k$, otherwise it is zero.
It will become apparent later that it is more convenient to denote $\frac{\partial C}{\partial z_{in}^l}$ simply as $d_{in}^l$ and accordingly define the \emph{local gradient vector} $\mathbf{d}_{i}^{l} =\left[d_{i1}^l, d_{i2}^l,\dots,d_{iN}^l\right]$ as follows:
\renewcommand{\arraystretch}{1}
\begin{equation}
\label{eq:delta_l}
\mathbf{d}_{i}^{l} 
\triangleq \frac{\partial C}{\partial \mathbf{z}_{i}^l}
=\left[
\displaystyle \frac{\partial C}{\partial z_{i1}^l}\textrm{,}~
\displaystyle \frac{\partial C}{\partial z_{i2}^l}\textrm{,}~
\dots \textrm{,}~
\displaystyle \frac{\partial C}{\partial z_{iN}^l}
\right]
\,.
\end{equation}
The terms $\left\{\frac{\partial z_{i1}^{l}}{\partial b_{ik}^{l}}\textrm{,}~\frac{\partial z_{i2}^{l}}{\partial b_{ik}^{l}}\textrm{,}~\dots\textrm{,}~\frac{\partial z_{iN}^{l}}{\partial b_{in}^{l}}\right\}$ ($1 \leq k \leq N$) can be obtained by differentiating Eq.~\ref{eq:feedforward_equation}:
:
\renewcommand{\arraystretch}{1.2}
\begin{equation}
\label{eq:partial_z_il_b_in}
\begin{split}
\frac{\partial \mathbf{z}_{i}^l}{\partial b_{in}^l}
&=
\begin{bmatrix}
\displaystyle
\frac{\partial z_{i1}^l}{\partial b_{ik}^l} \textrm{,}~ 
\frac{\partial z_{i2}^l}{\partial b_{ik}^l} \textrm{,}~\dots\textrm{,}~
\frac{\partial z_{iN}^l}{\partial b_{ik}^l}
\end{bmatrix}^T \\
&=
\begin{bmatrix}
\displaystyle
\frac{\partial b_{i1}^l}{\partial b_{ik}^l} \textrm{,}~ 
\frac{\partial b_{i2}^l}{\partial b_{ik}^l} \textrm{,}~\dots\textrm{,}~
\frac{\partial b_{iN}^l}{\partial b_{ik}^l}
\end{bmatrix}^T \\
&=
\underbrace{
\begin{bmatrix}
0 \textrm{,}~
\dots \textrm{,}~
1 \textrm{,}~
\dots \textrm{,}~
0 
\end{bmatrix}^T 
}_{\textrm{the }n \textrm{-th entry is }1}\textrm{.}
\end{split}
\end{equation}
It can be shown that $\frac{\partial C}{\partial b_{ik}^l}=d_{ik}^l$ and accordingly, 
%that the gradient in Eq.~(\ref{eq:derivative_b_fun}) is equal to the local gradient vector:
\renewcommand{\arraystretch}{2}
\begin{equation}
\label{eq:partial_C_over_partail_bi}
\frac{\partial C}{\partial \mathbf{b}_{i}^l}
= \left[\dots\,, \sum_{n=1}^N \frac{\partial C}{\partial z_{in}^l} \frac{\partial z_{in}^l}{\partial b_{ik}^l}\,, \dots \right]
=~ \mathbf{d}_{i}^{l}  \,.
\end{equation}
As a consequence, it is obvious that the derivative of $\mathbf{b}_{i}^l$ is equal to local gradient vector $\mathbf{d}_{i}^{l}$. Next, from the chain rule, we can calculate the derivative of $C$ w.r.t.~the weight vector $\mathbf{w}_{ij}^l$ similarly. With some algebra, it can be shown that
% \renewcommand{\arraystretch}{2}
% \label{eq:derivative_w_fun}
% \setlength{\arraycolsep}{2pt}
% \begin{equation}
% \begin{split}
% \frac{\partial C}{\partial \mathbf{w}_{ij}^l}
% &=\left[\dots,\frac{\partial C}{\partial w_{ijk}^l},\dots\right] \\
% &=\left[\dots,
% \sum_{n=1}^N \frac{\partial C}{\partial z_{in}^l} \frac{\partial z_{in}^l}{\partial w_{ijk}^l},\dots\right] \\
% &=\mathbf{d}_{i}^l\displaystyle \frac{\partial \mathbf{z}_{i}^l}{\partial \mathbf{w}_{ij}^l}
% = \mathbf{d}_{i}^{l}[\mathbf{a}_{j}^{l-1}]_{\ast}^T\,.
% \end{split}
% \end{equation}
\renewcommand{\arraystretch}{1.05}
\setlength{\arraycolsep}{2pt}
\begin{equation}
\label{eq:derivative_w_fun}
\begin{split}
\frac{\partial C}{\partial \mathbf{w}_{ij}^l}
&=
\begin{bmatrix}
\displaystyle
\frac{\partial C}{\partial w_{ij1}^l} \textrm{,}~ 
\frac{\partial C}{\partial w_{ij2}^l} \textrm{,}~\dots\textrm{,}~
\frac{\partial C}{\partial w_{ijN}^l}
\end{bmatrix} \\
&=
\begin{bmatrix}
\displaystyle \frac{\partial C}{\partial z_{i1}^l} \frac{\partial z_{i1}^l}{\partial w_{ij1}^l}+  
\displaystyle \frac{\partial C}{\partial z_{i2}^l} \frac{\partial z_{i2}^l}{\partial w_{ij1}^l}+ 
\displaystyle \dots+
\displaystyle \frac{\partial C}{\partial z_{iN}^l} \frac{\partial z_{iN}^l}{\partial w_{ij1}^l} \\ 
%------------------
\displaystyle \frac{\partial C}{\partial z_{i1}^l} \frac{\partial z_{i1}^l}{\partial w_{ij2}^l}+  
\displaystyle \frac{\partial C}{\partial z_{i2}^l} \frac{\partial z_{i2}^l}{\partial w_{ij2}^l}+ 
\displaystyle \dots+
\displaystyle \frac{\partial C}{\partial z_{iN}^l} \frac{\partial z_{iN}^l}{\partial w_{ij2}^l} \\
%------------------
\displaystyle \vdots   \\
%------------------
\displaystyle \frac{\partial C}{\partial z_{i1}^l} \frac{\partial z_{i1}^l}{\partial w_{ijN}^l}+  
\displaystyle \frac{\partial C}{\partial z_{i2}^l} \frac{\partial z_{i2}^l}{\partial w_{ijN}^l}+
\displaystyle \dots+
\displaystyle \frac{\partial C}{\partial z_{iN}^l} \frac{\partial z_{iN}^l}{\partial w_{ijN}^l} \\
\end{bmatrix}^T 
= \mathbf{d}_{i}^{l}
\begin{bmatrix}
  \displaystyle  \frac{\partial z_{i1}^l}{\partial w_{ij1}^l}
& \displaystyle  \frac{\partial z_{i1}^l}{\partial w_{ij2}^l}
& \displaystyle  \dots
& \displaystyle  \frac{\partial z_{iN}^l}{\partial w_{ijN}^l} \\ 
%------------------
  \displaystyle  \frac{\partial z_{i2}^l}{\partial w_{ij1}^l}  
& \displaystyle  \frac{\partial z_{i2}^l}{\partial w_{ij2}^l} 
& \displaystyle  \dots
& \displaystyle  \frac{\partial z_{iN}^l}{\partial w_{ijN}^l} \\
%------------------
\vdots & \vdots & \ddots & \vdots  \\
%------------------
  \displaystyle \frac{\partial z_{iN}^l}{\partial w_{ij1}^l}  
& \displaystyle \frac{\partial z_{iN}^l}{\partial w_{ij2}^l}
& \displaystyle \dots
& \displaystyle \frac{\partial z_{iN}^l}{\partial w_{ijN}^l} \\
\end{bmatrix} \\
&=\mathbf{d}_{i}^l\displaystyle \frac{\partial \mathbf{z}_{i}^l}{\partial \mathbf{w}_{ij}^l} 
= \mathbf{d}_{i}^{l}[\mathbf{a}_{j}^{l-1}]_{\ast}^T\,.
\end{split}
\end{equation}
From Eq.~(\ref{eq:feedforward_equation}), we see that the last equality holds because
\renewcommand{\arraystretch}{0.85}
\setlength{\arraycolsep}{1.5pt}
\begin{equation}
\label{eq:partial_c_over_w_ij_l}
\begin{split}
\displaystyle \frac{\partial \mathbf{z}_{i}^l}{\partial \mathbf{w}_{ij}^l}
=
\begin{bmatrix}
a_{j1}^{l-1}      & a_{jN}^{l-1}       & \dots     & a_{j2}^{l-1}   \\ 
a_{j2}^{l-1}      & a_{j1}^{l-1}       & \dots     & a_{j3}^{l-1}   \\
\vdots            & \vdots             & \ddots    & \vdots         \\
a_{jN}^{l-1}      & a_{j(N-1)}^{l-1}   & \dots     & a_{j1}^{l-1}   \\
\end{bmatrix}
=[\mathbf{a}_{j}^{l-1}]_{\ast}^T~\textrm{,}
\end{split}
\end{equation}
where $[\mathbf{a}_{j}^{l-1}]_{\ast}$ is the circulant matrix associated to $\mathbf{a}_{j}^{l-1}$. 
%The gradient vector of $\mathbf{w}_{ij}^{l}$ is formulated as a vector-matrix multiplication, where the matrix is still a circulant matrix.

We now show how to calculate $\mathbf{d}_{i}^{l}$, which requires backpropagation from the upper layer $l+1$ and is therefore more complicated. According to the chain rule, we have
%After derivative $\frac{\partial C}{\partial \mathbf{w}_{ij}^l}$ is defined, we derive the local gradient vector $\mathbf{d_{i}^l}$  from layer $l+1$, which is defined by: 
\renewcommand{\arraystretch}{2.0}
\begin{equation}
\label{eq:delta_propagation}
%\begin{split}
\mathbf{d}_{i}^{l}
\triangleq
\frac{\partial C}{\partial \mathbf{z}_{i}^l} 
=
\sum_{j=1}^{J^{l+1}} \frac{\partial C}{\partial \mathbf{z}_j^{l+1}}
\frac{\partial \mathbf{z}_j^{l+1}}{\partial \mathbf{a}_i^l}
\frac{\partial \mathbf{a}_i^l}{\partial \mathbf{z}_i^l} \,,
%\end{split}
\end{equation}
which involves vector-matrix-matrix multiplications. From Eq. (\ref{eq:partial_C_over_partail_bi}), we see that the first term in the summation %$\frac{\partial C}{\partial \mathbf{z}_j^{l+1}}$ 
is equal to $\mathbf{d}_{j}^{l+1}$. From Eq.~(\ref{eq:feedforward_equation}), and with some algebra, it can be shown that the middle term %$\frac{\partial \mathbf{z}_j^{l+1}}{\partial \mathbf{a}_i^l}$ 
is equal to $[\mathbf{w}_{ji}^{l+1}]_{\ast} \in \mathbb{R}^{N \times N}$, and the last term %$\frac{\partial \mathbf{a}_i^l}{\partial \mathbf{z}_i^l}$ 
is equal to the following $N \times N$ diagonal matrix: 
\renewcommand{\arraystretch}{1.5}
\setlength{\arraycolsep}{2pt}
\begin{equation}
\begin{split}
\Phi({\mathbf{z}_{i}^{l}})
\triangleq
\frac{\partial \mathbf{a}_{i}^{l}}{\partial \mathbf{z}_{i}^l}
=
\text{diag} \left(\left[\dot{\phi}({z_{i1}^{l}}), ~\dot{\phi}({z_{i1}^{2}}), ~\dots, ~\dot{\phi}({z_{i1}^{N}})\right]\right) \,,
\end{split}
\end{equation} 
where $\dot{\phi}(\cdot)$ denotes the derivative of the activation function. From the above description, we see that the update of $\mathbf{d}_{i}^{l}$ depends on $\mathbf{d}_{j}^{l+1}$ and we can do this layer-by-layer in the backward direction starting from the last layer. 
%where it is evident that the local gradient vector $\mathbf{d}_{i}^{l}$ is backpropagated from $\mathbf{d}_{k}^{l+1}$ in the layer $l+1$, which can be inferred from the output layer. 
In the output layer $L$, the local gradient vector $\mathbf{d}^{L}$ can be computed by:
\begin{equation}
\label{eq:delta_calculation_at_layer_L}
\begin{split}
\mathbf{d}^{L} 
&
\triangleq 
\frac{\partial C}{\partial \mathbf{z}^{L}}
= \frac{\partial C}{\partial \mathbf{y}}\frac{\partial\mathbf{y}}{\partial \mathbf{z}^{L}}\\
&
= \frac{\partial C}{\partial \mathbf{y}}\frac{\partial\mathbf{a}^L}{\partial \mathbf{z}^{L}}
= \frac{\partial C}{\partial \mathbf{y}}\Phi({\mathbf{z}^{L}}) \,,
\end{split}
\end{equation}
where the last term involves computing the derivative of the cost function $\epsilon(\cdot,\cdot)$.


\begin{algorithm}[t]
\caption{Backpropagation algorithm for updating deep cyclic-group networks}
\begin{algorithmic}[1]
\REQUIRE Training inputs $\mathbf{x}$ and related targets $\hat{\mathbf{y}}$; activation \\
		\hspace{0.45cm} function $\phi$; cost function $\epsilon$
\ENSURE  Parameters $\mathbf{\Theta}$ of weights and biases 

\WHILE{not converged}
\FOR{\textbf{each} minibatch  }
  \STATE Perform feedforward process with Eq.~(\ref{eq:feedforward_equation}) to get $\mathbf{z}$ 
  %--\ref{eq:activation_L} to find input vector $\mathbf{z}$ and output vector $\mathbf{a}$ of each neuron in each hidden layer
  \STATE Perform backpropagation with Eqs.~(\ref{eq:delta_propagation}) and (\ref{eq:delta_calculation_at_layer_L}) to get local gradient vector $\mathbf{d}$ per neuron per layer
  \STATE Update biases $\mathbf{b}$ with Eq.~(\ref{eq:partial_C_over_partail_bi})
  \STATE Update weights $\mathbf{w}$ with Eq.~(\ref{eq:derivative_w_fun})
  %with the gradients from
\ENDFOR 
\ENDWHILE 
\end{algorithmic}
\label{alg:cccnn_alg}
\end{algorithm}

We summarize the aforementioned process in Algorithm~\ref{alg:cccnn_alg}. We note that the proposed learning algorithm can be applied to vector neurons with arbitrary dimensions $N$. If the complexity of CNN is $\mathcal{O}(C)$, then the complexity of DCGN is $\mathcal{O}(NC)$. In other words, the time complexity of DCGN is $N$ times the complexity of CNN. However, if we view $N$ as a constant, the complexity of DCGN and CNN becomes the same.





%---------------------------------------------------------------
\section{Weight initialization}
\label{sec:ch3_weight_initialization}
Proper weight initialization can reduce the risk of gradient vanishing or exploding when the architecture of neural networks is deep. 
% It also helps the model parameters convergence. [Thomas deleted because convergence is mathematically unproven]
Below, we propose a weight initialization method for DCGN by extending the method proposed by \cite{xavier10aistats} and \cite{he15cv} to DCGNs. 
 
Initially, the weights are all scalars in the conventional CNNs but now all vectors in vector neurons. Since we know that variance is the expectation of the square minus square of the expectation, the equation of weight variance becomes:
\begin{equation}
\begin{split}
\textrm{Var}(\mathbf{w})
&=\mathbb{E}\left[\mathbf{w}\mathbf{w}^T\right]-\left(\mathbb{E}\left[\mathbf{w}\right]\right)^2 \\
&=\mathbb{E}\left[\vert \mathbf{w} \vert^2\right]-\left(\mathbb{E}\left[\mathbf{w}\right]\right)^2\textrm{,}
\end{split}
\end{equation}
where $\mathbb{E}\left[\mathbf{w}\right]$ is $0$ since the distribution of $\mathbf{w}$ is symmetric around $0$. 
However, estimating $\mathbb{E}\left[\vert \mathbf{w} \vert^2\right]$ is not easy to accomplish, so the main idea is to make the output variance equal to the input variance for each neuron. In doing so, we need to estimate the variance of each kernel map $\mathbf{w}$, which can be calculated as~\cite{gaudet17arxiv,parcollet18arxiv}:
\begin{equation}
\label{eq:w_variance}
\textrm{Var}(\mathbf{w})
=
\mathbb{E}\left[\vert \mathbf{w} \vert^2\right] \,.
\end{equation}

However, such an estimation is ill-defined for DCGN, since $\mathbf{w}$ is a vector not a scalar. That said, if we further assume $w_1$, $w_2$, $\ldots$, $w_N\sim\mathcal{N}\left(0,\sigma^2\right)$ i.i.d., then $\mathbb{E}\left[\vert\mathbf{w}\vert^2\right]=\mathbb{E}\left[w_1^2+w_2^2+\ldots+w_N^2\right]=\mathbb{E}\left[w_1^2\right]+\mathbb{E}\left[w_2^2\right]+\ldots+\mathbb{E}\left[w_N^2\right]=N\sigma^2$. Thus we can rewrite Eq.~(\ref{eq:w_variance}) as:
\begin{equation}
\label{eq:variance_equation}
\begin{split}
\textrm{Var}(\mathbf{w})
&=
\mathbb{E}\left[\vert \mathbf{w} \vert^2\right] \\
&=
\int_{-\infty}^{\infty}x^2f\left(x\right)dx \\
&=N\sigma^2\textrm{,}
\end{split}
\end{equation}
where 
$f(\cdot)$ stands for $N$-dimensional distribution with $N$ degree-of-freedom and $\sigma$ is an estimated parameter of the variance of $\mathbf{w}$. Then, if we use the initialization method proposed by \cite{xavier10aistats}, we have:
\begin{equation}
\textrm{Var}(\mathbf{w})=\frac{2}{\left(n_{in}+n_{out}\right)}\textrm{,}    
\end{equation}
where $n_{in}$ and $n_{out}$ denote the number of input and output units of that neuron, respectively. Accordingly, we have:
\begin{equation}
\begin{split}
\sigma^2=\frac{2}{N\left(n_{in}+n_{out}\right)}\textrm{.}
\end{split}
\end{equation}
On the other hand, if we use the initialization method proposed by \cite{he15cv}, which is designed for networks that use rectified linear units (ReLUs)~\cite{hinton10icml} as the activation function $\phi(\cdot)$, we would have $\textrm{Var}(\mathbf{w})=2/n_{in}$, and
$\sigma^2=\frac{2}{N\,n_{in}}$.
Then the variance is written as:
\begin{equation}
\begin{split}
\sigma=\sqrt{\frac{2}{N\,n_{in}}}\textrm{.}
\end{split}
\end{equation}


%---------------------------------------------------------------
\section{Batch normalization}
\label{sec:ch3_batch_normalization}
Batch normalization~\cite{ioffe15arxiv} is an important technique to stabilize and speed up the training process. The idea is to keep the input of each layer having zero mean and unit variance. 
The original version is designed for scalar neurons.
While \cite{gaudet17arxiv} have extended the method for DQNs, we present below how to extend it to general vector neurons with arbitrary dimension $N$. To ensure equal variance in the vector parts, we calculate the variance not only between different training data instances but also between the $N$ dimensions, leading to the following covariance matrix 
$\mathbf{V}$:
\begin{equation}
\begin{split}
\mathbf{V}=\left[\text{Cov}(u_i, v_j)\right]_{~i=1,\dots,N\,,~j=1,\dots,N} \,,
\end{split}
\end{equation}
where $\text{Cov}(\cdot)$ is the covariance function and $u_i$, $v_j$ denotes entries of two input vectors $\mathbf{u}$ and $\mathbf{v}$ for that layer.
Whitening approach~\cite{trabelsi17arxiv}, scaling the data by the square root of their variances along each of two-dimensional vectors, is used for complex-valued neurons. It can be completed by multiplying the inverse square root of the $2 \times 2$ covariance matrix. However, it is not easy to calculate the inverse square root of a matrix which dimension is more than $2 \times 2$.
Following \cite{gaudet17arxiv}, Cholesky decomposition is firstly proposed to for covariance in the deep quaternion networks, which is applied to four-dimensional data. In our proposed work, we use Cholesky decomposition on the covariance matrix to learn the shift vector $\mathbf{\beta} \in \mathbb{R}^N$ and the transformation matrix $\mathbf{\Gamma} \in \mathbb{R}^{N \times N}$.
%In real-valued batch normalization, two parameters are used to learned, $\beta$ and $\gamma$ respectively. The shift vector in our work $\beta$ becomes a $N$-dimensional vector. 
If the input vector $\tilde{\mathbf{x}}$ has been normalized to mean $0$ and variance $1$, then batch normalization is done by
\begin{equation}
\text{BN}(\tilde{\mathbf{x}})=\mathbf{\Gamma}\tilde{\mathbf{x}}+\mathbf{\beta}  \,.
\end{equation}
We note that $\mathbf{\Gamma}$ is a symmetric matrix, and that its diagonal can be initialized to $1/\sqrt{N}$ to obtain modulus of 1 for the variance of the normalized value. The other terms of $\mathbf{\Gamma}$ and all entries of $\mathbf{\beta}$ are initialized to $0$.

%---------------------------------------------------------------
\section{Pooling layers}
\label{sec:ch3_pooling_layers}
As shown in Figure \ref{fig:dcgn_arch}, we can optionally use pooling layers after the convolutional layer. A pooling layer is another type of layers simplifying or summarizing the information from convolved feature maps. It is also called subsampling or downsampling which reduces the dimensionality of each map and produce the downsampled one but retains the important information. 
The pooling layer is also used if necessary. The pooling methods we employ are max pooling and average pooling, which are also commonly used in CNNs. What the difference is that we propose to calculate the magnitude of each mode-$3$ fiber in the feature maps, and preserve the one (i.e., a vector) with the largest magnitude for max pooling. For mean pooling, we calculate the average of each element among the mode-$3$ fibers to generate new three-way feature maps.


%---------------------------------------------------------------
\section{Applications to regression problems}
\label{sec:ch3_applications_to_regression_problems}
DCGN can be applied to regression problems where we need to learn a nonlinear mapping between tensor inputs and tensor outputs. The inputs and outputs have the same mode-$3$ dimension, but not necessarily the same mode-$1$ and mode-$2$ dimensions. 
To fit the dimensionality of the output tensor, we need to use only one kernel map for the last convolutional layer of DCGN. For the activation function $\phi(\cdot)$, we can use rectified linear unit (ReLU) or leaky ReLU for all the layers except for the last layer, where we may want to use the sigmoid or tanh function. We can use squared error for the cost function $\epsilon(\mathbf{A},\mathbf{B})=\|\mathbf{A}-\mathbf{B}\|_2^2$, where $\mathbf{A}$ and $\mathbf{B}$ denote the groundtruth target and the output tensors of DCGN, respectively. Such a network is shown in Figure \ref{fig:dcgn_arch}.


%---------------------------------------------------------------
\section{Applications to classification problems}
\label{sec:ch3_applications_to_classification_problems}
When DCGN is applied to classification problems, the input is a three-way tensor and the corresponding output is a vector for one-hot representation (i.e., the class labels). We can also use only one kernel map for the last convolutional layer, but this is not mandatory. Moreover, we would add fully-connected layers to the end of DCGN for learning the classifier. We usually use the sigmoid function as the activation function for the last layer of the fully-connected layers, so that we can use the cross entropy as the cost function: $\epsilon(\mathbf{a},\mathbf{b})=-\sum [a_i \ln b_i + (1-a_i) \ln(1-b_i)]$, where $\mathbf{a}$ and $\mathbf{b}$ are the groundtruth label vector and the predicted one,
%probability distributions 
respectively. In such a case, DCGN performs the function of feature learning.
% $\epsilon(\mathbf{a},\mathbf{b})=-\sum [\mathbf{a} \ln\mathbf{b} + (1-\mathbf{a}) \ln(1-\mathbf{b})]$.
% $\epsilon(\mathbf{a},\mathbf{b})=-\sum_{x} \mathbf{a}(x)\log \mathbf{b}(x)$