## 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.

*H*,

_{i}*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*(

*H*) denote the prior probabilities,

_{i}*f*,

*g*the density functions, and

*f*,

*g*the distribution functions corresponding to the hypothesis

*H*,

_{i}*i*=0, 1, respectively. If

*P*(

*H*|

_{i}*x*) denotes the conditional (or posterior) probability of

*H*given

_{i}*X*=

*x*, then, using Bayes’ theorem, we have and Hence, that is, the logarithm of the likelihood ratio, defined as the negative difference between the logarithm of the odds in favor of

*H*

_{0}before and after the observation

*X*=

*x*, is the information in

*X*=

*x*for discrimination in favor of

*H*

_{0}against

*H*

_{1}(Kullback, 1959).

*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

*H*

_{0}against

*H*

_{1}is thus 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

*H*

_{1}against

*H*

_{0}and is given by 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).

*I*(

*f*,

*g*) and

*I*(

*g*,

*f*), one can define J-divergence

*J*(

*f*,

*g*) as which is a measure of the divergence or the difficulty of discriminating between the hypotheses

*H*

_{0}and

*H*

_{1}. 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.

*X*is said to be a member of the family of generalized inverse gaussian (GIG) distributions if its probability density function is given by 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 where and . The mean-variance relationship can be written as and thus indicates a quadratic dependence of variance on the mean. When ,

*f*(

*x*) reduces to an inverse gaussian distribution with density where and . The mean-variance relationship can be written as 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*, *KL ^{d}*, 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.

*f*(

*x*) and

*g*(

*x*), it can be shown that 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*), and 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*), and 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 *h _{ai}* of

*H*is the coefficient of time interval

*i*in synergy

*a*, and the entry

*w*of

_{ja}*W*is the expression level of synergy

*a*in muscle

*j*, where .

*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 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.

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 *R*^{2} 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 *R*^{2} (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.

Measure of Fit . | Behavior . | IP . | VI . | RA . | GA . | RI . |
---|---|---|---|---|---|---|

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 R^{2} | 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 R^{2} | 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 R^{2} | 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 Fit . | Behavior . | IP . | VI . | RA . | GA . | RI . |
---|---|---|---|---|---|---|

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 R^{2} | 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 R^{2} | 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 R^{2} | 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.

*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 -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 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 and MM updates are given by Similarly, when we use equation 2.13, the kernel of

*KL*divergence for the inverse gaussian model can be written as Heuristic updates for

*W*and

*H*are given by and MM updates are given by 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

*KL*,

^{H}_{G}*KL*,

^{MM}_{G}*KL*, and

^{H}_{IG}*KL*to represent the heuristic and MM algorithms for gamma and inverse gaussian models, respectively.

^{MM}_{IG}### 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 *KL ^{d}_{G}*,

*KL*, and

^{d}_{IG}*J*, respectively, to denote these three algorithms presented in theorems 1 to 3. Closed-form multiplicative update rules for

_{G}*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.

See the appendix.

#### 3.2.2. Gamma Model: Algorithm Based on *J* Divergence.

See the appendix.

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

*KL*divergence for the inverse gaussian model can be written as

See the appendix.

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.

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 *V _{ij}* in the denominator of various terms. Theoretically, this should not cause any numerical issues (such as division by zero) since

*V*>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

_{ij}*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.

*r*, the proportion of explained variation (or empirical uncertainty),

*R*

^{2}, is dependent on the particular algorithm and model used in the factorization. For the gaussian NMF algorithm,

*R*

^{2}is the well-known quantity given by 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, *R*^{2} is computed based on the corresponding minimum reconstruction error (*E*), as listed in Table 2. In the *R*^{2} 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 *R*^{2} 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 *R*^{2} in equation 4.1 for the gaussian model to nonlinear models such as the gamma and inverse gaussian. The algorithm-specific *R*^{2} 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).

### 4.2. Akaike Information Criterion (AIC).

*r*,

*AIC*is given by 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 *KL ^{H}_{G}*,

*KL*,

^{MM}_{G}*KL*, and

