Andreas Meinel
email: andreas.meinel@blbt.uni-freiburg.de

Sebastian Castaño-Candamil

Benjamin Blankertz

Fabien Lotte

Michael Tangermann
email: michael.tangermann@blbt.uni-freiburg.de

S Castaño-C

Characterizing Regularization Techniques for Spatial Filter Optimization in Oscillatory EEG Regression Problems Guidelines Derived from Simulation and Real-World Data

Keywords: EEG bandpower, subspace decomposition, brain-computer interface, single-trial analysis, source power comodulation, brain state decoding algorithm

We report on novel supervised algorithms for single-trial brain state decoding. Their reliability and robustness are essential to efficiently perform neurotechnological applications in closed-loop. When brain activity is assessed by multichannel recordings, spatial filters computed by the source power comodulation (SPoC) algorithm allow identifying oscillatory subspaces. They regress to a known continuous trial-wise variable reflecting, e.g. stimulus characteristics, cognitive processing or behavior. In small dataset scenarios, this supervised method tends to overfit to its training data as the involved recordings via electroencephalogram (EEG), magnetoencephalogram or local field potentials generally provide a low signal-to-noise ratio. To improve upon this, we propose and characterize three types of regularization techniques for SPoC: approaches using Tikhonov regularization (which requires model selection via cross-validation), combinations of Tikhonov regularization and covariance matrix normalization as well as strategies exploiting analytical covariance matrix shrinkage. All proposed techniques were evaluated both in a novel simulation framework and on real-world data. Based on the simulation findings, we saw our expectations fulfilled, that SPoC regularization generally reveals the largest benefit for small training sets and under severe label noise conditions. Relevant for practitioners, we derived operating ranges of regularization hyperparameters for cross-validation based approaches and offer open source code. Evaluating all methods additionally on real-world data, we observed an improved regression performance mainly for datasets from subjects with initially poor performance. With this proof-of-concept paper, we provided a gen-

Introduction

In many modern applications based on biomedical signals, machine learning software is extensively used to infer variables or states of interest [START_REF] Mahmud | Applications of Deep Learning and Reinforcement Learning to Biological Data[END_REF]. Examples are electrocardiographic imaging, where machine learning can be used to deduce cardiac activities or pathologies from multiple sensors [START_REF] Ramanathan | Noninvasive electrocardiographic imaging for cardiac electrophysiology and arrhythmia[END_REF], for controlling an upper-limb prosthesis from an amputee's electromyography recordings [START_REF] Farina | The extraction of neural information from the surface EMG for the control of upper-limb prostheses: Emerging avenues and challenges[END_REF], or for decoding users' mental states from their electroencephalographic (EEG) activity [START_REF] Clerc | SPoC: a novel framework for relating the amplitude of neuronal oscillations to behaviorally relevant parameters[END_REF][START_REF] Clerc | SPoC: a novel framework for relating the amplitude of neuronal oscillations to behaviorally relevant parameters[END_REF]. In all these applications, signals are recorded using multiple sensors, resulting in multivariate data that should be analyzed using robust classification or regression methods. Such machine learning problems are often very challenging.

Prominent neurotechnological systems are brain-computer interfaces (BCI) which typically utilize EEG recordings that enable users to interact with a computer or physical devices [START_REF] Jdr | Combining brain-computer interfaces and assistive technologies: State-of-the-art and challenges[END_REF]. Such practical closed-loop applications require the extraction of relevant and robust features [START_REF] Farquhar | Interactions between pre-processing and classification methods for event-related-potential classification[END_REF] from high-dimensional EEG data which unfortunately suffer from an inherently low signal-to-noise ratio [START_REF] Krusienski | Critical issues in state-of-the-art brain-computer interface signal processing[END_REF][START_REF] Makeig | Evolving signal processing for brain-computer interfaces[END_REF]. In addition, for most BCI applications only small calibration datasets are available to train the decoding algorithms -typically a few dozens or maximally a couple of hundreds of training samples -which further aggravates the situation [START_REF] Lotte | Signal processing approaches to minimize or suppress calibration time in oscillatory activity-based brain-computer interfaces[END_REF]. Thus, it is necessary to design robust decoding methods and training procedures, such that over-fitting to the training data is avoided [START_REF] Makeig | Evolving signal processing for brain-computer interfaces[END_REF].

A widely used approach for effective decoding of EEG signals are spatial filter methods. They learn a linear transformation to project multivariate EEG signals derived from several sensors to a lower dimensional subspace (de Cheveigné and [START_REF] De Cheveigné | Joint decorrelation, a versatile tool for multichannel data analysis[END_REF][START_REF] Blankertz | Optimizing spatial filters for robust EEG single-trial analysis[END_REF] for instance to remove artifacts [START_REF] Vos | Removal of muscle artifacts from EEG recordings of spoken language production[END_REF] or to extract task-related neural activity [START_REF] Makeig | Mining event-related brain dynamics[END_REF]). In the context of BCI, the most prominent algorithm for a supervised scenario is the common spatial pattern algorithm (CSP; [START_REF] Koles | The quantitative extraction and topographic mapping of the abnormal components in the clinical EEG[END_REF][START_REF] Ramoser | Optimal spatial filtering of single trial EEG during imagined hand movement[END_REF]). It is deployed for solving EEG classification tasks that are characterized by amplitude modulations of brain rhythms. Unfortunately, CSP is specifically sensitive towards noisy training data [START_REF] Reuderink | Robustness of the common spatial patterns algorithm in the BCI-pipeline[END_REF], non-stationarities [START_REF] Samek | Divergence-based framework for common spatial patterns algorithms[END_REF] and small datasets [START_REF] Grosse-Wentrup | Beamforming in non invasive brain-computer interfaces[END_REF][START_REF] Park | Filter bank regularized common spatial pattern ensemble for small sample motor imagery classification[END_REF]. To mitigate a subset of these limitations, regularization variants have been proposed for CSP [START_REF] Lotte | Regularizing common spatial patterns to improve BCI designs: Unified theory and new algorithms[END_REF][START_REF] Samek | Divergence-based framework for common spatial patterns algorithms[END_REF]. In general, regularization guides an optimization problem by adding prior information, thus limiting the space of possible solutions. Even though regularization is of specific importance for ill-posed problems such as source reconstruction [START_REF] Tian | EEG/MEG source reconstruction with spatial-temporal two-way regularized regression[END_REF], less underdetermined problems can also profit. For CSP, a broad bandwidth of regularization approaches has been published, such as L1-and L2-norm penalties [START_REF] Wang | Regularized filters for L1-norm-based common spatial patterns[END_REF][START_REF] Lotte | Regularizing common spatial patterns to improve BCI designs: Unified theory and new algorithms[END_REF][START_REF] Arvaneh | Optimizing the channel selection and classification accuracy in EEG-based BCI[END_REF][START_REF] Farquhar | Regularised CSP for sensor selection in BCI[END_REF], regularized transfer learning strategies that accumulate information across multiple sessions and subjects [START_REF] Cheng | Regularized common spatial patterns with subjectto-subject transfer of EEG signals[END_REF][START_REF] Devlaminck | Multisubject learning for common spatial patterns in motor-imagery BCI[END_REF][START_REF] Samek | Transferring subspaces between subjects in brain-computer interfacing[END_REF][START_REF] Kang | Composite common spatial pattern for subject-tosubject transfer[END_REF][START_REF] Lotte | Learning from other subjects helps reducing braincomputer interface calibration time[END_REF] and variants which favor invariant solutions across sessions/runs under EEG non-stationarities [START_REF] Arvaneh | Optimizing spatial filters by minimizing within-class dissimilarities in electroencephalogram-based braincomputer interface[END_REF][START_REF] Samek | Stationary common spatial patterns for brain-computer interfacing[END_REF][START_REF] Samek | Divergence-based framework for common spatial patterns algorithms[END_REF][START_REF] Cho | Increasing session-to-session transfer in a braincomputer interface with on-site background noise acquisition[END_REF].

