Robust extraction of basis functions for simultaneous and proportional myoelectric control via sparse non-negative matrix factorization

Objective. This paper proposes a novel simultaneous and proportional multiple degree of freedom (DOF) myoelectric control method for active prostheses. Approach. The approach is based on non-negative matrix factorization (NMF) of surface EMG signals with the inclusion of sparseness constraints. By applying a sparseness constraint to the control signal matrix, it is possible to extract the basis information from arbitrary movements (quasi-unsupervised approach) for multiple DOFs concurrently. Main Results. In online testing based on target hitting, able-bodied subjects reached a greater throughput (TP) when using sparse NMF (SNMF) than with classic NMF or with linear regression (LR). Accordingly, the completion time (CT) was shorter for SNMF than NMF or LR. The same observations were made in two patients with unilateral limb deficiencies. Significance. The addition of sparseness constraints to NMF allows for a quasi-unsupervised approach to myoelectric control with superior results with respect to previous methods for the simultaneous and proportional control of multi-DOF. The proposed factorization algorithm allows robust simultaneous and proportional control, is superior to previous supervised algorithms, and, because of minimal supervision, paves the way to online adaptation in myoelectric control.


Introduction
The surface electromyographic (EMG) signal has been used as the control source for multifunction prostheses for decades. In this study we will specifically focus on myoelectric control for wrist function replacement, given its clinical relevance.
Currently, state-of-the-art wrist-control prosthetic hands employ machine learning techniques to extract the neuromuscular control information from the surface EMG. Researchers have proposed several control schemes based on pattern recognition techniques [1], and have also made progresses in enhancing the control accuracy by spatial filtering [2], sensory feedback [3,4], and approaches insensitive to shifting and rotation of the electrodes [5]. Moreover, targeted muscle reinnervation is currently a successful clinical solution for prosthetic control in high-level amputations, that can be combined with pattern recognition [6][7][8]. However, with respect to the classic sequential control scheme of pattern recognition, natural limb movements consist in the simultaneous activations of multiple DOFs.
For simultaneous and proportional control, Jiang et al have analyzed the muscle synergy model and utilized the classic non-negative matrix factorization (NMF) algorithm [9] to blindly solve a multi-channel EMG factorization problem that provided continuous control signals. The NMF scheme adopts a DOF-wise training strategy, according to which simultaneous control of multi-DOF is obtained by linear summation of single-DOF activations. This allows simultaneous and proportional EMG control without multi-DOF training data [9,10]. Despite the promising results of this scheme, there are still shortcomings. The most relevant limitation of the NMF approach for myoelectric control is that it requires a DOF-wise calibration phase during which the user should activate single DOFs in a predefined order, e.g. first activate DOF1 for several trials, and then DOF2. This procedure may be inconvenient and, more importantly, the prosthetic user may not be able to exactly activate an individual DOF. If the activations in the training phase are not precisely of single DOFs, the DOF-wise based NMF scheme would generate an inadequate synergy matrix, leading to poor control. This behavior is due to the fact that the classic NMF cannot generate a sparse enough factorization [11,12].
In this study, we propose an alternative approach to the classic NMF scheme, our scheme can find out all the bases in the basis matrix robustly and concurrently in a quasiunsupervised way by factorizing the multi-channel surface EMG signals with the inclusion of sparseness constraints to the control signal matrix. In this scheme, in the calibration phase, the subject does not need to follow a predefined sequence to activate single DOFs, and can even activate more than one DOF simultaneously, which is theoretically impossible in the former classic NMF scheme. Moreover, with this approach, the bases in the synergy matrix can be extracted in one step, contrary to the classic NMF scheme that extracts the bases sequentially for each DOF. For this purpose, we add constraints to the factorization solution and specifically we require the solution to maximize the sparseness of the resulting control function. The sparseness constraint restricts the space of possible NMF solutions. Specifically, the solution with basis functions corresponding to single DOFs, the target solution, is the sparser one among the other infinite solutions. In this way, the factorization with constraints does not require preset calibration phases or activations of single DOFs and can be applied to any set of tasks performed by the user without labeling. In this paper, we present the mathematical background for this new approach and we test it in both able-bodied individuals and amputees. A preliminary version of this work has been reported as a conference abstract in reference [13].

