#  EEG Signal Processing using Wavelets for Accurate Seizure Detection through
Cost Sensitive Data Mining

Paul G. Grant<EMAIL_ADDRESS>
Md Zahidul Islam<EMAIL_ADDRESS>
School of Computing and Mathematics
Charles Sturt University
Bathurst, NSW 2795, Australia

###### Abstract

Epilepsy is one of the most common and yet diverse set of chronic neurological
disorders. This excessive or synchronous neuronal activity is termed
“seizure”. Electroencephalogram signal processing plays a significant role in
detection and prediction of epileptic seizures. In this paper we introduce an
approach that relies upon the properties of wavelets for seizure detection. We
utilise the Maximum Overlap Discrete Wavelet Transform which enables us to
reduce signal noise. Then from the variance exhibited in wavelet coefficients
we develop connectivity and communication efficiency between the electrodes as
these properties differ significantly during a seizure period in comparison to
a non-seizure period. We use basic statistical parameters derived from the
reconstructed noise reduced signal, electrode connectivity and the efficiency
of information transfer to build the attribute space. We have utilised data
that are publicly available to test our method that is found to be
significantly better than some existing approaches.

Keywords: EEG, Graph Theory, MODWT, Wavelet variance, Decision trees

##  1 Introduction

Epilepsy is described by recurring seizures caused by abnormal activity in the
brain. This maybe considered a chronic disorder of the central nervous system
that exposes individuals to experiencing recurrent seizures. A seizure is a
transient irregularity in the brain’s electrical processes that may also
produce disruptive physical symptoms (Zandi et al., 2010).
Electroencephalogram (EEG) is a signal which represents the electrical
activity of neurons within the brain. The signal is acquired from the surface
of the scalp and these signals are non-stationary, that is its statistical
properties vary over time. The important frequencies from the physiological
viewpoint lie in the range of 0.1 to 30 Hz. (Shrestha et al., 2019). A common
classification problem is seizure detection, where non-seizure and seizure EEG
records of patients need to be identified.

It is not an onerous task to find a number of papers and differing methods
that seek to determine the onset of seizure. The nearest-neighbour classifier
was used on EEG features extracted in both time and frequency domains to
detect the onset of epileptic seizures (Qu and Gotman, 1997). Also a method
using various derived statistical parameters as extracted features has been
proposed (Siddiqui et al., 2018). One problem with these methods is that all
the information in the signal is being considered. Some of that information
may simply be noise and may also contain specific frequency bands that are not
related to the event that is to be classified. In another study the
application of a method based on energy extraction at various frequency sub-
bands using the wavelet transform was implemented (Elsa Jacob et al., 2018).
This avoids processing redundant data by selecting the required sub-bands /
wavelet decomposition levels which contain relevant information. Statistical
parameters of the energy at each level were derived and used as attributes for
a Support Vector Machine (SVM) classifier. Using this method with statistical
parameters of energy, at different frequency levels, and in this case with 21
electrodes (arranged in the international “10-20” format, as an example see
Figure 1.), the amount of computation required increases quickly if we require
additional frequency bands to better segregate the signal or additional
statistical parameters to fully describe the energy feature. This results from
the need to compute the statistical parameters for the energy within each
additional sub-band/wavelet decomposition level required and of course for the
signal from each electrode.

Figure 1:  10 - 20 Electrode positioning; developed from Helwig (2018).

The EEG signal contains several spectral components. Wavelet analysis allows
frequency separation of such components (Percival and Walden, 2000). Use of
Wavelets in EEG analysis is quite often seen and is well described, with the
Discrete Wavelet transform (DWT) frequently encountered (Faust et al., 2015).
(For brief overview of DWT see Appendix A.3.1). The wavelet transform is
effective for representing different aspects of non-stationary signals such as
trends, discontinuities, and where signal patterns are recurring, here other
signal processing approaches are unlikely to be as effective (Daubechies,
1990). Some limitations of the DWT are the requirement of the signal to be of
length $2^{J}$ where $J\in\mathbb{Z}^{+}$ and is the upper level of wavelet
decomposition. Also as $J$ increases the number of wavelet coefficients within
each additional level systematically reduces, so to derive meaningful
statistical parameters at these upper levels of decomposition might not always
be possible.

Our proposed method minimises the use of non-required frequency ranges within
the signal, reconstructs the signal to its original length and overcomes some
of the possible shortfalls of the DWT. It also uses the interaction between
subsets of electrodes during an epoch. This interaction or connectivity can be
seen to vary significantly between the non-seizure and seizure periods. We
utilise statistical parameters from our reconstructed signal and electrode
connectivity as extracted features for classification. The organisation of
this paper is as follows; Section 2 outlines some related works to our
proposed method, Section 3 presents our method, Section 4 outlines
experimental results together with compared works. Concluding remarks are
presented in Section 5, followed by an Appendix which includes a Glossary.

##  2 Related Work

Time and frequency domains are commonly used in feature extraction methods for
EEG signals. To analyse the frequency domain features, the Discrete Fourier
Transform (DFT) has been applied, followed by a form of power spectral density
analysis of the EEG signals (Lee et al., 2014). Decomposing a signal in terms
of its frequency content, using the DFT which relies on sinusoids, results in
fine resolution within the frequency domain. However the Fourier series
representation is not very effective at time resolution (Vikash, 2015).
Another method of decomposing signals is the DWT, where in a wavelet
representation, we represent the signal in terms of functions that are
localised both in time and frequency. Here we have an adaptive time frequency
window, which provides good time resolution at high frequencies, and good
frequency resolution at low frequencies (Polikar, 1994).

Siddiqui et al. (2018) generated nine statistical parameters from EEG data to
extract time domain features and reduce dimension. The duration of data
(Epoch) was also altered, starting at 10 sec, and reducing in steps to 0.025
secs. The classifiers chosen there were decision forests, an ensemble of
decision trees. A concern with this methodology is that one is using all the
available information or frequency bandwidth within the signal. Such
techniques as DFT or DWT were not used to minimise the inclusion of some
sections of the signal’s bandwidth, it is quite possible that additional
noise111Here we define noise as frequency ranges unrelated to a seizure event.
may have been introduced. Hence the likelihood of introducing error into the
classification of the EEG signal.

EEG signals have also been decomposed into time–frequency representations
using the DWT (Parvez and Paul, 2013; Faust et al., 2015). This enables the
EEG signals to be decomposed into several sub-bands, this decomposition is
called Multi Resolution Analysis (MRA). By Parseval’s Theorem, the percentage
distribution of energy features in the EEG signal may be extracted at each of
these different sub-bands (Omerhodzic et al., 2013). The features extracted
from the wavelet coefficients at various levels or different frequency bands
can be used to determine the characteristics of the signal. Use of the DWT
continues to exhibit improving results on seizure detection (Faust et al.,
2015). One concern here is that with the DWT, as we increase the levels of
decomposition to extract the relevant frequency bandwidths we have fewer
wavelet detailed coefficients $\\{d_{i,j}\\}$, within these higher levels of
decomposition222Here $d_{i,j}$ represents the $i^{th}$ detail wavelet
coefficient at the $j^{th}$ scale or level and $\mathcal{D}_{j}$ represents
all the detail coefficients at level $j$. Similarly for the smooth wavelet
coefficients, $s_{i,J}$ and $\mathcal{S}_{J}$.. That is
$\\{d_{1,j},d_{2,j},\ldots,d_{k,j}\\}\in\mathcal{D}_{j}$ has fewer elements as
$j$ increases: $0\leq j\leq J$. Hence the difficulty in deriving meaningful
statistical parameters at the higher decomposition levels.

Chen et al. (2017) uses various different DWTs for seizure detection, also
choosing from nine statistical features derived from each of the wavelet
generated frequency bands then classifying (via SVM), the multi-channel EEG
recordings from Goldberger et al. (2000). The empirical evidence shown there
outlines that the choice of mother wavelet does not have a significant effect
on seizure detection, however it is very sensitive to decomposition level if
the features of non-seizure/seizure EEGs exhibit significant difference in the
different frequency bands. In this particular study each Epoch of EEG data
used was 20 seconds in duration. Of concern here is that many frequency bands
and the resulting DWT coefficient features could be redundant, causing
inaccuracy and unnecessary high computational cost. If we decompose to 6
levels:
$\mathcal{D}_{1},\mathcal{D}_{2},\ldots,\mathcal{D}_{6}\;\&\;\mathcal{S}_{6}$
and use all the wavelet coefficients within each of these levels to derive
statistical parameters, the result set of attributes to be used for
classification becomes quite large for multivariate EEG data.

While use of the DWT within EEG analysis is often encountered, however use of
the variance within the wavelet coefficients resulting from the transform,
mainly consider this wavelet variance as a selected feature only. This is
usually the standard deviation of the detail wavelet coefficients at each
decomposition scale or a linear combination of the variance across the
different wavelet decomposition scales (Janjarasjitt, 2010; Ashok and
Mahalakshmi, 2017). However little if any consideration is given to the
covariance between different signals/nodes, which may be calculated similarly
as the wavelet variance, this is called wavelet covariance, see definition 5
Appendix A.2, as well as Equations 9, 10 and 11. The covariance ( and wavelet
correlation, see Equation 12 ) can be used to derive electrode/node
connectivity during an epoch to provide additional features for
classification, especially if the connectivity between electrodes/nodes alters
significantly between the non-seizure and seizure states. As an example this
is shown in Figure 2, highlighting an epoch for both non-seizure and seizure
states, displaying an electrode’s number of connections to other electrodes.
It can be seen that during the non-seizure state, there exists numerous
connections between the electrodes, i.e. 6 electrodes have connections ranging
from 5 to 10 other electrodes and another 6 electrodes have connections
ranging from 15 to 20 other electrodes. Where as in comparison to non-seizure
here no electrode has more than eight connections to other electrodes and only
one electrode has a maximum of eight connections at the chosen wavelet
correlation threshold.

The discrimination of EEG events, using statistical parameters derived from
multivariate data as attributes are likely to include considerable information
unrelated to the event under study, noticed as noise and may induce error.
However using the DWT in an attempt to reduce noise by selecting the required
frequency bands (wavelet decomposition levels), problems may still arise. One
obvious issue is length of signal, if not a multiple of $2^{J}$ where
$J\in\mathbb{Z}^{+}$, then additional components would need to be added in as
extra coefficients or they may be simply omitted. Either option may introduce
bias, hence possibly confuse classification.

Figure 2:  Electrode Connections with correlation threshold = 0.125

##  3 Our Proposed Technique

As we are interested in the classification of neurological events, labelled as
either non-seizure or seizure states, where the data is presented in a
multivariate time series format. Here we present a novel approach for feature
extraction and classification using the inherent capabilities of the Maximum
Overlap discrete Wavelet Transform (MODWT, Appendix A.3.2) to:

$\displaystyle\blacklozenge$ decompose and then reconstruct the EEG signal
minimising the noise components $\displaystyle\blacklozenge$
$\displaystyle\left.\begin{array}[]{l}\text{construct a record of inter-
electrode connectivity and}\\\ \text{develop a metric for communication
efficiency of electrodes}\end{array}\right\\}\textit{ from wavelet
covariance.}$

From these wavelet derived methods, we build our attribute space and apply an
ensemble classifier to the transformed data.

The scenario we are considering is to have a database consisting of labelled
EEG signals from a number of previous or existing patients. We subset the
signal into smaller periods of time (Epochs) and label each epoch (non-seizure
or seizure), determined by where the epoch occurs in time, within the signal.
We construct a model from this training data, and apply the model on incoming
unlabelled EEG data, an epoch at a time, from a patient to classify/detect
non-seizure and seizure events.

When using raw data and deriving statistical parameters as attributes, then
our method if using same number of derived statistical parameters, we notice
slightly more attributes per record. However, we reduce the number of
attributes when compared to standard DWT decomposition methodology.

###  3.1 Main Steps

We outline the main steps of our method in this subsection, then explain each
of these steps in greater detail in the following subsections. With the raw
data chosen then after subsetting the data for each patient into the required
time duration or epoch, we begin a three step process to attain the required
attributes, then combine all into a data set followed by the 4th step,
Classification.

  1. Step 1.

Preparation of Attributes from statistical parameters of reconstructed signal.

  2. Step 2.

Preparation of Attributes from Connectivity of Electrodes

  3. Step 3.

Preparation of Attributes from Global Efficiency of Electrodes

  4. Step 4.

