# Correntropy-Based Logistic Regression
with Automatic Relevance Determination for
Robust Sparse Brain Activity Decoding

Yuanhao Li, Badong Chen,
Yuxi Shi, Natsue Yoshimura, and Yasuharu Koike This work was supported in part
by the Japan Society for the Promotion of Science (JSPS) KAKENHI under Grant
19H05728, in part by the Japan Science and Technology Agency (JST) Support for
the Pioneering Research Initiated by Next Generation (SPRING) under Grant
JPMJSP2106, and in part by the National Natural Science Foundation of China
under Grant U21A20485 and Grant 61976175. _(Corresponding author: Yuanhao
Li.)_ Yuanhao Li, Yuxi Shi, Natsue Yoshimura, and Yasuharu Koike are with the
Institute of Innovative Research, Tokyo Institute of Technology, Yokohama
226-8503, Japan. (correspondence e-mail: li.y.ay@m.titech.ac.jp)Badong Chen is
with the Institute of Artificial Intelligence and Robotics, Xi’an Jiaotong
University, Xi’an 710049, China.

###### Abstract

Recent studies have utilized sparse classifications to predict categorical
variables from high-dimensional brain activity signals to expose human’s
intentions and mental states, selecting the relevant features automatically in
the model training process. However, existing sparse classification models
will likely be prone to the performance degradation which is caused by noise
inherent in the brain recordings. To address this issue, we aim to propose a
new robust and sparse classification algorithm in this study. To this end, we
introduce the correntropy learning framework into the automatic relevance
determination based sparse classification model, proposing a new correntropy-
based robust sparse logistic regression algorithm. To demonstrate the superior
brain activity decoding performance of the proposed algorithm, we evaluate it
on a synthetic dataset, an electroencephalogram (EEG) dataset, and a
functional magnetic resonance imaging (fMRI) dataset. The extensive
experimental results confirm that not only the proposed method can achieve
higher classification accuracy in a noisy and high-dimensional classification
task, but also it would select those more informative features for the
decoding scenarios. Integrating the correntropy learning approach with the
automatic relevance determination technique will significantly improve the
robustness with respect to the noise, leading to more adequate robust sparse
brain decoding algorithm. It provides a more powerful approach in the real-
world brain activity decoding and the brain-computer interfaces.

###### Index Terms:

brain decoding, correntropy learning, automatic relevance determination,
electroencephalogram, functional magnetic resonance imaging.

##  I Introduction

Brain activity decoding technology has been increasingly driving the
developments of the brain-computer interface (BCI) systems, which aims to
reveal the human’s mental states [1, 2]. Non-invasive brain activity recording
approaches, such as multi-channel electroencephalogram (EEG) and functional
magnetic resonance imaging (fMRI), are prevalently employed to develop
advanced brain decoding paradigms, mainly due to their ease of operation and
surgery-free merit. Numerous brain decoding scenarios have been investigated
in the literature. For example, EEG-based motor imagery is a canonical task,
which aims to forecast the motion patterns imagined by the subjects [3]. The
P300 event-related potential, which could be detected from EEG signals, has
proved to be a promising way to control a speller BCI system [4]. Steady-state
visual evoked potentials were used to predict visual stimulus through
frequency domain analysis [5]. As for fMRI-based studies, the visual
perceptions of different edge orientations were successfully decoded in [6].
Further, a more challenging issue, constraint-free visual image
reconstruction, was solved in [7]. Moreover, memory detection has proved to be
practicable with fMRI signals [8].

To expose the information and patterns inherent in the brain recording
signals, in particular to predict categorical variables, such as whether the
subject is imaging left or right, one usually exploits multivariate
classification models [2, 9, 10], where the logistic regression has been
widely utilized in the literature [11]. Although the logistic regression model
performs well for a generic classification task, it is vulnerable to the
performance degradation caused by over-fitting, especially when the number of
features is larger than that of samples, i.e. high-dimensional problem. Brain
decoding tasks are considerably prone to high- dimensional problem which may
result from multiple channels and wide bandwidth of EEG signals, and thousands
of voxels in fMRI signals [12].

To address the over-fitting problem in brain decoding, many studies have
investigated various strategies. A primary method to avoid over-fitting is to
utilize $L_{2}$-norm regularization so that the model parameters are shrunk
towards zero and the variance of the classifier is decreased as a result [13].
Another approach is to select a subset of the features so that the dimension
would be reduced, which is commonly achieved by $L_{1}$-regularization to
realize a sparse model [14, 15]. The sparse model has been increasingly
attractive for brain decoding because the selected features could indicate the
spatio-temporal patterns relevant to the specific cognitive tasks [16, 17, 18,
19].

An alternative method to achieve a sparse classifier is to use a sparse prior
distribution on the model parameters and update the classification model from
Bayesian perspective. In fact, the $L_{1}$-regularization is equivalent to
employing a Laplacian prior distribution [14]. Compared to regularization, an
advantage of using a sparse prior distribution is that, one does not need to
adjust the regularization parameter manually [20]. Notably, the automatic
relevance determination (ARD) technique [21] was introduced into the logistic
regression classifier [16], which is a hierarchical sparse prior and has
proved to be more adequate than the Laplacian prior (equally
$L_{1}$-regularization) for feature selection [22]. The ARD-based sparse
logistic regression [16] (SLR) has been widely employed for brain decoding,
including EEG decoding [23, 24, 25, 26], fMRI decoding [7, 27, 28, 29], and
current source density analysis [30, 31, 32].

However, the decoding performance of SLR may be limited, since it assumes the
Bernoulli distribution on the data in binary classification, of which the
validity can be degenerated hugely in a noisy classification task [33]. It is
noted that non-invasive brain recordings are known to be prone to noises. For
example, blinks and eye movements deteriorate EEG signals [34], while fMRI
suffers both physiological noises and motions [35]. Even though many
preprocessing methods were developed to isolate noises, such as the blind
source separation by the independent component analysis [36], one can hardly
ensure that all noises are well removed by preprocessing.

In the present study, we aim to investigate a robust version of the SLR
algorithm, which can suppress the negative effects caused by noises inherent
in the brain measurement signals. To this end, we introduce the correntropy
learning framework into the ARD-based sparse model. Correntropy learning,
developed in information theoretic learning (ITL) [37], exhibits inspiring
robustness for various machine learning and signal processing tasks, including
regression [38, 39], classification [40, 41], principal component analysis
[42], and so on.

The remainder of this paper is organized as follows. Section II briefly
reviews the generic logistic regression algorithm, the hierarchical ARD prior
distribution, and the ARD-based SLR algorithm. In Section III, we introduce
the correntropy learning approach, and then propose a new robust sparse
classification algorithm by integrating the correntropy learning and the ARD
prior distribution. In Section IV, to demonstrate the robustness and feature
selection capability of the proposed algorithm, we present extensive
experimental results with a synthetic dataset, an EEG dataset, and an fMRI
dataset, respectively. In addition, we provide some discussions in Section V.
Finally, the paper is concluded in Section VI.

##  II Preliminaries

Before presenting the correntropy-based robust sparse classification for brain
activity decoding, we review in this section the logistic regression
classification model, the ARD hierarchical sparse prior, and the ARD-based SLR
algorithm.

###  II-A Logistic Regression

Logistic regression is a widely used probability-based classification model. A
linear discriminant function that recognizes two different classes can be
expressed by a weighted summation of each input feature

$f(\mathbf{x},\omega)=\sum_{d=1}^{D}{\omega_{d}x_{d}}+\omega_{0}$ (1)

where $\mathbf{x}=(x_{1},x_{2},\cdots,x_{D})\in\mathbb{R}^{D}$ denotes a
$D$-dimensional input sample to be classified, while
$\omega=(\omega_{0},\omega_{1},\omega_{2},\cdots,\omega_{D})$ represents the
logistic regression model parameter, with a bias term contained. In the
$\\{$0,1$\\}$-label context, logistic regression computes the probability that
an input sample belongs to class 1 through the _sigmoid_ function as

$y\triangleq P(t=1|\mathbf{x},\omega)=\frac{1}{1+\exp(-f(\mathbf{x},\omega))}$
(2)

where $t$ is the label. Assuming the Bernoulli distribution leads to the
opposite probability for class 0 as $P(t=0|\mathbf{x},\omega)=1-y$. Given a
finite dataset $\\{(\mathbf{x}_{n},t_{n})\\}_{n=1}^{N}$ the likelihood
function is