^{d}_{G}*J*), and three based on inverse gaussian (IG) noise (including

_{G}*KL*,

^{H}_{IG}*KL*, and

^{MM}_{IG}*KL*). 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 (

^{d}_{IG}*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

*R*

^{2}smaller than 10

^{−8}(with

*R*

^{2}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 *H*_{0}:*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.

Frog . | Intact Jump . | Deafferented Jump . | Intact Swim . | Deafferented Swim . |
---|---|---|---|---|

1 | 30,182 | 51,296 | 28,664 | 20,828 |

2 | 13,445 | 12,202 | 47,341 | 35,136 |

3 | 15,820 | 4,638 | 48,033 | 37,997 |

4 | 14,823 | 11,963 | 37,833 | 26,598 |

Frog . | Intact Jump . | Deafferented Jump . | Intact Swim . | Deafferented Swim . |
---|---|---|---|---|

1 | 30,182 | 51,296 | 28,664 | 20,828 |

2 | 13,445 | 12,202 | 47,341 | 35,136 |

3 | 15,820 | 4,638 | 48,033 | 37,997 |

4 | 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 *R*^{2} 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 *R*^{2} 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).

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 (*KL ^{H}_{IG}*,

*KL*, and

^{MM}_{IG}*KL*) consistently indicated a rank of 8 to 13 muscle synergies, gaussian,

^{d}_{IG}*KL*, and

^{d}_{G}*KL*suggested a rank of one to two synergies. The

^{H}_{G}*J*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

_{G}*J*algorithm are the most physiologically interpretable. We further verified that at these ranks, all algorithms returned synergies that explained the EMGs with an

_{G}*R*

^{2}of at least 85% (see Table 4). In addition, all three proposed algorithms returned synergies that explained the EMGs with an

*R*

^{2}of at least 93%, a significantly higher fraction compared to that of existing signal-dependent noise algorithms for which

*R*

^{2}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

*J*algorithm.

_{G}. | Intact Jump . | Deafferented Jump . | Intact Swim . | Deafferented Swim . |
---|---|---|---|---|

Algorithm . | r=3
. | r=3
. | r=4
. | r=4
. |

Gaussian | 87.69 1.21 | 87.38 1.04 | 85.43 1.83 | 87.89 4.51 |

KL ^{H}_{G} | 91.97 1.48 | 90.58 1.17 | 87.17 1.13 | 90.11 2.47 |

KL ^{MM}_{G} | 91.97 1.48 | 90.58 1.17 | 87.16 1.13 | 90.11 2.47 |

KL ^{d}_{G} | 98.93 0.35 | 98.74 0.11 | 96.46 0.41 | 97.62 0.77 |

J _{G} | 97.85 0.64 | 97.41 0.29 | 93.56 0.55 | 95.50 1.38 |

KL ^{H}_{IG} | 90.97 2.06 | 88.72 0.58 | 86.67 0.87 | 89.49 2.30 |

KL ^{MM}_{IG} | 90.89 2.07 | 88.68 0.52 | 86.50 0.83 | 89.43 2.25 |

KL ^{d}_{IG} | 99.79 0.11 | 99.77 0.05 | 98.94 0.30 | 99.37 0.22 |

. | Intact Jump . | Deafferented Jump . | Intact Swim . | Deafferented Swim . |
---|---|---|---|---|

Algorithm . | r=3
. | r=3
. | r=4
. | r=4
. |

Gaussian | 87.69 1.21 | 87.38 1.04 | 85.43 1.83 | 87.89 4.51 |

KL ^{H}_{G} | 91.97 1.48 | 90.58 1.17 | 87.17 1.13 | 90.11 2.47 |

KL ^{MM}_{G} | 91.97 1.48 | 90.58 1.17 | 87.16 1.13 | 90.11 2.47 |

KL ^{d}_{G} | 98.93 0.35 | 98.74 0.11 | 96.46 0.41 | 97.62 0.77 |

J _{G} | 97.85 0.64 | 97.41 0.29 | 93.56 0.55 | 95.50 1.38 |

KL ^{H}_{IG} | 90.97 2.06 | 88.72 0.58 | 86.67 0.87 | 89.49 2.30 |

KL ^{MM}_{IG} | 90.89 2.07 | 88.68 0.52 | 86.50 0.83 | 89.43 2.25 |

KL ^{d}_{IG} | 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 *J _{G}* 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

*J*(mean scalar product = 0.8698; );

_{G}*KL*(0.8695), and

^{d}_{IG}*KL*(0.8650). The worst-performing algorithm was gaussian (0.7648).

^{H}_{G}Table 5 lists the proportion of explained variance (*R*^{2}) 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 *R*^{2} values shown are averages across frogs (*N* = 4; mean SD).

. | . | Rank . | R^{2}
. |
---|---|---|---|

Algorithm . | Behavior . | (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 | |

KL ^{H}_{G} | 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 | |

KL ^{MM}_{G} | 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 | |

KL ^{d}_{G} | 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 | |

J _{G} | 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 | |

KL ^{H}_{IG} | 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 | |

KL ^{MM}_{IG} | 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 | |

KL ^{d}_{IG} | 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 |

. | . | Rank . | R^{2}
. |
---|---|---|---|

Algorithm . | Behavior . | (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 | |

KL ^{H}_{G} | 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 | |

KL ^{MM}_{G} | 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 | |

KL ^{d}_{G} | 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 | |

J _{G} | 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 | |

KL ^{H}_{IG} | 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 | |

KL ^{MM}_{IG} | 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 | |

KL ^{d}_{IG} | 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 (*R*^{2}), both at the ranks with minimum *AIC* and at the ranks determined by the *J _{G}* 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

*J*algorithm provided a much better balance between

_{G}*R*

^{2}and choice of rank based on minimum

*AIC*compared to any other algorithm. The ranks chosen by the

*KL*algorithm based on minimum

^{d}_{G}*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

*KL*algorithm explained the maximum variation (highest overall

^{d}_{IG}*R*

^{2}) among all algorithms, while the gaussian algorithm provided the smallest

*R*

^{2}and rank. Furthermore, the variability of

*R*

^{2}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 *J*_{G} Algorithm Were Physiologically Interpretable.

_{G}

In this section, we compare swim muscle synergies extracted using the standard gaussian algorithm with those identified by the *J _{G}* 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

*J*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

_{G}*J*algorithm, derived from a signal-dependent noise assumption, is better able to discover structures preserved after deafferentation than the traditional gaussian algorithm.

_{G}The muscular compositions of the synergies returned by the *J _{G}* 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

*J*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.

_{G}The activation pattern of synergy 4 identified by *J _{G}* (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 (*R*^{2}), 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 *KL ^{H}_{G}* and

*KL*as reference algorithms, the scalar product values (mean SD; over four frogs) between the synergies returned by gamma-based NMF algorithms and

^{H}_{IG}*KL*were (1) higher than those between the synergies returned by the gaussian NMF algorithm and

^{H}_{G}*KL*and (2) higher than those between inverse gaussian-based NMF algorithms and

^{H}_{G}*KL*(see Figure 6A). Similarly, the scalar product values (mean SD; over four frogs) between the synergies returned by inverse gaussian-based NMF algorithms and

^{H}_{G}*KL*were (1) higher than those between the synergies returned by the gaussian NMF algorithm and

^{H}_{IG}*KL*and (2) higher than those between gamma-based NMF algorithms and

^{H}_{IG}*KL*(see Figure 6B).

^{H}_{IG}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 *R*^{2} 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.

## 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*.

*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 where

*V*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,

_{ij}*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, *KL ^{H}_{G}*,

*KL*,

^{MM}_{G}*KL*,

^{d}_{G}*J*,

_{G}*KL*,

^{H}_{IG}*KL*, and

^{MM}_{IG}*KL*) 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

^{d}_{IG}*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

*R*

^{2}(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 *KL ^{d}_{G}*,

*J*, and all IG-based algorithms in the identification of both

_{G}*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,

*KL*, and

^{H}_{G}*KL*were comparable (see Figure 8A).

^{MM}_{G}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).

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, *).

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 *KL ^{MM}_{G}* 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

*KL*were not only much worse than gaussian but also matched the original quite poorly (see Figure 8B).

^{MM}_{G}## 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 *J _{G}* and

*KL*algorithms while

^{d}_{G}*KL*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.

^{d}_{IG}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.

*G*(

*H*,

_{aj}*H*)=

_{aj}*F*(

*H*). The minimizer of

_{aj}*F*(

*H*) is obtained by solving . Using equation A.2, we get Solving the above equation results in the update rule for

_{aj}*H*given in equation 3.13. Similarly, we can rewrite the cost function, equation 3.12, in terms of

*W*and obtain the update rule given in equation 3.14.

_{ia}### A.2. Proof of Theorem 2.

*G*(

*H*,

_{aj}*H*)=

_{aj}*F*(

*H*). The minimizer of

_{aj}*F*(

*H*) is obtained by solving . Using equation A.5, we get Solving the above equation results in the update rule for

_{aj}*H*given in equation 3.16. Similarly, we can rewrite the cost function, equation 3.15, in terms of

*W*and obtain the update rule given in equation 3.17.

_{ia}### A.3. Proof of Theorem 3.

*G*(

*H*,

_{aj}*H*)=

_{aj}*F*(

*H*). The minimizer of

_{aj}*F*(

*H*) is obtained by solving . Using equation A.8, we get Solving the above equation results in the update rule for

_{aj}*H*given in equation 3.19. Similarly, we can rewrite the cost function, equation 3.18 in terms of

*W*and obtain the update rule given in equation 3.20.

_{ia}## 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.