# WaveCatBoost for Probabilistic Forecasting of Regional Air Quality Data
††thanks: This work was supported in part by ASEAN-India Science and
Technology Collaboration (AISTIC) Funding from the Department of Science and
Technology, India (CRD/2020/000320) for providing incentives. ††thanks: J.
Borah (email<EMAIL_ADDRESS>& S. Majumdar (email<EMAIL_ADDRESS>is
with the Dept. of ECE, NIT Meghalaya, India ††thanks: T. Chakraborty (email:
<EMAIL_ADDRESS>is with the Sorbonne University, Abu Dhabi, UAE and
Sorbonne Center for Artificial Intelligence, Paris, France ††thanks: M.S.M.
Nadzir is with Dept. of Earth Sciences and Environment, Universiti Kebangsaan,
Selangor, Malaysia. (e-mail: shahrulnadzir@ukm.edu.my). ††thanks: M.G.
Cayetano is with the Institute of Environmental Science and Meteorology,
University of Philippines, Manila, Philippines (e-mail:
mcayetano@iesm.upd.edu.ph). ††thanks: J. Borah and T. Chakraborty (Joint first
authors) have equal contributions.

Jintu Borah, Tanujit Chakraborty, Md. Shahrul Md. Nadzir, Mylene G. Cayetano,
and Shubhankar Majumdar

###### Abstract

Accurate and reliable air quality forecasting is essential for protecting
public health, sustainable development, pollution control, and enhanced urban
planning. This letter presents a novel WaveCatBoost architecture designed to
forecast the real-time concentrations of air pollutants by combining the
maximal overlapping discrete wavelet transform (MODWT) with the CatBoost
model. This hybrid approach efficiently transforms time series into high-
frequency and low-frequency components, thereby extracting signal from noise
and improving prediction accuracy and robustness. Evaluation of two distinct
regional datasets, from the Central Air Pollution Control Board (CPCB) sensor
network and a low-cost air quality sensor system (LAQS), underscores the
superior performance of our proposed methodology in real-time forecasting
compared to the state-of-the-art statistical and deep learning architectures.
Moreover, we employ a conformal prediction strategy to provide probabilistic
bands with our forecasts.

###### Index Terms:

Air quality, CatBoost, Wavelet analysis, Conformal prediction, Real-time
forecasting.

##  I Introduction

Air pollution is a critical global concern, adversely affecting public health
and the environment. The rapid expansion of industrialization and urbanization
have surged the emission of air pollutants, resulting in significant
challenges to sustainable living. In response to these concerns, the World
Health Organization has issued guidelines for six major air pollutants and set
National Ambient Air Quality Standards. These pollutants, including nitrogen
dioxide ($NO_{2}$), ozone ($O_{3}$), carbon monoxide ($CO$), sulfur dioxide
($SO_{2}$), particulate matter with diameters of 2.5 mm or less ($PM_{2.5}$),
and 10 mm or less ($PM_{10}$), originate from various sources such as
transportation, industry, and natural processes. Despite considerable efforts
to reduce emissions and lower ambient concentrations of these pollutants, air
pollution remains a significant cause of mortality worldwide.

To mitigate these concerns related to polluted air, several countries have
implemented real-time air quality forecasting systems [1]. Accurate prediction
of air pollutant concentration levels is paramount for public awareness
campaigns, environmental management, and healthcare interventions, among many
others. By leveraging these forecasts, government regulations and public
policies can be designed to address pollution-related health concerns and
promote sustainable development initiatives. In recent years, machine learning
algorithms have emerged as powerful tools for enhancing the accuracy of air
quality predictions [2]. These computational techniques, driven by their
ability to decipher complex patterns and relationships within vast datasets,
offer a promising avenue for advancing pollutant forecasting models [3, 4].

However, existing methods encounter challenges adapting to diverse monitoring
environments and accurately predicting pollutant concentrations in real-time
[5]. Traditional models often struggle with the dynamic nature of air quality,
hindered by limitations in handling non-linear relationships and non-
stationary variations in pollutant dynamics. To address these challenges in
the forecasting methodologies, wavelet-based transformation is widely applied
in other applied domains, including epidemiology [6], economics [7], and
others. These wavelet-based forecasting approaches transform non-stationary
seasonal data into a series of uncorrelated distinct components. Subsequently,
the component series are modeled with various statistical and deep learning
[8] architectures to solve diverse problems in time series forecasting. While
the accuracies of these wavelet-based deep learning architectures surpass
traditional approaches for ultra-long time series, their training and
prediction times obstruct the generation of real-time forecasts. To mitigate
this challenge, decision tree-based boosting architectures are favored for
their faster convergence time in the presence of large datasets [9]. Several
boosting algorithms, including random forest, extreme gradient boosting, and
light gradient boosting machines, offer higher accuracy compared to deep
learning methods, but they overlook the sequential nature of the time series,
resulting in target leakage [10]. The recently developed CatBoost algorithm
employs an ordered boosting methodology to overcome this issue [11].

Motivated by this background, in this letter, we develop a novel wavelet-based
CatBoost (WaveCatBoost) model for generating accurate real-time forecasts of
air pollutant concentration levels to mitigate the overarching challenges of
reliable forecasting of air quality data. This integration aims to leverage
the strengths of both techniques, providing a robust and versatile model
capable of delivering reliable point and probabilistic forecasts across a
range of air pollutants via deploying a conformal prediction approach. By
conducting an extensive experimental evaluation and comparing existing models,
we demonstrate the superiority of our proposed approach in advancing air
quality forecasting and addressing the critical gaps present in the current
literature.

Figure 1: Pipeline of the proposed air quality forecasting system Figure 2:
Depiction of the forecast horizons considered in this study

##  II Material and Methodology

###  II-A Data Collection and Preprocessing

Real-time data on air pollutant concentration levels are gathered through two
wireless sensors, the Central Air Pollution Control Board (CPCB) sensor
network and a low-cost air quality sensor system named ID1, strategically
positioned in Meghalaya, India. These air quality sensors continuously monitor
the concentrations of NO2 (in ppb), O3 (in ppb), CO (in ppb), SO2 (in ppb),
PM2.5 (in $\mu$g/m$3$), and PM10 (in $\mu$g/m$3$) and transmit the data to the
server at one-minute intervals. The transmitted information is sequentially
stored in a tabular data format within the database for further analysis. We
extract the raw data from the server database for a year, from July 1, 2022,
to July 31, 2023. The preprocessing phase of our experiment involves adjusting
the missing data entries and computing hourly averages to generate near-real-
time data on the concentration of various air pollutants. Additionally, we
normalize the dataset using the min-max normalization strategy, which is
crucial for optimizing the convergence of data-driven forecasting algorithms
used in our analysis. A visual representation outlining the processes of data
collection, pre-processing, and analysis employed in this study is depicted in
Fig. 1.

###  II-B Proposed WaveCatboost

The key concept behind the wavelet-based CatBoost framework, which we named
WaveCatBoost, involves employing the scale-invariant maximal overlapping
discrete wavelet transformation (MODWT) [12] to denoise the air pollutant
concentrations series and modeling them using an ensemble of CatBoost models.
The MODWT decomposition is applied to filter the normalized air pollutant data
$(Y_{t})_{t=1}^{N}$ using a sequence of wavelet
$\left\\{p_{m};m=0,1,\ldots,\mathcal{M}-1\right\\}$ and scaling
$\left\\{q_{m};m=0,1,\ldots,\mathcal{M}-1\right\\}$ filters. These filters
adhere to the even-length scaling assumption,
$\sum_{m=0}^{\mathcal{M}-1}p_{m}p_{m+2n}=\sum_{m=0}^{M-1}q_{m}q_{m+2n}=0$ for
any non-zero integers $n$. Consequently, applying the pyramid algorithm [12]
generates wavelet and scaling coefficients for the transformed time series,
allowing the original series $Y_{t}$ to be represented as a set of equal-sized
uncorrelated time series:

$Y_{t}=\sum_{k=1}^{K}D_{k,t}+S_{K,t}\quad t=1,2,\ldots,N,$

where $K=\log_{e}N$ is the number of decomposition levels applied,
$D_{k,t}(k=1,2,\ldots,K)$ is the $k^{th}$ details series of $Y_{t}$, captures
the local fluctuations of the series, and the smooth series $S_{K,t}$
preserves the low-frequency trend of $Y_{t}$. This iterative decomposition
essentially helps filter the true signal from the noisy data and aids in
modeling long-term dependent, non-stationary, and seasonal data by breaking it
into several lower-resolution components. Subsequently, the task of generating
$h$-step ahead forecasts of $Y_{t}$ based on its historical data can be solved
by forecasting these low-resolution components based on their corresponding
prior observation. In our proposed WaveCatBoost architecture, we generate the
forecasts for each of the wavelet and scaling coefficients of $Y_{t}$ using
the gradient boosting CatBoost [11] algorithm and ensemble them using an
inverse MODWT (IMODWT) [12] approach as:

$\hat{Y}_{N+h}=\operatorname{IMODWT}\left(\hat{D}_{1,N+h},\ldots,\hat{D}_{K,N+h},\hat{S}_{K,N+h}\right),$
(1)

$\text{with
}\hat{D}_{k,N+h}=f\left(D_{k,1},D_{k,2},\ldots,D_{k,N}\right);k=1,2,\ldots,K$
and $\hat{S}_{K,N+h}=f\left(S_{K,1},S_{K,2},\ldots,S_{K,N}\right),$ where the
function $f$ indicates the CatBoost model applied on each of the $p$ lagged
values of details and smooth series. The CatBoost framework operates by
iteratively combining weak symmetric decision trees into a strong regressor.
Unlike traditional gradient boosting methods, it is prone to prediction shifts
from target leakage. CatBoost architecture employs an ordered boosting
mechanism. During the learning process, the training data undergoes a series
of time-dependent permutations, denoted as
$\left(\sigma_{0},\sigma_{1},\ldots,\sigma_{u}\right)$. The initial
permutation $\sigma_{0}$ determines the leaf values, while the subsequent $u$
permutations decide the structure of the weak learners. At each iteration $v$,
the algorithm considers a randomly chosen permutation $\sigma_{r}$. The
supporting models $M_{r,j}$ are then fitted to the preceding $j$ observations
of $\sigma_{r}$, generating the output $M_{r,j}(i)$ for the $i^{th}$ instance
in $\sigma_{r}$. To construct the tree $T_{v}$ at the current iteration,
gradients for the $i^{th}$ time-index are computed by considering previous
time steps in the $\sigma_{r}$ permutation, i.e.,
$\operatorname{grad}_{r,\sigma(i)-1}(i)$ where

$\operatorname{grad}_{r,j}(i)=\left.\frac{\partial
L\left(y_{i},A\right)}{\partial A}\right|_{A=M_{r,j}(i)}.$ (2)

The leaf value for the $i^{th}$ instance, used for evaluating candidate
splits, is obtained by averaging the corresponding gradients
$\operatorname{grad}_{r,\sigma_{r}(i)-1}$ of preceding examples belonging to
the same leaf node $\operatorname{Leaf}_{r\left(i\right)}$. Once the tree
$T_{v}$ is constructed, it is used for boosting all the supporting models
$M_{r^{\prime},j}$ with varied sets of leaf values depending on the first $j$
instances in the $\sigma_{r^{\prime}}$ permutation. When all the weak learners
are constructed, the leaf value of the CatBoost model is calculated by the
standard gradient boosting procedure. For generating the forecasts
$\hat{D}_{k,N+h};k=1,2,\ldots,K$ and $\hat{S}_{K,N+h}$, the leaf values are
computed based on target statistics learned from the CatBoost architecture
using lagged inputs. The algorithmic structure of the proposed WaveCatBoost
model is given in Algorithm 1.

###  II-C Probabilistic Forecasting using Conformal Pediction

Along with the point forecasts generated by the WaveCatBoost model, we focus
on quantifying the uncertainty associated with these forecasts using the
conformal prediction approach [13]. This non-parametric procedure generates a
probabilistic band around the point forecasts based on a conformal score
($CS_{t}$). To calculate $CS_{t}$ corresponding to $Y_{t}$ [14], the
architecture considers $p$ lagged values of the target series ($Y_{t-p}$) and
applies the WaveCatBoost and an uncertainty model $\mathbf{\hat{U}}$ such
that,
$CS_{t}=\frac{\left|Y_{t}-\operatorname{WaveCatBoost}\left(Y_{t-p}\right)\right|}{\mathbf{\hat{U}}\left(Y_{t-p}\right)}.$
Following the sequential nature of $Y_{t}$ and the $CS_{t}$, the conformal
quantile can be computed with a weighted aggregation method having a fixed
window $w_{t}=\mathbb{1}\left(t^{\prime}\geq
t-\kappa\right),\forall\;t^{\prime}<t$ of length $\kappa$ as

$CQ_{t}=\operatorname{inf}\left\\{q:\frac{1}{\operatorname{min\left(\kappa,t^{\prime}-1\right)}}\sum_{t^{\prime}=1}^{t-1}CS_{t^{\prime}}w_{t}\geq
1-\alpha\right\\}.$

The conformal prediction intervals based on these weighted conformal quantiles
can be estimated as

$\left[\operatorname{WaveCatBoost}\left(Y_{t-p}\right)\pm
CQ_{t}\mathbf{\hat{U}}\left(Y_{t-p}\right)\right].$ (3)

Algorithm 1 WaveCatBoost Model

0: Air pollutant concentration data $\left\\{Y_{t}\right\\}_{t=1}^{N}$

0: $h$-step ahead forecast
$\left\\{\hat{Y}_{N+1},\ldots,\hat{Y}_{N+h}\right\\}$

1: Initialize the MODWT algorithm with Haar filter and $K=\log_{e}N$ levels.

2: Transform the training data $\left\\{Y_{t}\right\\}$ using the MODWT
approach into $K$ uncorrelated details series
$\left\\{D_{k,t}\right\\}_{k=1}^{K}$ and a smooth series
$\left\\{S_{K,t}\right\\}$

3: for $k=1,2,\ldots,K$ do

4: Apply CatBoost model to the $k^{th}$ detail series $D_{k,t}$.

5: $\sigma_{0},\sigma_{1},\ldots,\sigma_{u}\leftarrow$ Time-dependent
permutations of $D_{k,t}$.

6: $\sigma_{0}$ decides the leaf nodes and the next $u$ permutations determine
the structure of the symmetric weak learners.

7: for $v=1,2,\ldots,V$ do

8: $\sigma_{r}\longleftarrow\;\text{random permutations
of}\;\left[\sigma_{0},\sigma_{1},\ldots,\sigma_{u}\right]$.

9: Fit $M_{r,j}$ to $j$ preceeding observations of $\sigma_{r}$ and generate
output $M_{r,j}(i)$ for $i^{th}$ instance of $\sigma_{r}$.

10: Calculate leaf values of the week learner $T_{v}$ by averaging the
gradients $\operatorname{grad}_{r,\sigma_{r}(i)-1}$ (Eq. 2) of previous
examples belonging to the same node $\operatorname{Leaf}_{r\left(i\right)}$.

11: Boost the models $M_{r^{\prime},j}$ of
$\sigma_{r^{\prime}};\;r^{\prime}\neq r$ using $T_{v}$.

12: end for

13: Calculate leaf value for the CatBoost model by applying gradient boosting
approach to $V$ week learners’ leaf nodes and generate the one-step ahead
forecast for $D_{k,t}$.

14: Iteratively generate multi-step ahead forecast.

15: end for

16: Similarly, utilize a CatBoost architecture on $S_{K,t}$ to generate the
corresponding one-step ahead and multi-step ahead forecast iteratively.

17: Apply Inverse MODWT transform on the component forecasts from the ensemble
of CatBoost models and generate the final forecasts of the desired horizon as
in Eq. 1.

##  III Experimental Results

In this section, we present experimental results showcasing the efficacy of
our proposed WaveCatBoost framework compared to benchmark forecasting methods.
We forecast hourly air pollution levels across different horizons through
subsequent experiments. It helps identify each model’s temporal sensitivity,
guiding the selection of the most suitable model for specific prediction
needs. Additionally, it ensures the models apply to a range of real-world
scenarios, from immediate health advisories to long-term urban planning. This
analysis informs resource allocation, enabling tailored responses for short-
term and long-term forecasts, and supports the construction of ensemble models
for enhanced accuracy. Furthermore, it deepens our scientific understanding of
air quality dynamics and aids in crafting effective policies and mitigation
strategies, addressing different timeframes, from rapid interventions to long-
term sustainability planning. Our analysis utilizes a rolling window
forecasting technique, assessing model performance across four distinct train-
test scenarios: with forecast horizons set at 1 day (24 test points), 7 days,
14 days, and 31 days. Fig. 2 presents a graphical representation of these
train-test configurations.

In our analysis we compare the performance of the proposed WaveCatBoost
architecture with state-of-the-art forecasting methodologies including
NLinear, DLinear, light gradient-boosting machines (LGBM), extreme gradient-
boosting (XGB), NBeats, NHiTS, Transformer, temporal convolution network
(TCN), recurrent neural network (RNN), gated recurrent unit (GRU), long-short-
term memory (LSTM), temporal fusion transformer (TFT), and CatBoost model [5].
We measure the forecast accuracy of various frameworks using the mean absolute
scaled error (MASE) metric, one of the most recommended performance measures
in forecasting literature, on the test data as follows:

$\text{MASE}=\frac{N-S}{h}\frac{\sum_{t=N+1}^{N+h}\left|\hat{Y}_{t}-Y_{t}\right|}{\sum_{t=S+1}^{N}\left|Y_{t}-Y_{t-S}\right|},$

where $h$ is the forecast horizon, $N$ is the length of training data,
$\hat{Y}_{t}$ is the forecast of $Y_{t}$ at time $t$ and $S$ is the
seasonality of the data.

(a)

(b)

Figure 3: MCB plots for MASE metric of the various forecasting models across
all forecast horizons utilizing two different sensors (a) ID1, (b) CPCB

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4: Actual vs Forecasts of the targeted pollutants for ID1 Sensor for 1
day forecast horizon (a) $NO_{2}$, (b) $O_{3}$, (c) $CO$, (d) $SO_{2}$, (e)
$PM_{2.5}$, and (f) $PM_{10}$

(a)

(b)

(c)

(d)

(e)

(f)

Figure 5: Actual vs Forecasts of the targeted pollutants for ID1 Sensor for 7
day forecast horizon (a) $NO_{2}$, (b) $O_{3}$, (c) $CO$, (d) $SO_{2}$, (e)
$PM_{2.5}$, and (f) $PM_{10}$ TABLE I: Performance of the proposed
WaveCatBoost model and selected baseline forecasters for CPCB and ID1 sensors,
evaluated based on the MASE metric for different forecast horizons.

Pollutant | Horizon | CPCB Sensor | ID1 Sensor
---|---|---|---
RNN | LGBM | TCN | TFT | Transformer | XGB | CatBoost | WaveCatBoost | RNN | LGBM | TCN | TFT | Transformer | XGB | CatBoost | WaveCatBoost
$NO_{2}$ | 1d | 2.45 | 0.79 | 1.45 | 1.41 | 3.28 | 0.91 | 0.82 | 0.57 | 2.05 | 1.27 | 2.15 | 1.62 | 2.07 | 1.40 | 1.20 | 0.91
7d | 1.35 | 0.66 | 0.67 | 2.65 | 0.93 | 0.66 | 1.04 | 0.84 | 1.66 | 0.82 | 1.64 | 1.76 | 1.76 | 0.95 | 0.84 | 0.52
14d | 0.81 | 0.71 | 1.79 | 4.27 | 2.05 | 0.72 | 0.73 | 0.72 | 1.57 | 0.85 | 1.56 | 2.06 | 1.63 | 0.92 | 1.69 | 0.42
31d | 1.98 | 0.75 | 1.87 | 1.18 | 1.35 | 0.76 | 1.17 | 0.47 | 1.47 | 1.15 | 1.48 | 1.95 | 1.86 | 1.14 | 1.36 | 0.40
$O_{3}$ | 1d | 5.19 | 1.87 | 1.85 | 2.43 | 11.47 | 1.52 | 2.69 | 0.98 | 3.65 | 1.61 | 3.32 | 2.79 | 3.75 | 1.72 | 1.49 | 1.42
7d | 5.46 | 1.68 | 4.66 | 6.16 | 5.82 | 1.51 | 4.84 | 0.76 | 3.57 | 1.64 | 2.55 | 2.84 | 3.18 | 1.96 | 1.79 | 0.72
14d | 6.24 | 1.23 | 4.38 | 14.68 | 9.54 | 1.27 | 1.77 | 1.11 | 3.85 | 1.80 | 3.67 | 3.74 | 3.12 | 2.06 | 2.32 | 0.64
31d | 3.56 | 2.46 | 4.26 | 2.71 | 3.56 | 2.44 | 5.37 | 0.25 | 3.59 | 2.72 | 3.63 | 3.56 | 3.61 | 2.87 | 2.66 | 1.11
$CO$ | 1d | 2.36 | 3.36 | 2.90 | 2.45 | 0.97 | 3.21 | 2.87 | 2.69 | 4.47 | 2.48 | 3.82 | 2.95 | 4.47 | 2.54 | 2.32 | 1.53
7d | 3.74 | 3.62 | 3.94 | 5.29 | 3.74 | 3.59 | 3.57 | 2.61 | 2.85 | 1.52 | 2.64 | 2.95 | 2.84 | 1.46 | 1.36 | 0.79
14d | 2.76 | 3.06 | 2.79 | 4.34 | 2.77 | 2.88 | 3.15 | 1.63 | 3.08 | 1.93 | 3.08 | 3.71 | 3.06 | 1.83 | 3.36 | 0.83
31d | 1.92 | 2.05 | 1.96 | 3.29 | 1.91 | 1.98 | 2.03 | 0.61 | 2.79 | 2.58 | 2.79 | 3.48 | 3.57 | 2.22 | 3.26 | 0.71
$SO_{2}$ | 1d | 34.35 | 11.46 | 28.11 | 7.76 | 40.15 | 10.62 | 9.35 | 23.93 | 1.81 | 0.45 | 0.85 | 0.49 | 1.79 | 0.43 | 0.44 | 1.31
7d | 16.83 | 15.43 | 14.62 | 20.72 | 11.59 | 14.39 | 5.84 | 7.30 | 2.36 | 0.34 | 2.23 | 0.61 | 2.21 | 0.33 | 0.30 | 0.61
14d | 10.70 | 7.35 | 9.37 | 24.75 | 9.57 | 8.14 | 11.31 | 6.17 | 2.15 | 0.22 | 2.82 | 1.66 | 3.18 | 0.26 | 0.24 | 0.58
31d | 9.06 | 11.22 | 10.34 | 11.44 | 10.79 | 9.13 | 9.81 | 5.33 | 2.35 | 0.41 | 2.95 | 0.81 | 1.88 | 0.26 | 0.31 | 1.31
$PM_{2.5}$ | 1d | 2.88 | 3.81 | 2.50 | 3.52 | 6.57 | 2.95 | 3.02 | 5.06 | 7.31 | 4.67 | 7.03 | 4.68 | 7.35 | 5.53 | 4.45 | 0.85
7d | 6.94 | 6.86 | 7.05 | 7.38 | 6.69 | 5.92 | 6.57 | 5.07 | 7.09 | 3.01 | 5.22 | 4.22 | 5.56 | 2.65 | 2.43 | 0.58
14d | 4.40 | 6.98 | 4.40 | 7.07 | 4.41 | 6.20 | 6.30 | 2.34 | 8.20 | 4.43 | 8.71 | 7.10 | 7.39 | 3.67 | 5.58 | 0.69
31d | 2.73 | 2.85 | 2.83 | 4.41 | 2.73 | 2.83 | 2.83 | 1.57 | 7.21 | 5.26 | 7.01 | 5.97 | 5.05 | 4.49 | 5.29 | 0.42
$PM_{10}$ | 1d | 1.01 | 6.69 | 1.79 | 6.37 | 4.16 | 3.81 | 5.77 | 3.12 | 7.31 | 4.34 | 7.06 | 4.24 | 7.35 | 5.59 | 4.22 | 0.86
7d | 5.20 | 5.30 | 5.32 | 6.19 | 5.11 | 4.99 | 5.24 | 3.70 | 7.37 | 2.83 | 5.42 | 4.34 | 5.81 | 2.46 | 2.27 | 0.57
14d | 3.53 | 4.59 | 3.55 | 5.90 | 3.54 | 4.36 | 4.89 | 1.48 | 8.48 | 4.29 | 9.08 | 7.08 | 7.65 | 3.29 | 5.44 | 0.67
31d | 2.39 | 2.40 | 2.51 | 3.64 | 2.37 | 2.43 | 2.40 | 0.98 | 7.57 | 5.39 | 7.35 | 6.09 | 5.36 | 4.50 | 5.47 | 0.41

Table I displays the performance of selected forecasting architectures,
incorporating only the best those that perform in at least one scenario. For
the CPCB sensor, the MASE metric values in the table indicate that the
WaveCatBoost model outperforms all the baseline architectures for long-range
forecasting periods of 14 days and 31 days for most air pollutants. This
improved performance is attributed to the hybrid training approach, where
MODWT aids in extracting signals from noise, and CatBoost accurately
extrapolates these transformed signals. In the 7-day forecasting scenario,
WaveCatBoost performs comparably to other boosting methods like LGBM, XGB, and
CatBoost across various pollutants. Furthermore, for short-term forecasts with
a 1-day horizon, WaveCatBoost demonstrates competitive performance with
transformer-based architectures. Additionally, the MASE metric for ID1 sensors
in Table I indicates that our proposed model outperforms all baseline models
across different horizons for all air pollutants except for SO2, where XGB and
CatBoost models offer competitive results.

Furthermore, we study the statistical significance of the improvement in model
performance using multiple comparisons with the best (MCB) test procedure
[15]. This distribution-free test ranks each model based on its MASE metric
across different air pollutants and computes their average rank. The model
with the lowest MASE metric is identified as the “best” method, determined by
achieving the minimum rank. The MCB plots depicted in Fig. 3 illustrate that
WaveCatBoost models exhibit superior performance for both the ID1 (Fig. 3a)
and CPCB (Fig. 3b) sensors. Additionally, the boundary of the critical
distance (blue line) of the best-performing model serves as the reference
value (shaded region) for the test. Given that most of the baseline
forecasters have non-overlapping critical distance (blue line) with the
reference value of the test, we can conclude that their performance is
significantly inferior to that of the proposed WaveCatBoost architecture. The
overall experimental analysis underscores the effectiveness of our proposed
model in outperforming a diverse range of benchmark models, emphasizing its
potential as a reliable tool for real-time air quality forecasting. To
quantify the uncertainty in the WaveCatBoost forecasts, we utilize the
conformal prediction approach as in Eq. 3. Figure 4 and 5 showcase the point
forecast generated by the WaveCatBoost model, ground truth, and the conformal
prediction intervals at $\alpha=0.05$ for the 1-day and 7-day horizon of ID1
sensors, respectively (as it gives a probabilistic band, therefore sometimes
the band consists of negative values which for this case should be considered
as zero). From the plot, it is evident that the proposed framework is
generalizable and robust, offering valuable insights for environmental
monitoring and public health interventions.

##  IV Conslusion

In this letter, we designed a WaveCatBoost model to forecast real-time air
pollutant concentration levels. These air pollutants, including gaseous and
particulate matter, lead to several carcinogenic and non-carcinogenic health
hazards. Our architecture integrates wavelet decomposition with the CatBoost
approach to effectively capture the non-stationary and long-term dependencies
inherent in pollutant time series data. Experimental evaluation using two
real-world datasets from Meghalaya, India, demonstrates that our proposal
outperforms baseline forecasting methods. Furthermore, we assess the
robustness of the model through statistical significance tests and provide
probabilistic bands for our forecasts based on a conformal prediction
approach. These findings are crucial for advancing the field of air quality
forecasts and guiding future research endeavors to enhance predictive models
for environmental sustainability. As a future direction, further analysis of
the WaveCatBoost architecture in forecasting air pollutant levels by
considering spatial dependencies emerges as a promising avenue in air quality
research.

## References

  * [1] K. B. Shaban, A. Kadri, and E. Rezk, “Urban air pollution monitoring system with forecasting models,” _IEEE Sensors Journal_ , vol. 16, no. 8, pp. 2598–2606, 2016.
  * [2] S. Du, T. Li, Y. Yang, and S.-J. Horng, “Deep air quality forecasting using hybrid deep learning framework,” _IEEE Transactions on Knowledge and Data Engineering_ , vol. 33, no. 6, pp. 2412–2424, 2019.
  * [3] J. Borah, S. Kumar, N. Kumar, M. S. M. Nadzir, M. G. Cayetano, H. Ghayvat, S. Majumdar, and N. Kumar, “Aicarebreath: Iot enabled location invariant novel unified model for predicting air pollutants to avoid related respiratory disease,” _IEEE Internet of Things Journal_ , 2023.
  * [4] K. Gu, J. Qiao, and W. Lin, “Recurrent air quality predictor based on meteorology-and pollution-related factors,” _IEEE Transactions on Industrial Informatics_ , vol. 14, no. 9, pp. 3946–3955, 2018.
  * [5] J. Herzen, F. Lässig, S. G. Piazzetta, T. Neuer, L. Tafti, G. Raille, T. Van Pottelbergh, M. Pasieka, A. Skrodzki, N. Huguenin _et al._ , “Darts: User-friendly modern machine learning for time series,” _Journal of Machine Learning Research_ , vol. 23, no. 124, pp. 1–6, 2022.
  * [6] M. Panja, T. Chakraborty, U. Kumar, and N. Liu, “Epicasting: an ensemble wavelet neural network for forecasting epidemics,” _Neural Networks_ , vol. 165, pp. 185–212, 2023.
  * [7] S. Sengupta, T. Chakraborty, and S. K. Singh, “Forecasting cpi inflation under economic policy and geo-political uncertainties,” _arXiv preprint arXiv:2401.00249_ , 2023.
  * [8] L. Sasal, T. Chakraborty, and A. Hadid, “W-transformers: a wavelet-based transformer framework for univariate time series forecasting,” in _2022 21st IEEE international conference on machine learning and applications (ICMLA)_. IEEE, 2022, pp. 671–676.
  * [9] L. Grinsztajn, E. Oyallon, and G. Varoquaux, “Why do tree-based models still outperform deep learning on typical tabular data?” _Advances in neural information processing systems_ , vol. 35, pp. 507–520, 2022.
  * [10] H. Liu and C. Chen, “Spatial air quality index prediction model based on decomposition, adaptive boosting, and three-stage feature selection: A case study in china,” _Journal of cleaner production_ , vol. 265, p. 121777, 2020.
  * [11] L. Prokhorenkova, G. Gusev, A. Vorobev, A. V. Dorogush, and A. Gulin, “Catboost: unbiased boosting with categorical features,” _Advances in neural information processing systems_ , vol. 31, 2018.
  * [12] D. B. Percival and A. T. Walden, _Wavelet methods for time series analysis_. Cambridge university press, 2000, vol. 4.
  * [13] V. Vovk, A. Gammerman, and G. Shafer, _Algorithmic learning in a random world_. Springer, 2005, vol. 29.
  * [14] A. N. Angelopoulos, S. Bates _et al._ , “Conformal prediction: A gentle introduction,” _Foundations and Trends® in Machine Learning_ , vol. 16, no. 4, pp. 494–591, 2023.
  * [15] A. J. Koning, P. H. Franses, M. Hibon, and H. O. Stekler, “The m3 competition: Statistical tests of the results,” _International journal of forecasting_ , vol. 21, no. 3, pp. 397–409, 2005.


# A SPA-based Manifold Learning Framework for Motor Imagery EEG Data
Classification11footnotemark: 1

Xiangyun Lia Peng Chenb<EMAIL_ADDRESS>Zhanpeng Baob West China
Biomedical Big Data Center, West China Hospital, Sichuan University,Chengdu
610041, PR China  School of Mechanical Engineering, Southwest Jiaotong
University, Chengdu 610031, PR China

###### Abstract

The electroencephalography (EEG) signal is a non-stationary, stochastic, and
highly non-linear bioelectric signal for which achieving high classification
accuracy is challenging, especially when the number of subjects is limited. As
frequently used solution, classifiers based on multilayer neural networks has
to be implemented without large training data sets and careful tuning. This
paper proposes a manifold learning framework to classify two types of EEG data
from motor imagery (MI) tasks by discovering lower dimensional geometric
structures. For feature extraction, it is implemented by Common Spatial
Pattern (CSP) from the preprocessed EEG signals. In the neighborhoods of the
features for classification, the local approximation to the support of the
data is obtained, and then the features are assigned to the classes with the
closest support. A spherical approximation (SPA) classifier is created using
spherelets for local approximation, and the extracted features are classified
with this manifold-based method. The SPA classifier achieves high accuracy in
the 2008 BCI competition data, and the analysis shows that this method can
significantly improve the decoding accuracy of MI tasks and exhibit strong
robustness for small sample datasets. It would be simple and efficient to tune
the two-parameters classifier for the online brain-computer
interface(BCI)system.

###### keywords:

Motor imagery electroencephalography , Brain-computer interface , Manifold
learning , Local manifold approximation , Classification algorithm

##  1 Introduction

Electroencephalography (EEG)-based brain-computer interface (BCI) technology
is a promising tool that enables users to use their brain activity to directly
control or interact with an external environment or device [1], especially for
persons with severe motor-related diseases [2]. Despite its low spatial
resolution, EEG is the most widely used technique in the field of BCI due to
its good temporal resolution [3], low equipment cost, and high compatibility
[4]. The BCI technique has been studied via EEG signals arising in different
physiological mechanisms, such as motor imagery (MI) [5], steady-state visual
evoked potentials (SSVEP) [6], and P300 [7]. Compared with BCI based on SSVEP
or P300, MI methods may have higher potential because they are independent of
external stimuli, thus allowing asynchronous control and communication, and
facilitating the application of BCI system [8]. Although there are slight
differences of EEG patterns between, recent advances in this field show that
patients can produce a cortical activation similar to healthy individuals
during the imagined movements [9]. Therefore, the MI-based BCI system is an
effective way to promote recovery during rehabilitation interventions, which
is particularly suitable for tetraplegic individuals to restore neuromuscular
connectivity and motor function [10]. The goal of MI-based system is to
provide a neurophysical communication channel between people and external
devices.For safety purpose of related technologies, noninvasive EEG-based BCIs
are widely used, in such applications as word spellers [11], wheelchair
control [12], and video games [13]. In addition, non-invasive BCI may also be
used in evaluating the brain activity of severely paralyzed patients to
predict the efficiency of invasive brain-computer interfaces [14]. For healthy
persons, BCI can also greatly enhance the experience of multimedia and video
games and thus have have vast commercial value [15].

The major challenges in a typical MI task are the effective extraction of EEG
features and and their proper classification. EEG signals is usually
represented as high dimensional data whose classification is prone to
overfitting and bias, especially for small training datasets [16].
Dimensionality reduction is an effective method to solve these problems for
obtaining a more compact representation of the high-dimensional space.
Manifold learning [17] is one of the most important nonlinear dimensionality
reduction techniques, and the goal of most manifold learning algorithms is to
obtain low-dimensional embedding of high-dimensional data so that the proximal
data points from the high-dimensional space remain close to each other as in
the high-dimensional space, so are the distant data points [18].

With the advantages promoting the classification performance with the limited
data settings due to its simple and direct mathematical representations, many
studies have recently proposed manifold learning algorithms for BCI
applications,and EEG classification on high-dimensional Riemannian manifold
has recently received increasing attention [19]. Barachant et al. [20]
proposed two classification algorithms, one of which introduces the concept of
Riemannian geodesic distance to compare the minimum Riemannian distance
between unlabeled data points and the Riemannian mean of labeled data points.
The other algorithm maps all data points in the Riemannian manifold into its
tangent space (TS),called the optimal hyperplane for classification [21],and
then applies the linear classification method on the tangent space. It is
different from the method of global sampling manifold information in reference
[20]. Xie et al used a simple but effective bilinear sub-meridian learning
(BSML) algorithm [18] in which the intrinsic sub-meridians are learned by
identifying a bilinear mapping that maximally preserves the local geometry and
overall structure of the original manifold for dimensionality reduction of the
SPD matrix space in the motor imagery BCIs. Li et al. [22] proposed a local
manifold approximation classification, which is mainly based on obtaining a
local approximation to the support of the data in each class within the
neighborhood of the feature to be classified, and assigning the feature to the
class with the closest support. The classifier performs a local approximation
using easily fitted spheres to result in a spherical approximation (SPA)
classifier.

Unlike general machine learning tasks, classification tasks for BCI often lack
sufficient labeled training data because of the limited number of subjects. In
practice, It is usually that case there is an insufficient number or even no
training samples at all, known as the small training data set problem. If the
amount of available data is small, the model is difficult to generalize, and
the classification accuracy tends to be lowered [16]. There is significant
variation in EEG signals between different individuals, characterized by the
inconsistent distribution of data in feature space, which prevents prevent the
model transferred from one person to other [23]. Therefore, designing an
algorithm with high classification accuracy constrained by small training data
sets is a challenging task in the current BCI field.

The innovations and contributions of this paper include the following three
aspects:

  1. 1.

The SPA classifier based on manifold learning is modified and applied to the
BCI system for decoding the features of MI-based EEG signals. In the case of
no geodesic information of the unknown manifold, SPA uses sphere for local
approximation, which is different from other manifold learning methods that
approximate geodesic distance from the global point of view. The sphere is a
simple geometric object that is easy to fit, and it also provides the
improvements of the hyperplane, so that the accuracy can be promoted by
approximating the local curvature. It is the first application of SPA in the
classification of EEG signals, especially for MI data.

  2. 2.

Compared with the usual algorithms for a large number of training samples,
This algorithm is suitable for the system such as BCI which lacks enough
training data, and can still present effective feature decoding performance of
EEG signals for a small number of samples with the SPA classifier, as
demonstrated by the experimental results. In clinical rehabilitation training,
such a classifier could work with the limited number of subjects.

  3. 3.

The design of the SPA classifier is simple, and there are only two parameters
needed to be tuned. The main parameters need to be adjusted are the size of
the local neighborhoods K and the dimension of the manifold approximating P.
After these two parameters are determined, the algorithm execution will become
very fast, and become suitable for online BCI systems requiring high real-time
performance.

The rest of this article is divided into the following sections. Section 2
describes the MI task experimental example and basic information on the EEG
data set for the 2008 BCI competition, as well as the preprocessing methods
for EEG signals, and presents the Common Space Pattern Feature Extraction
algorithm for MI-based classification problems. Section 3 describes the SPA
classifier based on manifold learning in detail. In Section 4, the
experimental results of this algorithm and other commonly used algorithms are
compared for evaluation and discussion.In the end, Section 5 provides
conclusions and outlines the future applications of SPA classifiers.

##  2 Motor imagery EEG signals classification preparation

###  2.1 Introduction to the BCI competition dataset

The dataset used in this paper is from the BCI Competition IV Dataset 2a,
which consisted of four MI tasks: motor imagery of the left hand, right hand,
both feet, and tongue for nine subjects. EEG signals from 22 Ag/AgCl
electrodes and 3 monopolar electrooculography (EOG) channels (with the left
mastoid as reference) were recorded at a sampling frequency of 250 Hz. The EEG
signal is band-pass filtered between 0.5 and 100 Hz, and the power line
interference is filtered by a 50 Hz notch filter. The timing scheme of the
experimental paradigm is shown in Figure 1. More detailed information about
the EEG experiment can be found in [24]. Taking into account the non-
stationary of the EEG data, the EEG data consists of two sessions, recorded on
different dates. Each session contains a total of 288 trials, and each type of
task is executed 72 times. It indicates that the dataset includes contaminated
data and such data contains some artifacts. However, these contaminated
samples are usually not removed to examine the robustness of the proposed
algorithm. In the experiment, only the data related to the motor imagery of
the left hand and right hand were extracted from the dataset, while the motor
imagery data of the tongue and feet were excluded.

Figure 1: BCI competition IV dataset 2a experimental paradigm

###  2.2 EEG data preprocessing process

####  2.2.1 Time Interception

The MI task requires subjects to repeatedly imagine limb movements for a
certain amount of time to generate stable and effective brain activity [25].
For the continuously recorded EEG signals, it only focuses on the motor
imagery part and then segments the band-pass filtered continuous EEG signals
through the time window. The segmented time length corresponds to a single
motor imagery trial. Besides, this study used signals between 500 ms and 2500
ms after the stimulus occurred and processed only the intercepted segments.

####  2.2.2 Band-pass filter

During motor imagery, event-related synchronization and desynchronization
(ERS/ERD) phenomena occur mainly in the 8-30 Hz band, covering both mu (8-13
Hz) and beta (14-28 Hz) rhythms [26]. Typically, a bandpass filter is used to
highlight the frequency components of interest while reducing some artifacts.
In this study, a fifth-order Butterworth bandpass filter was designed, which
eliminates most of the muscle activity, such as EMG (above 35 Hz) caused by
swallowing or facial movements. It also removes some ocular artifacts that
mainly affect the low-frequency components of the EEG signal, such as the EOG
(below 8 Hz) generated by blinking and motion.

###  2.3 Common Spatial Pattern (CSP) for Feature Extraction

At present, the spatial domain-based feature extraction method has become the
main feature extraction method for MI EEG signals, which is represented by the
common spatial pattern (CSP). CSP was first proposed by Fukunaga in 1970 and
then introduced into the analysis and processing of EEG signals by Koles [27],
and it has made a very large number of applications on MI EEG signals [28].
CSP is a supervised learning method, which requires the training data to be
labeled. CSP is mainly a feature extraction algorithm for two classification
tasks, capable of extracting the spatial distribution components of each class
from multi-channel brain-computer interface data.

Suppose $X_{1}$ and $X_{2}$ are the multi-channel evoked response time-space
signal matrices under the binary classification MI tasks, respectively. The
dimensions of $X_{1}$ and $X_{2}$ are all N*T, where N is the number of EEG
channels and T is the number of samples collected per channel.

To calculate its covariance matrix, assume that N<T. The normalized covariance
matrices $R_{1}$ and $R_{2}$ of $X_{1}$ and $X_{2}$ are:

$R_{1}=(X_{1}X_{1}^{T})/(trace(X_{1}X_{1}^{T}))R_{2}=(X_{2}X_{2}^{T})/(trace(X_{2}X_{2}^{T}))$
(1)

Then, the mixed space covariance matrix R is calculated by:

$R=\overline{R_{1}}+\overline{R_{2}}$ (2)

Perform the eigenvalue decomposition on the mixed space covariance matrix R:

$R=U\lambda U^{T}$ (3)

Arrange the eigenvalues in descending order, and calculate the whitening value
matrix P:

$P=\sqrt{\lambda^{-1}}U^{T}$ (4)

Then, the normalized covariance matrices $R_{1}$ and $R_{2}$ are transformed
as follows:

$S_{1}=PR_{1}P^{T}\quad S_{2}=PR_{2}P^{T}$ (5)

Next, do the principal component decomposition for $S_{1}$ and $S_{2}$:

$S_{1}=B_{1}\lambda_{1}B_{1}^{T}\quad S_{2}=B_{2}\lambda_{2}B_{2}^{T}$ (6)

The above formula proves that the eigenvector matrices of matrices $S_{1}$ and
$S_{2}$ are equal, i.e., $B_{1}$=$B_{2}$. Sort the eigenvalues in
$\lambda_{1}$ in descending order, and the corresponding eigenvalues in
$\lambda_{2}$ in ascending order. The transformation of the whitened EEG
signal with the eigenvectors corresponding to the largest eigenvalues in
$\lambda_{1}$ and $\lambda_{2}$ is optimal for separating the variances in two
signal matrices. The spatial filter corresponding to the projection matrix W
is:

$W=B^{T}P$ (7)

For the test data $X_{i}$, the feature vector f is extracted as follows:

$\left\\{\begin{array}[]{c}Z_{i}=W*X_{i}\\\
f_{i}=\log\left(\frac{{VAR}\left(Z_{i}\right)}{{sum}\left({VAR}\left(z_{i}\right)\right)}\right)\end{array}\right.$
(8)

The work flow of the system is shown as Figure 2, and the BCI competition IV
Dataset 2a with relatively few dimensions is utilized for classification.

Figure 2: Workflow of the system for EEG data classification

##  3 Principle and Implementation of SPA Classifier

The basic idea of manifold learning is to assume that the dataset under study
is uniformly sampled on a manifold structure in the high-dimensional data
space, then identify the low-dimensional manifold structure hidden in the
high-dimensional observation data space while keeping the neighborhood
relationship between the data sets unchanged, and the nonlinear mapping
relationship from the high-dimensional observation space to the low-
dimensional embedding space is constructed for dimension reduction or
visualization.

Manifold learning has the description [29]: given a dataset
$X=\left\\{x_{i},i=1,\ldots,N\right\\}\subset R^{m}$, assume that the samples
in X are generated from a low-dimensional dataset Y through an unknown
nonlinear transformation f, i.e.:

$x_{i}=f\left(y_{i}\right)+\varepsilon_{i}$ (9)

Among them, $\varepsilon_{i}$ denotes the noise, $y_{i}\in Y\subset R^{d}$,
$d\gg m$, and $f:R^{d}\rightarrow R^{m}$ is the embedding map of $C^{\infty}$.
Thus, the purpose of manifold learning is to get low-dimensional expressions
based on the given observed dataset X:

$Y=\left\\{y_{i},i=1,\ldots,N\right\\}\subset R^{d}$ (10)

Although many global and local methods have been proposed to identify low-
dimensional embeddings, few of them sample the manifold information from the
original data. Most global methods perform low-dimensional embedding by
approximating the geodesic distance without the geodesic information of the
unknown manifold, which can lead to bias [17]. The SPA method introduced in
this paper mainly uses local manifold approximation and easy-to-fit spheres
for local manifold approximation to classify the MI data. The following
describes the details of the SPA classifier.

For the binary classification problem of MI, suppose there are two types of MI
data with different types of motions, labeled as 1 and 2, respectively. The
features of the MI data for the two different types of motions tend to
approach two different manifolds $M_{1}$ and $M_{2}$, and both of which are
embedded in $R^{D}$ with intrinsic dimension $p<D$. The two manifolds of the
EEG signals might be highly nonlinear and complex, even have different
curvatures or even gaps.

For a given test sample x, the distances between the sample and the two
manifolds, denoted by $d_{1}$ and $d_{2}$, are first calculated, and then the
sample is assigned to the group with a shorter distance. Actually, the two
manifolds $M_{1}$ and $M_{2}$ are unknown, but the training data contains n
different sample data and the corresponding labels. With this known
information, it is possible to obtain an exact local approximation of the
manifolds $M_{1}$ and $M_{2}$ in a high-order domain of a feature x to be
processed. In this way, a wide variety of local approximations can be
considered, and the limited training data can be used completely and
efficiently.

Algorithm 1: SPCA(Spherical Principal Components Analysis)algorithm
---
1 Input: $X=\left\\{x_{i}\right\\}_{i=1}^{n}$, Dimension of manifold p
2 Output: Sphere $S_{P}(V,c,r)$
3 $V=\left(v_{1},\ldots,v_{p+1}\right)$, Eigenvalues are in decreasing order,
where $v_{j}$ is the j-th eigenvector of the covariance matrix
4 $\xi_{i}=\bar{X}+V\bar{V}\left(X_{i}-\bar{X}\right)$
5
$c=-\frac{1}{2}\left\\{\sum_{\mathrm{i}=1}^{n}\left(\xi_{i}-\bar{\xi}\right)\left(\xi_{i}-\bar{\xi}\right)^{\top}\right\\}^{-1}\left\\{\sum_{\mathrm{i}=1}^{n}\left(\xi_{i}^{\top}\xi_{i}-\frac{1}{n}\sum_{j=1}^{n}\xi_{j}^{\top}\xi_{j}\right)\left(\xi_{i}-\bar{\xi}\right)\right\\}$
6 $r=\frac{1}{n}\sum_{\mathrm{i}=1}^{n}\left\|\xi_{i}-c\right\|$

The algorithm regards the local spherical approximation as the sphere which is
a simple geometric object easy to fit, and the sphere provides a hyperplane
generalization that can significantly improve accuracy by approximating local
curvature. First, the center, radius and size of each sphere are optimized to
provide the best local approximation. Then, the distances $d_{1}$ and $d_{2}$
between the test sample x and the two manifolds can be easily calculated based
on the local spherical approximation. The sample to be classified by the
distance to different manifolds embed in 3D space is shown as Figure 3.

Figure 3: Sample to be classified by the distance to different manifolds

Let $X_{[k]}^{l}$ be the k-nearest neighbors between x samples with label $l$.
The points are fitted to a sphere using the spherical principal component
analysis (SPCA) algorithm to obtain a local manifold approximation
$\widehat{M_{l}}$ around x in class l. Then, the value of $\mathrm{d}_{l}$ is
approximated by $\widehat{\mathrm{d}_{l}}:=d\left(x,\widehat{M_{l}}\right)$,
and the label y is selected as the value of $l$. The SPCA algorithm produces
an estimated p-dimensional sphere $S_{P}(V,c,r)$ with center c and radius r in
the subspace V when applied to data $X=\left\\{x_{i}\right\\}_{i=1}^{n}$.

Algorithm 2: SPA(SPherical Approximation)algorithm
---
1 Input: training dataset$\left\\{x_{i},y_{i}\right\\}_{i=1}^{n}$, parameters
k and p,test data x
2 Output: Prediction y of test data x
3 Let L be the number of array
4 for l=1:L
5 Find the K-nearest neighbors of x under label l, denoted as
$X_{[k]}^{l}\subset\left\\{x_{i}\mid y_{i}=l\right\\}$
6 Compute the p-dimensional spherical approximation using Algorithm 1
7 Compute the projection of x onto the sphere field S by :
$({\mathord{\buildrel{\lower
3.0pt\hbox{$\scriptscriptstyle\frown$}}\over{x}}^{l}})=Pro{j^{l}}(x)={c_{l}}+\frac{{{r_{l}}}}{{\left\|{{V_{l}}{V_{l}}^{\top}(x-{c_{l}})}\right\|}}{V_{l}}{V_{l}}^{\top}(x-{c_{l}})$
8 Calculate the distance $\widehat{\mathrm{d}}_{l}$ between x and the sphere
$S_{p}^{l}\left(V_{l},c_{l},r_{l}\right):\widehat{\mathrm{d}}_{l}=\left\|x-\widehat{x}^{\imath}\right\|$
9 End
10 Assign x to the group with the smallest distance:
$y=\mathop{argmin}\limits_{l=1,\ldots,L}\widehat{{{\rm{d}}_{l}}}$

The overall process of the SPA classifier is to first obtain the parameters V,
c and r of the spherical classifier by Algorithm 1 (SPCA algorithm), and then
implement the classification by Algorithm 2 (SPA algorithm).

The SPA classifier is simple to implement and efficient to perform, which
could detect non-linear support differences between groups. The parameters to
be tuned for the SPA classifier are the size of the local neighborhood k, and
the dimension p that approximates the data denoising support. When these two
parameters are fixed, the algorithm becomes fast.

For the selection of parameters k and p, k and p are brought in from a certain
range, and the highest accuracy is the accuracy of the test. For the selection
of p, the interval is 1,2,3,4; For the selection of k, the minimum value is
set to 8, increasing in the order of 1, and the maximum value is set to the
number of the small number of the two types of samples. If the number is too
large, the maximum value is set to 46. If p is greater than the number of a
small number of samples, the algorithm will fail.

The SPA classifier theoretically guarantees the impact of the growth of the
number of training samples n on the classification performance. Theorem 1,
which corresponds to the data in the presence of noise, is given below.

Theorem 1: Let $M_{1}$ and $M_{2}$ correspond to two Compact Riemannian
Manifolds. Assume that
$\left\\{z_{i}\right\\}_{i=1}^{n}{}_{\sim}{}^{iid}\rho$, $supp(\rho)=M_{1}\cup
M_{2}$, $z_{i}\in M_{yi}$ and $x_{i}=z_{i}+\epsilon_{i}$, where
$\epsilon_{i}\sim N\left(0,\sigma^{2}I_{D}\right)$. Given a test sample x with
label y, let the SPA classifier obtain the predicated label $\widehat{y_{n}}$.
Let $B_{\delta}(M):=\\{x\mid\mathrm{d}(x-n)<\delta\\}$ and $\delta>0$, then:

$\lim_{n\rightarrow\infty}{P}\left(y\neq\widehat{y_{n}}\right)\leq\rho\left(B_{\delta}\left(M_{1}\right)\cap
B_{\delta}\left(M_{2}\right)\right)+\exp\left\\{-\frac{\delta^{2}}{8\sigma^{2}}+\frac{D}{2}\log\left(\frac{\delta^{2}}{4\sigma^{2}}\right)-\frac{D}{2}(\log(D)-1)\right\\}$
(11)

In theorem 1, the data in class $l$ are distributed around $M_{l}$ as Gaussian
noise. This theorem shows that as the size of training sample n increases, the
probability of the algorithm producing an error class label is gradually
bounded by:

$\rho\left(B_{\delta}\left(M_{1}\right)\cap
B_{\delta}\left(M_{2}\right)\right)+\exp\left\\{-\frac{\delta^{2}}{8\sigma^{2}}+\frac{D}{2}\log\left(\frac{\delta^{2}}{4\sigma^{2}}\right)-\frac{D}{2}(\log(D)-1)\right\\}$
(12)

When the noise level attenuates to 0, i.e., $\sigma\rightarrow 0$. Let
$\delta=\sqrt{\sigma}\rightarrow 0$, so the first term converges to the
$\rho\left(M_{1}\cap M_{2}\right)$, which is the intersection region between
$M_{1}$ and $M_{2}$ are assigned by the probability p. Since
$\frac{\delta^{2}}{\sigma^{2}}\rightarrow\infty$, the second term converges to
0, which is the case that the noise is zero. The probability of correct
classification by the algorithm is ultimately greater than
$1-\rho\left(M_{1}\cap M_{2}\right)$. As the deduction, the limit is 1 when
$\rho\left(M_{1}\cap M_{2}\right)=0$, which means the SPA classifier has
perfect classification performance when it has enough training samples, as
long as the classes are geometrically separable. Theorem 1 takes into account
the noise, which is more practical in most applications.
$\frac{\delta^{2}}{\sigma^{2}}$ can be considered as the signal-noise ratio.
The larger the $\frac{\delta^{2}}{\sigma^{2}}$ , the better the performance of
the SPA classifier.

##  4 Results Evaluation and Discussion

The Accuracy [30] denotes the percentage of the number of correctly classified
samples to the total number of samples, and is calculated as:

$Accuracy=\frac{TP+TN}{TP+TN+FP+FN}\times 100\%$ (13)

In the above formula, True Positive (TP) indicates the number of samples with
positive category and positive classification prediction. True Negative (TN)
indicates the number of samples with negative category and negative
prediction. False Negative (FN) indicates the number of samples with positive
category and negative prediction. And False Positive (FP) indicates the number
of samples with negative category and positive prediction.

Algorithm evaluation: The algorithm proposed in this paper is evaluated by
comparing the CSP+SPA algorithm with the following five algorithms.

  1. 1.

CSP+LDA: Motion imagery classification using CSP and LDA algorithm [31].

  2. 2.

CSP+SVM: Motion imagery classification using CSP and SVM algorithm [32].

  3. 3.

MDRM: The minimum distance to the Riemannian mean is used for classification
on high-dimensional Riemannian manifolds [20].

  4. 4.

TS+LDA: LDA was applied to the high-dimensional tangent space for
classification [20].

  5. 5.

TS+SVM: SVM was applied to the high-dimensional tangent space for
classification [33].

Parameter Settings: In the TS+LDA algorithm, since the tangent space is
$m=n\times(n+1)/2$-dimensional space, there may be a problem that the sample
dimension exceeds the number of trials per class, and the regularized
classification algorithm is usually used to solve this problem. In the TS+SVM
algorithm, the radial basis function kernel (RFB) is selected as the SVM
kernel[32].

The following two experiments are designed for the BCI competition data,
Experiment 1 for the normal number of samples and Experiment 2 for the small
training dataset.

Since cross-validation allows testing the model in the training phase, in this
paper, the performance of the proposed algorithm is first tested on the BCI IV
Dataset 2a using a 10-fold cross-validation procedure. The BCI IV Dataset 2a
is randomly divided into 10 subsets of equal size. In each run, 9 subsets are
elected for training and 1 subset is left for testing.

Table 1: 10-fold cross-validation results (%) Subject | Method
---|---
TS+LDA | TS+SVM | MDRM | CSP+SPA | CSP+SVM | CSP+LDA
S01 | 84.29 | 86.07 | 80.71 | 88.93 | 87.86 | 82.14
S02 | 59.64 | 60.36 | 54.29 | 70.36 | 58.57 | 59.29
S03 | 96.43 | 96.07 | 91.07 | 97.50 | 93.93 | 94.29
S04 | 81.07 | 82.50 | 76.07 | 83.93 | 76.79 | 75.71
S05 | 59.64 | 60.00 | 57.50 | 71.43 | 61.07 | 60.71
S06 | 69.29 | 68.93 | 68.21 | 73.93 | 68.93 | 68.21
S07 | 80.00 | 81.07 | 75.36 | 90.36 | 83.57 | 82.14
S08 | 95.00 | 94.64 | 92.14 | 98.57 | 95.36 | 94.64
S09 | 86.07 | 85.36 | 77.14 | 86.43 | 82.50 | 81.79
Mean | 79.05 | 79.44 | 74.72 | 84.60 | 78.73 | 77.66
Std |  ($\pm$)12.85 |  ($\pm$)12.70 |  ($\pm$)12.34 | ($\pm$)10.04 |  ($\pm$)12.68 |  ($\pm$)12.20

BCI IV Dataset 2a is a four-category dataset, and this paper focuses on
extracting the motion imagery of the left hand and right hand in the dataset.
The classification accuracies of all algorithms are given in Table 1, and the
classification accuracy of CSP+SPA is the highest on each subject. CSP+SPA,
TS+LDA and TS+SVM not only have higher average classification accuracy but
also have a lower standard deviation of classification accuracy. Among them,
CSP+SPA has the lowest standard deviation, which indicates that the proposed
SPA classifier is more robust to the variance of subjects than the other
algorithms. The above results show that the SPA classifier works well for
classifying the data of BCI IV Dataset 2a, which reflects the effectiveness of
the algorithm.

For many reasons, the training set available in BCI applications is usually
quite small. Reducing the number of training trials required for a specific
task is an important goal for BCI, so the experiment is set up to evaluate the
performance of the SPA classifier on small datasets. In the experiment, the
training and testing samples of each subject are aggregated into a total
sample set for that subject, from which different proportions of samples are
selected for the experiment. 1/2, 1/3, 1/6, 1/12 of the BCI IV Dataset 2a
(i.e., 144, 96, 48, and 24 trials randomly selected) were used as training
samples in the evaluation, and 20 replicate experiments were conducted each
time. Table 2 records the average classification accuracy with 1/2 and 1/6 of
the total data as the training dataset. As the training sample size decreases
from 1/2 to 1/6, the performance of all algorithms decreases. However,
compared to other methods, the proposed method demonstrates less performance
degradation and has the highest classification accuracy in most subjects.

Table 2: Accuracy in 1/2 and 1/6 training sets results (%) Sample Number | Method | Mean | subject
---|---|---|---
S01 | S02 | S03 | S04 | S05 | S06 | S07 | S08 | S09
1/2 Training Sample (144 trials) | TS+LDA | 75.65 | 82.08 | 56.84 | 92.74 | 76.11 | 57.08 | 64.72 | 76.49 | 92.15 | 82.60
TS+SVM | 74.57 | 80.20 | 55.83 | 92.66 | 73.84 | 56.36 | 63.64 | 74.88 | 91.67 | 82.02
MDRM | 71.32 | 81.25 | 47.78 | 90.38 | 66.18 | 52.22 | 60.07 | 75.94 | 91.32 | 76.74
CSP+SPA | 77.94 | 85.69 | 60.59 | 94.83 | 75.73 | 61.28 | 65.49 | 81.77 | 95.76 | 80.35
CSP+SVM | 75.25 | 84.39 | 54.86 | 93.96 | 72.36 | 53.91 | 63.47 | 79.67 | 93.91 | 80.73
CSP+LDA | 74.15 | 78.92 | 55.52 | 92.99 | 72.19 | 55.03 | 63.54 | 78.33 | 91.46 | 79.38
1/6 Training Sample (48 trials) | TS+LDA | 70.26 | 77.21 | 53.21 | 89.40 | 66.10 | 53.00 | 60.13 | 66.02 | 88.19 | 79.10
TS+SVM | 67.90 | 72.62 | 55.27 | 85.60 | 60.51 | 54.46 | 58.31 | 60.01 | 85.33 | 78.94
MDRM | 68.64 | 79.15 | 49.40 | 89.90 | 58.54 | 49.67 | 57.56 | 67.77 | 88.50 | 77.31
CSP+SPA | 72.79 | 80.25 | 56.96 | 91.50 | 68.52 | 56.42 | 59.23 | 74.27 | 91.94 | 76.00
CSP+SVM | 71.35 | 79.16 | 54.75 | 90.77 | 64.07 | 51.61 | 58.80 | 71.35 | 91.34 | 80.28
CSP+LDA | 65.62 | 67.71 | 51.67 | 84.04 | 61.52 | 51.48 | 55.48 | 66.23 | 83.79 | 68.67

To further analyze the effect of the number of training samples on the
performance of different subjects, the average classification accuracy of
different subjects in the BCI IV Dataset 2a with the training set percentage
of 1/2, 1/3, 1/6 and 1/12 is given. The line charts of the well-performing
subject S08 and the poor-performing subject S06 are given below.

Figure 4: Classification accuracy of S08 Figure 5: Classification accuracy of
S05

As can be seen in Figure 2, the classification accuracy of the SPA classifier
and the other algorithms on the well-performing subject S08 increases as the
proportion of the training set increases. It is also obvious from the figure
that the SPA classifier has the best performance in different proportions of
the datasets. In Figure 3, there are cases where the classification accuracy
increases despite a small proportion, which is due to the randomness of the
poor performance of the subjects.

##  5 Conclusion and Future Research

MI does not rely on any external visual or auditory stimuli, which facilitates
the implementation of portable BCI systems. It has broad application prospects
in the fields of neural engineering, robotics, and rehabilitation engineering.
Our work proposed a new classifier based on the spherical local manifold
approximation, and optimized it for the best size of the local neighborhoods
and the dimension of manifold facing the application of BCI systems.

To evaluate the performance of the proposed classification method, it runs on
the 2008 BCI competition dataset. The experimental results demonstrate the
effectiveness of the method in decoding and classifying two types of MI tasks,
and the proposed method has higher decoding accuracy and shows strong
robustness for small training datasets compared to other algorithms. The
future work will focus on applying and improving the SPA classifier on the
online rehabilitation devices for clinical experiments.

## References

  * [1] U. Chaudhary, N. Birbaumer, A. Ramos-Murguialday, Brain–computer interfaces for communication and rehabilitation, Nature Reviews Neurology 12 (9) (2016) 513\.
  * [2] I. Lazarou, S. Nikolopoulos, P. C. Petrantonakis, I. Kompatsiaris, M. Tsolaki, Eeg-based brain–computer interfaces for communication and rehabilitation of people with motor impairment: a novel approach of the 21st century, Frontiers in human neuroscience 12 (2018) 14.
  * [3] M. J. Khan, M. J. Hong, K.-S. Hong, Decoding of four movement directions using hybrid nirs-eeg brain-computer interface, Frontiers in human neuroscience 8 (2014) 244.
  * [4] Y. Wang, S. Gao, X. Gao, Common spatial pattern method for channel selelction in motor imagery based brain-computer interface, in: 2005 IEEE engineering in medicine and biology 27th annual conference, IEEE, 2006, pp. 5392–5395.
  * [5] C. Brunner, M. Naeem, R. Leeb, B. Graimann, G. Pfurtscheller, Spatial filtering and selection of optimized components in four class motor imagery eeg data using independent components analysis, Pattern recognition letters 28 (8) (2007) 957–964.
  * [6] Z. Lin, C. Zhang, W. Wu, X. Gao, Frequency recognition based on canonical correlation analysis for ssvep-based bcis, IEEE transactions on biomedical engineering 53 (12) (2006) 2610–2614.
  * [7] D. J. Krusienski, E. W. Sellers, D. J. McFarland, T. M. Vaughan, J. R. Wolpaw, Toward enhanced p300 speller performance, Journal of neuroscience methods 167 (1) (2008) 15–21.
  * [8] N. Naseer, K.-S. Hong, Decoding answers to four-choice questions using functional near infrared spectroscopy, Journal of Near Infrared Spectroscopy 23 (1) (2015) 23–31.
  * [9] S. T. Foldes, D. J. Weber, J. L. Collinger, Altered modulation of sensorimotor rhythms with chronic paralysis, Journal of neurophysiology 118 (4) (2017) 2412–2420.
  * [10] F. Pichiorri, G. Morone, M. Petti, J. Toppi, I. Pisotta, M. Molinari, S. Paolucci, M. Inghilleri, L. Astolfi, F. Cincotti, et al., Brain–computer interface boosts motor imagery practice during stroke recovery, Annals of neurology 77 (5) (2015) 851–865.
  * [11] S. Martens, J. Leiva, A generative model approach for decoding in the visual event-related potential-based brain–computer interface speller, Journal of Neural Engineering 7 (2) (2010) 026003.
  * [12] F. Galán, M. Nuttin, E. Lew, P. W. Ferrez, G. Vanacker, J. Philips, J. d. R. Millán, A brain-actuated wheelchair: asynchronous and non-invasive brain–computer interfaces for continuous control of robots, Clinical neurophysiology 119 (9) (2008) 2159–2169.
  * [13] M. Tangermann, M. Krauledat, K. Grzeska, M. Sagebaum, B. Blankertz, C. Vidaurre, K.-R. Müller, Playing pinball with non-invasive bci., in: NIPS, 2008, pp. 1641–1648.
  * [14] R. Fukuma, T. Yanagisawa, S. Yorifuji, R. Kato, H. Yokoi, M. Hirata, Y. Saitoh, H. Kishima, Y. Kamitani, T. Yoshimine, Closed-loop control of a neuroprosthetic hand by magnetoencephalographic signals, PLoS One 10 (7) (2015) e0131547.
  * [15] A. Nijholt, D. P.-O. Bos, B. Reuderink, Turning shortcomings into challenges: Brain–computer interfaces for games, Entertainment computing 1 (2) (2009) 85–94.
  * [16] Q. Jackson, D. A. Landgrebe, An adaptive classifier design for high-dimensional data analysis with a limited training data set, IEEE Transactions on Geoscience and Remote Sensing 39 (12) (2001) 2664–2679.
  * [17] D. Freedman, Efficient simplicial reconstructions of manifolds from their samples, IEEE transactions on pattern analysis and machine intelligence 24 (10) (2002) 1349–1357.
  * [18] X. Xie, Z. L. Yu, H. Lu, Z. Gu, Y. Li, Motor imagery classification based on bilinear sub-manifold learning of symmetric positive-definite matrices, IEEE Transactions on Neural Systems and Rehabilitation Engineering 25 (6) (2016) 504–516.
  * [19] A. Barachant, S. Bonnet, M. Congedo, C. Jutten, Riemannian geometry applied to bci classification, in: International conference on latent variable analysis and signal separation, Springer, 2010, pp. 629–636.
  * [20] A. Barachant, S. Bonnet, M. Congedo, C. Jutten, Multiclass brain–computer interface classification by riemannian geometry, IEEE Transactions on Biomedical Engineering 59 (4) (2011) 920–928.
  * [21] O. Tuzel, F. Porikli, P. Meer, Pedestrian detection via classification on riemannian manifolds, IEEE transactions on pattern analysis and machine intelligence 30 (10) (2008) 1713–1727.
  * [22] D. Li, D. B. Dunson, Classification via local manifold approximation, arXiv preprint arXiv:1903.00985.
  * [23] A. Craik, A. Kilicarslan, J. L. Contreras-Vidal, Classification and transfer learning of eeg during a kinesthetic motor imagery task using deep convolutional neural networks, in: 2019 41st Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), IEEE, 2019, pp. 3046–3049.
  * [24] M. Tangermann, K.-R. Müller, A. Aertsen, N. Birbaumer, C. Braun, C. Brunner, R. Leeb, C. Mehring, K. J. Miller, G. Mueller-Putz, et al., Review of the bci competition iv, Frontiers in neuroscience 6 (2012) 55.
  * [25] J. Asensio-Cubero, J. Q. Gan, R. Palaniappan, Extracting optimal tempo-spatial features using local discriminant bases and common spatial patterns for brain computer interfacing, Biomedical Signal Processing and Control 8 (6) (2013) 772–778.
  * [26] G. Pfurtscheller, C. Brunner, A. Schlögl, F. L. Da Silva, Mu rhythm (de) synchronization and eeg single-trial classification of different motor imagery tasks, NeuroImage 31 (1) (2006) 153–159.
  * [27] Z. J. Koles, M. S. Lazar, S. Z. Zhou, Spatial patterns underlying population differences in the background eeg, Brain topography 2 (4) (1990) 275–284.
  * [28] H. Ramoser, J. Muller-Gerking, G. Pfurtscheller, Optimal spatial filtering of single trial eeg during imagined hand movement, IEEE transactions on rehabilitation engineering 8 (4) (2000) 441–446.
  * [29] V. De Silva, J. B. Tenenbaum, Global versus local methods in nonlinear dimensionality reduction, Advances in neural information processing systems (2003) 721–728.
  * [30] G. Xu, X. Shen, S. Chen, Y. Zong, C. Zhang, H. Yue, M. Liu, F. Chen, W. Che, A deep transfer convolutional neural network framework for eeg signal classification, IEEE Access 7 (2019) 112767–112776.
  * [31] A. M. Martinez, A. C. Kak, Pca versus lda, IEEE transactions on pattern analysis and machine intelligence 23 (2) (2001) 228–233.
  * [32] H. Sun, Y. Xiang, Y. Sun, H. Zhu, J. Zeng, On-line eeg classification for brain-computer interface based on csp and svm, in: 2010 3rd International Congress on Image and Signal Processing, Vol. 9, IEEE, 2010, pp. 4105–4108.
  * [33] Y. Liang, Y. Ma, Calibrating eeg features in motor imagery classification tasks with a small amount of current data using multisource fusion transfer learning, Biomedical Signal Processing and Control 62 (2020) 102101.