$\begin{split}&P(t_{1},t_{2},\cdots,t_{N}|\mathbf{x}_{1},\mathbf{x}_{2},\cdots,\mathbf{x}_{N},\omega)\;\\\
=&\prod_{n=1}^{N}{P(t_{n}|\mathbf{x}_{n},\omega)}=\prod_{n=1}^{N}{y_{n}^{t_{n}}(1-y_{n})^{1-t_{n}}}\;\\\
\end{split}$ (3)

in which $y_{n}$ represents the probability of the _n_ -th sample for class 1
by

$y_{n}\triangleq
P(t_{n}=1|\mathbf{x}_{n},\omega)=\frac{1}{1+\exp(-f(\mathbf{x}_{n},\omega))}$
(4)

Maximizing the likelihood probability in (LABEL:equ:likelihood_function) with
an observed dataset is equal to maximizing the logarithmic form, yielding

$\omega^{*}=\text{argmax}\sum_{n=1}^{N}{(t_{n}\log
y_{n}+(1-t_{n})\log(1-y_{n}))}$ (5)

in which $\omega^{*}$ indicates the optimal model parameter. After one trains
the logistic regression classifier, for a new testing sample $\mathbf{x}$, the
label is predicted as class 1 if $f(\mathbf{x},\omega^{*})>0$ (or equally
$P(t=1|\mathbf{x},\omega^{*})>0.5$), otherwise as class 0.

###  II-B Automatic Relevance Determination

The maximum likelihood estimation (MLE), as is mentioned above, only takes
into account the distribution for the data. On the other hand, one could
assume a prior distribution for model parameter moreover.

ARD is a sophisticated sparse prior with a hierarchical form, proved to
exhibit superior feature selection capability than the Laplacian prior [22].
ARD employs the anisotropic Gaussian distribution for each model parameter,
that each parameter has an individual variance

$P(\omega_{d}|\alpha_{d})=\mathcal{N}(0,\alpha_{d}^{-1})\qquad d=0,1,\cdots,D$
(6)

in which $\mathcal{N}$ is the univariate Gaussian distribution, and
$\alpha_{d}$ is called the relevance parameter, controlling the possible range
of the corresponding model parameter. Each hyper-parameter is further assumed
with the non-informative prior distribution

$P_{0}(\alpha_{d})=\alpha_{d}^{-1}\qquad d=0,1,\cdots,D$ (7)

in a full-Bayesian formulation.

###  II-C ARD-Based Sparse Logistic Regression

The ARD-based SLR algorithm estimates the model parameters and the relevance
parameters by a maximum a posteriori (MAP) method. The posterior distribution
for $\omega$ is calculated

$P(\omega|{T},{X})=\frac{\int
P({T}|\omega,{X})P(\omega|\alpha)P_{0}(\alpha)\text{d}\alpha}{\int
P({T}|\omega,{X})P(\omega|\alpha)P_{0}(\alpha)\text{d}\alpha\text{d}\omega}$
(8)

where ${T}$ and ${X}$ denote the labels and attributes of the observed
dataset, respectively, and $\alpha=(\alpha_{1},\alpha_{2},\cdots,\alpha_{D})$
is a collection of the relevance parameters. The integration in (8) is,
however, difficult to calculate analytically. Hence, variational inference
[43] was introduced for approximation, defining the following free energy on
$Q(\omega,\alpha)$ to approximate the true joint posterior distribution
$P(\omega,\alpha|{T},{X})$ by

$F(Q(\omega,\alpha))=\int
Q(\omega,\alpha)\log\frac{P({T},\omega,\alpha)}{Q(\omega,\alpha)}\text{d}\alpha\text{d}\omega$
(9)

where the dependency on ${X}$ is omitted for simplicity. The free energy
$F(Q(\omega,\alpha))$ would be maximized when and only when $Q(\omega,\alpha)$
is equal to $P(\omega,\alpha|{T},{X})$. To maximize $F(Q(\omega,\alpha))$, one
further assumes the conditional independence between $\omega$ and $\alpha$ as
$Q(\omega,\alpha)=Q(\omega)Q(\alpha)$, leading to

$F(Q(\omega)Q(\alpha))=\int
Q(\omega)Q(\alpha)\log\frac{P({T},\omega,\alpha)}{Q(\omega)Q(\alpha)}\text{d}\alpha\text{d}\omega$
(10)

Thus one can maximize the free energy alternately with respect to $Q(\omega)$
and $Q(\alpha)$, respectively, with the following updates

$\begin{split}&\omega\text{-step:}\quad\log Q(\omega)=\left<\log
P({T},\omega,\alpha)\right>_{Q(\alpha)}+const\;\\\
&\alpha\text{-step:}\quad\log Q(\alpha)=\left<\log
P({T},\omega,\alpha)\right>_{Q(\omega)}+const\;\\\ \end{split}$ (11)

in which $\left<x\right>_{Q(x)}$ denotes the expectation of $x$ with respect
to the probability distribution $Q(x)$.

Despite the exceptional capability for feature selection, SLR suffers the
significant performance degradation resulting from noises in practice, because
its foundation is the assumption of Bernoulli distribution for the likelihood,
which was proved to exhibit poor robustness with respect to noises [33].
Hence, the optimization process for $\omega$-step will be degenerated in a
noisy classification task. The deteriorated solution for $\omega$-step would
further destabilize the calculation for the relevance parameters in
$\alpha$-step, which provides imprecise $\alpha_{d}$ for the next
$\omega$-step in turn. Thus, both model parameter $\omega$ and relevance
parameter $\alpha$ are deviated from the desired values.

##  III Method

###  III-A Correntropy Learning

Correntropy, developed as a generalized correlation function for stochastic
process, has been further generalized to measure the similarity between two
arbitrary variables in a reproducing kernel Hilbert space [38]. Correntropy
has a close relation with the information potential of the quadratic Renyi’s
entropy [37], in which the Parzen’s window method is employed to estimate the
probability density of the data distribution. Concerning two arbitrary random
variables $A$ and $B$, the correntropy is defined

$\mathcal{V}(A,B)=E[k(A-B)]$ (12)

in which $E$ is the mathematical expectation while $k$ is a kernel function.
One usually uses the Gaussian kernel function which is denoted by
$k_{\sigma}(x)\triangleq\exp(-x^{2}/2\sigma^{2})$ with kernel bandwidth
$\sigma$. Thus, $N$ observations of the variables $A$ and $B$, in practice,
bring about the following empirical estimation of correntropy

$\begin{split}\hat{\mathcal{V}}(A,B)&=\frac{1}{N}\sum_{n=1}^{N}{k_{\sigma}(a_{n}-b_{n})}\;\\\
&=\frac{1}{N}\sum_{n=1}^{N}{\exp(-\frac{(a_{n}-b_{n})^{2}}{2\sigma^{2}})}\;\\\
\end{split}$ (13)

If one desires to maximize the similarity between two variables in a machine
learning task, such as between the predicted and the target output, one can
train a learning machine to maximize their correntropy, which is called as the
maximum correntropy criterion (MCC) [38]. Employing a correntropy-based
learning criterion exhibits many benefits. First, the value of correntropy is
mainly determined by the Gaussian kernel function $k_{\sigma}$ along $A=B$,
which means it is a local similarity measure and could alleviate the negative
effect of large deviation caused by noise. It was also proved to extract more
information from the dataset because it involves all the even-numbered
statistical moments. In addition, it has a close relation with the
$m$-estimation. One could find a comprehensive explanation of correntropy in
[38].

###  III-B Correntropy-Based Robust Sparse Classification

The MAP calculation in (8), when one omits the dependency on ${X}$, could be
rewritten by integration as

$P(\omega|{T})=\frac{P({T}|\omega)P_{0}(\omega)}{P({T})}\propto
P({T}|\omega)P_{0}(\omega)$ (14)

where $P({T})$ is a constant that is usually called evidence, and
$P_{0}(\omega)=\int P(\omega|\alpha)P_{0}(\alpha)\text{d}\alpha$ indicates the
prior distribution for $\omega$ by integrating out $\alpha$. Since the
logarithm is a monotonically increasing function, MAP is equal to

$\begin{split}\max\ &P(\omega|{T})\;\\\ \Leftrightarrow\max\ &\log
P(\omega|{T})\;\\\ \Leftrightarrow\max\ &\log P({T}|\omega)+\log
P_{0}(\omega)\;\\\ \end{split}$ (15)