Sparse non-negative matrix factorization
We propose to blindly estimate the control functions through SNMF of multi-channel EMG signals. Before providing an overview of the sparse NMF (SNMF) algorithm and showing how it can be adapted for its use in myoelectric control, we first provide a brief introduction to NMF.
The NMF method consists in finding two non-negative matrices whose product is a good approximation of the recorded signal matrix. At the same time, the NMF representation should be able to learn the parts of the objects (the recorded signal matrices), i.e. it should be able to discover the local signal characteristics in an additive way. This is in contrast to other methods, such as PCA (principal components analysis) and VQ (vector quantization), that learn holistic, not partsbased, representations. For face representation, for example, the basis face in NMF, PCA, and VQ is localized feature [11,12], eigenface [14], and whole face [15], respectively. The non-negative constraint used in NMF allows only additive, not subtractive, combinations, which leads to the parts-based representation [11,12]. The two factorized non-negative matrices are called basis matrix and coefficient matrix (control signal matrix), respectively. NMF has recently received considerable attention as a framework to solve a variety of problems, e.g. document clustering [16], face recognition [17], graph clustering [18,19], DNA gene expression analysis [20], and myoelectric control [9]. However, the NMF [12] algorithm has some limitations, which has been discussed extensively [17,18,21]. One of the notable issues of NMF is that it does not result in sparse enough solutions [21][22][23]. Researchers tried to solve this problem by incorporating sparseness constraints [20,21,[24][25][26][27][28]. These approaches extended the NMF framework to include a sparseness parameter to learn more localized features.
In myoeletric control, a movement of more than one DOF can be decomposed as the linear combination of basic DOFs, and the latent synergy basis and the control signal are nonnegative [29,30].
Let us denote the root mean square (RMS) values of the N-channel T-length surface EMG signal as Z, in which the tth column is the EMG signal at time t and denote m the number of DOFs. Then, the multi-channel observation can be approximated by the product of a N × 2m latent non-negative synergy basis matrix W and a 2m × T latent non-negative control signal matrix F: The above factorization problem can be solved in a semisupervised way. Specifically, by using a DOF-wise training strategy, online experiments showed that NMF achieved good performance [9,29,30]. The DOF-wise training strategy, according to which the user is asked to perform single-DOF movements in a predefined order, is needed to obtain a sparse factorization with NMF. The sparseness of the solution is obtained by imposing a specific single-DOF activation sequence for calibration. Myoelectric control based on NMF is then based on the summation property of DOFs in the muscle signal domain. The control signals are identified for the single DOFs and then used in combination for simultaneous control. If the training is not performed using predefined single-DOF activations, the NMF will not provides a factorization with basis functions that are linear combinations of single DOFs. In this case, the DOFs would not be separated for the control, which would not be intuitive. Unfortunately, single-DOF movements, as required by NMF, are not always possible because the user may not be able to exactly separate the activations of DOFs. In this study we propose a mathematical way to generate a sparse solution without the need for the specific set of calibration data to be generated by single DOF activations. We achieve this by imposing a sparseness constraint to the control signals, which leads to the SNMF [20] scheme. The sparseness constraint that we used determines a much sparser solution than with NMF that corresponds to single DOFs as basis functions.
Mathematically, the degree of sparseness can be controlled by l1-norm or l0-norm. For computational convenience, we selected the l1-norm based SNMF method, SNMF [20]. The objective function of the SNMF method is: where F(:, t) is the tth column vector of F (equation (1)), 'Fro' is the Frobenius norm, and λ > 0 is a regularization parameter to balance the trade-off between the accuracy of the factorization and the degree of sparseness of F. In the experiment, we chose the optimal λ through offline cross-validation (the details are presented in section 3.1). In the current study, m = 2 (see section 2.4 for the experimental setup). The superscript '+' and '−' represent the positive and negative direction of each DOF, respectively. The above formulation can then be rewritten as where e 1×2m is a row vector with all entries equal to 1 and 0 1×T equal to 0. Equation (3) can be efficiently solved via the alternating non-negative least squares (ANLS) method. In ANLS, W and F are iteratively updated by fixing the other one: Equations (4a) and (4b) are the classic least square problem, and each has a closed-form solution. Moreover, the detailed derivation in [20] shows that SNMF can converge to a stationary point. SNMF can extract all the basis functions in one step from recordings generated by arbitrary combinations of DOFs by the user (concurrent way). Only the ordering of the control signals (basis functions) will be undetermined with this approach but this ambiguity can be simply resolved in practice. In practice, we need to record the label of the first two single DOFs to decide the sequence of the basis in the control signal matrix F, which is the reason we refer to our method as quasi-unsupervised instead of completely unsupervised.
In addition to avoid the need for a pre-defined sequence of single-DOF activations, it is expected that the SNMF would outperform the NMF practically since the users (especially patients but even able-bodied subjects) are likely not able to precisely articulate single DOFs individually. If there are coactivations of DOFs in the set of signals used for factorization, the solutions achieved by NMF and SNMF would be different. In this case, the solution of SNMF is expected to be superior because it would impose the sparseness and would be closer to represent single DOFs in the columns of the synergy matrix.

Online estimation of control signals
After obtaining the synergy basis matrix W from a set of EMG signals, we assume it to be constant (at least over the course of the test experiments). To estimate control signals related to the intended activation levels of the DOFs, we apply the pseudo inverse of W (denoted by W + ), and multiply it by newly recorded surface EMG signals Z. The estimated control signals F is therefore: To ensure that no components are overshadowed by the others, each component of F in (6) is normalized with respect to its maximum value. The estimated control signals in (6) are further scaled by the scalar correction factors τ ij , which are used to account for the indetermination of the signal power (range of control signals) in the factorization process (see section 2.4 for details on these factors): where we name F 1 , F 2 the control signals for DOF1 and DOF2, respectively. The multiplicative factors τ ij are determined for each subject so that the final control signals, F = [ F 1 ; F 2 ], match the range of joint angles in the respective DOFs. These factors were manually set for each subject before operating online tasks. The obtained control signals, low-pass filtered at 6 Hz (kinematic bandwidth), are then used by the subjects for controlling the DOFs.

Subjects
Eight able-bodied subjects (six males, two females, 20-45 years old) and two subjects with unilateral limb deficiencies (25 and 46 years old) participated in the study. The experimental protocol was approved by the local ethic committee, and all subjects read and signed the informed consent form prior to their participation.

Experimental protocol
During the experimental session, each subject was instructed to sit comfortably in a neutral position, with the dominant upper limb pointing downwards at the side of the body. We evenly placed 16 monopolar electrodes in two rings around the forearm, at a distance of 1/3 of the forearm length in the distal direction. The distance between the two rings of electrodes (inter-electrode distance in the proximal-distal direction) was 20 mm, and along the arm circumference direction it was in the range 25-35 mm, depending on the size of the forearm. The experimental protocol consisted of two phases. First, a set of surface EMG signals were collected for the factorization (needed for applying NMF that was used as control) and subsequently an online validation test was performed. During the validation phase, the subject performed goal-directed target hitting tasks controlling two DOFs. The detail of the two phases is described in the following.  1). The raw EMG signals collected in this phase were used for training NMF and SNMF, as well as linear regression (LR) that is another method proposed for simultaneous and proportional control of multi-DOF [31,32]. LR estimates the matrix W in equation (5) in the following closed-form (supervised approach): where I is the identity matrix and the regularization constant β is optimized by grid-search in the nested cross-validation [32].
For LR and NMF, we separately trained the data of DOF1 and DOF2 with a fully known trial order. The SNMF was trained either by single-DOF data (without any knowledge of the DOF or trial order) or by combination of DOFs. After the algorithms were calibrated, they were tested online. The training procedure for NMF consists in solving equation (1); Z is factorized to W and F, W is the basis matrix, and F is the control signal matrix. The testing procedure for NMF consists in calculating equation (5); F is the control signal that we obtain. The training procedure for SNMF is similar to that for NMF with the difference that the sparseness constraint is added to F when factorizing Z. The training procedure for LR consists in calculating the LR model by solving the liner system of equations.

