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.

1 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-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-7].

The expectation-maximization (EM) algorithm has also been applied to solve the sparse state estimate problem in dynamic systems [8-12]. In particular, in [8-10], the EM algorithm is employed to update the parameters of the Bernoulli-Gaussian prior and the measurement noise. These parameters are then used in the generalized approximate message passing algorithm [8-10]. In [12,13], the EM algorithm is used to iteratively estimate the parameters that describe the prior distribution and noise variances. Moreover, in [14], the EM algorithm is used for blind identification, where the sparse state is explored. Note that in the above works, only linear dynamic systems are considered.

In this paper, we focus on the sparse parameter estimation problem instead of the sparse state estimation problem. We consider a general nonlinear dynamic system, where both the state equation and the measurement equation are parameterized by some unknown parameters which are assumed to be sparse. One particular application is the inference of gene regulatory networks. The gene regulatory network can be modeled by the state-space model [15], in which the gene regulations are represented by the unknown parameters. The gene regulatory network is known to be sparse due to the fact that a gene directly regulates or is regularized by a small number of genes [16-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 sparsity-enforced 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.

2 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:

where and are the state vector and the observation vector at time k, respectively; is the unknown parameter vector; (·) and (·) are two nonlinear functions; is the process noise, and is the measurement noise. It is assumed that both { } and { } are independent noise processes and they are mutually independent. Note that the nonlinear functions and are assumed to be differentiable.

Define the notation . The problem considered in this paper is to estimate the unknown system parameter vector from the length-K measurement data . 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, 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-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 and . We assume that the initial state 1 is independent of the parameter . Hence, with the prior given in (3), the optimization in (5) can be rewritten as

where =[λ1,λ2,⋯,λ ] , and ‘ ∘’ denotes the point-wise 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( , | , ′) and p( | , ′), the integrals, and the minimization in (8).

3 The E-step: computing the Q-function

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

3.1 Forward recursion

The forward recursion is composed of two steps: prediction and filtering. Specifically, given the prior probability density function (PDF) p( | ) at time k−1, we need to compute the predicted conditional PDF p( | ); then, given the measurement at time k, we update the filtered PDF p( | ). 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: (1) Given , has a Gaussian distribution, i.e., ; and (2) ( , ) are jointly Gaussian, given .

It then follows that the predictive PDF is Gaussian, i.e., , with [24,26,27]

where )d , and denotes the multivariate Gaussian PDF with mean and covariance .

Moreover, the filtered PDF is also Gaussian, i.e., [24,26,27], where

with

3.2 Backward recursion

In the backward recursion, we compute the smoothed PDFs p( , | ) and p( | ). Here, the approximate assumption made is that conditioned on , and are jointly Gaussian [30], i.e.,

Due to the Markov property of the state-space model, we have p( | , )=p( | , ). Therefore, we can write [30]

Now, assume that

It then follows from (17) and (19) that [30]

where

3.3 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 {( ,w ),i=1,…,N} can be used to approximate the integral

then the general Gaussian-type integral can be approximated by

where = and 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 ( ) 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 ( ) is a highly nonlinear function. The GHQ rule, in contrast, can capture arbitrary degree information of ( ) 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 . The quadrature points and the corresponding weights are given, respectively, by

and

where κ is a tunable parameter, and 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].

4 The M-step: solving the optimization problem

Solving the ℓ1 optimization problems in (8) is not trivial since |θ | is nondifferentiable at θ =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 sparse-simplex 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.

4.1 Iterative thresholding algorithm

Denote as the two summation terms in (8). We consider the optimization problem in (8)

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 the around the solution at the tth iteration is given by

where is the gradient of the negative Q-function and α is such that α mimics the Hessian . Then, is given by

where 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 α mimics the Hessian , i.e., α ≈ , solving α in the least-squares sense, we have

The solution to (32) is given by , where

is the soft thresholding function with sign( ) and max{| |− , } 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.

4.2 Adaptive selection of

So far, the parameters λ 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 ℓ1-minimization problems. λ used for the next iteration are computed from the value of the current solution. A good choice of λ 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 λ =1,∀i. Then, we run the iterative thresholding algorithm to obtain . Next, we update λ as , 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 the maximization step completes. Note that for the iterative re-weighted thresholding algorithm, the assumption that θ has a Laplacian prior no longer holds.

5 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 =[x1,,…,x ] , where x , denotes the gene expression level of the ith gene at time k.

In this case, the nonlinear function ( ) in the general dynamic Equation (1) is given by [15]

with

and

In (41), is an n×n regulatory coefficient matrix with the element a denoting the regulation coefficient from gene j to gene i. A positive coefficient a indicates that gene j activates gene i, and a negative θ indicates that gene j represses gene i. The parameter to be estimated is = which is sparse.

For the measurement model, we have

5.1 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 ( ) is randomly generated from a Gaussian distribution with mean 0 and variance 2. The EM, rEM, and rEM , as well as the basis pursuit de-noising dynamic filtering (BPDN-DF) method and the ℓ1 optimization method, are tested. Here, rEM denotes the version of the rEM algorithm with the iterative re-weighted thresholding discussed in Section 4.2.

As a performance metrics, the receiver operating characteristic (ROC) curve is frequently used. However, for this specific example, with the increasing of the false-positive 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 . 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 using the true parameters (denoted by ) and the support of using the estimated parameters (denoted by ). Note that we assign [ ] =1 if A ≠0 and [ ] =0 if A =0. Similarly, we assign if and if . (2) Compute the number of zero elements of as the matched elements. 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.