The poor robustness of the original SLR algorithm results from the log
likelihood item $\log P({T}|\omega)$. Since likelihood measures the
probability that the prediction is equal to the target output, comparably, we
could employ the local correntropy instead to evaluate the similarity between
prediction and desired output, motivated by the admirable robustness of
correntropy learning criterion. Therefore, we use the correntropy similarity
between the prediction and target $\mathcal{V}({T},{Y})$, in which $Y$ is the
collection of the predicted probability $y$, to substitute the non-robust log
likelihood $\log P({T}|\omega)$. The superior robustness of correntropy
$\mathcal{V}({T},{Y})$, compared to log likelihood $\log P({T}|\omega)$, was
verified on ordinary logistic regression model and radial basis function
network with extensive experimental results [33, 40].

To coordinate the correntropy item $\mathcal{V}({T},{Y})$ to the Bayesian
derivation of SLR described in Section II-C, we could rewrite the free energy
maximization in (10) as

$\begin{split}\max&\int Q(\omega)Q(\alpha)\times\;\\\ &(\log
P({T},\omega,\alpha)-\log Q(\omega)-\log
Q(\alpha))\text{d}\alpha\text{d}\omega\;\\\ \end{split}$ (16)

where the log joint distribution $\log P({T},\omega,\alpha)$ is decomposed

$\log P({T},\omega,\alpha)=\log P({T}|\omega)+\log P(\omega|\alpha)+\log
P_{0}(\alpha)$ (17)

from which one can also find the non-robust $\log P({T}|\omega)$. Here we
propose the correntropy-based pseudo log joint distribution $\log
P_{c}({T},\omega,\alpha)$, which is defined by

$\log P_{c}({T},\omega,\alpha)\triangleq\mathcal{V}({T},{Y})+\log
P(\omega|\alpha)+\log P_{0}(\alpha)$ (18)

where $P_{c}({T},\omega,\alpha)$ is called pseudo joint distribution since one
finds its integration over all values is not equal to one.

By defining the novel pseudo joint distribution $P_{c}({T},\omega,\alpha)$, we
reformulate the free energy maximization as

$\begin{split}\max&\int Q(\omega)Q(\alpha)\times\;\\\ &(\log
P_{c}({T},\omega,\alpha)-\log Q(\omega)-\log
Q(\alpha))\text{d}\alpha\text{d}\omega\;\\\ \end{split}$ (19)

Similarly, one can acquire the following alternate optimization

$\begin{split}&\omega\text{-step:}\quad\log Q_{c}(\omega)=\left<\log
P_{c}({T},\omega,\alpha)\right>_{Q_{c}(\alpha)}+const\;\\\
&\alpha\text{-step:}\quad\log Q_{c}(\alpha)=\left<\log
P_{c}({T},\omega,\alpha)\right>_{Q_{c}(\omega)}+const\;\\\ \end{split}$ (20)

Computing the expectation and omitting the constant in $\omega$-step, one
obtains

$\log
Q_{c}(\omega)=\frac{1}{N}\sum_{n=1}^{N}{\exp(-\frac{(t_{n}-y_{n})^{2}}{2\sigma^{2}})}-\frac{1}{2}\omega^{t}\varLambda\omega$
(21)

where $\varLambda=\text{diag}(\alpha_{1},\cdots,\alpha_{D})$. One finds that,
however, $\omega$-step is analytically intractable since the likelihood and
the prior are not conjugate. We further employ the Laplacian approximation to
$\log Q_{c}(\omega)$, a quadratic approximation at the MAP estimate, denoted
by $\omega^{*}$, yielding

$\log Q_{c}(\omega)\approx\log
Q_{c}(\omega^{*})-\frac{1}{2}(\omega-\omega^{*})^{t}H_{c}(\omega^{*})(\omega-\omega^{*})$
(22)

where $H_{c}(\omega^{*})$ is the negative Hessian matrix of $\log
Q_{c}(\omega)$ at $\omega^{*}$. Thus, $Q_{c}(\omega)$ is approximated with a
Gaussian distribution $\mathcal{N}(\omega^{*},S_{c})$, where $S_{c}\triangleq
H_{c}(\omega^{*})^{-1}$. The gradient of $\log Q_{c}(\omega)$ with respect to
model parameter is given by

$\begin{split}&\frac{\partial\log Q_{c}(\omega)}{\partial\omega}\;\\\
&=\frac{1}{N\sigma^{2}}\sum_{n=1}^{N}{\exp(-\frac{e_{n}^{2}}{2\sigma^{2}})e_{n}y_{n}(1-y_{n})\mathbf{x}_{n}}-\varLambda\omega\;\\\
\end{split}$ (23)

where $e_{n}=t_{n}-y_{n}$ is the prediction error. The Hessian matrix of $\log
Q_{c}(\omega)$ is given in (28).

In $\alpha$-step, given $Q_{c}(\omega)=\mathcal{N}(\omega^{*},S_{c})$, one can
obtain

$\log
Q_{c}(\alpha)=-\frac{1}{2}\sum_{d=1}^{D}{(\alpha_{d}(\omega_{d}^{*2}+s_{d}^{2})+\log\alpha_{d})}$
(24)

where $\omega^{*}$ is the posterior mean acquired from the $\omega$-step, and
$s_{d}^{2}$ is the _d_ -th diagonal element in $S_{c}$. $Q_{c}(\alpha)$ could
be regarded to obey the following Gamma distribution

$Q_{c}(\alpha)=\prod_{d=1}^{D}{Q_{c}(\alpha_{d})}=\prod_{d=1}^{D}{\varGamma(\alpha_{d}^{*},\frac{1}{2})}$
(25)

in which $\varGamma(\alpha_{d}^{*},\frac{1}{2})$ is the Gamma distribution
with the degree of freedom being $\frac{1}{2}$ and the expectation being
$\alpha_{d}^{*}$ that is

$\alpha_{d}^{*}=\frac{1}{\omega_{d}^{*2}+s_{d}^{2}}$ (26)

The reformulated ARD-based sparse logistic regression with the correntropy
learning framework, proposed in this paper, is named as CSLR. We then consider
the optimization for $\omega$-step and $\alpha$-step, respectively. In
$\omega$-step, $\log Q_{c}(\omega)$ is equivalent to an $L_{2}$-regularized
MCC-based logistic regression. Although it is non-convex as a result of the
integration of sigmoid function and Gaussian kernel function, it is acceptable
to obtain a local optimum for $\omega$-step, because it has been proved that
any local optimums of regularized $m$-estimation are sufficiently close to the
global optimum [44] while correntropy is exactly a robust formulation of the
Welsch $m$-estimator [38]. To acquire a local optimum for $\omega$-step, one
could utilize the half-quadratic (HQ) technique [33]. We provide the
mathematical derivation of HQ optimization for $\omega$ in Appendix A.
Considering the update for $\alpha$, one can use the following rule to
accelerate the convergence

$\alpha_{d}^{*}=\frac{1-\alpha_{d}^{*}s_{d}^{2}}{\omega_{d}^{*2}}$ (27)

which is motivated by the effective number of parameters [21].

$\frac{\partial^{2}\log
Q_{c}(\omega)}{\partial\omega\partial\omega^{t}}=\frac{1}{N\sigma^{2}}\sum_{n=1}^{N}{\mathbf{x}_{n}^{t}\left\\{\exp(-\frac{e_{n}^{2}}{2\sigma^{2}})\left[(\frac{e_{n}^{2}}{\sigma^{2}}-1)y_{n}^{2}(1-y_{n})^{2}+e_{n}y_{n}(1-y_{n})(1-2y_{n})\right]\right\\}\mathbf{x}_{n}}-\varLambda$
(28)

CSLR executes $\omega$-step and $\alpha$-step alternately, updating the model
parameters and relevance parameters. During the model training, the relevance
parameters of the irrelevant features can diverge to infinity, that the
probability density of corresponding model parameters is distributed at zero,
thus pruning irrelevant features and obtaining a sparse classifier. In
practice, one could initialize each $\alpha_{d}$ as 1 and set an upper limit,
such as $10^{8}$. If $\alpha_{d}$ exceeds the upper limit, the corresponding
features will be pruned in the subsequent model training. The proposed CSLR
for robust sparse classification is summarized in Algorithm 1.

