\section{Introduction}\label{sec:intro}

Takagi-Sugeno-Kang (TSK) fuzzy systems~\cite{nguyen2019fuzzy} have achieved great success in numerous applications, including both classification and regression problems. Many optimization approaches have been proposed for them.

There are generally three strategies for fine-tuning the TSK fuzzy system parameters after initialization: 1) evolutionary algorithms \cite{shi1999implementation,wu2006genetic}; 2) gradient descent (GD) based algorithms \cite{Wang1992b}; and, 3) GD plus least squares estimation (LSE), represented by the popular adaptive-network-based fuzzy inference system (ANFIS) \cite{jang1993anfis}. However, these approaches may have challenges when the size and/or the dimensionality of the data increase. Evolutionary algorithms need to keep a large population of candidate solutions, and evaluate the fitness of each, which result in high computational cost and heavy memory requirement for big data. Traditional GD needs to compute the gradients from the entire dataset to iteratively update the model parameters, which may be very slow, or even impossible, when the data size is very large. The memory requirement and computational cost of LSE also increase rapidly when the data size and/or dimensionality increase. Additionally, as shown in \cite{drwuGD2019}, ANFIS may result in significant overfitting in regression problems.

Many efforts have been spent to tackling the difficulty in optimizing the TSK fuzzy systems on big and/or high-dimensional data~\cite{jin2000fuzzy, deng2016hierarchical,chung2008minimum}. Dimensionality reduction and/or feature selection are usually used to reduce the number of fuzzy partitions (rules). Traditional dimensionality reduction techniques such as principal component analysis (PCA) has been used for TSK fuzzy system optimization~\cite{nilashi2015multi, lau2013fault}. There are also methods focusing on learning a sparse subspace of the original feature space to reduce the number of antecedents in each rule~\cite{deng2016survey,deng2010enhanced}. Once the number of antecedents is determined, different optimization approaches can be used to tune the TSK fuzzy system on large datasets. For example, Chung \emph{et al.} \cite{chung2008minimum} utilized the equivalence between minimum enclosing ball and the Mamdani-Larsen fuzzy inference system to train the latter using the former. Gacto \emph{et al.} \cite{gacto2014metsk} proposed a multi-objective evolutionary algorithm to optimize TSK fuzzy systems for high-dimensional large-scale regression problems.

Mini-batch gradient descent (MBGD)~\cite{goodfellow2016deep,Ruder2016} based optimization, which is particularly popular in deep learning, can also be a solution to training TSK fuzzy systems on large and high-dimensional datasets. In each iteration, MBGD computes the gradients from a randomly selected small batch of data, instead of the entire dataset~\cite{bottou2010large}. Different batch sizes can be used, according to the trade-off among the available memory, the training speed, and the expected generalization performance. The original MBGD used a constant learning rate to update the model's parameters~\cite{bottou2010large}. Later, Sutskever \emph{et al.}~\cite{sutskever2013importance} found that adding a momentum to MBGD can improve the final training performance. However, it still needs to manually select a learning rate, and the convergence may be very slow at the beginning. Kingma and Ba~\cite{kingma2015adam} proposed the well-known Adam algorithm to automatically rescale the gradients to achieve adaptive and individualized learning rate for each parameter, which leads to faster convergence. However, the generalization performance of Adam may not be as good as the momentum~\cite{wilson2017marginal}; so, Keskar and Socher~\cite{keskar2017improving} also tried to combine the advantages of momentum and Adam to achieve both fast convergence and good generalization. Recently, Luo \emph{et al.}~\cite{luo2019adaptive} also proposed AdaBound to improve Adam. AdaBound uses an adaptive bound for the learning rate of each parameter to force the optimizer to behave like Adam at the beginning and like stochastic GD at the end. Our very recent research \cite{drwuGD2019} has found that TSK fuzzy systems can achieve better performance with AdaBound than Adam for regression problems.

Although MBGD-based optimization has many advantages, it may be easily trapped into a local-minimum, and may face the gradient vanishing problem. Many other techniques have been proposed to complement MBGD for better performance. In 2015, Ioffe and Szegedy~\cite{ioffe2015batch} proposed the well-known batch normalization (BN) approach to accelerate the training of deep neural networks by reducing the internal covariate shift\footnote{Recently some researchers had different opinions on why BN works. For example, Santurkar \emph{et al.}~\cite{santurkar2018does} argued that BN may not reduce the internal covariate shift; instead, it helps improve the Lipschitzness of both the loss and the gradients, and also reduces the dependency on the training hyper-parameters, such as the learning rate and the regularization weights.}. BN normalizes the input distribution of each layer, so it also alleviates the gradient vanishing problem. It has been used almost ubiquitously in deep learning, and many variants \cite{ba2016layer,fan2017revisit,clevert2015fast,wu2018group} have also been proposed.

This paper, following our previous research \cite{drwuGD2019} on MBGD-based optimization of TSK fuzzy systems for regression problems, considers classification problems. We use AdaBound, as in \cite{drwuGD2019}, to adjust the learning rates. Additionally, we propose two novel techniques for training TSK fuzzy systems for classification problems, namely, uniform regularization (UR) and BN. Our main contributions are:
\begin{enumerate}
\item We introduce a novel UR term to the cross-entropy loss function in training TSK fuzzy classifiers, which forces all rules to have similar average firing levels on the entire dataset. Experiments show that UR can improve the generalization performance of TSK fuzzy classifiers.
\item We extend BN from the training of deep neural networks to the training of TSK fuzzy classifiers, and show that it can speed up the convergence in training and improve the generalization performance in testing.
\item We further integrate UR and BN, and show that the combined approach outperforms each individual ones.
\end{enumerate}

The remainder of this paper is organized as follows: Section~\ref{sec:methods} introduces the proposed UR and BN approaches. Section~\ref{sec:res} presents the experimental results to validate the performances of UR and BN. Section~\ref{sec:conclusion} draws conclusions and points out some future research directions.

\section{UR and BN}\label{sec:methods}