Online control phase.
The subject controlled the same arrow as in the training phase using the control signals obtained as in equation (5) for NMF and SNMF or by LR. The subject controlled the same arrow as in the training phase using the control signals obtained as in equation (5) for NMF and SNMF or by LR ( figure 2). The online control was in position mode, i.e. without surface EMG activity, the arrow returned to its original position. Being in position mode the control needed to be proportional since the position depended on the intensity of activation, i.e. the spatial displacement of the arrow on the respective DOF is directly related to the control signals in equation (7). Considering that the relative magnitudes of the 2 directions of the two DOFs were undetermined and the range of motion did not necessarily match the intended range, we scaled the factors to all four bases so that the subjects could reach the maximal displacements in all directions. The scaling was subject-specific and was made in a simple practical way. Each subject was asked to activate each DoF maximally and the maximum activation was assigned to the maximum range of motion of the corresponding DoF. This was done before testing and was not adapted during testing. The online control performance of each algorithm was evaluated through a series of target hitting tests. Targets were displayed sequentially at different locations on the screen in the form of empty circles. The radius of each circle was 8 density-independent pixels (dp) with the entire working interface of 360 × 60 dp. The subjects controlled the arrow using their forearm muscle contractions to place its tip within the circle for at least 300 ms, within a 20 s time interval (figure 1). If this was accomplished, the trial was considered successful. Otherwise, the trial was considered as a failure. Three types of targets were presented to the subjects, with 20 targets for each type. The targets of the first type were presented randomly on the unit-circle with different rotation angles. In this case, the subjects could complete the task by activating DOF2 only, and thus rotating the arrow on the screen. The targets of the second type were randomly displaced on the left and right side. Therefore, in this case, the tasks could be accomplished by moving the arrow to the left or right by activating DOF1 only. The targets of the third type were placed on the entire working interface so that the subjects had to simultaneously control both DOF1 and DOF2 (position control mode). Figure 1 illustrates these three types of targets. To fairly assess the performance of the methods, the order of the targets was presented randomly for each subject. In addition, the order of the tested algorithms was also randomly selected and the subject was blind to the algorithm used. The comparison between LR, NMF and the proposed SNMF was done for able-bodied subjects only whereas only SNMF and NMF were compared in amputees. This choice was due to time constraints in the experiments with amputees.