Algorithm 1 CSLR for Robust Sparse Classification

1:input: training samples $\\{(\mathbf{x}_{n},t_{n})\\}_{n=1}^{N}$; Gaussian
kernel bandwidth $\sigma$; threshold for relevance parameter $\alpha_{\max}$;

2:initialize: model parameters $\omega_{d}$ ($d=0,1,\cdots,D$); relevance
parameters $\alpha_{d}=1$ ($d=0,1,\cdots,D$);

3:output: model parameters $\omega_{d}$ ($d=0,1,\cdots,D$)

4:repeat

5: $\omega$-step: update $\omega_{d}$ according to Appendix A;

6: $\alpha$-step: update $\alpha_{d}$ according to Equation (27);

7: if $\alpha_{d}\geqslant\alpha_{\max}$ then

8: adjust the corresponding model parameters to zero and prune the
corresponding features from the samples in the following iterations

9: end if

10:until the parameter change is small enough or the number of iterations
exceeds a predetermined limit

##  IV Experiments

In this section, the proposed CSLR algorithm was evaluated with a synthetic
dataset, an EEG-based motor imagery dataset, and an fMRI-based visual
reconstruction dataset, respectively, and was compared to the baseline, the
original SLR algorithm. For the kernel-based CSLR algorithm, we utilized the
five-fold cross-validation method to select the optimal kernel bandwidth from
20 values: 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.2, 1.4, 1.6,
1.8, 2, 4, 7, 10, 30, and 100. The thresholds for the relevance parameters
were specified as $10^{8}$ for both CSLR and SLR algorithms. The attributes
were normalized such that each feature was of mean $0$ and variance $1$ before
utilizing the classification algorithms. The maximum number for iterations is
set as 300 for both CSLR and SLR.

Note that, the crucial purpose of this study is to demonstrate the superior
robustness of the correntropy learning method for sparse Bayesian problem,
especially with the hierarchical ARD prior, which has not been investigated in
the literature. Hence, the principal performance comparison shall be
considered with the original ARD-based sparse model, i.e. SLR. We provide a
detailed discussion for the relations between CSLR and other existing methods
in Section V-B.

###  IV-A Synthetic Dataset

####  IV-A1 Dataset Description

First, we considered a noisy and high-dimensional synthetic dataset, with
which CSLR and SLR algorithms were evaluated in respect of classification
accuracy and feature selection. This example was created with a similar method
as in [33]. We first randomly generated 300 i.i.d. 111independent and
identically distributed training samples and 300 i.i.d. testing samples, which
obeyed a 500-dimensional multivariate standard normal distribution. For the
label of each sample, we generated a sparse true solution. Specifically, the
true solution $\omega^{*}$ was a 500-dimensional vector, where only five
components were relevant to the label while the other 495 components were
equal to zero:

$\omega^{*}=[\overset{\text{500-dimensional}}{\overbrace{\omega^{*}_{1},\omega^{*}_{2},\omega^{*}_{3},\omega^{*}_{4},\omega^{*}_{5},\underset{\text{495
components}}{\underbrace{0,0,\cdots,0}}}}]$ (29)

where the non-zero components were separately subject to the univariate
standard normal distribution. For each sample, the label was assigned 1 if the
product between the corresponding attribute and $\omega^{*}$ was larger than
0, otherwise assigned 0. Thus, we would train the classifiers with 300
training samples by 500 features, and evaluate them on the other 300 testing
samples.

Figure 1: Two corruption models were employed for the synthetic dataset: (a)
sample contamination (b) arbitrary contamination. Figure 2: Classification
accuracy on the noisy and high-dimensional synthetic example under two
different contamination models: (a) sample contamination (b) arbitrary
contamination. The results are averaged across $100$ Monte-Carlo repetitions,
where the error bar indicates the corresponding standard deviation. Figure 3:
Number of selected features and F1-score for feature selection: (a) sample
contamination (b) arbitrary contamination. The error bar denotes the standard
deviation. The average results and the corresponding standard deviations are
obtained from $100$ Monte-Carlo repetitions.

Considering the contamination for this synthetic dataset, we utilized two
different corruption models according to [45], as shown in Fig. 1. The sample
contamination means undivided samples are corrupted, and the arbitrary
contamination means any arbitrary elements in the attribute matrix may be
corrupted. To contaminate the data, we randomly selected a certain ratio of
samples or elements and replaced their attributes with zero- mean Gaussian
distributed noises with the following different standard deviations: 0.1, 0.3,
0.7, 1.0, 2.0, and 3.0. To assess the robustness of different classifiers, as
was suggested in [46], one is supposed to only contaminate the training data.
Hence, only the 300 training samples were corrupted in this synthetic dataset.
For both sample and arbitrary contamination, the ratio of the corrupted
samples/elements was increased from 0 to 1.0 with a step 0.05.

####  IV-A2 Results

We evaluated CSLR and SLR algorithms with 100 Monte-Carlo repetitions on this
synthetic dataset. First, we considered the classification accuracy, which
means the proportion of the correctly classified samples in the testing
dataset. The average classification accuracy for sample contamination and
arbitrary contamination is illustrated in Fig. 2 (a) and (b), respectively. As
one could observe in Fig. 2, for both sample contamination and arbitrary
contamination, the proposed CSLR outperformed the original SLR algorithm
significantly when the training data suffered corruption under each noise
standard deviation.

In addition, we assessed the capability for feature selection. The feature
selection can be regarded as an unbalanced binary classification, in which
there were 5 relevant features and 495 irrelevant features. The sparse
classifiers would select features in the model training. Hence, we could
evaluate the validity of the selected features with the true
‘relevant’/‘irrelevant’ labels. We used a comprehensive performance indicator,
F1-score 222Denoting ‘relevant’ as positive and ‘irrelevant’ as negative,
F1-score is the harmonic mean of the precision (TP/(TP+FP)) and the recall
(TP/(TP+FN)), which is computed by
2$\times$precision$\times$recall/(precision+recall)., to evaluate this
unbalanced binary classification. Fig. 3 shows the number of selected features
and F1-score for each algorithm in sample and arbitrary contamination,
respectively. The noise standard deviation was set as 1.0. As is shown in Fig.
3, CSLR always selected fewer features than SLR for both sample and arbitrary
contamination. More importantly, CSLR achieved a higher F1-score for the
feature selection than SLR under most noise proportions.

###  IV-B EEG-Based Motor Imagery Dataset

In addition, we evaluated the proposed CSLR algorithm on an EEG-based motor
imagery decoding task with the galvanic vestibular stimulation (GVS). This
dataset was proposed in our previous work [26].

####  IV-B1 Dataset Description

Ten healthy subjects were involved in this experiment. Their brain activities
during the experiment were recorded by a 64-channel EEG recording system at
2048 Hz. A customized GVS instrument was used in parallel to induce the
sensory feedback. In the experiment, the subjects were required to keep their
eyes closed and to imagine the motions (forward/backward), with random voice-
based cues. After 3-second cue period, the GVS instrument started to stimulate
the subjects with four directions for 0.5 second: forward/backward/left/right.
Then, the subjects rested for 3 seconds with a beep cue. A schematic diagram
of one trial is illustrated in Fig. 4. Each subject participated in 6
sessions, where each session consisted of 60 trials. Thus, the data of each
subject contained 360 trials. For each subject, the directions of the GVS-
induced sensory feedback were identical with the imagined motion directions
for 180 trials, while were inconsistent for the other 180 trials. One is
supposed to predict whether the direction of motor imagery is consistent with
the GVS-induced sensory feedback: match or mismatch. One can find a
comprehensive description of this dataset in [26].

Figure 4: A schematic diagram of the EEG-based motor imagery experiment with
GVS-induced sensory feedback.

####  IV-B2 Decoding Paradigm

To achieve an application for real-time BCI system, the EEG data was used for
decoding analysis with a rather raw state in this dataset. The decoding
paradigm is listed in what follows:

  * 1.

The EEG signals were down-sampled to 512 Hz.

  * 2.

The 360 trials were reordered randomly with their labels (match/mismatch).

  * 3.

The reordered dataset was separated, the first 80$\%$ defined as training
while the remaining 20$\%$ as testing.

  * 4.

The absolute magnitude of the EEG recording during the GVS period was used as
feature for each 100 ms duration (0–0.1 s, 0.1–0.2 s, 0.2-0.3 s, 0.3–0.4 s,
0.4–0.5 s).

  * 5.

