Optimal Motor Unit Subset Selection for Accurate Motor Intention Decoding: Towards Dexterous Real-Time Interfacing

Objective: Motor unit (MU) discharge timings encode human motor intentions to the finest degree. Whilst tapping into such information can bring significant gains to a range of applications, current approaches to MU decoding from surface signals do not scale well with the demands of dexterous human-machine interfacing (HMI). To optimize the forward estimation accuracy and time-efficiency of such systems, we propose the inclusion of task-wise initialization and MU subset selection. Methods: Offline analyses were conducted on data recorded from 11 non-disabled subjects. Task-wise decomposition was applied to identify MUs from high-density surface electromyography (HD-sEMG) pertaining to 18 wrist/forearm motor tasks. The activities of a selected subset of MUs were extracted from test data and used for forward estimation of intended motor tasks and joint kinematics. To that end, various combinations of subset selection and estimation algorithms (both regression and classification-based) were tested for a range of subset sizes. Results: The mutual information-based minimum Redundancy Maximum Relevance (mRMR-MI) criterion retained MUs with the highest predicative power. When the portion of tracked MUs was reduced down to 25%, the regression performance decreased only by 3% (R2=0.79) while classification accuracy dropped by 2.7% (accuracy = 74%) when kernel-based estimators were considered. Conclusion and Significance: Careful selection of tracked MUs can optimize the efficiency of MU-driven interfacing. In particular, prioritization of MUs exhibiting strong nonlinear relationships with target motions is best leveraged by kernel-based estimators. Hence, this frees resources for more robust and adaptive MU decoding techniques to be implemented in future.