Statistical analysis
In all ANOVA tests (2-sided tests), the full model was conducted. The main factors of the analysis were the target type (single DOF activations and two DOFs activation), and the algorithm (LR, NMF and SNMF). The subject was a random factor. First, the significant interaction between the two factors (target type and algorithm) was tested. When a significant interaction was detected, focused one-way ANOVA were performed by fixing the level of one of the interacting factors. In case of no significant interactions, the effects of the interactions were pooled and only the effects of the main factors were analyzed. When statistical significance was detected, the Tukey-Kramer multiple comparison with Bonferroni correction was performed to identify the order of the levels of that factor. The significance level was set to 0.05 for all tests.

Evaluation of the sparseness of the control signals
We used an average signal-to-noise ratio (ASNR) to characterize the factorization solutions obtained with NMF and SNMF. During trials with activation of single DOFs in the test set, the single trial SNR was the ratio between the intended control signal (integrated over time) and the unintended control signal: where i indicates the trial index, [T i + 1, T i+1 ] is the time interval of the ith trial, and the subscript 'int' and 'unint' of F(t) respectively represent the intended and unintended control signal at time t. When a single DOF activation is required, the intended control signal is the one corresponding to the DOF to be activated and the unintended control signal is the signal activation the other DOF. Then, the ASNR over I trials was obtained by averaging the single trial SNR: Throughput (TP, (bit s −1 )) Ratio between the task difficulty index of each target and the completion time Completion time (CT, (s)) The time it took the subject to complete the successful attempt Efficiency coefficient (EC, (%)) Ratio between the length of the optimal path from the initial point to the target and the actual trajectory realized (100% indicates a perfect execution) Overshoots (O) Number of occurrences that the tip of the arrow passes through the target before the dwelling time was reached

