\chapter{Arbitrary Bilinear Product Neural Networks}
\label{chp:abpinn_method}

The idea of ABIPNN is to extend the data type of each neuron from scalars to vectors by replacing multiplications with arbitrary bilinear products. Such products can be used to learn the associations between vector elements in the input (see Fig.~\ref{fig:vnn_explain}). 
In the following, tensors are represented by bold uppercase calligraphic letters, matrices by bold uppercase letters, and vectors by bold lowercase letters. Matrix slices of tensors are represented as matrices and vector slices of matrices are represented as vectors. The notational conventions used in this paper are summarized in Table~\ref{tab:convention_table}. 


%---------------------------------------------------------------
\section{Generalizing N-D vector neurons with bilinear Products}
\label{sec:ch2_bilinear_product_operation}
To begin with, we identify a generalization of all the products used in previous vector-valued NN literature, including vector product~\cite{nitta93ijcnn}, quaternion multiplication~\cite{arena94iscas}, octonion multiplication~\cite{popa16icann}, and so on. We have been able to confirm that all of these products are bilinear (defined below). Motivated by this important observation, we will extend the vector-valued NN~\cite{nitta07ijcai} to use any kind of bilinear product between two $N$-dimensional vectors.


From the inspiration of Chan's work~\cite{chan16tsp}, we could assume an $N$-dimensional space with standard basis $\left\{\mathbf{e}_1,~\mathbf{e}_2,~\ldots,~\mathbf{e}_N\right\}$ as column vectors. Let $\mathbf{p}=\sum_{n=1}^N p_n\mathbf{e}_n=\left[p_1,~p_2,~\ldots,~p_N\right]^T$ and $\mathbf{q}=\sum_{n=1}^N q_n\mathbf{e}_n=\left[q_1,~q_2,~\ldots,~q_N\right]^T$ be two vectors in this space. A product $\bullet\colon \mathbb{R}^N\times \mathbb{R}^N\to \mathbb{R}^N$ is bilinear if and only if $\mathbf{p}\bullet \mathbf{q}$ is linear when we hold one of $\mathbf{p}$ or $\mathbf{q}$ fixed. If $\bullet$ is bilinear, we have:
\begin{equation}
\mathbf{p}\bullet \mathbf{q}=\mathbf{p}\bullet\left(\sum_{n=1}^N q_n\mathbf{e}_n\right)\textrm{,}
\end{equation}
where we have expanded the second term. Due to bilinearity,
\begin{equation}
\renewcommand{\arraystretch}{0.60}
\label{eq:biliner_product_second_term}
\begin{split}
\mathbf{p}\bullet \mathbf{q}&=\sum_{n=1}^N \left(\mathbf{p}\bullet \mathbf{e}_n\right)q_n\\
&=\left[\mathbf{p}\bullet 
\begin{bmatrix}
1 \\
0 \\
\vdots \\
0 
\end{bmatrix},
~\mathbf{p}\bullet 
\begin{bmatrix}
0 \\
1 \\
\vdots \\
0
\end{bmatrix},
~\ldots,~\mathbf{p}\bullet 
\begin{bmatrix}
0 \\
0 \\
\vdots \\
1 
\end{bmatrix}
\right]
\begin{bmatrix}
q_1 \\
q_2 \\
\vdots \\
q_N 
\end{bmatrix}\\
&=\left[\mathbf{p}\bullet \mathbf{e}_1,~\mathbf{p}\bullet \mathbf{e}_2,~\ldots,~\mathbf{p}\bullet \mathbf{e}_N\right]\mathbf{q}.
\end{split}
\end{equation}
If we let $\left[\mathbf{p}\right]_\bullet=\left[\mathbf{p}\bullet \mathbf{e}_1,~\mathbf{p}\bullet \mathbf{e}_2,~\ldots,~\mathbf{p}\bullet \mathbf{e}_N\right]$, then:
\begin{equation}
\label{eq:biliner_product_second_term_result}
\mathbf{p}\bullet \mathbf{q}=\left[\mathbf{p}\right]_\bullet \mathbf{q}\textrm{.}
\end{equation}
Similarly, if we expand the first term, we can write:
\begin{equation}
\label{eq:bilinear_product_first_term}
\begin{split}
\mathbf{p}\bullet\mathbf{q}
&=\left(\sum_{n=1}^{N}p_n\mathbf{e}_n\right)\bullet\mathbf{q} =\sum_{n=1}^{N}p_n\left(\mathbf{e}_{n}\bullet\mathbf{q}\right)\\
&=\left[\mathbf{e}_1\bullet\mathbf{q},~\mathbf{e}_2\bullet\mathbf{q},~\dots,~ \mathbf{e}_N\bullet\mathbf{q}\right]\mathbf{p}\textrm{.}
\end{split}
\end{equation}
Let $\left[\mathbf{q}\right]_\bullet^\dagger=\left[\mathbf{e}_1\bullet \mathbf{q},~\mathbf{e}_2\bullet \mathbf{q},~\ldots,~\mathbf{e}_N\bullet \mathbf{q}\right]$, we also have:
\begin{equation}
\label{eq:bilinear_product_first_term_result}
\mathbf{p}\bullet\mathbf{q}=\left[\mathbf{q}\right]_{\bullet}^{\dagger}\mathbf{p}.
\end{equation}
Eqs.~\ref{eq:biliner_product_second_term_result} and \ref{eq:bilinear_product_first_term_result} are the two possible matrix representations of the bilinear product $\bullet$. Here we call $\left[\mathbf{p}\right]_\bullet$ its matrix representation and $\left[\mathbf{q}\right]_\bullet^\dagger$ its transmuted representation (the term \emph{transmuted} is originally used for quaternions~\cite{ell07arxiv} but here we extend it for general bilinear products). In what follows, the feedforward and backpropagation processes will be described by using the above notations for bilinear products and their representations. For convenience, in the rest of this paper, if $\mathbf{p}$ or $\mathbf{q}$ are not column vectors, then $\mathbf{p}\bullet\mathbf{q}$ will implicitly reshape them so that the above equations make sense.

\newpage