As a result, there were 288 training samples (80$\%$) and 72 testing samples
(20$\%$) with 3,264 features (64 channels $\times$ 51 samplings in 100 ms).
The classifiers were trained on the training set and were evaluated on the
testing set.

  * 6.

Steps 2.–5. were implemented with 20 repetitions. Then we obtained the results
for each subject.

####  IV-B3 Results

The average classification accuracy in each 100 ms decoding window is shown in
Fig. 5 for each subject from 20 repetitions. The whiskers represent the
corresponding standard deviations. We also examined whether there existed
statistically significant difference between the accuracy obtained by CSLR and
SLR, respectively, according to a paired _t_ -test with $p<0.01$. Among the 50
conditions in total (10 subjects $\times$ 5 decoding windows), the proposed
CSLR achieved statistically higher classification accuracy in 44 conditions.
Though we failed to find significant difference in the remaining 6 conditions,
CSLR realized higher average classification accuracy than SLR as well.
Additionally, the average classification accuracy across a total of 10
subjects in each decoding window and across all decoding windows is shown in
TABLE I. The higher accuracy under each decoding window is marked in bold. One
can observe that the proposed CSLR achieved higher average accuracy than SLR
under each decoding window and across all the decoding windows.

Figure 5: Classification accuracy on the EEG-based motor imagery dataset with GVS-induced sensory feedback. The results are averaged across 20 repetitions, where the whiskers denote the standard deviations. ‘$\ast$’ indicates statistically significant difference according to a paired _t_ -test ($p<0.01$). TABLE I: Average classification accuracy across ten subjects in each decoding window and across all decoding windows with different classification algorithms. Decoding Window (s) | 0-0.1 | 0.1-0.2 | 0.2-0.3 | 0.3-0.4 | 0.4-0.5 | average
---|---|---|---|---|---|---
SLR | 75.75$\pm$5.12 | 79.38$\pm$4.70 | 81.03$\pm$4.82 | 80.74$\pm$4.70 | 79.18$\pm$5.35 | 79.22$\pm$5.28
CSLR | 84.15$\pm$3.90 | 84.32$\pm$3.65 | 83.62$\pm$3.79 | 84.08$\pm$3.72 | 83.66$\pm$3.77 | 83.97$\pm$3.77
Figure 6: Spatial patterns obtained by SLR and CSLR in EEG-based motor imagery dataset with GVS-induced sensory feedback. The contribution of each EEG channel $W(\text{ch})$ is shown in (a), sorted in descending order according to the spatial contribution by SLR. The top 5 EEG channels with the maximal spatial contributions are plotted in black circles, while the other channels with the top 16 spatial contributions are plotted in gray circles for SLR in (b) and CSLR in (c), respectively. TABLE II: Average classification accuracy across ten subjects in each decoding window and across all decoding windows with the generic logistic regression algorithm and the top 5 EEG channels selected by different sparse classifiers. Decoding Window (s) | 0-0.1 | 0.1-0.2 | 0.2-0.3 | 0.3-0.4 | 0.4-0.5 | average
---|---|---|---|---|---|---
top 5 channels selected by SLR | 79.86$\pm$4.94 | 84.04$\pm$3.72 | 83.78$\pm$3.90 | 83.64$\pm$3.75 | 83.37$\pm$4.15 | 82.94$\pm$4.37
top 5 channels selected by CSLR | 81.10$\pm$4.07 | 84.24$\pm$3.82 | 83.84$\pm$3.89 | 83.64$\pm$3.64 | 83.69$\pm$3.91 | 83.30$\pm$4.04

We further studied the spatial patterns for the classification models acquired
by CSLR and SLR, respectively, calculating how much each channel contributed
to the whole classification model. The element of the trained model parameter
is denoted by $\omega_{\text{ch,temp}}$ which corresponds to the EEG channel
‘ch’ and the sampling time ‘temp’. The spatial contribution for the channel
‘ch’, denoted by $W(\text{ch})$, is calculated by the ratio between the
summation of the absolute values of model parameters, which correspond to the
current channel, and the whole classification model

$W(\text{ch})=\frac{\sum_{\text{temp}}{\|\omega_{\text{ch,temp}}\|}}{\sum_{\text{ch}}\sum_{\text{temp}}{\|\omega_{\text{ch,temp}}\|}}$
(30)

The spatial contribution for each EEG channel is illustrated in Fig. 6 (a) for
SLR and CSLR, respectively, averaged across a total of ten subjects and the
five decoding windows. Moreover, the top 5 EEG channels with the maximal
spatial contributions are shown in black circles while the other top 16 EEG
channels are presented in gray circles in Fig. 6 (b) and (c) for SLR and CSLR,
respectively. To demonstrate the superior capability of feature selection for
CSLR, we utilized the EEG data from the respective top 5 channels and
implemented the same decoding paradigm by the MLE-based generic logistic
regression model. The average classification accuracy for each decoding window
across all subjects is listed in TABLE II. One can observe that the generic
logistic regression model exhibited higher average accuracy with the top 5
channels that were selected by CSLR in four decoding windows and higher
average accuracy across all decoding windows.

###  IV-C fMRI-Based Visual Reconstruction Dataset

Furthermore, the proposed CSLR was assessed on an fMRI-based visual image
reconstruction task, which was originally proposed in [7]. This dataset is
public available at Brainliner, a neuroscience data sharing platform developed
by Advanced Telecommunications Research Institute International,
Japan.333brainliner.jp/data/brainliner/Visual$\\_$Image$\\_$Reconstruction

Figure 7: Diagrammatic results on the fMRI-based visual stimulus
reconstruction dataset: (a) The comparison between the original and the
reconstructed visual stimulus by CSLR and SLR, respectively, for three
different categories in the figure image session. A total of 40 images were
presented to the subject and then reconstructed by the fMRI signals for each
category. The bottom rows illustrate the average reconstructed visual images
for each kind of the presented figure images. (b) The spatial correlation
(upper) and mean squared error (bottom) between the original and the
reconstructed visual stimulus for each category, averaged across the
corresponding 40 stimulus blocks. The error bars indicate the standard
deviations. ‘$\ast$’ means statistically significant difference according to a
paired _t_ -test with $p<0.01$. Figure 8: Number of selected features for the
100 local decoders that were ultimately trained by 440 random images. ‘$\ast$’
means CSLR selected fewer features than SLR with statistically significant
difference according to a paired _t_ -test with $p<0.01$.

####  IV-C1 Dataset Description

The subject was watching visual images consisting of 10 $\times$ 10 square
patches. Each patch was either a homogeneous gray area or flickering at 6 Hz.
Brain activities of the subject were recorded simultaneously by fMRI signals.
The dataset consists of 2 sessions: one random image session and one figure
image session. In the random image session, the shown images were formed in
stochastic patterns. 440 different random images in total were presented to
the subject. Each stimulus block lasted 6 s, followed with a 6 s rest period.
In the figure image session, there were 3 types of figure images: geometric,
alphabet letter layout 1, and alphabet letter layout 2. Each type had 40
blocks totally. Each stimulus block lasted 12 s, followed by 12 s rest. For
geometric images, five shapes were presented 8 times. For alphabet letter
layout 1, five letters were presented 8 times. For alphabet letter layout 2,
ten letters were presented 4 times. The preprocessing for fMRI data was
identical to the original study. In the following analysis, the V1 and V2
regions were utilized to reconstruct the images, in which 1,698 voxels in
total were involved. More details of the experiment can be found in [7].

####  IV-C2 Decoding Paradigm

The random image session was used to train reconstruction models with
different classifiers while the figure image session was utilized to evaluate
the reconstruction performance. As in [7], we used a linear combination of
local image decoders to reconstruct the 10 $\times$ 10 images by a local image
basis of size 1 $\times$ 1\. Accordingly, for a 10 $\times$ 10 image, we would
train 100 individual binary classifiers for each pixel, which predicted if the
pixel was flickering or a gray area. Afterward, the predicted contrast values
of each pixel were combined to obtain the final reconstructed images with the
combination coefficients, which were acquired by 10-fold cross-validation in
the random image session. A total of 440 random images were separated equally
into nine training groups and one validation group. The local decoders with
100 classifiers were trained on the nine training groups. Then, we found the
optimal non-negative combination coefficients by minimizing the sum of the
square error between the true and the predicted validation group by the
trained local decoders. The eventual combination coefficients were averaged by
the cross validations. Then, we retrained the local decoders with all the 440
random images and integrated them with the combination coefficients to
reconstruct the visual stimulus for the figure image session.