Assigning Class labels and Training Classifiers

###  3.2 Step 1: Preparation of Attributes from statistical parameters of
reconstructed signal.

With data, originating from an EEG skull cap sampling method, see Figure 1 for
an example of the standard EEG “10-20” skull cap electrode configuration. The
data appears as a time series or signal from each electrode. The non-seizure
and seizure periods are previously identified. We segment the signal into
epochs of suitable duration and attach class labels. Any epoch with part of
its duration within the identified seizure period is defined as seizure,
elsewhere non-seizure. As an example see Figure 3. Each segment along the
X-Axis within a left and right arrow represents a epoch. The segment on the
left hand end, entirely falls during a non seizure event and hence is labelled
as non-seizure. The 2nd and 3rd epochs from the left, either partially or
fully occur during a seizure event and hence are labelled as seizure. We apply
to each epoch the wavelet transform.

The decomposition function of the wavelet transform may be represented as a
tree of low and high pass filters (LPF & HPF), with each step further
decomposing the low pass filter, Figure 4 displays such decomposition to 3
levels. In this example the resultant transform would consist of the levels
$\\{H1,H2,H3\;\&\;L3\\}$ and the usual nomenclature of the MODWT for these
levels is
$\\{\widetilde{\mathcal{D}_{1}},\widetilde{\mathcal{D}_{2}},\widetilde{\mathcal{D}_{3}}\;\&\;\widetilde{\mathcal{S}_{3}}\\}$
respectively. Each of these MODWT levels consists of wavelet coefficients,
with $\\{d_{i,j}\\}\in\widetilde{\mathcal{D}_{j}}$ as detail coefficients and
$\\{s_{i,J}\\}\in\widetilde{\mathcal{S}_{J}}$ for smooth
coefficients333Similar to Section 2 for the DWT here
$\widetilde{\mathcal{D}_{j}}$ and $\widetilde{\mathcal{S}_{J}}$ represent the
MODWT levels $j$ or $J$ : $0\leq j\leq J$. Figure 5 shows the wavelet
coefficients resulting from a 3 level MODWT, at each of the detail and smooth
levels, across time for an epoch.

EEG waves are conventionally grouped into frequency ranges and named: delta,
theta, alpha, beta and gamma waves, where each set of these waves occupies a
different non overlapping frequency bandwidth. Then dependent upon the
sampling frequency of the initial data, we select our level of wavelet
decomposition to align as close as possible with the frequency bandwidths
required, to analyse the EEG signal (Elsa Jacob et al., 2018). However, some
frequency bands are unlikely to be related to the EEG event we seek to
classify, as it has been noted that EEG signals do not have any useful
frequency components above 30 Hz (Subasi, 2007). Therefore after wavelet
decomposition, we omit or smooth the non-related frequency bands and
reconstruct our signal using the inverse wavelet transform444The inverse
wavelet transform is fully documented (Percival and Walden, 2000, Chapter 5)..
Figure 6 shows an original signal compared to the reconstructed signal where
MODWT level $\widetilde{\mathcal{D}_{1}}$, has been omitted prior to applying
the inverse wavelet transform.

For each electrode (say, 1 to n), after decomposition and reduction of
unwanted noise, the inverse transform produces a reconstructed signal of same
length as the original data. We then take statistical parameters being:
minimum, maximum, mean, standard deviation and normalised energy555Normalised
energy ($NormE$) here is the sum of the squares of the time series values
divided by the number of elements in the time series, i.e. Norm.Energy =
$\frac{\sum_{i}^{k}X_{i}^{2}}{k}$ from each electrode’s reconstructed signal
and form the attributes. We have for each epoch; say epoch $i$
$min_{i1},max_{i1},mean_{i1},std_{i1},normE_{i1},min_{i2},max_{i2},mean_{i2},\ldots,mean_{in},std_{in},NormE_{in}$.
We use the statistical parameters to form a list of attributes representing a
single epoch. Table 1 represents the data set produced in this step. The
following steps add additional attributes to the data set as explained later.
A data set is generated from the EEG signal of each patient in order to create
patient specific training data set. One may also create a training data set
combining all patients. Algorithm 1 in Section 4 outlines this Step 1 in a
structured procedure format.

Figure 3: Data from a single electrode, with seizure duration defined and
Epochs labelled.  Figure 4: Wavelet Decomposition to 3 levels, from Grant and
Islam (2019) Figure 5: Wavelet Coefficients from 3 level decomposition, for an
epoch Figure 6: FP1 and Reconstructed FP1, omitting level
$\widetilde{\mathcal{D}_{1}}$, for an epoch Table 1: Example: Data set
resulting from Step 1

| Attributes derived from each electrode’s signal
---|---
$Individal_{j}$ | Elec-1 | Elec-2 | … | Elec-n
Epoch 1 | $min_{11},max_{11},\ldots,normE_{11}$ | $min_{12},max_{12},\ldots,normE_{12}$ | … | $min_{1n},max_{1n},\ldots,normE_{1n}$
Epoch 2 | $min_{21},max_{21},\ldots,normE_{21}$ | $min_{22},max_{22},\ldots,normE_{22}$ | … | $min_{2n},max_{2n},\ldots,normE_{2n}$
⋮ | ⋮ |  |  | ⋮
Epoch m | $min_{m1},max_{m1},\ldots,normE_{m1}$ | $min_{m2},max_{m2},\ldots,normE_{m2}$ | … | $min_{mn},max_{mn},\ldots,normE_{mn}$

###  3.3 Step 2: Preparation of Attributes from Connectivity of Electrodes

Using the MODWT to decompose the time series from each electrode per epoch, we
derive wavelet variance, and similarly wavelet covariance and hence wavelet
correlation between each of the transformed signals from each electrode.
Wavelet variance is calculated at each level of wavelet decomposition and is
equal to the variance of the wavelet coefficients at that level. The sum of
the wavelet variances from all decomposition levels equals the sample variance
of the time series, similarly for wavelet covariance between each of the
transformed signals. Covariance provides a measure of joint variation for two
random sequences. The covariance of a random sequence with itself would be the
variance. see Appendix A.3.2, A.3.3 and (Percival and Walden, 2000, Chapter
8).

Wavelet correlation which is derived from wavelet covariance see Equation 12,
provides us with a metric allowing us to determine association (joint
variability) between the transformed signals from each electrode to another,
during an epoch. By setting a user defined lower bound on the correlation
value between two electrodes, we may define if a connection (or association)
exists between these electrodes. If the correlation value between two
transformed signals is greater than or equal to the user defined lower
bound666The value of this lower bound may be derived by the number of
connections available determined by that specific value, given that we need to
derive statistical parameters from these connection values. then we conclude
that a connection exists between the two during the epoch.

We then represent those connections in an adjacency matrix where a connection
is represented by the value 1. See Table 2 as an example of an adjacency
matrix. Here, Electrode 1 (Elec-1) and Electrode 2 (Elec-2) are connected
since the correlation between the signals, as calculated by Equation 12 is
higher than the user defined lower bound. From the derived adjacency matrix we
take the sum of the rows.
In Table 2, Electrode 1 is directly connected to two other electrodes i.e
Elec-2 & Elec-3, therefore sum of the row is 2. That is the degree of
connection is equal to 2.
Electrode 3 is directly connected to five other electrodes, therefore the sum
of row is 5. Hence the degree of connection is equal to 5.
For each electrode in our network we now have a degree of connection or number
of direct connections during the epoch. Figure 2 is a representative example
of this. From these row sums (or degree of connection) we derive the
statistical parameters: mean, maximum, 1st Quartile, 3rd Quartile and
skewness. We represent these as: $meanC,maxC,\ldots,skewC$. This provides us
with additional attributes for the epoch in question. These additional
attributes for each epoch are combined with the data set outlined in Step 1
(as shown in Table 1) increasing its dimension. The updated data set is shown
in Table 3. A procedural format of this Step 2 is shown in Algorithm 1.

Table 2: Example: Adjacency matrix

| Elec-1 | Elec-2 | Elec-3 | Elec-4 | Elec-5 | Elec-6 | Elec-7 | … | Elec-n | Sum
---|---|---|---|---|---|---|---|---|---|---
Elec-1 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | … | 0 | 2
Elec-2 | 1 | 0 | 1 | 0 | 0 | 1 | 1 |  | 0 | 4
Elec-3 | 1 | 1 | 0 | 1 | 1 | 1 | 0 |  | 0 | 5
Elec-4 | 0 | 0 | 1 | 0 | 0 | 0 | 0 |  | 0 | 1
Elec-5 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | … | 0 | 1
Elec-6 | 0 | 1 | 1 | 0 | 0 | 0 | 1 |  | 0 | 3
Elec-7 | 0 | 1 | 0 | 0 | 0 | 1 | 0 |  | 0 | 2
⋮ |  | ⋮ |  |  |  |  | ⋮ |  |  | ⋮
Elec-n | 0 | 0 | 0 | 0 | 0 | 0 | 0 | … | 0 | 0

Table 3: Example: Data set resulting after Step 2

| Attributes derived by Step 1 and Step 2
---|---
| Step 1 | Step 2
$Individal_{j}$ | Elec-1 | … | Elec-n | Connectivity
Epoch 1 | $min_{11},max_{11},\ldots,normE_{11}$ | … |  $min_{1n},max_{1n},\ldots,normE_{1n}$, | $meanC_{1},maxC_{1},\ldots,skewC_{1}$
Epoch 2 | $min_{21},max_{21},\ldots,normE_{21}$ | … |  $min_{2n},max_{2n},\ldots,normE_{2n}$, | $meanC_{2},maxC_{2},\ldots,skewC_{2}$
⋮ | ⋮ |  | ⋮ | ⋮
Epoch m | $min_{m1},max_{m1},\ldots,normE_{m1}$ | … |  $min_{mn},max_{mn},\ldots,normE_{mn}$, | $meanC_{m},maxC_{m},\ldots,skewC_{m}$

###  3.4 Step 3. Preparation of Attributes from Global Efficiency of
Electrodes.

Section 3.3 shows that from the wavelet coefficients we have a methodology for
the construction of an adjacency matrix derived from the correlation between
the electrodes. From that adjacency matrix we are able to construct a visual
representation of the network or Graph. Figure 7 shows a sub-graph derived
from only the first seven electrodes of Table 2. In Figure 7 an electrode is
represented as a numbered node and a pair of electrodes/nodes having a
correlation above the threshold are connected through an edge between them.
Figure 7, has seven nodes/vertices (i.e. N = 7) and nine edges. Each single
edge has an assumed unit length of 1. From this constructed network diagram
(and hence the adjacency matrix) we derive how “efficiently”777Efficiency of
$\text{node}_{i}$ is related to the sum of the inverses of the minimum
distances between $\text{node}_{i}$ and all other nodes in the network. the
electrodes are communicating across their network of connections. We calculate
the global efficiency ($Eglob_{i}$) for each node.

For every electrode $i$ we compute the reciprocal of shortest path length
($R_{i,j}$) to each other electrode $j$. For example the shortest path length
from electrode 5 to Electrode 7 is three, hence $R_{5,7}$ is $\frac{1}{3}$.
This is shown in the corresponding cell (encircled) in Table 4. For each
electrode $i$ we then then compute $\frac{\sum\limits_{j=1}^{N}R_{i,j}}{N-1}$
where $N$ is the number of nodes and $R_{i,i}=0$. The last row of Table 4
shows these values. We call these values global efficiency ($Eglob_{i}$)
values of an $\text{electrode}_{i}$ for an epoch, see definition 7.

Figure 7:  A Visual representation of Table 2 Table 4: Example: Global Efficiency | Elec1 | Elec2 | Elec3 | Elec4 | Elec5 | Elec6 | Elec7
---|---|---|---|---|---|---|---
Elec1 | 0 | 1 | 1 | 1/2 | 1/2 | 1/2 | 1/2
Elec2 | 1 | 0 | 1 | 1/2 | 1/2 | 1 | 1
Elec3 | 1 | 1 | 0 | 1 | 1 | 1 | 1/2
Elec4 | 1/2 | 1/2 | 1 | 0 | 1/2 | 1/2 | 1/3
Elec5 | 1/2 | 1/2 | 1 | 1/2 | 0 | 1/2 | 1/3
Elec6 | 1/2 | 1 | 1 | 1/2 | 1/2 | 0 | 1
Elec7 | 1/2 | 1 | 1/2 | 1/3 | 1/3 | 1 | 0
Sum/(N-1) | 2/3 | 5/6 | 11/12 | 5/9 | 5/9 | 3/4 | 11/18