I. INTRODUCTION
S URFACE electromyography (sEMG) is one of the most prominent means for noninvasively accessing user intention in human-machine interfacing (HMI).Specifically, in the control of upper-limb prosthetics, sEMG-operated or myoelectric systems remain the only class of powered devices that are clinically available.Here, state-of-the-art approaches have typically relied on pattern recognition and regression algorithms to decode user intention from the time-domain features of multiple recorded channels [1], [2].While suitable for devices offering basic functionality, the maximal throughput using such global features is limited by the poor specificity of the sEMG signal [3].As sEMG is an interferent signal comprised of superimposed action potentials originating from spatially distributed sources inside the muscle [4], the effects of spatial filtering (eg.volume conductor) and amplitude cancellation can confound the relationship between global amplitude characteristics of the signal and the underlying neural drive [5], [6], [7], [8].Conversely, knowledge on the firing times of motor neurons innervating these muscles would provide access to the most basic unit of neural drive responsible for instigating force generation, potentially facilitating more intuitive and dexterous interfacing [9], [10].
While the discharge times of individual MUs can be automatically obtained via spike sorting of intramuscular recordings [11], [12], the invasive nature of the signal acquisition presents an added obstacle to clinical feasibility.Alternatively, the use of high-density sEMG (HD-sEMG) recordings in conjunction with convolutive blind source separation (BSS) [12], [13], [14] have been shown to accurately extract MU activity during voluntary contractions.Such decomposition methods applied to HD-sEMG recordings are capable of identifying larger populations of active MUs [10] and are applicable to a variety of musculatures with components that serve both shared and specific motor functions [15], [16].
Online extensions of such decompositions techniques have since been developed to operate in real-time, extracting MU activity from smaller time windows of HD-sEMG (50-250 ms) [17], [18], [19], [20], [21].The majority of such methods share a common process in the pre-identification of MUs and their respective separation vectors (MU filters) from training data via batch (offline) decomposition.Here, MU filters refer to spatiotemporal projections that compensate for the action potentials of the identified MUs.Their direct application in decomposing newly acquired HD-sEMG has received some validation regarding isometric contractions [19], [22], [23].More recently, the incorporation of deep learning has also been proposed [24], [25].Coined as hybrid approaches [10], BSSbased decomposition is still relied upon for initial extraction of MU activities subsequently used in the training of a deep network.Such an approach has been shown to be more robust in online decomposition [25].For addressing the non-stationary conditions of non-isometric contractions, means of adapting the MU filters in real-time have also been proposed [17], [20].
In the context of HMI, online decomposition is used to extract motoneuronal activity to support the estimation of motor intention.In [18], real-time control of a hand prostheses over 4 motion classes was demonstrated using the decomposed activity of five MUs.In [21], control over two degrees-of-freedom (DoFs) driven by MU activity extracted from the forearm was shown to outperform interfacing with conventional global sEMG features.Online experimentation has also shown subjects to be capable of learning novel MU recruitment orders given appropriate feedback, expanding the potential control bandwidth that one muscle can support [26].
To facilitate HMIs with greater functionality, it is beneficial to track large populations of MUs.However, a single decomposition step to process HD-sEMG comprised of multiple motor tasks tends to identify inadequate sources for accurate estimation of motor intention.Thus, recent studies have highlighted the need to perform decomposition in a task-wise manner where segments of the HD-sEMG signal pertaining to distinct motor tasks are decomposed separately [27], [28].To address the redundancy of identifying a single MU multiple times, sources that share a significant portion of discharges (20-30%) over the full repertoire of training motions are deemed to be equal.The removal of duplicate MU extractions, however, does not ensure optimal efficiency.Given that the hardware capabilities of wearable HMIs tend to be more constrictive than the machines employed in laboratory-based interfacing experiments, it is practical to consider the prioritization of MUs that benefit estimation performance the most.
To that end, feature selection algorithms may be considered.These have typically been used to improve the efficiency of machine learning in applications involving high data dimensionality and can be categorized as either wrapper, embedded, or filter-based methods [29], [30], [31].Whereas wrapper and embedded methods require the output of the learned model to inform the search of features, filter methods are model-agnostic and rely on the ranking of candidates based on some measure of dependency with the target variable.The simplest measures of dependency include Pearson's or Spearman's correlation coefficients [32].Here, since the target kinematics are multivariate, Canonical Correlation Analysis (CCA) can be utilized instead [33], [34].Another popular measure is mutual information (MI) which is capable of detecting nonlinear dependencies [35], [36], [37].The union of the most informative features, however, does not yield the most informative subset [38].Thus, the goal to maximize joint dependency through an exhaustive search of all possible feature combinations is computationally intractable for most cases.Therefore, various forward-searching algorithms have been developed [32], [39], [40], amongst which the minimal Redundancy Maximal Relevance (mRMR) framework [41], [42] has received widespread usage [36], [43], [44], [45].
This work, therefore, investigates the use of feature selection techniques in finding the optimal MU subsets to be tracked.In the context of wearable HMIs, limited computational resources have to be shared between the tasks of online decomposition and motor intention estimation.Given that the machine learning algorithms employed to accomplish the latter can have wildly varying complexities, it is worthwhile to also investigate compromising the quantity of tracked MUs in order to accommodate more computationally expensive estimation algorithms.The incorporation of subset selection to an MU-driven estimation pipeline was analyzed in cross-validation format using HD-sEMG from 18 motor tasks pertaining to the single and pair-wise combined activations of three wrist/forearm DoFs.From the train data set, MUs were identified via task-wise batch decomposition and the subsequent MU subset selection step was performed with different selection criteria.The mRMR scheme was employed using CCA and MI as measures of dependency (mRMR-CCA and mRMR-MI, respectively).For comparison, their Maximal Relavance counterparts (MR-CCA and MR-MI, respectively) along with two arbitrary selection schemes based on random selection and maximal MU activity were also tested.Assessment of the selection criterion was made based on the changes to open-loop accuracy with different estimation algorithms as MU subset sizes were reduced.The estimation algorithms were driven by the spiking activities of the subset MUs, extracted using an online decomposition approach applied to the test set data.Results using three regression algorithms: linear regression (LR), multilayer perceptron (MLP) and kernel ridge regression (KRR), and three classification algorithms: linear discriminate analysis (LDA), k-Nearest Neighbors, (kNN) and kernel support vector machine (KSVM) were obtained.

