Abstract

Nonnegative matrix factorization (NMF) by the multiplicative updates algorithm is a powerful machine learning method for decomposing a high-dimensional nonnegative matrix V into two nonnegative matrices, W and H, where . It has been successfully applied in the analysis and interpretation of large-scale data arising in neuroscience, computational biology, and natural language processing, among other areas. A distinctive feature of NMF is its nonnegativity constraints that allow only additive linear combinations of the data, thus enabling it to learn parts that have distinct physical representations in reality. In this letter, we describe an information-theoretic approach to NMF for signal-dependent noise based on the generalized inverse gaussian model. Specifically, we propose three novel algorithms in this setting, each based on multiplicative updates, and prove monotonicity of updates using the EM algorithm. In addition, we develop algorithm-specific measures to evaluate their goodness of fit on data. Our methods are demonstrated using experimental data from electromyography studies, as well as simulated data in the extraction of muscle synergies, and compared with existing algorithms for signal-dependent noise.

1.  Introduction

Nonnegative matrix factorization (NMF) was introduced as an unsupervised, parts-based learning paradigm, in which a high-dimensional nonnegative matrix V is decomposed into two matrices, W and H, each with nonnegative entries, , by a multiplicative updates algorithm (Lee & Seung, 1999, 2001). In the past decade, NMF has been increasingly applied in a variety of areas involving large-scale data. These include neuroscience, computational biology, natural language processing, information retrieval, biomedical signal processing and image analysis. (For a review of its applications, see Devarajan, 2008.)

Lee and Seung (2001) outlined algorithms for NMF based on the gaussian and Poisson likelihoods. Since their seminal work, numerous variants, extensions, and generalizations of the original NMF algorithm have been proposed in the literature. For example, Hoyer (2004), Shahnaz, Berry, Pauca, and Plemmons (2006), Pascual-Montano, Carazo, Kochi, Lehmann, and Pascual-Marqui (2006), and Berry, Browne, Langville, Pauca, and Plemmons (2007) extended NMF to include sparseness constraints. Wang, Kossenkov, and Ochs (2006) introduced LS-NMF that incorporated variability in the data. Cheung and Tresch (2005) and Devarajan and Cheung (2012) extended the NMF algorithm to include members of the exponential family of distributions while Devarajan and Ebrahimi (2005, 2008), Devarajan (2006), and Devarajan, Wang, and Ebrahimi (2011) formulated a generalized approach to NMF based on the Poisson likelihood that included various well-known distance measures as special cases. Dhillon and Sra (2005) and Kompass (2007) have also proposed generalized divergence measures for NMF. Cichocki, Zdunek, and Amari (2006), Cichocki, Lee, Kim, and Choi (2008), Cichocki, Zdunek, Phan, and Amari (2009), and Cichocki, Cruces, and Amari (2011) extensively developed a series of generalized algorithms for NMF based on - and -divergences while Févotte and Idier (2011) recently extended it by proposing some novel algorithms. The work of Cichocki et al. (2009) provides a detailed reference on this subject.

The main focus of this letter is on NMF algorithms for signal-dependent noise with particular emphasis on the generalized inverse gaussian family of distributions. This family includes the well-known gamma model for signal-dependent noise as a special case. It also includes the inverse gaussian model as a special case, among others. Each model incorporates signal dependence in noise in structurally different ways based on the mean-variance relationship, as evidenced in the following sections. These models are embedded within the framework of the exponential family of models outlined in Cheung and Tresch (2005) and can be obtained as special cases of -divergence proposed in Cichocki et al. (2006, 2009). In each case, the NMF algorithm is based on maximizing the likelihood or, equivalently, minimizing a cost function defined by the Kullback-Leibler divergence between the input matrix V and the reconstructed matrix WH.

We describe an approach to NMF for signal-dependent noise by extending the standard likelihood approach to include two well-known alternative cost functions from information theory for quantifying this divergence: the dual Kullback-Leibler divergence and the J-divergence. Based on these measures, we propose three NMF algorithms applicable when the data exhibit signal-dependent noise. For each algorithm, we provide a rigorous proof of monotonicity of updates using the EM algorithm. We describe a principled method for selecting the appropriate rank of the factorization and develop algorithm-specific measures to quantify the variation explained by the chosen model. We demonstrate the applicability of our methods using experimental data from electromyography (EMG) studies, as well as simulated data in extracting muscle synergies, and compare the performance of our proposed methods with existing algorithms for signal-dependent noise.

The remainder of the letter is organized as follows. Section 2 provides the necessary background required for the information-theoretic approach we describe. It is intended to serve as a brief tutorial on fundamental concepts, terminology, basic quantities of interest for our problem, and their interpretations. Section 3 provides a detailed overview of existing NMF algorithms for signal-dependent noise and places them within the broader context of NMF algorithms based on generalized divergence measures. Furthermore, it describes our proposed NMF algorithms for signal-dependent noise and provides multiplicative update rules. Section 4 outlines methods for model selection and evaluation, while section 5 presents an application of our methods to EMG data and a comparison to existing methods. Section 6 provides a summary and conclusions. Detailed proofs of monotonicity of updates for the proposed algorithms are relegated to the appendix.

2.  Background and Preliminaries

2.1.  Directed Divergence and Divergence.

Suppose we are interested in testing a set of hypotheses denoted by Hi, i=0, 1 that a random variable X is from population i with probability measure Assume that and are absolutely continuous with respect to each other and that X takes values on the entire real line. Let P(Hi) denote the prior probabilities, f, g the density functions, and f, g the distribution functions corresponding to the hypothesis Hi, i=0, 1, respectively. If P(Hi|x) denotes the conditional (or posterior) probability of Hi given X=x, then, using Bayes’ theorem, we have
formula
and
formula
Hence,
formula
that is, the logarithm of the likelihood ratio, defined as the negative difference between the logarithm of the odds in favor of H0 before and after the observation X=x, is the information in X=x for discrimination in favor of H0 against H1 (Kullback, 1959).
Suppose that x is not given and there is not specific information on the whereabouts of x other than . The mean information per observation, averaged over all the values x of X, for discrimination in favor of H0 against H1 is thus
formula
2.1
This quantity is known as Kullback-Leibler divergence between f and g, the negative log likelihood or empirical entropy. Similarly, the measure I(g, f) is defined as the mean information per observation, averaged over all the values x of X, for discrimination in favor of H1 against H0 and is given by
formula
2.2
This quantity is known as the dual Kullback-Leibler divergence between f and g or the negative empirical log likelihood. In light of these definitions, I(f, g) and I(g, f) are also referred to as directed divergences. These quantities are nonnegative definite and are zero if and only if f(x)=g(x) almost everywhere (Kullback, 1959; Owen, 2001).
Using directed divergences I(f, g) and I(g, f), one can define J-divergence J(f, g) as
formula
2.3
which is a measure of the divergence or the difficulty of discriminating between the hypotheses H0 and H1. A key feature of J(f, g) is symmetry with respect to the measures and . It has all the properties of a distance measure (metric) except the triangle inequality, is nonnegative definite, and is zero if and only if f(x)=g(x) almost everywhere (Kullback, 1959).

2.2.  Motivating NMF for Signal-Dependent Noise.

2.2.1.  The Generalized Inverse Gaussian Distribution.

A nonnegative random variable X is said to be a member of the family of generalized inverse gaussian (GIG) distributions if its probability density function is given by
formula
2.4
where , , , and, is the modified Bessel function of the third kind with index . In the limiting case when , f(x) reduces to a gamma distribution with density
formula
2.5
where and . The mean-variance relationship can be written as
formula
2.6
and thus indicates a quadratic dependence of variance on the mean. When , f(x) reduces to an inverse gaussian distribution with density
formula
2.7
where and . The mean-variance relationship can be written as
formula
2.8
and thus indicates a cubic dependence of variance on the mean.

Other special cases of the GIG family of distributions include the inverse gamma and hyperbolic distributions (Eberlein & Hammerstein, 2004). In this letter, we focus on the gamma and inverse gaussian distributions as data-generating models for signal-dependent noise in the context of NMF.

2.2.2.  Divergence Measures for Signal-Dependent Noise.

The discrimination information functions defined above were introduced by Kullback and Leibler (1951) and serve as divergence measures for comparing two distributions or probability models. Using appropriate densities for f(x) and g(x) in equations 2.1 to 2.3 based on the gaussian (normal), gamma or inverse gaussian models, one can obtain various divergence measures for NMF based on the empirical entropy, empirical likelihood, or a combination of these. Throughout the presentation, we use KL, KLd, and J to denote Kullback-Leibler, dual Kullback-Leibler, and J-divergence, respectively. In each case, the subscripts N, G, and IG are used to refer to the gaussian, gamma, and inverse gaussian models, respectively. The term KL divergence has been used in the literature to refer to that based on the Poisson model. However, it is important to note that the terms KL, dual KL, and J-divergence used in this letter refer to generic divergence measures between any two densities f and g as defined in equations 2.1 to 2.3 and are not model specific. Each divergence measure can be defined specifically for a particular choice of model, as outlined below.

For two normal random variables with means and (and equal variance ) and corresponding probability density functions f(x) and g(x), it can be shown that
formula
2.9
Using equation 2.5, it can be shown that for two gamma random variables with means and (and common shape parameter ) and corresponding probability density functions f(x) and g(x),
formula
2.10
formula
2.11
and
formula
2.12
Similarly, using equation 2.7, it can be shown that for two inverse gaussian random variables with means and (and common shape parameter ) and corresponding probability density functions f(x) and g(x),
formula
2.13
and
formula
2.14
It will become clear in section 3 that none of the parameter coefficients in the divergence equations 2.9 to 2.14 play any role in the derivation of the NMF algorithms. Hence, we assume that without loss of generality. It should be noted that this assumption is consistent with the basic formulation in NMF and that numerous other divergence measures available in the literature for NMF implicitly make such assumptions (Cichocki et al., 2006, 2008, 2009; Févotte & Idier, 2011; Kompass, 2007; Devarajan & Cheung, 2012).

