Regularized EM algorithm for sparse parameter estimation in nonlinear dynamic systems with application to gene regulatory network inference

Parameter estimation in dynamic systems finds applications in various disciplines, including system biology. The well-known expectation-maximization (EM) algorithm is a popular method and has been widely used to solve system identification and parameter estimation problems. However, the conventional EM algorithm cannot exploit the sparsity. On the other hand, in gene regulatory network inference problems, the parameters to be estimated often exhibit sparse structure. In this paper, a regularized expectation-maximization (rEM) algorithm for sparse parameter estimation in nonlinear dynamic systems is proposed that is based on the maximum a posteriori (MAP) estimation and can incorporate the sparse prior. The expectation step involves the forward Gaussian approximation filtering and the backward Gaussian approximation smoothing. The maximization step employs a re-weighted iterative thresholding method. The proposed algorithm is then applied to gene regulatory network inference. Results based on both synthetic and real data show the effectiveness of the proposed algorithm.


Introduction
The dynamic system is a widely used modeling tool that finds applications in many engineering disciplines. Techniques for state estimation in dynamic systems have been well established. Recently, the problem of sparse state estimate has received significant interest. For example, various approaches to static sparse state estimation have been developed in [1][2][3][4], where the problem is essentially an underdetermined inverse problem, i.e., the number of measurements is small compared to the number of states. Extensions of these methods for dynamic sparse state estimation have been addressed in [5][6][7].
The expectation-maximization (EM) algorithm has also been applied to solve the sparse state estimate problem in dynamic systems [8][9][10][11][12]. In particular, in [8][9][10], the http://bsb.eurasipjournals.com/content/2014/1/5 is known to be sparse due to the fact that a gene directly regulates or is regularized by a small number of genes [16][17][18][19]. The EM algorithm has been applied to parameter estimation in dynamic systems [20]. However, the EM algorithm cannot exploit the sparsity of the parameters. Here, we propose a regularized expectation-maximization (rEM) algorithm for sparse parameter estimation in nonlinear dynamic systems. Specifically, the sparsity of the parameters is imposed by a Laplace prior and we consider the approximate maximum a posteriori (MAP) estimate of the parameters. It should be emphasized that the proposed method is an approximate MAP-EM algorithm based on various Gaussian assumptions and quadrature procedures for approximating Gaussian integrals. Note that the MAP-EM algorithm may get stuck at local minima or saddle points. Similar to the conventional EM algorithm, the rEM algorithm also consists of an expectation step and a maximization step. The expectation step involves the forward Gaussian approximation filtering and the backward Gaussian approximation smoothing. The maximization step involves solving an 1 minimization problem for which a re-weighted iterative thresholding algorithm is employed. To illustrate the proposed sparse parameter estimation method in dynamic systems, we consider the gene-regulatory network inference based on gene expression data.
The unscented Kalman filter has been used in the inference of gene regulatory network [15,21,22]. However, the methods proposed in [15,21,22] are fundamentally different with the method proposed in this paper. Firstly, the unscented Kalman filter is only used once in [15,21,22], while it is used in each iteration of the rEM algorithm in this paper. Secondly, not only the unscented Kalman filter but also the unscented Kalman smoother is used in our proposed rEM algorithm. In essence, only the observation before time k is used to the estimation at time k in the unscented Kalman filter. However, in our rEM algorithm, all observation data is used to the estimation at time k (by the unscented Kalman smoother). The fundamental difference between the proposed work and that of [9] is that the proposed work is for the sparse parameter estimation problem of the dynamic system, while that of [9] is only for the sparse parameter estimation of the static problem. In addition, a general nonlinear dynamic system is involved in our work and only linear system is involved in the work of [9]. The main difference between the proposed work and that of [23] is that the sparsity constraint is enforced. The main contribution of this paper is to use the sparsityenforced EM algorithm to solve the sparse parameter estimation problem. In addition, the reweighted iterative threshold algorithm is proposed to solve the 1 optimization algorithm. To the best knowledge of the authors, the proposed rEM with the reweighted iterative threshold optimization algorithm is innovative. Furthermore, we have systematically investigated the performance of the proposed algorithm and compared the results with other conventional algorithms.
The remainder of this paper is organized as follows. In Section 2, the problem of the sparse parameter estimation in dynamic systems is introduced and the regularized EM algorithm is formulated. In Section 3, the E-step of rEM that involves forward-backward recursions and Gaussian approximations is discussed. Section 4 discusses the 1 optimization problem involved in the maximization step. Application of the proposed rEM algorithm to gene regulatory network inference is discussed in Section 5. Concluding remarks are given in Section 6.