Taking a closer look into the BCI decoding literature, a variety of methods for oscillatory EEG classification problems can be found, but for the regression case the choice still is extremely limited [START_REF] Wu | Spatial filtering for EEGbased regression problems in brain-computer interface (BCI)[END_REF] even though regression methods allow tackling highly interesting problems. Examples are the estimation of continuous mental workload levels [START_REF] Frey | Framework for electroencephalography-based evaluation of user experience[END_REF][START_REF] Schultze-Kraft | Unsupervised classification of operator workload from brain signals[END_REF], decoding the depth of cognitive processing [START_REF] Nicolae | Assessing the depth of cognitive processing as the basis for potential user-state adaptation[END_REF], predicting singletrial motor performance [START_REF] Meinel | Pre-trial EEGbased single-trial motor performance prediction to enhance neuroergonomics for a hand force task[END_REF] or continuous decoding of movement trajectories ( Úbeda et al, 2017). A spatial filtering solution, which solves an EEG regression problem, was provided by Dähne et al (2014) with the source power comodulation algorithm (SPoC). It optimizes spatial filters that describe oscillatory subspace components, whose bandpower co-modulate with a given continuous univariate target variable. Comparing the formulations of the objective functions of SPoC and CSP, both can be translated into a Rayleigh quotient. Thus, the known limitations of CSP regarding noise, non-stationarities and limited data were found to also apply to SPoC [START_REF] Castaño-Candamil | Probing meaningfulness of oscillatory EEG components with bootstrapping, label noise and reduced training sets[END_REF]. However, regularization approaches for SPoC have not yet been explored.

The main goal of this paper is to close this gap. Thus, we present generally applicable regularization variants for SPoC to improve the algorithm's robustness. Therefore, we first evaluate regularized SPoC variants on a very recent simulation approach based on post-hoc labeling of arbitrary EEG recordings. This allows probing the stability of the regularized variants under reduced training datasets, varying label noise conditions and different strengths of oscillatory sources. In a second step, we evaluate the methods on a real-world dataset and compare the findings of both scenarios. As regularization introduces additional hyperparameters, we will compare model selection via cross-validation with an analytical solution. Finally, we provide the practitioner with information on how to determine suitable parameters for SPoC regularization and provide open source software for regularized SPoC. This paper builds upon preliminary results reported in a conference paper by [START_REF] Meinel | Tikhonov regularization enhances EEGbased spatial filtering for single-trial regression[END_REF].

2 Supervised Spatial Filter Optimization for Single-Trial EEG Regression Let x(t) ∈ R N c describe multivariate EEG signals acquired from N c channels at time sample t. Moreover, the signals are assumed to be bandpass filtered. A spatial filter w ∈ R N c can be interpreted as a projection of x(t) from the N c -dimensional sensor space to a one-dimensional source component ŝ(t) = w x(t) according to the linear model of the EEG Parra et al (2005).

Prerequisites for the SPoC algorithm

Based upon the course of events in the conducted experimental paradigm, an EEG recording x(t) can be translated into N e single epochs X(e) ∈ R N c ×N s with N s sample points per epoch. Hereafter, • refers to an average across epochs. For each epoch e, we assume to have access to a continuous epoch-wise target variable z(e) which is required to be standardized to zero mean and unit variance. This variable z provides label information about the experimental paradigm and could represent e.g. stimulus intensity, behavioral responses or cognitive measures. Subsequently, our goal is to search for a source component ŝ within the full EEG signal, whose epoch-wise power Φ ŝ(e) is linearly related to the target variable z(e).

A spatial filter w defines an estimated source ŝ when applied to the EEG signal. Due to the preceding bandpass filtering, its power can be assessed by its variance within the epoch e such that Φ ŝ(e) = Var[ŝ(t)](e) = Var[w x(t)](e).

The central principle of SPoC is to search an optimal spatial filter w * such that the epoch-wise power of the resulting estimated source component ŝ maximally co-modulates with the known target variable z. Once this optimal filter has been found, it can then be employed to estimate the target variable z(e) from the bandpassed signal X(t, e)1 using the variance:

z est (e) ≈ Φ ŝ(e) = Var[ŝ(t)](e) = w Σ(e)w (1)
where Σ(e) = (N s -1) -1 X(e) X(e) denotes the epoch-wise spatial covariance matrix.

Striving to find a formulation of the overall optimization function requires two ingredients. The first one is expressed by the (Euclidean) mean of the power Φ ŝ(e) across epochs by:

Φ ŝ(e) Eq. (1) = w Σ(e) w = w Σ avg w (2)
with Σ avg := Σ(e) defining the averaged covariance matrix across all N e epochs.

The second ingredient is expressed as the covariance between the epoch-wise power of the source component and the target variable z. One can show the following relation [START_REF] Clerc | SPoC: a novel framework for relating the amplitude of neuronal oscillations to behaviorally relevant parameters[END_REF]:

J 1 (w) := Cov[Φ ŝ(e), z(e)] = w Σ z w (3)
where Σ z := Σ(e) z(e) defines the label-weighted covariance matrix averaged over epochs. While the original SPoC formulation comprised two different optimization strategies, we will restrict further analysis to SPoC λ which optimizes covariance instead of correlation, but allows deriving a closed-form solution of the spatial filters.

Optimizing covariance -SPoC λ algorithm

As the covariance is directly affected by the scaling of its arguments, it requires a constraint upon possible solutions. In SPoC λ this is tackled by a filter norm constraint J 2 (w) := Var[ŝ](e) = w Σ avg w = 1 which translates into the following Rayleigh quotient:

J λ (w) = J 1 J 2 = w Σ z w w Σ avg w (4)
The optimization task w * = argmax w J λ (w) can be transfered into a generalized eigenvalue problem (de Cheveigné and [START_REF] De Cheveigné | Joint decorrelation, a versatile tool for multichannel data analysis[END_REF] and thus delivers a closed-form solution. Overall, the approach returns a full set {w * j } j=1,..,N c of N c spatial filters with j indexing the rank. It is determined in descending order of the eigenvalues and thereby according to the covariance. Throughout the remaining paper, the SPoC λ algorithm is used. It will be referred to by the term SPoC.

Regularization for regression based spatial filtering

In most BCI scenarios small training datasets of a high dimensionality are encountered [START_REF] Makeig | Evolving signal processing for brain-computer interfaces[END_REF]. In this setting, SPoC shows an impairing sensitivity and thus might be prone to overfit the training data [START_REF] Castaño-Candamil | Probing meaningfulness of oscillatory EEG components with bootstrapping, label noise and reduced training sets[END_REF]. A common machine learning strategy in such situations is to add prior information and thus regularize the objective function of an algorithm.

Similar to the regularization strategies proposed by [START_REF] Lotte | Regularizing common spatial patterns to improve BCI designs: Unified theory and new algorithms[END_REF] for CSP, there are two possible branches of regularization strategies for the SPoC algorithm: The first is to directly add prior information on the level of the objective function in Eq. 4. This leads to a restriction of the solution space of possible filters. The second one directly addresses the involved empirical covariance matrices which suffer from small training sets and noisy data. Poorly estimated covariance matrices will not characterize the intended neural activity well. Therefore, regularization on the level of covariance matrices intends to improve their estimation and thus enhance the spatial filtering optimization. In the following, we will propose two regularization approaches, one from each branch of strategies.

Additional penalty on the objective function

Introducing a regularization to the objective function of SPoC can be achieved by adding a penalty term P (w) to the denominator of the Rayleigh quotient stated in Eq. 4:

J λ (w) = w Σ z w (1 -α)w Σ avg w + αP (w) (5)
where α ≥ 0 is the regularization parameter that modulates the strength of the penalty. In general, the term P (w) penalizes solutions of w that do not fulfill a specified prior. Thereby it increases the algorithm's robustness to outliers and small training sets.

In this paper, we select a simple quadratic penalty of the form:

P (w) = w Iw = w 2 (6)
using the identity matrix I ∈ R N c ×N c . This penalty is known as Tikhonov regularization [START_REF] Tikhonov | Regularization of incorrectly posed problems[END_REF] and has similarly been introduced for CSP [START_REF] Lotte | Regularizing common spatial patterns to improve BCI designs: Unified theory and new algorithms[END_REF]. As the penalty P (w) scales with the spatial filter norm, solutions w with small weights are preferred. Regarding utmost regularization strength in Eq. 5 expressed by α = 1, the Rayleigh quotient simplifies to the one of the principal component analysis (PCA, (De Bie et al, 2005)) meaning that a PCA on the z-weighted covariance matrix is computed. For the introduced Tikhonov regularization of SPoC, model selection wrt. α can be done via cross-validation (CV).

Trace normalization

SPoC λ optimizes covariance which is not scale-invariant. This drawback might be mitigated by the norm constraint, but to directly compensate for the relative scaling of the covariance matrices in Eq. 5, a normalization of all covariance elements by the trace tr[•] might also be a suitable strategy as already proposed for CSP [START_REF] Ramoser | Optimal spatial filtering of single trial EEG during imagined hand movement[END_REF][START_REF] Lu | Regularized common spatial pattern with aggregation for EEG classification in small-sample setting[END_REF]:

Σ(e) = Σ(e) tr[Σ(e)]
;

Σ avg = Σ avg tr[Σ avg ] (7)
Here, we investigate the effect of applying trace norm to Σ(e) and Σ avg entering Eq. 5, but not upon Σ z as its label-weighting shall be maintained.

Equivalence to covariance shrinkage

Inserting the given Tikhonov penalty P (w) of Eq. 6 into the objective function in Eq. 5, enables to factorize the denominator to a shrinkage of the averaged covariance matrix Σ avg towards the identity matrix I ∈ R N c ×N c :

Σ avg = (1 -α)Σ avg + αI (8)
By that we have shown that substituting Σ avg by the shrinked version Σ avg in the objective function of SPoC (see Eq. 4) is equivalent to the Tikhonov formulation stated in Eq. 5 and 6.

Regularization of covariance matrices

In parallel to the proposed Tikhonov regularization which builds upon a CV procedure for model selection, there are faster ways of determining a suitable regularization strength. We will focus on two strategies for covariance shrinkage which allow to use an analytic solution to determine the regularization parameter.

Automatic shrinkage of sample covariance matrices

When estimating a sample covariance matrix S ∈ R N c ×N c based on N train training data samples, there was a systematic bias reported in the setting of N c > N train : large eigenvalues get overestimated while small eigenvalues tend to be underestimated [START_REF] Bai | Spectral Analysis of Large Dimensional Random Matrices[END_REF]. The situation can be improved by shrinking the covariance matrix S towards the identity matrix I ( [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF][START_REF] Schäfer | A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics[END_REF]:

S = (1 -α * )S + α * νI (9)
Under the assumption of i.i.d. data and thus in the absence of outliers, Ledoit & Wolf derived a closed-form solution for the optimal shrinkage parameter α * and the optimal scaling parameter ν by minimizing the expected mean squared error.

For the exact closed-form solution of α * and ν, we refer the reader to [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF]; [START_REF] Schäfer | A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics[END_REF]; [START_REF] Bartz | Covariance shrinkage for autocorrelated data[END_REF]. This closedform solution holds the advantage of directly computing an estimate of α * without cross-validation. Note that the additional scaling factor ν takes a similar role as the trace normalization introduced for Tikhonov regularization with the difference that it only takes diagonal terms into account.

Automatic shrinkage of averaged covariance matrix

As shown in section 3.1, the Tikhonov penalty introduced for the SPoC objective function can be rewritten as a shrinkage of the covariance matrix Σ avg , which was gained by averaging across the epoch-wise covariances Σ(e). Thus, one can directly apply the closed-form solution for α * and ν, but it first requires to estimate the averaged covariance matrix as Σ avg = (N s • N e -1) -1 x cat x cat using a concatenated data matrix of all N e epochs, namely

x cat = [x(1), ..., x(N e )] ∈ R N c ×(N s •N e ) .
To compensate for signal non-stationarities, each data epoch X(e) should be corrected to channel-wise zero mean prior to concatenation.

Automatic shrinkage of epoch-wise covariance matrix

SPoC includes the label-weighted covariance matrix Σ z which holds all the available label information. A direct covariance shrinkage for Σ z was tested in pilot experiments (data not shown), but this turned out not to be beneficial -probably because adding a regularization term would diminish the contained label information. As both Σ z and Σ avg require a computation of the epoch-wise covariance Σ(e), we propose to choose this matrix as regularization target using the previously mentioned closed-form solution for α * and ν in order to derive a shrinked estimate Σ(e).

Overview on evaluated SPoC regularization variants

In section 3.1 and 3.2 different regularization strategies were introduced. In Table 3.3 an overview over all proposed approaches is given. The first three rows summarize Tikhonov regularization variants which all require an estimation of α by means of cross-validation. Among them, 'Tik-SPoC' comprises Tikhonov regularization only according to Eq. 5, while 'NTik-SPoC' considers an additional trace norm both for Σ(e) and Σ avg . The largest extent of regularization is realized by 'ASNTik-SPoC' which uses the same strategy as 'NTik-SPoC' with additional automatic shrinkage on Σ(e) for the computation of Σ z . As this term applies to the numerator (N) of the objective function, this is marked accordingly in Table 3.3. The last two rows in Table 3.3 summarize automatic shrinkage approaches using the closed-form solution by Ledoit & Wolf (LW). Applying automatic shrinkage to the averaged covariance matrix will be referred to as automatic Tikhonov regularization 'aTik-SPoC'. In contrast, using automatic shrinkage directly upon Σ(e) in the numerator (E) and denominator (D) of the objective function will be referred to as 'AS-SPoC'.

4 Experiments and validation procedure

Simulation data

In this work, we aim to characterize and benchmark the introduced regularization techniques to the SPoC λ algorithm. However, in the majority of real-world EEG experiments there is no ground truth s available which severely challenges the validation procedure. To compensate for this, a novel data-driven simulation approach for labeling datasets was utilized [START_REF] Castaño-Candamil | Post-hoc labeling of arbitrary EEG recordings for data-efficient evaluation of neural decoding methods[END_REF]. It generates ground-truth label information based on known sources from arbitrary pre-recorded EEG measurements. This post-hoc data labeling allowed obtaining noiseless labels from a relatively large amount of EEG data (here up to 1000 epochs) while conserving the real statistics of the neural activity including nonstationarities of the signal. Furthermore, the approach provided full control over label noise and allows studying its influence upon the decoding performance. In the following, a detailed description for the dataset generation is given.

Preprocessing

In total, 40 datasets of a single motor imagery session per subject formed the basis for the simulation. The experimental design of the motor imagery paradigm is described in detail in [START_REF] Blankertz | Neurophysiological predictor of SMR-based BCI performance[END_REF]. From the recorded EEG, we utilized the signals of 63 passive EEG channels placed according to the extended 10-20 system. The preprocessing of each raw EEG dataset consisted of a high-pass filtering at 0.2 Hz, low-pass filtering at 48 Hz and sub-sampling to 120 Hz. For each dataset, the continuous EEG recordings of active task periods (from the task cue to the end of the imagery interval) were segmented into non-overlapping epochs of 1000 ms duration. Artifact epochs were identified by a min-max threshold and by a variance criterion. The latter was additionally applied to detect and remove outlier channels. Details about the artifact preprocessing are described in [START_REF] Meinel | Pre-trial EEGbased single-trial motor performance prediction to enhance neuroergonomics for a hand force task[END_REF].

As SPoC's stability depending on the number of training epochs N e shall be studied in this paper, we discarded datasets with N e < 1000. Similarly, datasets where more than 10% of the original EEG channels had to be rejected, have been removed from further analysis. Applying these criteria, the data of 12 out of 40 subjects remained.

Post-hoc labeling of pre-recorded EEG data

As illustrated in Fig. 1 (A), the following steps were applied to generate continuous labels z true from pre-recorded EEG datasets in a data-driven way:

1. Bandpass filtering of the data to a frequency band of interest. For our analysis, we choose the alpha-band frequency range of [8,12] Hz. 2. Based on the bandpass filtered data, an ICA decomposition (fastICA, [START_REF] Hyvarinen | Fast and robust fixed-point algorithms for independent component analysis[END_REF]) into N in = 20 independent components (ICs) was computed. 3. To identify and remove artifactual components in an automatic way, the datadriven classification approach MARA [START_REF] Winkler | Robust artifactual independent component classification for bci practitioners[END_REF] for the identification of artifactual components was applied. A posterior probability threshold (p art = 10 -8 ) describing the probability of an artifact feature was applied for discarding components of non-neural origin resulting in N sel ≤ N in selected ICs. 4. The log-bandpower for each selected component j, with j = {1, .., N sel }, was computed by the Hilbert transform and averaged in each 1 s time interval which defined the epoch-wise known target variable z j (e) as sketched in Fig. 1 (B).

Overall, the preprocessed data of 12 subjects resulted in 145 oscillatory components (≈ 12 per subject) which survived MARA. For each selected IC, the logbandpower activation was sampled across N e = 1000 epochs and thus delivered continuous epoch-wise labels z true (e) to the respective epoched EEG signals X(e).

We expect the SPoC decoding accuracy to be sensitive to the strength of envelope changes of an oscillatory component. The simulation design enables to empirically study this influence by extracting the absolute width in bandpower fluctuations of a single selected IC across a full session. Therefore, we define the fluctuation width of the jth IC as σ z := Var[z j (e)] calculated across the N e = 1000 epochs as illustrated in Fig. 1 (C).

Probing the algorithms under reduced datasets and label noise

In an offline analysis using the generated 145 labeled datasets, all five introduced SPoC regularization variants and the standard SPoC approach were evaluated in a 10-fold chronological CV. For each epoch e, an estimate of the target variable z est (e) was derived according to Eq. 1 by selecting the highest ranked spatial filter obtained from the training data.

To analyze the benefit of regularization under different dataset sizes, we evaluated the algorithms' stability by systematically reducing each dataset with originally 1000 epochs to smaller data chunks. Therefore, epochs from the session end were removed. For each of the 145 labeled datasets, 22 discrete, logarithmically scaled dataset sizes N e ∈ [20, 1000] (respectively training set sizes N train ) were tested. Similarly, we probed the stability of our approaches under varying label noise conditions. Therefore, each sample of the target variable distribution z true was modified by adding normally distributed label noise, resulting in a noisy label set z noisy which was used for the CV procedure. According to the label noise model proposed by [START_REF] Castaño-Candamil | Post-hoc labeling of arbitrary EEG recordings for data-efficient evaluation of neural decoding methods[END_REF], the correlation between the undistorted and the noisy labels ρ n = Corr(z true , z noisy ) can be controlled via the label noise parameter ξ n := 1 -ρ n . A value ξ n = 1 refers to maximal label noise, while ξ n = 0 indicates that the labels are completely noise free. Five fixed levels for ξ n were evaluated.

For the CV-based regularized SPoC variants (see Table 3.3), the regularization strength α was varied in a range α ∈ {0; [10 -8 , 10 0 ]}. Overall, 20 discrete, logarithmically scaled α levels were analyzed. To summarize, we tested all algorithms on different hyperparameter sets Ω = {(N train , ξ n , α)}.

Real-world scenario

Dataset for evaluation

In order to examine the regularization methods in a real-world decoding scenario, we utilized data of 18 subjects who participated in a repetitive visuo-motor hand force task with 400 trials per session. The task enables to extract a single-trial motor performance metric such as reaction time (RT) or a cursor path length (CPL). During the full session, EEG from 63 passive Ag/AgCl electrodes (Easy-Cap) placed according to the extended 10-20 system was recorded by multichannel EEG amplifiers (BrainAmp DC, Brain Products) with a sampling rate of 1 kHz. In each trial, a "get-ready" interval preceded a visually presented "go-cue", the latter indicating the start of a motor execution phase. In an offline analysis by [START_REF] Meinel | Pre-trial EEGbased single-trial motor performance prediction to enhance neuroergonomics for a hand force task[END_REF], we found that oscillatory bandpower features recorded during the get-ready interval can partially explain upcoming single-trial motor performance. For further details, please see [START_REF] Meinel | Pre-trial EEGbased single-trial motor performance prediction to enhance neuroergonomics for a hand force task[END_REF].

Building upon these findings, in this paper the EEG signals were segmented into epochs along the time interval [-500, +500] ms relative to the go-cue in each trial to decode RT of the upcoming motor task. After data preprocessing and outlier rejection following the workflow described in [START_REF] Meinel | Pre-trial EEGbased single-trial motor performance prediction to enhance neuroergonomics for a hand force task[END_REF], we now restricted any further analysis to oscillatory features within the alpha-band frequency range of [8,12] Hz. The bandpass filter was realized applying a zero-phase butterworth filter of 6th order. The number of epochs N e surviving the preprocessing varied across subjects from 142 to 352.

Evaluation scheme

All algorithms were evaluated within a (nested) 10-fold chronological CV. The three CV-based regularization variants demanded an additional inner CV to estimate the individually optimal regularization parameter α * . It was chosen among 15 discrete, logarithmically scaled values in the range α ∈ [10 -8 , 1]. The α-value maximizing the z-AUC evaluation score (details see section 4.3) was selected and applied to the outer CV in order to train the respective spatial filtering algorithm and linear regression model. The methods 'aTik-SPoC' and 'AS-SPoC' allow for an analytical estimate of α * and hence did not require an inner CV. In contrast to the simulation scenario, the total number of ground truth neural source(s) which might (partially) explain the target variable z true is not known a priori. By applying a regression model, we assume that several sources might contribute to explain the labels z true .

For each α in the inner or outer CV, the following scheme was applied: a spatial filter set {w (i) } i=1,..,N c was computed on training data x tr . The first N f eat = 4 highest ranked components were selected as input to train a linear regression model with coefficients {β j } j=0,..,N f eat . The model was trained upon the log-bandpower features Φ j,tr = log(Var[w (j) tr x tr ]). On each feature Φ j,tr , the mean µ j,tr and the variance σ j,tr was estimated in order to standardize the data to zero mean and unit variance before entering the regression model. Given unseen test data x te , the logbandpower features Φ j,te (e) = log(Var[w (j) tr x te ])(e) for each selected spatial filter w (j) tr were first standardized by µ j,tr and σ j,tr . Subsequently, the corresponding coefficients β j of the trained linear regression model enabled to estimate the target variable z est (e) via:

z est (e) = β 0 + N f eat j=1 β j Φ j,te (e) (10)

Evaluation scores

To compare the estimated labels z est with the known or measured labels z true in the simulation and real-world scenarios across the proposed regularization variants, different evaluation scores can be considered [START_REF] Meinel | Pre-trial EEGbased single-trial motor performance prediction to enhance neuroergonomics for a hand force task[END_REF]. In general, the Pearson correlation coefficient could be utilized but has the drawback, that it is very sensitive to the number of samples [START_REF] Kenney | Mathematics of statistics[END_REF]. Therefore, we instead decided to utilize the following three scores in this paper:

-Angle θ between spatial filters: The design of the simulation scenario gives access to each ground truth spatial filter w true . As all proposed SPoC variants directly optimize for a spatial filter estimate w with arbitrary sign and amplitude (this characteristic is inherited from the formulation as an eigenvalue problem), the angle θ between the spatial filters can directly serve as an evaluation metric:

θ r = arccos w w true w w true θ = θ r , θ r ≤ π/2 π -θ r , θ r > π/2 (11)
with 0 < θ < π. A perfect decoding will be expressed by an angle θ = 0. Please note that the angle θ can only be estimated within the simulation scenario. -Separability z-AUC of labels: Another possibility is to transfer the continuous labels z true into a two-class scenario according to the median of z true . This enables the utilization of the receiver operating characteristics (ROC) curve which is calculated upon the estimated target variable z est given the true two-class labels [START_REF] Fawcett | An introduction to ROC analysis[END_REF]. As ROC performance can be reduced to a scalar value by calculating the area under the ROC curve (AUC), we will name this metric z-AUC as it characterizes the separability of the estimated target variable z est . The z-AUC score can be directly evaluated in both scenarios [START_REF] Meinel | Tikhonov regularization enhances EEGbased spatial filtering for single-trial regression[END_REF]. A perfect decoding corresponds to z-AUC = 1 while chance level correspondents to a value of 0.5. -Relative z-AUC performance: The score z-AUC ref corresponds to the baseline performance of SPoC without any regularization. In this paper, we will compare it to performances obtained by the proposed regularized variants (see Table 3.3). Given a hyperparameter set Ω, the target variable obtained under these hyperparameters z est (w(Ω)) can be estimated using Eq. 1 and the corresponding z-AUC can be computed. For fixed Ω, the performance of a regularized SPoC variant z-AUC reg (Ω) can be assessed as the relative change of z-AUC to the baseline SPoC performance:

rel. z-AUC(Ω) := z-AUC reg (Ω) -z-AUC ref (Ω) z-AUC ref (Ω) (12)
If rel. z-AUC > 0, this directly corresponds to a relative performance increase compared to SPoC and vice versa.

Results

First, we studied the characteristics of the regularization algorithms on 145 analysis problems within the simulation framework. It allows assessing the influence of (hyper)parameters such as regularization strength, dataset size or label noise under controlled conditions. Second, the approaches were tested on real-world data to verify the transferability of the findings and to provide rules of thumb for the practitioner.

Simulation data

Labeling according to bandpower fluctuation width

The SPoC algorithm optimizes for oscillatory components that co-modulate in their bandpower with a given target variable. In Fig. 2, the relation between the fluctuation width σ z and the baseline SPoC performance z-AUC ref on the full dataset N e = 1000 is shown for each of the 145 ICs (correlation R = 0.31 with p = 2.20 • 10 -4 ). The results indicate that the decoding quality of SPoC depends on the fluctuation width σ z of the underlying neural component, with stronger fluctuation width being related to higher decoding quality. For further analysis of the simulation data, all 145 ICs were labeled according to their bandpower fluctuation width σ z into three classes determined by the lower and upper quartile according to the distribution of σ z across all components (see color coding in Fig. 2). In the following, we will show the decoding performances z-AUC GA and θ GA as grand average for each corresponding fluctuation width class.

Sensitivity to regularization parameter

Regarding the CV-based regularized SPoC versions, their sensitivity to the regularization parameter α is reported in Fig. 3 exemplarily for the 'high σ z ' class. It reflects the grand average (GA) of all components contained in this class and provides different evaluation scores. The first row reports the z-AUC GA while the second row summarizes the angle θ GA between filters. A regularization benefit is expressed via an increasing z-AUC GA or a decreasing θ GA relative to the performance level at α = 10 -8 . A few observations can be summarized from Fig. 3: First, the two evaluation scores z-AUC and θ are highly (anti-)correlated across the shown dataset scenarios and SPoC regularization variants. As in real-world data the ground truth will not be known a priori, further analysis will need to be restricted to the metric z-AUC. Second, an increase of the training set size N train (left to right column) leads to a lower sensitivity wrt. α. Third, a comparison of α sensitivity ranges across the three regularization variants yields that 'NTik-SPoC' and 'ASNTik-SPoC' are sensitive in the interval 10 -6 ≤ α ≤ 1 while 'Tik-SPoC' is only sensitive within 10 -3 ≤ α ≤ 1. Fourth, 'NTik-SPoC' and 'ASNTik-SPoC' behave highly similar, while 'Tik-SPoC' shows a qualitative different behavior.

Based on these observations, further analysis will focus on differences between 'NTik-SPoC' and 'Tik-SPoC'. Fifth, extreme regularization with α = 1 leads to a drop of decoding performance regardless of the approach, while in the absence of regularization (α = 0) a slight improvement due to trace normalization can be reported for 'NTik-SPoC'.

Influence of reduced datasets and fluctuation width

The simulation scenario grants access to test the stability of different regularized SPoC variants under reduced datasets. For the CV-based methods, a sensitivity analysis for the regularization strength α under 22 training set sizes N train is shown in Fig. 4 for 'Tik-SPoC' (first row) and 'NTik-SPoC' (second row). The two columns in Fig. 4 reveal the influence of the components' fluctuation width σ z (left: low, right: high). We observed, that regularization has the strongest effects for components with large σ z and for small training sets. With increasing training set size N train , the sensitivity range for α shifts towards smaller α values. Comparing the depicted methods, 'NTik-SPoC' shows a higher sensitivity to regularization strength α compared to 'Tik-SPoC'. Interestingly, for all subplots (A)-(D) the curves along different N train values converge at α = 1, as for this value the SPoC methods collapse to a PCA on the z-weighted covariance. Even for this extreme choice of α, data characterized by higher σ z reaches a better decoding performance than data with lower σ z .

To quantify the decoding performances across methods, the maximum GA performance z-AUC max := z-AUC GA (α * ) is reported in Fig. 5 (A) and(B) in the absence of label noise. Therefore, the optimal regularization strength α * = arg max α z-AUC GA (α) is selected for fixed N train and σ z class. For variants using the LW estimate, this selection is not necessary as there is an analytic solution for α * such that z-AUC max = z-AUC GA . Accordingly, the relative performance change rel. z-AUC(α * ) is reported on the GA level in Fig. 5 and (F) report the statistical significance of the findings. Therefore, a one-sided Wilcoxon signed rank test was applied to test if the median of performance differences (z-AUC max,ref (Ω) -z-AUC max,reg (Ω)) is smaller or equal to zero for fixed N train and σ z . If a p-value p < 0.05 was found (not corrected for multiple testing), the configuration Ω reveals a significant difference among the two methods, indicated by a colored data point in Fig. 5 (E) and (F). The following observations can be reported: First, the absolute decoding performance strongly depends on N train regardless of the regularization method and σ z class. Second, there is a relative performance increase of all introduced regularization methods up to training sets of size N train ≈ 60 on the grand average level. For larger datasets the regularization does not reveal an additional benefit on the grand average. Third, our results indicate that regularization is beneficial for various methods in the 'high σ z ' class, while this is not the case for 'low σ z '. Here, a noticeable case is reported by the performance of 'AS-SPoC' which drastically looses performance for N train 50.

Stability under label noise and reduced data

As in most real-world scenarios label noise challenges the decoding performance of subspace methods like SPoC [START_REF] Castaño-Candamil | Probing meaningfulness of oscillatory EEG components with bootstrapping, label noise and reduced training sets[END_REF]. Thus, we studied its influence for reduced datasets within the simulation data. Fig. 6 exemplary shows the degrading decoding performance under label noise conditions for 'aTik-SPoC' and 'AS-SPoC' for 'high σ z '. Both methods have in common, that performance estimates are very noisy under small dataset size and increasing label noise. Regarding the maximally achievable decoding performance for both meth- 7 Simulation results: landscape of the grand average relative performance changes in z-AUC dependent on the training set size N train and regularization strength α for ICs of high fluctuation width σz. The isolines of relative performance changes were interpolated along a grid search. No label noise was applied to generate maps (A) and (B) for methods 'Tik-SPoC' and 'NTik-SPoC', respectively. The second row reports the landscapes including a label noise level of ξn = 0.4 for both methods. Additional diamond markers in subplots (B) and (D) depict the grand average of α * for 'aTik-SPoC', which is independent of label noise. This method utilizes analytically derived values of α * and may serve as a reference for the CV-based 'NTik-SPoC'.

D NTik-SPoC aTik-SPoC B ξ n =0 ξ n =0.4 ξ n =0.4 -0.2 -0.1 0 0.1 0.2 rel. z-AUC GA Fig.
ods at N train = 900, the absolute performance z-AUC GA scales almost linearly with the amount of label noise ξ n . Referring to the relative performance change shown in (C) and (D) as well as the statistical tests in (E) and (F), they reveal that under increased levels of label noise ξ n even larger training sets can profit from regularization when compared to the unregularized SPoC. While for ξ n = 0 a relative performance increase on the GA can be found up to N train ≈ 60, for ξ n = 0.6 it increases up to N train ≈ 800. This effect is stronger for 'AS-SPoC' than for 'aTik-SPoC'. Despite not shown here, we would like to mention, that under increased label noise the performance gain of the regularized variants with larger N train can be observed also for the 'low σ z ' case, but with a lower overall decoding performance.

Optimal regularization parameter ranges

To identify suitable ranges of the regularization parameter for the CV-based methods, color-coded contour maps of relative performance changes are provided in Fig. 7. The maps show the grand average rel. z-AUC GA within the (N train , α) hyperparameter space separately for the two methods 'Tik-SPoC' (first column) and 'NTik-SPoC' (second column). Maps in the upper row summarize the performance changes in the absence of label noise (ξ n = 0) while the lower one provides these results under systematic label noise (ξ n = 0.4). The blue areas in each map mark ranges in the hyperparameter space, where a relative performance increase is In each subplot, a marker represents one of 18 subjects. Above each scatter plot, the p-value of a one-sided Wilcoxon rank sum test is given as well as the percentage of subjects for which the regularized variant outperforms baseline SPoC. Additional percentage values in brackets exclude data points located inside the grey shaded area. The latter marks a threshold criterion on z-AUC for meaningful predictions.

obtained, while "no-go" areas in red associate with a decrease of decoding quality. When comparing Fig. 7 (A) and (B), we observe that the trace norm in 'NTik-SPoC' induces a reduction of optimal α values by a few orders of magnitude as well as a larger sensitivity range compared to 'Tik-SPoC'. Both plots reveal consistently a "no-go" area towards the top right corner, which indicates, that on the grand average strong regularization is detrimental, when large training datasets without label noise are available. With additional label noise in Fig. 7 (C) and (D), the heterogeneity of the relative performance landscape increases and the "no-go" areas at the top right shift towards larger N train . In accordance with the automatic shrinkage based methods visualized in Fig. 6 we find, that the inclusion of label noise ξ n into the simulation has the effect that regularization might even be beneficial for large training sets.

For different training set sizes N train , we now compare the CV-based estimates of α with those of 'aTik-SPoC', which makes use of an analytical solution α * . The grand average of α * is plotted into Fig. 7 (B) and (D). As the analytical solution for α * [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF][START_REF] Schäfer | A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics[END_REF] is proportional to N -2 train , it should scale anti-proportional with log 10 (N train ), which in fact was observed in Fig. 7 (B). It is worth to mention, that the analytic choices of α * are not influenced by label noise -compare maps (B) and (D) -as the involved covariance shrinkage (see Eq. 9) does not make use of the label information. 5.2 Real-world data scenario

Comparison of regularized SPoC variants

In Fig. 8, the subject-wise performance comparison of all regularized SPoC variants to standard SPoC is depicted. To compare each regularized variant to its baseline, we report two different group statistics. First, the overall ratio of subjects for which the regularized variant outperforms standard SPoC is provided. Second, the values in brackets consider only those individual performances which cross a threshold of minimum meaningful performance z-AUC th = 0.59. For details on how this chance level has been determined via group analysis of predictors, we refer to [START_REF] Meinel | Pre-trial EEGbased single-trial motor performance prediction to enhance neuroergonomics for a hand force task[END_REF]. To verify if a regularized variant reaches a statistically significantly higher performance compared to standard SPoC, a one-sided Wilcoxon rank sum test was evaluated on the group level. The corresponding p-values are reported in the plot headers of Fig. 8 (A)-(E).

The following observations for the RT decoding on real-world data were made: First, in contrast to all other regularization variants, the performance changes induced by 'aTik-SPoC' are negligible small. Second, across the remaining regularization approaches we observed a tendency towards larger benefits for initially poorly performing subjects. On the group level, all regularization methods except 'Tik-SPoC' registered the majority of data points above the bisectrix. The CVbased 'NTik-SPoC' and 'ASNTik-SPoC' behave very similarly, which have been observed before on simulation data (see Fig. 5). Both approaches significantly outperform the baseline SPoC performance.

Selected regularization strengths

The regularization parameter values α f obtained on real data by the nested CVbased regularization variants across folds f are evaluated in Fig. 9 (A) and(B). Its plots should be compared with the maps for simulations depicted in Fig. 7. The median Med[α f ] across folds is shown for each subject. It's color encodes the associated z-AUC performance. The results indicate that 'NTik-SPoC' operates in smaller α ranges than 'Tik-SPoC' does, which is in accordance with the observations from the simulation in Fig. 7. For the majority of subjects, the regularization strength is outside the "no-go" areas of the simulation as α was selected by nested CV from the interval [10 -8 , 1]. For a few subjects, a large α was chosen.

As expected from simulations, this strong regularization is linked with a low absolute decoding level. The median of the analytically computed α * f across folds for 'aTik-SPoC' are presented in Fig. 9 (C). For most subjects a way smaller median regularization strength is chosen compared to the CV-based 'NTik-SPoC' method, while we observe that the analytical solution does not elicit a significant decoding improvement (see Fig. 8 (D)).

Discussion

In summary, we have proposed a set of novel regularization techniques for SPoC. We investigated their effectiveness by evaluating their performance both on simulated and on real-world datasets. Overall, 'NTik-SPoC' based on Tikhonov regularization and additional covariance normalization turned out to be the most beneficial technique.

Simulation scenario

A closer look upon the simulation results clearly shows that the regularization benefit for SPoC strongly depends on the dataset size, prevalent label noise conditions as well as on the fluctuation width of the underlying component. As a strong absolute performance variability across datasets was present in the simulation, the reported grand average performance provides a way less optimistic view than single dataset results do. Largest regularization benefit was reported for low amount of data and components with large fluctuation widths. The latter observation might be explained by the intrinsic difficulty of SPoC to recover sources of small bandpower changes.

Intuitively, additional label noise reduces the information content per data point such that the estimation of Σ z gets more demanding. Theoretically, this disadvantage could be compensated by either enlarging the training set or by adding regularization. Using the large amount of simulation data, we were able to show that under label noise conditions even larger datasets profit from regularization.

Surprisingly, in the simulation we found that 'AS-SPoC' looses performance for large datasets (especially for 'low σ z ' components) while it outperformed standard SPoC on small datasets and revealed a good performance on real-world data as well. This observation might be explained as follows: In the simulation data, the target variable is directly estimated from the EEG (IC) epoch. As such, there should be enough samples in each epoch to estimate reliably the target variable, since it was created this way. Epoch regularization might thus not be necessary here. However, for real data this might not be the case, as the target variable is not directly dependent on the EEG epoch and contains an even unknown label noise level. As such, epoch regularization might be much more useful in that case.

The direct transferability of the simulation results to real-world data is limited by three major differences: First, in real-world experiments the number of neural sources is not known a priori. Thus, a good decoding of source power typically requires the use of several components and of a regression model. Second, in realworld experiments both, label noise and the components' fluctuation widths act as latent variables and cannot directly be estimated. Third, while in the simulation we can almost perfectly recover the label information given sufficient amount of data (z-AUC > 0.9), in real-world experiments we clearly expect a decreased upper limit of the decoding performance. This is a strong indicator for the assumption, that solely bandpower information may not suffice to fully explain the labels.

Real-world scenario

Based on the real-world data, we could show that predominantly the decoding performance of initially poorly performing subjects was improved by almost all regularized SPoC approaches (except 'aTik-SPoC'). However, we cannot report a single regularization variant that systematically performed best on all subjects.

Two important aspects can be transfered from the simulation to the real-world data. First, the simulation allowed deriving an operating range of the regularization hyperparameter α for each CV-based regularization variant. When comparing these findings with the real-word data, we found that the optimal choice of the regularization intensity α for the CV-based techniques is in good accordance with the derived "no-go" areas obtained from our simulations.

Second, according to the simulation under label noise in Fig. 6, we could gain an estimate of the label noise conditions ξ n of any real-world dataset directly by comparing the absolute achievable decoding levels with the real-world decoding performances in Fig. 8. As an example, for the best performing subject of Fig. 8 (E) with z-AUC ≈ 0.78 on N train = 310 data points, the label noise level can be estimated as ξ n ≈ 0.2 according to Fig. 6 (B). Despite such estimates may not perfectly represent the ground truth, they might be beneficial for comparing data from multiple experimental paradigms e.g. in order to choose most suitable regularization strategies.

CV-based vs. analytical model selection

Overall, we introduced three CV-based Tikhonov regularization methods for SPoC (see overview in Table 3.3) and compared their performance against two variants based on automatic covariance shrinkage. Although the decoding performance of all three Tikhonov variants are on comparable levels, they strongly differ in terms of their sensitivity range for the regularization parameter. This information, however, is of great importance when it comes to choosing parameters by cross-validation. Interestingly, we found that 'NTik-SPoC' and 'ASNTik-SPoC' profit from a logarithmically scaled search space wrt. regularization parameter α while 'Tik-SPoC' could also cope with a linear scaling. We conclude that this behavior is introduced by the additional trace normalization. When comparing 'NTik-SPoC' and 'ASNTik-SPoC', the inclusion of additional LW-based shrinkage for the numerator regularization realized by 'ASNTik-SPoC' does not boost performance significantly. Accordingly, 'NTik-SPoC' seems preferable in a direct comparison due to its lower computational effort. In future work, an alternative data-driven estimation of the regularization parameter without cross-validation might be achieved e.g. by utilizing a Bayesian framework which estimates the regularization strength via expectation maximization [START_REF] Mattout | MEG source localization under multiple constraints: An extended Bayesian framework[END_REF].

Comparing both LW-based covariance shrinkage based approaches, 'AS-SPoC' seems to be the better choice compared to 'aTik-SPoC'. Three arguments support this view. First, referring to the label noise challenged simulation in Fig. 6 we found that 'AS-SPoC' profits from regularization under high label noise even for larger training sets (N train 300) while this effect was less pronounced for 'aTik-SPoC'. Second, we found that the analytically derived regularization parameter for 'aTik-SPoC' across subjects is chosen way smaller compared to values chosen by CV for 'NTik-SPoC'. For 'aTik-SPoC', the concatenation of epochs results in N s • N train sample points to estimate Σ avg . As the LW-based regularization parameter is anti-proportional to the number of samples [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF][START_REF] Schäfer | A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics[END_REF], an overly small regularization parameter is chosen, irrespectively of whether the covariance estimate did improve. Third, the analytic approach makes an i.i.d. assumption about the data. A violation thereof due to outliers might be compensated with a CV-based strategy but not by 'aTik-SPoC'. The i.i.d. assumption might also be violated for 'AS-SPoC' when the LW-based analytical solution for the trial-wise covariance matrix is challenged by autocorrelated data of a single epoch. A potential mitigation may be provided by alternative covariance shrinkage estimators that accounts for autocorrelated data as proposed by [START_REF] Bartz | Covariance shrinkage for autocorrelated data[END_REF]. Alternatively closed-form solutions for covariance shrinkage assuming elliptical distributions could also prove superior to the LW-based solution [START_REF] Chen | Robust shrinkage estimation of highdimensional covariance matrices[END_REF].

Guidance for the practitioner

Both, simulation and real-world data results strongly indicate that there is not one single regularization variant that outperforms all others. Different global parameters, such as dataset size, the noise conditions or non-stationarity in the data influence the achievable decoding accuracy.

The work by [START_REF] Engemann | Automated model selection in covariance estimation and spatial whitening of MEG and EEG signals[END_REF] reported the superiority of CVbased compared to analytical model selection in the context of spatial whitening of M/EEG data. This supports our proposal to prefer the CV-based approaches 'Tik-SPoC' or 'NTik-SPoC' over the LW-based 'AS-SPoC' method. All three methods, however, are analytically solvable by an eigenvalue decomposition and require relatively low computational effort. As they may come up with partially disjunct components, we thus propose in practice to evaluate all three variants in parallel. The final feature set should be selected in a data-driven strategy to deduce the overall most relevant oscillatory components for a given application scenario.

Conclusion

We investigated novel regularization variants for SPoC and reported their characteristics in a simulation and real-world data scenario. Initially, we applied a novel data-driven simulation framework that by design enables to generate labeled EEG datasets. The simulation delivered two main results:

First, it allowed comparing and explaining characteristics of the regularized SPoC algorithms. We could study the influence of varying training set sizes, label noise and of the bandpower fluctuation width of the neural sources of interest. On the one hand, we found that the achievable overall decoding performance decays under increased label noise conditions and smaller datasets. On the other hand, small datasets and label noise were the settings under which several regularized SPoC variants could outperform the original unregularized algorithm. As most real-world experiments come with an unknown amount of label noise, we expect that the benefits of regularization would transfer into real-world problems. Second, the simulation outcomes offered a guideline for practitioners. It proposes to tune the search for a suitable regularization parameter to a log-scaled search space. Furthermore, it indicates that the number of training data points and label noise present in the data should guide the choice of this parameter.

As an additional validation, we tested the regularized SPoC algorithms on realworld EEG data. Its outcome supported the guidelines obtained by simulation concerning the choice of regularization parameters and achievable performance improvements. We found, that individual datasets could profit strongly from single forms of regularization. As a consequence, we recommend testing several versions of regularization if decoding performance is to be optimized in practice.

While we have chosen to compare relatively simple and general regularization techniques, this work could be expanded to more sophisticated regularization strategies e.g. to realize session-to-session or subject-to-subject transfer scenarios. The presented regularization framework and the evaluation strategy using simulated and real-world datasets may pave this way.

Information Sharing Statement

The Matlab code for the proposed SPoC regularizations is accessible on GitHub under https://github.com/ameinel/regularized_SPoC. The datasets are available upon request from the authors.

Fig. 1

 1 Fig. 1 Procedure for data-driven post-hoc labeling of arbitrary pre-recorded EEG signals. (A) Processing pipeline to extract independent components (ICs). (B) For each IC and epoch e, the log-bandpower average of the epoch serves as a ground truth label ztrue(e). (C) Distribution of ztrue over all epochs of an exemplary IC. Its bandpower fluctuation width is described by σz.

Fig. 2

 2 Fig. 2 Simulation data: scatter plot relating the fluctuation width σz of each selected ICA component to their baseline SPoC performance z-AUC ref for non-reduced datasets with 1000 epochs. Based on its σz-distribution, the dataset of each IC was labeled into one of three classes, defined by the quartile thresholds Q 25 and Q 75 .

Fig. 3 Fig. 4

 34 Fig.3Simulation results: influence of regularization strength α onto the decoding accuracy of three SPoC variants regularized via CV. The grand average performance z-AUC GA is reported in the top row (subplots (A) and (B)), while subplots (C) and (D) in the lower row report the angle θ GA between the estimated highest ranked and the ground truth filter as evaluation score. Subplots in the left and right columns differ in the number of training data points (epochs) used for SPoC decoding. Results are reported for the class 'σz -high'.

Fig. 5

 5 Fig. 5 Simulation results: influence of training set size and fluctuation width upon decoding performance of optimal regularization strength α * . The top row depicts the grand average absolute performance of five regularized SPoC variants for ICs that either have low (A) or high (B) bandpower fluctuation width. The middle row depicts performance increase or decrease of the five regularized methods relative to the baseline SPoC method without any regularization and again separately for IC's of low (C) and high (D) fluctuation width. Subplots (E) and (F) reveal color-coded points for each training set size where the regularized variant significantly outperformed the baseline method (Wilcoxon signed rank test with p < 0.05).

Fig. 6

 6 Fig.6Simulation results: interaction between label noise level ξn and dataset size N train . A level of ξn = 0 states the absence of label noise. All curves report the grand average results for ICs belonging to the 'high σz' class. Subplots (A) and (B) in the top row provide the absolute grand average performances for 'aTik-SPoC' and 'AS-SPoC', while the middle row depicts relative performance changes. The dots in (E) and (F) indicate configurations, for which the regularized variant significantly outperformed the baseline method (Wilcoxon signed rank test with p < 0.05).

Fig. 8

 8 Fig.8Real data: scatter plots (A)-(E) compare the performance of different regularized SPoC variants with the unregularized baseline method SPoC. In each subplot, a marker represents one of 18 subjects. Above each scatter plot, the p-value of a one-sided Wilcoxon rank sum test is given as well as the percentage of subjects for which the regularized variant outperforms baseline SPoC. Additional percentage values in brackets exclude data points located inside the grey shaded area. The latter marks a threshold criterion on z-AUC for meaningful predictions.

Fig. 9

 9 Fig. 9 Real data: median regularization strength across the 10-fold chronological crossvalidation for each dataset as a function of the training set size N train color coded by the achieved z-AUC decoding performance.

Table 1

 1 Overview on introduced regularization variants for SPoC using two model selection schemes: cross-validation (CV) and based on the Ledoit & Wolf shrinkage estimator (LW). For Σ(e) it is explicitly marked if regularization is applied to the numerator (N) and/or to the denominator (D) of the objective function. The checkmarks refer to covariance normalization.

	Method	Regularization Target Normalization
		Σ(e)	Σavg	Σ(e), Σavg
	Tik-SPoC	-	CV	-
	NTik-SPoC	-	CV
	ASNTik-SPoC	LW: N	CV
	AS-SPoC	LW: N,D	-	-
	aTik-SPoC	-	LW	-

To simplify the notation, epoched data X(t, e) will further on be written as a matrix X(e).

Acknowledgements This work was fully supported by BrainLinks-BrainTools Cluster of Excellence funded by the German Research Foundation (DFG), grant number EXC 1086. For the data analysis, the authors acknowledge support by the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) by grant no. INST 39/963-1 FUGG. Fabien Lotte received research support from the French National Research Agency with the REBEL project (grant ANR-15-CE23-0013-01) and the European Research Council with the BrainConquest project (grant ERC-2016-STG-714567). For parts of the data analysis, the Matlab-based BBCI toolbox was utilized (Blankertz et al, 2016). The authors declare that they have no conflict of interest. This is a post-peer-review, pre-copyedit version of an article published in Neuroinformatics (Springer). The final authenticated version is available online at: http://dx.doi.org/10.1007/s12021-018-9396-7