####  IV-C3 Results

We assessed CSLR and SLR algorithms separately with the above-mentioned
decoding method. The reconstructed images of the figure image session are
illustrated in Fig. 7 (a) by each type, in comparison to the original visual
stimulus. Further, we quantified the spatial correlation and mean squared
error (mse) between the original and the reconstructed visual stimulus. The
average result across 40 images for each figure image category is shown in
Fig. 7 (b). For each category, the proposed CSLR achieved higher spatial
correlation and lower mse than SLR with statistically significant difference
according to a paired _t_ \- test with $p<0.01$.

In addition, we considered the feature selection in this fMRI dataset as well.
For the ultimate 100 local decoders which were trained by 440 random images to
reconstruct the figure image, we counted the number of non-zero model
parameters for each decoder (classifier), i.e. the number of selected features
in each classification model. The result is illustrated by boxplot in Fig. 8
for SLR and CSLR, respectively. One perceives that CSLR selected statistically
fewer features (voxels) in fMRI decoding than SLR under a paired _t_ -test
with $p<0.01$.

##  V Discussion

###  V-A Classification Accuracy and Feature Selection

First, we would like to give a more detailed discussion about the experimental
results. There are two critical aspects to view the classification model in a
high-dimensional problem, i.e. the classification accuracy and the feature
selection.

The proposed CSLR was demonstrated by the experimental consequences to achieve
higher classification accuracy in noisy and high-dimensional decoding tasks.
In the synthetic dataset, CSLR realized nearly identical classification
accuracy as SLR when there was no contamination in the data. After we added
artificial contamination in the attribute matrix, especially when the noise
standard deviation is larger than 0.4 for sample contamination and all the
noise standard deviations for arbitrary contamination, one sees in Fig. 2 that
significant performance degradation happened to SLR even though the noise
proportion is equal to 0.05. In contrast, the proposed CSLR realized much less
performance degradation than SLR when the dataset was corrupted. Overall, in
this noisy high-dimensional toy dataset, except for when the noise proportion
is equal to zero or close to 1.0, the proposed CSLR almost always showed
better results than SLR. Next, for the EEG data, CSLR achieved statistically
higher accuracy than SLR for 44 conditions among a total of 50 scenarios (10
subjects $\times$ 5 decoding windows), while also realized higher average
accuracy in the remaining 6 conditions though without significant difference
(Fig. 5). In summary, the average classification accuracy for all subjects and
decoding windows was improved by 4.75$\%$ (TABLE I). Finally, one can observe
from the fMRI-based visual reconstruction results in Fig. 7 that, the
reconstructed images by CSLR are commonly more legible and closer to the
original stimulus. Quantitatively, CSLR achieved higher spatial correlation
and lower mse than SLR with statistically significant difference.

On the other hand, CSLR can select a more informative set of features. In Fig.
3, the number of selected features by SLR is 12.66$\pm$2.25 without
contamination, while is 4.40$\pm$0.96 for CSLR, which is closer to the true
number of the five relevant features. Meanwhile, CSLR achieved considerably
higher F1-score in feature selection (0.708$\pm$0.101 vs 0.443$\pm$0.105) than
SLR. When the data was contaminated, the number of selected features and
F1-score were obviously affected for SLR, while CSLR effectively suppressed
the negative effects of corruption by comparison. For the EEG dataset, since
it could be difficult to subjectively judge which spatial pattern is more
informative as shown in Fig. 6, we used the top 5 EEG channels selected by SLR
and CSLR, respectively, and assessed them by a generic logistic regression
model. Even by the same classifier, the top 5 channels of CSLR led to higher
accuracy than SLR (TABLE II), indicating that the spatial pattern of CSLR
could be more informative for the decoding task. For the fMRI dataset, CSLR
selected much fewer voxels than SLR, see Fig. 8. Notably we only employed the
voxels in the V1 and V2 regions, which are the most informative areas for
visual reconstruction, according to [7]. Despite such feature selection in
advance, CSLR gave better visual reconstruction results through much fewer
voxels.

Note that, in the experiments of this study, we did not make any artificial
contamination on the brain measurement signals, which is considerably
different from the related literature such as [47, 48]. We only evaluated the
proposed CSLR under the original brain decoding scenarios, in which it
outperformed the existing non-robust method. On the one hand, this verifies
that many adverse noises could exist in the brain recording signals. On the
other hand, in turn, these results further demonstrated the legitimacy to
employ the correntropy learning approach in brain decoding tasks.

###  V-B Extension to Other Machine Learning Models

CSLR algorithm could be viewed as a preliminary proposal to integrate the
correntropy learning framework with the sparse ARD technique. More generally
speaking, it was proposed by combining a robust information theoretic learning
criterion and sparse Bayesian learning means. To the best of our knowledge,
this is the first time such a proposal has been investigated. The existing
robust sparse machine learning and signal processing algorithms motivated by
information theoretic learning mainly employed a sparse regularization item to
a robust information theoretic learning based objective function. For example,
a new robust version of the sparse representation classifier (SRC) for face
recognition was proposed by using the $L_{1}$-regularization for the
correntropy-based objective function [49, 50]. In [51], a robust sparse
subspace learning algorithm was proposed for unsupervised feature selection,
with a regularized correntropy-based model. Moreover, a robust sparse adaptive
filtering was proposed in [52] by regularizing the MCC with the correntropy
induced metric (CIM). Different from the existing studies, the proposed CSLR
realizes a sparse model from the perspective of Bayesian learning with the ARD
prior. A crucial advantage, as mentioned above, is one does not need to adjust
the hyper-parameter manually that controls the strength of regularization
term. Thus, we would like to argue that the derivation of CSLR could be
intuitively generalized to the other machine learning models. For instance,
the $L_{1}$-regularized MCC-based models can be optimized with the Bayesian
learning framework, since $L_{1}$-regularization is equal to utilizing a
Laplacian prior. More importantly, one could derive new robust sparse
algorithms for other machine learning problems by employing the correntropy
learning criterion and the sparse ARD prior, because the real-world data could
suffer noise and numerous irrelevant features frequently.

###  V-C Extension to Multi-Class Classification

In this study, we proposed a new robust sparse classification model by
integrating the correntropy learning method with the ARD-based sparse logistic
regression algorithm which focuses on the binary classification problems. The
proposed CSLR can be easily generalized to the multi-class classification
problems. If there are $C$ classes to be classified, each class will have an
individual linear discriminant function

$f_{c}(\mathbf{x},\omega^{c})=\sum_{d=1}^{D}{\omega_{d}^{c}x_{d}}+\omega_{0}^{c}\quad(c=1,\cdots,C)$
(31)

where $\omega^{c}$ denotes the model parameter for _c_ -th class. Thus we
could calculate the probability that the _n_ -th sample belongs to _c_ -th
class through the _softmax_ function as

$\begin{split}y_{n}^{c}\triangleq
P(t_{n}=&c)=\frac{\exp(f_{c}(\mathbf{x},\omega^{c}))}{\sum_{k=1}^{C}\exp(f_{k}(\mathbf{x},\omega^{k}))}\\\
&(c=1,\cdots,C)\\\ \end{split}$ (32)

In multi-class problem, one usually utilizes the _one-hot_ coding for label
expression, such as $t_{n}=(1,0,\cdots,0)$ if _n_ -th sample belongs to the
first class. Similarly, one obtains the prediction by
$y_{n}=(y_{n}^{1},\cdots,y_{n}^{C})$. Hence, the prediction error would be a
$C$-dimensional vector by subtraction $e_{n}=t_{n}-y_{n}\in\mathbb{R}^{C}$.
One can then derive the multi-class CSLR algorithm naturally with the multi-
dimensional errors, which would be investigated in detail in our future works.
Multi-class robust sparse classifiers could further improve the brain activity
decoding performance. For example, multiscale local image decoders were proved
to show better visual reconstruction results in [7] than only using the
decoder of 1 $\times$ 1 size, by combining 1 $\times$ 1, 1 $\times$ 2, 2
$\times$ 1, and 2 $\times$ 2 decoders. Since the decoders for the other scales
necessitate multi-class classification, we merely used the 1 $\times$ 1
decoder in the evaluation. For our future works, we will assess the multi-
class CSLR algorithm with the multiscale decoders.

