GRASP is a three-stage algorithm for learning the structure of a BN. In the first (pruning) stage, we designed a Double Filtering (DF) method to find the cover of the skeleton of the BN, where the skeleton of a BN is defined as the BN structure after removing the direction of all the edges, and the cover is defined as a superset of undirected edges containing all the edges of the skeleton. In the second (structure searching) stage, we developed an adaptive sequential Monte Carlo (adSMC) method to search the BN structure on the undirected network found in the first stage based on Bayesian information criterion (BIC) score. To reclaim the potentially missed edges, we designed a Random Order Hill Climbing (ROHC) method as the third stage.

The first stage, namely Double Filtering (DF) method, contains two filtering processes. The first filtering was done by unconditioned tests, filtering out the nodes that are not ancestors or descendants of a given node X i . The second filtering was built on conditioned tests, further filtering out the nodes that are not direct neighbors (parents or children) of X i .

Suppose we have p nodes. For a given node X i , let nbr ( X i ) be the set of nodes that have an undirected edge with X i . Initialize nbr ( X i ) = ∅ the empty set for all i . The procedure of the DF method is as follows:

1. First filtering. For each pair of nodes ( X i and X j , i ≠ j ), record the p-value of their Mutual Information (MI) test as p ij . If p ij < α , α being a predetermined significance level for the test, add X j into n b r ( X i )   and X i into n b r ( X j ) . The MI test is also used by other BN structure learning methods (Campos 1998). Let N i be the number of elements in nbr ( X i ) after first filtering. Sort nbr ( X i ) using their p -values in ascending order and denote these nodes as X 1 ( i ) ,   X 2 ( i ) ,   … , X N i ( i ) .

After applying the DF method on all p nodes, the collection of nbr C ( X i ) ,   i = 1,2 , ... , p gives us the skeleton of BN.

In the pruned space, we designed an adaptive sequential Monte Carlo (adSMC) method to search the structure of the Bayesian network. In a traditional sequential Monte Carlo, the random variable of all p nodes (or features) X ∈ ℝ p is decomposed into M blocks ( x 1 ,   x 2 ,   … ,   x M ) each with p i features x i ∈ ℝ p i and ∑ i = 1 M p i = p , and the decomposition is predefined and fixed throughout the whole sampling procedure. One usually samples   x 1 at first, then x 2 , and so on. However, the sequence each variable is sampled (namely sampling sequence in this study) based on any prior decomposition may not be the most efficient one. The optimal sampling sequence may need to be decided dynamically. For example, when x 1 ,   x 2 , … ,   x m − 1 have been sampled for some 0 < m ≤ M , the conditional distribution f (   x m | x 1 ,   x 2 , … ,   x m − 1 ) may have a small set of candidate decompositions (to satisfy the acyclic condition) which limits the diversity of the SMC samples. Therefore, we designed our sampling block x m conditioning on the current sampled structure x 1 ,   x 2 , … ,   x m − 1 to increase the diversity and quality of obtained samples (see Supplementary Figure S1 in Supplementary Materials for an example).

For each SMC sample, we start with all possible fully connected triplets (three nodes connected by three undirected edges) discovered in the first stage. We sample one such triplet having the least outside connection, e.g., the one having the least undirected edges connected to its nodes (Figure 1A). These triplets are likely to be restricted to certain configuration by the sampled structure; therefore, to sample them earlier allows more variety in their configurations. When all fully connected triplets are sampled, partially connected triplets (two undirected edges among three nodes) are considered (Figure 1B). Lastly, we consider pairs (the remaining undirected edges, Figure 1C). For partially connected triplets and pairs, the configurations with the least outside connections are sampled first.