From these global efficiency values for an epoch as shown in the last row of
Table 4 we derive statistical parameters: mean, maximum, 1st Quartile, 3rd
Quartile & skewness. We represent them as $MeanG_{i},MaxG_{i},\dots.SkewG_{i}$
for the $i$th epoch. These attributes are again added to the data set as shown
in Table 5. Again this step in procedural form is seen in Algorithm 1.

Table 5: Example: Data set resulting after Step 3.

| Attributes derived by Step 1, Step 2 & Step 3
---|---
| Step 1 |  | Step 2 | Step 3
$Individal_{j}$ | Elec-1 to n |  | Connectivity | Global Efficiency
Epoch 1 |  $min_{11},max_{11},\ldots,normE_{1n}$, |  |  $meanC_{1},maxC_{1},\ldots,skewC_{1}$, | $MeanG_{1},MaxG_{1},\dots SkewG_{1}$
Epoch 2 |  $min_{21},max_{21},\ldots,normE_{2n}$, |  |  $meanC_{2},maxC_{2},\ldots,skewC_{2}$, | $MeanG_{2},MaxG_{2},\dots SkewG_{2}$
⋮ | ⋮ |  | ⋮ | ⋮
Epoch m |  $min_{m1},max_{m1},\ldots,normE_{mn}$, |  |  $meanC_{m},maxC_{m},\ldots,skewC_{m}$, | $MeanG_{m},MaxG_{m},\dots SkewG_{m}$

Table 6: Example: Data set resulting from the 3 Steps and attaching label

| Attributes derived from each Step
---|---
$Individual_{j}$ | Step 1 | Step 2 | Step 3 | Label
Epoch 1 |  $min_{11},max_{11},\ldots,normE_{1n}$, |  $meanC_{1},maxC_{1},\ldots,skewC_{1}$, |  $MeanG_{1},MaxG_{1},\dots SkewG_{1}$, | non-seizure
Epoch 2 |  $min_{21},max_{21},\ldots,normE_{2n}$, |  $meanC_{2},maxC_{2},\ldots,skewC_{2}$, |  $MeanG_{2},MaxG_{2},\dots SkewG_{2}$, | seizure
⋮ | ⋮ | ⋮ | ⋮ |  ⋮
Epoch m |  $min_{m1},max_{m1},\ldots,normE_{mn}$, |  $meanC_{m},maxC_{m},\ldots,skewC_{m}$, |  $MeanG_{m},MaxG_{m},\dots SkewG_{m}$, | non-seizure

###  3.5 Step 4. Assigning Class labels and Training Classifiers

In this step we first assign the class labels to the epochs. That is for every
row/epoch in Table 5 we assign an additional column with a class label, see
Table 6. If the epoch entirely falls within a non-seizure period then the
epoch is labelled as non-seizure (as also explained in Section 3.2 and Figure
3) otherwise the epoch is labelled as seizure. As the consequence of
misclassification of a seizure record as a non-seizure record can be much
higher than a non-seizure record with a seizure prediction, we use a cost
sensitive classifier. Thus we aim to reduce the misclassification of seizure
records with non-seizure records (i.e. False Negative predictions). Once we
train our cost sensitive classifier on the training data set, being Table 6,
we are ready to predict/detect seizure in the new unlabelled signal. When we
receive new unlabelled signal we first convert into a record as shown in Table
5. The trained classifier is then used to label the record.

input : A Signal input $DB_{I}$, of individuals’ EEG recordings

output : $A_{i,j}$ : Output for individual per epoch from the Cost sensitive
Classifier.

Initialization: Build transformed database

$I_{j}\leftarrow DB_{I}$ //Download initial data file per individual

for  _All individuals $I_{j}$_ do

$I_{i,j}\leftarrow I_{j}$ //Segment individual’s data into 1 second Epochs

for  _Each Individual by Epoch $I_{i,j}$_ do

Step 1.

for  _Each Electrode $I_{i,j,e}$_ do

$\mathbf{\widetilde{W}_{i}}\leftarrow I_{i,j,e}$ // Apply MODWT, decompose to
required levels

$I_{i,j,e,R}\leftarrow\mathbf{\widetilde{W}_{i}}$// Omit noise levels, use
Inverse MODWT to reconstruct reduced signal

$Sp_{i}\leftarrow I_{i,j,e,R,}$ // Derive statistical parameters from new
signal

end for

$\text{RW}_{i,j,e}\leftarrow Sp_{1},Sp_{2},\dots,Sp_{n}$ //Form attributes
from statistical parameters

end for

for  _Each Individual by epoch $I_{i,j}$_ do

Step 2.

for  _Each Electrode $I_{i,j,e}$_ do

$\mathbf{\widetilde{W}_{i}}\leftarrow I_{i,j,e}$ // Apply MODWT, decompose to
first wavelet detail level

end for

for  _Each electrode pair $\\{I_{i,j,\tau},I_{i,j,\zeta}\\}$ _ do

$cov\\{\widetilde{W}_{i,\tau},\widetilde{W}_{i,\zeta}\\}\leftarrow\widetilde{W}_{i,\tau},\widetilde{W}_{i,\zeta}$
// derive wavelet covariance between electrodes

$\rho_{\tau\zeta}\leftarrow
cov\\{\widetilde{W}_{i,\tau},\widetilde{W}_{i,\zeta}\\}$//derive correlation
from standardised covariance

end for

$Thresh\leftarrow 0<t_{h}\leq 1$ // determine correlation threshold level
$t_{h}$

$Adj_{ee}\leftarrow 1\text{where}\,\rho\geq Thresh$ //construct adjacency
matrix

$Sp_{i}\leftarrow\sum r_{i}$ //derive statistical parameters from row sums of
adj matrix

$C_{i,j}\leftarrow Sp_{1},Sp_{2},\dots,Sp_{k}$ // Use statistical parameters
to form attributes

$\text{RW}_{i,j,e},C_{i,j}\leftarrow C_{i,j}$ // combine with attributes from
Step 1.

Step 3.

$Eglob_{i}\leftarrow Adj_{ee}$ //From Adjacency matrix derive Global
efficiency

$Sp_{i}\leftarrow Eglob_{i}$ // derive statistical parameters

$G_{i,j}\leftarrow Sp_{1},Sp_{2},\dots,Sp_{k}$ // Use statistical parameters
to form attributes

$\text{RW}_{i,j,e},C_{i,j},G_{i,j}\leftarrow G_{i,j}$ // combine with
attributes from Step 2

end for

Step4. $DB_{T}\leftarrow\text{RW}_{i,j,e},C_{i,j},G_{i,j},label$ // add label
and build database (for training classifier upon)

end for

$A_{i,j}\leftarrow DB_{T}$ //Apply Cost sensitive Classifier
$\qquad\blacksquare$

Algorithm 1 EEG Analysis via MODWT decomposition and Wavelet Correlation

##  4 Experimental Results

###  4.1 Data

We used data from Goldberger et al. (2000), which was collected at the
Children’s Hospital Boston, Massachusetts (MIT). This consists of EEG
recordings from paediatric subjects with intractable seizures. There are
recordings from 24 different investigations (cases) from 23 subjects (5 males,
17 females and 1 with gender not specified). The International 10-20 system of
EEG electrode positions (for 23 electrodes) and nomenclature was used for
these recordings with the data sampled at 256Hz. More details about the data
set may be found from https://archive.physionet.org/pn6/chbmit/. For this
study we used a single file from each of 12 individuals, these individuals
were identified in the MIT data by the case numbers; 1, 2, 3, 4, 5, 6, 7, 8,
10, 11, 23 and 24. The cases were selected on the basis of

  * •

the data across cases was reported with same identical contiguous electrode
sequence,
i.e.$\\{FP1,F7,T7,\ldots,FT10,P8\\}$.

  * •

each case was a different individual

  * •

limited selection to files up to two hours in duration

  * •

no seizure event within the first 10 minutes.

Most of the selected files were approx one hour in duration. As the files have
256 data points per second for each electrode signal, therefore one hour
equates to $3600\times 256=921600$ data points, per electrode.

A few files were two hours in length, for these files we selected either the
first or second hour depending upon which hour had the first seizure event.
This raw data is presented with indicated periods of seizure, we subset the
data into one second intervals (i.e. an epoch length of one second) and
attached labels indicating the class of that epoch, i.e. non-seizure or
seizure. We have approx 1 hour of recording from each patient. For two of the
patients there were no seizures recorded within the hour sample, we still
included these two into our training set. Of these two, one patient simply had
no seizures recorded in files of up to two hours duration, the other patient’s
files had the sequence in which the electrodes were reported, altered before
first file with a seizure occurred.

###  4.2  Synopsis of results and comparisons

In this study True Positive (TP) refers to events where epochs labelled
seizure were classified as seizure. False Negative (FN) refers to events where
seizure labelled epochs were classed as non-seizure. True negative (TN) refers
to epochs labelled non-seizure and classified non-seizure and False Positive
(FP) refers to a an epoch labelled non-seizure and classified seizure. We use
the following metrics or indicators: Accuracy, Recall, Precision and f-score.
They are defined as,

$\displaystyle Accuracy\,\%$ $\displaystyle=$
$\displaystyle\frac{TP+TN}{TP+FN+TN+FP}\times 100$ $\displaystyle Recall$
$\displaystyle=$ $\displaystyle\frac{TP}{TP+FN}$ $\displaystyle Precision$
$\displaystyle=$ $\displaystyle\frac{TP}{TP+FP}$ $\displaystyle f-score$
$\displaystyle=$
$\displaystyle\frac{2}{(Recall)^{-1}+(Precision)^{-1}}=2\times\frac{Precision\times
Recall}{Precision+Recall}$

We use these indicators to compare our methodology to the other methods. As we
are interested in correctly classifying seizure we use a cost sensitive/class
imbalanced algorithm. Here we choose CSForest (Siers and Islam, 2015) and set
the cost penalty structure as shown in Table 7.

Table 7:  Penalty Cost Structure for CSForest.  | Actual Value |
---|---|---
| Positive | Negative |
Predicted | Positive | $TP=1$ | $FP=1$ |
Value | Negative | $FN=15$ | $TN=0$ |

We arrived at this cost structure by training on all individuals in the
transformed data set, except individual 1. From this initial reduced data set,
we implemented a 66% training 34% test split on this data set for validation
of the CSForest model (i.e we repeated this training using 66% of data for
training and testing on 34%, altering the penalty cost, with the objective of
increasing the f-score while maintaining a high Recall value). As this penalty
cost structure shown in Table 7 provided us with a close to optimal result, in
regards to f-score888The f-score is the harmonic mean of precision and recall,
sometimes called a balanced f-score., we left the CSForest cost structure
unchanged throughout the rest of experiment across all individuals in our data
set.

We also compared our results to the two other methods outlined in Section 4.3
and used leave-one-out cross validation for all methods. By using the macro
average method, see Definition 6 Appendix A.2, for determining an average
f-score across individuals per method, our method returned a higher average
f-score than the other two comparative methods, as shown in Table 8. For the
f-score results at an individual level see Table 9. Our technique (labelled
Method No.3 in Tables 8 & 9) usually returns better results in respect to this
balanced f-score across these individuals. We choose the f-score as it weights
both Precision and Recall equally.

Table 8: Average f-score results via Method. Method No. | Method origin | Average f-score
---|---|---
1 | Siddiqui et al. (2018) | 0.37
2 | Chen et al. (2017) | 0.42
3 | Our Proposed Technique | 0.58
Table 9: f-scores via Individual across Methods Comparative results
---
Method Number Method Number
f-scores f-scores
Ind. No. | 1 | 2 | 3 | Ind. No. | 1 | 2 | 3
1 | 0.136 | 0.178 | 0.646 | 5 | 0.664 | 0.550 | 0.735
2 | 0.081 | 0.062 | 0.326 | 7 | 0.328 | 0.724 | 0.840
3 | 0.213 | 0.042 | 0.385 | 8 | 0.32 | 0.037 | 0.455
4 | 0.039 | 0.0 | 0.185 | 10 | 0.556 | 0.87 | 0.881
23 | 0.018 | 0.0 | 0.395 | 24 | 0.562 | 0.584 | 0.574

###### Remark 1