%%-------------------------- Notation Table --------------------------%%
\begin{table}[h]
\renewcommand{\arraystretch}{1.2}
\caption{\textsc{Notational Conventions}}
\label{tab:convention_table}
% \vspace{-15pt}
\centering  
\begin{tabular}{|l l|}  
% \hline 
\hline
$\bm{\mathcal{X}}$                                    & Input tensor  \\
$\mathbf{X}_{m}$=$\bm{\mathcal{X}}_{:m:}$             & Lateral slice from input tensor $\bm{\mathcal{X}}$ \\
$\mathbf{x}$                                          & Vector from input tensor $\bm{\mathcal{X}}$  \\
$\bm{\mathcal{Y}}$                                    & Output tensor  \\
$\mathbf{Y}_{m}$=$\bm{\mathcal{Y}}_{:m:}$             & Lateral slice from output tensor $\bm{\mathcal{X}}$ \\
$\bm{\mathcal{W}}^{l}$                                & Weight tensor in layer $l$        \\ 
$\bm{\mathcal{W}}^{l+1}$                              & Weight tensor in layer $l+1$      \\  
$\mathbf{w}_{ij}^{l}=\bm{\mathcal{W}}_{ij:}^{l}$      & Vector from tensor $\bm{\mathcal{W}}^{l}$  \\
$\mathbf{w}_{ki}^{l+1}=\bm{\mathcal{W}}_{ki:}^{l+1}$  & Vector from tensor $\bm{\mathcal{W}}^{l+1}$  \\ 
$\bm{\mathcal{B}}_{}^{l}$                             & Bias tensor in layer $l$        \\ 
$\bm{\mathcal{B}}_{}^{l+1}$                           & Bias tensor in layer $l+1$      \\  
$\mathbf{b}_{i}^{l}=\bm{\mathcal{B}}_{im:}^{l}$       & Vector from tensor $\bm{\mathcal{B}}^{l}$  \\
$\mathbf{b}_{k}^{l+1}=\bm{\mathcal{B}}_{km:}^{l+1}$   & Vector from tensor $\bm{\mathcal{B}}^{l+1}$  \\ 
$\bm{\mathcal{Z}}^{L}$                                & Input tensor in the layer $L$      \\ 
$\bm{\mathcal{Z}}^{l}$                                & Input hidden tensor in layer $l$        \\ 
$\bm{\mathcal{Z}}^{l+1}$                              & Input hidden tensor in layer $l+1$        \\ 
$\mathbf{z}_{g}^{L}=\bm{\mathcal{Z}}_{gm:}^{L}$       & Vector from tensor $\bm{\mathcal{Z}}^{L}$  \\
$\mathbf{z}_{i}^{l}=\bm{\mathcal{Z}}_{im:}^{l}$       & Vector from tensor $\bm{\mathcal{Z}}^{l}$  \\
$\mathbf{z}_{k}^{l+1}=\bm{\mathcal{Z}}_{km:}^{l+1}$   & Vector from tensor $\bm{\mathcal{Z}}^{l+1}$  \\
$\bm{\mathcal{A}}^{l-1}$                              & Output hidden tensor in layer $l-1$      \\ 
$\bm{\mathcal{A}}^{l}$                                & Output hidden tensor in layer $l$        \\ 
$\mathbf{a}_{j}^{l-1}=\bm{\mathcal{A}}_{jm:}^{l-1}$   & Vector from tensor $\bm{\mathcal{A}}^{l-1}$  \\
$\mathbf{a}_{i}^{l}=\bm{\mathcal{A}}_{im:}^{l}$       & Vector from tensor $\bm{\mathcal{A}}^{l}$  \\ 
$\bm{\mathcal{D}}^{l+1}$                              & Local gradient tensor in layer $l+1$      \\ 
$\bm{\mathcal{D}}^{l}$                                & Local gradient tensor in layer $l$        \\ 
$\mathbf{d}_{j}^{l+1}=\bm{\mathcal{D}}_{im:}^{l+1}$   & Vector from tensor $\bm{\mathcal{D}}^{l+1}$  \\
$\mathbf{d}_{i}^{l}=\bm{\mathcal{D}}_{km:}^{l}$       & Vector from tensor $\bm{\mathcal{D}}^{l}$  \\ 
$\bullet\colon \mathbb{R}^N\times \mathbb{R}^N\to \mathbb{R}^N$ & Arbitrary bilinear product \\
$\left[\cdot\right]_\bullet$                          & Matrix representation of $\bullet$ (Eq.~\ref{eq:biliner_product_second_term_result})\\
$\left[\cdot\right]_\bullet^\dagger$                  & Transmuted representation of $\bullet$ (Eq.~\ref{eq:bilinear_product_first_term_result})\\
% \hline
\hline
\end{tabular}
\end{table}
%-------------------------- Notation Table --------------------------%
\newpage

\section{Feedforward process for ABIPNN}
\label{sec:ch2_feedforward}
As before, we assume the inputs, outputs, weights, and biases are $N$-dimensional vectors. Suppose we have an $L$-layer vector-valued NN. For 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}
\mathbf{z}_i^l=\sum_{j=1}^{R_{l-1}}\mathbf{w}_{ij}^{l}\bullet\mathbf{a}_j^{l-1}+\mathbf{b}_{i}^{l}~\textrm{,}
\end{equation}
where $\bullet$ denotes an arbitrary bilinear product, $\mathbf{w}_{ij}^l$ stands for the weight vector connecting a neuron $j$ ($1 \leq j \leq R_{l-1}$) in layer $l-1$ to a neuron $i$ ($1 \leq i \leq R_l$) in layer $l$, and
$\mathbf{w}_{ij}^l=\sum_{n=1}^N w_{ijn}^l\mathbf{e}_n=\left[w_{ij1}^l\textrm{,}~w_{ij2}^l\textrm{,}\dots\textrm{,}w_{ijN}^l\right]^T\in \mathbb{R}^N$. 
The output vector of neuron $j$ in layer $l-1$ is $\mathbf{a}_{j}^{l-1}=\sum_{n=1}^N a_{jn}^{l-1}\mathbf{e}_n=\left[a_{j1}^{l-1}\textrm{,}~a_{j2}^{l-1}\textrm{,}\dots\textrm{,}~a_{jN}^{l-1}\right]^T\in \mathbb{R}^N$,
the bias vector of neuron $i$ in layer $l$ is $\mathbf{b}_{i}^l=\sum_{n=1}^N b_{in}^{l}\mathbf{e}_n=\left[b_{i1}^l\textrm{,}~b_{i2}^l\textrm{,}\dots\textrm{,}~b_{iN}^l\right]^T \in \mathbb{R}^N$. 
As a result, $\mathbf{z}_i^l=\sum_{n=1}^N z_{in}^{l}\mathbf{e}_n=\left[z_{i1}^l\textrm{,}~z_{i2}^l\textrm{,}\dots\textrm{,}~z_{iN}^l\right]^T \in \mathbb{R}^N$. 
When $l$ equals $1$, the vector $\mathbf{z}$ is simply the input vector $\mathbf{x}$. After $\mathbf{z}_i^l$ is calculated, the output vector $\mathbf{a}_{i}^{l}$ of the hidden neuron $i$ is as follows:
\renewcommand{\arraystretch}{0.75}
\begin{equation}
\label{eq:feedforeard_at_layer_l}
\mathbf{a}_i^l=
\begin{bmatrix}
a_{i1}^l \\
a_{i2}^l \\
\vdots \\
a_{iN}^l
\end{bmatrix}
=\phi(\mathbf{z}_i^l)=
\begin{bmatrix}
\phi(z_{i1}^l) \\
\phi(z_{i2}^l) \\
\vdots           \\
\phi(z_{iN}^l)
\end{bmatrix} \textrm{,}
\end{equation}
where $\phi$ could be any differentiable activation function. Then the output vector of neuron $g$ $\left(1 \leq g \leq R_{L}\right)$ in output layer $L$ is defined as:
\renewcommand{\arraystretch}{0.75}
\begin{equation}
\label{eq:feedforeard_at_layer_L}
\mathbf{y}_g^L=
\begin{bmatrix}
y_{g1}^{L} \\
y_{g2}^{L} \\
\vdots \\
y_{gN}^{L}
\end{bmatrix}
=\phi(\mathbf{z}_g^L)=
\begin{bmatrix}
\phi(z_{g1}^{L}) \\
\phi(z_{g2}^{L}) \\
\vdots           \\
\phi(z_{gN}^{L})
\end{bmatrix} \textrm{,}
\end{equation}
where $\mathbf{z}_g^L$ is the input vector and $\mathbf{y}_g^L$ is the output vector of a neuron $g$. The objective of the training process is to estimate the parameters that minimize the cost function, defined as:  
\begin{equation}
\label{eq:cost_function}
C(\bm{\Theta})=\sum_{m=1}^M loss(\mathbf{Y}_m, f(\mathbf{X}_m;\bm{\Theta}))~\textrm{,}
\end{equation}
where $M$ stands for the number of training data patterns, $\mathbf{\Theta}$ is the set of all training parameters (weights and biases), and $\mathbf{Y}_m$ is the training label related to the input $\bm{X}_m$. The output $f(\mathbf{X}_m;\mathbf{\Theta})$ is the predicted version of $\mathbf{Y}_m$, made up of $\mathbf{y}_g^L$ mentioned in Eq.~\ref{eq:cost_function}. The $loss(\mathbf{Y}_m, f(\mathbf{X}_m;\bm{\Theta}))$
function measures the difference between the predicted results and the training labels. 