A. Subjects
Eleven healthy subjects, seven male and four female, all right-handed, aged 26-34, participated in the experiment.The study was approved by the local ethical board of Aalto University.In accordance with the Declaration of Helsinki, all participants gave their written informed consent prior to the experiments.

B. Experimental Protocol
HD-sEMG was recorded from each subject's dominant side with three 8 × 8 electrode matrices spaced evenly around the bulk of the forearm.All 192 channels were sampled at 2048 Hz by a benchtop bioamplifier (OT Bioelettronica, IT).Wrist joint angles, along with pronation and supination, were recorded at 80 Hz with three wireless Inertial Measurement Units (IMUs) (Xsens Technologies B.V, NL) attached to the posterior sides of the upper-arm, mid-forearm and hand.
Throughout the experiment, subjects were seated upright with their recorded limb relaxed by their side.Three repetitions of single-DoF motions pertaining to wrist flexion/extension (FL/EX), radial/ulnar deviation (RD/UD) and forearm pronation/supination (PR/SU) were recorded.Each repetition followed a trapezoidal activation profile of 2 s ramp time and 10 s plateau time.In addition, three repetitions of each pair-wise combination of the previously recorded motions were obtained, resulting in a dataset of 18 individual motor tasks.The recordings and subsequent processing were carried out using a custom Matlab R2021b (MathWorks Inc, MA, USA) framework.
Three-fold cross-validation was conducted where 2 repetitions from each motor task were aggregated to form the training data from which MUs were initially identified, subset selection was conducted, and estimators were trained.For estimators with tunable hyper-parameters, an inner two-fold optimization step was conducted within the training set.The remaining repetition from each motor task formed the test set used to perform simulated real-time decoding where the HD-sEMG was processed in sliding windows of 100 ms length which advanced in increments of 50 ms.