We motivate nonnegative matrix factorizations in the context of dimension reduction of high-dimensional electromyography (EMG) data. EMG data are typically presented as a matrix in which the rows correspond to different muscles, the columns to disjoint, sequentially sampled time intervals, and each entry to the EMG signal of a given muscle in a given time interval. In EMG studies, the number of muscles, p, is typically, fewer than fifty; the number of time intervals, n, is typically in the tens of thousands; and the matrix of EMG signal intensities V is of size so that each column of V represents an activation vector in the muscle space at one time instance. Our goal is to find a small number of muscle synergies, each defined as a nonnegative, time-invariant activation balance profile in the p-dimensional muscle space. This is accomplished with a decomposition of the matrix V into two matrices with nonnegative entries, , where W has a size , so that each column is a time-invariant muscle synergy in the p-dimensional muscle space and the matrix H has size , so that each column contains the activation coefficients for the r synergies in W for one time instance. The number of synergies r is chosen so that (n+p)r<np. The entry hai of H is the coefficient of time interval i in synergy a, and the entry wja of W is the expression level of synergy a in muscle j, where .

The first step in obtaining an approximate factorization for V is to define cost functions that measure the divergence between the observed matrix V and the product of the factored matrices WH. We can express this in the form of a linear model as
formula
2.15
where represents noise. NMF algorithms for signal-dependent noise based on KL divergence (see equations 2.10 and 2.13, respectively) for gamma and inverse gaussian models exist in the literature (Cheung & Tresch, 2005; Cichocki et al, 2009; Févotte & Idier, 2011). In this letter, we propose three novel NMF algorithms for handling signal-dependent noise, specifically based on dual KL and J-divergence for gamma and inverse gaussian models (see equations 2.11, 2.12, and 2.14). The appropriate cost function for each model is obtained by simply substituting and in equations 2.9 to 2.14 with V and WH, respectively.

2.2.3.  Application of NMF to EMG Data.

We present an application involving the analysis of EMG data, electrical signals recorded from muscles that reflect how they are activated by the nervous system for a particular posture or movement. It is well known in the literature that EMG data exhibit signal-dependent noise (Harris & Wolpert, 1998; Cheung, d'Avella, Tresch, & Bizzi, 2005). One long-standing question in neuroscience concerns how the motor system coordinates the activations of hundreds of skeletal muscles, representing hundreds of degrees of freedom to be controlled (Bernstein, 1967). They define an immense volume of possible motor commands that the central nervous system (CNS) must search through for the execution of even an apparently simple movement. It has been argued that the CNS simplifies the complexity of movement and postural control arising from high dimensionality by activating groups of muscles as individual units, known as muscle synergies (Tresch, Saltiel, d'Avella, & Bizzi, 2002; Giszter, Patil, & Hart, 2007; Bizzi, Cheung, d'Avella, Saltiel, & Tresch, 2008; Ting, Chvatal, Safavynia, & McKay, 2012). As modules of motor control, muscle synergies serve to reduce the search space of motor commands, reduce potential redundancy of motor commands for a given movement, and facilitate learning of new motor skills (Poggio & Bizzi, 2004). Numerous laboratories have utilized linear factorization algorithms to extract muscle synergies from multichannel EMGs recorded from humans and animals (reviewed in Bizzi & Cheung, 2013). In particular, several studies have demonstrated that the muscle synergies returned by the gaussian NMF could be neurophysiological entities utilized by the CNS for the production of natural motor behaviors (Saltiel, Wyler-Duda, d'Avella, Tresch, & Bizzi, 2001; Tresch, Cheung, & d'Avella, 2006; Overduin, d'Avella, Carmena, & Bizzi, 2012).

It has been further posited that muscle synergies for locomotion are basic units of the so-called central pattern generators (CPGs) whose organizations are independent of the pattern of sensory feedback. Cheung et al. (2005) sought to demonstrate this possibility by recording hind-limb EMGs from bullfrogs during jumping and swimming, before and after deafferentation or the surgical procedure of eliminating sensory inflow into the spinal cord by cutting the dorsal nerve roots. By applying a manipulated version of the gaussian NMF to the data matrix that pooled the intact and deafferented EMGs together, they found three to six, out of four to six, muscle synergies were preserved after deafferentation. The preserved synergies were then interpreted as basic components of the CPGs. Here, we ask whether the NMF algorithms based on gamma or inverse gaussian noise can better identify CPG components by discovering more muscle synergies shared between the pre- and postdeafferentation data sets than the traditional gaussian NMF. We hypothesize that the NMF algorithms derived from signal-dependent noise outperform the gaussian NMF in their ability to discover shared muscle synergies, because signal-dependent noise formulations should better model the noise properties of EMG signals than a gaussian formulation (Harris & Wolpert, 1998).

The data analyzed here were previously described in Cheung et al. (2005). EMGs during unrestrained jumping and swimming were collected from four adult bullfrogs (Rana catesbeiana) before and after a complete hind limb deafferentation was achieved by severing dorsal roots 7 to 9. Intramuscular EMG electrodes were surgically implanted into the following muscles in the right hind limb: rectus internus major (RI), adductor magnus (AD), semimembranosus (SM), semitendinosus (ST), iliopsoas (IP), vastus internus (VI), rectus femoris anticus (RA), gastrocnemius (GA), tibialis anticus (TA), peroneus (PE), biceps (BI), sartorius (SA), and vastus externus (VE). The collected EMG signals were amplified (gain of 10,000) and bandpass-filtered (10–1000 Hz) through differential alternating-current amplifiers, then digitized at 1000 Hz. Using custom software written in Matlab (R2010b; Math-Works, Natick, MA), the EMG signals were further high-pass-filtered with a window-based finite impulse response (FIR) filter (50th order; cutoff of 50 Hz) to remove any motion artifacts, then rectified, low-pass-filtered (FIR; 50th order; 20 Hz), and finally integrated over 10 ms intervals. The preprocessed data of each muscle were then normalized to the maximum EMG value of that muscle attained in the entire experiment.