Significance of sparseness for control signals
First we illustrate the effect of sparseness constraint on the control signals when using NMF and SNMF. For this purpose we show the control signals obtained by the two methods when using the same set of signals for the factorization. Specifically, only signals generated by single-DOF activations of one of the subjects were used for factorization, according to the needs for NMF. In all the experiments, there were 9 trials of DOF1 (including Extension and Flexion), and 9 trials of DOF2 (including Supination and Pronation). As mentioned previously, for SNMF, we need to choose the optimal param eter λ through cross-validation. The detailed selection procedure of λ opt is as follows. We first collected offline trials (for each subject) and randomly split them into K-folds (K = 3 in this paper). Then, we took each 2 of the 3 folds (2/3 trials) to form the training set and the remaining fold (1/3 trials) to form the validation set. By rotating the validation, we obtained a total of 3 training sets and the respective 3 validation sets for each subject. After that, we predefined a series of values of λ (e.g. from 0.001 to 10). The range of λ can be arbitrary. λ determines the degree of sparseness (a smaller λ generates less sparseness). For each λ, we used each of the training set to learn the basis function via equation (2), and applied them to the respective validation dataset to obtain the control signal and the ASNR. Thus, we had 3 ASNRs for each λ and average them as the final ASNR. On evaluating all λs, we selected the optimal λ as that associated with the minimal ASNR. Finally, we used the learnt optimal λ (among the predefined scope) for online testing. Figure 3 shows the differential control signals (as described in equation (7)) of one subject generated by NMF and SNMF during testing with the user attempting to control individual DOFs. For this example, NMF and SNMF were based on factorizing the EMG generated during single DOF activations, as requested by NMF. When using SNMF, the control signal corresponding to the unintended DOF was substantially more suppressed than when using NMF. With NMF, in some cases, the control signal of the unintended DOF had similar magnitude as that of the intended DOF. According to the representative results in figure 3, the resulting ASNR for NMF was 2.30 while for SNMF it was 7.16. Therefore, even when single-DOF activations were used for deriving the matrix W, the sparseness constraint led to a superior solution with respect to NMF by reducing the undesirable activations of the unintended DOF.