There were two other individuals included in our data set and neither of these
had a seizure state in the file selected. Since there is no seizure record,
all three classification methods all provided near perfect results in regards
to accuracy and f-score for these two individuals, (not shown).

####  4.2.1  Metrics chosen

In our situation we wish maximise True Positives, i.e. those with seizure
classified as seizure as well as minimise False Negatives, i.e. a seizure
classified as a non-seizure. The metric Recall is important here. Recall
provides a measure of the True Positives, divided by the sum of True Positives
and False Negatives. Alternatively, Precision which is equal to True positives
divided by the sum of True Positives and False Positives. Here False Positives
are those in a non-seizure state but classified as seizure. Therefore while we
aim to detect all seizures as seizures, classifying non-seizures as seizures
is at worst, a waste of resources. However classifying seizures as non-
seizures may lead to catastrophic outcomes.

Table 10: Statistical results via Individual across Methods

Comparative results
---
Method Number Method Number
Ind. No. |  Metric | 1 | 2 | 3 | Ind. No. |  Metric | 1 | 2 | 3
1 |  Accuracy% | 98.94 | 98.87 | 99.03 | 5 |  Accuracy% | 97.47 | 98.00 | 98.56
|  Recall | 0.073 | 0.098 | 0.780 |  |  Recall | 0.776 | 0 .379 | 0.621
|  Precision | 1 | 1 | 0.552 |  |  Precision | 0.581 | 1 | 0.900
|  f-score | 0.136 | 0.178 | 0.646 |  |  f-score | 0.664 | 0.550 | 0.735
2 |  Accuracy% | 74.77 | 97.46 | 94.35 | 7 |  Accuracy% | 96.48 | 98.87 | 99.28
|  Recall | 0.488 | 0.037 | 0.598 |  |  Recall | 0.330 | 0.567 | 0.732
|  Precision | 0.044 | 0.200 | 0.224 |  |  Precision | 0.327 | 1 | 0.986
|  f-score | 0.081 | 0.062 | 0.326 |  |  f-score | 0.328 | 0.724 | 0.840
3 |  Accuracy% | 97.94 | 98.72 | 96.89 | 8 |  Accuracy% | 95.39 | 95.61 | 92.08
|  Recall | 0.208 | 0.021 | 0.729 |  |  Recall | 0.241 | 0.018 | 0.735
|  Precision | 0.217 | 1 | 0.261 |  |  Precision | 0.476 | 1 | 0.33
|  f-score | 0.213 | 0.042 | 0.385 |  |  f-score | 0.32 | 0.037 | 0.455
4 |  Accuracy% | 98.67 | 98.36 | 94.83 | 10 |  Accuracy% | 97.77 | 99.61 | 99.61
|  Recall | 0.02 | 0.0 | 0.429 |  |  Recall | 0.82 | 0.77 | 0.813
|  Precision | 0.50 | 0.0 | 0.118 |  |  Precision | 0.420 | 1 | 0.963
|  f-score | 0.039 | 0.0 | 0.185 |  |  f-score | 0.556 | 0.87 | 0.881
23 |  Accuracy% | 96.89 | 96.86 | 96.25 | 24 |  Accuracy% | 98.91 | 98.97 | 98.81
|  Recall | 0.009 | 0.0 | 0.389 |  |  Recall | 0.50 | 0.520 | 0.580
|  Precision | 1 | 0.0 | 0.40 |  |  Precision | 0.641 | 0.667 | 0.569
|  f-score | 0.018 | 0.0 | 0.395 |  |  f-score | 0.562 | 0.584 | 0.574

Hence we consider Recall slightly more relevant than Precision for this
exercise as Precision does not take into account False Negatives. Precision
does however consider False Positives. Table 10 which highlights all the
chosen indicators at an individual level, demonstrates that our proposed
method (Method Number 3), returns higher values of Recall across nearly all
individuals. In some instances Precision values from the comparative methods
are higher than from our method, i.e. Individual No 1. This could simply be a
situation of a small proportion of the overall number of seizures classified
as seizures but with no non-seizures classified as seizures. There is no
indication in this Precision metric of how many False Negatives were generated
in the classification process. This is why we have chosen f-score as our main
indicator of overall performance, to balance out Recall and Precision.
Accuracy, shown in Tables 10 & 12 is a misleading metric for such imbalanced
data. In our final data set of all individuals we have approx 42816 records of
which 822 are labelled seizure. i.e. 1.92%. Hence if we simply classify all
records as non-seizure then our Accuracy will be 98.08%.

####  4.2.2 Alternative Cost Structure within Classifier

If a major concern was False Positives, we could alter our Recall and
Precision values by changing the penalty cost structure within the CSForest
classifier. In comparison to those values shown in Table 7, by reducing the
penalty cost on False Negatives we may reduce the Recall value, conversely
while increasing the False Positive cost we may increase the Precision value.
As an example of such penalty costs see Table 11.

Table 11:  Alternative Penalty Cost Structure for CSForest.  | Actual Value |
---|---|---
| Positive | Negative |
Predicted | Positive | $TP=1$ | $FP=2$ |
Value | Negative | $FN=2$ | $TN=0$ |

Applying the CSForest classifier to the transformed data for the individuals
No. 1 and No. 8, but using the cost structure as shown in Table 11, we note
from the results in Table 12 that the Recall, Precision and f-score have
altered when compared to Table 10.

Table 12: Statistical results via Individual with Alternative Cost Structure
Our Proposed method with Costs Altered
---
Ind. No. | Metric | Result | Ind. No. | Metric | Result
1 | Accuracy% | 98.89 | 8 | Accuracy% | 96.36
| Recall | 0.024 |  | Recall | 0.191
| Precision | 1.0 |  | Precision | 1.0
| f-score | 0.048 |  | f-score | 0.321

###  4.3 Comparative methods

We compared our classification results to two other published methods, using
the exact same raw data we had used, i.e one second epochs from the MIT data
set, with results mentioned in Section 4.2, Tables 8, 9 & 10. The two works we
compared to are as follows:

  * Method 1.

Siddiqui et al. (2018) used nine statistical parameters, i.e. (Min, Max, …,
Skewness) derived from each electrode signal to build the attribute space,
resulting in $9\times 23=207$ attributes plus label per record, then applied
an ensemble classifier to data, in this case being SysFor (Islam and Giggins,
2012).

  * Method 2.

Chen et al. (2017), utilised the DWT to decompose the electrode signal into
different frequency bands, selected a subset of these frequency bands and
derived statistical parameters from the wavelet coefficients in the chosen
frequency bands (wavelet decomposition scales) then applied the classifier
SVM, (via a leave one out methodology), to their transformed data. This was
repeated for different mother wavelets, (i.e. Daubechies, Haar, Symlets, etc
…). For comparison we used this publication’s methodology for the Haar
wavelet, being six frequency bands together with three statistical parameters:
Minimum, Standard deviation and Skewness. This resulted in $6\times 3\times
23=414$ attributes plus label per record.

###  4.4 Application of our proposed Method in the The Experiments

For this experiment, we used open source software, being R (R Core Team, 2019)
for all Graph Theory methods, Wavelet and Statistical analysis. WEKA (Hall et
al., 2009) was utilised for the Classification step.

####  4.4.1 Step 1: Preparation of Attributes from statistical parameters of
reconstructed signal.

From the initial data that had been segmented into one second epochs, we
initially applied the wavelet transform999Using wavelet filter “s8” which had
been previously used successfully in Achard (2012) & Chen et al. (2017).. As
the original data had been sampled at 256Hz, using the MODWT decomposition we
could decompose the signal into six frequency bands to match to the advised
bandwidths, our decomposition/frequency bandwidth range shown in Table 13.

We observed that after such wavelet decomposition, only minor signal energy
was apparent in the first two wavelet decomposition bands101010These sub-bands
are also called Crystals $\widetilde{\mathcal{D}_{1}}$ &
$\widetilde{\mathcal{D}_{2}}$, which represent 32Hz to 128Hz. This was noticed
across many electrodes, both in non-seizure and seizure states, as an example
from one individual for one electrode see Figure 8.

Table 13: Frequency Band decomposition Level | Frequency | Band
---|---|---
$\widetilde{\mathcal{D}_{1}}$ | 64-128Hz | High Gamma
$\widetilde{\mathcal{D}_{2}}$ | 32-64Hz | Gamma
$\widetilde{\mathcal{D}_{3}}$ | 16-32Hz | Beta
$\widetilde{\mathcal{D}_{4}}$ | 8-16Hz | Alpha
$\widetilde{\mathcal{D}_{5}}$ | 4-8Hz | Theta
$\widetilde{\mathcal{S}_{5}}$ | 0-4Hz | Delta

Signal energy at each sub-band (Crystal), e.g.
$\widetilde{\mathcal{D}_{4}}\sim$ Alpha, is defined similarly as in total
energy in signal (time series), see Definition 2, Appendix A.2

As an example: Energy of Alpha sub-band:

$Energy|\widetilde{\mathcal{D}_{4}}|=\sum_{i=1}^{n}|d_{i,4}|^{2}\text{ ,where
}d_{1,4},d_{2,4},\ldots,d_{n,4}\in\widetilde{\mathcal{D}_{4}}$ (1)

(a) Individual 1 Electrode T7, Energy across Wavelet levels, Non-Seizure state

(b) Individual 1 Electrode T7, Energy across Wavelet levels, Seizure state

Figure 8: MODWT Decomposition of 1 second Epochs

After eliminating the the first two wavelet levels (Crystals), we then applied
the inverse wavelet transform to reconstruct the signal, where the non-
relevant frequency ranges had been omitted.

#####  feature extraction.

After application of the inverse wavelet transform on the reduced set of
wavelet coefficients, for the epoch, this resulted in a time series of length
256 points (1 second epoch). For each of these reconstructed signals, per
epoch we derived the statistical parameters and followed the procedure as
outlined in in Section 3.2. This provides in total for the 23 electrodes, 115
attributes for each epoch. Table 1 is an example of the output derived here..

####  4.4.2 Step 2: Preparation of Attributes from Connectivity of Electrodes

By applying the MODWT to our data, using same wavelet filter as in Section
4.4.1, for the time series from each electrode, per epoch we derived wavelet
covariance and hence wavelet correlation between each time series. We then
constructed an adjacency matrix highlighting the direct connections between
the electrodes. This was performed using only the wavelet coefficients from
the first level of MODWT decomposition, i.e.
$d_{i,1}\in\widetilde{\mathcal{D}_{1}}:i\leq n$ and setting a lower bound on
the correlation value of 0.125. Using the lower bound value we determined if a
connection exists between the electrodes if the correlation value is greater
than or equal to the lower bound111111The value 0.125 was chosen as with
wavelet decomposition level 1, this provided over 150 electrode connections
during the non-seizure state and over 30 connections during the seizure state,
per epoch.. We represented those connections in an adjacency matrix, as
described in Section 3.3.

##### feature extraction

From the derived adjacency matrix we calculate the sum of the rows. For each
electrode in our network we now have a degree of connections or number of
direct connections during the epoch. From the row sums (or degree of
connection) we derive the statistical parameters and follow the procedure as
described in Section 3.3. For each epoch this provides us with another five
attributes for each record. Table 3 is an example of this.

The data initially chosen consisted of results from 23 electrodes or channels.
As one channel was the reverse of the other, ie channel 3: T7-P7, channel 19:
P7-T7, we dropped channel 19 when calculating the wavelet covariances and
hence from the adjacency matrix construction, we had 22 row sums to derive the
statistical parameters from, per epoch.

####  4.4.3 Step 3. Preparation of Attributes from Global Efficiency of
Electrodes

From Section 3.3 where wavelet covariance was shown to provide a methodology
for the construction of an adjacency matrix derived from the wavelet
correlation between the electrodes, we derive how “efficiently” the electrodes
are at communicating across their network. This provides us with the global
efficiency per node, as described in Section 3.4

We performed this calculation using the adjacency matrix derived in Step 2 of
our experiment, (Section 4.4.2), for each one second epoch and 22 electrodes
(nodes).

##### feature extraction

For each electrode in our network we now have an efficiency value,
$Eglob_{j}:0\leq Eglob_{j}\leq 1$ or how well that electrode communicates or
transfers information within the network during the epoch. As we have these
values, 22 in total, for each epoch, we derive from them the statistical
parameters as described in Section 3.4 . This provides us with an additional
five parameters. We then place these attributes directly following the
attributes derived in Step 1 and Step 2 of our method, see Table 5 as an
example.