C. Decomposition Algorithms 1) sEMG Mixture Model:
The decomposition methodology employed in this work is based on a convolutive mixture model for sEMG generation shown in (1) [12]: where z(k) is the column vector of M multichannel HD-sEMG measurements (observations) at sample k, τ (k −l), the column vector of delta train values from N identifiable sources and ε(k), additive noise which also encompasses signal contributions from non-identifiable sources.A(l) is the matrix of action potentials where each column encodes the action potential of a MU at l samples after an impulse.While Equation (1) describes a convolutive mixture of finite-length impulse response filters, the observation vector can be artificially extended with its time-delayed versions.This transforms the de-mixing problem to that of instantaneous spatiotemporal filters.Accounting for this extension, the mixture model becomes: where the tilde accent indicates extension by a factor of E such that z(k) is now a column of M • E elements while Ã now encompasses MU action potentials up to a duration of E samples.
2) Initialization With Task-Wise Batch Decomposition: To enable the estimation of MU activities in real-time, initial batch decomposition of training data is necessary for pre-identification of active sources and their corresponding MU filters.First, the HD-sEMG signals allocated for system training are partitioned by motor task such that separate batch decomposition steps can be conducted on each data block, Z c .This has been shown to maximize the number of task-specific MUs [27].Here, variables derived from a specific motor task are denoted by subscript c ∈ {1, . . ., C} with C = 18 in this work.
The batch decomposition algorithm employed follows that described in [12].Here, a brief overview is provided while the interested reader is referred to the original publication for further exposition.The algorithm can split into the Sphering, Extraction, Refinement and Acceptance steps.In the Sphering step, the HD-sEMGs signals are extended, zero-meaned and a zero-phase component analysis (ZCA) whitening transform is calculated such that transformed observations will be decorrelated.This preprocessing creates the necessary conditioning for the following steps to converge to valid MU filters and sources.Candidate MU filters are first sequentially identified by execution of Extraction and Refinement steps while the final set of filters are chosen in the Acceptance step.In the Extraction step, a separation vector is obtained by iterating a fixed-point algorithm which optimizes source signal sparsity.These iterations are interlaced with Gram-Schmidt orthogonalization against past solutions to prevent convergence to the same sources.The second Refinement step iteratively improves the quality of the projection and source signal in terms of a silhouette measure (SIL) which is a measure of accuracy of the separation [12].Previous studies have shown that sources with a relatively high SIL value are often physiologically accurate in both isometric and isotonic conditions [46].To convert the continuous source signal to a delta train, the signal is first rectified by element-wise power of 2 and peak detection algorithm is applied.Resultant peaks are then clustered by the kmeans++ algorithm to differentiate MU spikes from noise peaks.Finally, the Acceptance step validates sources that satisfy a minimum SIL threshold (0.90 as per [12]), discarding duplicate sources that have been repeatedly extracted as timeshifted versions.These are identified as sources with over 30% common discharges (±0. 3) Online Decomposition Algorithm: The online decomposition algorithm, which is based on the main Extraction step of [12], as recently implemented by [18], operates on far smaller windows of HD-sEMG, under strict time constraints.By reapplying preprocessing transforms and MU filters retained from batch decomposition (Section II-C.2),MU activities can be more efficiently estimated: where Zwin (t) is the t th window of extended HD-sEMG signals, S c (t) = [s c 1 (t), s c 2 (t), . . ., s c Nc (t)] ⊺ is the corresponding matrix of estimated activities for MUs identified from motor task c and 1 is a vector of ones of appropriate length.Fig. 1A gives an overview of how these retained preprocessing transforms and MU filters support the online algorithm when executed on a newly acquired window of HD-sEMG.
Spike detection is then performed on each of the rectified source signals where detected peaks are classified as either spikes or noise based on distances from the stored centroids.If the SIL of the spike train meets the minimum threshold of 0.93, the decomposition quality is deemed to be satisfactory and the spike events are accepted, otherwise, all of that source's spike events within the process window are discarded.The output of the online algorithm is a vector of non-negative integers representing decomposed spike counts, x c (t) = x c 1 (t), x c 2 (t), . . ., x c Nc (t) ⊺ .Each element corresponds to the number of discharges detected for a MU (identified from the c th motor task) within the observations stored in Zwin (t).The decomposed spike count feature has been used in past studies linking MU discharge patterns and motor intention estimation [28], [47].

D. MU Subset Selection
A full feature matrix is first constructed by extracting the activities of all identified MUs over the full training data set.This is achieved by applying the online decomposition method such that the activities of MUs initially identified from individual motor tasks are extracted over the entire repertoire of training movements (Fig. 1(b)).
For formulating the subset selection methods, it is convenient to treat the full collection of MUs as a set of random variables, }, which can be paired with recorded kinematics or class labels.Here, the notation of '(t)' has been omitted from x i for convenience.
The subset selection step now identifies a subset, S, based on some optimality criterion or selection scheme where the cardinality of the subset is bounded, |S| ≤ ξ |F|, 0 ≤ ξ ≤ 1. Future deployment of the online decomposition algorithm would then only need to extract the activities of MUs within S. In this study, ξ = 1, 0.5 and 0.25 was tested.
1) Canonical Correlation Analysis: A basic approach is to prioritize MUs whose activities exhibit the strongest dependency with the recorded kinematics, defined as a random vector y = y 1 , . . ., y D ⊺ , where D is the number of DoFs to be estimated.Given two random vectors q and v with Q and V elements, respectively, CCA finds respective projections α and β such that the Pearson correlation between the projected samples (canoncical variables) are maximized: where qv is the cross-covariance matrix of q and v while qq and vv are their respective auto-covariance matrices [33].The optimal subset under the MR-CCA criterion is thus found by satisfying the following condition: max 2) Mutual Information: An alternative criterion is to consider the mutual information between the individual MU activities and the class label, ℓ, which annotates the motion type at a particular time and takes values from {0, 1, 2, . . ., C}.Given paired discrete random variables, q and v, with values occupying Q and V, respectively, their mutual information is calculated as the reduction in uncertainty that knowledge of one affords to the value of the other: where H(•) is the Shannon entropy of the argument variable and p(•) is the event probability.For this work, the marginal and joint distributions can be exactly represented by histograms with integer-centered bins.The optimal subset under the MR-MI criterion is thus found by satisfying (7) [41]: 3) Maximal Relevance and Minimal Redundancy: To reduce the redundancy within S, an incremental approach can be taken where, at each step, candidate MUs are also penalized by their degree of dependency with the MUs already selected.Incorporating this, the step-wise optimization criterion for mRMR-CCA and mRMR-MI are expressed as ( 8) and ( 9), respectively, max max Here, ( 9) is equivalent the Mutual Information Quotient (MIQ) scheme presented in [42] while (8) is the variant which considers relevancy and redundancy in terms of CCA [34].4) Selection by Random and by Maximum Activity: For comparative purposes, two naïve selection schemes were also tested.One selects MUs by random while the other method prioritized MUs that had the highest spike counts during the training movements.