Online control performance
The NMF, SNMF and LR were compared in the online performance using the same set of signals for the extraction of the matrix W. Since NMF requires single-DOF activations, these signals were first obtained from the single-DOF trials (ideal conditions for NMF). The results in this condition will be shown first and will be followed by those obtained when using a combination of DOF as input for the factorization.
For single-DOF trials in the factorization, the statistical results of ANOVA for LR, NMF, and SNMF for able-bodied subjects are summarized in figure 4. In all ANOVA tests, the full model was tested first. The full model ANOVA for the 4 performance indices found no significance interaction between the two main factors (target type and algorithm) (p = 0.464, p = 0.153, p = 0.207, and p = 0.082 for TP, CT, EC, and O, respectively). Therefore, the reduced ANOVA was tested only for the main factors (target type and algorithm), and Bonferroni comparison was used to show the statistical significance. The TP was significantly greater for SNMF (3.01 ± 1.99 bit s −1 ) than NMF (2.58 ± 1.51 bit s −1 ; p < 0.001) and LR (2.47 ± 1.61 bit s −1 ; p < 0.001), whereas it was similar for NMF and LR (p = 0.794). The task completion   The statistical results for the amputee patients are shown in figure 5. The TP was 1.92 ± 1.99 bit s −1 (NMF) and 2.53 ± 1.83 bit s −1 (SNMF); task completion time (CT) was 7.44 ± 6.28 s (NMF) and 5.33 ± 5.58 s (SNMF); efficiency coefficient (EC) was 33.67 ± 27.37% (NMF) and 47.62 ± 28.27% (SNMF); and overshoot (O) was 2.27 ± 4.31 (NMF), and 0.98 ± 1.56 (SNMF). There was no interaction between the target type and algorithm (p = 0.941, p = 0.665, p = 0.062, and p = 0.172 for TP, CT, EC, and O, respectively). The online performance of SNMF was significantly better than that of NMF according to all 4 performance indices (p = 0.013 for TP, p = 0.004 for CT, p < 0.001 for EC, and p = 0.006 for O).
Finally, we tested the robustness of the solution provided by SNMF with mixed multi-DOF signals as input. For this test, we again compared SNMF and NMF. When NMF was trained (for the factorization) with multi-DOF data, the subjects were not able to control the interface and could not complete any task. Conversely, the online control based on SNMF was possible. The results obtained by SNMF when applied to multi-DOF data for the factorization were compared with the results obtained by NMF when the factorization was performed using single-DOF data (as in the previous tests), which was the only condition for which control with NMF was possible. This comparison was conducted on two able-bodied subjects from the initial subject sample and its results are shown in figure 6. There was no significant difference between the control by SNMF and NMF even though the first was trained by multi-DOF data (arbitrary movements of the user) and the second by single-DOF constrained activations. Therefore, benefiting from the sparseness constraint in the factorization stage, SNMF allowed simultaneous and proportional myoelectric control independent on the data used for extracting W.

Discussion
The proposed method provides a solution for simultaneous and proportional myoelectric control given EMG signals without labeling or any other constraint. Therefore, the method is  basically unsupervised (only the ordering of the basis of the control signals is undetermined, but this ambiguity can be simply resolved in practice, as discussed in Methods). It is based on the assumption that in the muscle signal domain, activations can be obtained as linear combinations of basis functions. The best basis functions are those corresponding to single-DOF since these are naturally associated to the physiological muscle activation determining the corresponding tasks. However, extracting a set of basis functions from EMG factorization poses indeterminacy when no further conditions are imposed. Indeed, the same space of activations in the muscle signal domain can be covered by any linear combination of the activation signals corresponding to the single-DOFs. For these reasons, NMF applied in myoelectric control requires that the basis functions are extracted one by one from signals generated during single-DOF tasks of the user (DOF-wise training strategy). The proposed SNMF approach eliminates this need since, from a set of EMG signals generated by an arbitrary task, it selects the factorization with maximum sparseness, which is the one corresponding to basis functions associated to the single-DOFs.
In addition to the advantage of a better set of basis functions even when the signals used are from single-DOF activations, the SNMF method can be used for factorizing signals without the constraint of being generated from single-DOF data (since the activations cannot be in practice ideal single-DOF anyway). Therefore, the same solution will be achieved for a large repertoire of signals, in an unsupervised way.
Sparse representation algorithms have several applications in signal processing. Sparse representation classification (SRC) methods can resist the presence of noisy data and preserve high classification accuracy for corrupted data [33], as also shown for EEG and EMG processing [34,35].

Conclusion
We proposed a SNMF scheme for robust and concurrent extraction of basis functions for myoelectric control in upper limb prostheses with an unsupervised scheme, which not only simplifies the training stage for the users, but also improves the quality of the control. This factorization scheme has advantages over the classic NMF in that it selects the sparse solution among infinite possible solutions. In this way, it does not need a preset calibration phase of constraint tasks by the user (the DOF-wise based training scheme) but rather the constraint on the solution is mathematically imposed. This new scheme has been compared to NMF and LR and proved to be superior to these methods when using single-DOF activations for the factorization (ideal condition for NMF). In addition, the method can be used by factorizing EMG signals from arbitrary tasks made of combinations of DOFs.
In conclusion, we have proposed and validated a factorization algorithm for myoelectric control that has superior performance than previously tested methods for the simultaneous and proportional control of multiple DOFs and may pave the way to adaptive myoelectric control systems.