These three steps combined together provides a total 115 + 5 + 5 = 125
attributes.

####  4.4.4 Step 4. Assigning Class labels and Training Classifiers

From the initial three steps in our method, we have constructed a new data set
from the original data, where we have transformed our multivariate EEG data
into a univariate list of attributes per Epoch and attached labels as
described in Section 3.5, where Table 6 is an example of this. Our constructed
database of 12 patients which equates to approx 42,816 records, each record of
125 attributes plus label. We applied the CSForest classifier to our data,
with the cost structure as outlined in Section 4.2, testing each individual
against the training data set constructed from the other 11 individuals in
this study, i.e. leave one out cross validation method. This classification
was undertaken via WEKA. The results of the three different methods per
individual are shown in Table 10.

##  5 Conclusion

In this study, we investigated the use of wavelets and their inherent ability
to provide analysis into the time frequency space. We utilised a wavelet
transform, that is defined for all lengths of sample, able to provide
selection of relative frequencies as well as deriving correlation between
signals to construct different indicators. Together with some basic
methodology from graph theory, we applied these to an open source data set and
compared our results with two other published methods which we applied to
exactly the same data. Our method compared well, providing overall better
results, in regards to f-score and Recall.

Our methodology of utilising three steps (consisting of Wavelet applications
and Graph theory) to construct the attributes to which the classifier will be
applied, does require some additional computation, but returns a similar
number of attributes when compared to usual methods that construct statistical
parameters for attributes derived from raw data. However this method returns
considerably less attributes when compared to methods deriving statistical
parameters from wavelet coefficients across different decomposition levels.

During classification we kept the initial penalty cost structure within the
cost sensitive classifier unchanged. This could be altered by determining the
optimal cost structure each time the training data set is altered or after
inclusion of additional data. The chosen data had 22 electrode positions that
could be used to derive node correlation and connectivity. This limited the
number of connections available at higher levels of wavelet decomposition
without considerably lowering the correlation threshold. These higher levels
of wavelet decomposition may have been better suited for seizure
discrimination however we did not use owing to the limit on derived
connections. For electrode connectivity calculation we used the first Crystal
$\widetilde{\mathcal{D}_{1}}$ as this decomposition level provided a larger
number of connections across the electrodes, especially during seizure, in
comparison to the other wavelet levels.

Also we did not weight the distance between the nodes, simply gave a derived
connection a distance value or path length of one. Further investigation using
data sets with many more nodes together with spatial coordinates which would
enable weighting for distance apart, may prove better at determining electrode
connectivity especially in the surrounding locality of seizure, to provide
enhanced discrimination between the different EEG states.

Conditional upon acceptance, the associated R code will be made available via
Github.

## A

###  A.1 Glossary

  * $A_{i,j}$

Output from classifier, per epoch

  * $Adj_{ee}$

Adjacency Matrix

  * $C_{i,j}$

List of attributes, Step 2

  * $\text{DB}_{T}$

Database of transformed data

  * $d_{i,j}$

$i^{th}$ detail wavelet coefficient at $j^{th}$ level

  * $\mathcal{D}_{j}$

$j^{th}$ level of detail coefficients, DWT

  * $\widetilde{\mathcal{D}_{j}}$

$j^{th}$ level of detail coefficients, MODWT

  * DWT

Discrete Wavelet transform

  * EEG

Electroencephalogram

  * $Eglob_{i}$

Global Efficiency value, $\text{electrode}_{i}$

  * FN

False Negative

  * FP

False Positive

  * $G_{i,j}$

List of attributes, Step 3

  * HPF

High Pass filter

  * Hz

Hertz, One cycle per second

  * $I_{j}$

Individual’s data

  * $I_{i,j}$

Individual’s data by epoch

  * $I_{i,j,e}$

Individual’s data by epoch by electrode

  * $I_{i,j,e,R}$

Individual’s reconstructed data

  * LPF

Low Pass filter

  * MODWT

Maximum Overlap Discrete Wavelet Transform

  * MRA

Multi Resolution Analysis

  * $\rho_{\tau\zeta}$

Correlation between the two electrodes ($\tau\,,\,\zeta$) per epoch

  * Quartiles

values that divide a list of numbers into quarters

  * $r_{i}$

a row in Adjacency matrix

  * $R_{i,j}$

Reciprocal of shortest path between two electrodes

  * $\text{RW}_{i,j,e}$

List of attributes, Step 1

  * $s_{iJ}$

$i^{th}$ smooth wavelet coefficient at $J^{th}$ level

  * $\mathcal{S}_{J}$

$J^{th}$ level of smooth coefficients, DWT

  * $\widetilde{\mathcal{S}_{J}}$

$J^{th}$ level of smooth coefficients, MODWT

  * $Sp_{i}$

a set of statistical parameters

  * Skewness

Measure of asymmetry of a probability distribution.

  * SVM

Support Vector Machine, a classifier

  * $t_{h}$

derived threshold level

  * TN

True Negative

  * TP

True Positive

  * $\mu$V

MicroVolt, one millionth of a volt

  * W

A vector of DWT wavelet coefficients

  * $W_{j}$

vector of detail DWT coefficients, $j^{th}$ level

  * $V_{J}$

vector of smooth DWT coefficients, $J^{th}$ level

###  A.2 Definitions

###### Definition 1

Let a sequence $X_{1},X_{2},\ldots,X_{n}$ represent a time series of $n$
elements, denoted as $\\{X_{t}:t=1,\ldots,n\\},\;X_{t}\in\mathbb{R}$.

###### Definition 2

The energy within a time series $X_{t}$ is the sum of the squared norm of
$X_{t}$,

$\parallel X_{t}\parallel^{2}\equiv\langle
X,X\rangle=\sum\limits_{t=1}^{n}X^{2}_{t}$

###### Definition 3

let $\overline{X}$ denote the sample mean of $X_{t}$

$\overline{X}\equiv\frac{1}{n}\sum\limits_{t=1}^{n}X_{t}$

###### Definition 4

let $\hat{\sigma}^{2}_{X}$ denote the sample variance of $X_{t}$

$\hat{\sigma}^{2}_{X}\equiv\frac{1}{n}\sum\limits_{t=1}^{n}(X_{t}-\overline{X})^{2}$

###### Definition 5

let cov{$X_{t},Y_{t}$} denote the covariance of two random variate sequences
$X_{t}\text{and}\;Y_{t}$

$cov\\{X_{t},Y_{t}\\}=\sum_{t=1}^{n}\frac{(X_{t}-\overline{X})(Y_{t}-\overline{Y})}{n}$

###### Definition 6

macro average Recall (MAR) $\displaystyle=$
$\displaystyle\frac{R_{1}+R_{2}+\dots+R_{N}}{N}$ macro average Precision (MAP)
$\displaystyle=$ $\displaystyle\frac{P_{1}+P_{2}+\dots+P_{N}}{N}$ macro
average f-score $\displaystyle=$
$\displaystyle\frac{2}{\frac{1}{MAR}+\frac{1}{MAP}}=2\times\frac{MAR\times
MAP}{MAR+MAP}$

see Torgo (2014)

###### Definition 7

In a Graph G (or network) the global efficiency measures how the information
is propagating in the entire network. Efficiency = the inverse of the harmonic
mean of the minimum path length $L_{jk}$ between a node j and all the other
nodes k in the graph.

$\text{Eglob}_{j}=\frac{1}{N-1}\sum_{k\in G}\frac{1}{L_{jk}}$

as shown in Achard (2012)

###  A.3 Wavelets

####  A.3.1 DWT

The discrete wavelet transform (DWT) returns a data vector of the same length
as that of the input. Usually, in this new vector many data points are almost
zero. This transform provides an additive decomposition of $J$+1 levels.

$X_{t}\equiv\sum^{J}_{j=1}\mathcal{D}_{j}+\mathcal{S}_{J}$ (2)

where

  * •

$\mathcal{D}_{j}$ $j^{th}$ level of wavelet detail coefficients, related to
changes at scale $j$

  * •

$\mathcal{S}_{J}$ $J^{th}$ level of wavelet smooth coefficients, related to
averages at scale $J$

  * •

$\\{X_{t}:t=1,\ldots,n\\}$ where $n=2^{J}$ :
$J\in\mathbb{Z^{+}},\;X_{t}\in\mathbb{R}$,

The DWT is a linear transform which decomposes $X_{t}$, into $J$ levels giving
$n$ DWT coefficients; the wavelet coefficients are obtained by premultiplying
$X$ by $\mathcal{W}$

$\textbf{W}=\mathcal{W}X$ (3)

  * •

W is a vector of DWT coefficients ($j$th component is $W_{j}$)

  * •

$\mathcal{W}$ is $n\times n$ orthonormal transform matrix; i.e.,
$\mathcal{W}^{T}\mathcal{W}=I_{n}$, where $I_{n}$ is $n\times n$ identity
matrix

  * •

inverse of $\mathcal{W}$ is its transpose,
$\implies\mathcal{W}\mathcal{W}^{T}=I_{n}$

  * •

$\therefore\mathcal{W}^{T}$W = $\mathcal{W}^{T}\mathcal{W}X=X$

W is partitioned into $J+1$ subvectors

  * •

$\textbf{W}_{j}$ has $n/2^{j}$ elements 121212 note:
$\sum^{J}_{j=1}\frac{n}{2^{j}}=\frac{n}{2}+\frac{n}{4}+\dots+2+1=2^{J}-1=n-1$

  * •

$\textbf{V}_{J}$ has one element 131313If we decompose $X_{t}$ to level $I\leq
J$ then $\textbf{V}_{I}$ has $n/2^{I}$ elements.

the synthesis equation for the DWT is:

$X=\mathcal{W}^{T}\textbf{W}=\left[\mathcal{W}^{T}_{1},\mathcal{W}^{T}_{2},\ldots,\mathcal{W}^{T}_{J},\mathcal{V}^{T}_{J}\right]\left[\begin{tabular}[]{c}$\textbf{W}_{1}$\\\
$\textbf{W}_{2}$\\\ {\scriptsize$\vdots$}\\\ $\textbf{W}_{J}$\\\
$\textbf{V}_{J}$\\\ \end{tabular}\right]$ (4)

Equation 4 leads to additive decomposition which expresses $X$ as the sum of
$J+1$ vectors, each of which is associated with a particular scale
$\mathcal{T}_{j}$

$=\sum^{J}_{j=1}\mathcal{W}^{T}_{j}\textbf{W}_{j}+\mathcal{V}^{T}_{J}\textbf{V}_{J}\equiv\sum^{J}_{j=1}\mathcal{D}_{j}+\mathcal{S}_{J}$
(5)

  * •

$\mathcal{D}_{j}\equiv\mathcal{W}^{T}_{j}\textbf{W}_{j}$ is portion of
synthesis due to scale $\mathcal{T}_{j}$, called the jth ’detail’

  * •

$\mathcal{S}_{J}\equiv\,\mathcal{V}^{T}_{J}\textbf{V}_{j}$ is a vector called
the ’smooth’ of the $J$th order

####  A.3.2 MODWT

However another variant of the wavelet transform exists called the Maximum
Overlap Discrete Wavelet Transform (MODWT), where the number of coefficients
returned by the transform at each level of decomposition is the same as that
of the input.

  * •

The MODWT is not orthonormal and is highly redundant.

  * •

The MODWT is defined for $X_{t}$ on all sample sizes, $n$ does not need to be
a power of two

  * •

The MODWT to level $J$ has $(J+1)\times n$ coefficients , the DWT has a total
of $n$ coefficients for any chosen $J$.

  * •

Similarly to the DWT , the MODWT permits an MRA

  * •

Similarly to the DWT , the MODWT permits an analysis of variance to be based
on MODWT coefficients.

As with the DWT, from the MODWT we may also obtain a scale based additive
decomposition. Where the $\,\mathbf{\widetilde{}}\,$ represents the
coefficients are derived from the MODWT.

$X_{t}=\mathcal{\widetilde{W}}^{T}_{j}\widetilde{\textbf{W}}_{j}+\mathcal{\widetilde{V}}^{T}_{J}\widetilde{\textbf{V}}_{J}\equiv\sum^{J}_{j=1}\mathcal{\widetilde{D}}_{j}+\mathcal{\widetilde{S}}_{J}$
(6)