5.1.1 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 . 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 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 provides more matched elements than the EM algorithm.

Figure 1

RMSE of rEM with different and rEM .

Figure 2

SF of rEM with different and rEM .

Figure 3

The number of matched elements of rEM with different and rEM .

5.1.2 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 Figures 4, 5, 6. It can be seen that the RMSE of rEM with , =0.01 is smaller than that with , =0.1 . In addition, the rEM with , =0.01 provides a larger number of matched elements than that with , =0.1 as shown in Figure 6. Hence, the estimation accuracy is better when the process noise and measure noise are small.

Figure 4

RMSE of rEM with different noise levels.

Figure 5

SF of rEM with different noise levels.

Figure 6

The number of matched elements of rEM with different noise levels.

5.1.3 The effect of the number of observations

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

Figure 7

RMSE of rEM with different lengths of observations.

Figure 8

SF of rEM with different lengths of observations.

Figure 9

The number of matched elements of rEM with different lengths of observations.

5.1.4 The effect of

In order to test the performance of κ, the rEM 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 with different κ is close. Specifically, it can be seen that the RMSE of rEM using κ=−1 and rEM using κ=−3 is less than that of rEM using κ=0. The sparsity factor of rEM using κ=−1 is more close to 1 than that of rEM using κ=−3 and κ=0. Moreover, the number of matched elements of rEM using κ=−1 is more than that of rEM using κ=−3 and κ=0. Hence, the performance of rEM using κ=−1 is the best in this case.

Figure 10

RMSE of rEM with different .

Figure 11

SR of rEM with different .

Figure 12

The number of matched elements of rEM with different .

5.1.5 Effect of sparsity level

The performance comparison of the rEM and the conventional EM with different sparsity levels of is shown in Figures 13, 14, 15. In this subsection, another which is denser than the previously used is given by

Figure 13

RMSE of the EM and rEM for normal and denser .

Figure 14

SF of the EM and rEM for normal and denser .

Figure 15

The number of matched elements of the EM and rEM for normal and denser .

Note that ‘(Denser)’ is used to denote the result using shown in Equation 48. It can be seen that the RMSE of rEM (Denser) is comparable to that of the EM(Denser). However, the sparsity factor of rEM (Denser) is closer to 1 than that of the EM(Denser) which means that the rEM (Denser) is better. In addition, the number of matched elements of the rEM (Denser) is large than that of the EM(Denser), which means that the rEM (Denser) is better than the EM(Denser) in finding the support of the unknown parameters. The performance of the rEM (Denser), however, is worse than that of the rEM 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.

5.1.6 Comparison with 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:

where and .

We also compare the ℓ1 optimization method with the proposed rEMw 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 optimization-based methods cannot achieve better performance than that of using the rEM .

Figure 16

RMSE of the optimization with different , reweighted optimization, EM, and rEM .

Figure 17

SF of the optimization with different , reweighted optimization, EM, and rEM .

Figure 18

The number of matched elements of the optimization with different , reweighted optimization, EM, and rEM .

5.1.7 Comparison with BPDN-DF

To solve the problem using BPDN-DF, the model in (41) and (44) are modified as

and

respectively. Note that . Then, is given by [36]

where =[λ1,⋯,λ20] with λ =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 , the RMSE of the BPDN-DF is much larger than that of the rEMw. In addition, as shown in Figure 21, the rEMw 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 matrix is an ill-conditioned matrix. In the simulation, λ =0.1, j=5,⋯,20. Based on our tests by using other values of λ, there is no obvious improvement.

Figure 19

RMSE of BPDN-DF and rEM .

Figure 20

SF of BPDN-DF and rEM .

Figure 21

The number of matched elements of BPDN-DF and rEM .

5.2 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 , 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 and covariances =0.01 8 and =0.01 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

where the number of true positives (TP #) denotes the number of links correctly predicted by the inference algorithm, the number of false positives (FP #) denotes the number of incorrectly predicted links. The number of true negatives (TN #) denotes the number of correctly predicted non-links, and the number of false negatives (FN #) denotes the number of missed links by the inference algorithm [15].

The ROC curves of the EM, the rEM, and the rEM 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 performs better than the rEM and the convectional EM algorithms.

Figure 22

ROCs of the EM, the rEM, and the rEM .

In addition, the sparse solution is obtained by using rEM and rEM while it cannot be obtained by using the EM algorithm. The sparsity factor of rEM and rEM is shown in Figure 23; the sparsity of the solution given by rEM is closer to the ground truth than that given by the EM algorithm.

Figure 23

Sparsity factor of the rEM and the rEM .

In Figure 24, the ROC curves of the rEM , ℓ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 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.

Figure 24

ROCs of the rEM , optimization method, and BPDN-DF.

5.3 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 0=0.5 . The process noise and measurement noise are assumed to be Gaussian noise with zero mean and covariance 0.32 and 0.42 , 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 by the EM algorithm is

The inferred by the rEM with λ=1 is

The inferred by the rEM is

The state estimation provided by the UKF based on the model using the inferred parameters of the EM, the rEM, the rEM , 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 is close to the true gene expression data. In addition, The rEM algorithm provide sparser solution than the rEM algorithm. Both the rEM and the rEM algorithms give sparser solutions than the EM algorithm which validates the effectiveness of the proposed method.

Figure 25

True malaria gene expression and estimated gene expression by different algorithms.

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

Competing interests

The authors declare that they have no competing interests.