E. Motor Intention Estimation Algorithms
Three regression-based kinematics estimation methods (LR, MLP, KRR) and three classification-based motor-task estimation methods (LDA, kNN, KSVM) were tested.Prior to being passed to the estimation algorithms, the decomposed spike counts of selected MUs were smoothened by a 4th order moving average filter.Estimated joint angles from regression algorithms were further smoothened by a 6th order moving average filter.
1) Regression: In LR, a linear least squares mapping between y and the smoothed neural signals of MUs in S is established by the Moore-Penrose pseudoinverse method [48].
For MLP-based estimation, single hidden-layer feedforward networks using the tanh activation function are trained via the Levenberg-Marquardt backpropagation algorithm.Similar to [49], each DoF is estimated by a dedicated network while the optimal hidden-layer node counts are obtained via grid search.To compensate for randomized initialization, training processes are rerun 10 times and the model with the lowest train error is used.
With KRR, a mapping is formed by the inner products between samples projected to a higher dimensional kernel feature space.The radial basis function (also referred to as the Gaussian kernel) is employed.Two hyperparameters, the ridge regularization scale and the spread of the Gaussian kernel, are optimized via grid search.
2) Classification: In LDA-based classification, the predicted class is given by the maximum likelihood estimate (MLE) where classes-conditional densities are assumed as Gaussians with equal covariance matrices.Despite the simplistic modelling, the method has been found to yield practical results in past studies and is often used as a benchmarking tool for gauging more sophisticated techniques [2], [50].
kNN incorporates all training samples as the classification model.Given a queried input, the majority class amongst proximal (by Euclidean distance) training samples is given as the estimate.The optimal numbers of proximal training samples to consider are found via grid search.
KSVM forms a decision boundary between two clusters by maximizing the region between two marginal hyperplanes identified by support vectors.As KSVM is strictly a binary classifier, extending its use to a multi-class scenario can be done via a One vs.One topology where C(C − 1)/2 classifiers are formed, one for each pairwise comparison, and the estimated class is given by majority vote.As with KRR, the Gaussian kernel is used to permit nonlinear decision boundaries.Two hyperparameters, the penalization weight of training misclassifications, and the spread of the Gaussian kernel, are optimized by grid search.

F. Performance Metrics and Statistical Analyses
The performance of regression-based estimation algorithms were measured in terms of the R 2 metric [49] while categorical prediction of motor tasks was gauged by classification accuracy.To detect significant differences between the use of full MU sets and subsets of various sizes and composition, repeated-measures ANOVA (RM-ANOVA) was conducted.Data from each estimation algorithm was analyzed separately with samples taken from the averaged cross-validation results of each subject.Normality of datasets were verified by Shapiro-Wilks tests.Detection of statistically significant differences was followed by Bonferroni-corrected pairwise comparisons.For direct comparisons between subset selection schemes, focused RM-ANOVA was also conducted at each subset size.Similarly, detection of significant effects were followed by Bonferroni-corrected pairwise comparisons.Statistical analyses were conducted with SPSS Statistics 28 (IBM, Armonk, NY, USA) with significance thresholds set at P < 0.01.