##  VI Conclusion

In this paper, we proposed a novel robust and sparse logistic regression
algorithm, by designing a correntropy-based pseudo joint distribution with the
sparse automatic relevance determination technique, which is optimized by
alternate updates from Bayesian perspective. The correntropy-motivated robust
sparse logistic regression algorithm, named as CSLR, was evaluated on various
noisy and high-dimensional classification tasks. The experimental results
demonstrated that CSLR can realize better classification accuracy and feature
selection for brain activity decoding tasks. We will generalize the proposed
CSLR for the multi-class situations and other machine learning algorithms in
our future studies. In addition, we are working on deriving the correntropy-
based likelihood function to incorporate it into the Bayesian framework with
more solid theoretical guarantees.

##  Appendix A Model Parameter Optimization

We aim to address the optimization issue in $\omega$-step with the half-
quadratic technique. Denoting the parameter space as $\mathcal{W}$, we can
rewrite $\omega$-step as

$\omega^{*}=\mathop{arg\max}_{\omega\in\mathcal{W}}\frac{1}{N}\sum_{n=1}^{N}{\exp(-\frac{e_{n}^{2}}{2\sigma^{2}})}-\frac{1}{2}\omega^{t}\varLambda\omega$
(33)

We define a convex function $\varphi(v)=-v\log(-v)+v$, in which $v<0$. The
Gaussian kernel function can be expressed by

$\exp(-\frac{e_{n}^{2}}{2\sigma^{2}})=\mathop{\sup}_{v<0}\\{v\frac{e_{n}^{2}}{2\sigma^{2}}-\varphi(v)\\}$
(34)

in which the supremum is achieved when and only when
$v=-\exp(-e_{n}^{2}/2\sigma^{2})<0$. Thus, we could rewrite (33) as

$\omega^{*}=\mathop{arg\max}_{\omega\in\mathcal{W},v_{n}<0}\frac{1}{N}\sum_{n=1}^{N}{(v_{n}\frac{e_{n}^{2}}{2\sigma^{2}}-\varphi(v_{n}))}-\frac{1}{2}\omega^{t}\varLambda\omega$
(35)

by introducing auxiliary variables $\\{v_{n}\\}_{n=1}^{N}$. We then optimize
(35) by alternately updating $\omega$ and $v_{n}$. Given the current model
parameters, one acquires $\\{e_{n}\\}_{n=1}^{N}$ with the logistic regression
model. The auxiliary variables are updated by

$v_{n}=-\exp(-\frac{e_{n}^{2}}{2\sigma^{2}})\quad(n=1,\cdots,N)$ (36)

Then one could fix the auxiliary variables and optimize $\omega$ with

$\omega^{*}=\mathop{arg\max}_{\omega\in\mathcal{W}}\frac{1}{N}\sum_{n=1}^{N}{v_{n}\frac{e_{n}^{2}}{2\sigma^{2}}}-\frac{1}{2}\omega^{t}\varLambda\omega$
(37)

for which one could use gradient-based optimization. Although (37) remains a
non-convex problem, the iterative half-quadratic approach is profitable, since
we can easily guarantee the model parameters to converge to a local optimum,
as proved in what follows.

Denoting the objective function in (35) as $E(\omega,v_{n})$ we have

$E(\omega,v_{n})=\frac{1}{N}\sum_{n=1}^{N}{(v_{n}\frac{e_{n}^{2}}{2\sigma^{2}}-\varphi(v_{n}))}-\frac{1}{2}\omega^{t}\varLambda\omega$
(38)

By comparing the values of $E(\omega,v_{n})$ for the (_k_ -1)-th and _k_ -th
half-quadratic iterations, we could prove the convergence to a local optimum
as

$E(\omega^{k-1},v_{n}^{k-1})\leqslant E(\omega^{k-1},v_{n}^{k})\leqslant
E(\omega^{k},v_{n}^{k})$ (39)

in which the first inequality is guaranteed by (34)(36). We can also establish
the second inequality provided we have a larger objective function value for
(37) with fixing the auxiliaries $v_{n}^{k}$. Thus, for each half-quadratic
iteration, $E(\omega,v_{n})$ is guaranteed to increase, until the model
parameters reach a local optimum.