The probabilities of possible configurations of triplets and pairs are proportional to their BIC (Bayesian Information Criterion) score defined as BIC = ln ( n ) ⋅ k − 2 ⋅ ln ( L ^ ) , where   n represents sample size, k represents number of parameters of the partial BN model, and L ^ = p ( x | θ ^ , Model ) is the maximized value of the likelihood function of the model, with estimated parameters θ ^ from observed data   x . The probability is set to be 0 if certain configuration fails to satisfy the acyclicity condition. In summary, for triplets (pairs), we calculate the BIC for all possible configurations of this triplet (pair) and sample one configuration with probability P ( configuration   i ) ∝ { exp ( BIC i T ) ,  if    configuration   i   does   not   result   in   a   loop 0 , if    configuration   i   result   in   a   loop , (1) Where exp ( ⋅ ) is the exponential function, T is temperature controlling how greedy we want the searching to be, and ∝ means proportional to.

The main algorithm used in the second step is as follows.

Step 1

Sample one fully connected triplet { X I , X J , X K } with the least outside connection; Choose a configuration between these three nodes with probability described in Eq. 1; Then remove the connection between X I , X J and X K from the skeleton.

Step 2

Repeat step 1 until all fully connected triplets are sampled.

Step 3

Sample one partially connected triplet { X I , X J , X K } with the least outside connection. Choose a configuration between these three nodes with probability described in Eq. 1. Then remove the connection between X I , X J and X K (if applicable) from the skeleton.

Step 4

Repeat step 3 until all partially connected triplets are sampled.

Step 5

Sample one pair { X I , X J } with the least outside connection. Choose a configuration between them (either   X I → X J   or   X I ← X J ) with probability described in Eq. 1. Then remove the connection between X I and X J from the skeleton.

Step 6

Repeat step 5 until all pairs are sampled and no more unsampled edges in the skeleton.

Since each SMC sample is generated independently, we can run our algorithm in parallel on multiple CPUs/GPU cores to speed up the sampling process.

We mentioned earlier that one disadvantage of the traditional two-stage method was that the edges missed in the first stage will never be recovered. Therefore, in the third stage we designed a Random Order Hill Climbing (ROHC) method to identify the possible missed edges and refine the network. The general idea is described as follows:

1. Generate a permutation of 1,2 , … ,   p for each network sampled by adSMC in stage 2, suppose m 1 , … ,   m p is such a permutation.

One could also view this stage as a further ascent to the local optima to ensure we achieve the best possible BIC score.

In `general, generating more SMC samples gives a higher chance to reach the optimum. However, more samples also require more computation time; therefore, a balance between running time and sample sizes must be made. In most of our simulation study and practical problems, we found that around 20,000 samples were often good enough for finding a network with a satisfactory BIC score.

To measure the effectiveness of edge screening methods, we employed the precision, recall and f-score measurements. Precision is defined as TP/(TP + FP), recall is defined as TP/(TP + FN), and f-score is the harmonic mean of precision and recall, 2 (precision × × recall)/(precision + recall), where TP means true positive (number of true undirected edges identified), FP false positive (number of non-edges identified as undirected edges), and FN false negative (number of undirected edges not identified).

In our study, recall measures the percentage of true edges (irrespective of their directions) identified; therefore, it is the most important measurement in edge screening stage, since as we discussed earlier, any missed edges in stage one may never be reclaimed in a traditional two stage approach. Besides the recall, f-score is also important since it measures a balanced performance in terms of both precision and recall. It is obvious that if we propose all possible edges, we will always identify all true edges, but that will not do any pruning to the searching space. Thus, a high f-score is desired for a decent edge screening strategy.

We used Bayesian Information Criterion (BIC) as the score function in both second stage and third stage. BIC has the score-equivalent property (Appendix definition 10), which can reduce the searching space, since if we could find one network in the equivalent class, we found the true network. And the consistency property of BIC score guarantees that the true network has the highest score asymptotically.

The networks used to generate simulated data (Table 1) are from actual decision making processes of a wide range of real applications including risk management, tech support, and disease diagnosis. All networks are obtained from Bayesian Network Repository maintained by M. Scutari http://www.bnlearn.com/bnrepository/.

We randomly generated data with 1,000, 2,000, and 5,000 observations, and we generated 10 datasets for each size of observations. All results reported in this section are based on averages of 10 datasets. Observation size in this article refers to the number of data points, and shall not be confused with number of sequential Monte Carlo samples. The datasets were generated using R package bnlearn (Scutari 2009; Nagarajan, Scutari and Lèbre 2013).

The principal of the edge screening stage is pruning the searching space as much as possible while the remaining edges in the pruned space still possess as many true edges as possible. We compare our method to five other methods including max-min parent-child (mmpc) (Tsamardinos, et al., 2006), grow-shrink (gs) (Margaritis 2003), incremental association (iamb) (Tsamardinos, et al., 2003), fast iamb, and inter iamb (Yaramakala, et al., 2005). For all methods, we fixed the significance level (α) to 0.01.

The simulation study results (Figure 2 and Figure S2) showed that our double filtering (DF) method was able to identify the most edges (highest recall) for each of the observation size we tested. In some cases we observed that with even 1,000 observations, our method achieved a higher recall than the other methods using 5,000 observations and the f-scores are still comparable (e.g., Alarm, Hepar2 and etc.). For some networks (Child, Insurance), not only the recalls were higher but also the f-scores were higher for DF. The results confirmed that DF identifies true edges more accurately than other methods and it often requires fewer observations. Higher recall is desired in the first stage (the edge screening stage) since any missed edges will not be sampled in the second stage.