Problem statement and the MAP-EM algorithm
We consider a general discrete-time nonlinear dynamic system with unknown parameters, given by the following state and measurement equations: and where x k and y k are the state vector and the observation vector at time k, respectively; θ is the unknown parameter vector; f (·) and h(·) are two nonlinear functions; It is assumed that both {u k } and {v k } are independent noise processes and they are mutually independent. Note that the nonlinear functions f and h are assumed to be differentiable. Define the notation Y k [ y 1 , · · · , y k ]. The problem considered in this paper is to estimate the unknown system parameter vector θ from the length-K measurement data Y K . We assume that θ is sparse. In particular, it has a Laplacian prior distribution which is commonly used as a sparse prior, In the EM algorithm and the MAP-EM algorithm [23], given an estimate θ , a new estimate θ is given by and respectively. Note that the regularized EM can be viewed as a special MAP-EM. To differentiate the sparsity-enforced EM algorithm from the general MAP-EM algorithm, rEM is used. In this paper, the following assumptions are made.
(1) The probability density function of the state is assumed to be Gaussian. The Bayesian filter is optimal; however, http://bsb.eurasipjournals.com/content/2014/1/5 exact finite-dimensional solutions do not exist. Hence, numerical approximation has to be made. The Gaussian approximation is frequently assumed due to the relatively low complexity and high accuracy [24][25][26]. (2) The integrals are approximated by various quadrature methods. Many numerical rules, such as Gauss-Hermite quadrature [25], unscented transformation [27], cubature rule [24], and the sparse grid quadrature [26], as well as the Monte Carlo method [28], can be used to approximate the integral. However, the quadrature rule is the best when computational complexity and accuracy are both considered [29].
We next consider the expression of the Q-function in (5). Due to the Markovian structure of the state-space model (1) to (2), we have Therefore, where c k . We assume that the initial state x 1 is independent of the parameter θ . Hence, with the prior given in (3), the optimization in (5) can be rewritten as where λ = [λ 1 , λ 2 , · · · , λ m ] T , and '•' denotes the pointwise multiplication. Note that in many applications, the unknown parameters θ are only related to the state equation, but not to the measurement equation. Therefore, the second term in (8) can be removed. In the next section, we discuss the procedures for computing the densities p(x k , x k−1 |Y K , θ ) and p(x k |Y K , θ ), the integrals, and the minimization in (8).

The E-step: computing the Q-function
We first discuss the calculation of the probability density functions of the states p(x k , x k−1 |Y K , θ ) and p(x k |Y K , θ ) in (8), which involves a forward recursion of a point-based Gaussian approximation filter to compute p(x k |Y k , θ ) and p(x k+1 |Y k , θ ), k = 1, 2, ..., K , and a backward recursion of a point-based Gaussian approximation smoother to compute p(x k , x k−1 |Y K , θ ) and p(x k |Y K , θ ), k = K , K − 1, ..., 1. For notational simplicity, in the remainder of this section, we drop the parameter θ .

Forward recursion
The forward recursion is composed of two steps: prediction and filtering. Specifically, given the prior probability density function (PDF) p(x k−1 |Y k−1 ) at time k − 1, we need to compute the predicted conditional PDF p(x k |Y k−1 ); then, given the measurement y k at time k, we update the filtered PDF p(x k |Y k ). These PDF recursions are in general computationally intractable unless the system is linear and Gaussian. The Gaussian approximation filters are based on the following two assumptions: It then follows that the predictive PDF is Gaussian, i.e., [24,26,27] where and φ x;x, P denotes the multivariate Gaussian PDF with meanx and covariance P.