Also a scale based energy decomposition, (can be seen as a basis for ANOVA)

$\parallel
X_{t}\parallel^{2}=\sum^{J}_{j=1}\parallel\widetilde{\textbf{W}}_{j}\parallel^{2}+\parallel\widetilde{\textbf{V}}_{J}\parallel^{2}$
(7)

We may obtain a decomposition of the sample variance, from Definition 4, it
follows:

$\hat{\sigma}^{2}_{X}=\frac{1}{n}\sum_{j=1}^{J}\parallel\widetilde{\textbf{W}}_{j}\parallel^{2}+\frac{1}{n}\parallel\widetilde{\textbf{V}}_{J}\parallel^{2}-\overline{X})^{2}\equiv\sum_{j=1}^{J}\nu_{X}^{2}(\mathcal{T}_{j})$
(8)

We may consider a wavelet variance $\nu_{X}^{2}(\mathcal{T}_{j})$ as a portion
of $\hat{\sigma}^{2}_{X}$ related to changes in averages over scale
$\mathcal{T}_{j}$ i.e., scale by scale analysis of variance. Or the wavelet
variance for scale $\mathcal{T}_{j}$ may be defined as the variance of
$\widetilde{W}_{j}$

$\nu^{2}_{X}(\mathcal{T}_{j})\equiv\text{var}\\{\widetilde{W}_{j}\\}$ (9)

####  A.3.3 Wavelet correlation

If we have two sequences (or signals) $X_{t}\,\&\,Y_{t}$ as per Definition 1,
then with wavelet coefficients $\\{\widetilde{W}_{X,j}\\}$ and
$\\{\widetilde{W}_{Y,j}\\}$, we have: wavelet covariance

$\nu_{XY}(\mathcal{T}_{j})\equiv
cov\\{\widetilde{W}_{X,j},\widetilde{W}_{Y,j}\\}$ (10)

which similar to wavelet variance, provides a decomposition at scale
$\mathcal{T}_{j}$, of the covariance between $\\{X_{t}\\}$ and $\\{Y_{t}\\}$:

$\sum_{j=1}^{\infty}\nu_{XY}(\mathcal{T}_{j})=cov\\{X_{t},Y_{t}\\}$ (11)

Wavelet covariance may be standardised. which yields, wavelet correlation:

$\rho_{XY}\equiv\frac{cov\\{\widetilde{W}_{X,j},\widetilde{W}_{Y,j}\\}}{(var\\{\widetilde{W}_{X,j}\\},var\\{\widetilde{W}_{Y,j}\\})^{1/2}}=\frac{\nu_{XY}(\mathcal{T}_{j})}{\nu_{X}(\mathcal{T}_{j})\nu_{Y}(\mathcal{T}_{j})}$
(12)

where $\nu_{X}(\mathcal{T}_{j})\;\text{and}\;\nu_{Y}(\mathcal{T}_{j})$ are the
wavelet variances for $\\{X_{t}\\}\;\text{and}\;\\{Y_{t}\\}$, (see Percival
and Walden, 2000, Chap 8).

## References

  * Achard (2012) S. Achard.  _brainwaver: Basic wavelet analysis of multivariate time series with a visualisation and parametrisation using graph theory_ , 2012.  URL https://CRAN.R-project.org/package=brainwaver.  R package version 1.6.
  * Ashok and Mahalakshmi (2017) S. Ashok and P. Mahalakshmi.  Wavelet-based feature extraction for classification of epileptic seizure EEG signal.  _Journal of Medical Engineering and Technology_ , 41:1–11, 11 2017.  doi: 10.1080/03091902.2017.1394388.
  * Chen et al. (2017) D. Chen, S. Wan, J. Xiang, and F. Bao.  A high-performance seizure detection algorithm based on discrete wavelet transform (DWT) and EEG.  _PLOS ONE_ , 12:e0173138, 03 2017.  doi: 10.1371/journal.pone.0173138.
  * Daubechies (1990) I. Daubechies.  The wavelet transform, time-frequency localization and signal analysis.  _IEEE Trans. Inf. Theory_ , 36:961–1005, 1990.
  * Elsa Jacob et al. (2018) J. Elsa Jacob, G. Nair, T. Iype, and A. Cherian.  Diagnosis of encephalopathy based on energies of EEG subbands using discrete wavelet transform and support vector machine.  _Neurology Research International_ , 2018:1–9, 07 2018\.  doi: 10.1155/2018/1613456.
  * Faust et al. (2015) O. Faust, U. R. Acharya, H. Adeli, and A. Adeli.  Wavelet-based EEG processing for computer-aided seizure detection and epilepsy diagnosis.  _Seizure_ , 26, 01 2015.  doi: 10.1016/j.seizure.2015.01.012.
  * Goldberger et al. (2000) A. L. Goldberger, L. A. N. Amaral, L. Glass, J. M. Hausdorff, P. Ch. Ivanov, R. G. Mark, J. E. Mietus, G. B. Moody, C.-K. Peng, and H. E. Stanley.  PhysioBank, PhysioToolkit, and PhysioNet: Components of a new research resource for complex physiologic signals.  _Circulation_ , 101(23):e215–e220, 2000.  Circulation Electronic Pages: http://circ.ahajournals.org/content/101/23/e215.full PMID:1085218; doi: 10.1161/01.CIR.101.23.e215.
  * Grant and Islam (2019) P. Grant and M. Z. Islam.  Clustering noisy temporal data.  In _15th International Conference on Advanced Data Mining and Applications (ADMA 2019)_ , Lecture Notes in Computer Science, pages 185–194. Springer, November 2019.  ISBN 9783030352301.  doi: https://doi.org/10.1007/978-3-030-35231-8˙13.
  * Hall et al. (2009) M. Hall, E. Frank, G. Holmes, B. Pfahringer, P. Reutemann, and I. H. Witten.  The WEKA data mining software: an update.  _SIGKDD Explorations_ , 11(1):10–18, 2009.
  * Helwig (2018) N. E. Helwig.  _eegkit: Toolkit for Electroencephalography Data_ , 2018.  URL https://CRAN.R-project.org/package=eegkit.  R package version 1.0-4.
  * Islam and Giggins (2012) M. Z. Islam and H. Giggins.  Knowledge discovery through SysFor: A systematically developed forest of multiple decision trees.  _Conferences in Research and Practice in Information Technology Series_ , 121, 01 2012.
  * Janjarasjitt (2010) S. Janjarasjitt.  Classification of the epileptic EEGs using the wavelet-based scale variance feature.  _Int J Appl Biomed Eng_ , 3:19–25, 01 2010.
  * Lee et al. (2014) S Lee, J Lim, J Kim, J Yang, and Y Lee.  Classification of normal and epileptic seizure eeg signals using wavelet transform, phase-space reconstruction, and euclidean distance.  _Comput Methods Programs Biomed_ , 116(1):10–25, 08 2014\.  doi: doi:10.1016/j.cmpb.2014.04.012.
  * Omerhodzic et al. (2013) I. Omerhodzic, S. Avdakovic, A. Nuhanovic, and K. Dizdarevic.  Energy distribution of EEG signals: EEG signal wavelet-neural network classifier.  _Int. J. Biol. Life Sci._ , 6, 07 2013.
  * Parvez and Paul (2013) Md Z. Parvez and M. Paul.  EEG signal classification using frequency band analysis towards epileptic seizure prediction.  _16th Int’l Conf. Computer and Information Technology, ICCIT 2013_ , 12 2013.  doi: 10.1109/ICCITechn.2014.6997315.
  * Percival and Walden (2000) D. Percival and A. Walden.  _Wavelet Methods for Time Series Analysis_.  Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, 2000.  doi: 10.1017/CBO9780511841040.
  * Polikar (1994) R. Polikar.  The Wavelet Tutorial, The engineer’s ultimate guide to wavelet analysis, 1994.  URL users.rowan.edu/~polikar/WTtutorial.html.
  * Qu and Gotman (1997) H. Qu and J. A. Gotman.  A patient-specific algorithm for the detection of seizure onset in long-term eeg monitoring: Possible use as a warning devices.  _IEEE Trans Biomed Eng_ , IT-44(2):115, 1997.
  * R Core Team (2019) R Core Team.  _R: A Language and Environment for Statistical Computing_.  R Foundation for Statistical Computing, Vienna, Austria, 2019.  URL https://www.R-project.org/.
  * Shrestha et al. (2019) S. Shrestha, R. Shrestha, and B. Thapa.  Implementing neural network and multi resolution analysis in EEG signal for early detection of epilepsy.  _SCITECH Nepal_ , 14:8–16, 09 2019.  doi: 10.3126/scitech.v14i1.25528.
  * Siddiqui et al. (2018) Md K. Siddiqui, M. Z. Islam, and A. Kabir.  A novel quick seizure detection and localization through brain data mining on ECoG dataset.  _Neural Computing and Applications_ , pages 1–14, 03 2018.  doi: 10.1007/s00521-018-3381-9.
  * Siers and Islam (2015) M. Siers and M. Z. Islam.  Software defect prediction using a cost sensitive decision forest and voting, and a potential solution to the class imbalance problem.  _Information Systems_ , 51:62–71, 07 2015.  doi: 10.1016/j.is.2015.02.006.
  * Subasi (2007) A. Subasi.  EEG signal classification using wavelet feature extraction and a mixture of expert model. expert systems with applications 32, 1084-1093.  _Expert Systems with Applications_ , 32:1084–1093, 05 2007\.  doi: 10.1016/j.eswa.2006.02.005.
  * Torgo (2014) L. Torgo.  An infra-structure for performance estimation and experimental comparison of predictive models in R.  _CoRR_ , abs/1412.0436, 2014.  URL http://arxiv.org/abs/1412.0436.
  * Vikash (2015) R. Vikash.  Comparison of DFT and wavelet based image modification techniques.  _International Journal of Computer Science and Mobile Computing_ , 4:61–65, 8 2015.  URL https://www.ijcsmc.com/docs/papers/August2015/V4I8201524.pdf.
  * Zandi et al. (2010) A. Zandi, M. Javidan, G. Dumont, and R. Tafrershi.  Automated real-time epileptic seizure detection in scalp EEG recordings using an algorithm based on wavelet packet transform.  _IEEE transactions on bio-medical engineering_ , 57:1639–51, 07 2010.  doi: 10.1109/TBME.2010.2046417.


# Stochastic Local Winner-Takes-All Networks Enable Profound Adversarial
Robustness

Konstantinos P. Panousis
Cyprus University of Technology
Limassol, Cyprus
<EMAIL_ADDRESS>
&Sotirios Chatzis
Cyprus University of Technology
Limassol, Cyprus
<EMAIL_ADDRESS>
Sergios Theodoridis
National and Kapodistrian University of Athens, Greece
Aalborg University, Aalborg, Denmark
<EMAIL_ADDRESS>
###### Abstract

This work explores the potency of stochastic competition-based activations,
namely Stochastic Local Winner-Takes-All (LWTA), against powerful (gradient-
based) white-box and black-box adversarial attacks; we especially focus on
Adversarial Training settings. In our work, we replace the conventional ReLU-
based nonlinearities with blocks comprising locally and stochastically
competing linear units. The output of each network layer now yields a sparse
output, depending on the outcome of winner sampling in each block. We rely on
the Variational Bayesian framework for training and inference; we incorporate
conventional PGD-based adversarial training arguments to increase the overall
adversarial robustness. As we experimentally show, the arising networks yield
state-of-the-art robustness against powerful adversarial attacks while
retaining very high classification rate in the benign case.

##  1 Introduction

In this paper, we revisit the novel stochastic formulation of deep networks
with LWTA activations (Panousis et al.,, 2019, 2021), and delve deeper into
its potency against adversarial attacks in the context of PGD-based
Adversarial Training (AT) (Madry et al.,, 2017). We evaluate the robustness of
the emerging networks against powerful gradient-based white-box as well as
black-box adversarial attacks using the well-known AutoAttack (AA) framework
(Croce and Hein,, 2020). We provide the related source code at:
https://github.com/konpanousis/Adversarial-LWTA-AutoAttack. We experimentally
show that Stochastic LWTA-based networks not only yield state-of-the-art
accuracy against all the considered attacks, but do so while retaining very
high classification accuracy in the benign case.

##  2 Stochastic Local Winner-Takes-All (LWTA) Networks