III. RESULTS
The number of MUs decomposed from each motor task is given in Table I.On average, 20.3 ± 8.8 viable MUs were extracted via batch decomposition from two training repetitions of each motor task.On the other hand, application of the batch decomposition algorithm on HD-sEMG data encompassing all motor tasks yielded only 4.8 ± 2.8 MUs which is insufficient for prediction of all recorded motor tasks.Table II gives the total number of MUs extracted from each subject in each cross-validation fold.On average, 365.9±103.0MUs were extracted.This also corresponds to the average size of F.
Fig. 3(a-d) shows the baseline estimation performances from using the full sets of extracted MUs.A weak positive correlation between the number of MUs in F and performance is shown by linear regression analyses that give R-squared

TABLE II TOTAL NUMBER OF MUS DECOMPOSED PER CROSS-VALIDATION FOLD FROM EACH SUBJECT
values of 0.329 and 0.275 for regression and classificationbased estimators, respectively.The cross-validation average R 2 values for each DoF are 0.79 ± 0.08, 0.72 ± 0.08 and 0.76±0.09for DoF 1, 2 and 3, respectively.On average, the motor task with the highest classification accuracy is FL+UD (0.90±0.08) while the lowest is RD (0.62±0.18).Estimation performances from different subset selection criteria and subset size combinations are displayed in Fig. 4(a).Shapiro-Wilk's tests confirmed the normality of all results.Statistically significant differences were detected amongst the subset size/selection criteria combinations for all estimation algorithms.From the pair-wise comparisons, dependencybased selection criteria (MR/mRMR-CCA/MI) prevented statistically significant performance drops at 50% subset size  across all estimation algorithms except LR.At 25% subset size, only the MI-based selection schemes (MR/mRMR-MI) prevented significant performance reductions with the same estimation algorithms.From the focused RM-ANOVA, significant differences between the selection criteria were detected at every subset size across all estimators.Fig. 4(b) shows the results from the Bonferonni-corrected pairwise comparisons.Across the 12 separate analyses, no significant differences were detected amongst the dependency-based selection criteria.However, mRMR-MI accrued the highest counts of statistically significant performance gains over the naïve selection schemes of maximum activity (8/12) and random selection (12/12).The average R 2 and classification accuracies obtained at 25% subset sizes are shown in Tables III and IV, respectively.Overall, subsets selected by mRMR-MI provided the lowest reduction in performance (−3.5% R 2 and −2.1% classification accuracy) while randomized selection performed the worst (−14.8%R 2 and −11.6% classification accuracy).On average, consideration of redundancy (mRMR) improved upon MR, while MI, as a measure of dependency, proved to be more appropriate than CCA for this application.

IV. DISCUSSION
In this work, we have proposed an approach to MU-based HMI designed to scale efficiently with the number of supported target functions.The main components of the method are the initial task-wise extraction of MUs and the subset selection step.Pseudo-online testing was conducted to gauge the feasibility of the approach and to identify an ideal selection criterion based on the responses of various estimation algorithms to the different subset compositions.