## References

  * [1] J. R. Wolpaw, N. Birbaumer, D. J. McFarland, G. Pfurtscheller, and T. M. Vaughan, “Brain–computer interfaces for communication and control,” _Clinical neurophysiology_ , vol. 113, no. 6, pp. 767–791, 2002.
  * [2] J.-D. Haynes and G. Rees, “Decoding mental states from brain activity in humans,” _Nature reviews neuroscience_ , vol. 7, no. 7, pp. 523–534, 2006\.
  * [3] N. Padfield, J. Zabalza, H. Zhao, V. Masero, and J. Ren, “Eeg-based brain-computer interfaces using motor-imagery: Techniques and challenges,” _Sensors_ , vol. 19, no. 6, p. 1423, 2019.
  * [4] D. J. Krusienski, E. W. Sellers, F. Cabestaing, S. Bayoudh, D. J. McFarland, T. M. Vaughan, and J. R. Wolpaw, “A comparison of classification techniques for the p300 speller,” _Journal of neural engineering_ , vol. 3, no. 4, p. 299, 2006.
  * [5] C. Jia, X. Gao, B. Hong, and S. Gao, “Frequency and phase mixed coding in ssvep-based brain–computer interface,” _IEEE Transactions on Biomedical Engineering_ , vol. 58, no. 1, pp. 200–206, 2010.
  * [6] Y. Kamitani and F. Tong, “Decoding the visual and subjective contents of the human brain,” _Nature neuroscience_ , vol. 8, no. 5, pp. 679–685, 2005.
  * [7] Y. Miyawaki, H. Uchida, O. Yamashita, M.-a. Sato, Y. Morito, H. C. Tanabe, N. Sadato, and Y. Kamitani, “Visual image reconstruction from human brain activity using a combination of multiscale local image decoders,” _Neuron_ , vol. 60, no. 5, pp. 915–929, 2008.
  * [8] J. Peth, T. Sommer, M. N. Hebart, G. Vossel, C. Büchel, and M. Gamer, “Memory detection using fmri—does the encoding context matter?” _NeuroImage_ , vol. 113, pp. 164–174, 2015.
  * [9] K. A. Norman, S. M. Polyn, G. J. Detre, and J. V. Haxby, “Beyond mind-reading: multi-voxel pattern analysis of fmri data,” _Trends in cognitive sciences_ , vol. 10, no. 9, pp. 424–430, 2006.
  * [10] F. Pereira, T. Mitchell, and M. Botvinick, “Machine learning classifiers and fmri: a tutorial overview,” _Neuroimage_ , vol. 45, no. 1, pp. S199–S209, 2009.
  * [11] S. Menard, _Applied logistic regression analysis_. Sage, 2002, no. 106.
  * [12] K.-R. Müller, M. Krauledat, G. Dornhege, G. Curio, and B. Blankertz, “Machine learning techniques for brain-computer interfaces,” _Biomed. Tech_ , vol. 49, no. 1, pp. 11–22, 2004.
  * [13] F. Salehi, E. Abbasi, and B. Hassibi, “The impact of regularization on high-dimensional logistic regression,” _Advances in Neural Information Processing Systems_ , vol. 32, 2019.
  * [14] R. Tibshirani, “Regression shrinkage and selection via the lasso,” _Journal of the Royal Statistical Society: Series B (Methodological)_ , vol. 58, no. 1, pp. 267–288, 1996.
  * [15] A. Y. Ng, “Feature selection, l 1 vs. l 2 regularization, and rotational invariance,” in _Proceedings of the twenty-first international conference on Machine learning_ , 2004, p. 78.
  * [16] O. Yamashita, M.-a. Sato, T. Yoshioka, F. Tong, and Y. Kamitani, “Sparse estimation automatically selects voxels relevant for the decoding of fmri activity patterns,” _NeuroImage_ , vol. 42, no. 4, pp. 1414–1429, 2008.
  * [17] S. Ryali, K. Supekar, D. A. Abrams, and V. Menon, “Sparse logistic regression for whole-brain classification of fmri data,” _NeuroImage_ , vol. 51, no. 2, pp. 752–764, 2010.
  * [18] E. Satake, K. Majima, S. C. Aoki, and Y. Kamitani, “Sparse ordinal logistic regression and its application to brain decoding,” _Frontiers in neuroinformatics_ , vol. 12, p. 51, 2018.
  * [19] G. Ganesh, E. Burdet, M. Haruno, and M. Kawato, “Sparse linear regression for reconstructing muscle activity from human cortical fmri,” _Neuroimage_ , vol. 42, no. 4, pp. 1463–1472, 2008.
  * [20] M. A. Figueiredo, “Adaptive sparseness for supervised learning,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 25, no. 9, pp. 1150–1159, 2003.
  * [21] D. J. MacKay, “Bayesian interpolation,” _Neural computation_ , vol. 4, no. 3, pp. 415–447, 1992.
  * [22] D. Wipf and S. Nagarajan, “A new view of automatic relevance determination,” _Advances in neural information processing systems_ , vol. 20, 2007.
  * [23] G. Lisi, T. Noda, and J. Morimoto, “Decoding the erd/ers: influence of afferent input induced by a leg assistive robot,” _Frontiers in systems neuroscience_ , vol. 8, p. 85, 2014.
  * [24] G. Lisi and J. Morimoto, “Eeg single-trial detection of gait speed changes during treadmill walk,” _PloS one_ , vol. 10, no. 5, p. e0125479, 2015.
  * [25] G. Ganesh, K. Nakamura, S. Saetia, A. M. Tobar, E. Yoshida, H. Ando, N. Yoshimura, and Y. Koike, “Utilizing sensory prediction errors for movement intention decoding: a new methodology,” _Science Advances_ , vol. 4, no. 5, p. eaaq0183, 2018.
  * [26] Y. Shi, G. Ganesh, H. Ando, Y. Koike, E. Yoshida, and N. Yoshimura, “Galvanic vestibular stimulation-based prediction error decoding and channel optimization,” _International Journal of Neural Systems_ , vol. 31, no. 11, p. 2150034, 2021.
  * [27] K. Shibata, T. Watanabe, Y. Sasaki, and M. Kawato, “Perceptual learning incepted by decoded fmri neurofeedback without stimulus presentation,” _science_ , vol. 334, no. 6061, pp. 1413–1415, 2011.
  * [28] N. Yahata, J. Morimoto, R. Hashimoto, G. Lisi, K. Shibata, Y. Kawakubo, H. Kuwabara, M. Kuroda, T. Yamada, F. Megumi _et al._ , “A small number of abnormal brain connections predicts adult autism spectrum disorder,” _Nature communications_ , vol. 7, no. 1, pp. 1–12, 2016.
  * [29] T. Horikawa and Y. Kamitani, “Generic decoding of seen and imagined objects using hierarchical visual features,” _Nature communications_ , vol. 8, no. 1, pp. 1–15, 2017.
  * [30] H. Morioka, A. Kanemura, S. Morimoto, T. Yoshioka, S. Oba, M. Kawanabe, and S. Ishii, “Decoding spatial attention by using cortical currents estimated from electroencephalography with near-infrared spectroscopy prior information,” _Neuroimage_ , vol. 90, pp. 128–139, 2014.
  * [31] N. Yoshimura, A. Nishimoto, A. N. Belkacem, D. Shin, H. Kambara, T. Hanakawa, and Y. Koike, “Decoding of covert vowel articulation using electroencephalography cortical currents,” _Frontiers in neuroscience_ , vol. 10, p. 175, 2016.
  * [32] A. Mejia Tobar, R. Hyoudou, K. Kita, T. Nakamura, H. Kambara, Y. Ogata, T. Hanakawa, Y. Koike, and N. Yoshimura, “Decoding of ankle flexion and extension from cortical current sources estimated from non-invasive brain activity recording methods,” _Frontiers in neuroscience_ , vol. 11, p. 733, 2018.
  * [33] Y. Li, B. Chen, N. Yoshimura, and Y. Koike, “Restricted minimum error entropy criterion for robust classification,” _IEEE Transactions on Neural Networks and Learning Systems_ , 2021.
  * [34] T. Ball, M. Kern, I. Mutschler, A. Aertsen, and A. Schulze-Bonhage, “Signal quality of simultaneously recorded invasive and non-invasive eeg,” _Neuroimage_ , vol. 46, no. 3, pp. 708–716, 2009.
  * [35] T. T. Liu, “Noise contributions to the fmri signal: An overview,” _NeuroImage_ , vol. 143, pp. 141–151, 2016.
  * [36] T.-P. Jung, S. Makeig, C. Humphries, T.-W. Lee, M. J. Mckeown, V. Iragui, and T. J. Sejnowski, “Removing electroencephalographic artifacts by blind source separation,” _Psychophysiology_ , vol. 37, no. 2, pp. 163–178, 2000.
  * [37] J. C. Principe, _Information theoretic learning: Renyi’s entropy and kernel perspectives_. Springer Science & Business Media, 2010.
  * [38] W. Liu, P. P. Pokharel, and J. C. Principe, “Correntropy: Properties and applications in non-gaussian signal processing,” _IEEE Transactions on signal processing_ , vol. 55, no. 11, pp. 5286–5298, 2007.
  * [39] Y. Feng, X. Huang, L. Shi, Y. Yang, J. A. Suykens _et al._ , “Learning with the maximum correntropy criterion induced losses for regression.” _J. Mach. Learn. Res._ , vol. 16, no. 30, pp. 993–1034, 2015.
  * [40] A. Singh, R. Pokharel, and J. Principe, “The c-loss function for pattern classification,” _Pattern Recognition_ , vol. 47, no. 1, pp. 441–453, 2014\.
  * [41] G. Xu, B.-G. Hu, and J. C. Principe, “Robust c-loss kernel classifiers,” _IEEE transactions on neural networks and learning systems_ , vol. 29, no. 3, pp. 510–522, 2016.
  * [42] R. He, B.-G. Hu, W.-S. Zheng, and X.-W. Kong, “Robust principal component analysis based on maximum correntropy criterion,” _IEEE Transactions on Image Processing_ , vol. 20, no. 6, pp. 1485–1494, 2011.
  * [43] D. M. Blei, A. Kucukelbir, and J. D. McAuliffe, “Variational inference: A review for statisticians,” _Journal of the American statistical Association_ , vol. 112, no. 518, pp. 859–877, 2017.
  * [44] P.-L. Loh and M. J. Wainwright, “Regularized m-estimators with nonconvexity: Statistical and algorithmic theory for local optima,” _Advances in Neural Information Processing Systems_ , vol. 26, 2013.
  * [45] Y. Chen, C. Caramanis, and S. Mannor, “Robust sparse regression under adversarial corruption,” in _International Conference on Machine Learning_. PMLR, 2013, pp. 774–782.
  * [46] X. Zhu and X. Wu, “Class noise vs. attribute noise: A quantitative study,” _Artificial intelligence review_ , vol. 22, no. 3, pp. 177–210, 2004.
  * [47] H. Wang, Q. Tang, and W. Zheng, “L1-norm-based common spatial patterns,” _IEEE Transactions on Biomedical Engineering_ , vol. 59, no. 3, pp. 653–662, 2011.
  * [48] B. Chen, Y. Li, J. Dong, N. Lu, and J. Qin, “Common spatial patterns based on the quantized minimum error entropy criterion,” _IEEE Transactions on Systems, Man, and Cybernetics: Systems_ , vol. 50, no. 11, pp. 4557–4568, 2018\.
  * [49] R. He, W.-S. Zheng, and B.-G. Hu, “Maximum correntropy criterion for robust face recognition,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 33, no. 8, pp. 1561–1576, 2010.
  * [50] R. He, W.-S. Zheng, B.-G. Hu, and X.-W. Kong, “A regularized correntropy framework for robust pattern recognition,” _Neural computation_ , vol. 23, no. 8, pp. 2074–2100, 2011.
  * [51] N. Zhou, Y. Xu, H. Cheng, Z. Yuan, and B. Chen, “Maximum correntropy criterion-based sparse subspace learning for unsupervised feature selection,” _IEEE Transactions on Circuits and Systems for Video Technology_ , vol. 29, no. 2, pp. 404–417, 2017.
  * [52] W. Ma, H. Qu, G. Gui, L. Xu, J. Zhao, and B. Chen, “Maximum correntropy criterion based sparse adaptive filtering algorithms for robust channel estimation under non-gaussian environments,” _Journal of the Franklin Institute_ , vol. 352, no. 7, pp. 2708–2727, 2015.