This section introduces the details of the TSK fuzzy classifier under consideration, our proposed UR for regularizing the loss function, and BN for more efficient and effective training of the TSK fuzzy classifier. Python implementation of our algorithm can be downloaded at \emph{https://github.com/YuqiCui/TSK\_BN\_UR}.

\subsection{The TSK Fuzzy Classifier}

Let the training dataset be $\mathcal{D} = \{\bm{x}_n, y_n\}_{n=1}^{N}$, in which $\bm{x}_n=[x_{n,1},...,x_{n,D}]^T\in \mathbb{R}^{D \times 1}$ is a $D$-dimensional feature vector, and $y_n \in \{1,2,...,C\}$ the corresponding class label for a $C$-class classification problem.

Suppose the TSK fuzzy classifier has $R$ rules, in the following form:
\begin{align}
	\begin{split}
		\textup{Rule}_r:~&\textup{IF}~x_1~\textup{is}~X_{r, 1}~\textup{and}~ \cdots ~\textup{and}~x_D~\textup{is}~X_{r, D}, \\
		&\textup{THEN}~y_r^1(\bm{x}) = b_{r,0}^1+\sum_{d=1}^{D}b_{r,d}^1\cdot x_d~\textup{and} \cdots \\
&\qquad\qquad \textup{and}~y_r^C(\bm{x}) = b_{r,0}^C+\sum_{d=1}^{D}b_{r,d}^C\cdot x_d
	\end{split} \label{eq:rule}
\end{align}
where $X_{r, d}$ ($r=1,...,R$; $d=1,...,D$) is the membership function (MF) for the $d$-th antecedent in the $r$-th rule, and $b_{r,0}^c$ and $b_{r,d}^c$ ($c=1,...,C$) are the consequent parameters for the $c$-th class.

Different types of MFs can be used in our algorithm, as long as they are differentiable. For simplicity, Gaussian MFs are considered in this paper, and the membership grade of $x_d$ on $X_{r,d}$ is:
\begin{align}
	\mu_{X_{r, d}}(x_d) &= \exp\left(-\frac{(x_d-m_{r,d})^2}{2\sigma_{r,d}^2}\right), \label{eq:MG}
\end{align}
where $m_{r,d}$ and $\sigma_{r,d}$ are the center and the standard deviation of the Gaussian MF, respectively.

The output of the TSK fuzzy classifier for the $c$-th class is:
\begin{align}
y^c(\bm{x})=\frac{\sum_{r=1}^Rf_r(\bm{x})y_r^c(\bm{x})}{\sum_{r=1}^Rf_r(\bm{x})}, \label{eq:yx}
\end{align}
where
\begin{align}
f_r(\bm{x})=\prod_{d=1}^D \mu_{X_{r,d}}(x_d) =\textup{exp}\left(-\sum_{d=1}^D\frac{(x_d-m_{r,d})^2}{2\sigma_{r,d}^2}\right) \label{eq:prod}
\end{align}
is the firing level of Rule~$r$. We can also re-write (\ref{eq:yx}) as:
\begin{align}
y^c(\bm{x})=\sum_{r=1}^R\overline{f}_r(\bm{x})y_r^c(\bm{x}),
\end{align}
where
\begin{align}
\overline{f}_r(\bm{x})=\frac{f_r(\bm{x})}{\sum_{i=1}^R f_i(\bm{x})}
\end{align}
is the normalized firing level of Rule~$r$.

Once the output vector $\bm{y}(\bm{x})=[y^1(\bm{x}),...,y^C(\bm{x})]^T$ is obtained, the input $\bm{x}$ is assigned to the class with the largest $y^c(\bm{x})$.

To optimize the TSK fuzzy classifier, we need to fine-tune the antecedent MF parameters $m_{r,d}$ and $\sigma_{r,d}$, and the consequent parameters $b_{r,0}^c$ and $b_{r,d}^c$, where $r=1,...,R$, $d=1,...,D$, and $c=1,...,C$.

\subsection{Uniform Regularization (UR)}\label{subsec:ur}

Mixture of experts (MoE)~\cite{Jacobs1991}, which is functionally equivalent to TSK fuzzy systems~\cite{Bersini1997,Andersen1998,wu2019functional}, is a popular machine learning algorithm. Its model is shown in Fig.~\ref{fig:MoE}. It trains multiple local experts, each taking care of only a small local region of the input space. For a new input, the gating network determines the activations (weights) of the local experts, and the final output is a weighted average of the local expert outputs.

\begin{figure}[htbp]\centering
\includegraphics[width=.85\linewidth,clip]{MoE}
\caption{Mixture of experts (MoE)~\cite{Jacobs1991}.} \label{fig:MoE}
\end{figure}

Although MoE has been used successfully in many applications, it may suffer from the ``rich get richer" effect \cite{shen2019mixture,shazeer2017outrageously}: once an expert is slightly better than others, it is always picked by the gating network, whereas other experts starve and are rarely used. This is bad for the generalization performance of the overall model.

Since MoE and TSK fuzzy systems are functionally equivalent~\cite{wu2019functional}, TSK fuzzy systems may also suffer from the ``rich get richer" effect, i.e., only a few rules are always activated with large firing levels, whereas others have very small firing levels, and hence not adequately tuned in training. A remedy to the ``rich get richer" effect in TSK fuzzy systems is to force the rules to be fired at similar degrees in the input space, so that each rule contributes about equally to the output.

Next, we propose UR to achieve this goal.

UR forces the rules to have similar average firing levels, by minimizing the following loss:
\begin{align}
	\ell_{UR} = \sum_{r=1}^R\left(\frac{1}{N}\sum_{n=1}^N \overline{f}_r(\bm{x}_n) - \tau\right)^2,
\end{align}
where $N$ is the number of training examples, and $\tau$ the expected firing level of each rule, which is set to $1/C$ in this paper (recall that $C$ is the number of classes).

$\ell_{UR}$ can then be added to the original loss function in MBGD-based training of TSK fuzzy classifiers, i.e., for each mini-batch with $N$ training samples,
\begin{align}
	\mathcal{L} = \ell + \alpha \ell_2 + \lambda \sum_{r=1}^R\left(\frac{1}{N}\sum_{n=1}^N\overline{f}_r(\bm{x}_n)-\frac{1}{R}\right)^2, \label{eq:ur}
\end{align}
where $\ell$ is the cross-entropy loss between the estimated class probabilities [obtained by applying \emph{softmax} to $\bm{y}(\bm{x})$] and the true class probabilities, $\ell_2$ the L2 regularization of the rule consequent parameters, and $\alpha$ and $\lambda$ the trade-off parameters.

\subsection{Batch Normalization (BN)}\label{subsec:bn}

BN \cite{ioffe2015batch} is a very powerful technique in optimizing deep neural networks~\cite{he2016deep,zagoruyko2016wide,huang2017densely}. It normalizes the data distribution in each mini-batch to accelerate the training. For a mini-batch $\mathcal{B}=\{\bm{x}_n\}_{n=1}^N$, the output of BN is \cite{ioffe2015batch}:
\begin{align}
\bm{x}'_n = BN(\bm{x}_n) = \gamma \frac{\bm{x}_n-\bm{m}_\mathcal{B}}{\sqrt{\bm{\sigma}_\mathcal{B}^2+\epsilon}}+\beta, \label{eq:bn}
\end{align}
where $\bm{m}_\mathcal{B}$ and $\bm{\sigma}_\mathcal{B}$ are the mean and the standard deviation of the samples in the mini-batch, respectively, $\gamma$ and $\beta$ are parameters to be learned during training, and $\epsilon$ is usually set to $1e-8$ to avoid being divided by zero. During training, exponential weighted averages of $\bm{m}_\mathcal{B}$ and $\bm{\sigma}_\mathcal{B}$ are recorded so that they can be used in the test phase.

Since TSK fuzzy systems and neural networks share lots of similarity \cite{wu2019functional}, we can extend BN to the optimization of TSK fuzzy classifiers, as shown in Fig.~\ref{fig:bn}. In the training phase, we first compute the firing level of each rule using the unmodified inputs, as in traditional TSK fuzzy systems. Then, we use BN to normalize the inputs, according to their mean and standard deviation in the current mini-batch. The normalized inputs are then used to compute the rule consequents. The final output is a weighted average of the rule consequents, the weights being the corresponding rule firing levels.

\begin{figure}[htpb] \centering
\includegraphics[width=.85\linewidth,clip]{BN1}
	\caption{BN in training a TSK fuzzy classifier. All rule consequents share the same BN layer.}
	\label{fig:bn}
\end{figure}

At the testing phase, the BN operation can be merged into the consequent layer. Assume that after training, we obtain a BN layer with learned $\bm{m}=(m_1,...,m_D)^T$, $\bm{\sigma}=(\sigma_1,...,\sigma_D)^T$, $\gamma$ and $\beta$. Then, the output $y_r$ of the $r$-th rule with BN is:
\begin{align}
	y_r(BN(\bm{x}_n)) &= b_{r,0}+ \gamma\sum_{d=1}^D b_{r,d}\frac{x_{n,d}-m_d}{\sqrt{\sigma_d^2+\epsilon}} + \beta D,
\end{align}
which can be re-written as:
\begin{align}
	y_r(BN(\bm{x}_n)) = b'_{r,0}+\sum_{d=1}^D b'_{r,d}x_{n,d}, \label{eq:bnrew}
\end{align}
where
\begin{align}
b'_{r,0} &= b_{r,0}+\beta D
-\gamma\sum_{d=1}^D\frac{m_db_{r,d}}{\sqrt{\sigma_d^2+\epsilon}},\\
b'_{r,d}& = \gamma \frac{b_{r,d}}{\sqrt{\sigma_d^2+\epsilon}}.
\end{align}
By doing this, the original architecture of the TSK fuzzy classifier is kept unchanged.

We also tested two variants of BN, as shown in Fig.~\ref{fig:bnv}. The TSK with global BN (\texttt{TSK-MBGD-UR-GBN}) approach in Fig.~\ref{fig:bnv2} uses the BN normalized inputs in both antecedents and consequents to compute the final output. In this case, the output of \texttt{TSK-MBGD-UR-GBN} for Class~$c$ is:
\begin{align}
	y^c(\bm{x})=\sum_{r=1}^R\overline{f}_r(BN(\bm{x}))y_r^c(BN(\bm{x})).
\end{align}
The TSK with rule-specific BN (\texttt{TSK-MBGD-UR-RBN}) approach in Fig.~\ref{fig:bnv3} uses the raw inputs to compute the antecedents, and rule-specific BN to compute each consequent individually. The output of \texttt{TSK-MBGD-UR-RBN} for Class~$c$ is:
\begin{align}
	y^c(\bm{x})=\sum_{r=1}^R\overline{f}_r(\bm{x})y_r^c(BN_r(\bm{x})),
\end{align}
where $BN_r$ represents the BN operation for the $r$-th rule.

\texttt{TSK-MBGD-UR-GBN} has the same computational cost as \texttt{TSK-MBGD-UR-BN}, but \texttt{TSK-MBGD-UR-RBN} has $R$ times more BN parameters, and hence higher computational cost. Both of them can be re-expressed in the original TSK architecture. We also evaluate their performances in Section~\ref{subsec:bnres}.

\begin{figure}[htpb] \centering
	\subfigure[]{\includegraphics[width=0.9\columnwidth,clip]{BN2}\label{fig:bnv2}}
	\subfigure[]{\includegraphics[width=0.9\columnwidth,clip]{BN3}\label{fig:bnv3}}
	\caption{(a) TSK fuzzy system with global BN (\texttt{TSK-MBGD-UR-GBN}); and, (b) TSK fuzzy system with rule-specific BN (\texttt{TSK-MBGD-UR-RBN}).} 	\label{fig:bnv}
\end{figure}

\section{Experiments and Results}\label{sec:res}

This section validates the performances of our proposed UR and BN on multiple datasets from various application domains, with varying size and feature dimensionality.

\subsection{Datasets}

We evaluated our proposed algorithms on 12 classification datasets from the UCI Machine Learning Repository\footnote{http://archive.ics.uci.edu/ml/index.php}. Their characteristics are summarized in Table~\ref{tab:data}. For each dataset, we randomly selected 70\% samples as the training set and the remaining 30\% as the test set for 30 times to get 30 different data splits. We ran each algorithm on these 30 data splits and report the average performance.

\begin{table}[htpb]
\begin{threeparttable} \setlength{\tabcolsep}{1.5mm}{
\caption{Summary of the 12 datasets.} \label{tab:data}
\begin{tabular}{ccccc} \hline
Index & Dataset                     & No. of Samples & No. of Features & No. of Classes \\ \hline
1 & Vehicle\tnote{1}                     & 846            & 18              & 4              \\
2 & Biodeg\tnote{2}                      & 1,055           & 41              & 2              \\
3 & DRD\tnote{3} & 1151           & 19              & 2              \\
4 & Yeast\tnote{4}                       & 1,484           & 8               & 10             \\
5 & Steel\tnote{5}                       & 1,941           & 27              & 7              \\
6 & IS\tnote{6}           & 2,310           & 19              & 7              \\
7 & Abalone\tnote{7}                     & 4,177           & 10              & 3              \\
8 & Waveform21\tnote{8}                  & 5,000           & 21              & 3              \\
9 & Page-blocks\tnote{9}                & 5,473           & 10              & 5              \\
10 & Satellite\tnote{10}                   & 6,435           & 36              & 6              \\
11 & Clave\tnote{11}                       & 10,798          & 16              & 4              \\
12 & MAGIC\tnote{12}                       & 19,020          & 10              & 2              \\ \hline
\end{tabular}}
\begin{tablenotes}
\item[1] https://archive.ics.uci.edu/ml/datasets/Statlog+\%28Vehicle+Silhouettes\%29
\item[2] https://archive.ics.uci.edu/ml/datasets/QSAR+biodegradation
\item[3] https://archive.ics.uci.edu/ml/datasets/Diabetic+Retinopathy+Debrecen +Data+Set
\item[4] https://archive.ics.uci.edu/ml/datasets/Yeast
\item[5] https://archive.ics.uci.edu/ml/datasets/Steel+Plates+Faults
\item[6] https://archive.ics.uci.edu/ml/datasets/Image+Segmentation
\item[7] https://archive.ics.uci.edu/ml/datasets/Abalone
\item[8] https://archive.ics.uci.edu/ml/datasets/Waveform+Database+Generator +(Version+1)
\item[9] https://archive.ics.uci.edu/ml/datasets/Page+Blocks+Classification
\item[10] https://archive.ics.uci.edu/ml/datasets/Statlog+(Landsat+Satellite)
\item[11] https://archive.ics.uci.edu/ml/datasets/Firm-Teacher\_Clave-Direction\_Classification
\item[12] https://archive.ics.uci.edu/ml/datasets/MAGIC+Gamma+Telescope
\end{tablenotes}
\end{threeparttable}
\end{table}

Some datasets contain both numerical features and categorical features. The categorical features were converted into numerical ones by one-hot coding. We $z$-normalized each feature using the mean and standard deviation computed from the training set.

\subsection{Algorithms}

We compared nine algorithms to validate our proposed approaches. Among them, four were tree based approaches (\texttt{DT}, \texttt{RF}, \texttt{PART}, and \texttt{JRip}), one was a TSK fuzzy system optimized by a traditional approach (\texttt{TSK-FCM-LSE}), and the remaining four were TSK fuzzy systems optimized by MBGD based approaches (\texttt{TSK-MBGD}, \texttt{TSK-MBGD-BN}, \texttt{TSK-MBGD-UR}, \texttt{TSK-MBGD-UR-BN}).

The details of these nine algorithms are as follows:
\begin{enumerate}
\item \texttt{DT}: Decision tree implemented in scikit-learn\footnote{https://scikit-learn.org/stable/modules/generated/sklearn.tree.DecisionTree Classifier.html} in Python. We used 5-fold cross-validation to select the maximum depth of the tree from $\{3, 4, 5, 6, 7\}$ on the training set. Other parameters were set by default.
\item \texttt{RF}: Random forest implemented in scikit-learn\footnote{https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.Random ForestClassifier.html} in Python. We set the number of trees to 20 and used 5-fold cross-validation to select the maximum depth of the trees from $\{3, 4, 5, 6, 7\}$ on the training set. Other parameters were set by default.
\item \texttt{PART}~\cite{frank1998generating}: The PART (partial decision tree) classifier implemented in RWeka\footnote{https://cran.r-project.org/web/packages/RWeka/index.html}. All parameters were set by default.
\item \texttt{JRip}~\cite{cohen1995repeated}: The RIPPER (Repeated Incremental Pruning to Produce Error Reduction) classifier implemented in RWeka. All parameters were set by default.
\item \texttt{TSK-FCM-LSE}~\cite{jang1997neuro}: We used fuzzy $c$-means (FCM) clustering to estimate the antecedent parameters, and LSE with L2 regularization to estimate the consequent parameters.
\item \texttt{TSK-MBGD}: We used MBGD and AdaBound~\cite{luo2019adaptive} to optimize both the antecedent and the consequent parameters.
\item \texttt{TSK-MBGD-UR}: We used MBGD, AdaBound and UR (Section~\ref{subsec:ur}) to optimize both the antecedent and the consequent parameters. The UR weight $\lambda$ in (\ref{eq:ur}) was selected from $\{0.1, 1, 10, 20, 50\}$ by cross-validation on the training set.
\item \texttt{TSK-MBGD-BN}: We used MBGD, AdaBound and BN (Section~\ref{subsec:bn}) to optimize both the antecedent and the consequent parameters.
\item \texttt{TSK-MBGD-UR-BN}: We used MBGD, AdaBound, BN and UR to optimize both the antecedent and the consequent parameters. The UR weight $\lambda$ in (\ref{eq:ur}) was selected from $\{0.1, 1, 10, 20, 50\}$ by cross-validation on the training set.
\end{enumerate}

For \texttt{TSK-FCM-LSE}, \texttt{TSK-MBGD}, \texttt{TSK-MBGD-BN}, \texttt{TSK-MBGD-UR} and \texttt{TSK-MBGD-UR-BN}, we set the L2 regularization weight $\alpha=0.05$, and the number of rules $R=20$. For \texttt{TSK-MBGD}, \texttt{TSK-MBGD-BN}, \texttt{TSK-MBGD-UR} and \texttt{TSK-MBGD-UR-BN}, we set the learning rate of AdaBound to 0.01, following our previous work \cite{drwuGD2019}. In order to make use of all data in the training set and to reduce overfitting simultaneously, we randomly sampled 20\% data from the training set and trained the TSK model with early stopping five times. The maximum epoch number was 2,000, and the patience of early stopping 40. We recorded the number of epochs at stopping in each run, and trained the final model with the average stopping epoch number on the entire training set.

$k$-mean clustering was used in the MBGD-based algorithms (\texttt{TSK-MBGD}, \texttt{TSK-MBGD-BN}, \texttt{TSK-MBGD-UR}, and \texttt{TSK-MBGD-UR-BN}) to initialize the antecedent parameters. We performed $k$-means clustering on the training set, where $k$ equaled $R$, the number of rules. We then initialized the rule centers to the cluster centers, and randomly initialized the standard deviation $\sigma_{r,d}$ from a Gaussian distribution $\mathcal{N}(1, 0.2)$. For the consequent parameters, we set the initial bias of each rule to zero, and the attribute weight $b_{r,d}$ ($r=1, ..., R$; $d=1, ..., D$) randomly from a uniform distribution $U(-1, 1)$.

\subsection{Performance Measures}

The raw classification accuracy (RCA), which is the total number of correctly classified test samples divided by the total number of test samples, was used as our primary performance measure.

Since some datasets have significant class imbalance, in addition to the RCA, we also computed the balanced classification accuracy (BCA), which is the mean of the per-class RCAs, as our second performance measure.

\subsection{Experimental Results}

The average test RCAs and BCAs are shown in Tables~\ref{tab:acc} and \ref{tab:bca}, respectively. The largest value (best performance) on each dataset is marked in bold. To facilitate the comparison, we also show the ranks of the RCAs and BCAs in Tables~\ref{tab:accrank} and \ref{tab:bcarank}, respectively.

The following observations can be made from the above four tables:
\begin{enumerate}
	\item \emph{Generally, UR improved both RCA and BCA.} Comparing \texttt{TSK-MBGD} with \texttt{TSK-MBGD-UR}, and \texttt{TSK-MBGD-BN} with \texttt{TSK-MBGD-UR-BN}, we can conclude that generally UR improved the classification performance, regardless of whether BN was used or not. The average ranks in the last row of Tables~\ref{tab:accrank} and \ref{tab:bcarank} demonstrate this more clearly: the average rank of \texttt{TSK-MBGD-UR} (\texttt{TSK-MBGD-UR-BN}) was smaller than that of \texttt{TSK-MBGD} (\texttt{TSK-MBGD-BN}).

	\item \emph{Generally, BN improved both RCA and BCA.} Comparing \texttt{TSK-MBGD} with \texttt{TSK-MBGD-BN}, and \texttt{TSK-MBGD-UR} with \texttt{TSK-MBGD-UR-BN}, we can conclude that generally BN improved the classification performance, regardless of whether UR was used or not. The average ranks in the last row of Tables~\ref{tab:accrank} and \ref{tab:bcarank} demonstrate this more clearly: the average rank of \texttt{TSK-MBGD-BN} (\texttt{TSK-MBGD-UR-BN}) was smaller than that of \texttt{TSK-MBGD} (\texttt{TSK-MBGD-UR}).

	\item \emph{Generally, integrating BN and UR achieved further RCA and BCA improvements.} Comparing \texttt{TSK-MBGD-UR-BN} with \texttt{TSK-MBGD}, \texttt{TSK-MBGD-UR} and \texttt{TSK-MBGD-BN}, we can conclude that \texttt{TSK-MBGD-UR-BN} almost always performed the best on both RCA and BCA, as shown in Fig.~\ref{fig:improv}. This indicated that BN and UR are somehow complementary, and hence integrating them may achieve better performance than using each one alone.

	\item \emph{Overall, \texttt{TSK-MBGD-UR-BN} achieved the best performance among the nine algorithms.}  The last row of Table~\ref{tab:bcarank} shows that \texttt{TSK-MBGD-UR-BN} achieved the best average BCA performance, and the last row of Table~\ref{tab:accrank} shows that \texttt{TSK-MBGD-UR-BN} achieved the second best average RCA performance. Interestingly, \texttt{RF} had the best average rank on RCA, but only ranked the fifth on BCA, suggesting that \texttt{RF} may tend to overlook the minority classes. On the contrary, \texttt{TSK-MBGD-UR-BN} performed well on both RCA and BCA.
\end{enumerate}

\begin{table*}[htpb]\centering \setlength{\tabcolsep}{1.6mm}
\caption{Average RCAs of the nine algorithms on the 12 datasets.} \label{tab:acc}
\begin{tabular}{cccccccccc}\hline
Dataset                     & \texttt{CART}   & \texttt{RF}              & \texttt{JRip}   & \texttt{PART}            & \texttt{TSK-FCM-LSE}             & \texttt{TSK-MBGD}    & \texttt{TSK-MBGD-BN}          & \texttt{TSK-MBGD-UR} & \texttt{TSK-MBGD-UR-BN}       \\ \hline
Vehicle                     & 0.6907 & 0.7407          & 0.6892 & 0.7110          & 0.7411          & 0.6970 & 0.7354          & 0.7089 & \textbf{0.7907} \\
Biodeg                      & 0.8202 & 0.8572          & 0.8222 & 0.8362          & 0.8377          & 0.8523 & 0.8531          & 0.8539 & \textbf{0.8609} \\
DRD & 0.6283 & 0.6589          & 0.6240 & 0.6364          & \textbf{0.6824} & 0.6623 & 0.6618          & 0.6713 & 0.6720          \\
Yeast                       & 0.5564 & \textbf{0.5963} & 0.5731 & 0.5340          & 0.5851          & 0.5673 & 0.5770          & 0.5722 & 0.5725          \\
Steel                       & 0.7017 & 0.7328          & 0.7135 & 0.7120          & 0.6527          & 0.5864 & 0.7110          & 0.7248 & \textbf{0.7350} \\
IS           & 0.9320 & 0.9529          & 0.9481 & \textbf{0.9608} & 0.9571          & 0.5762 & 0.7557          & 0.8559 & 0.9501          \\
Abalone                     & 0.7170 & 0.7314          & 0.7254 & 0.7104          & \textbf{0.7323} & 0.5821 & 0.7129          & 0.6238 & 0.7306          \\
Waveform21                  & 0.7641 & 0.8369          & 0.7908 & 0.7843          & \textbf{0.8647} & 0.6779 & 0.8002          & 0.8363 & 0.8234          \\
Page-blocks                 & 0.9651 & \textbf{0.9688} & 0.9681 & 0.9677          & 0.9499          & 0.9375 & 0.9419          & 0.9515 & 0.9580          \\
Satellite                   & 0.8524 & 0.8863          & 0.8587 & 0.8592          & 0.8864          & 0.4890 & 0.8001          & 0.8929 & \textbf{0.8943} \\
Clave                       & 0.7103 & 0.7600          & 0.7344 & 0.7779          & 0.7690          & 0.8223 & \textbf{0.8427} & 0.8187 & 0.8192          \\
MAGIC                       & 0.8427 & 0.8531          & 0.8455 & 0.8488          & 0.8319          & 0.7347 & 0.7861          & 0.8574 & \textbf{0.8392} \\ \hline
Average                        & 0.7651 & 0.7979          & 0.7744 & 0.7782          & 0.7909          & 0.6821 & 0.7648          & 0.7806 & \textbf{0.8038} \\\hline
\end{tabular}
\end{table*}

\begin{table*}[htpb]\centering \setlength{\tabcolsep}{1.6mm}
\caption{Average BCAs of the nine algorithms on the 12 datasets.} \label{tab:bca}
\begin{tabular}{cccccccccc}
\hline
Dataset                     & \texttt{CART}   & \texttt{RF}              & \texttt{JRip}   & \texttt{PART}            & \texttt{TSK-FCM-LSE}             & \texttt{TSK-MBGD}    & \texttt{TSK-MBGD-BN}          & \texttt{TSK-MBGD-UR} & \texttt{TSK-MBGD-UR-BN}       \\ \hline
Vehicle           & 0.6936 & 0.744  & 0.6939          & 0.7131          & 0.7443          & 0.7010 & 0.7380 & 0.7127          & \textbf{0.7930} \\
Biodeg            & 0.7973 & 0.8306 & 0.7899          & 0.8122          & 0.8205          & 0.8368 & 0.8318 & 0.8390          & \textbf{0.8439} \\
DRD               & 0.634  & 0.6624 & 0.6227          & 0.6422          & \textbf{0.6845} & 0.6642 & 0.6634 & 0.6717          & 0.6729          \\
Yeast             & 0.3998 & 0.4867 & 0.5203          & 0.4889          & 0.5102          & 0.4951 & 0.5184 & 0.4946          & \textbf{0.5332} \\
Steel             & 0.7005 & 0.6937 & 0.7129          & 0.7267          & 0.6319          & 0.5933 & 0.7258 & 0.7245          & \textbf{0.7515} \\
IS & 0.932  & 0.9529 & 0.9481          & \textbf{0.9607} & 0.9571          & 0.5762 & 0.7557 & 0.8559          & 0.9501          \\
Abalone           & 0.5319 & 0.5362 & 0.5371          & 0.5280          & 0.5402          & 0.4567 & 0.5236 & 0.4791          & \textbf{0.5402} \\
Waveform21        & 0.7637 & 0.8365 & 0.7905          & 0.7844          & \textbf{0.8645} & 0.6784 & 0.8003 & 0.8362          & 0.8233          \\
Page-blocks       & 0.7986 & 0.7385 & \textbf{0.8192} & 0.8162          & 0.6003          & 0.5129 & 0.5609 & 0.6033          & 0.671           \\
Satellite         & 0.8204 & 0.8480 & 0.8308          & 0.834           & 0.8558          & 0.4337 & 0.7651 & 0.8679          & \textbf{0.8700} \\
Clave             & 0.4701 & 0.4878 & 0.4985          & \textbf{0.6507} & 0.4825          & 0.5876 & 0.6468 & 0.6374          & 0.6421          \\
MAGIC             & 0.8058 & 0.8108 & 0.8052          & 0.8135          & 0.7886          & 0.6325 & 0.7128 & \textbf{0.8225} & 0.7934          \\ \hline
Average               & 0.6956 & 0.7190 & 0.714           & 0.7309          & 0.7067          & 0.5974 & 0.6869 & 0.7120          & \textbf{0.7404} \\ \hline
\end{tabular}
\end{table*}

\begin{table*}[htpb]\centering
\caption{RCA ranks of the nine algorithms on the 12 datasets.}\label{tab:accrank}
\begin{tabular}{cccccccccc}\hline
Dataset                     & \texttt{CART}   & \texttt{RF}              & \texttt{JRip}   & \texttt{PART}            & \texttt{TSK-FCM-LSE}             & \texttt{TSK-MBGD}    & \texttt{TSK-MBGD-BN}          & \texttt{TSK-MBGD-UR} & \texttt{TSK-MBGD-UR-BN}       \\ \hline
Vehicle                     & 8    & 3   & 9    & 5    & 2   & 7   & 4      & 6      & 1         \\
Biodeg                      & 9    & 2   & 8    & 7    & 6   & 5   & 4      & 3      & 1         \\
DRD & 8    & 6   & 9    & 7    & 1   & 4   & 5      & 3      & 2         \\
Yeast                       & 8    & 1   & 4    & 9    & 2   & 7   & 3      & 6      & 5         \\
Steel                       & 7    & 2   & 4    & 5    & 8   & 9   & 6      & 3      & 1         \\
IS           & 6    & 3   & 5    & 1    & 2   & 9   & 8      & 7      & 4         \\
Abalone                     & 5    & 2   & 4    & 7    & 1   & 9   & 6      & 8      & 3         \\
Waveform21                  & 8    & 2   & 6    & 7    & 1   & 9   & 5      & 3      & 4         \\
Page-blocks                 & 4    & 1   & 2    & 3    & 7   & 9   & 8      & 6      & 5         \\
Satellite                   & 7    & 4   & 6    & 5    & 3   & 9   & 8      & 2      & 1         \\
Clave                       & 9    & 7   & 8    & 5    & 6   & 2   & 1      & 4      & 3         \\
MAGIC                       & 5    & 2   & 4    & 3    & 7   & 9   & 8      & 1      & 6         \\\hline
Average                        & 7.0  & 2.9 & 5.8  & 5.3  & 3.8 & 7.3 & 5.5    & 4.3    & 3.0      \\\hline
\end{tabular}
\end{table*}

\begin{table*}[htbp]\centering
\caption{BCA ranks of the nine algorithms on the 12 datasets.}\label{tab:bcarank}
\begin{tabular}{cccccccccc}
\hline
Dataset                     & \texttt{CART}   & \texttt{RF}              & \texttt{JRip}   & \texttt{PART}            & \texttt{TSK-FCM-LSE}             & \texttt{TSK-MBGD}    & \texttt{TSK-MBGD-BN}          & \texttt{TSK-MBGD-UR} & \texttt{TSK-MBGD-UR-BN}       \\ \hline
Vehicle                     & 9    & 3   & 8    & 5    & 2   & 7   & 4      & 6      & 1         \\
Biodeg                      & 8    & 5   & 9    & 7    & 6   & 3   & 4      & 2      & 1         \\
DRD & 8    & 6   & 9    & 7    & 1   & 4   & 5      & 3      & 2         \\
Yeast                       & 9    & 8   & 2    & 7    & 4   & 5   & 3      & 6      & 1         \\
Steel                       & 6    & 7   & 5    & 2    & 8   & 9   & 3      & 4      & 1         \\
IS           & 6    & 3   & 5    & 1    & 2   & 9   & 8      & 7      & 4         \\
Abalone                     & 5    & 4   & 3    & 6    & 1   & 9   & 7      & 8      & 2         \\
Waveform21                  & 8    & 2   & 6    & 7    & 1   & 9   & 5      & 3      & 4         \\
Page-blocks                 & 3    & 4   & 1    & 2    & 7   & 9   & 8      & 6      & 5         \\
Satellite                   & 7    & 4   & 6    & 5    & 3   & 9   & 8      & 2      & 1         \\
Clave                       & 9    & 7   & 6    & 1    & 8   & 5   & 2      & 4      & 3         \\
MAGIC                       & 4    & 3   & 5    & 2    & 7   & 9   & 8      & 1      & 6         \\ \hline
Average                        & 6.8  & 4.7 & 5.4  & 4.3  & 4.2 & 7.3 & 5.4    & 4.3    & 2.6       \\ \hline
\end{tabular}
\end{table*}

\begin{figure}
	\subfigure[]{\includegraphics[width=0.9\columnwidth,clip]{acc_improv}}
	\subfigure[]{\includegraphics[width=0.9\columnwidth,clip]{bca_improv}}
	\caption{(a) RCAs and (b) BCAs of the four MBGD-based TSK fuzzy classifiers on the 12 datasets. Datasets were sorted according to the RCAs of the \texttt{TSK-MBGD} model. The indices along the horizontal axis denote the dataset indices in Table~\ref{tab:data}.} \label{fig:improv}
\end{figure}

\subsection{Statistical Analysis}

To further evaluate the performance improvement of our proposed \texttt{TSK-MBGD-UR-BN} over others, we also performed non-parametric multiple comparison tests on the RCAs and BCAs using Dunn's procedure~\cite{dunn1964multiple}, with a $p$-value correction using the False Discovery Rate method~\cite{benjamini1995controlling}. The results are shown in Table~\ref{tab:pvalue}, where the statistically significant ones are marked in bold.

Table~\ref{tab:pvalue} demonstrates that our proposed BN and UR can significantly improve the generalization performance of the traditional MBGD optimization for TSK fuzzy classifiers. \texttt{TSK-MBGD-UR-BN} statistically significantly outperformed \texttt{CART}, \texttt{JRip}, \texttt{PART}, \texttt{TSK-MBGD} and \texttt{TSK-MBGD-BN} on RCA, and also statistically significantly outperformed \texttt{CART}, \texttt{JRip}, \texttt{TSK-FCM-LSE}, \texttt{TSK-MBGD} and \texttt{TSK-MBGD-BN} on BCA. Although the performance improvement of \texttt{TSK-MBGD-UR-BN} over \texttt{RF} and \texttt{TSK-MBGD-UR} were not statistically significant, they were quite close to the threshold, especially for the BCA.

\begin{table*}[htpb]\centering
\caption{$p$-values of non-parametric multiple comparisons on the RCAs and BCAs.} \label{tab:pvalue}
\begin{tabular}{cccccccccc}\hline
                           & Metric & \texttt{CART}            & \texttt{RF}              & \texttt{JRip}            & \texttt{PART}            & \texttt{TSK-FCM-LSE}         & \texttt{TSK-MBGD}             & \texttt{TSK-MBGD-BN}          & \texttt{TSK-MBGD-UR} \\\hline
\multirow{2}{*}{\texttt{TSK-MBGD-BN}}    & RCA    & \textbf{0.0097} & 0.0628 & 0.1368          & 0.3239          & 0.2723          & \textbf{0.0000} & -               & -      \\
                           & BCA    & 0.1547          & 0.1627 & 0.4508          & 0.0981          & 0.4518          & \textbf{0.0000} & -               & -      \\\hline
\multirow{2}{*}{\texttt{TSK-MBGD-UR}}    & RCA    & \textbf{0.0001} & 0.3740 & \textbf{0.0090} & 0.0460          & 0.2452          & \textbf{0.0000} & 0.1146          & -      \\
                           & BCA    & \textbf{0.0036} & 0.2900 & 0.0912          & 0.4420          & 0.0921          & \textbf{0.0000} & 0.0731          & -      \\\hline
\multirow{2}{*}{\texttt{TSK-MBGD-UR-BN}} & RCA    & \textbf{0.0000} & 0.2113 & \textbf{0.0002} & \textbf{0.0025} & 0.0404          & \textbf{0.0000} & \textbf{0.0094} & 0.1409 \\
                           & BCA    & \textbf{0.0000} & 0.0291 & \textbf{0.0021} & 0.0730          & \textbf{0.0022} & \textbf{0.0000} & \textbf{0.0013} & 0.0986
\\\hline
\end{tabular}
\end{table*}

\subsection{Effect of UR}

As mentioned in Section~\ref{subsec:ur}, using MBGD to optimize the TSK fuzzy system may face the ``rich get richer" problem. To demonstrate this, Fig.~\ref{fig:urfiring} shows the average normalized firing levels of the rules on the entire dataset after the four MBGD-based TSK models were trained, on three representative datasets. For \texttt{TSK-MBGD}, a few ``richest" rules had much larger average firing levels than others, and hence the rules contributed significantly differently to the output. BN may help alleviate this problem a little bit, as the average normalized rule firing levels in \texttt{TSK-MBGD-BN} were more uniform than those in \texttt{TSK-MBGD}, which also resulted in better classification performances, as demonstrated in the previous subsection. However, UR had the most direct effect on alleviating the ``rich get richer" problem, as \texttt{TSK-MBGD-UR} (\texttt{TSK-MBGD-UR-BN}) had much more uniform average normalized rule firing levels than \texttt{TSK-MBGD} (\texttt{TSK-MBGD-BN}), and hence also better classification performance.

\begin{figure}[htpb]\centering
	\subfigure[]{\includegraphics[width=0.9\columnwidth,clip]{satellite_URfiring}}
	\subfigure[]{\includegraphics[width=0.9\columnwidth,clip]{vehicle_URfiring}}
	\subfigure[]{\includegraphics[width=0.9\columnwidth,clip]{biodeg_URfiring}}
	\caption{Average normalized rule firing levels of \texttt{TSK-MBGD}, \texttt{TSK-MBGD-BN}, \texttt{TSK-MBGD-UR} and \texttt{TSK-MBGD-UR-BN} on (a) Satellite, (b) Vehicle, and (c) Biodeg datasets. } 	\label{fig:urfiring}
\end{figure}

Note that we set $\tau=1/C$ in (\ref{eq:ur}), where $C=6$ for Satellite, $C=4$ for Vehicle, and $C=2$ for Biodeg. However, the actual average normalized rule firing levels were not exactly $\tau$ on these datasets. Our experiments showed that although UR cannot guarantee the average normalized rule firing levels to be around $\tau$, it can indeed make the rules fired more uniformly.

Why may making the rules fired more uniformly help improve the generalization performance? In \cite{wu2019functional} we pointed out that a TSK fuzzy system may be functionally equivalent to an adaptive stacking ensemble model, in which each rule can be viewed as a base learner, and the aggregation weights equal the corresponding rule firing levels. When the rule firing levels are more uniform, generally more rules are utilized in computing the output, i.e., more base learners are used in the stacking ensemble model, which may help improve the generalization performance.

To demonstrate this, we computed the entropy of the normalized rule firing levels for each input example:
\begin{align}
	E = -\sum_r^R \overline{f}_r~\log\overline{f}_r, 	\label{eq:entropy}
\end{align}
where $\overline{f}_r$ is the normalized firing level of the $r$-th rule. Generally, a larger entropy means more rules were fired.

Fig.~\ref{fig:entro} shows the histogram of the entropy distributions on the Satellite dataset. When training TSK fuzzy systems without UR, many samples had close to zero $E$, i.e., all except one rule had firing levels close to zero. When UR was added, the number of examples with close to zero $E$ decreased significantly, i.e., more rules with larger firing levels were used in computing the output.

\begin{figure}[htpb]\centering
	\subfigure[]{\includegraphics[width=0.95\columnwidth,clip]{TSKFiringEntropy}}
	\subfigure[]{\includegraphics[width=0.95\columnwidth,clip]{TSKBNFiringEntropy}}
	\caption{Histogram of the normalized rule firing level entropy $E$ of (a) \texttt{TSK-MBGD} and \texttt{TSK-MBGD-UR}, and, (b) \texttt{TSK-MBGD-BN} and \texttt{TSK-MBGD-UR-BN}, on the Satellite dataset.} 	\label{fig:entro}
\end{figure}

\subsection{Effect of BN}\label{subsec:bnres}

We also used the Satellite dataset to analyze the effect of BN.

We set the UR weight $\lambda=1$ and recorded the training loss and test BCA in the first 20 training epochs. This process was repeated 10 times, and the average results are shown in Figs.~\ref{subfig:histb} and \ref{subfig:hista}, respectively. BN resulted in smaller training losses and better generalization performances in testing.

There is still no agreement on theoretically why BN is helpful in optimizing deep neural networks~\cite{santurkar2018does}; thus, it is also challenging to analyze theoretically why BN can help the optimization of TSK fuzzy systems. Nevertheless, we performed an empirical study to peek into this, by recording the L1 norm of the antecedent parameters' gradients and the L1 norm of the consequent parameters' gradients in the first 20 training epochs on the Satellite dataset. The results are shown in Figs.~\ref{subfig:grada} and \ref{subfig:gradb}, respectively. BN significantly increased the gradients of both antecedent and consequent parameters. With the same learning rate, this can expedite the convergence.

\begin{figure}[htpb]\centering
	\subfigure[]{\includegraphics[width=0.9\columnwidth,clip]{bn_train_loss}\label{subfig:histb}}
	\subfigure[]{\includegraphics[width=0.9\columnwidth,clip]{bn_test_bca}\label{subfig:hista}}
	\subfigure[]{\includegraphics[width=0.9\columnwidth,clip]{ante_grad}\label{subfig:grada}}
	\subfigure[]{\includegraphics[width=0.9\columnwidth,clip]{cons_grad}\label{subfig:gradb}}
		\caption{(a) Training loss, (b) test BCA, (c) L1 norm of the antecedent parameters' gradients, and (d) L1 norm of the consequent parameters' gradients, in the first 20 training epochs on the Satellite dataset. The horizontal axis starts from 3 epochs so that the differences among the curves can be more clearly visualized.} 	\label{fig:bn_hist}
\end{figure}

We also evaluated the performances of the two BN variants introduced in Section~\ref{subsec:bn}. The BCAs of \texttt{TSK-MBGD-UR}, \texttt{TSK-MBGD-UR-BN}, \texttt{TSK-MBGD-UR-GBN} and \texttt{TSK-MBGD-UR-RBN} are shown in Table~\ref{tab:bcabnv}. \texttt{TSK-MBGD-UR-BN} performed the best, and \texttt{TSK-MBGD-UR-GBN} the worst. Since \texttt{TSK-MBGD-UR-RBN} had more parameters to optimize, its training was not as stable as \texttt{TSK-MBGD-UR-BN} and \texttt{TSK-MBGD-UR-GBN}. Therefore, \texttt{TSK-MBGD-UR-BN} is the best choice.

\begin{table}[htpb]\centering
\caption{Average BCAs of the three BN variants on the 12 datasets.} \label{tab:bcabnv}
\setlength{\tabcolsep}{1mm}
\begin{tabular}{ccccc}\hline
\multirow{2}{*}{Dataset}        & \texttt{TSK-MBGD}          & \texttt{TSK-MBGD}      & \texttt{TSK-MBGD}       & \texttt{TSK-MBGD} \\
& \texttt{-UR}          & \texttt{-UR-BN}      & \texttt{-UR-GBN}       & \texttt{-UR-RBN} \\\hline
Vehicle           & 0.7127          & \textbf{0.7930}  & 0.7261 & 0.7679              \\
Biodeg            & 0.8390          & 0.8439          & 0.8422 & \textbf{0.8440}      \\
DRD               & 0.6717          & \textbf{0.6729} & 0.6636 & 0.6650               \\
Yeast             & 0.4946          & 0.5332          & 0.4352 & \textbf{0.5339}      \\
Steel             & 0.7245          & \textbf{0.7515} & 0.7332  & 0.7219              \\
IS & 0.8559          & \textbf{0.9501} & 0.9115 & 0.8938               \\
Abalone           & 0.4791          & \textbf{0.5402} & 0.4924 & 0.5275               \\
Waveform21        & \textbf{0.8362} & 0.8233          & 0.8232 & 0.8334               \\
Page-blocks       & 0.6033          & \textbf{0.6710} & 0.5912  & 0.6333              \\
Satellite         & 0.8679          & \textbf{0.8700} & 0.8679  & 0.8216              \\
Clave             & 0.6374          & 0.6421          & 0.6090   & \textbf{0.6442}    \\
MAGIC             & 0.8225          & 0.7934    & \textbf{0.8319}       & 0.8318     \\\hline
Average & 0.7121          &  \textbf{0.7404}       & 0.7106  & 0.7265  \\\hline
\end{tabular}
\end{table}

\subsection{Effect of the Batch Size}

The batch size is an important hyper-parameter in MBGD-based optimization. It determines the memory requirement and the convergence speed in training. A larger batch size leads to faster convergence but also requires more memory. In \cite{keskar2016large}, the authors analyzed the effect of the batch size on the generalization performance. Their results showed that using a larger batch size causes degradation in the model generalization performance, because it tends to converge to a shaper minimum, which makes the model sensitive to noise. A similar finding was presented in \cite{masters2018revisiting} that a smaller batch size leads to more stable and reliable training. However, since we used the mean, standard deviation and mean firing level of each batch to compute the losses, too small batch size may also lead to poor performance.

We validated our model on the Satellite dataset with batch size varying from 16 to 2,048. The test RCAs and BCAs averaged over 30 runs are shown in Fig.~\ref{fig:batch}. The test performance decreased with too small or too large batch sizes. For \texttt{TSK-MBGD-UR-BN}, it seems that a batch size within [64, 256] is a good choice.

\begin{figure}[htpb]\centering
\includegraphics[width=\columnwidth, clip]{batch_size}
\caption{Average RCAs and BCAs of \texttt{TSK-MBGD-UR-BN} on the Satellite dataset, using different batch sizes.} \label{fig:batch}
\end{figure}

\section{Conclusions and Future Research}\label{sec:conclusion}

TSK fuzzy systems are powerful and frequently used machine learning models, for both regression and classification. However, they may not be easily applicable to large and/or high-dimensional datasets. Our very recent research \cite{drwuGD2019} proposed an MBGD-based efficient and effective training algorithm (MBGD-RDA) for TSK fuzzy systems for regression problems. This paper has proposed an MBGD-based algorithm, \texttt{TSK-MBGD-UR-BN}, to train TSK fuzzy systems for classification problems. It can deal with both small and big data with different dimensionalities, and may be the only algorithm that can train a TSK fuzzy classifier on big and high-dimensional datasets. \texttt{TSK-MBGD-UR-BN} integrates two novel techniques, which are also first proposed in this paper:
\begin{enumerate}
\item UR, which is a regularization term in the loss function to ensure that all rules are fired similarly on average, and hence to improve the generalization performance.
\item BN, which normalizes the inputs in computing the rule consequents to speedup the convergence and to improve the generalization.
\end{enumerate}

Experiments on 12 UCI datasets from various domains, with varying size and feature dimensionality, demonstrated that each of UR and BN has its own unique advantages, and integrating them can achieve the best classification performance. \texttt{TSK-MBGD-UR-BN}, together with MBGD-RDA proposed in \cite{drwuGD2019}, shall greatly promote the applications of TSK fuzzy systems in both classification and regression, especially for big data problems.

The proposed \texttt{TSK-MBGD-UR-BN} also has some limitations, which will be addressed in our future research. First, for very high dimensional data, fuzzy partitions of the input space become very complicated, and numeric underflow may happen when the product $t$-norm is used. Further research shall consider rules that automatically select the most relevant attributes as the antecedents. Second, we shall investigate how to improve the interpretability of data-driven TSK fuzzy systems. This is also partially linked to the first problem, as reducing the number of antecedents can improve the interpretability of the rules.


\section{Introduction} \label{sec:intro}

The existence of the continuous and supersonic solar wind predicted by \citet{Parker1958} is a
consequence of the supersonic expansion of the solar corona. Fast solar wind
streams originate from coronal holes, and as they interact with the slow solar wind in front of it,
Corotating Interaction Regions (CIRs) are formed. The shocks associated with CIRs can accelerate
energetic ions far beyond 1 AU, and these particles stream into the inner heliosphere
\citep{VanHollebeke1978, Mason1999}. Unlike in the case of Coronal Mass Ejections (CMEs), the
probability of CIR shocks forming within 1 AU is only $\sim$30\% \citep{Jian2006}, and whenever a
solar energetic particle (SEP) event is observed at 1 AU, particle flux enhancements associated with
these CIRs at larger heliocentric distances show a higher particle flux intensity \citep{Lario2000b}. The physical processes leading to the
acceleration of these solar energetic particles
are described in detail by \citet{Lee1982} and \citet{Desai2016}. The CIR associated energetic
particles in the inner heliosphere were examined using observations by various spacecraft from
1 AU to $\sim$9 AU \citep{VanHollebeke1978}. However, most of the vantage points for these
observations, including Ulysses were beyond 3 AU, except for two   observations of Pioneer 11 at
$\sim$1.5 AU, of which one of them was ambiguous \citep{VanHollebeke1978}, and a few  observations
in mid\textendash 1992 and late 1994 when Ulysses moved from $\sim$5 AU near the ecliptic to a
latitude of 80$^\circ$ S at $\sim$2.5 AU \citep{Richardson2004}. Apart from these isolated
observations, most of the data points were from 1 AU \citep{Burlaga1985, Jian2006} and
between 3 and 5
AU \citep{VanHollebeke1978,  Burlaga1984b, Lario2000b,Crooker1999}.  \citet{Morgan2010} while studying the radar absorption events
at Mars caused by
energetic particles impinging on the Martian ionosphere mentioned that the forward and reverse
shocks associated with a CIR evolve between 0.5 and 2 AU \citep{Crooker1999} and therefore the
orbit of Mars, between 1.4 and 1.7  AU, is a key region to observe the acceleration of these
energetic ions.\\

The Mars Atmosphere and Volatile EvolutioN (MAVEN) spacecraft was launched in November 2013 primarily to understand the atmospheric
escape from Mars \citep{Jakosky2015a}. MAVEN has a unique suite of instruments
to measure the solar wind velocity, solar wind particle density, interplanetary magnetic field
(IMF),
solar wind electron pitch angle distributions, solar EUV irradiance, and solar energetic electron
and ion fluxes along with Martian ionosphere-thermosphere measurements. The  SEP
flux measurements from MAVEN  provide a new vantage point near $\sim$1.5 AU to study the
acceleration of energetic particles  by CIR and CME driven shocks.\\

In this study, we compare the SEP fluxes during the CIR  event of June 2015 and CME events of March
2015 as observed near Earth (1 AU) and near Mars ($\sim$1.5 AU) to understand the acceleration of
the energetic particles beyond 1 AU, especially for the CIR driven event. The uniqueness of this
study is that for the first time, we have a vantage point to study the energetic particle fluxes
near 1.5 AU associated with CIRs. The CME events are shown to demonstrate the differences between
the SEP fluxes in CIRs and other impulsive events at these vantage points. The direct measurements
presented here, from these two key regions in space would be useful to understand the evolution of
energetic particles in CIR and the changes in their spectral characteristics, within short
heliocentric distances.\\

\section{Data}\label{sec:data}

The energetic particle measurements at 1 AU for the selected events are from the Electron Proton
Alpha Monitor (EPAM) sensor onboard Advanced Composition Explorer (ACE) satellite. The EPAM
provides particle fluxes in the energy range from 47 keV to 4.75 MeV, in 8 channels. We have used
ion fluxes measured in 4 channels in the energy ranges, 0.31-0.58 MeV, 0.58-1.05 MeV,  1.05-1.89
MeV, and 1.89-4.75 MeV by the Low Energy Magnetic Spectrometer of the EPAM sensor. These data are
obtained from the ACE data center (\url{http://www.srl.caltech.edu/ACE/ASC/level2/}). The solar wind
velocity,  the IMF near 1 AU and the Geostationary Operational Environmental Satellite (GOES)
energetic particle flux variation ($>$ 10 MeV and $>$ 30 MeV) are obtained from the NASA Space
Physics Data Facility (SPDF) OMNIWeb data center (\url{https://omniweb.gsfc.nasa.gov/}).\\

The solar wind velocity and IMF values at 1.5 AU are obtained from the Solar Wind Ion Analyzer
(SWIA) and Magnetometer (MAG) instruments aboard MAVEN spacecraft. The SWIA and MAG data are used to
compute the upstream solar wind parameters following the method by \citet{Halekas2016}. The SEP
fluxes are obtained from the Solar Energetic Particle (SEP) instrument aboard MAVEN. This instrument
consists of two identical sensors, SEP 1 and SEP 2, each consisting of a pair of double\textendash
ended solid\textendash state telescopes  to measure 20 keV\textendash 1 MeV electrons and 30
keV\textendash 6 MeV ions in four orthogonal view directions \citep{Larson2015}. The data used in
this study  are the ion data  in the form of energy fluxes measured by the SEP 1 sensor in the 1F
direction  that typically views the Parker spiral direction \citep{Larson2015}. The gaps in the SEP
data are when the instrument attenuators close at altitudes below 500 km \citep{Lee2017}. These
datasets from the MAVEN instruments are downloaded from the the Planetary Data System
(\url{https://pds.nasa.gov/}). The level 2, version 01, revision 01 (V01\_R01) data of SWIA, level
2, version 01, revision 01/02 (V01\_R01/R02) data of MAG, and level 2, version 04, revision 02
(V04\_R02) data of SEP are used for the analysis. The Solar Dynamics Observatory (SDO) images
are taken from \url{https://solarmonitor.org/}. The ENLIL simulations during June and
March 2015
are taken from ENLIL Solar Wind Prediction (\url{http://helioweather.net/},
\url{https://iswa.ccmc.gsfc.nasa.gov/}).\\

\begin{figure}[h]
\includegraphics[width=1.0\linewidth]{figure1_sep_ace-cir.pdf}\\
\caption{(a) The variation of IMF (left Y-axis) and solar wind velocity (right Y-axis) for
the
period 5\textemdash17 June 2015 observed by ACE at 1 AU. (b) The SEP particle fluxes at different
energy ranges as observed by the EPAM sensor onboard ACE satellite
for the period 5\textemdash17 June 2015.}
\end{figure}

\begin{figure}[h]
\includegraphics[width=1.0\linewidth]{figure2_june_maven_rev.pdf}\\
\caption{(a) The variation of IMF  (left Y-axis) and solar wind velocity
(right Y-axis) for the period 17\textemdash24 June 2015. (b)The energy spectrum of ions in the 30
KeV to $\sim$6 MeV range observed by the MAVEN SEP instrument for the
period 17\textemdash24 June 2015. The Y-axis shows the energy range and the color axis shows the
logarithmic particle flux. The initial period shown in white had very low fluxes ($<$10) whereas the
gaps on 24 June 2015 corresponds to absence of good quality data. (c) The integrated SEP particle
fluxes in different energy ranges (comparable to ACE measurements) estimated from the MAVEN SEP
data.}
\end{figure}

\begin{figure}[h]
\includegraphics[width=1.0\linewidth]{figure3_power_law_ace-maven.pdf}\\
\caption{ (a) The energy spectrum during peak of ion flux enhancement in the CIR event on 8 June
2015 detected by ACE at $\sim$1 AU. The  power laws fitted with the observations are also shown.
The observations by ACE are extrapolated to 30 MeV, to show that at both 10 MeV and at 30 MeV, the
expected flux levels are below the detection limits \citep{Rodriguez2014} of GOES.
(b) The energy spectrum obtained from MAVEN.}
\end{figure}

\begin{figure}[h]
\includegraphics[width=1.0\linewidth]{figure4_sep_aceMar2015.pdf}\\
\caption{Same as Figure 1, but for the period 5\textemdash24 March 2015.}
\end{figure}

\begin{figure}[h]
\includegraphics[width=1.0\linewidth]{figure5_march_maven_rev.pdf}\\
\caption{Same as Figure 2, but for the period 1\textemdash13 March 2015. }
\end{figure}

\section{Overview of events}\label{sec:event}

A  comprehensive overview of the solar cycle 24 space weather conditions observed during the first
$\sim$1.9 years of the MAVEN science mission had been provided by \citet{Lee2017}. The events
described in this study are selected from \citet{Lee2017}. During June\textendash July 2015, a
number of low\textendash latitude coronal holes appeared and generated recurrent high\textendash
speed streams. The SDO/Atmospheric Imaging Assembly (AIA) composite images of the Earth-facing disk
of the solar corona shows a persistent low-latitude coronal hole that remain more or less similar
even after a 27 day solar rotation period. The images of the coronal hole are given in
\citet{Lee2017}. The ENLIL predictions (see Appendix \ref{appendix}) reveal that the CIR
hits Earth on 7 June 2015. The ENLIL snapshots further confirm that the same coronal hole
with negative IMF polarity is rotated towards Mars and arrive at Mars on 22 June 2015.  Assuming to
a first approximation that the rotation period for the coronal hole source is about 27.3 days, with
Mars located at 1.56 AU and in solar conjunction with Earth the solar wind stream would arrive
$\sim$15.5 days later (half a Carrington rotation $+$ $\sim$2.5 days to propagate from 1 AU to 1.56
AU) after the observation at Earth on 7 June \citep{Lee2017}. The ENLIL simulations also
show that at STEREO location the CIR stream would hit on 19 June 2015 (see \ref{appendix}). These
simulations confirm that these streams are originated from the same coronal hole. Further, the study
by \citet{Fisk1980} reports that typically, the high-speed stream is relatively steady in time, the
CIR and the shocks associated with it are stationary in the frame corotating with the Sun.
Therefore, the energetic ions from a single CIR event can be  seen for a period of 17 days and
$\sim$225$^\circ$ in solar longitude. \citet{Reames1997} also find clear evidence of energetic ions
from CIR shocks over a vast spatial region corresponding to far more than half a solar rotation.\\

End of February 2015 through early March 2015 was a period of intense space weather disturbances
at Mars \citep{Lee2017}.
Bright CME loops were observed in the
southeast quadrant of the Solar
and Heliospheric Observatory Large Angle and Spectrometric Coronagraph (LASCO)
C2 field of view on 6 March 2015, and peak enhancements in solar wind speed, density, dynamic
pressure, and IMF were observed at Mars on 8 March
2015 (\citet{Thampi2018} and references therein). The solar source,
for these events, Active Region 12297 was  on the backside of the Sun (as viewed from Earth), which
rotated over the east solar limb and produced a
M8.9 class flare that was observed by both MAVEN and
GOES on 7 March 2015 \citep{Lee2017}. As this active region rotated toward the direction of Earth,
it produced a number of flares and CMEs including the 17 March 2015 superstorm at Earth, popularly known as the
``St. Patrick's Day Storm''.
ENLIL simulations confirm the CME arrivals and shows that these events at Earth and Mars
are
produced by different eruptions from the same active region (see \ref{appendix}). \\

\section{Observations}\label{sec:obs}

Figure 1a shows the variation of IMF and solar wind velocity as observed near 1 AU, during
5\textendash 17 June 2015. The CIR reached the vicinity of Earth on 7 June 2015, marked by an
enhancement in velocity. This enhancement in velocity continued to remain near the highest value
almost for 4 days and then remained at a slightly lower value (which is still higher than the
velocities observed during 5\textendash 7 June 2015) upto  17 June 2015. The Z\textendash component
of the IMF showed largest fluctuations on 7 and 8 June, followed by fluctuations of smaller
magnitude, a typical feature of a CIR driven storm \citep{Borovsky2006}. Figure 1b shows the
variation of the energetic particle fluxes of energy range from $\sim$300 keV to $\sim$4.7 MeV as
observed by the EPAM sensor onboard ACE satellite located near 1 AU. The signature of the CIR
arrival is seen as $\sim$2 orders of magnitude increase in the SEP fluxes. The arrival of the CIR is
preceded by a calm period with almost no SEP fluxes. This is another typical feature of
CIR events \citep{Borovsky2006}. The flux of particles in the energy range $\sim$0.3 MeV to
$\sim$0.5 MeV increased to $>$100 particles cm$^{-2}$s$^{-1}$sr$^{-1}$MeV$^{-1}$, whereas the
particle flux with energies $>=$1 MeV remained $\le$ 10 particles
cm$^{-2}$s$^{-1}$sr$^{-1}$MeV$^{-1}$.\\

Figure 2a shows the variation of IMF and solar wind velocity as observed
by the MAG and SWIA instruments aboard MAVEN located at 1.56 AU. The total magnetic
field increased to $\sim$9 nT coincident to the arrival of CIR. The solar wind velocity
increased from the quiescent value of $\sim$300 km s$^{-1}$ to $\sim$700 km s$^{-1}$ and remained to
be $>$ 600
km s$^{-1}$ for the next day as well and remained high upto the end of June
\citep{Lee2017}.\\

Figure 2b shows the differential energy flux spectrum of the solar energetic ions observed by the
MAVEN SEP instrument. During the CIR event (on 22 and 23 June 2015), there is an enhancement in the
flux of energetic particles. To see the actual number of accelerated energetic ions using the MAVEN
SEP instrument and compare them with the observations at 1 AU by ACE, we estimated the integrated
flux at energy bins comparable to that of EPAM/ACE observations (Figure 2c). The onset of disturbed
fluxes from 8 June at Earth shown in Figure 1b (coinciding with the enhancements in velocity and
magnetic field)  correspond to the enhancement seen from 22 to 23 June at Mars. The time differences
seen in the SEP enhancement at these locations is consistent with the estimated arrival time based
on the solar rotation period for the coronal hole source. It must be noted that the SEP fluxes
before 19 June at Mars had very low values and the source of the slightly enhanced SEPs seen on
19-20 June could be related to the streaming particles from field lines connected to a CME that
erupted over the west solar limb, although the CME has not directly hit Mars \citep{Lee2017}. The
SEP fluxes return to  lower values (below 100 particles cm$^{-2}$s$^{-1}$sr$^{-1}$MeV$^{-1}$) before
the arrival of CIR on 22 June. The arrival of the CIR is marked by the enhancements in solar wind
velocity and magnetic field, and coincident enhancement in SEP fluxes by more than 2 orders from
the pre-CIR levels.  In comparison to the observations
from 1 AU, fluxes enhanced in the similar energy ranges at 1.5 AU. The peak ion flux in the
$\sim$1\textendash 2 MeV bin enhanced by $\sim$300 times at 1.5 AU. This indicates that
acceleration took place beyond 1 AU, but within 1.5 AU. Similarly, $\sim$140 times enhancement is
observed for the flux in the $\sim$600 keV\textendash 1 MeV energy bin. Therefore, we can infer that
the acceleration by the shocks associated with CIR could cause an enhancement of around two orders
of magnitude in the SEP fluxes when observed near 1.5 AU.\\

Figure 3a shows the energy spectrum corresponding to the peak of ion flux enhancement for the CIR
event on 8 June 2015 detected by ACE. Double power laws are fitted, one with
E$^{-1.88}$ for flux below 0.2 MeV and the second one with E$^{-2.78}$ for the flux between 0.2 to
3 MeV. These spectral indices are similar to those reported by \citet{Bucik2009} for the May 2007
events. To see whether we can expect higher energy particles at the energy ranges of GOES,
the ACE spectrum is extrapolated to 30 MeV with the same power law. With this approximation
itself, it can be seen that at both 10 MeV and at 30 MeV the expected flux levels are below the
detection limits of GOES \citep{Rodriguez2014}. A double
power law or a power law with exponential rollover at a few to tens of MeV has been reported in many
events (e.g. \citet{Mason1999}).  If so,  the expected flux at the higher energy ranges will be
still lower than
what we estimated from a simple power law. Therefore, for this CIR event we do not expect that GOES
would detect the SEP fluxes, and the observations confirm that GOES did not see any particles above its detection threshold.\\

Figure 3b shows the  energy spectrum during peak of ion flux enhancement in the CIR event on 23
June 2015 detected by MAVEN, with power law fits. The first
portion (below 60 keV) may have contributions from the oxygen pickup ions \citep{Larson2015}. The flux in the energy range 0.5\textendash3 MeV
shows a spectral index of -2.82.   At both
these
locations, the observed energy spectra continued to extend as power laws right down to the lowest
energy measured by the instruments, which corroborates with the results of \citet{Mason1997} and
\citet{Chotoo2000}.\\

To demonstrate the difference between CIRs and other impulsive events, we
show in Figure 4a (Figure 5a) the variation of IMF and solar wind velocity
observed near 1 AU (1.5 AU), during 5\textendash 24 March 2015 (1\textemdash 13 March 2015). Figure
4a shows that the IMF showed large fluctuations
maximizing on 17 March. The solar wind velocity also showed concurrent
enhancement and fluctuations afterwards. Figure 4b shows the SEP fluxes observed by the  ACE satellite for the same period. The fluxes showed
enhancement in all
the energy levels during the CME event. Compared with the CIR related energetic particle fluxes,
these observations are more than 2 orders of magnitude higher in all the energy ranges observed by
EPAM/ACE. This is the typical characteristic of CME driven shock acceleration of energetic particles
and one of the major difference between CME and CIR driven storms at 1 AU. The GOES observations
(not illustrated) also indicate the presence of particles in the $>$ 10 MeV energy range, that is,
a peak flux of $\sim$7 cm$^{-2}$s$^{-1}$sr$^{-1}$ on 17 March.\\

Figure 5a shows the variation of IMF and solar wind velocity observed by  MAVEN for the period
1\textemdash 13 March 2015 and Figure 5b shows the corresponding SEP energy flux spectrum.  The
first  velocity enhancement occurred on 3 March and the highest velocity was observed on 8 March.
IMF variations were observed during both the events with  higher magnitudes during the second event.
The SEP spectrum showed flux enhancements with almost similar levels on both the occasions. MAVEN
observed the CME shock only after the arrival of the highest energy SEP ions and these high energy
($>$1 MeV) ions arrived first followed by the lower energy particles, during 7\textendash8 March.
Figure 5c shows the integrated particle fluxes at Mars in the energy ranges comparable to EPAM/ACE.
The flux levels showed increase at $\sim$1.5 AU at all the energy bins. In the energy range of
$\sim$500 keV to 1 MeV, the enhancement is only  $\sim$5 times. \\

\section{Discussion and Conclusions}
Most of the previous studies concerning CIR driven SEP flux enhancements were
confined to data from observations beyond 3 AU. The radial gradient of CIR associated energetic
particle streams was measured by instruments on multiple probes to the inner heliosphere
(Helios 1 and 2,  Mariner 10, Ulysses), and the outer solar system (Pioneer 10 and 11) with
reference spacecraft at 1 AU such as IMP 7 and 8, ACE, WIND, and STEREO, and most of these
measurements pertained to total ion flux in the 0.9\textendash 2.2 MeV energy range. The Ulysses observations
at $\sim$1.4 AU, during its rapid latitude scan  phase  had shown the presence of reverse waves
that had not yet steepened into shocks \citep{Gosling1995c}. Apart from this, there was one
unambiguous observation of CIR induced acceleration of ions  near 1.5 AU, by Pioneer 11, which was
total ion fluxes in the energy range 0.9\textendash 2.2 MeV \citep{VanHollebeke1978}. The
enhancement in the fluxes reported here is almost an order of magnitude higher than that reported
by \citet{VanHollebeke1978}. \\

The  observations presented here clearly illustrate the following aspects:\\

(a) During the CME events in March 2015, the SEP observations from both 1 and 1.5 AU show
enhancement of high energy particle fluxes. The observations suggests that shock acceleration
takes place at 1 AU itself during CME events \citep{Lario2000b}. The comparison facilitates to
demonstrate the differences between SEP acceleration in CIRs and other impulsive events like
CMEs.\\

(b) During the CIR event in June 2015, at 1 AU, GOES did not observe any particle flux above its
detection threshold, and the particle flux $>=$1 MeV detected by ACE remained $\le$ 10 particles
cm$^{-2}$s$^{-1}$sr$^{-1}$MeV$^{-1}$. Compared to this, the SEP instrument aboard MAVEN showed high
SEP fluxes upto $\sim$3 MeV. Comparison shows that the acceleration by the shocks associated with
CIR can produce an enhancement of around two orders of magnitude in the SEP fluxes when observed
near 1.5 AU. Interestingly, the peak high energy particle fluxes at Mars are of almost similar
levels during CME driven event indicating similar levels of acceleration during these events at 1.5
AU.  In contrast, these flux levels are drastically different in the case of ACE observations at 1
AU ($>$ one order of magnitude higher in the $\sim$1\textendash 2 MeV range during CME), consistent
with the expectation that strong shock acceleration takes place inside 1 AU during CME events
\citep{Li2005, Desai2008}.  The comparison of energy fluxes at different ranges may point toward the
change of energy with distance as particles travel several AUs in the heliosphere under
cross\textendash field diffusion \citep{Zhao2016}. These observations are unique because of the
addition of a vantage point  at 1.5 AU, and as such observations were previously very sparse between
1 and 3 AU.\\

\newpage


\section{Introduction}

Our capabilities of both perceiving and reacting to the affective behavior of other persons are fine-tuned based on the observed social response of our interaction peers. We usually perceive how others are behaving towards us by reading their affective behavior through the processing of audio/visual cues \cite{lang1997motivated}. We are then able to provide an affective response, either by modulating our voice, using facial expressions, producing body movements, or realizing a specific action within the interaction context. The response is usually modulated by a contagion effect \cite{hatfield1993emotional}, where we are impacted by and replicate the perceived affective state of others, or by reacting with a proper social response \cite{hao2017enhancing}, which is accepted and expected by the interaction peers. Deciding the proper response is a learning process modulated by the interpretation of the impact that our behavior causes within the interaction context \cite{ochsner2008cognitive}. The learning association between perception and action configures a long-term adaptation mechanism which updates how we perceive and respond to others over a lifelong span. Adapting such a mechanism into a robotic platform would allow it to provide socially acceptable responses when in interactive social scenarios, and would benefit greatly the acceptability and human-likeliness of companion and assistive robots \cite{crossman2018influence}.

Adapting this perception-action mechanism into a social robot demands the development and integration of two mechanisms: an adaptive emotion appraisal, which is updated over time towards novel persons, and a behavior generation which learns based on the perception of the social impact that a certain action caused within the interaction.

The state-of-the-art in emotion recognition for social robots was much benefited with the recent advents of deep neural networks \cite{chen2018softmax}. Such models can learn, implicitly, how to represent high-level stimuli, such as facial expressions and prosody in speech. Such capability allows these models to represent multimodal emotion expressions efficiently and robustly, and they currently post as the state-of-the-art solutions in different emotion recognition tasks. When deployed to real-world scenarios, in particular, involving social robots, however, these solutions tend to fail, mostly due to their inability to adapt quickly to how specific persons express emotions \cite{sariyanidi2014automatic}. Deep learning models have millions of different parameters which need to be updated to depict different affective characteristics. Given the subjectivity on perceiving and labeling emotion expressions, to update these models to learn the specific characteristics of how a person expresses emotions without a large number of training samples of that specific person posts as a machine learning challenge.

Generating behavior based on affective perception has been deeply studied by the social robots area \cite{hirokawa2018adaptive}. Most of the proposed models, however, employ an emotion contagion strategy which usually reproduces the perceived facial expressions or body movement \cite{van2018generic}, or by using a simple decision tree on how to behave when a certain emotion was perceived \cite{tuyen2018emotional}. Such solutions are ineffective when dealing with complex interaction scenarios, such as interactions which are based on empathetic responses, very common to conversational robots and robots are driven by social responses \cite{leite2013social}. As soon as the robot has to provide a response which was not scripted by pre-defined decision mechanisms, it will usually fail to maintain an human-level of engagement with the users \cite{sidner2004look}, which was observed to reduce its applicability when interacting with children or elderly citizens in therapy or health-care scenarios \cite{chandrasekaran2015human}. One of the most common machine learning solutions to generate robot behavior is to provide a learning mechanism based on reinforcement learning. Usually, in such learning scheme, a robot reads the feedback from the environment and use it in a reward function to modulate its learning. The advantage of using reinforcement learning is that the robot can learn by itself, with little or no human supervision, by measuring the effect that its actions caused in the environment. When deployed in Human-Robot Interaction (HRI) scenarios, however, reinforcement learning strategies are usually not effective as they demand several thousands of examples to learn associations between perception and actions. As in most social interaction scenarios the robot interacts with a human, thousands of interactions would be necessary before the robot learns meaningful responses \cite{churamani2018learning}, reducing its applicability in real-world scenarios.

This position paper focuses on a multi-subject interaction scenario where different humans play a competitive card game, known as UNO, with a robot. The robot will be an active part of the interaction and will play against humans. As such, the robot will perceive the individual and group affective behavior and the game status and derive an action which is composed of behavioral response (facial expression and body movement), and a game-related decision (e.g. play a card or collect a card).

\begin{figure}[t]
\begin{center}
   \includegraphics[width=0.55\linewidth]{images/icubRobot.png}
\end{center}
   \caption{The iCub humanoid robot. We will use it in our proposed scenario due to its social capabilities, both for the perception of audio/visual cues and the generation of facial expressions and body movement.}
\label{fig:icubRobot}
\end{figure}

We also describe in this paper the first steps towards a hybrid neural architecture to be embedded into an iCub robot\cite{metta2008icub}, illustrated in Figure \ref{fig:icubRobot}. We drive the development of the proposed framework based on two directions: an adaptive perception neural network which uses an unsupervised affective memory mechanism to overcome the problems of online learning in deep neural networks to adapt towards how specific persons express emotions, and to model an intrinsic affective state; and an interactive reinforcement learning strategy based on learning the social and contextual impact of a certain action by measuring the affective responses of the players, avoiding the necessity of active human interaction in the learning loop.

This paper details the road-map on how the individual models that compose the framework will be developed and integrated. We also describe, and justify, in details the evaluation scenario and protocol to obtain a set of objective evaluation measures to evaluate the impact of our proposed framework.

\section{Playing UNO with a Robot}

The iCub robot is one of the most advanced platforms for Human-Robot Interaction (HRI). It was evaluated in different interaction scenarios and it provides a high-level of engagement, acceptability, and reliability when interacting with different persons. These are important characteristics to improve the perception of the robot as an active social member in our scenario, which we assume it will be an important characteristic to elicit natural affective behavior.

In our proposed interaction scenario, the iCub is playing the card game named UNO with three other persons. In this game, each player receives a set of cards, with different colors, numbers, and possible actions. Starting from a pre-selected card, the players have to discard a card which matches the color or number of the pre-selected card. The game happens in turns. The first player which discard all the cards at hand wins the game. The action cards help to improve the competitiveness of the game by providing specific game-change events, such as ``make the next player collect 4 cards", or `` invert the playing order". In our scenario, the iCub knows the rules of the game and can do a set of pre-defined actions such as check its cards, discarding and collecting cards, calling ``UNO", and passing its turn. The robot can also display certain affective responses, using different body movements and facial expressions.

While playing, the robot will observe the individual affective behavior of the other humans and will decide which game-related action to take and which affective expression it will display. The UNO card game was chosen as it is a relatively controllable scenario, where each player has its turn to make an action, a highly competitive game in which trying to identify, or even manipulate, the other players' affective behavior is an important game strategy, and, most importantly, the game provides several turns of context-dependent interaction as illustrated in Figure \ref{fig:scenario}.

\begin{figure}[t]
\begin{center}
   \includegraphics[width=0.80\linewidth]{images/Scenario.png}
\end{center}
   \caption{Our competitive game scenario involving three persons and the iCub robot. The robot will read audio/visual cues of the persons and the game context. The iCub will produce an action composed of a game-related action, and affective behavior. The impact of the produced action on the game and the persons is read and used to identify if the action was proper, given the situation, and to modulate the learning of novel actions.}
\label{fig:scenario}
\end{figure}

The iCub must take into consideration how the three persons are behaving while playing to be able to react specifically and in an appropriate manner to each of them. The appropriate manner will be selected based on different desired targets: the iCub must win the game (maximizing the game strategy), it must improve the engagement of the other players (maximizing the generation of affective behavior), or it must provide a playing behavior which is accepted by the subjects as human-like (maximizing both affective behavior and game strategies).

As each person expresses emotions in a highly individual manner, the iCub must adapt its perception to each of them. But, as they pose as a group, it also must take into consideration how the affective behavior of each person impacts the reading of the emotional state of the entire group. The iCub must take into consideration the gradual evolution of the affective states of each player within the game, and associate it with the temporal change within the game context.

We derived this scenario from common real-world applications for companion and tutoring robots used in elderly care centers and interaction with children. Our scenario, however, increases the complexity of the interaction by providing an unrestricted turn-taking initiative, where the persons will be instructed to act as natural as possible.  To provide a real-time response with the iCub is of much importance, as it must be perceived as an active and autonomous player. Also, the interplay between the individual perception and the group affective behavior is very important for grounding the produced actions and for the developing of different game-related strategies.

\section{Towards a Cognitive Framework}

\subsection{Affective Modeling}

Recently, we proposed a novel hybrid neural framework to deal with the problem of adapting deep neural models towards specific characteristics of perceiving emotion expressions in an online manner \cite{barros2016developing}. A series of pre-trained convolution channels were used to represent facial expressions, body movement, and prosody on the speech from one person. The GWR network received the multimodal representation as input, and learning, online, how to cluster the sensory stimuli into different affective concepts. The GWR adapts to the perceived stimuli by creating or removing neurons and it can learn, quickly, how to model novel information. Each neuron represents a prototype emotional concept, and are used to represent affective expressions. We explored the use of individual GWRs to learn prototype emotional concepts from individual persons and used them as specific affective memories which represent how that specific person expresses emotions \cite{barros2017self}. The affective memories were evaluated on an emotion recognition task and were shown to be a state-of-the-art model on recognizing spontaneous emotion expressions.

To deal with the proposed scenario, we will take inspirations on the affective memory networks. We must devise a novel affective memory model which can deal with the temporal changes within the interaction, in particular, able to integrate the perceived stimuli asynchronously. The model will process facial expressions, body movement, and speech signals from all persons in the scenario,  which are all expressed within different contextual times. To deal with this difference we will investigate the use of gated recurrent units (GRUs) \cite{cho2014learning} which will contribute to a time-compressed representation of an expression. A general attention mechanism will also important to avoid that intense emotional changes are neutralized by long periods of displaying the same affective behavior. To address this problem we will integrate the concepts of local and global attention to the convolution channels of the model. Also, we will investigate the use of recurrent pooling layers will be used to highlight specific periods of the input stimuli which contribute more to emotion recognition.

To provide a general affective analysis of the situation, we will integrate the individual affective memories into one emotional read. We will investigate the use of recurrent gamma connections \cite{Parisi2018} within a GWR to model the impact that each person will have in the general emotional representation. The GWR, once again, will provide the model with an online association learning mechanism, which will be able to depict, in an online manner, how both persons are behaving based on their affective information.

 It is important to ground the proposed model on existing literature on emotion processing. This allows a direct comparison with humans not only on the performance level but also on behavioral aspects. Thus, the design of the novel learning rules, the topological and conceptual aspects of the model will be inspired by behavioral knowledge on how humans perceived emotion expressions.

\subsection{Modulating Game Strategies with Social Reactions}

To learn how to provide proper behavior, we will develop and evaluate a reinforcement learning model based on actor-critic neural networks. The actor-network will learn how to associate the current group perception state of the robot and the affective memory related to the person that the robot is currently aiming to play with, with a set of game and affective behavior actions.

By perceiving the reaction of the players to the selected action, and the change in the game context, the iCub will maximize the production of a proper reaction. The critic network will learn how to maximize the actor-network output based on the desired target: a reward function based only on maximizing the game score of the iCub will produce a set of behaviors which will drive the iCub to win the game; by providing a set of reaction that maintain the group affective state positive, will direct the iCub to produce behaviors which seem as enjoyable by the other players; and creating a reward function based on both improving the game score of the iCub and maintaining a positive affective behavior on the group, will develop a certain game strategy that we assume it will be perceived as more human-like.

To maximize the learning strategy based on the few data points collected in each interaction, we will investigate the use of predictive stimulation \cite{nagai2019predictive} to provide a continual learning mechanism for the actor-critic network. We will collect videos of persons playing the game, and use it to create a replay memory pool, which integrates perception-action-reaction associations.  Different modulation techniques, in particular, based on the interplay between exploration and exploitation will be devised to avoid that the model overfits towards the memory replay.

\section{Evaluation protocol}

We will evaluate the proposed model using two strategies: one offline, which will support the topological development of each model, and one with real persons, which will enforce the real-time processing properties demanded by our scenario.

The offline training will be important for the development of novel learning rules and hyperparameter optimization of the model. We will simulate the game using the strategies collected from the recorded data, and the iCub simulator with a set of pre-defined game strategies. We will evaluate the individual modules of the proposed model, as well as different integration strategies. For this evaluation, objective measures such as the accuracy on recognizing emotion expressions, the concordance between the model and external human observers and the processing time will be maximized.

The real-world evaluation will involve the described scenario. Ideally, we will have multiple sessions with the same participants happening over a month to capture the important characteristics of the affective behavior of each participant. We will shuffle each group of participants during the sessions to avoid perception bias associated always with the same group of participants. The expected behavior is that the robot learns how to play the game based on the three different goals, and devise specific playing characteristics for each person. After each interaction, the participants will fill out different HRI-related questionnaires. We are in particular interested in the Asch`s personality impression questionnaire \cite{asch1946forming}, which will give us the persons´ assessment of how the robot is behaving, and the Godspeed questionnaire \cite{bartneck2009measurement}, which will help us to understand the persons’ general opinion about the impact of the robot in the interaction. Combining these results with the objective measures of the offline evaluation will give us a general perspective on how the proposed model impacts the general acceptance of the robot in a competitive scenario.

\section{Conclusion}

The action-perception cycle of perceiving emotion expressions, and responding with socially accepted actions is one of the most important aspects of human-human natural communication. In this paper, we describe our intentions to develop a hybrid neural framework for learning proper affective responses with an iCub robot in a competitive game scenario. The framework will comprise of two modules, the first one implementing deep and unsupervised neural networks for learning how to represent the group affective behavior from the specific characteristics of individual persons. The second module uses an actor-critic reinforcement learning neural network to provide actions, based on the perceived affective behavior. The impact of a chosen actions, based on three different desired targets, will be used as learning modulation for the model. We also present a series of evaluation strategies on how to provide a balance between objective and subjective analysis to improve the robustness of the model and to assess its impact on our HRI scenario.

\section*{Acknowledgment}
The authors gratefully acknowledge partial support from the German Research Foundation DFG under project CML (TRR 169). Alessanda Sciutti was supported by a Starting Grant from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation
programme, G.A. No 804388, wHiSPER.

\bibliographystyle{plain}