A. Motor Task-Wise MU Identification
The task-wise batch decomposition scheme employed here is similar to the 'segment-wise decomposition' scheme in [27].There, the approach was verified, on both synthetic and experimental data, to extract more MUs from HD-sEMG activity pertaining to multiple motor tasks.This is in agreement with the results obtained here.However, a number of past studies on MU-based interfacing [18], [21], [47] utilize a different approach.Namely, the entire population of MUs used for motion estimation were extracted from a single decomposition process conducted on the full aggregation of training motions.Of note is that the number of supported motions in such studies were relatively low, not exceeding four motion classes.Indeed, decomposition of multi-task signals tend to yield less sources which limits the performance of motor estimation [27].Conversely, in [28] and [51], separate decompositions were conducted to extract MUs from different motions.In both cases, the rationale for decomposing MUs from each task separately was to maximize the number of identified MUs.The former study extracted MUs from the antagonistic activations of three DoFs (6 motor tasks) while the latter extracted MUs from single and combined finger flexions (10 motor tasks).Here, we further extend the application of this approach by identifying and tracking MUs from 18 motor tasks which can be used to facilitate to a highly dexterous HMI.

B. MU Subset Selection
By identifying an optimal subset of MUs, the computational demand of online decomposition is reduced as the activities of redundant MUs are no longer tracked.Moreover, the model complexity of the downstream estimation algorithms are also reduced due to the lower input feature dimensionality.In that regard, it is important to distinguish the suitability of feature selection over dimensionality reduction techniques as the latter operates by learning projections to subspaces which would still require the extraction of the full set of MU activities.
The mRMR criterion was shown to be effective in selecting the best MUs given the large performance disparities compared to the naïve selection methods.This difference became even more pronounced at smaller subset sizes (Fig. 4(a)).The use of MI as a measure of dependency was also shown to outperform CCA in almost all cases (Table III and Table IV).This may be attributed to the capability of MI in capturing both linear and nonlinear relationships.As such, nonlinear kernelized CCA may then offer improved feature ranking over CCA [33], [45].While the subset portions tested here are arbitrary, subset sizes in real-world use cases would be informed by computation resource constraints.
One by-product of the task-wise decomposition scheme is that the same MU can be identified multiple times from the separate decomposition processes.This is expected as many motor tasks supported by MU-based interfacing share recruitment of the same muscles.In other works, duplicate extractions have been identified by either shared spike events or similarity in action potential profiles [27], [28], [51].Here, the consideration of redundancy inherently remedies repeated MU extractions as duplicate source activities will exhibit high dependencies with one another.Hence, the mRMR criterion automatically addresses the issue of duplicate extractions while also ensuring that the optimal copy of the MU filter for the application is preserved.

C. MU-Driven Estimation Algorithms
In this study, the proposed initialization pipeline for MU-based HMI was paired with six different estimation algorithms from two myocontrol paradigms: classification and regression-based control.The latter natively supports 'Simultaneous and Proportional Control' which is considered a more intuitive mode of interfacing as it permits fluid translation between motor tasks and activation intensities [52].Regardless, recent studies in MU-driven interfacing have incorporated estimators from both paradigms.For instance, [18], [53], [54] utilized linear SVM while [16] used LDA, [21], [51] used LR and in [28], MLPs with two hidden layers were used.Of note is that kernel methods have not been employed despite their favorable performance in conventional myocontrol studies [49].
Here, we directly compare these estimation algorithms, where again, kernel-based methods (KRR and KSVM) yielded the highest baseline estimation accuracies.One major draw-back of kernel methods is that they require calculation of the kernel function between the input and all training data points.This can be computationally cumbersome and hard to execute under the strict time constraints for HMI applications.Hence, MU subset selection may be used to make the deployment of these estimation algorithms more feasible as online decomposition and kernel function overheads are reduced.The results show this is indeed worthwhile as KRR trained with 25% of MUs selected by mRMR-MI outperformed LDA and MLP when those algorithms had access to the full set of MUs (Table III).Likewise, KSVM trained with 25% of MUs selected by mRMR-MI performed equally with LDA and outperformed kNN when trained with the full set of MUs (Table IV).