\newpage
%---------------------------------------------------------------
\section{Backpropagation algorithm for ABIPNN}
\label{sec:ch2_backpropagate}
The bilinear product is also employed in the process of backpropagation learning. 
Before we present the error backpropagation process (Algorithm~\ref{alg:ABIPalg}), we need to first derive:
\begin{enumerate}
\item The gradients of the biases $\mathbf{b}_{i}^{l}$.
\item The gradients of the weights $\mathbf{w}_{ij}^{l}$.
\item The backpropagation of local gradients $\mathbf{d}_{k}^{l+1}\rightarrow\mathbf{d}_{i}^{l}$.
\end{enumerate}
Following Eq.~\ref{eq:biliner_product_second_term_result}, the bilinear product can be viewed as matrix-vector multiplication. Hence, Eq.~\ref{eq:feedforward} can be rewritten as
\renewcommand{\arraystretch}{1}
\begin{equation}
\begin{split}
\label{eq:feedforward_extend_part1}
\mathbf{z}_i^l
&=\sum_{j=1}^{R_{l-1}}\mathbf{w}_{ij}^{l}\bullet\mathbf{a}_{j}^{l-1}+\mathbf{b}_{i}^{l}  =\sum_{j=1}^{R_{l-1}}\left[\mathbf{w}_{ij}^{l}\right]_{\bullet}\mathbf{a}_{j}^{l-1}+\mathbf{b}_{i}^{l}\\
&=\sum_{j=1}^{R_{l-1}}\left[
\mathbf{w}_{ij}^l\bullet\mathbf{e}_1\textrm{,}~
\mathbf{w}_{ij}^l\bullet\mathbf{e}_2\textrm{,}~
\dots \textrm{,}
\mathbf{w}_{ij}^l\bullet\mathbf{e}_N\right]\mathbf{a}_{j}^{l-1}+\mathbf{b}_{i}^{l}\\
&=\sum_{j=1}^{R_{l-1}}\left(
a_{j1}^{l-1}\mathbf{w}_{ij}^l\bullet\mathbf{e}_1 + 
\dots +
a_{jN}^{l-1}\mathbf{w}_{ij}^l\bullet\mathbf{e}_N\right) +\mathbf{b}_{i}^{l}\textrm{.}\\
\end{split}
\end{equation}
Here, $\mathbf{w}_{ij}^{l}$ can be formulated as the summation of scalar-vector multiplications:
\begin{equation}
\begin{split}
\label{eq:feedforward_extend_part2}
\mathbf{z}_i^l=
&\sum_{j=1}^{R_{l-1}}\left(
a_{j1}^{l-1}
\left(
w_{ij1}^l \mathbf{e}_1+
%w_{ij2}^l \mathbf{e}_2+
\dots+w_{ijN}^l \mathbf{e}_N \right)\right. \bullet\mathbf{e}_1  \\
&\quad\thinspace+
a_{j2}^{l-1}
\left(
w_{ij1}^l \mathbf{e}_1+
%w_{ij2}^l \mathbf{e}_2+
\dots+w_{ijN}^l \mathbf{e}_N \right) \bullet\mathbf{e}_2 + \dots\\
&\quad\thinspace+
a_{jN}^{l-1}\left.
\left(
w_{ij1}^l \mathbf{e}_1+
%w_{ij2}^l \mathbf{e}_2+
\dots+w_{ijN}^l \mathbf{e}_N \right) \bullet\mathbf{e}_N\right) \thinspace+\mathbf{b}_i^l\textrm{.}
% ~\textrm{,}
\end{split}
\end{equation}
To estimate the first-order partial derivatives of a vector-valued function, we apply the concept of Jacobian matrix, which is a matrix containing all the partial derivatives. For the vector-valued function $f:\mathbb{R}^N\to\mathbb{R}^M$, the Jacobian matrix is defined as $\left[\frac{\partial f}{\partial \bm{x}}\right]_{ij}= \frac{\partial}{\partial x_{j}}f(\bm{x})_i\in\mathbb{R}^{N\times M}$.
The partial derivative of C with respect to $\mathbf{b}_{i}^{l}$ is as follows:
\renewcommand{\arraystretch}{1.1}
\begin{equation}
\label{eq:derivative_b_fun}
\begin{split}
\frac{\partial C}{\partial \mathbf{b}_{i}^l}
&=
\begin{bmatrix}
\displaystyle
\frac{\partial C}{\partial b_{i1}^l} \textrm{,}~ 
\frac{\partial C}{\partial b_{i2}^l} \textrm{,}~\dots\textrm{,}~
\frac{\partial C}{\partial b_{iN}^l}
\end{bmatrix} \\
&=
\begin{bmatrix}
\displaystyle \frac{\partial C}{\partial z_{i1}^l} \frac{\partial z_{i1}^l}{\partial b_{i1}^l}+  
\displaystyle \frac{\partial C}{\partial z_{i2}^l} \frac{\partial z_{i2}^l}{\partial b_{i1}^l}+ 
\displaystyle \dots+
\displaystyle \frac{\partial C}{\partial z_{iN}^l} \frac{\partial z_{iN}^l}{\partial b_{i1}^l} \\ 
%------------------
\displaystyle \frac{\partial C}{\partial z_{i1}^l} \frac{\partial z_{i1}^l}{\partial b_{i2}^l}+  
\displaystyle \frac{\partial C}{\partial z_{i2}^l} \frac{\partial z_{i2}^l}{\partial b_{i2}^l}+ 
\displaystyle \dots+
\displaystyle \frac{\partial C}{\partial z_{iN}^l} \frac{\partial z_{iN}^l}{\partial b_{i2}^l} \\
% ------------------
\displaystyle \frac{\partial C}{\partial z_{i1}^l} \frac{\partial z_{i1}^l}{\partial b_{i3}^l}+  
\displaystyle \frac{\partial C}{\partial z_{i2}^l} \frac{\partial z_{i2}^l}{\partial b_{i3}^l}+
\displaystyle \dots+
\displaystyle \frac{\partial C}{\partial z_{iN}^l} \frac{\partial z_{iN}^l}{\partial b_{i3}^l} \\
% ------------------
\displaystyle \vdots   \\
%------------------
\displaystyle \frac{\partial C}{\partial z_{i1}^l} \frac{\partial z_{i1}^l}{\partial b_{iN}^l}+  
\displaystyle \frac{\partial C}{\partial z_{i2}^l} \frac{\partial z_{i2}^l}{\partial b_{iN}^l}+
\displaystyle \dots+
\displaystyle \frac{\partial C}{\partial z_{iN}^l} \frac{\partial z_{iN}^l}{\partial b_{iN}^l} \\
\end{bmatrix}^T \textrm{,}
\end{split}
\end{equation}
where $\frac{\partial C}{\partial \mathbf{b}_{i}^l}$ stands for the gradient vector of $\mathbf{b}_{i}^l$ belonging to neuron $i$ in layer $l$. The terms $\frac{\partial C}{\partial z_{in}^l}$ ($1 \leq n \leq N$) are extracted into an $N$-dimensional local gradient vector $\frac{\partial C}{\partial \mathbf{z}_{i}^l}$ and we call it $\mathbf{d}_{i}^{l}$ for convenience:
\begin{equation}
\label{eq:delta_l}
\renewcommand{\arraystretch}{1}
\begin{split}
\frac{\partial C}{\partial \mathbf{z}_{i}^l}
&=\begin{bmatrix}
\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}
\end{bmatrix}\\
&\triangleq\begin{bmatrix}
d_{i1}^l\textrm{,}~
d_{i2}^l\textrm{,}~
\dots\textrm{,}~
d_{iN}^l
\end{bmatrix}
=\mathbf{d}_{i}^{l} \textrm{.}
\end{split}
\end{equation}
The terms $\left\{\frac{\partial z_{i1}^{l}}{\partial b_{in}^{l}}\textrm{,}~\frac{\partial z_{i2}^{l}}{\partial b_{in}^{l}}\textrm{,}~\dots\textrm{,}~\frac{\partial z_{iN}^{l}}{\partial b_{in}^{l}}\right\}$ ($1 \leq n \leq N$) can be obtained by differentiating Eq.~\ref{eq:feedforward_extend_part2}:
\renewcommand{\arraystretch}{1.0}
\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_{in}^l} \textrm{,}~ 
\frac{\partial z_{i2}^l}{\partial b_{in}^l} \textrm{,}~\dots\textrm{,}~
\frac{\partial z_{iN}^l}{\partial b_{in}^l}
\end{bmatrix}^T \\
&=
\begin{bmatrix}
\displaystyle
\frac{\partial b_{i1}^l}{\partial b_{in}^l} \textrm{,}~ 
\frac{\partial b_{i2}^l}{\partial b_{in}^l} \textrm{,}~\dots\textrm{,}~
\frac{\partial b_{iN}^l}{\partial b_{in}^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}
By combining Eqs.~\ref{eq:delta_l} and \ref{eq:partial_z_il_b_in}, the derivative $\frac{\partial C}{\partial b_{in}^l}$ can be formulated as a dot product:
\renewcommand{\arraystretch}{0.7}
\begin{equation}
\label{eq:partial_C_bin_result}
\begin{split}
\frac{\partial C}{\partial b_{in}^l}
&=
\mathbf{d}_i^l\left[\frac{\partial \mathbf{z}_i^l}{\partial b_{in}^l}\right] =
\underbrace{
\begin{bmatrix}
d_{i1}^l\textrm{,}~
d_{i2}^l\textrm{,}~ 
\dots\textrm{,}~ 
d_{iN}^l
\end{bmatrix} 
}_{1 \times N} 
\underbrace{
\begin{bmatrix}
0 \\
\vdots \\
1 \\
\vdots \\ 
0 
\end{bmatrix}
}_{N \times 1} 
=
d_{in}^l \textrm{.}
\end{split}
\end{equation}
Then, each entry of $\frac{\partial C}{\partial \mathbf{b}_{i}^l}$ can be estimated by the same procedure from Eqs.~\ref{eq:partial_z_il_b_in} and \ref{eq:partial_C_bin_result}. After each entry is estimated, Eq.~\ref{eq:derivative_b_fun} can be derived as follows:
\renewcommand{\arraystretch}{0.8}
\begin{equation}
\label{eq:partial_C_over_partail_bi}
\begin{split}
\frac{\partial C}{\partial \mathbf{b}_{i}^l}
&=
\begin{bmatrix}
\displaystyle
\frac{\partial C}{\partial b_{i1}^l} \textrm{,}~ 
%\frac{\partial C}{\partial b_{i2}^l} \textrm{,}~
\dots\textrm{,} ~
\frac{\partial C}{\partial b_{iN}^l}
\end{bmatrix} =
\begin{bmatrix}
d_{i1}^{l} \textrm{,}~
%d_{i2}^{l} \textrm{,}~
\dots \textrm{,}~
d_{iN}^{l}
\end{bmatrix}
= \mathbf{d}_{i}^{l}  \textrm{.}
\end{split}
\end{equation}
As a result, it is evident that the derivative of $\mathbf{b}_{i}^l$ is completely equal to local gradient vector $\mathbf{d}_{i}^{l}$. Next, we derive the partial derivative of $C$ with respect to $\mathbf{w}_{ij}^l$ as follows: 
\renewcommand{\arraystretch}{1.1}
\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 \frac{\partial C}{\partial z_{i1}^l} \frac{\partial z_{i1}^l}{\partial w_{ij3}^l}+  
\displaystyle \frac{\partial C}{\partial z_{i2}^l} \frac{\partial z_{i2}^l}{\partial w_{ij3}^l}+
\displaystyle \dots+
\displaystyle \frac{\partial C}{\partial z_{iN}^l} \frac{\partial z_{iN}^l}{\partial w_{ij3}^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 
\end{split}
\end{equation}
where $\frac{\partial C}{\partial \mathbf{w}_{ij}^l}$ stands for the gradient vector of $\mathbf{w}_{ij}^l$ connecting neuron $j$ to neuron $i$ in layer $l$ and it consists of $N$ elements. 
From Eq.~\ref{eq:delta_l}, here the terms $\left\{\frac{\partial C}{\partial z_{i1}^{l}}\textrm{,}~\frac{\partial C}{\partial z_{i2}^{l}}\textrm{,}~\dots\textrm{,}~\frac{\partial C}{\partial z_{iN}^{l}}\right\}$
are extracted as the local gradient vector $\mathbf{d}_{i}^{l}$.
Likewise, we can obtain the terms 
$\left\{\frac{\partial z_{i1}^l}{\partial w_{ijn}^{l}}\textrm{,}~\frac{\partial z_{i2}^l}{\partial w_{ijn}^{l}}\textrm{,}~\dots\textrm{,}~\frac{\partial z_{iN}^l}{\partial w_{ijn}^{l}}\right\}$ ($1 \leq n \leq N$)
by differentiating Eq.~\ref{eq:feedforward_extend_part2}:
\renewcommand{\arraystretch}{1.0}
\begin{equation}
\label{eq:partial_z_il_w_ijn}
\begin{split}
&\frac{\partial \mathbf{z}_{i}^l}{\partial w_{ijn}^l}
=
\begin{bmatrix}
\displaystyle
\frac{\partial z_{i1}^l}{\partial w_{ijn}^l} \textrm{,}~ 
\frac{\partial z_{i2}^l}{\partial w_{ijn}^l} \textrm{,}~\dots\textrm{,}~
\frac{\partial z_{iN}^l}{\partial w_{ijn}^l}
\end{bmatrix}^T \\
&=
a_{j1}^{l-1}\mathbf{e}_n\bullet\mathbf{e}_1+
a_{j2}^{l-1}\mathbf{e}_n\bullet\mathbf{e}_2+\dots+
a_{jN}^{l-1}\mathbf{e}_n\bullet\mathbf{e}_N\\
&=
\begin{bmatrix}
\mathbf{e}_n\bullet\mathbf{e}_1 \textrm{,}~
\mathbf{e}_n\bullet\mathbf{e}_2 \textrm{,} \dots \textrm{,}~
\mathbf{e}_n\bullet\mathbf{e}_N 
\end{bmatrix}
\begin{bmatrix}
a_{j1}^{l-1},
\dots 
a_{jN}^{l-1}
\end{bmatrix}^T \\
% \begin{bmatrix}
% a_{j1}^{l-1}\\
% a_{j2}^{l-1}\\
% \vdots \\
% a_{jN}^{l-1}\\
% \end{bmatrix} \\
&=
\begin{bmatrix}
\mathbf{e}_n\bullet\mathbf{e}_1 \textrm{,}~
\mathbf{e}_n\bullet\mathbf{e}_2 \textrm{,} \dots \textrm{,}~
\mathbf{e}_n\bullet\mathbf{e}_N 
\end{bmatrix}
\mathbf{a}_{j}^{l-1} \textrm{,}
\end{split}
\end{equation}
which is formulated as a matrix-vector multiplication. Each column in the matrix is the bilinear product of two standard bases. Different bilinear products contribute to different results. Then, the derivative $\frac{\partial C}{\partial w_{ijn}^l}$ can be formulated as a dot product by combining Eqs.~\ref{eq:delta_l} and \ref{eq:partial_z_il_w_ijn}:
\renewcommand{\arraystretch}{1.0}
\begin{equation}
\label{eq:partial_C_over_wijn_result}
\begin{split}
\frac{\partial C}{\partial w_{ijn}^l}
&=
\mathbf{d}_i^l\left[\frac{\partial \mathbf{z}_i^l}{\partial w_{ijn}^l}\right]\\
&=
\underbrace{
\begin{bmatrix}
d_{i1}^l\textrm{,}~
d_{i2}^l\textrm{,}~ 
\dots\textrm{,}~ 
d_{iN}^l
\end{bmatrix} 
}_{1 \times N}
\\
&\quad\quad\quad\quad\quad
\left[
\underbrace{
\begin{bmatrix}
\mathbf{e}_n\bullet\mathbf{e}_1 \textrm{,}~
\mathbf{e}_n\bullet\mathbf{e}_2 \textrm{,}~
\dots \textrm{,}~
\mathbf{e}_n\bullet\mathbf{e}_N 
\end{bmatrix}
}_{N \times N}
\underbrace{
\mathbf{a}_{j}^{l-1} 
}_{N \times 1} 
\right]\\
&=
\mathbf{d}_i^l\left[\sum_{h=1}^{N}a_{jh}^{l-1}\mathbf{e}_n\bullet\mathbf{e}_h\right]\textrm{.}
\end{split}
\end{equation}
Each entry of $\frac{\partial C}{\partial \mathbf{w}_{ij}^l}$ can be estimated by the same procedure from Eqs.~\ref{eq:partial_z_il_w_ijn} and \ref{eq:partial_C_over_wijn_result}. After each entry is estimated, Eq.~\ref{eq:derivative_w_fun} can be rewritten as:
\renewcommand{\arraystretch}{1.0}
\begin{equation}
\label{eq:partial_C_over_partail_wij}
\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
\mathbf{d}_i^l\left[\sum_{h=1}^{N}a_{jh}^{l-1}\mathbf{e}_1\bullet\mathbf{e}_h\right] %\\
%\mathbf{d}_i^l\left[\sum_{h=1}^{N}a_{jh}^{l-1}\mathbf{e}_2\bullet\mathbf{e}_h\right] %\\
,\dots, %\\
\mathbf{d}_i^l\left[\sum_{h=1}^{N}a_{jh}^{l-1}\mathbf{e}_N\bullet\mathbf{e}_h\right]
\end{bmatrix}\\%^T\\ 
&=
\mathbf{d}_{i}^{l}
\left[
\mathbf{e}_1\bullet\left[\sum_{h=1}^{N}a_{jh}^{l-1}\mathbf{e}_h\right]\textrm{,}~
%\mathbf{e}_2\bullet\left[\sum_{h=1}^{N}a_{jh}^{l-1}\mathbf{e}_h\right]\textrm{,} 
\right.
%\\
%&\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad
\dots\textrm{,}
\left.\mathbf{e}_N\bullet\left[\sum_{h=1}^{N}a_{jh}^{l-1}\mathbf{e}_h\right]\right]\textrm{,}
\end{split}
\end{equation}
and referring to Eq.~\ref{eq:bilinear_product_first_term}, we can rewrite Eq.~\ref{eq:partial_C_over_partail_wij} into:
\begin{equation}
\label{eq:partial_C_over_partial_wij_result}
\begin{split}
\frac{\partial C}{\partial \mathbf{w}_{ij}^l}
&=
\mathbf{d}_{i}^{l}
\left[ 
\mathbf{e}_{1}\bullet\mathbf{a}_{j}^{l-1}\textrm{,}~
\mathbf{e}_{2}\bullet\mathbf{a}_{j}^{l-1}\textrm{,}~
\dots\textrm{,}~
\mathbf{e}_{N}\bullet\mathbf{a}_{j}^{l-1}
\right] \\
&=
\mathbf{d}_{i}^{l}\left[\mathbf{a}_{j}^{l-1}\right]_{\bullet}^{\dagger}\textrm{,}
\end{split}
\end{equation}
which is a vector-matrix multiplication with the transmuted representation. After the derivative $\frac{\partial C}{\partial \mathbf{w}_{ij}^l}$ is calculated, we derive the local gradient vector $\mathbf{d}_{i}^l$ from layer $l+1$ by: 
\renewcommand{\arraystretch}{1}
\begin{equation}
\label{eq:delta_propagation}
\begin{split}
\mathbf{d}_{i}^{l}
=
\frac{\partial C}{\partial \mathbf{z}_{i}^l} 
&=
\sum_{k=1}^{R_{l+1}} \frac{\partial C}{\partial \mathbf{z}_{k}^{l+1}}
\frac{\partial \mathbf{z}_{k}^{l+1}}{\partial \mathbf{a}_{i}^{l}}
\frac{\partial \mathbf{a}_{i}^{l}}{\partial \mathbf{z}_{i}^l} \\
&=
\sum_{k=1}^{R_{l+1}} \frac{\partial C}{\partial \mathbf{z}_{k}^{l+1}} 
\frac{\partial \mathbf{z}_{k}^{l+1}}{\partial \mathbf{a}_{i}^{l}}
\dot{\phi}({\mathbf{z}_{i}^{l}}) \textrm{,}
\end{split}
\end{equation}
which is a vector-matrix-matrix multiplication, where $\mathbf{z}_{k}^{l+1}$ stands for an $N$-dimensional input vector of neuron $k$ in layer $l+1$ and $\mathbf{a}_{i}^{l}$ stands for an $N$-dimensional output vector of neuron $i$ in layer $l$. The vector $\frac{\partial C}{\partial \mathbf{z}_{k}^{l+1}}$ is regarded as an $N$-dimensional local gradient vector, defined as $\mathbf{d}_{k}^{l+1}$ in layer $l+1$. Referring to Eqs. \ref{eq:derivative_b_fun} to \ref{eq:partial_C_over_partail_bi}, we can see that $\mathbf{d}_{k}^{l+1}$ can be regarded as the derivative of $\mathbf{b}_{k}^{l+1}$. 
The matrix $\dot{\phi}({\mathbf{z}_{i}^{l}})$ represents the derivative of the activation function, which is an $N \times N$ diagonal matrix:
\renewcommand{\arraystretch}{1.0}
\begin{equation}
\begin{split}
\dot{\phi}({\mathbf{z}_{i}^{l}})
=
\frac{\partial \mathbf{a}_{i}^{l}}{\partial \mathbf{z}_{i}^l} 
=
\begin{bmatrix}
\displaystyle \frac{\partial a_{i1}^{l}}{\partial z_{i1}^{l}}  & 0          & \dots	     & 0            \\ 
0          & \displaystyle \frac{\partial a_{i2}^{l}}{\partial z_{i2}^{l}}  & \dots        & 0          \\
\vdots     & \vdots     & \ddots       & \vdots                                                         \\
0          & 0          & \dots        & \displaystyle \frac{\partial a_{iN}^{l}}{\partial z_{iN}^{l}}  \\
\end{bmatrix} 
=
\begin{bmatrix}
\dot{\phi}({z_{i1}^{l}})  & 0                         & \dots	           & 0                         \\ 
0                         & \dot{\phi}({z_{i2}^{l}})  & \dots              & 0                         \\
\vdots                    & \vdots                    & \ddots             & \vdots                    \\
0                    	  & 0                    	  & \dots              & \dot{\phi}({z_{iN}^{l}})  \\
\end{bmatrix} \textrm{,}
\end{split}
\end{equation} 
where $\dot{\phi}$ can be an arbitrary differentiable activation function.
Next, the matrix $\frac{\partial \mathbf{z}_{k}^{l+1}}{\partial \mathbf{a}_{i}^{l}}$ is also a $N\times N$ square matrix:
\renewcommand{\arraystretch}{1.1}
\begin{equation}
\label{eq:partial_zk_over_partial_ai}
\begin{split}
\frac{\partial \mathbf{z}_{k}^{l+1}}{\partial \mathbf{a}_{i}^{l}} 
&=
\begin{bmatrix}
& \displaystyle \frac{\partial z_{k1}^{l+1}}{\partial a_{i1}^{l}}   
& \displaystyle \frac{\partial z_{k1}^{l+1}}{\partial a_{i2}^{l}} 
& \displaystyle \dots
& \displaystyle \frac{\partial z_{k1}^{l+1}}{\partial a_{iN}^{l}}  \\ 
%------------------
& \displaystyle \frac{\partial z_{k2}^{l+1}}{\partial a_{i1}^{l}} 
& \displaystyle \frac{\partial z_{k2}^{l+1}}{\partial a_{i2}^{l}}  
& \displaystyle \dots
& \displaystyle \frac{\partial z_{k2}^{l+1}}{\partial a_{iN}^{l}}  \\
%------------------
& \vdots   
& \vdots 
& \ddots
& \vdots  \\
%------------------
& \displaystyle \frac{\partial z_{kN}^{l+1}}{\partial a_{i1}^{l}}  
& \displaystyle \frac{\partial z_{kN}^{l+1}}{\partial a_{i2}^{l}} 
& \displaystyle \dots
& \displaystyle \frac{\partial z_{kN}^{l+1}}{\partial a_{iN}^{l}}  \\
\end{bmatrix} \textrm{.}
\end{split}
\end{equation}
We calculate the above derivatives columnwise. The calculations are based on the feedforward process of $\mathbf{z}_k^{l+1}$:
\renewcommand{\arraystretch}{1.2}
\begin{equation}
\begin{split}
\label{eq:feedforward_extend_l+1}
\mathbf{z}_k^{l+1}
&=\sum_{i=1}^{R_l}\mathbf{w}_{ki}^{l+1}\bullet\mathbf{a}_{i}^{l}+\mathbf{b}_{k}^{l+1}
=\sum_{i=1}^{R_l}\mathbf{w}_{ki}^{l+1}\bullet\left(\sum_{n=1}^{N}a_{in}^{l}\mathbf{e}_{n}\right)+\mathbf{b}_{k}^{l+1}\\
&=\sum_{i=1}^{R_l}\mathbf{w}_{ki}^{l+1}\bullet\left[a_{i1}^{l}\mathbf{e}_{1}+a_{i2}^{l}\mathbf{e}_{2}+\dots +a_{iN}^{l}\mathbf{e}_{N}\right]+\mathbf{b}_{k}^{l+1} \textrm{.}\\
\end{split}
\end{equation}
Then, the derivative of the $n$-th column of Eq.~\ref{eq:partial_zk_over_partial_ai} is derived by differentiating Eq.~\ref{eq:feedforward_extend_l+1}:
\renewcommand{\arraystretch}{1.3}
\begin{equation}
\begin{split}
\frac{\partial \mathbf{z}_k^{l+1}}{\partial a_{in}^l}
&=
\begin{bmatrix}
\displaystyle
\frac{\partial z_{k1}^{l+1}}{\partial a_{in}^l} \textrm{,} 
\frac{\partial z_{k2}^{l+1}}{\partial a_{in}^l} \textrm{,}\dots\textrm{,} 
\frac{\partial z_{kN}^{l+1}}{\partial a_{in}^l}
\end{bmatrix}^T \\
&=\mathbf{w}_{ki}^{l+1}\bullet\mathbf{e}_n \\
&=
\left(
w_{ki1}^{l+1}\mathbf{e}_1+w_{ki2}^{l+1}\mathbf{e}_2+\dots+w_{kiN}^{l+1}\mathbf{e}_N
\right)\bullet\mathbf{e}_n\\
&=
\sum_{h=1}^{N}w_{kih}^{l+1}\mathbf{e}_h\bullet\mathbf{e}_n \\
&=
\underbrace{
\left[
\mathbf{e}_1\bullet\mathbf{e}_n\textrm{,}~
\mathbf{e}_2\bullet\mathbf{e}_n\textrm{,}~
\dots~\textrm{,}
\mathbf{e}_N\bullet\mathbf{e}_n
\right]}_{N \times N}
\underbrace{
\mathbf{w}_{ki}^{l+1} 
}_{N \times 1} \textrm{,}
\end{split}
\end{equation}
resulting in a column vector. Since each column of the matrix $\frac{\partial \mathbf{z}_{k}^{l+1}}{\partial \mathbf{a}_{i}^l}$ can be estimated by the same procedure mentioned above, the derivatives of the matrix can be derived as:
\begin{equation}
\label{eq:partial_zk_over_partial_ai_result}
\begin{split}
\frac{\partial \mathbf{z}_{k}^{l+1}}{\partial \mathbf{a}_{i}^l}
&=
\left[
\left[\sum_{h=1}^{N}w_{kih}^{l+1}\mathbf{e}_h\bullet\mathbf{e}_1\right]\textrm{,}~
\left[\sum_{h=1}^{N}w_{kih}^{l+1}\mathbf{e}_h\bullet\mathbf{e}_2\right]\textrm{,} 
\right.\\
&\quad\quad\quad\quad\quad\quad\quad\quad\thinspace
\dots
\textrm{,}
\left.
\left[\sum_{h=1}^{N}w_{kih}^{l+1}\mathbf{e}_h\bullet\mathbf{e}_N\right]\right] \\
&=
\left[
\mathbf{w}_{ki}^{l+1}\bullet\mathbf{e}_1 \textrm{,}~
\mathbf{w}_{ki}^{l+1}\bullet\mathbf{e}_2 \textrm{,}~
\dots \textrm{,}~
\mathbf{w}_{ki}^{l+1}\bullet\mathbf{e}_N
\right] \\
&=
\left[\mathbf{w}_{ki}^{l+1}\right]_{\bullet}\textrm{,}
\end{split}
\end{equation}
in which the matrix is made up of $N$ vectors. After estimating the matrix $\frac{\partial \mathbf{z}_{k}^{l+1}}{\partial \mathbf{a}_{i}^l}$, Eq.~\ref{eq:delta_propagation} can be rewritten as follows:
\begin{equation}
\label{eq:delta_propagation_new}
\begin{split}
\mathbf{d}_{i}^{l}
&=
\sum_{k=1}^{R_{l+1}}\mathbf{d}_{k}^{l+1}\left[\mathbf{w}_{ki}^{l+1}\right]_{\bullet}
\dot{\phi}({\mathbf{z}_{i}^{l}})\textrm{,}
\end{split}
\end{equation}
which becomes a vector-matrix-vector multiplication. It is evident that local gradient vector $\mathbf{d}_{i}^{l}$ in layer $l$ is backpropagated from $\mathbf{d}_{k}^{l+1}$ in layer $l+1$. From this we see that the local gradient can be inferred from the output layer. In the output layer $L$, the local gradient vector $\mathbf{d}_{g}^{L}$ is defined as:
\renewcommand{\arraystretch}{1.1}
\begin{equation}
\label{eq:delta_calculation_at_layer_L}
\begin{split}
\frac{\partial C}{\partial \mathbf{z}_{g}^{L}}
&=\begin{bmatrix}
\displaystyle
\frac{\partial C}{\partial z_{g1}^{L}} \textrm{,}~
\frac{\partial C}{\partial z_{g2}^{L}} \textrm{,}~
\dots\textrm{,}~
\frac{\partial C}{\partial z_{gN}^{L}} 
\end{bmatrix} \\
&=\begin{bmatrix}
\displaystyle
\frac{\partial C}{\partial y_{g1}^{L}} \frac{\partial y_{g1}^{L}}{\partial z_{g1}^{L}} \textrm{,}~
\frac{\partial C}{\partial y_{g2}^{L}} \frac{\partial y_{g2}^{L}}{\partial z_{g2}^{L}} \textrm{,}~
\dots\textrm{,}~
\frac{\partial C}{\partial y_{gN}^{L}} \frac{\partial y_{gN}^{L}}{\partial z_{gN}^{L}}
\end{bmatrix} \\
&=\begin{bmatrix}
\displaystyle
\frac{\partial C}{\partial y_{g1}^{L}} \dot{\phi}({z_{g1}^{L}}) \textrm{,}~
\frac{\partial C}{\partial y_{g2}^{L}} \dot{\phi}({z_{g2}^{L}}) \textrm{,}~
\dots\textrm{,}~
\frac{\partial C}{\partial y_{gN}^{L}} \dot{\phi}({z_{gN}^{L}})
\end{bmatrix} \textrm{,}
\end{split}
\end{equation}
where $\frac{\partial C}{\partial y_{gn}^{L}}$ ($1 \leq n \leq N$) is the derivative of the loss function used in the feedforward process, and $\displaystyle \dot{\phi}({z_{gn}^{L}})$ is the derivative of the activation function. The above process is summarized in Algorithm~\ref{alg:ABIPalg}.


\begin{algorithm}[t]
\caption{Arbitrary bilinear product backpropagation}
\begin{algorithmic}[1]
\REQUIRE Training inputs $\left\{\mathbf{X}_{m}\right\}_{m=1}^{M}$ \\
		\hspace{0.42cm} Training targets $\left\{\mathbf{Y}_{m}\right\}_{m=1}^{M}$ \\
		\hspace{0.42cm} Bilinear product $\bullet$ \\
		\hspace{0.42cm} Activation function $\phi$
\ENSURE  Parameters $\mathbf{\Theta}=\left\{\mathbfcal{W}^{l}, \mathbfcal{B}^{l}\right\}_{l=1}^{L}$ 

\WHILE{not converged}
\FOR{\textbf{each} minibatch  $\in \left\{\mathbf{X}\right\}_{m=1}^{M}$}
    \STATE Compute $\left\{\mathbfcal{Z}^{l},\mathbfcal{A}^{l}\right\}_{l=1}^{L}$ with Eqs.~\ref{eq:feedforward}--\ref{eq:feedforeard_at_layer_L} (feedforward) 
  \STATE Compute $\left\{\mathbfcal{D}^{l}\right\}^{L}_{l=1}$ with Eqs.~\ref{eq:delta_propagation_new} and \ref{eq:delta_calculation_at_layer_L} (backprop)
  \STATE Update $\left\{\mathbfcal{W}^{l}\right\}^{L}_{l=1}$ with the gradients from Eq.~\ref{eq:partial_C_over_partial_wij_result}
  \STATE Update $\left\{\mathbfcal{B}^{l}\right\}^{L}_{l=1}$ with the gradients from Eq.~\ref{eq:partial_C_over_partail_bi}
\ENDFOR 
\ENDWHILE 
\end{algorithmic}
\label{alg:ABIPalg}
\end{algorithm}

\newpage

%---------------------------------------------------------------
\section{Example: Circular convolution as the bilinear product}
\label{sec:ch2_circular_convolution_induction}
In the above, we generalized the vector NNs using arbitrary bilinear products.  This algorithm can be realized to train model parameters for scalar neurons or vector neurons from training data based on matrices or tensors. In this subsection, we derive $\left[\mathbf{a}_{j}^{l-1}\right]_{\bullet}^{\dagger}$ and $\left[\mathbf{w}_{ki}^{l+1}\right]_{\bullet}$ for circular convolution. The derivation is described between layer $l-1$ and layer $l$. 

Suppose $\mathbf{w}_{ij}^{l}\bullet\mathbf{a}_{ij}^{l-1}$ stands for the usual circular convolution, $\mathbf{w}_{ij}^{l}=\left[w_{ij1}^l\textrm{,}~w_{ij2}^l\textrm{,}\dots\textrm{,}~w_{ijN}^l\right]^T\in \mathbb{R}^N$ and $\mathbf{a}_{j}^{l-1}=\left[a_{j1}^{l-1}\textrm{,}~a_{j2}^{l-1}\textrm{,}\dots\textrm{,}~a_{jN}^{l-1}\right]^T\in \mathbb{R}^N$ are both $N$-dimensional vectors, and the matrix representation of $\mathbf{w}_{ij}^{l}\bullet\mathbf{a}_{ij}^{l-1}$ is:
\renewcommand{\arraystretch}{0.95}
\begin{equation}
\label{eq:circular_convolution}
\begin{split}
\left[\mathbf{w}_{ij}^{l}\right]_{\bullet}=
\begin{bmatrix}
w_{ij1}^l  & w_{ijN}^l      & w_{ij(N-1)}^l  & \dots     & w_{ij2}^l   \\ 
w_{ij2}^l  & w_{ij1}^l      & w_{ijN}^l      & \dots     & w_{ij3}^l   \\
w_{ij3}^l  & w_{ij2}^l      & w_{ij1}^l      & \dots     & w_{ij4}^l   \\
\vdots     & \vdots         & \vdots         & \ddots    & \vdots      \\
w_{ijN}^l  & w_{ij(N-1)}^l  & w_{ij(N-2)}^l  & \dots     & w_{ij1}^l   \\
\end{bmatrix} \textrm{,}
\end{split}
\end{equation}
where the weight vector is formulated as an $N \times N$ square matrix with $w_{ij1}^l$ on the main diagonal. The matrix $\left[\mathbf{a}_{j}^{l-1}\right]_{\bullet}^{\dagger}$ for Eq.~\ref{eq:partial_C_over_partial_wij_result} is extended as follows:
\renewcommand{\arraystretch}{0.95}
\begin{equation}
\begin{split}
\left[\mathbf{a}_{j}^{l-1}\right]_{\bullet}^{\dagger} =
\begin{bmatrix}
a_{j1}^{l-1}  & a_{jN}^{l-1}      & a_{j(N-1)}^{l-1}  & \dots     & a_{j2}^{l-1}  \\ 
a_{j2}^{l-1}  & a_{j1}^{l-1}      & a_{jN}^{l-1}      & \dots     & a_{j3}^{l-1}  \\
a_{j3}^{l-1}  & a_{j2}^{l-1}      & a_{j1}^{l-1}      & \dots     & a_{j4}^{l-1}  \\
\vdots        & \vdots            & \vdots            & \ddots    & \vdots        \\
a_{jN}^{l-1}  & a_{j(N-1)}^{l-1}  & a_{j(N-2)}^{l-1}  & \dots     & a_{j1}^{l-1}  \\
\end{bmatrix} \textrm{,}
\end{split}
\end{equation}
where the permutation of elements in the matrix $\left[\mathbf{a}_{j}^{l-1}\right]_{\bullet}^{\dagger}$ is identical to the matrix $\left[\mathbf{w}_{ij}^{l}\right]_{\bullet}$. The matrix $\left[\mathbf{w}_{ki}^{l+1}\right]_{\bullet}$ for circular convolution is extended as:
\renewcommand{\arraystretch}{0.95}
\begin{equation}
\begin{split}
\left[\mathbf{w}_{ki}^{l+1}\right]_{\bullet}=
\begin{bmatrix}
w_{ki1}^{l+1}  & w_{kiN}^{l+1}      & w_{ki(N-1)}^{l+1}  & \dots     & w_{ki2}^{l+1}   \\ 
w_{ki2}^{l+1}  & w_{ki1}^{l+1}      & w_{kiN}^{l+1}      & \dots     & w_{ki3}^{l+1}   \\
w_{ki3}^{l+1}  & w_{ki2}^{l+1}      & w_{ki1}^{l+1}      & \dots     & w_{ki4}^{l+1}   \\
\vdots         & \vdots             & \vdots             & \ddots    & \vdots          \\
w_{kiN}^{l+1}  & w_{ki(N-1)}^{l+1}  & w_{ki(N-2)}^{l+1}  & \dots     & w_{ki1}^{l+1}   \\
\end{bmatrix} \textrm{.}
\end{split}
\end{equation}
Note that when $N$ equals two, this architecture becomes the hyperbolic NN \cite{buchholz00ijcnn}. When $N$ is larger than two, this product is also known as polar complex multiplication \cite{olariu02elsevier,chan16tsp}. Detail investigations of more bilinear products are described in Appendix~\ref{appendix:bilinear_product_induction}.