Backward recursion
In the backward recursion, we compute the smoothed PDFs p(x k , x k+1 |Y K ) and p(x k |Y K ). Here, the approximate assumption made is that conditioned on y k , x k and x k+1 are jointly Gaussian [30], i.e., with Due to the Markov property of the state-space model, Now, assume that It then follows from (17) and (19) that [30] x k wherẽ

Approximating the integrals
The integrals associated with the expectations in the forward-backward recursions for computing the approximate state PDFs, i.e., (9), (10), (13), (15), (16), and (18), as well as the integrals involved in computing the function Q(θ , θ ) in (8), are integrals of Gaussian type that can be efficiently approximated by various quadrature methods. Specifically, if a set of weighted points (γ i , w i ), i = 1, . . . , N can be used to approximate the integral then the general Gaussian-type integral can be approximated by where P = SS T and S can be obtained by Cholesky decomposition or singular value decomposition (SVD).
By using different point sets, different Gaussian approximation filters and smoothers can be obtained, such as the Gauss-Hermite quadrature (GHQ) [25], the unscented transform (UT) [27], the spherical-radial cubature rule (CR) [24], the sparse grid quadrature rule (SGQ) [26], and the quasi Monte Carlo method (QMC) [28]. Both the UT and the CR are the third-degree numerical rules which means the integration can be exactly calculated when g(x) is a polynomial with the degree up to three. In addition, the form of the CR is identical to the UT with a specific parameter. The main advantage of the UT and the CR is that the number of points required by the rule increases linearly with the dimension. However, one problem of the UT and the CR is that the high-order information of the nonlinear function is difficult to capture so that the accuracy may be low when g(x) is a highly nonlinear function. The GHQ rule, in contrast, can capture arbitrary degree information of g(x) by using more points. It has been proven that GHQ can provide more accurate results than the UT or the CR [25,26]. Similarly, the QMC method can also obtain more accurate results than the UT. However, both the GHQ rule and the QMC method require a large number of points for the high-dimensional problem. Specifically, the number of points required by the GHQ rule increases exponentially with the dimension. To achieve a similar performance of the GHQ with a small number of points, the SGQ is proposed [26], where the number of points increases only polynomially with the dimension.
For the numerical results in this paper, the UT is used in the Gaussian approximation filter and smoother, where we have N = 2n + 1, with n being the dimension of the state vector x k . The quadrature points and the corresponding weights are given, respectively, by http://bsb.eurasipjournals.com/content/2014/1/5 and where κ is a tunable parameter, and e i is the ith n dimensional unit vector. Note that κ = 0 is used as the default value in this paper, as in the cubature Kalman filter [24]. In addition, κ = −3 can also be used as in the unscented Kalman filter [27].

The M-step: solving the 1 optimization problem
Solving the 1 optimization problems in (8) is not trivial since |θ i | is nondifferentiable at θ i = 0. The 1 optimization is a useful tool to obtain sparse solutions. Methods for solving linear inverse problems with sparse constraints are reviewed in [1]. Some more recent developments include the projected scaled subgradient [31] method, the gradient support pursuit method [32], and the greedy sparsesimplex method [33]. In this paper, for the maximization step in the proposed rEM algorithm, due to the simplicity of implementation, we will employ a modified version of the iterative thresholding algorithm.

Iterative thresholding algorithm
DenoteQ(θ , θ ) as the two summation terms in (8). We consider the optimization problem in (8) arg min The solution to (29) can be iteratively obtained by solving a sequence of optimization problems [34]. As in the Newton's method, the Taylor series expansion of thẽ Q(θ , θ ) around the solution θ t at the tth iteration is given byQ where ∇Q is the gradient of the negative Q-function and α t is such that α t I mimics the Hessian ∇ 2Q . Then, θ t+1 is given by where z denotes the variable to be optimized in the objective function.
The equivalent form of (31) is given by Note that Equation 34 is derived as follows. Because we require that α t I mimics the Hessian ∇ 2Q , i.e., α t s t ≈ r t , solving α t in the least-squares sense, we have The solution to (32) is given by is the soft thresholding function with sign(u t ) and max |u t | − a, 0 being component-wise operators. Finally, the iterative procedure for solving (29) is given by And the iteration stops when the following condition is met: where is a given small number.