Let us consider an input $\boldsymbol{x}\in\mathbb{R}^{J}$ presented to a
conventional deep neural network layer comprising $K$ hidden units and weight
matrix $\boldsymbol{W}\in\mathbb{R}^{J\times K}$. Each hidden unit $k$ within
the layer performs an inner product computation
$h_{k}=\boldsymbol{w}_{k}^{T}\boldsymbol{x}=\sum_{j=1}^{J}w_{jk}\cdot
x_{j}\in\mathbb{R}$; this (usually) passes through a non-linear function
$\sigma(\cdot)$. Thus, the final layer output vector
$\boldsymbol{y}\in\mathbb{R}^{K}$ is formed via the concatenation of the non-
linear activation of each individual hidden unit, such that
$\boldsymbol{y}=[y_{1},\dots,y_{K}]$, where $y_{k}=\sigma(h_{k})$.

In contrast, in a fully connected LWTA-based network layer, singular nonlinear
units are replaced by $U$ linear competing units aggregated together in a
(LWTA) block; in the following, we denote with $B$ the number of LWTA blocks
in an LWTA-based layer. The associated weights are now arranged as a three
dimensional matrix $\boldsymbol{W}\in\mathbb{R}^{J\times B\times U}$
signifying that the input $\boldsymbol{x}$ is presented to each block and each
unit therein. Specifically, in the LWTA-based framework, the $u$th competing
unit within the $b$th block computes its activation $h_{b,u}$ via the standard
inner product computation
$h_{b,u}=\boldsymbol{w}_{b,u}^{T}\boldsymbol{x}=\sum_{j=1}^{J}w_{j,b,u}\cdot
x_{j}\in\mathbb{R}$; then, competition takes place among the units in the
block.

The underlying principle is that out of the $U$ units in an LWTA block, only
one can be the winner; this unit gets to pass its linear activation to the
next layer, while the rest output zero values. Thus, the output of an LWTA
layer $\boldsymbol{y}\in\mathbb{R}^{B\cdot U}$ is composed of $B$ subvectors
$\boldsymbol{y}_{b}\in\mathbb{R}^{U}$, one for each LWTA block and each with a
single non-zero entry. It is apparent that this competition process results in
a sparse representation, since all units, except one in each block, produce a
zero output. In related literature, the competition procedure is
deterministic, i.e., the winner unit is the one with the highest activation.
However, stochastic competition principles have been recently proposed in
Panousis et al., (2019, 2021); Voskou et al., (2021).

Figure 1: A detailed bisection of the $b$th Stochastic LWTA block in an LWTA
layer. Presented with an input $\boldsymbol{x}\in\mathbb{R}^{J}$, each unit
$u=1,\dots,U$ computes its activation $h_{b,u}$ via different weights
$\boldsymbol{w}_{b,u}\in\mathbb{R}^{J}$, i.e.,
$h_{b,u}=\boldsymbol{w}_{b,u}^{T}\boldsymbol{x}$. The linear responses of the
units are concatenated, such that
$\boldsymbol{h}_{b}=[h_{b,1},\dots,h_{b,U}]$, and transformed into
probabilities via the softmax operation. Then, a Discrete sample
$\boldsymbol{\xi}_{b}=[\xi_{b,1},\dots,\xi_{b,U}]$ is drawn; this constitutes
an one-hot vector with a single non-zero entry at position $u^{\prime}$,
denoting the winner unit in the block. This winner unit, $u^{\prime}$, passes
its linear response to the next layer; the rest pass zero values.

In this context, to encode the winner unit in each of the $B$ LWTA blocks that
constitute a stochastic LWTA layer, we introduce an appropriate set of
discrete latent indicator vectors
$\boldsymbol{\xi}\in\mathrm{one\\_hot}(U)^{B}$. This vector comprises $B$
component subvectors, where each component entails exactly one non-zero value
at the index position that corresponds to the winner unit in each respective
LWTA block.

Thus, the output $\boldsymbol{y}$ of a stochastic LWTA layer’s $(b,u)$th
component $y_{b,u}$ is defined as:

$\displaystyle y_{b,u}=\xi_{b,u}\sum_{j=1}^{J}w_{j,b,u}\cdot
x_{j}\in\mathbb{R}$ (1)

where $\xi_{b,u}$ denotes the $u$th component of $\boldsymbol{\xi}_{b}$, and
$\boldsymbol{\xi}_{b}\in\mathrm{one\\_hot}(U)$ holds the $b$th subvector of
$\boldsymbol{\xi}$.

We postulate that the latent winner indicator variables
$\boldsymbol{\xi}_{b},\forall b$ in Eq.(1) are drawn from a data-driven
Categorical distribution with probabilities proportional to the intermediate
linear computations that each unit performs. Therefore, the higher the linear
response of a particular unit in a particular block, the higher its
probability of it being the winner in said block; this yields:

$\displaystyle
q(\boldsymbol{\xi}_{b})=\mathrm{Categorical}\Big{(}\boldsymbol{\xi}_{b}\Big{|}\mathrm{softmax}\Big{(}\sum_{j=1}^{J}[w_{j,b,u}]_{u=1}^{U}\cdot
x_{j}\Big{)}$ (2)

where $[w_{j,b,u}]_{u=1}^{U}$ denotes the vector concatenation of the set
$\\{w_{j,b,u}\\}_{u=1}^{U}$. A graphical illustration of the proposed
stochastic LWTA block is depicted in Fig. 1. Each stochastic LWTA layer
comprises multiple such LWTA blocks, as illustrated in Fig. 2(a). Note that
stochastically selecting the winner of each LWTA block in each layer
introduces stochasticity to the network activations. Presented with the same
input, different subnetworks may be activated and a different subpath is
followed to the output, as a result of winner sampling.

$x_{1}$$x_{J}$$\dots$$\dots$$\dots$$\dots$ $\xi$= 1 $\xi$= 0 $\xi$= 0 $\xi$= 1
LWTA layer $1$ $1$ $B$ $B$ LWTA layer Input layer Output layer

(a)

(b)

Figure 2:  (a) A graphical representation of our competition-based modeling
approach. Rectangles denote LWTA blocks, and circles the competing units
therein. The winner units are denoted with bold contours ($\xi=1$). (b) The
convolutional LWTA variant. Competition takes place position-wise among the
feature maps comprising a kernel. For each position, only the winner feature
map contains a non-zero entry; for the rest feature maps in the kernel, the
value at said position is zero.

Further, to also accommodate in our study architectures based on the
convolutional operation, we adopt the corresponding variant of the Stochastic
LWTA activation in Panousis et al., (2019). We assume an input tensor
$\boldsymbol{X}\in\mathbb{R}^{H\times L\times C}$ and define a set of kernels,
each with weights $\boldsymbol{W}_{b}\in\mathbb{R}^{h\times l\times C\times
U}$, where $h,l,C,U$ are the kernel height, length, channels and competing
feature maps, and $b=1,\dots,B$. Analogously to the grouping of linear units
in dense layers, in this case, local competition is performed among feature
maps on a _position-wise_ basis. Each kernel is treated as an LWTA block with
competing feature maps; each layer comprises $B$ kernels. Specifically, each
feature map $u=1,\dots,U$ in the $b$th LWTA block of a convolutional LWTA
layer computes:

$\displaystyle\boldsymbol{H}_{b,u}=\boldsymbol{W}_{b,u}\star\boldsymbol{X}\in\mathbb{R}^{H\times
L}$ (3)

Then, competition takes place among competing feature maps on a position-wise
basis. The competitive random sampling procedure reads:

$\displaystyle
q(\boldsymbol{\xi}_{b,h^{\prime},l^{\prime}})=\mathrm{Categorical}\left(\boldsymbol{\xi}_{b,h^{\prime},l^{\prime}}\
\Big{|}\mathrm{softmax}\left(\left[\boldsymbol{H}_{b,1,h^{\prime},l^{\prime}},\dots,\boldsymbol{H}_{b,U,h^{\prime},l^{\prime}}\right]\right)\right),\
\forall h^{\prime},l^{\prime}$ (4)

In each kernel, $b=1,\dots,B$, for each position, only the winner feature map
contains a non-zero entry; all the rest feature maps contain zero values at
this position. This yields sparse feature maps with mutually exclusive active
positions.

Thus, at a given layer of the proposed convolutional variant, the output
$\boldsymbol{Y}\in\mathbb{R}^{H\times L\times B\cdot U}$ is obtained via
concatenation of the subtensors $\boldsymbol{Y}_{b,u}$ that read:

$\displaystyle\boldsymbol{Y}_{b,u}=\boldsymbol{\Xi}_{b,u}\Big{(}\boldsymbol{W}_{b,u}\star\boldsymbol{X}\Big{)},\
\forall b,u$ (5)

where
$\boldsymbol{\Xi}_{b,u}=[\xi_{b,u,h^{\prime},l^{\prime}}]_{h^{\prime},l^{\prime}=1}^{H,L}$.

Figure 3: A detailed bisection of the $b$th convolutional stochastic LWTA
block. Presented with an input $\boldsymbol{X}\in\mathbb{R}^{H\times L\times
C}$, competition now takes place among feature maps on a position-specific
basis. Only the winner feature map contains a non-zero entry in a specific
position. This leads to sparse feature maps, each comprising uniquely
position-wise activated pixels.

The corresponding illustration of the proposed stochastic convolutional LWTA
block is depicted in Fig. 3. Convolutional Stochastic LWTA-based layers
comprise multiple such blocks as shown in Fig. 2(b).

###  2.1 Training

Since a network composed of such stochastic LWTA layers entails latent
variables $\boldsymbol{\xi}$, we resort to a Bayesian treatment to perform
effective parameter estimation. To this end, we turn to a stochastic gradient
variational Bayes treatment (Kingma and Welling,, 2014) for scalability. The
resulting objective takes the form of an evidence lower-bound (ELBO) as
described next.

Considering data $\mathcal{D}=\\{X_{i},Y_{i}\\}_{i=1}^{N}$, we define the
categorical cross-entropy between the data labels $Y_{i}$ and the class
probabilities $f(X_{i};\hat{\boldsymbol{\xi}})$, generated by the penultimate
Softmax layer of a Stochastic LWTA-based network, as
$\mathrm{CE}(Y_{i},f(X_{i};\hat{\boldsymbol{\xi}}))$. Here,
$\hat{\boldsymbol{\xi}}$ denotes sample instances of all the latent variables,
$\boldsymbol{\xi}$, in all layers. We stress that the output class
probabilities depend on the winner selection process in each layer, which is
stochastic.

This way, the Evidence Lower Bound (ELBO) reads:

$\displaystyle\mathcal{L}=-\sum_{X_{i},Y_{i}\in
D}\mathrm{CE}(Y_{i},f(X_{i};\hat{\boldsymbol{\xi}}))-\mathrm{KL}[q(\boldsymbol{\xi})||p(\boldsymbol{\xi})])$
(6)

For simplicity, and without loss of generality, we consider a symmetric
Categorical distribution for the latent variable indicators
$\boldsymbol{\xi}$; hence, $p(\boldsymbol{\xi}_{b})=\mathrm{Categorical}(1/U)\
\forall b$ for dense layers, and
$p(\boldsymbol{\xi}_{b,h^{\prime},l^{\prime}})=\mathrm{Categorical}(1/U)\
\forall b,h^{\prime},l^{\prime}$ for convolutional ones. In our work, we
perform Monte-Carlo sampling using a single reparameterized sample for each of
the corresponding latent variables. These are obtained via the
reparameterization trick of the continuous relaxation of the Categorical
distribution (Maddison et al.,, 2017; Jang et al.,, 2017) as described next.
We focus on the reparameterization trick for the dense case; the convolutional
case is analogous.

Let $\tilde{\xi}$ denote the probabilities of $q(\boldsymbol{\xi})$ (Eqs. (2)
and (4)). Then, the samples $\hat{\boldsymbol{\xi}}$ are expressed as:

$\displaystyle\hat{\xi}_{b,u}=\mathrm{Softmax}((\log\tilde{\xi}_{b,u}+g_{b,u})/\tau),\
\forall b=1,\dots,B,\ u=1,\dots,U$ (7)

where $g_{b,u}=-\log(-\log V_{b,u}),\ V_{b,u}\sim\mathrm{Uniform}(0,1)$, and
$\tau\in(0,\infty)$ is a temperature factor, controlling how “closely” the
continuous relaxation approximates the Categorical distribution.

On this basis, we can write the KL divergence term present in Eq. (6) as:

$\displaystyle\begin{split}\mathrm{KL}[q(\boldsymbol{\xi}_{b})||p(\boldsymbol{\xi}_{b})]&=\mathbb{E}_{q(\boldsymbol{\xi}_{b})}[\log
q(\boldsymbol{\xi}_{b})-\log p(\boldsymbol{\xi}_{b})]\\\ &\approx\log
q(\hat{\boldsymbol{\xi}}_{b})-\log p(\hat{\boldsymbol{\xi}}_{b}),\ \forall
b\end{split}$ (8)

Hence, the final ELBO expression yields:

$\displaystyle\mathcal{L}=-\sum_{X_{i},Y_{i}\in
D}\mathrm{CE}(Y_{i},f(X_{i};\hat{\boldsymbol{\xi}}))-\sum_{b}\left(\log
q(\hat{\boldsymbol{\xi}}_{b})-\log p(\hat{\boldsymbol{\xi}}_{b})\right)$ (9)

###  2.2 Prediction

At prediction time, we directly draw $L$ samples from the trained posteriors
$q(\boldsymbol{\xi})$ in order to determine the winning units in each block of
the network. As mentioned before, each time we sample for the same input, a
different subpath is followed. This is a key aspect that stochastically alters
the information flow in the network and obstructs an adversary from attacking
the model.

The sampling process results in a set of $L$ output logits of the network,
which we can average to obtain the final prediction:

$\displaystyle
f(X_{i};\hat{\boldsymbol{\xi}})=\frac{1}{L}\sum_{l=1}^{L}f(X_{i};\hat{\boldsymbol{\xi}}^{l})$
(10)

where $\hat{\boldsymbol{\xi}}^{l}$ denotes a sample drawn from
$q(\boldsymbol{\xi})$.

##  3 Experimental Results

We investigate the potency of LWTA-based networks against adversarial attacks
under an Adversarial Training regime; we employ a PGD adversary (Madry et
al.,, 2017). To this end, we use the well-known WideResNet-34 (Zagoruyko and
Komodakis,, 2016) architecture, considering three different widen factors: 1,
5, and 10; we focus on the CIFAR-10 dataset and adopt experimental settings
similar to Wu et al., (2021). We use a batch-size of 128 and an initial
learning rate of 0.1; we halve the learning rate at every epoch after the
$75^{th}$ epoch. We use a single sample $L=1$ for prediction. All experiments
were performed using a single NVIDIA Quadro P6000.

For evaluating the robustness of the proposed structure, we initially consider
the conventional PGD attack with 20 steps, step size $0.007$ and
$\epsilon=8/255$. In Table 1, we compare the robustness of LWTA-based
WideResNet networks against the baseline results of Wu et al., (2021). As we
observe, our Stochastic LWTA-based networks yield significant improvements in
robustness under a traditional PGD attack; they retain extremely high natural
accuracy (up to $\approx 13\%$ better), while exhibiting a staggering, up to
$\approx 32.6\%$, difference in robust accuracy compared to the exact same
architectures employing the conventional ReLU-based nonlinearities and trained
in the exact same fashion.

Table 1: Natural and Robust accuracy under a conventional PGD attack with 20
steps and $0.007$ step-size using WideResNet-34 models with different widen
factors. We use the same PGD-based Adversarial Training scheme for all models
(Madry et al.,, 2017).

Adversarial Training-PGD
---
| Natural Accuracy ($\%$) | Robust Accuracy ($\%$)
Widen Factor | Baseline | Stochastic LWTA | Baseline | Stochastic LWTA
1 | 74.04 | 87.0 | 49.24 | 81.87
5 | 83.95 | 91.88 | 54.36 | 83.4
10 | 85.41 | 92.26 | 55.78 | 84.3

Table 2: Robust Accuracy $(\%)$ comparison under the AutoAttack framework.
$\dagger$ denotes models that are trained with additional unlabeled data. The
AutoAttack performance corresponds to the final robust accuracy after
employing all the attacks in AA. Results directly from the AA leaderboard.

Method | AutoAttack
---|---
TRADES(Zhang et al.,, 2019) | 53.08
Early-Stop (Rice et al.,, 2020) | 53.42
FAT (Zhang et al.,, 2020) | 53.51
HE (Pang et al.,, 2020) | 53.74
WAR (Wu et al.,, 2021) | 54.73
Pre-training (Hendrycks et al.,, 2019)$\dagger$ | 54.92
MART (Wang et al.,, 2020)$\dagger$ | 56.29
HYDRA (Sehwag et al.,, 2020)$\dagger$ | 57.14
RST (Carmon et al.,, 2019)$\dagger$ | 59.53
Gowal et al., (2021)$\dagger$ | 65.88
WAR (Wu et al.,, 2021)$\dagger$ | 61.84
Ours (Stochastic-LWTA/PGD/WideResNet-34-1) | 74.71
Ours (Stochastic-LWTA/PGD/WideResNet-34-5) | 81.22
Ours (Stochastic-LWTA/PGD/WideResNet-34-10) | 82.60

Further, and to ensure that our approach does not cause the well-known
obfuscated gradient problem (Athalye et al.,, 2018), we resort to stronger
parameter-free attacks using the newly introduced AutoAttack (AA) framework
(Croce and Hein,, 2020). AA comprises an ensemble of four powerful white-box
and black-box attacks, e.g., the commonly employed A-PGD attack; this is a
step-free variant of the standard PGD attack (Madry et al.,, 2017), which
avoids the complexity and ambiguity of step-size selection. In addition, for
the entailed $L_{\infty}$ attack, we use the common $\epsilon=8/255$ value.
Thus, in Table 2, we compare the LWTA-based networks to several recent state-
of-the-art approaches evaluated on AA111https://github.com/fra31/auto-attack.
The reported accuracies correspond to the final reported robust accuracy of
the methods after sequentially performing all the considered AA attacks. Once
again, we observe that the proposed networks yield state-of-the-art robustness
against all SOTA methods, with an improvement of $\approx 16.72\%$, even when
compared with methods that employ substantial data augmentation to increase
robustness, e.g. Gowal et al., (2021). These results vouch for the potency of
Stochastic LWTA networks in adversarial settings.

Finally, since our considered networks consist of stochastic components, i.e.
the competitive random sampling procedure to determine the winner in each LWTA
block, the output of the classifier might change at each iteration; this
obstructs the attacker from successfully altering the final decision. To
counter such randomness in the involved computations, Croce and Hein, (2020)
combine the APGD attack with an averaging procedure of 20 computations of the
gradient at the same point. This technique is known as Expectation over
Transformation (EoT) (Athalye et al.,, 2018). Thus, we use AA combined with
EoT for further performance evaluation of the proposed LWTA-based networks.
The corresponding results are presented in Table 3. As we observe, all of the
considered networks retain state-of-the-art robustness against the powerful AA
& EoT attacks. This conspicuously supports the usefulness of Stochastic LWTA
activations towards adversarial robustness.

Table 3: Robustness against AA combined with 20 iterations of EoT. APGD-DLR corresponds to the APGD attack, using a different loss, i.e., the Difference of Logits Ratio (Croce and Hein,, 2020). Widen Factor | Nat. Acc. | APGD | APGD-DLR
---|---|---|---
1 | 87.00 | 79.67 | 76.15
5 | 91.88 | 81.67 | 77.65
10 | 92.26 | 82.55 | 79.00

##  4 Conclusions

In this work, we explored the potency of Stochastic LWTA-based networks
against powerful white-box and black-box attacks. The experimental results
vouch for the efficacy of the arising networks, yielding state-of-the-art
robustness in all the experimental settings. We obtained an immense
improvement in robustness compared to the second best performing alternative,
which notably relies on substantial data augmentation. A potentially key
principle towards adversarial robustness may be the stochastic alteration of
the information flow in Stochastic LWTA-based networks; this, arises from the
considered data-driven winner selection mechanism in each LWTA block.
Different subpaths, stochastically emerging even for the same input,
essentially obstruct the adversary from successfully attacking the model.
Further evaluation against stochasticity countermeasures, i.e., Expectation
over Transformation (EoT), further validate our findings, as the induced
decrease in the final robustness was negligible.

## Acknowledgements

This work has received funding from the European Union’s Horizon 2020 research
and innovation program under grant agreement No 872139, project aiD.

## References

  * Athalye et al., (2018) Athalye, A., Carlini, N., and Wagner, D. A. (2018).  Obfuscated gradients give a false sense of security: Circumventing defenses to adversarial examples.  In Proc. ICML.
  * Carmon et al., (2019) Carmon, Y., Raghunathan, A., Schmidt, L., Duchi, J. C., and Liang, P. S. (2019).  Unlabeled data improves adversarial robustness.  In Proc. NIPS.
  * Croce and Hein, (2020) Croce, F. and Hein, M. (2020).  Reliable evaluation of adversarial robustness with an ensemble of diverse parameter-free attacks.  In Proc. ICML.
  * Gowal et al., (2021) Gowal, S., Qin, C., Uesato, J., Mann, T., and Kohli, P. (2021).  Uncovering the limits of adversarial training against norm-bounded adversarial examples.  arXiv preprint arXiv:2010.03593.
  * Hendrycks et al., (2019) Hendrycks, D., Lee, K., and Mazeika, M. (2019).  Using pre-training can improve model robustness and uncertainty.  In Proc. ICML.
  * Jang et al., (2017) Jang, E., Gu, S., and Poole, B. (2017).  Categorical reparametrization with gumbel-softmax.  In Proc. ICLR.
  * Kingma and Welling, (2014) Kingma, D. P. and Welling, M. (2014).  Auto-encoding variational bayes.  In Proc. ICLR.
  * Maddison et al., (2017) Maddison, C. J., Mnih, A., and Teh, Y. W. (2017).  The concrete distribution: A continuous relaxation of discrete random variables.  In Proc. ICLR.
  * Madry et al., (2017) Madry, A., Makelov, A., Schmidt, L., Tsipras, D., and Vladu, A. (2017).  Towards deep learning models resistant to adversarial attacks.  arXiv preprint arXiv:1706.06083.
  * Pang et al., (2020) Pang, T., Yang, X., Dong, Y., Xu, K., Zhu, J., and Su, H. (2020).  Boosting adversarial training with hypersphere embedding.  In Proc. NIPS.
  * Panousis et al., (2021) Panousis, K., Chatzis, S., Alexos, A., and Theodoridis, S. (2021).  Local competition and stochasticity for adversarial robustness in deep learning.  In Proc. AISTATS.
  * Panousis et al., (2019) Panousis, K., Chatzis, S., and Theodoridis, S. (2019).  Nonparametric Bayesian deep networks with local competition.  In Proc. ICML.
  * Rice et al., (2020) Rice, L., Wong, E., and Kolter, J. Z. (2020).  Overfitting in adversarially robust deep learning.  arXiv preprint arXiv:2002.11569.
  * Sehwag et al., (2020) Sehwag, V., Wang, S., Mittal, P., and Jana, S. (2020).  Hydra: Pruning adversarially robust neural networks.  In Proc. NIPS.
  * Voskou et al., (2021) Voskou, A., Panousis, K. P., Kosmopoulos, D., Metaxas, D. N., and Chatzis, S. (2021).  Stochastic transformer networks with linear competing units: Application to end-to-end sl translation.  In In Proc. ICCV.
  * Wang et al., (2020) Wang, Y., Zou, D., Yi, J., Bailey, J., Ma, X., and Gu, Q. (2020).  Improving adversarial robustness requires revisiting misclassified examples.  In Proc. ICLR.
  * Wu et al., (2021) Wu, B., Chen, J., Cai, D., He, X., and Gu, Q. (2021).  Do wider neural networks really help adversarial robustness?  In In Proc. NIPS.
  * Zagoruyko and Komodakis, (2016) Zagoruyko, S. and Komodakis, N. (2016).  Wide residual networks.  In Proc. BMVC.
  * Zhang et al., (2019) Zhang, H., Chen, H., Song, Z., Boning, D., Dhillon, I. S., and Hsieh, C.-J. (2019).  The limitations of adversarial training and the blind-spot attack.  arXiv preprint arXiv:1901.04684.
  * Zhang et al., (2020) Zhang, J., Xu, X., Han, B., Niu, G., Cui, L., Sugiyama, M., and Kankanhalli, M. (2020).  Attacks which do not kill training make adversarial learning stronger.  In Proc. ICML.