D. Limitations and Future Work
In this work, only pseudo-online analyses was conducted though it is well known that the relative real-time control performances between estimation algorithms can differ from their open-loop accuracies [50], [55].This is largely due to feedback-enabled user adaptation.Thus, one may expect the relative performances between algorithms to change.Moreover, it is still unknown if an interaction exists between user-learning and different MU subset compositions and sizes.
Batch and online decomposition algorithms were applied to HD-sEMG recorded over dynamic contractions even though the decomposition algorithms rely upon an assumption of stationarity.While past works have similarly applied such methods over dynamic contraction data [51], no quantitative data regarding the decomposition accuracy of such cases currently exist.Here, a higher SIL threshold (0.93) was employed in the acceptance of spikes to compensate for the greater variances between the data used for initialization of MU filters and the data presented in subsequent testing phases.Meanwhile, recent works have focused on adaptive MU filters that can compensate for the changes in action potentials due to dynamic contractions [20].Given that any MU filter adaptation increases the complexity of online decomposition, the MU subset optimization work here becomes more relevant.

V. CONCLUSION
We have further demonstrated the feasibility of MU-based HMI by extending the task-wise decomposition scheme to extract MUs that facilitate estimation of 18 motor tasks.Given the strict time constraints of HMI applications and the limited computational capabilities available on wearable devices, optimization of the online decomposition and motor estimation steps is warranted.To that end, mRMR-MI was shown to be a highly effective criterion for the selection of MU subsets that exhibit a high collective predictive power.In turn, the activities of these MU subsets were best leveraged by kernel-based algorithms for the estimation of motor intention.
5 ms timing tolerance) and the version with the lowest coefficient of variance in inter-spike intervals (COV isi ) is retained.The accepted MU filters, B c = b c 1 , b c 2 , . . ., b c Nc , where N c is the number of MUs extracted from motor task c, are stored for subsequent online extraction of source signals.For online spike detection, the spike and noise centroids for each source, hi c n and lo c n , respectively, are also stored: c = {(hi c n , lo c n ), n = 1, . . ., N c }. Preprocessing transforms, such as the subtractive observation means, E[z c (k)], and whitening matrix, W c , are also retained for online application.
1. (a) Schematic for batch and online decomposition techniques showing the parameters that are transferred.(b) Initialization process of the proposed MU-based interfacing.

Fig. 2 .
Fig.2.Example of subset MU activity extracted over unseen HD-sEMG data pertaining to all 18 motor tasks which include single and combined DoF activations.The subset was identified via mRMR-MI with a size of 25%.The corresponding estimated kinematics obtained using KRR are shown at the bottom.Estimated DoFs pertaining to wrist flexion/extension, radial/ulnar deviation and forearm pronation/supination are shown as solid blue, red and yellow lines, respectively.The recorded ground-truth kinematics are shown as dotted lines.This particular example has an R 2 of 0.86.

Fig. 3 .
Fig. 3. Baseline performances from using the full sets of extracted MUs.(a) Kinematics estimation performances plotted against the number of MUs in F. Results from all estimators are aggregated and displayed with each cross representing the result from one cross-validation fold.The mean result from each subject is shown as a circle.Horizontal and vertical error bars indicate one standard deviations in total MU count and R 2 , respectively.The line-of-best fit is shown as a gray dashed line.(b) Classification accuracy displayed in the same format.(c) The average R 2 across each DoF.(d) The average classification accuracy across each motor task.

Fig. 4 .
Fig. 4. (a) Violin plots of estimation performance with different estimation algorithms, subset selection criteria and subset sizes.Each dot represents one sample (subject average).Light shaded areas represent probability density functions estimated by kernel density estimation while darker shaded blocks show the 1st-3rd quartile range and medians are indicated by black notches.Statistically significant differences with corresponding full set is indicated by asterisks.(b) Pairwise-comparisons focused at each subset selection size.The color shading of each off-diagonal cell corresponds to performance differences calculated as the subtraction of the column scheme from the row scheme.Significant differences are marked by asterisks.

TABLE I NUMBER
OF DECOMPOSED MUS FROM EACH MOTOR TASK

TABLE IV CLASSIFICATION
ACCURACY WITH FULL SET OF MUS AND 25% SUBSET SIZE.(RELATIVE DROP SHOWN IN BRACKETS)