Adaptive selection of λ
So far, the parameters λ i in the Laplace prior are fixed.
Here, we propose to adaptively tune them based on the output of the iterative thresholding algorithm. The algorithm consists of solving a sequence of weighted 1minimization problems. λ i used for the next iteration are computed from the value of the current solution. A good choice of λ i is to make them counteract the influence of the magnitude of the 1 penalty function [35]. Following this idea, we propose an iterative re-weighted thresholding algorithm. At the beginning of the maximization step, we set λ i = 1, ∀i. Then, we run the iterative thresholding algorithm to obtain θ . Next, we update λ i as λ i = 1 |θ i |+ , ∀i, where is a small positive number, and run the iterative thresholding algorithm again using the new λ. The above process is repeated until it converges at the point where http://bsb.eurasipjournals.com/content/2014/1/5 the maximization step completes. Note that for the iterative re-weighted thresholding algorithm, the assumption that θ has a Laplacian prior no longer holds.

Application to gene regulatory network inference
The gene regulatory network can be described by a graph in which genes are viewed as nodes and edges depict causal relations between genes. By analyzing collected gene expression levels over a period of time, one can find some regulatory relations between different genes. Under the discrete-time state-space modeling, for a gene regulatory network with n genes, the state vector x k = [ x 1,k , . . . , x n,k ] T , where x i,k , denotes the gene expression level of the ith gene at time k.
In this case, the nonlinear function f (x) in the general dynamic Equation (1) is given by [15] f with and In (41), A is an n × n regulatory coefficient matrix with the element a ij denoting the regulation coefficient from gene j to gene i. A positive coefficient a ij indicates that gene j activates gene i, and a negative θ ij indicates that gene j represses gene i. The parameter to be estimated is θ = A which is sparse.
For the measurement model, we have

