We propose a new method for detecting changes in Markov network structure between two sets of samples. Instead of naively fitting two Markov network models separately to the two data sets and figuring out their difference, we directly learn the network structure change by estimating the ratio of Markov network models. This density-ratio formulation naturally allows us to introduce sparsity in the network structure change, which highly contributes to enhancing interpretability. Furthermore, computation of the normalization term, a critical bottleneck of the naive approach, can be remarkably mitigated. We also give the dual formulation of the optimization problem, which further reduces the computation cost for large-scale Markov networks. Through experiments, we demonstrate the usefulness of our method.
Changes in interactions between random variables are interesting in many real-world phenomena. For example, genes may interact with each other in different ways when external stimuli change, co-occurrence between words may appear or disappear when the domains of text corpora shift, and correlation among pixels may change when a surveillance camera captures anomalous activities. Discovering such changes in interactions is a task of great interest in machine learning and data mining because it provides useful insights into underlying mechanisms in many real-world applications.
In this letter, we consider the problem of detecting changes in conditional independence among random variables between two sets of data. Such conditional independence structure can be expressed using an undirected graphical model called a Markov network (MN) (Bishop, 2006; Wainwright & Jordan, 2008; Koller & Friedman, 2009), where nodes and edges represent variables and their conditional dependencies, respectively. As a simple and widely applicable case, the pairwise MN model has been thoroughly studied (Ravikumar, Wainwright, & Lafferty, 2010; Lee, Ganapathi, & Koller, 2007). Following this line, we also focus on the pairwise MN model as a representative example.
A naive approach to change detection in MNs is the two-step procedure of first estimating two MNs separately from two sets of data by maximum likelihood estimation (MLE) and then comparing the structure of the learned MNs. However, MLE is often computationally intractable due to the normalization factor included in the density model. Therefore, gaussianity is often assumed in practice for computing the normalization factor analytically (Hastie, Tibshirani, & Friedman, 2001), though this gaussian assumption is highly restrictive in practice. We may utilize importance sampling (Robert & Casella, 2005) to numerically compute the normalization factor, but an inappropriate choice of the instrumental distribution may lead to an estimate with high variance (Wasserman, 2010); for more discussions on sampling techniques (see Gelman, 1995, and Hinton, 2002). Hyvärinen (2005) and Gutmann and Hyvärinen (2012) have explored an alternative approach to avoid computing the normalization factor that is not based on MLE. However, the two-step procedure has the conceptual weakness that structure change is not directly learned. This indirect nature causes a crucial problem. Suppose that we want to learn a sparse structure change. For learning sparse changes, we may utilize -regularized MLE (Banerjee, El Ghaoui, & d'Aspremont, 2008; Friedman, Hastie, & Tibshirani, 2008; Lee et al., 2007), which produces sparse MNs, and thus the change between MNs also becomes sparse. However, this approach does not work if each MN is dense but only change is sparse.
To mitigate this indirect nature, the fused lasso (Tibshirani, Saunders, Rosset, Zhu, & Knight, 2005) is helpfull, where two MNs are simultaneously learned with a sparsity-inducing penalty on the difference between two MN parameters (Zhang & Wang, 2010; Danaher, Wang, & Witten, 2013). Although this fused-lasso approach allows us to learn sparse structure change naturally, the restrictive gaussian assumption is still necessary to obtain the solution in a computationally tractable way.
The nonparanormal assumption (Liu, Lafferty, & Wasserman, 2009; Liu, Han, Yuan, Lafferty, & Wasserman, 2012) is a useful generalization of the gaussian assumption. A nonparanormal distribution is a semiparametric gaussian copula where each gaussian variable is transformed by a monotone nonlinear function. Nonparanormal distributions are much more flexible than gaussian distributions thanks to the feature-wise nonlinear transformation, while the normalization factors can still be computed analytically. Thus, the fused lasso method combined with nonparanormal models would be one of the state-of-the-art approaches to change detection in MNs. However, the method is still based on separate modeling of two MNs, and its computation for more general nongaussian distributions is challenging.
In this letter, we propose a more direct approach to structural change learning in MNs based on density ratio estimation (DRE) (Sugiyama, Suzuki, & Kanamori, 2012a). Our method does not separately model two MNs, but directly models the change in two MNs. This idea follows Vapnik's principle (Vapnik, 1998): “If you possess a restricted amount of information for solving some problem, try to solve the problem directly and never solve a more general problem as an intermediate step. It is possible that the available information is sufficient for a direct solution but is insufficient for solving a more general intermediate problem.” This principle was used in the development of support vector machines (SVMs): rather than modeling two classes of samples, an SVM directly learns a decision boundary that is sufficient for performing pattern recognition. In the current context, estimating two MNs is more general than detecting changes in MNs (see Figure 1). By directly detecting changes in MNs, we can also halve the number of parameters, from two MNs to one MN difference.
Another important advantage of our DRE-based method is that the normalization factor can be approximated efficiently, because the normalization term in a density ratio function takes the form of the expectation over a data distribution and thus can be simply approximated by the sample average without additional sampling. Through experiments on gene expression and Twitter data analysis, we demonstrate the usefulness of our proposed approach.
The remainder of this letter is structured as follows. In section 2, we formulate the problem of detecting structural changes and review currently available approaches. We propose our DRE-based structural change detection method in section 3. Results of illustrative and real-world experiments are reported in sections 4 and 5, respectively. Finally, we conclude our work and show the future direction in section 6.
2. Problem Formulation and Related Methods
In this section, we formulate the problem of change detection in Markov network structure and review existing approaches.
2.1. Problem Formulation.
Given two densities that can be parameterized using and , our goal is to discover the changes in parameters from P to Q, that is, .
2.2. Sparse Maximum Likelihood Estimation and Graphical Lasso.
Sparse changes in conditional independence structure between P and Q can be detected by comparing two MNs estimated separately using sparse MLE. However, this approach implicitly assumes that two MNs are sparse, which is not necessarily true even if the change is sparse.
2.3. Fused-Lasso Method.
To more naturally handle sparse changes in conditional independence structure between P and Q, a method based on fused-lasso (Tibshirani et al., 2005) has been developed (Zhang & Wang, 2010). This method directly sparsifies the difference between parameters.
Since the Flasso-based method directly sparsifies the change in MN structure, it can work well even when each MN is not sparse. However, using models other than gaussian is difficult because of the normalization issue described in section 2.2.
2.4. Nonparanormal Extensions.
In the above methods, gaussianity is required in practice to compute the normalization factor efficiently, a highly restrictive assumption. To overcome this restriction, it has become popular to perform structure learning under the nonparanormal settings (Liu et al., 2009, 2012), where the gaussian distribution is replaced by a semiparametric gaussian copula.
A random vector is said to follow a nonparanormal distribution if there exists a set of monotone and differentiable functions, , such that follows the gaussian distribution. Nonparanormal distributions are much more flexible than gaussian distributions thanks to the nonlinear transformation , while the normalization factors can still be computed in an analytical way.
However, the nonparanormal transformation is restricted to be element-wise, which is still restrictive to express complex distributions.
2.5. Maximum Likelihood Estimation for Nongaussian Models by Importance Sampling.
A numerical way to obtain the MLE solution under general nongaussian distributions is importance sampling.
However, importance sampling tends to produce an estimate with large variance if the instrumental distribution is not carefully chosen. Although it is often suggested to use a density whose shape is similar to the function to be integrated but with thicker tails as , it is not straightforward in practice to decide which to choose, especially when the dimensionality of x is high (Wasserman, 2010).
3. Direct Learning of Structural Changes via Density Ratio Estimation
The Flasso method can more naturally handle sparse changes in MNs than separate sparse MLE. However, the method is still based on separate modeling of two MNs, and its computation for general high-dimensional nongaussian distributions is challenging. In this section, we propose to directly learn structural changes based on density ratio estimation (Sugiyama et al., 2012a). Our approach does not involve separate modeling of each MN and allows us to approximate the normalization term efficiently for any distributions.
3.1. Density Ratio Formulation for Structural Change Detection.
3.2. Direct Density-Ratio Estimation.
Density ratio estimation has been recently introduced to the machine learning community and has proven to be useful in a wide range of applications (Sugiyama et al., 2012a). Here, we concentrate on the density ratio estimator called the Kullback-Leibler importance estimation procedure (KLIEP) for log-linear models (Sugiyama et al., 2008; Tsuboi, Kashima, Hido, Bickel, & Sugiyama, 2009).
3.3. Sparsity-Inducing Norm.
To find a sparse change between P and Q, we propose regularizing the KLIEP solution with a sparsity-inducing norm . Note that the MLE approach sparsifies both and so that the difference is also sparsified, while we directly sparsify the difference ; thus our method can still work well even if and are dense.
3.4. Dual Formulation for High-Dimensional Data.
The solution of the optimization problem, equation 3.3, can be easily obtained by standard sparse optimization methods. However, in the case where the input dimensionality d is high (often the case in our setup), the dimensionality of parameter vector is large, and thus obtaining the solution can be computationally expensive. Here, we derive a dual optimization problem (Boyd & Vandenberghe, 2004), which can be solved more efficiently for high-dimensional (see Figure 2).
Note that the dimensionality of the dual variable is equal to nQ, while that of is quadratic with respect to the input dimensionality d because we are considering pairwise factors. Thus, if d is not small and nQ is not very large (often the case in our experiments shown later), solving the dual optimization problem would be computationally more efficient. Furthermore, the dual objective (and its gradient) can be computed efficiently in parallel for each (u, v), a useful property when handling large-scale MNs. Note that the dual objective is differentiable everywhere, while the primal objective is not.
4. Numerical Experiments
In this section, we compare the performance of the proposed KLIEP-based method, the Flasso method, and the Glasso method for gaussian models, nonparanormal models, and nongaussian models. Results are reported on data sets with three different underlying distributions: multivariate gaussian, nonparanormal, and nongaussian “diamond” distributions. We also investigate the computation time of the primal and dual formulations as a function of the input dimensionality. The Matlab implementation of the primal and dual methods are available online at http://sugiyama-www.cs.titech.ac.jp/~song/SCD.html.
4.1. Gaussian Distribution.
We compare the performance of the KLIEP, Flasso, and Glasso methods. Because all methods use the same gaussian model, the difference in performance is caused only by the difference in estimation methods. We repeat the experiments 20 times with randomly generated data sets and report the results in Figure 3.
The top six graphs are examples of regularization paths.6 The dashed lines represent changed edges in the ground truth, while the solid lines represent unchanged edges. The top row is for n = 100, while the middle row is for n = 50. The bottom three graphs are the data-generating distribution and averaged precision-recall (P-R) curves with standard error over 20 runs. The P-R curves are plotted by varying the group-sparsity control parameter with in KLIEP and Flasso, and by varying the sparsity control parameters as in Glasso.
In the regularization path plots, solid vertical lines show the regularization parameter values picked based on holdout data and as follows:
When n=100, KLIEP and Flasso clearly distinguish changed (dashed lines) and unchanged (solid lines) edges in terms of parameter magnitude. However, when the sample size is halved to n=50, the separation is visually rather unclear in the case of Flasso. In contrast, the paths of changed and unchanged edges are still almost disjoint in the case of KLIEP. The Glasso method performs rather poorly in both cases. A similar tendency can also be observed in the P-R curve plot: when the sample size is n=100, KLIEP and Flasso work equally well, but KLIEP gains its lead when the sample size is reduced to n=50. Glasso does not perform well in both cases.
4.2. Nonparanormal Distribution.
The experiments are conducted on 20 randomly generated data sets with n = 50 and 100, respectively. The regularization paths, data-generating distribution, and averaged P-R curves are plotted in Figure 4. The results show that Flasso clearly suffers from performance degradation compared with the gaussian case, perhaps because the number of samples is too small for the complicated nonparanormal distribution. Due to the two-step estimation scheme, the performance of Glasso is poor. In contrast, KLIEP separates changed and unchanged edges clearly for both n = 50 and n = 100. The P-R curves also show the same tendency.
4.3. “Diamond” Distribution with No Pearson Correlation.
In the experiments in section 4.2, though samples are nongaussian, the Pearson correlation is not zero. Therefore, methods assuming gaussianity can still capture some linear correlation between random variables. Here, we consider a more challenging case with a diamond-shaped distribution within the exponential family that has zero Pearson correlation between variables. Thus, the methods assuming gaussianity cannot extract any information in principle from this data set.
We set d = 9 and nP=nQ=5000. AP is randomly generated with 35% sparsity, while AQ is created by randomly removing edges in AP so that the sparsity level is dropped to 15%. Samples from the above distribution are drawn by using a slice sampling method (Neal, 2003). Since generating samples from high-dimensional distributions is nontrivial and time-consuming, we focus on a relatively low-dimensional case. To avoid sampling error, which may mislead the experimental evaluation, we also increase the sample size, so that the erratic points generated by accident will not affect the overall population.
The averaged P-R curves over 20 data sets are shown in Figure 5e. KLIEP with the polynomial model significantly outperforms all the other methods, while the IS-Glasso and especially IS-Flasso give better results than the KLIEP, Flasso, and Glasso methods with the gaussian and nonparanormal models. This means that the polynomial basis function is indeed helpful in handling completely nongaussian data. However, as discussed in section 2.2, it is difficult to use such a basis function in Glasso and Flasso because of the computational intractability of the normalization term. Although IS-Glasso can approximate integrals, the result shows that such approximation of integrals does not lead to very good performance. In comparison, the result of the IS-Flasso method is much improved thanks to the coupled sparsity regularization, but it is still not comparable to KLIEP.
The regularization paths of KLIEP with the polynomial model illustrated in Figure 5b show the usefulness of the proposed method in change detection under nongaussianity. We also give regularization paths obtained by the IS-Flasso and IS-Glasso methods on the same data set in Figures 5c and 5d, respectively. The graphs show that both methods do not separate changed and unchanged edges well, though the IS-Flasso method works slightly better.
4.4. Computation Time: Dual Versus Primal Optimization Problems.
Finally, we compare the computation time of the proposed KLIEP method when solving the dual optimization problem, equation 3.4, and the primal optimization problem, equation 3.3. Both the optimization problems are solved by using the same convex optimizer minFunc.7 The data sets are generated from two gaussian distributions constructed in the same way as in section 4.1. we draw 150 samples separately from two distributions with dimension d = 40, 50, 60, 70, 80. We then perform change detection by computing the regularization paths using 20 choices of ranging from 10−4 to 100 and fix . The results are plotted in Figure 6.
It can be seen from the graph that as the dimensionality increases, the computation time for solving the primal optimization problem is sharply increased, while that for solving the dual optimization problem grows only moderately: when d=80, the computation time for obtaining the primal solution is almost 10 times more than that required for obtaining the dual solution. Thus, the dual formulation is computationally much more efficient than the primal formulation.
In this section, we report the experimental results on a synthetic gene expression data set and a Twitter data set.
5.1. Synthetic Gene Expression Data Set.
A gene regulatory network encodes interactions between DNA segments. However, the way genes interact may change due to environmental or biological stimuli. In this experiment, we focus on detecting such changes. We use SynTReN, a generator of gene regulatory networks used for benchmark validation of bioinformatics algorithms (Van den Bulcke et al., 2006).
We first choose a subnetwork containing 13 nodes from an existing signaling network in Saccharomyces cerevisiae (shown in Figure 7a). Three types of interactions are modeled: activation (ac), deactivation (re), and dual (du). Fifty samples are generated in the first stage, after which we change the types of interactions in six edges and generate 50 samples again. Four types of changes are considered: ac re, re ac, du ac, and du re.
We use KLIEP and IS-Flasso with the polynomial transform function for . The regularization parameter in KLIEP and Flasso is tested with choices . We set the instrumental distribution as the standard normal and use sample for approximating integrals in IS-Flasso.
The regularization paths on one example data set for KLIEP, IS-Flasso, and the plain Flasso with the gaussian model are plotted in Figures 7b, 7c, and 7d, respectively. Averaged P-R curves over 20 simulation runs are shown in Figure 7e. We can see clearly from the KLIEP regularization paths shown in Figure 7b that the magnitude of estimated parameters on the changed pairwise interactions is much higher than that of the unchanged edges. IS-Flasso also achieves rather clear separation between changed and unchanged interactions, though a few unchanged interactions drop to zero at the final stage. Flasso gives many false alarms by assigning nonzero values to the unchanged edges, even after some changed edges hit zeros.
Reflecting a similar pattern, the P-R curves plotted in Figure 7e show that the proposed KLIEP method has the best performance among all three methods. We can also see that the IS-Flasso method achieves significant improvement over the plain Flasso method with the gaussian model. The improvement from Flasso to IS-Flasso shows that the use of the polynomial basis is helpful on this data set, and the improvement from IS-Flasso to KLIEP shows that the direct estimation can further boost the performance.
5.2. Twitter Storytelling.
Finally, we use KLIEP with the polynomial transform function for and Flasso as event detectors from Twitter. More specifically, we choose the Deepwater Horizon oil spill as the target event and hope that our method can recover some story lines from Twitter as the news events develop.8 Counting the frequencies of 10 keywords (BP, oil, spill, Mexico, gulf, coast, Hayward, Halliburton, Transocean, and Obama), we obtain a data set by sampling four times per day from February 1, 2010, to October 15, 2010, resulting in 1061 data samples.
We segment the data into two parts: the first 300 samples collected before the day of the oil spill (April 20, 2010) are regarded as conforming to a 10-dimensional joint distribution Q, while the second set of samples that are in an arbitrary 50-day window after the oil spill accident happened is regarded as following distribution P. Thus, the MN of Q encodes the original conditional independence of frequencies between 10 keywords, while the underlying MN of P has changed since an event occurred. We expect that unveiling changes in MNs between P and Q can recover the drift of popular topic trends on Twitter in terms of the dependency among keywords.
The detected change graphs (i.e., the graphs with only detected changing edges) on 10 keywords are illustrated in Figure 8. The edges are selected at a certain value of indicated by the maximal cross-validated log likelihood (CVLL). Since the edge set that is picked by CVLL may not be sparse in general, we sparsify the graph based on the permutation test as follows. We randomly shuffle the samples between P and Q and repeatedly run change detection algorithms 100 times; then we observe detected edges by CVLL. Finally, we select the edges that are detected using the original nonshuffled data set and remove those that were detected in the shuffled data sets more than five times (i.e., the significance level 5%). For KLIEP, k is also tuned by using CVLL. In Figure 8, we plot detected change graphs generated using samples of P starting from April 17, July 6, and July 26, respectively.
The initial explosion happened on April 20, 2010. Both methods discover dependency changes between keywords. Generally, KLIEP captures more conditional independence changes between keywords than the Flasso method, especially when comparing Figures 8c and 8f. At the first two stages (see Figures 8a, 8b, 8d, and 8e) the keyword Obama is very well connected with other keywords in the results given by both methods. Indeed, at the early development of this event, he the president lies in the center of the news stories, and his media exposure peaks after his visit to the Louisiana coast (May 2, May 28, and June 5) and his meeting with BP CEO Tony Hayward on June 16. Notably, both methods highlight the “gulf-obama-coast” triangle in Figures 8a and 8d and the “bp-obama-hayward” chain in Figures 8b and 8e.
However, there are some important differences worth mentioning. First, the Flasso method misses the “transocean-hayward-obama” triangle in Figures 8d and 8e. Transocean is the contracted operator in the Deepwater Horizon platform, where the initial explosion happened. On Figure 8c, the chain “bp-spill-oil” may indicate that the phrase “bp spill” or “oil spill” has been publicly recognized by the Twitter community since then, while the “hayward-bp-mexico” triangle, although relatively weak, may link to the event that Hayward stepped down as CEO on July 27.
It is also noted that Flasso cannot find any changed edges in Figure 8f, perhaps due to the gaussian restriction.
6. Discussion, Conclusion, and Future Work
In this letter, we proposed a direct approach to learning sparse changes in MNs by density ratio estimation. Rather than fitting two MNs separately to data and comparing them to detect a change, we estimated the ratio of the probability densities of two MNs where changes can be naturally encoded as sparsity patterns in estimated parameters. This direct modeling allows us to halve the number of parameters and approximate the normalization term in the density ratio model by a sample average without sampling. We also showed that the number of parameters to be optimized can be further reduced with the dual formulation, which is highly useful when the dimensionality is high. Through experiments on artificial and real-world data sets, we demonstrated the usefulness of the proposed method over state-of-the-art methods, including nonparanormal-based methods, and sampling-based methods.
Our important future work is to theoretically elucidate the advantage of the proposed method, beyond Vapnik's principle of solving the target problem directly. The relation to score matching (Hyvärinen, 2005), which avoids computing the normalization term in density estimation, is also an interesting issue to be further investigated. Considering higher-order MN models such as the hierarchical log-linear model (Schmidt & Murphy, 2010) is a promising direction for extension.
In the context of change detection, we are mainly interested in the situation where p and q are close to each other (if p and q are completely different, it is straightforward to detect changes). When p and q are similar, density ratio estimation for p(x)/q(x) or q(x)/p(x) performs similarly. However, given the asymmetry of density ratios, the solutions for p(x)/q(x) or q(x)/p(x) are generally different. The choice of the numerator and denominator in the ratio is left for future investigation.
Detecting changes in MNs is the main target of this letter. Estimating the difference or divergence between two probability distributions has been studied under a more general context in the statistics and machine learning communities (Amari & Nagaoka, 2000; Eguchi & Copas, 2006; Wang, Kulkarni, & Verdú, 2009; Sugiyama, Suzuki, & Kanamori, 2012b; Sugiyama, Liu, et al., 2013). In fact, the estimation of the Kullback-Leibler divergence (Kullback & Leibler, 1951) is related to the KLIEP-type density ratio estimation method (Nguyen, Wainwright, & Jordan, 2010), and the estimation of the Pearson divergence (Pearson, 1900) is related to the squared-loss density ratio estimation method (Kanamori, Hido, & Sugiyama, 2009). However, the density-ratio-based divergences tend to be sensitive to outliers. To overcome this problem, a divergence measure based on relative density ratios was introduced, and its direct estimation method was developed (Yamada, Suzuki, Kanamori, Hachiya, & Sugiyama, 2013). L2-distance is another popular difference measure between probability density functions. L2-distance is symmetric, unlike the Kullback-Leibler divergence and the Pearson divergence, and its direct estimation method has been investigated recently (Sugiyama, Suzuki, et al., 2013; Kim & Scott, 2010).
Change detection in time series is a related topic. A straightforward approach is to evaluate the difference (dissimilarity) between two consecutive segments of time-series signals. Various methods have been developed to identify the difference by fitting two models to two segments of time series separately, for example, the singular spectrum transform (Moskvina & Zhigljavsky, 2003; Ide & Tsuda, 2007), subspace identification (Kawahara, Yairi, & Machida, 2007), and the method based on the one-class support vector machine (Desobry, Davy, & Doncarli, 2005). In the same way as this letter, direct modeling of the change has also been explored for change detection in time-series (Kawahara & Sugiyama, 2012; Liu, Yamada, Collier, & Sugiyama, 2013; Sugiyama, Suzuki, et al., 2013).
Appendix: Derivation of the Dual Optimization Problem
S.L. is supported by the JST PRESTO program and the JSPS fellowship. J.Q. is supported by the JST PRESTO program. M.U.G. is supported by the Finnish Centre-of-Excellence in Computational Inference Research COIN (251170). T.S. is partially supported by MEXT Kakenhi 25730013, and the Aihara Project, the FIRST program from JSPS, initiated by CSTP. M.S. is supported by the JST CREST program and AOARD.
Note that the proposed algorithm itself can be applied to any MNs containing more than two elements in each factor.
From here on, we simplify as .
For implementation simplicity, we maximize the joint likelihood of p and q instead of its feature-wise conditional likelihood. We also switch the first penalty term from to .
We set for not breaking the symmetry of the precision matrix.
Paths of univariate factors are omitted for clear visibility.
An earlier version of this work was presented at the European Conference on Machine Learning and Principles, and Practice of Knowledge Discovery in Databases (ECML/PKDD2013), September 23–27, 2013.