The EMG signal is a spatiotemporal summation of the motor action potentials traveling along the muscle fibers of the thousands of motor units in the recorded muscle. The high-frequency components of the EMG reflect, in addition to noise, the contribution of these action potentials. In motor neuroscience, it is customary to perform low-pass filtering on the recorded EMGs to obtain an “envelope” of muscle activation, which should reflect the higher-level control signals originating from the brain and spinal cord that specify the degree of muscle contraction for generating the desired force (with the force magnitude dictated by the muscle's force-length and force-velocity relationships). Since we are interested in discovering structures at the level of control signals for muscular contraction (i.e., muscle synergies) and since signal dependent noise is thought to occur at this control-signal level, it is appropriate to apply NMF to filtered EMG data. There is a sizable literature on using the gaussian NMF algorithm for extracting muscle synergies from filtered EMG data (reviewed in Bizzi & Cheung, 2013). Moreover, filtering the EMGs before NMF extraction would allow an easier comparison of our results with those in the literature.

One approach to visualizing signal dependence in these data is to plot the variability as a function of the mean for EMG signals from each muscle separately. Any observed trend in the mean-variance relationship would indicate some form of signal dependence, the exact nature of the relationship being dependent on the data-generating mechanism. For each muscle, the mean and variance of the EMG signal were computed for moving windows across time. Several window sizes ranging from 3 to 50 were explored, and it was observed that mean-variance relationship was not sensitive to the choice of window size. Our choice of window size was based on a physiological justification. There has been an earlier result suggesting that bursts with 275 msec duration could be a fundamental pulse unit in the frog spinal cord (Hart & Giszter, 2004). Since our integration time interval is 10 ms, we used a window size of 28 so that each window corresponded to the duration of this fundamental drive.

Using equation 2.6, we can rewrite the mean-variance relationship for a gamma model in terms of the standard deviation as
formula
Taking the logarithm on both sides, we obtain
formula

Figure 1A shows a plot of the logarithm of the estimated standard deviation against the logarithm of the estimated mean for moving windows across time for the intact jump of one frog for the muscle TA. The black solid line represents a linear fit to the data. The proximity of the estimated slope of this line to unity provides strong evidence that a gamma model adequately represents the mean-variance relationship in these data. The logarithmic transformation provides variance stabilization and aids in interpreting the slope of the fit. Panels A to D in Figure 1 display such plots for each of the four behaviors of selected muscles and frogs. At the top of each panel in this figure, the estimate of the slope and goodness-of-fit measures such as the root mean squared error (RSE) and adjusted R2 are listed along with frog behavior and name of muscle. This figure is representative of the mean-variance relationship typically observed in our frog EMG data. Table 1 lists the estimates of slope and adjusted R2 (mean SD for N=4 frogs) from the least squares fit for each behavior and muscle. It is evident from these results that the gamma model provides an overall good fit of the EMG data.

Figure 1:

Illustration of the mean-variance relationship for the frog EMG data. Plot of the logarithm of the estimated standard deviation against the logarithm of the estimated mean for moving windows across time for each behavior of selected muscles and frogs. Each panel displays the mean-variance relationship for a particular behavior. (A) Intact jump. (B) deafferented jump. (C) intact swim. (D) deafferented swim. In each panel, the black solid line represents a linear fit to the data and estimates of the slope, root mean squared error (RSE), and adjusted R2 are listed at the top of each panel.

Figure 1:

Illustration of the mean-variance relationship for the frog EMG data. Plot of the logarithm of the estimated standard deviation against the logarithm of the estimated mean for moving windows across time for each behavior of selected muscles and frogs. Each panel displays the mean-variance relationship for a particular behavior. (A) Intact jump. (B) deafferented jump. (C) intact swim. (D) deafferented swim. In each panel, the black solid line represents a linear fit to the data and estimates of the slope, root mean squared error (RSE), and adjusted R2 are listed at the top of each panel.

Table 1:
Empirical Assessment of Mean-Variance Relationship: Estimated Slope () and Adjusted R2 from Least Squares Fit of Standard Deviation versus Mean (Log Scale) (Mean SD for N=4 Frogs).
Measure of FitBehaviorIPVIRAGARI
 Intact jump 1.03 0.09 1.25 0.09 1.21 0.09 1.21 0.05 1.82 0.11 
 Deafferented jump 1.02 0.04 1.18 0.10 1.13 0.07 1.23 0.04 1.52 0.09 
 Intact swim 1.09 0.06 1.32 0.12 1.13 0.06 1.21 0.14 2.37 0.91 
 Deafferented swim 1.12 0.12 1.30 0.15 1.11 0.14 1.21 0.08 2.38 0.63 
Adjusted R2 Intact jump 0.92 0.04 0.94 0.05 0.92 0.02 0.94 0.02 0.85 0.02 
 Deafferented jump 0.88 0.02 0.92 0.03 0.91 0.02 0.94 0.01 0.90 0.15 
 Intact swim 0.82 0.03 0.90 0.02 0.87 0.05 0.90 0.06 0.74 0.11 
 Deafferented swim 0.84 0.03 0.88 0.04 0.85 0.04 0.90 0.03 0.78 0.06 
  AD SM ST TA PE 
 Intact jump 1.35 0.08 1.32 0.08 0.90 0.05 1.04 0.04 1.16 0.09 
 Deafferented jump 1.31 0.12 1.27 0.04 0.90 0.05 0.94 0.07 1.08 0.12 
 Intact swim 1.52 0.11 1.50 0.10 1.16 0.10 1.04 0.05 1.02 0.09 
 Deafferented swim 1.48 0.29 1.56 0.22 1.33 0.16 1.09 0.13 1.03 0.13 
Adjusted R2 Intact jump 0.93 0.01 0.91 0.03 0.82 0.08 0.90 0.03 0.92 0.04 
 Deafferented jump 0.92 0.01 0.90 0.02 0.76 0.08 0.86 0.03 0.89 0.05 
 Intact swim 0.89 0.02 0.90 0.01 0.86 0.01 0.80 0.05 0.82 0.05 
 Deafferented swim 0.88 0.03 0.90 0.03 0.88 0.03 0.88 0.09 0.81 0.06 
  BI SA VE   
 Intact jump 1.03 0.03 1.17 0.08 1.34 0.12   
 Deafferented jump 1.03 0.06 1.06 0.03 1.17 0.04   
 Intact swim 1.15 0.05 1.24 0.18 1.31 0.17   
 Deafferented swim 1.19 0.10 1.34 0.15 1.21 0.18   
Adjusted R2 Intact jump 0.88 0.04 0.88 0.04 0.92 0.02   
 Deafferented jump 0.89 0.02 0.88 0.03 0.90 0.03   
 Intact swim 0.86 0.03 0.86 0.03 0.88 0.04   
 Deafferented swim 0.88 0.04 0.86 0.06 0.86 0.04   
Measure of FitBehaviorIPVIRAGARI
 Intact jump 1.03 0.09 1.25 0.09 1.21 0.09 1.21 0.05 1.82 0.11 
 Deafferented jump 1.02 0.04 1.18 0.10 1.13 0.07 1.23 0.04 1.52 0.09 
 Intact swim 1.09 0.06 1.32 0.12 1.13 0.06 1.21 0.14 2.37 0.91 
 Deafferented swim 1.12 0.12 1.30 0.15 1.11 0.14 1.21 0.08 2.38 0.63 
Adjusted R2 Intact jump 0.92 0.04 0.94 0.05 0.92 0.02 0.94 0.02 0.85 0.02 
 Deafferented jump 0.88 0.02 0.92 0.03 0.91 0.02 0.94 0.01 0.90 0.15 
 Intact swim 0.82 0.03 0.90 0.02 0.87 0.05 0.90 0.06 0.74 0.11 
 Deafferented swim 0.84 0.03 0.88 0.04 0.85 0.04 0.90 0.03 0.78 0.06 
  AD SM ST TA PE 
 Intact jump 1.35 0.08 1.32 0.08 0.90 0.05 1.04 0.04 1.16 0.09 
 Deafferented jump 1.31 0.12 1.27 0.04 0.90 0.05 0.94 0.07 1.08 0.12 
 Intact swim 1.52 0.11 1.50 0.10 1.16 0.10 1.04 0.05 1.02 0.09 
 Deafferented swim 1.48 0.29 1.56 0.22 1.33 0.16 1.09 0.13 1.03 0.13 
Adjusted R2 Intact jump 0.93 0.01 0.91 0.03 0.82 0.08 0.90 0.03 0.92 0.04 
 Deafferented jump 0.92 0.01 0.90 0.02 0.76 0.08 0.86 0.03 0.89 0.05 
 Intact swim 0.89 0.02 0.90 0.01 0.86 0.01 0.80 0.05 0.82 0.05 
 Deafferented swim 0.88 0.03 0.90 0.03 0.88 0.03 0.88 0.09 0.81 0.06 
  BI SA VE   
 Intact jump 1.03 0.03 1.17 0.08 1.34 0.12   
 Deafferented jump 1.03 0.06 1.06 0.03 1.17 0.04   
 Intact swim 1.15 0.05 1.24 0.18 1.31 0.17   
 Deafferented swim 1.19 0.10 1.34 0.15 1.21 0.18   
Adjusted R2 Intact jump 0.88 0.04 0.88 0.04 0.92 0.02   
 Deafferented jump 0.89 0.02 0.88 0.03 0.90 0.03   
 Intact swim 0.86 0.03 0.86 0.03 0.88 0.04   
 Deafferented swim 0.88 0.04 0.86 0.06 0.86 0.04   

In the following section, we discuss NMF algorithms for signal-dependent noise. We begin with a survey of existing work in this area before describing three novel algorithms for this problem. A detailed analysis of the data sets described here, including a comparison of existing approaches to our proposed methods, is presented in section 5.

3.  NMF Algorithms for Signal-Dependent Noise

3.1.  Existing Work.

Cheung and Tresch (2005) proposed a heuristic NMF algorithm for the exponential family of distributions that embeds the gaussian, Poisson, gamma, and inverse gaussian models. They provided generalized multiplicative update rules for W and H by modifying the step size in the gradient based on the negative log likelihood (or, equivalently, KL divergence). In independent work, Cichocki et al. (2006) also proposed a similar heuristic algorithm based on the generalized -divergence and provided multiplicative update rules for W and H. -divergence between the input matrix V and reconstructed matrix WH is given by
formula
3.1
-divergence includes the gaussian , Poisson , gamma , and inverse gaussian models as special cases. It should be noted that this generalized divergence is related to the other divergence measures independently described in the literature (such as those in Kompass, 2007; Cichocki et al., 2008; Devarajan et al., 2011; and Devarajan & Ebrahimi, 2005) via transformations. For NMF algorithms based on gamma and inverse gaussian models stemming from the work of Cheung and Tresch (2005) and Cichocki et al. (2006), monotonicity of updates cannot be established, and they remain heuristic. Recently, however, Févotte and Idier (2011) proposed a rigorous majorization-maximization (MM) algorithm based on -divergence that enables monotonicity of updates for W and H to be theoretically established. Moreover, they provided generalized multiplicative update rules for W and H that were seen to be different from heuristic updates. We refer to these as the heuristic and MM algorithms for gamma and inverse gaussian models. In each algorithm, it is straightforward to see that the divergence measure for gamma and inverse gaussian models is that based on KL divergence. For the NMF problem , when we used equation 2.10 the kernel of KL divergence for the gamma model can be written as
formula
3.2
This is commonly referred to as the Itakuro-Saito divergence (Cichocki et al., 2009; Févotte & Idier, 2011). Heuristic updates for W and H are given by
formula
3.3
formula
3.4
and MM updates are given by
formula
3.5
formula
3.6
Similarly, when we use equation 2.13, the kernel of KL divergence for the inverse gaussian model can be written as
formula
3.7
Heuristic updates for W and H are given by
formula
3.8
formula
3.9
and MM updates are given by
formula
3.10
formula
3.11
Using results from Cheung and Tresch (2005), Cichocki et al. (2006), and Févotte and Idier (2011), it is straightforward to obtain the above update rules in each specific case. For consistency, we use the notation KLHG, KLMMG, KLHIG, and KLMMIG to represent the heuristic and MM algorithms for gamma and inverse gaussian models, respectively.

3.2.  Proposed Algorithms.

In this section, we propose two novel NMF algorithms based on dual KL divergence, one each for the gamma and inverse gaussian models, and one algorithm based on J-divergence for the gamma model. We use the notation KLdG, KLdIG, and JG, respectively, to denote these three algorithms presented in theorems 1 to 3. Closed-form multiplicative update rules for W and H are provided for each, while proofs of monotonicity of updates are detailed in the appendix.

3.2.1.  Gamma Model: Algorithm Based on Dual KL Divergence.

Using equation 2.11, the kernel of dual KL divergence for the gamma model can be written as
formula
3.12
Theorem 1. 

The divergence KLdG(V, WH) in equation 3.12 is nonincreasing under the multiplicative update rules for W and H given by equations 3.13 and 3.14. It is also invariant under these updates if and only if W and H are at a stationary point of the divergence.

Proof. 

See the appendix.

Update rules for H and W are
formula
3.13
formula
3.14

3.2.2.  Gamma Model: Algorithm Based on J Divergence.

When equation 2.12 is used, the kernel of J divergence for the gamma model can be written as
formula
3.15
Theorem 2. 

The divergence JG(V, WH) defined in equation 3.15 is nonincreasing under the multiplicative update rules for W and H given by equations 3.16 and 3.17. It is also invariant under these updates if and only if W and H are at a stationary point of the divergence.

Proof. 

See the appendix.

Update rules for H and W are
formula
3.16
formula
3.17

3.2.3.  Inverse Gaussian Model: Algorithm Based on Dual KL Divergence.

When equation 2.14 is used, the kernel of dual KL divergence for the inverse gaussian model can be written as
formula
3.18
Theorem 3. 

The divergence KLdIG(V, WH) defined in equation 3.18 is nonincreasing under the multiplicative update rules for W and H given by equations 3.19 and 3.20. It is also invariant under these updates if and only if W and H are at a stationary point of the divergence.

Proof. 

See the appendix.

Update rules for H and W are
formula
3.19
formula
3.20

When equations 2.13 and 2.14 are used, J-divergence for the inverse gaussian model can be written in terms of V and WH. However, we note that closed-form multiplicative updates cannot be obtained using the EM approach, and the monotonicity of updates cannot be theoretically established.

Remark. 

The divergences 3.12, 3.15, and 3.18 and their corresponding update rules for W and H (equations 3.13 and 3.14, 3.16 and 3.17, and 3.19 and 3.20, respectively) contain Vij in the denominator of various terms. Theoretically, this should not cause any numerical issues (such as division by zero) since Vij>0 for both gamma and inverse gaussian models (i.e., zero is not in their domain). However, in a practical setting, due to data preprocessing, a few zero entries may sometimes occur in the input matrix V. In such cases, it is reasonable to set the zero entries to the smallest nonzero entry in V.

4.  Model Selection and Measuring Goodness of Fit

Starting with random initial values for W and H, the multiplicative update rules for any given NMF algorithm outlined in section 3 ensure monotonicity of updates for that run; however, the algorithm may not necessarily converge to the same solution on each run. In general, NMF algorithms are prone to this problem of local minima. For a given NMF algorithm and a prespecified rank r factorization, the corresponding divergence (or reconstruction error) computed at the final converged values of W and H for a set of random initial values can be used directly in model selection and to measure goodness of fit. One solution is to use the factorization from the run that results in the best reconstruction (quantified by minimum reconstruction error across multiple runs) for evaluation using different quantities. Below, we define two quantities of interest for this purpose based on algorithm-specific minimum reconstruction error E.

4.1.  Proportion of Explained Variation.

We propose several new measures to quantify the variation explained by the various algorithms for signal-dependent noise discussed in this letter. For each prespecified rank r, the proportion of explained variation (or empirical uncertainty), R2, is dependent on the particular algorithm and model used in the factorization. For the gaussian NMF algorithm, R2 is the well-known quantity given by
formula
4.1
where RSS is the residual sum of squares, SST is the total sum of squares, and is the minimum reconstruction error (E), calculated based on the kernel of the gaussian likelihood . The gaussian likelihood for NMF is obtained using equation 2.9 and was first proposed by Lee and Seung (2001).

For each algorithm, R2 is computed based on the corresponding minimum reconstruction error (E), as listed in Table 2. In the R2 column of this table, the numerator of each quantity within parentheses (other than the gaussian) is the minimum reconstruction error (E) calculated using equations 3.2, 3.7, 3.12, 3.15 and 3.18, as appropriate. The quantity (WH)ij in each numerator is the (i, j)th entry of the reconstructed matrix WH (also obtained as for a given rank r). In the corresponding denominator of each quantity, each entry of the reconstructed matrix WH is replaced by the grand mean of all entries of the input matrix V, . The underlying principle in the calculation of R2 is that the algorithm-specific reconstruction error E quantifies the performance of the model as determined by the entries (WH)ij, while in the absence of the model , the best approximation of (WH)ij is provided simply by the grand mean of all observations in the data. This is a direct extension of the definition of R2 in equation 4.1 for the gaussian model to nonlinear models such as the gamma and inverse gaussian. The algorithm-specific R2 measures the proportionate reduction in uncertainty due to the inclusion of W and H and, therefore, can be interpreted in terms of information content of the data (see Cameron & Windmeijer, 1996, 1997, for more details).

Table 2:
Algorithm-Specific Proportion of Explained Variation R2.
AlgorithmR2
Gaussian  
KLHG, KLMMG  
KLdG  
JG  
KLHIG, KLMMIG  
KLdIG  
AlgorithmR2
Gaussian  
KLHG, KLMMG  
KLdG  
JG  
KLHIG, KLMMIG  
KLdIG  

4.2.  Akaike Information Criterion (AIC).

For a particular algorithm and a prespecified rank r, AIC is given by
formula
4.2
where E is the corresponding minimum reconstruction error, is the total number of parameters estimated in the model for a input matrix V, and and for the gaussian, gamma, and inverse gaussian models, respectively. The model (rank r factorization) that results in the smallest AIC is chosen as the optimal model. The calculation of algorithm-specific E is detailed in section 4.1, and the determination of is outlined in section 5.

5.  Implementation of Algorithms on EMG Data

In this section, we present a detailed application of NMF algorithms based on signal-dependent noise in the analysis of the EMG data described in section 2. Time-invariant muscle synergies were extracted from each of the intact and deafferented EMG data sets of each frog using each of the eight NMF algorithms described earlier, including one based on normally distributed noise (gaussian), four based on gamma noise (including KLHG, KLMMG, KLdG, and JG), and three based on inverse gaussian (IG) noise (including KLHIG, KLMMIG, and KLdIG). The NMF update rules were implemented using Matlab. It should be noted that none of the preprocessed data sets contained zero entries. For every extraction, the muscle synergies (W) and their associated time-varying activation coefficients (H) were initialized with random matrices whose components were uniformly distributed between 0 and 1. Convergence was defined as having 20 consecutive iterations with a change of R2 smaller than 10−8 (with R2 for each algorithm defined in Table 2), but if convergence was not reached within 500 iterations, the extraction was terminated. The number of muscle synergies r extracted from each data set was successively increased from 1 to 13; at each number, extraction was repeated 20 times, each time with different random initial matrices.

AIC was calculated as follows. Let c denote the parameter , or depending on the model. In the specification of the divergence for each algorithm (see section 2.2.2, equations 2.9 to 2.14), we assumed that c=1 without loss of generality. In order to ensure that the EMG data fit this assumption, a global test of the null hypothesis H0:c=1 against the two-sided alternative was performed for each model. The mean-variance relationship for the gamma and inverse gaussian models can be written using equations 2.6 and 2.8, respectively, and is described in detail in section 2.2.3. This relationship was used to obtain an estimate of or in these models. For the gaussian model, was estimated using the approach described in Morup and Hansen (2009). In addition, these parameters were estimated using standard maximum likelihood methods. In each case, the estimate of c was approximately 1, and the 95% confidence interval for this estimate included 1, (p-values for these tests ranged from 0.15 to 0.77) thereby providing strong evidence in favor of the null hypothesis c=1. Based on this empirical evidence, c was taken to be 1 and the appropriate value of was used in equation 4.2 for each model. The best model order was selected by identifying the rank r giving the minimum AIC for each data set and algorithm. Table 3 lists the dimensionality of the data set, that is, number of columns n in ψ in equation 4.2, for each frog and behavior.

Table 3:
Dimensionality of the Data Set for Each Frog and Behavior: Number of Columns n in for Calculating AIC in Equation 4.2.
FrogIntact JumpDeafferented JumpIntact SwimDeafferented Swim
30,182 51,296 28,664 20,828 
13,445 12,202 47,341 35,136 
15,820 4,638 48,033 37,997 
14,823 11,963 37,833 26,598 
FrogIntact JumpDeafferented JumpIntact SwimDeafferented Swim
30,182 51,296 28,664 20,828 
13,445 12,202 47,341 35,136 
15,820 4,638 48,033 37,997 
14,823 11,963 37,833 26,598 

Since for this application we are primarily interested in the ability of each algorithm to identify features shared between the intact and deafferented EMG data sets (or features interpretable as units of CPGs), the performance of each algorithm was assessed by the similarity between the intact and deafferented muscle synergies, quantified with two measures. The first measure used was the scalar product between best-matching pairs of intact and deafferented synergies, calculated after the synergies were normalized to unit vectors. The second measure used was the cosine of the principal angles between the subspaces spanned by the intact and deafferented synergy sets (Golub & Van Loan, 1983). Both measures were used in Cheung et al. (2005).

5.1.  NMF Algorithms Based on Signal-Dependent Noise Outperformed Gaussian NMF.

In analysis of motor patterns from natural behaviors, it has remained difficult to determine, a priori, the number of muscle synergies composing the data set. Most previous studies on muscle synergies have relied on ad hoc measures to determine this number either by locating the cusp of the R2 curve plotted against the rank r (d'Avella, Saltiel, & Bizzi, 2003; Cheung et al., 2005; Tresch, Cheung, & d'Avella, 2006) or by finding the minimum number of synergies that produced an R2 greater than a certain arbitrary threshold (Cheung et al., 2012). Here, we explored using the AIC as a principled, objective measure of selecting the model order that best described the data without over-fitting (Akaike, 1987). For every algorithm and data set, we plotted the AIC against the number of synergies extracted r, and the preferred model order was indicated by the number at which the curve attained a minimum AIC (indicated in Figure 2 by the asterisks).

Figure 2:

Selecting the number of muscle synergies for the JG algorithm using AIC. To determine the model order, the number of muscle synergies extracted was successively increased from 1 to 13. At each number of synergies, the AIC was calculated using equation 4.2. (A) Plot of AIC against the number of muscle synergies extracted for both the intact (black solid) and deafferented (dotted) jump (4 frogs; mean SD). (B) Plot of AIC for intact (solid black) and deafferented (dotted) swim. The model order with minimum AIC was found to be 3 for jump and 4 for swim (*).

Figure 2:

Selecting the number of muscle synergies for the JG algorithm using AIC. To determine the model order, the number of muscle synergies extracted was successively increased from 1 to 13. At each number of synergies, the AIC was calculated using equation 4.2. (A) Plot of AIC against the number of muscle synergies extracted for both the intact (black solid) and deafferented (dotted) jump (4 frogs; mean SD). (B) Plot of AIC for intact (solid black) and deafferented (dotted) swim. The model order with minimum AIC was found to be 3 for jump and 4 for swim (*).

The model order selected using the AIC was in general consistent across animals for all four behaviors (intact jump, deafferented jump, intact swim, and deafferented swim) and all algorithms. For the gaussian and all gamma-based algorithms, the selected model order in each behavior differed at most by only 1 across frogs; for the IG-based algorithms, the selected order differed across frogs by 2 to 3 in most instances and by 4 only in one instance (see Figure 3). While IG algorithms (KLHIG, KLMMIG, and KLdIG) consistently indicated a rank of 8 to 13 muscle synergies, gaussian, KLdG, and KLHG suggested a rank of one to two synergies. The JG algorithm indicated a rank of three synergies for intact and deafferented jump (see Figures 2A and 3), and four synergies for intact and deafferented swim (see Figures 2B and 3). A previous study (d'Avella et al., 2003) has argued, based on a detailed kinematic analysis, that there are at least three muscle synergies underlying frog hind-limb kicking, with two synergies controlling kick direction during limb extension and one for executing limb flexion. It thus appears that the model orders selected for the JG algorithm are the most physiologically interpretable. We further verified that at these ranks, all algorithms returned synergies that explained the EMGs with an R2 of at least 85% (see Table 4). In addition, all three proposed algorithms returned synergies that explained the EMGs with an R2 of at least 93%, a significantly higher fraction compared to that of existing signal-dependent noise algorithms for which R2 values ranged from 85% to 93%. In the ensuing analysis, we will measure the performances of all algorithms at the ranks determined from the results of the JG algorithm.

Figure 3:

Number of muscle synergies selected for the different NMF algorithms. For each behavior (A, intact and deafferented jump; B, intact and deafferented swim) and each algorithm, the number of muscle synergies selected for each frog was determined by selecting the rank with minimum AIC. Note that the selected numbers for all behaviors and algorithms were quite consistent across animals.

Figure 3:

Number of muscle synergies selected for the different NMF algorithms. For each behavior (A, intact and deafferented jump; B, intact and deafferented swim) and each algorithm, the number of muscle synergies selected for each frog was determined by selecting the rank with minimum AIC. Note that the selected numbers for all behaviors and algorithms were quite consistent across animals.

Table 4:
The Proportion of Explained Variation (R2) Achieved by the NMF Algorithms in Four Frog Behaviors at the Rank Determined by the JG Algorithm.
Intact JumpDeafferented JumpIntact SwimDeafferented Swim
Algorithmr=3r=3r=4r=4
Gaussian 87.69 1.21 87.38 1.04 85.43 1.83 87.89 4.51 
KLHG 91.97 1.48 90.58 1.17 87.17 1.13 90.11 2.47 
KLMMG 91.97 1.48 90.58 1.17 87.16 1.13 90.11 2.47 
KLdG 98.93 0.35 98.74 0.11 96.46 0.41 97.62 0.77 
JG 97.85 0.64 97.41 0.29 93.56 0.55 95.50 1.38 
KLHIG 90.97 2.06 88.72 0.58 86.67 0.87 89.49 2.30 
KLMMIG 90.89 2.07 88.68 0.52 86.50 0.83 89.43 2.25 
KLdIG 99.79 0.11 99.77 0.05 98.94 0.30 99.37 0.22 
Intact JumpDeafferented JumpIntact SwimDeafferented Swim
Algorithmr=3r=3r=4r=4
Gaussian 87.69 1.21 87.38 1.04 85.43 1.83 87.89 4.51 
KLHG 91.97 1.48 90.58 1.17 87.17 1.13 90.11 2.47 
KLMMG 91.97 1.48 90.58 1.17 87.16 1.13 90.11 2.47 
KLdG 98.93 0.35 98.74 0.11 96.46 0.41 97.62 0.77 
JG 97.85 0.64 97.41 0.29 93.56 0.55 95.50 1.38 
KLHIG 90.97 2.06 88.72 0.58 86.67 0.87 89.49 2.30 
KLMMIG 90.89 2.07 88.68 0.52 86.50 0.83 89.43 2.25 
KLdIG 99.79 0.11 99.77 0.05 98.94 0.30 99.37 0.22 

The seven NMF algorithms based on signal-dependent noise outperformed the gaussian NMF in their ability to identify features shared between the intact and deafferented EMG data sets. This is indicated by the generally higher similarity between the intact and deafferented muscle synergy sets, measured by both the scalar product (see Figures 4A and 4C) and the cosine of principal angle (see Figures 4B and 4D), when the nongaussian NMFs were applied. In particular, for both similarity measures, performance of the JG algorithm exceeded that of the gaussian algorithm in three of four frogs in the jump data sets (frogs 2, 3, 4; see Figures 4A and 4B) and also in three of four frogs in the swim data sets (frogs 1, 2, 3; see Figures 4C and 4D). Overall, the three best-performing algorithms in terms of these measures were JG (mean scalar product = 0.8698; ); KLdIG (0.8695), and KLHG (0.8650). The worst-performing algorithm was gaussian (0.7648).

Figure 4:

Signal-dependent noise NMFs outperformed the gaussian NMF. In this application, we are primarily interested in each algorithm's ability to identify structures shared between the intact and deafferented data sets; thus, our measures of algorithm performance are based on quantifying the similarity between the intact and deafferented muscle synergies. For both the scalar-product (A and C) and principal angle (B and D) measures, overall the seven NMFs based on signal-dependent noise outperformed the gaussian NMF in their ability to extract features shared between data sets. In each graph, the level of similarity achieved by the gaussian algorithm (black) is marked by a horizontal black dotted line for ease of visual inspection.

Figure 4:

Signal-dependent noise NMFs outperformed the gaussian NMF. In this application, we are primarily interested in each algorithm's ability to identify structures shared between the intact and deafferented data sets; thus, our measures of algorithm performance are based on quantifying the similarity between the intact and deafferented muscle synergies. For both the scalar-product (A and C) and principal angle (B and D) measures, overall the seven NMFs based on signal-dependent noise outperformed the gaussian NMF in their ability to extract features shared between data sets. In each graph, the level of similarity achieved by the gaussian algorithm (black) is marked by a horizontal black dotted line for ease of visual inspection.

Table 5 lists the proportion of explained variance (R2) achieved by the NMF algorithms at the ranks with minimum AIC. For every algorithm, the number of muscle synergies underlying each behavior of each frog was determined by selecting the rank that resulted in the smallest AIC values. All R2 values shown are averages across frogs (N = 4; mean SD).

Table 5:
The Proportion of Explained Variation (R2) Achieved by the NMF Algorithms at the Ranks with Minimum AIC.
RankR2
AlgorithmBehavior(N=4; median SD)(N=4; mean SD)
Gaussian Intact jump 3 0.58 85.16 4.08 
 Deafferented jump 2 0.50 82.91 3.52 
 Intact swim 1 0.00 56.78 4.89 
 Deafferented swim 1 0.00 56.58 12.33 
KLHG Intact jump 2 0.50 86.66 3.64 
 Deafferented jump 2 0.50 85.22 1.09 
 Intact swim 2 0.00 76.66 3.68 
 Deafferented swim 2 0.00 80.00 5.07 
KLMMG Intact jump 2 0.50 86.66 3.64 
 Deafferented jump 2 0.50 85.22 1.09 
 Intact swim 2 0.00 76.66 3.69 
 Deafferented swim 2 0.00 80.00 5.07 
KLdG Intact jump 2 0.50 98.21 0.77 
 Deafferented jump 2 0.50 98.05 0.17 
 Intact swim 2 0.50 92.81 2.80 
 Deafferented swim 2 0.50 94.78 2.70 
JG Intact jump 3 0.50 98.05 0.79 
 Deafferented jump 3 0.50 97.61 0.47 
 Intact swim 4 0.50 93.03 1.28 
 Deafferented swim 4 0.50 95.18 1.05 
KLHIG Intact jump 10 0.96 99.07 0.38 
 Deafferented jump 11 1.41 98.84 0.54 
 Intact swim 11 0.82 99.25 0.34 
 Deafferented swim 12 1.29 99.61 0.30 
KLMMIG Intact jump 9 0.82 98.41 0.45 
 Deafferented jump 10 1.73 98.67 0.49 
 Intact swim 11 0.96 98.50 0.46 
 Deafferented swim 10 1.50 98.95 0.43 
KLdIG Intact jump 10 1.50 99.98 0.02 
 Deafferented jump 10 1.41 99.98 0.01 
 Intact swim 11 0.58 99.92 0.03 
 Deafferented swim 11 1.41 99.97 0.02 
RankR2
AlgorithmBehavior(N=4; median SD)(N=4; mean SD)
Gaussian Intact jump 3 0.58 85.16 4.08 
 Deafferented jump 2 0.50 82.91 3.52 
 Intact swim 1 0.00 56.78 4.89 
 Deafferented swim 1 0.00 56.58 12.33 
KLHG Intact jump 2 0.50 86.66 3.64 
 Deafferented jump 2 0.50 85.22 1.09 
 Intact swim 2 0.00 76.66 3.68 
 Deafferented swim 2 0.00 80.00 5.07 
KLMMG Intact jump 2 0.50 86.66 3.64 
 Deafferented jump 2 0.50 85.22 1.09 
 Intact swim 2 0.00 76.66 3.69 
 Deafferented swim 2 0.00 80.00 5.07 
KLdG Intact jump 2 0.50 98.21 0.77 
 Deafferented jump 2 0.50 98.05 0.17 
 Intact swim 2 0.50 92.81 2.80 
 Deafferented swim 2 0.50 94.78 2.70 
JG Intact jump 3 0.50 98.05 0.79 
 Deafferented jump 3 0.50 97.61 0.47 
 Intact swim 4 0.50 93.03 1.28 
 Deafferented swim 4 0.50 95.18 1.05 
KLHIG Intact jump 10 0.96 99.07 0.38 
 Deafferented jump 11 1.41 98.84 0.54 
 Intact swim 11 0.82 99.25 0.34 
 Deafferented swim 12 1.29 99.61 0.30 
KLMMIG Intact jump 9 0.82 98.41 0.45 
 Deafferented jump 10 1.73 98.67 0.49 
 Intact swim 11 0.96 98.50 0.46 
 Deafferented swim 10 1.50 98.95 0.43 
KLdIG Intact jump 10 1.50 99.98 0.02 
 Deafferented jump 10 1.41 99.98 0.01 
 Intact swim 11 0.58 99.92 0.03 
 Deafferented swim 11 1.41 99.97 0.02 

It is clear from the results shown in Tables 4 and 5 that all three proposed algorithms outperformed existing algorithms in terms of fraction of explained variation (R2), both at the ranks with minimum AIC and at the ranks determined by the JG algorithm. A closer look also revealed that the variability of this fraction (estimated by the standard deviation) was significantly lower for the proposed algorithms relative to existing methods, indicating a higher overall confidence level in the variation explained by these methods. The JG algorithm provided a much better balance between R2 and choice of rank based on minimum AIC compared to any other algorithm. The ranks chosen by the KLdG algorithm based on minimum AIC were similar to those of existing gamma-based algorithms; however, this algorithm was able to explain a much higher fraction of variation in the data. The KLdIG algorithm explained the maximum variation (highest overall R2) among all algorithms, while the gaussian algorithm provided the smallest R2 and rank. Furthermore, the variability of R2 was also the highest for the gaussian algorithm, suggesting an overall lower confidence level in the variation explained by this algorithm.

5.2.  Muscle Synergies Extracted by the JG Algorithm Were Physiologically Interpretable.

In this section, we compare swim muscle synergies extracted using the standard gaussian algorithm with those identified by the JG algorithm in one specific individual (frog 2) and illustrate how the latter set could be more physiologically interpretable. In the extraction results returned by the gaussian formulation, a very high similarity between the pre- and postdeafferentation synergies was observed in two of the synergy pairs (scalar product >0.90; see Figure 5A, synergies 1 to 2), a moderate similarity in one synergy pair (scalar product = 0.90; see Figure 5A, synergy 3), and total dissimilarity in the last pair (scalar product = 0.06; see Figure 5A, synergy 4). By contrast, the JG algorithm found three synergy pairs with high similarity (scalar product >0.90; see Figure 5B, synergies 1 to 3); in the last pair, the similarity was modest (scalar product = 0.62; see Figure 5B, synergy 4), but the muscles active in both the intact and deafferented synergy vectors were the same (RI, AD, SM, and ST). Overall, the synergy extraction results from this frog demonstrate that the JG algorithm, derived from a signal-dependent noise assumption, is better able to discover structures preserved after deafferentation than the traditional gaussian algorithm.

The muscular compositions of the synergies returned by the JG algorithm could also be biomechanically interpreted, Synergy 3 (see Figure 5B), for instance, was composed of the hip extensor SM; knee extensors VI, RA, and VE; and the ankle extensor GA. Examination of the time-varying coefficients associated with this synergy revealed that it was active only during the extension phase of every swim cycle; thus, it is likely that muscle synergy 3 functions to propel the animal forward through extension of the hip, knee, and ankle joints. Synergy 1 (see Figures 5A and 5B), discovered by both the gaussian and JG algorithms, consisted of the hip flexors IP and BI; synergy 2 (see Figures 5A and 5B), on the other hand, consisted primarily of the ankle flexors TA and PE, and the hip/knee flexor SA. It is no surprise that both of these synergies were active during the flexion phase of every swim cycle.

Figure 5:

Muscle synergies extracted by the JG algorithm were physiologically interpretable. (A) Intact (black) and deafferented (white) muscle synergies for swimming (frog 2) returned by the gaussian algorithm. The scalar product similarity between each synergy pair is indicated above the pair. The intact and deafferented synergies for pair 4 were totally dissimilar. (B) Intact (black) and deafferented (white) muscle synergies for swimming (frog 2) returned by the JG algorithm. Here, even in the least similar pair (pair 4, scalar product = 0.62), the sets of muscles found to be active in the intact and deafferented synergies were still identical. (C) The correlation coefficient between the activation of muscle synergy 3 (the extension synergy) and those of muscle synergies 1, 2, and 4, respectively, before (black) and after (white) deafferentation. Note that the correlation between synergies 3 and 4 increased five-fold after deafferentation. This suggests that sensory feedback is essential in triggering or maintaining the activation of synergy 4 during the flexion phase of the swim cycle.

Figure 5:

Muscle synergies extracted by the JG algorithm were physiologically interpretable. (A) Intact (black) and deafferented (white) muscle synergies for swimming (frog 2) returned by the gaussian algorithm. The scalar product similarity between each synergy pair is indicated above the pair. The intact and deafferented synergies for pair 4 were totally dissimilar. (B) Intact (black) and deafferented (white) muscle synergies for swimming (frog 2) returned by the JG algorithm. Here, even in the least similar pair (pair 4, scalar product = 0.62), the sets of muscles found to be active in the intact and deafferented synergies were still identical. (C) The correlation coefficient between the activation of muscle synergy 3 (the extension synergy) and those of muscle synergies 1, 2, and 4, respectively, before (black) and after (white) deafferentation. Note that the correlation between synergies 3 and 4 increased five-fold after deafferentation. This suggests that sensory feedback is essential in triggering or maintaining the activation of synergy 4 during the flexion phase of the swim cycle.

The activation pattern of synergy 4 identified by JG (see Figure 5B) was more complex. During the intact state, it was primarily active during limb flexion; after deafferentation it was activated only during limb extension. Consistent with this switch of activation phase for this synergy after the loss of sensory feedback, the correlation coefficient between the activation coefficients of synergy 4 and those of the extension synergy 3 increased five-fold after deafferentation (see Figure 5C). Since three muscles in this synergy—RI, SM, and ST—have both hip extension and knee flexion actions, it is possible that before deafferentation, this synergy executes knee flexion, while after deafferentation, it aids limb extension. It thus appears that sensory feedback functions both to inhibit its activation during extension and facilitates or triggers its activation during flexion. Such an inference about the contribution of afferents to the inhibition and activation of this muscle synergy would be difficult with the synergy sets obtained by the gaussian NMF (see Figure 5A) given that the gaussian algorithm failed to discover this synergy from the deafferented data set.

5.3.  Comparison of Results Using Algorithms Derived from the Same Noise Distribution.

In the preceding sections, we compared the performance of various algorithms in extracting muscle synergies based on AIC, the fraction of explained variation (R2), their ability to identify features shared between the deafferented and intact EMG data (measured by the scalar dot product and cosine of the principal angle), and their physiological interpretability. In this section, we perform a comparison of muscle synergies extracted by different NMF algorithms from the same EMG data set in order to understand how the underlying noise assumption and cost function used in the NMF algorithm may affect the muscular compositions of the extracted synergies. Algorithms based on the same noise distribution but different cost functions tended to return similar muscle synergies. For instance, when we use KLHG and KLHIG as reference algorithms, the scalar product values (mean SD; over four frogs) between the synergies returned by gamma-based NMF algorithms and KLHG were (1) higher than those between the synergies returned by the gaussian NMF algorithm and KLHG and (2) higher than those between inverse gaussian-based NMF algorithms and KLHG (see Figure 6A). Similarly, the scalar product values (mean SD; over four frogs) between the synergies returned by inverse gaussian-based NMF algorithms and KLHIG were (1) higher than those between the synergies returned by the gaussian NMF algorithm and KLHIG and (2) higher than those between gamma-based NMF algorithms and KLHIG (see Figure 6B).

Figure 6:

Comparison of results using NMF algorithms derived from the same noise distribution. We performed a comparison of the muscle synergies extracted by different NMF algorithms from the same EMG data set in order to understand the effects of the NMF-noise distribution and the cost function employed on the muscular compositions of the extracted muscle synergies. (A) In each frog, the set of muscle synergies extracted by each algorithm was matched to the set returned by the gamma-based KLHG algorithm (*), and their similarity was quantified by the scalar product values averaged across the synergy set. Shown in the plot are values averaged across frogs (N = 4; mean SD). Values for the KLHG were 1.0 by definition. In this comparison, scalar product values from the gamma algorithms tended to be higher than those from the gaussian or IG-based algorithms. This difference is especially obvious for the intact jump and deafferented jump data sets. (B) Same as panel A, except that the comparison was performed by matching synergies of each algorithm to synergies returned by the IG-based KLHIG algorithm (*). In this comparison, scalar product values from IG-based algorithms tended to be higher than those from the gaussian or gammabased algorithms. Again, this difference is especially obvious for the intact jump and deafferented jump data sets.

Figure 6:

Comparison of results using NMF algorithms derived from the same noise distribution. We performed a comparison of the muscle synergies extracted by different NMF algorithms from the same EMG data set in order to understand the effects of the NMF-noise distribution and the cost function employed on the muscular compositions of the extracted muscle synergies. (A) In each frog, the set of muscle synergies extracted by each algorithm was matched to the set returned by the gamma-based KLHG algorithm (*), and their similarity was quantified by the scalar product values averaged across the synergy set. Shown in the plot are values averaged across frogs (N = 4; mean SD). Values for the KLHG were 1.0 by definition. In this comparison, scalar product values from the gamma algorithms tended to be higher than those from the gaussian or IG-based algorithms. This difference is especially obvious for the intact jump and deafferented jump data sets. (B) Same as panel A, except that the comparison was performed by matching synergies of each algorithm to synergies returned by the IG-based KLHIG algorithm (*). In this comparison, scalar product values from IG-based algorithms tended to be higher than those from the gaussian or gammabased algorithms. Again, this difference is especially obvious for the intact jump and deafferented jump data sets.

Our analysis shows that in our frog EMGs, algorithms derived from the same noise distributions tended to return similar muscle synergies. The noise distribution appears to play a critical role in determining the muscular compositions of the synergies extracted from the data. Similarly, the cost function (divergence measure) employed for formulating the update rules exerts its own influence on the extraction results. Indeed, the best rank (rank with minimal AIC) and R2 values from algorithms assuming the same noise distribution but employing different cost functions were still different (see Table 5). This is because these algorithms derived from different cost functions returned different activation coefficients (H). As an illustrative example, we present in Figure 7 the muscle synergies extracted by all eight algorithms from the deafferented jump EMGs of frog 2. In this case, the results produced by the four gamma-based NMF algorithms are nearly identical. However, it is important to note that the synergies extracted by the three IG-based NMF algorithms are quite different, exposing the activation of different muscles as determined by the choice of cost function.

Figure 7:

Both the noise distribution and the cost function employed for formulating the NMF update rules could influence the muscular compositions of the extracted muscle synergies. Here we show the muscle synergies extracted from one particular data set (frog 2, deafferented jump) by different NMF algorithms. The results returned by the four gamma-based algorithms were almost identical (as suggested by Figure 6). However, the gamma synergies were clearly different from the gaussian and IG-based synergies. Also, the muscle synergies returned by the three IG-based synergies were also somewhat different from each other. Thus, both the noise distribution and the cost function used for deriving the NMF update rules could influence the structures of the basis vectors extracted.

Figure 7:

Both the noise distribution and the cost function employed for formulating the NMF update rules could influence the muscular compositions of the extracted muscle synergies. Here we show the muscle synergies extracted from one particular data set (frog 2, deafferented jump) by different NMF algorithms. The results returned by the four gamma-based algorithms were almost identical (as suggested by Figure 6). However, the gamma synergies were clearly different from the gaussian and IG-based synergies. Also, the muscle synergies returned by the three IG-based synergies were also somewhat different from each other. Thus, both the noise distribution and the cost function used for deriving the NMF update rules could influence the structures of the basis vectors extracted.

6.  Evaluating NMF Algorithms on Simulated Data Sets

In this section, we present a detailed application of the proposed NMF algorithms to the analysis of simulated data. We implemented the algorithms on simulated data sets generated by known muscle synergies (W) and time-varying activation coefficients (H), so that the performance of each NMF algorithm can be evaluated by comparing the extracted results with the original W and H.

In our simulations, we are interested in how well each algorithm performs as a function of noise distribution and noise level in the data. For every distribution and noise amplitude tested, 10 simulated data sets were generated. Each data set, consisting of 15 muscles and 5000 time points, was produced by linearly combining five muscle synergies. The components of both W and H were drawn from a uniform distribution defined over (0, 1). The simulated data were then corrupted by one of the three noise types—gaussian, gamma, and inverse gaussian—at different noise magnitudes quantified by the signal-to-noise ratio (SNR), defined as
formula
where Vij is the original, uncorrupted data point and is the noise-corrupted data point. For gaussian noise with mean and variance , noise for each data point was generated by the Matlab function, normrnd, with and set to 0.04, 0.05, 0.07, 0.1, 0.2, 0.3, 0.5, 1.0, 1.5, and 2.0, respectively. These choices of produced data with an SNR ranging from 0.17 to 225. For gamma noise with mean (see equation 2.5), the Matlab function gamrnd was used, with and set to 0.1, 0.5, 0.25, 1.0, 2.5, 5.0, 10, 50, 100, 250, and 500, respectively (SNR of 0.1 to 500). For inverse gaussian noise with mean (see equation 2.7), noise for each data point was generated by combining the Matlab functions makedist and random, with and set to 0.1, 0.2, 0.5, 0.75, 1, 2, 5, 7, 10, 12, 14, 16, 18, 20, 30, 40, and 100, respectively (SNR of 0.14 to 139).

The eight NMF algorithms described in this letter (gaussian, KLHG, KLMMG, KLdG, JG, KLHIG, KLMMIG, and KLdIG) were then applied to each of the simulated data sets for extracting five muscle synergies. In every extraction, the NMF update rules were implemented using Matlab (R2013b). The W and H matrices were initialized with random components drawn from a uniform distribution over (0, 1). Convergence was defined as having 20 consecutive iterations with a change of algorithm-specific R2 (see Table 2) smaller than 10−8, but if convergence was not achieved within 500 iterations, the extraction was terminated. Extraction was repeated 20 times for each data set, each time with different initial random matrices. The extraction repetition with the smallest reconstruction error among the 20 repetitions was then selected for performance evaluation. The ability of each algorithm in identifying the muscle synergies was quantified by the scalar product between the original and extracted synergy vectors (after the synergies were normalized to unit vectors), averaged over the five synergies. For the activation coefficients, performance was assessed by the Pearson's correlation coefficient () between the components in the original H and those in the extracted H, again averaged over the five synergies.

For gaussian-noise data sets, the gaussian algorithm outperformed KLdG, JG, and all IG-based algorithms in the identification of both W and H (see Figure 8; *, p<0.05; Student's t-test). The superiority in performance of the gaussian NMF over all other algorithms was especially obvious for the extraction of H (see Figure 8B). For the extraction of W, however, the performances of gaussian, KLHG, and KLMMG were comparable (see Figure 8A).

Figure 8:

Gaussian NMF outperformed the signal-dependent noise NMFs in data sets corrupted by gaussian noise. We evaluated the performance of each algorithm in simulated data sets (N = 10) generated by known W (15 × 5 matrix) and H (5 × 5000 matrix), but corrupted by random gaussian noise at different signal-to-noise ratios (SNR). (A) Performance of NMF algorithms in identifying the basis vectors (W). Performance of each algorithm in each data set was quantified by the scalar product between the extracted vectors and the original vectors, averaged across the five basis vectors in the W matrix. Shown in the plot are mean scalar product values, defined as above, averaged across 10 simulated data sets. The gaussian NMF algorithm outperformed all IG-based NMF algorithms and two of the gamma-based NMF algorithms (KLdG and JG) over a wide range of SNR (*; Student's t-test; p<0.05). (B) Performance of NMF algorithms in identifying the coefficients (H). Performance of each algorithm in each data set was quantified by the Pearson's correlation coefficient () between the extracted coefficients and the original coefficients (over a total of 5 × 5000 = 25,000 values). Shown in the plot are values averaged across the 10 simulated data sets. The gaussian NMF algorithm outperformed all of the gamma- and IG-based NMF algorithms over almost all tested SNR (*; p<0.05).

Figure 8:

Gaussian NMF outperformed the signal-dependent noise NMFs in data sets corrupted by gaussian noise. We evaluated the performance of each algorithm in simulated data sets (N = 10) generated by known W (15 × 5 matrix) and H (5 × 5000 matrix), but corrupted by random gaussian noise at different signal-to-noise ratios (SNR). (A) Performance of NMF algorithms in identifying the basis vectors (W). Performance of each algorithm in each data set was quantified by the scalar product between the extracted vectors and the original vectors, averaged across the five basis vectors in the W matrix. Shown in the plot are mean scalar product values, defined as above, averaged across 10 simulated data sets. The gaussian NMF algorithm outperformed all IG-based NMF algorithms and two of the gamma-based NMF algorithms (KLdG and JG) over a wide range of SNR (*; Student's t-test; p<0.05). (B) Performance of NMF algorithms in identifying the coefficients (H). Performance of each algorithm in each data set was quantified by the Pearson's correlation coefficient () between the extracted coefficients and the original coefficients (over a total of 5 × 5000 = 25,000 values). Shown in the plot are values averaged across the 10 simulated data sets. The gaussian NMF algorithm outperformed all of the gamma- and IG-based NMF algorithms over almost all tested SNR (*; p<0.05).

In data sets with simulated gamma noise, the gamma- and IG-based algorithms performed equally well, and better than gaussian, in the identification of H over all tested noise magnitudes (see Figure 9B). For W identification, the gamma- and IG-based algorithms were similar in performance when the SNR was above (see Figure 9A). The gamma algorithms outperformed all other algorithms when noise magnitude was very high (see Figure 9A; +, p<0.05).

Figure 9:

Gamma-based NMF algorithms outperformed the gaussian NMF algorithm in data sets corrupted by gamma noise. We evaluated the performance of each algorithm in simulated data sets (N = 10) generated by known W (15 × 5 matrix) and H (5 × 5000 matrix), but corrupted by random gamma noise at different signal-to-noise ratios (SNR). (A) Performance of NMF algorithms in identifying the basis vectors (W). Performance of each algorithm in each data set was quantified by the scalar product between the extracted vectors and the original vectors, averaged across the five basis vectors in the W matrix. Shown in the plot are mean scalar product values, defined as above, averaged across 10 simulated data sets. Gamma-based algorithms outperformed the gaussian algorithm (but not the IG-based algorithms) at moderate noise magnitude (*; Student's t-test; p<0.05); but at high noise magnitudes, the gamma algorithms performed better than both gaussian- and IG-NMF algorithms (+; p<0.05). (B) Performance of the NMF algorithms in identifying the coefficients (H). Performance of each algorithm in each data set was quantified by the Pearson's correlation coefficient () between the extracted coefficients and the original coefficients (over a total of 5 × 5000 = 25,000 values). Shown in the plot are values averaged across the 10 simulated data sets. Gamma-based algorithms outperformed the gaussian algorithm, but not the IG-based algorithms, at moderate noise levels (*; p<0.05).

Figure 9:

Gamma-based NMF algorithms outperformed the gaussian NMF algorithm in data sets corrupted by gamma noise. We evaluated the performance of each algorithm in simulated data sets (N = 10) generated by known W (15 × 5 matrix) and H (5 × 5000 matrix), but corrupted by random gamma noise at different signal-to-noise ratios (SNR). (A) Performance of NMF algorithms in identifying the basis vectors (W). Performance of each algorithm in each data set was quantified by the scalar product between the extracted vectors and the original vectors, averaged across the five basis vectors in the W matrix. Shown in the plot are mean scalar product values, defined as above, averaged across 10 simulated data sets. Gamma-based algorithms outperformed the gaussian algorithm (but not the IG-based algorithms) at moderate noise magnitude (*; Student's t-test; p<0.05); but at high noise magnitudes, the gamma algorithms performed better than both gaussian- and IG-NMF algorithms (+; p<0.05). (B) Performance of the NMF algorithms in identifying the coefficients (H). Performance of each algorithm in each data set was quantified by the Pearson's correlation coefficient () between the extracted coefficients and the original coefficients (over a total of 5 × 5000 = 25,000 values). Shown in the plot are values averaged across the 10 simulated data sets. Gamma-based algorithms outperformed the gaussian algorithm, but not the IG-based algorithms, at moderate noise levels (*; p<0.05).

In data sets corrupted by inverse gaussian noise, for W identification, not surprisingly the IG-based algorithms outperformed the gamma-based algorithms (see Figure 10A; *), which in turn outperformed the gaussian (see Figure 10A; +, *). For H identification, while the performances of all signal-dependent noise NMFs were almost indistinguishable, they clearly did much better than the gaussian (see Figure 10B, *).

Figure 10:

The inverse gaussian NMF algorithms outperformed the gaussian- and gamma-based NMF algorithms in data sets corrupted by inverse gaussian noise. We evaluated the performance of each algorithm in simulated data sets (N = 10) generated by known W (15 × 5 matrix) and H (5 × 5000 matrix), but corrupted by random inverse gaussian (IG) noise at different signal-to-noise ratios (SNR). (A) Performance of NMF algorithms in identifying the basis vectors (W). Performance of each algorithm in each data set was quantified by the scalar product between the extracted vectors and the original vectors, averaged across the five basis vectors in the W matrix. Shown in the plot are mean scalar product values, defined as above, averaged across 10 simulated data sets. At moderate noise levels, IG-based algorithms clearly outperformed both the gaussian and gamma algorithms (*; Student's t-test; p<0.05); at high and low noise levels, IG-based algorithms still performed better than the gaussian (but not the gamma) algorithm (+; p<0.05). (B) Performance of the NMF algorithms in identifying the coefficients (H). Performance of each algorithm in each data set was quantified by the Pearson's correlation coefficient () between the extracted coefficients and the original coefficients (over a total of 5 × 5000 = 25,000 values). Shown in the plot are values averaged across the 10 simulated data sets. IG-based algorithms outperformed the gaussian NMF over a wide range of SNR (*; p<0.05).

Figure 10:

The inverse gaussian NMF algorithms outperformed the gaussian- and gamma-based NMF algorithms in data sets corrupted by inverse gaussian noise. We evaluated the performance of each algorithm in simulated data sets (N = 10) generated by known W (15 × 5 matrix) and H (5 × 5000 matrix), but corrupted by random inverse gaussian (IG) noise at different signal-to-noise ratios (SNR). (A) Performance of NMF algorithms in identifying the basis vectors (W). Performance of each algorithm in each data set was quantified by the scalar product between the extracted vectors and the original vectors, averaged across the five basis vectors in the W matrix. Shown in the plot are mean scalar product values, defined as above, averaged across 10 simulated data sets. At moderate noise levels, IG-based algorithms clearly outperformed both the gaussian and gamma algorithms (*; Student's t-test; p<0.05); at high and low noise levels, IG-based algorithms still performed better than the gaussian (but not the gamma) algorithm (+; p<0.05). (B) Performance of the NMF algorithms in identifying the coefficients (H). Performance of each algorithm in each data set was quantified by the Pearson's correlation coefficient () between the extracted coefficients and the original coefficients (over a total of 5 × 5000 = 25,000 values). Shown in the plot are values averaged across the 10 simulated data sets. IG-based algorithms outperformed the gaussian NMF over a wide range of SNR (*; p<0.05).

Overall, the simulation results highlight the need for using the NMF algorithm derived from a noise distribution that matches the noise type of the data for the most accurate identification of both W and H. However, under certain conditions, even when the noise assumed by the NMF algorithm and the data noise type do not completely agree, the extracted results may still contain substantial information about the underlying data structure. We have seen, for instance, that in data with gamma noise, IG-based NMF algorithms could identify muscle synergies as well as gamma-based NMF algorithms could. It should be noted that even when the identified W is reasonably close to the original generating bases, the H identified by the same algorithm may not be as accurate (and vice versa). For example, in data with gaussian noise, at the gamma-based KLMMG algorithm performed at the same level as the gaussian algorithm for W identification and returned muscle synergies that matched the originals with scalar product >0.8 (see Figure 8A); however, for H, the extraction results from KLMMG were not only much worse than gaussian but also matched the original quite poorly (see Figure 8B).

7.  Some Recommendations

As argued by our extraction results from simulated data (see Figures 8, 9, and 10), for a data set with known noise properties, using NMF algorithms derived from a noise distribution that matches that of the data should yield the most accurate estimations of both W and H. The noise distribution of the data can be determined using the exploratory approach outlined in section 2.2.3 for the EMG data presented in this letter. However, if the noise characteristics of the data are not known or cannot be reasonably determined, it is preferable to first evaluate the algorithms on a data set with a clear prediction, based on the biology of the processes generating the data, of what the underlying W or H could be, and see which algorithm produces results that best match such predictions. The best-performing algorithm can then be used in other data sets of a similar nature for a further understanding of the biology. The noise distribution assumed by this best-performing algorithm would allow an understanding of the noise structure of the data set. For the frog EMG data presented here, two of the three best-performing algorithms (i.e., algorithms that discovered the most number of synergies shared between the intact and deafferented data, or features interpretable as CPG components, while explaining most of the variation in the data) were derived from the gamma noise distribution. These are the JG and KLdG algorithms while KLdIG is the third best-performing algorithm. Thus, it is evident that signal-dependent noise in the EMGs is closer to the gamma than to the inverse gaussian distribution. This also agrees with the empirical observation and conclusion in section 2.2.3 that the gamma model provides a good fit to the mean-SD relationship of the EMG data.

Although we focused primarily on NMF algorithms for handling EMG data with signal-dependent noise in this letter, other matrix factorization algorithms have been used for extracting muscle synergies in addition to NMF (Tresch et al., 2006). In a variety of applications, NMF has been shown to provide a parts-based local representation of the data in contrast to the holistic representation provided by vector quantization and the distributed representation provided by principal component analysis (PCA) (Devarajan, 2008). PCA is based on the gaussian model and requires nonoverlapping, orthogonal components with mixed signs. On the other hand, independent component analysis (ICA) seeks a linear representation of nongaussian data such that the resulting components are statistically independent (Hyvärinen & Oja, 2000; Devarajan, 2011). The representation provided by ICA has been shown to capture the essential structure of the data in various applications involving blind source separation. Independence implies uncorrelatedness, and in the case of the gaussian distribution, they are equivalent, implying independent principal components. Thus, PCA and ICA require different but stronger assumptions, particularly with regard to application to EMG data. In general, PCA provides dimensionality reduction, while ICA results in perceptually relevant components. NMF provides interpretable components; however, it is limited by the nonnegativity requirement on the input data and the resulting components. From an exploratory data analysis perspective, it is important to note that each of these methods comes with its owns merits and demerits and that the extent of its usefulness depends on the specific application at hand. When data occur naturally on the nonnegative scale such as the EMG signals presented in this letter, it appears more intuitive and reasonable to apply a factorization that retains the nonnegativity requirement on the resulting components (muscle synergies). These nonnegativity constraints are compatible with the intuitive notion of combining parts to form a whole. In NMF, these components are additive, linear combinations of the parts that are overlapping and nonorthogonal. The resulting “parts” extracted by NMF from the EMG data can naturally be interpreted as representations of motor primitives—or basic modules of motor control—whose existence has been demonstrated in physiological experiments (Bizzi & Cheung, 2013). Moreover, the extension of this approach to nongaussian models described in this letter is particularly relevant for applications involving signal dependent noise. Such modeling flexibility is not provided by other methods.

8.  Summary and Conclusion

In this letter, we proposed a comprehensive extension of methods for handling data with signal-dependent noise in NMF. We outlined three novel algorithms based on dual KL and J-divergence for the gamma and inverse gaussian models. A rigorous proof of monotonicity of updates has been provided for each algorithm. In addition, algorithm-specific measures for quantifying the variation explained by the chosen model have been proposed. Using EMG as well as simulated data, we demonstrated superior performance of these algorithms in delineating muscle synergies by systematically comparing them with existing approaches for signal-dependent noise. It is evident from the methods and results presented that there is a need for more general models for data in which the variance of the signal depends on its mean. It is not entirely surprising that among all algorithms considered, those based on signal-dependent noise clearly outperformed the gaussian model. However, among all algorithms considered based on signal-dependent noise, those based on dual KL and J-divergence showed the best overall performance in terms of both selecting the appropriate model for a given data set and the fraction of variation in the data that was explained by the chosen model. For each data set considered, all three proposed algorithms explained the variation in the data better than existing methods. The variability in the explained variation was also observed to be the smallest for the proposed algorithms. In particular, muscle synergies extracted by J-divergence were the most physiologically interpretable and corroborated with previous findings. The proposed methods therefore provide useful alternatives to current approaches in handling signal-dependent noise and would augment the literature on this topic.

Appendix: Proofs of Theorems

We present detailed proofs of the theorems stated in section 3. In the proof of each theorem, we make use of an auxiliary function similar to the one used in the expectation-maximization (EM) algorithm (Dempster, Laird, & Rubin, 1977; Lee & Seung, 2001). Note that for h real, is an auxiliary function for F(h) if and G(h, h)=F(h) where G and F are scalar valued functions. Also, if G is an auxiliary function, then F is nonincreasing under the update .

A.1.  Proof of Theorem 1.

The cost function, equation 3.12, can be rewritten as
formula
A.1
Its auxiliary function is
formula
A.2
where such that .
Note that . Therefore, and G(Haj, Haj)=F(Haj). The minimizer of F(Haj) is obtained by solving . Using equation A.2, we get
formula
A.3
Solving the above equation results in the update rule for H given in equation 3.13. Similarly, we can rewrite the cost function, equation 3.12, in terms of Wia and obtain the update rule given in equation 3.14.

A.2.  Proof of Theorem 2.

The cost function, equation 3.15, can be rewritten as
formula
A.4
Its auxiliary function is
formula
A.5
where is as defined in the proof of theorem 2.
Note that . Therefore, and G(Haj, Haj)=F(Haj). The minimizer of F(Haj) is obtained by solving . Using equation A.5, we get
formula
A.6
Solving the above equation results in the update rule for H given in equation 3.16. Similarly, we can rewrite the cost function, equation 3.15, in terms of Wia and obtain the update rule given in equation 3.17.

A.3.  Proof of Theorem 3.

The cost function, equation 3.18, can be re-written as
formula
A.7
Its auxiliary function is
formula
A.8
where is as defined in the proof of theorem 2.
Note that . Therefore, and G(Haj, Haj)=F(Haj). The minimizer of F(Haj) is obtained by solving . Using equation A.8, we get
formula
A.9
Solving the above equation results in the update rule for H given in equation 3.19. Similarly, we can rewrite the cost function, equation 3.18 in terms of Wia and obtain the update rule given in equation 3.20.

Acknowledgments

K.D.’s research was supported in part by NIH grant PA CA06297. V.C.K.C. was supported by NIH grants NS44393 and RC1-NS068103-01 to Emilio Bizzi. We thank Emilio Bizzi of MIT for his enthusiastic support of this study.

References

Akaike
,
H.
(
1987
).
Factor analysis and AIC
.
Psychometrika
,
52
(
3
),
317
332
.
Bernstein
,
N.
(
1967
).
The co-ordination and regulation of movements
.
Oxford
:
Pergamon
.
Berry
,
M. W.
,
Browne
,
M.
,
Langville
,
A. N.
,
Pauca
,
V. P.
, &
Plemmons
,
R. J.
(
2007
).
Algorithms and applications for approximate nonnegative matrix factorization
.
Computational Statistics and Data Analysis
,
52
,
155
173
.
Bizzi
,
E.
, &
Cheung
,
V. C. K.
(
2013
).
The neural origin of muscle synergies
.
Frontiers in Computational Neuroscience
,
7
:
51
.
doi:10.3389/fncom.2013.00051
Bizzi
,
E.
,
Cheung
,
V. C. K.
,
d'Avella
,
A.
,
Saltiel
,
P.
, &
Tresch
,
M.
(
2008
).
Combining modules for movement
.
Brain Research Reviews
,
57
,
125
133
.
Cameron
,
A. C.
, &
Windmeijer
,
F. A. G.
(
1996
).
R-squared measures for count data regression models with applications to health care utilization
.
Journal of Business and Economic Statistics
,
14
(
2
),
209
220
.
Cameron
,
A. C.
, &
Windmeijer
,
F. A. G.
(
1997
).
An R-squared measure of goodness of fit for some common nonlinear regression models
.
Journal of Econometrics
,
77
(
2
),
329
342
.
Cheung
,
V. C. K.
,
d'Avella
,
A.
,
Tresch
,
M. C.
, &
Bizzi
,
E.
(
2005
).
Central and sensory contributions to the activation and organization of muscle synergies during natural motor behaviors
.
Journal of Neuroscience
,
25
(
27
),
6419
6434
.
Cheung
,
V. C. K.
, &
Tresch
,
M. C.
(
2005
).
Nonnegative matrix factorization algorithms modeling noise distributions within the exponential family
. In
Proceedings of the 2005 IEEE Engineering in Medicine and Biology 27th Annual Conference
(pp.
4990
4993
).
Piscataway NJ
:
IEEE
.
Cheung
,
V. C. K.
,
Turolla
,
A.
,
Agostini
,
M.
,
Silvoni
,
S.
,
Bennis
,
C.
,
Kasi
,
P.
, …
Bizzi
, E.
(
2012
).
Muscle synergy patterns as physiological markers of motor cortical damage
.
Proceedings of the National Academy of the Sciences USA
,
109
(
36
),
14652
14656
.
Cichocki
,
A.
,
Cruces
,
S.
, &
Amari
,
S.
(
2011
).
Generalized alpha-beta divergences and their application to robust nonnegative matrix factorization
.
Entropy
,
13
,
134
170
.
Cichocki
,
A.
,
Lee
,
H.
,
Kim
,
Y.-D.
, &
Choi
,
S.
(
2008
).
Non-negative matrix factorization with