Inference of gene regulatory network with four genes
In the simulations, we consider a network with four genes. The true gene regulatory coefficients matrix is given by To compare the EM algorithm with the proposed rEM algorithm, the simulation was conduced ten times. In each time, the initial value of A(θ ) is randomly generated from a Gaussian distribution with mean 0 and variance 2. The EM, rEM, and rEM w , as well as the basis pursuit de-noising dynamic filtering (BPDN-DF) method and the As a performance metrics, the receiver operating characteristic (ROC) curve is frequently used. However, for this specific example, with the increasing of the falsepositive rate, the true-positive rate given by rEM and EM is always high (close to 1) which makes the distinguishment of the performance of rEM algorithm and EM algorithm difficult. Hence, the root mean-squared error (RMSE) and the sparsity factor (SF) are used in this section. The RMSE is defined by whereĀ denotes the estimated A. The SF is given by where φ 0 and φ are the number of zero values of the estimated parameter and the number of zero values of true parameter, respectively. It can be seen that the estimation is over sparse if the sparsity factor is greater than 1.
In addition, to test the effectiveness of the proposed method at finding the support of the unknown parameters, the number of matched elements is used and can be obtained by the following procedures: (1) Compute the support of A using the true parameters (denoted by A s ) and the support of A using the estimated parameters (denoted byĀ s ). Note that we assign [ It is easy to see that it is effective at finding the support of the unknown parameters when the number of matched elements is large.

The effect of different λ
The performance of rEM using different λ (10, 5, 1, 0.5, 0.1) is compared with the EM algorithm and the rEM w . The RMSE and SF are shown in Figures 1 and 2, respectively. The RMSE does not increase monotonously with the decreasing of parameter λ. It can be seen that the rEM with λ = 5 has better performance than that using other λ. In addition, the rEM with all λ except λ = 10 outperforms the EM algorithm. It provides smaller RMSE and sparser result. The rEM w provides the smallest RMSE and sparsest parameter estimation. The number of matched elements of test algorithms with different λ is given in Figure 3. It can also be seen that rEM w provides more matched elements than the EM algorithm.

The effect of noise
Two different cases are tested. In the first case, the covariance of the process noise and measurement noise are chosen to be 0.01. In the second case, they are chosen to be 0.1. The performance of two test cases is shown in http://bsb.eurasipjournals.com/content/2014/1/5  It can be seen that the RMSE of rEM w with U, R = 0.01I is smaller than that with U, R = 0.1I. In addition, the rEM w with U, R = 0.01 provides a larger number of matched elements than that with U, R = 0.1 as shown in Figure 6. Hence, the estimation accuracy is better when the process noise and measure noise are small.

The effect of the number of observations
In order to test the effect of the number of observations, the rEM w algorithm with 10 and 20 observations are tested. The simulation results are shown in Figures 7, 8, 9. It can be seen that rEM w with more observations gives less RMSE. In addition, as shown in Figure 9, rEM w with more observations gives slightly better result for finding the support of the unknown parameters.

The effect of κ
In order to test the performance of κ, the rEM w algorithm with different κ (0,-1,-3) is tested. The performance results are shown in Figures 10, 11, 12. Note that the cubature rule corresponds to κ = 0, and the unscented transformation corresponds to κ = −1. Roughly speaking, the performance of rEM w with different κ is close. Specifically, it can be seen that the RMSE of rEM w using κ = −1 and rEM w using κ = −3 is less than that of rEM w using κ = 0. The sparsity factor of rEM w using κ = −1 is more close to 1 than that of rEM w using κ = −3 and κ = 0. Moreover, the number of matched elements of rEM w using κ = −1 is more than that of rEM w using κ = −3 and κ = 0. Hence, the performance of rEM using κ = −1 is the best in this case.

Effect of sparsity level
The performance comparison of the rEM w and the conventional EM with different sparsity levels of A is shown in Figures 13, 14, 15. In this subsection, another A which is denser than the previously used A is given by Note that '(Denser)' is used to denote the result using A shown in Equation 48. It can be seen that the RMSE of rEM w (Denser) is comparable to that of the EM(Denser). However, the sparsity factor of rEM w (Denser) is closer to 1 than that of the EM(Denser) which means that  the rEM w (Denser) is better. In addition, the number of matched elements of the rEM w (Denser) is large than that of the EM(Denser), which means that the rEM w (Denser) is better than the EM(Denser) in finding the support of the unknown parameters. The performance of the rEM w (Denser), however, is worse than that of the rEM w in terms of the improvement of the RMSE. Hence, the rEM algorithm may have close performance with the EM algorithm when the sparsity is not obvious.

Comparison with 1 optimization
We compare the proposed rEM algorithm and the 1 optimization-based method, as well as the conventional EM algorithm. The 1 optimization is a popular approach to obtain the sparse solution. For the problem under consideration, it obtains an estimate of θ by solving the following optimization problem: We also compare the 1 optimization method with the proposed rEM w algorithm, and the results are shown in Figures 16, 17, 18. Seven different λ (5, 2, 1, 0.5, 0.1, 0.05, and 0.01) are used in the 1 optimization method. The RMSE does not decrease monotonously with the decreasing of the parameter λ. Among all tested values, the 1 optimization method with λ = 0.1 gives the smallest RMSE. However, the sparsity factor of the 1 optimization with λ = 0.1 is far from the ideal value 1. The 1 optimization with λ = 5 gives the best support detection as shown  in Figure 18. The re-weighted 1 optimization algorithm is also used in the simulation. However, all 1 optimizationbased methods cannot achieve better performance than that of using the rEM w .

Comparison with BPDN-DF
To solve the problem using BPDN-DF, the model in (41) and (44) are modified as and Then,x k is given by [36] x k = arg miñ where λ = [λ 1 , · · · , λ 20 ] with λ i = 0, i = 1, 2, 3, 4 since our objective is to explore the sparsity of the parameter θ . The exact same initial values used in testing EM and rEM are used to test the performance of the BPDN-DF.
The simulation results are shown in Figures 19, 20, 21. It can be seen that although the sparsity factor of BPDN-DF is comparable with that of the rEM w , the RMSE of the BPDN-DF is much larger than that of the rEM w . In addition, as shown in Figure 21, the rEM w is better than the BPDN-DF in finding the support of the unknown parameters. The possible reason is that the BPDN-DF does not  consider the noise in the dynamic system, and the measurement matrixH k is an ill-conditioned matrix. In the simulation, λ j = 0.1, j = 5, · · · , 20. Based on our tests by using other values of λ, there is no obvious improvement.

Inference of gene regulatory network with eight genes
In this section, we test the proposed algorithm using a larger gene regulatory network which includes eight genes; the performances of the EM, the rEM, the rEM w , the 1 optimization method, and BPDN-DF are given. Forty data points are collected to infer the structure of the network. The system noise and measurement noise are assumed to be Gaussian-distributed with means 0 and covariances U k = 0.01I 8 and R k = 0.01I 8 , respectively. The connection coefficient matrix is given by  For testing, each coefficient inÂ is initialized from a Gaussian distribution with mean 0 and variance 1. The system state is initialized using the first measurement.
The metric used to evaluate the inferred GRN is the ROC curve, in which the true-positive rate (TPR) and the false-positive rate (FPR) are involved. They are given by (FN#) denotes the number of missed links by the inference algorithm [15].
The ROC curves of the EM, the rEM, and the rEM w are compared in Figure 22. The rEM with different λ is tested. In Figure 22, the curves of rEM with four typical values of λ are shown. There is no obvious improvement by using other λ. From the figure, it can be seen that the rEM w performs better than the rEM and the convectional EM algorithms.
In addition, the sparse solution is obtained by using rEM and rEM w while it cannot be obtained by using the EM algorithm. The sparsity factor of rEM and rEM w is shown in Figure 23; the sparsity of the solution given by rEM w is closer to the ground truth than that given by the EM algorithm.  In Figure 24, the ROC curves of the rEM w , 1 optimization method, and BPDN-DF are compared. Similarly, the 1 optimization method with different λ is tested, and only four curves are shown in the figure. By using other values, there is no obvious improvement. The BPDN-DF with different λ has no obvious difference in the test. From Figure 24, it can be seen that the rEM w performs much better than the 1 optimization method and BPDN-DF algorithm. Hence, the sparsity factor of 1 optimization method and BPDN-DF is not shown.

Inference of gene regulatory network from malaria expression data
The dataset with the first six gene expression data of malaria is given in reference [37] and is used in this section. The initial covariance for the algorithm is P 0 = 0.5I. The process noise and measurement noise are assumed to be Gaussian noise with zero mean and covariance 0.3 2 I and 0.4 2 I, respectively. In the following, we show the inference results of the parameter and the state estimation provided by the unscented Kalman filter (UKF) based on the model using the inferred parameters.
The inferred A by the EM algorithm is The inferred A by the rEM with λ = 1 is The state estimation provided by the UKF based on the model using the inferred parameters of the EM, the rEM, the rEM w , and the true gene expression is shown in Figure 25. The left top and right top panels are the expression of the first gene and the second gene, respectively. The left center and right center panels are the expression of the third gene and the fourth gene, respectively. The left bottom and right bottom panels are the expression of the fifth gene and the sixth gene, respectively. It can be seen that the estimate gene expression using the UKF and parameters given by the EM, the rEM, and the rEM w is close to the true gene expression data. In addition, The rEM w algorithm provide sparser   Figure 25 True malaria gene expression and estimated gene expression by different algorithms. http://bsb.eurasipjournals.com/content/2014/1/5 solution than the rEM algorithm. Both the rEM and the rEM w algorithms give sparser solutions than the EM algorithm which validates the effectiveness of the proposed method.

Conclusions
In this paper, we have considered the problem of sparse parameter estimation in a general nonlinear dynamic system, and we have proposed an approximate MAP-EM solution, called the rEM algorithm. The expectation step involves the forward Gaussian approximation filtering and the backward Gaussian approximation smoothing. The maximization step employs a re-weighted iterative thresholding method. We have provided examples of the inference of gene regulatory network based on expression data. Comparisons with the traditional EM algorithm as well as with the existing approach to solving sparse problems such as the 1 optimization and the BPDN-DF show that the proposed rEM algorithm provides both more accurate estimation result and sparser solutions.