Uncovering Hidden Structure through Parallel Problem Decomposition for the Set Basis Problem: Application to Materials Discovery

1 Introduction

Exploiting parallelism and multi-core architectures is a natural way to speed up computations in many domains. Recently, there has been great success in parallel computation in fields such as scientific computing and information retrieval [Dean and Ghemawat, 2008; Chu et al., 2007].

Parallelism has also been taken into account as a promising way to solve hard combinatorial problems. However, it remains challenging to exploit parallelism to speed up combinatorial search because of the intricate non-local nature of the interactions between variables in hard problems [Hamadi and Wintersteiger, 2013]. One class of approaches in this domain is divide-and-conquer, which dynamically splits the Carla P. Gomes Cornell University, USA gomes@cs.cornell.edu Bart Selman Cornell University, USA selman@cs.cornell.edu search space into sub-spaces, and allocates each sub-space to a parallel node [Chrabakh and Wolski, 2003; Chu et al., 2008; Rao and Kumar, 1993; Regin et al., 2013; Moisan et al., 2013; Fischetti et al., 2014]. A key challenge in this approach is that the solution time for subproblems can vary by several orders of magnitude and is highly unpredictable. Frequent load re-balancing is required to keep all processors busy, but the load re-balancing process can result in a substantial overhead cost. Another class of approaches harnesses portfolio strategies, which runs a portfolio of solvers (of different type or with different randomization) in parallel, and terminates as soon as one of the algorithms completes. [Xu et al., 2008; Leyton-Brown et al., 2003; Malitsky et al., 2011; Kadioglu et al., 2011; Hamadi and Sais, 2009; Biere, 2010; Kottler and Kaufmann, 2011; Schubert et al., 2010; O'Mahony et al., 2008]. Parallel portfolio approaches can be highly effective. They do require however the use of a collection of effective solvers that each excel at different types of problem instances. In certain areas, such as SAT/SMT solving, we have such collections of solvers but for other combinatorial tasks, we do not have many different solvers available.

In this paper, we exploit parallelism to boost combinatorial search in a novel way. Our framework complements the two parallel approaches discussed before. In our approach parallelism is used as a preprossessing step to identify a promising portion of the search space to be explored by a complete sequential solver. In our scheme, a set of parallel processes are first deployed to solve a series of related subproblems. Next, the solutions to these subproblems are aggregated to obtain an initial guess for a candidate solution to the original problem. The aggregation is based on a key empirical observation that solutions to the subproblems, when properly aggregated, provide information about solutions for the original problem. Lastly, a global sequential solver searches for a solution in an iterative deepening manner, starting from the promising portion of the search space identified by the previous aggregation step. At a high level, the initial guess obtained by aggregating solutions to subproblems provides the so-called backdoor information to the sequential solver, by forcing it to start from the most promising portion of the search space. A backdoor set is a set of variables, such that once their values are set correctly, the remaining problem can be solved in polynomial time [Williams et al., 2003; Dilkina et al., 2009; Hamadi et al., 2011].

We empirically show that a global solver, when initialized with proper information obtained by solving the subproblems, can solve a set of instances in seconds, while it takes for the same solver hours to days to find the solution without initialization. The strategy also outperforms state-ofthe-art incomplete solvers in terms of solution quality.

We apply the parallel scheme to an NP-complete problem called the Set Basis Problem, in which we are given a collection of subsets C of a finite set U. The task is to find another, hopefully smaller, collection of subsets of U, called a "set basis", such that each subset in C can be represented exactly by a union of sets from the set basis. Intuitively, the set basis provides a compact representation of the original collection of sets. The set basis problem occurs in a range of applications, most prominently in machine learning, e.g., used as a special type of matrix factorization technique [Miettinen et al., 2008]. It also has applications in system security and protection, where it is referred to as the role mining problem in access control [Vaidya et al., 2007]. It also has applications in secure broadcasting [Shu et al., 2006] and computational biology [Nau et al., 1978].

While having many natural applications, our work is motivated by a novel application in the field of computational sustainability [Gomes, Winter 2009], concerning the discovery of new materials for renewable energy sources such as improved fuel cell catalysts [Le Bras et al., 2011]. In this domain, the set basis problem is used to find a succinct explanation of a large set of measurements (X-ray diffraction patterns) that are represented in a discrete way as sets. Mathematically, this corresponds to a generalized version of the set basis problem with extra constraints. Our parallel solver can be applied to this generalized version as well, and we demonstrate significant speedups on a set of challenging benchmarks. Our work opens up a novel angle for using parallelism to solve a set of hard combinatorial problems.

2 Set Basis Problem

Throughout the paper, sets will be denoted by uppercase letters, while members of a set will be denoted by lowercase letters. A collection of sets will be denoted using calligraphic letters.

The classic Set Basis Problem is defined as follows:

• Given: a collection C of subsets of a finite universe U, C = {C1, C2, . . . , Cm} and a positive integer K;
• Find: a collection B = {B1, . . . , BK} where each Bi is a subset of U, and for each Ci ∈ C, there exists a subcollection Bi ⊆ B, such that Ci = ∪B∈Bi B. In this case, we say Bi covers Ci, and we say C is collectively covered by B. Following common notations, B is referred to as a basis, and we call each Bj ∈ B a basis set. If Bj ∈ Bi and Bi covers Ci, we call Bj a contributor of Ci, and call Ci a sponsor of Bj . C1, C2, . . . , Cm are referred to as original sets.

Intuitively, similar to the basis vectors in linear algebra, which provides a succinct representation of a linear space, a set basis with smallest cardinality K plays the role of a compact representation of a collection of sets. The Set Basis Problem is shown to be NP-hard in [Stockmeyer, 1975]. We use I(C) to denote an instance of the set basis problem which finds the basis for C. A simple instance and its solution is reported in Table 1.

Most algorithms used in solving set basis problems are incomplete algorithms. These algorithms are based on heuristics that work well in certain domains, but often fail at covering sets exactly. For a survey, see Molloy et al [Molloy et al., 2009]. The authors of [Ene et al., 2008] implement the only complete solver we are aware of. The idea is to translate the set basis problem as a graph coloring problem, and then use existing graph coloring solvers. They also develop a useful prepossessing technique, which can significantly reduce the problem complexity.

The Set Basis Problem has a useful dual property, which has been implicitly used by previous researchers [Vaidya et al., 2006; Ene et al., 2008]. We formalize the idea by introducing Theorem 2.1.

Definition (Closure) For a collection of sets C, define the closure of C, denoted as C, which includes the collection of all possible intersections of sets in C:

• ∀Ci ∈ C, Ci ∈ C.
• For A ∈ C and B ∈ C, A∩B ∈ C.

Theorem 2.1. For original sets C = {C1, C2, . . . , Cn}, suppose {B1, . . . , BK} is a basis that collectively covers C. Define Ci = {Cj ∈ C|Bi ⊆ Cj}. Then B′i = ∩C∈Ci C (i = 1 . . .K) collectively covers C as well. Note for every B′i (i = 1 . . .K), B′i ∈ C.

One can check Theorem 2.1 by examining the example in Table 1. The full proof is available in the supplementary materials [Xue et al., 2015]. From the theorem, any set basis problem has a solution of minimum cardinality, where each basis set is in C. Therefore, it is sufficient to only search for basis within the closure C. Hence throughout this paper, we assume all basis sets are within its closure for any solutions to set basis problems. Theorem 2.1 also implies that each basis set Bi ∈ C is an intersection of all its sponsor sets. One can observe this fact in Table 1. It motivates our dual approach to solve the set basis problem, in which we search for possible sponsors for each basis set.

3 Parallel Scheme

The main intuition of our parallel scheme comes from an empirical observation on the structure of the solutions of the benchmark problems we considered. For each benchmark, we solve a series of related simplified subproblems, where we restrict ourselves to finding basis for a subset of the original collection of sets C. Interestingly, the solutions found by solving these subproblems are connected to the basis in the global solution. Although strictly speaking, the basis found for one sub-problem can only be expected to be a solution for that particular sub-problem, we observe empirically that almost all basis sets from sub-problems are supersets of one or more basis sets for the original, global problem. One intuitive explanation is as follows: Recall that from Theorem 2.1 each basis set can be obtained as the intersection of its sponsors. This fact applies both to the original global problem and its relaxed versions (subproblems). Since there are fewer sets to be covered in the subproblems, basis sets for the subproblems are likely to have fewer sponsors, compared to the ones for the global problem. When we take the intersection of fewer sets, we get a larger intersection. Hence we observe that it is often the case that a basis set for a subproblem is a superset of a basis set for the global problem.

Now suppose two subproblem basis sets A and B are both supersets of one basis set C in the global solution. If we intersect A with B, then the elements of C will remain in the intersection, but other elements from A or B will likely be removed. In practice, we can often obtain a basis set in the global solution by intersecting only a few basis sets from the solutions to subproblems.

Let us walk through the example in Table 1. First consider the subproblem consisting of the first 5 original sets, C0...C4. It can be shown that a minimum set basis is B1,1 = {x0, x1}, B1,2 = {x0, x3}, B1,3 = {x2, x3}, B1,4 = {x4}. As another subproblem we consider the collection of all original sets except for C0 and C2. We obtain a minimum basis B2,1 = {x0, x3}, B2,2 = {x2, x3}, B2,3 = {x3, x4}, B2,4 = {x0, x1, x3}. We see that each basis set of these two subproblems contains at least one of the basis sets of the original, full set basis problem. For example, B2 = {x0, x1} ⊆ B1,1 and B2 ⊆ B2,4. Moreover, one can obtain all basis sets except B0 for the original problem by intersecting these basis sets. For example, B3 = {x3} = B1,2 ∩ B2,2.

Given this observation, we design a parallel scheme that works in two phases - an exploration phase, followed by an aggregation phase. The whole process is shown in Figure 1, and the two phases are detailed in subsequent subsections.

Exploration Phase: we use a set of parallel processes. Each one solves a sub-problem obtained by restricting the global problem to finding the minimum basis for a subset of the original collection of sets C.

Aggregation Phase: we first identify an initial solution can-

1

1

1

2

2

2

s

s

s

didate by looking at all the possible intersections among basis sets found by solving sub-problems in the exploration phase. We then use this candidate to initialize a complete solver, which expands the search in an iterative deepening way to achieve completeness, iteratively adding portions of the search space that are "close" to the initial candidate.

3.1 Exploration Phase

The exploration phase utilizes parallel processes to solve a series of sub-problems. Recall the global problem is to find a basis of size K for the collection C. Let {C1, C2, . . . , Cs} be a decomposition of C, which satisfies C = C1 ∪ C2 ∪ . . . ∪ Cs. The sub-problem I(Ci) restricted on Ci is defined as:

• Given: Ci ⊆ C;
• Find: a basis Bi = {Bi,1, Bi,2, . . . , Bi,ki } with smallest cardinality, such that every set C ′ ∈ Ci is covered by the union of a sub-collection of Bi.

The sub-problem is similar to the global problem, however, with one key difference: we are solving an optimization problem where we look for a minimum basis, as opposed to the global problem, which is the decision version of the Set Basis Problem. In practice, the optimization is done by repeatedly solving the decision problem, with increasing values of K. We observe empirically that the optimization is crucial for us to get meaningful basis sets to be used in later aggregation steps. If we do not enforce the minimum cardinality constraint, the problem becomes under-constrained and there could be redundant basis sets found in this phase, which have no connections with the ones in the global solution.

Sets C1, C2, . . . , Cs need not be mutually exclusive in the decomposition. We group similar sets into one subproblem in our algorithm, so the resulting subproblem will have a small basis. To obtain collection Ci for the i-th subproblem, we start from an initial collection of a singleton Ci = {Ci1 }, where Ci1 is a randomly picked element from C. We then add T -1 sets most similar to Ci1 , using the Jaccard similarity coefficient |A∩B| |A∪B| . This results in a collection of T sets which look similar. Notice that this is our method to find a collection of similar sets. We expect other approach can work equally well.

3.2 Aggregation Phase

In the aggregation phase, a centralized process searches for the exact solution, starting from basis sets that are "close" to a candidate solution selected from the closure of basis sets found by solving sub-problems, and then expands its search space iteratively to achieve completeness.

To obtain a good initial candidate solution, we begin with a pre-solving step, in which we restrict ourselves to find a good global solution only within the closure of basis sets found by solving sub-problems. This is of course an incomplete procedure, because the solution might lie outside the closure. However, due to the empirical connections between the basis sets found by parallel subproblem solving and the ones in the final solution, we often find the global solution at this step.

If we cannot find a solution in the pre-solving step, the algorithm continues with a re-solving step, in which an iterative deepening process is applied. It starts with the best K basis sets found in the pre-solving step 1 , and iteratively expands the search space until it finds a global solution. The two steps are detailed as follows.

Pre-solving Step

Suppose Bi is the basis found for sub-problem I(Ci), let B0 = ∪si=1Bi and B0 be the closure of B0. The algorithm solves the following problem in the pre-solving step:

• Given: B0;
• Find: Basis B∗ = {B∗ 1 , . . . ,B∗K } from B0, such that B∗ 1 , . . . ,B∗K minimizes the total number of uncovered elements and falsely covered elements in C 2 .

In practice B0 is still a huge space, so this optimization problem is hard to solve. We thus apply an incomplete algorithm, which only samples a small subset U ⊆ B0 and then select the best K basis sets from U. It does not affect the later re-solving step, since it can start the iterative deepening process from any B∗ , whether optimal in B0 or not.

The incomplete algorithm first forms U by conducting multiple random walks in the space of B0. Each random walk starts with a random basis set B ∈ B0, and randomly intersects it with other basis sets in B0 to obtain a new member in B0. All these sets are collected to form U. With probability p, the algorithm chooses to intersect with the basis which maximizes the cardinality of the intersection. With probability (1 - p), the algorithm intersects with a random set. In our experiment, p is set to 0.95, and we repeat this random walk several times with different initial sets to make U large enough. Next the algorithm selects the optimal basis of size K from U which maximizes the coverage of the initial set collection, using a Mixed Integer Programming (MIP) formulation. The pseudocode of the incomplete algorithm and the MIP formulation are both in the supplementary materials [Xue et al., 2015].

1 Best in terms of coverage of the initial set collection.

2 An uncovered element of set Cj is one element contained in Cj , but is not covered by any basis set that are contributors to Cj . A falsely covered element of set Cj is one element that is in one basis set that is a contributor to set Cj, but is not contained in Cj.

Re-solving Step

The final step is the re-solving step. It takes as input the basis B∗ = {B∗ 1 , B∗ 2 , . . . ,B∗K } from the pre-solving step, and searches for a complete solution to I(C) in an iterative deepening manner. The algorithm starts from a highly restricted space D1, which is a small space close to B∗ . If the algorithm can find a global solution in D1, then it terminates and returns the solution. Otherwise, it expands its search space to D2, and searches again in this expanded space, and so on. At the last step, searching in Dn is equivalent to searching in the original unconstrained space C, which is equivalent to solving the global set-basis problem without initialization at all. However, this situation is rarely seen in our experiments.

In practice, D1, . . . , Dn are specified by adding extra constraints to the original MIP formulation for the global problem, then iteratively removing them. Dn corresponds to the case where all extra constraints are removed.

The actual design of D1, . . . , Dn relies on the MIP formulation. In our MIP formulation (which is detailed in the supplementary materials [Xue et al., 2015]), there are indicator variables yi,k (1 ≤ i ≤ n and 1 ≤ k ≤ K), where yi,k = 1 if and only if the i-th element is contained in the k-th basis set Bk. We also have indicator variables zk,j, where zk,j is one if and only if the basis set Bk is a contributor of the original set Cj (or equivalently, Cj is a sponsor set for Bk).

Because we empirically observe that B∗ 1 , B∗ 2 , . . . ,B∗K are often super-sets of the basis sets in the exact solution, we construct the constrained space D1, . . . ,Dn by enforcing the sponsor sets of certain basis sets. Notice that this is a straightforward step in the MIP formulation, since we only need to fix the corresponding indicator variables zk,j to 1 to enforce Cj as a sponsor set for Bk. The hope is that these clamped variables will include a subset of backdoor variables for the original search problem [Williams et al., 2003; Dilkina et al., 2009; Hamadi et al., 2011]. The runtime of the sequential solver is dramatically reduced when the aggregation phase is successful in identifying a promising portion of the search space.

As pointed out by Theorem 2.1, we can represent B∗ 1 , B∗ 2 , . . . ,B∗K in terms of their sponsors:

in which C11,C12, . . . ,C21, . . . ,CKsK are all original sets in collection C. For the first restricted search space D1, we enforce the constraint that the sponsors for the i-th basis set Bi must contain all the sponsors of B∗ i for all i ∈ {1, . . . ,K}. Notice this implies Bi ⊆ B∗ i .

In later steps, we gradually relax these extra constraints, by freeing some of the indicator variables zk,j's which were clamped to 1 in previous steps. Dn denotes the search space when all these constraints are removed, which is equivalent to searching the entire space. The last thing is to decide the order used to remove these sponsor constraints. Intuitively, if one particular set is discovered many times as a sponsor set in the solutions to subproblems, then it should have a high chance to be the sponsor set in the global solution, because it fits in the solutions to many subproblems. Given this intuition, we associate each sponsor set with a confidence score, and define n thresholds: 0 = p1 < . . . < pn = +∞. In the k-th round (search in Dk), we remove all the sponsor sets whose confidence score is lower than pk. We define the confidence score of a particular set as the number of times it appears as a sponsor of a basis set in subproblem solutions, which can be aggregated from the solutions to subproblems.

4 Experiments

We test the performance of our parallel scheme on both the classic set basis problem, and on a novel application in materials science.

4.1 Classic Set Basis Problem

Setup We test the parallel scheme on synthetic instances. We use a random ensemble similar to Molloy et al [Molloy et al., 2009], where every synthetic instance is characterized by n, m, k, e, p. To generate one synthetic instance, we first generate k basis sets. Every set contains [n p 100 ] objects, uniformly sampled from a finite universe of n elements. We then generate m sets. Each set is a union of e randomly chosen basis sets from those initially generated.

We develop a Mixed Integer Programming (MIP) model to solve the set basis problem (detailed in the supplementary materials [Xue et al., 2015]). The MIP model takes the original sets C and an integer K, and either returns a basis of size K that covers C exactly, or reports failure. We compare the performance of the MIP formulation with and without the initialization obtained using the parallel scheme described in the previous section.

We empirically observe high variability in the running times of the sub-problems solved in the exploration phase, as commonly observed for NP-hard problems. Luckily, our parallel scheme can still be used without waiting for every sub-problem to complete. Specifically, we set up a cut-off threshold of 90%, such that the central process waits until 90% of sub-problems are solved before carrying out the aggregation phase. We also run 5 instances of the aggregation phase in parallel, each with a different random initialization, and terminate as soon as the fastest one finds a solution. In our experiment, n = m = 100. All MIPs are solved using IBM CPLEX 12.6, on a cluster of Intel x5690 3.46GHz core processors with 4 gigabytes of memory. We let each subproblem contain T = 15 sets for all instances.

Results Results obtained with and without initialization from parallel sub-problem solving are reported in Table 2. First, we see it takes much less wall-clock time (typically, by several orders of magnitude) for the complete solver to find the exact solution if it is initialized with the information collected from the sub-problems. The improvements are significant even when taking into account the time required for solving sub-problems in the exploration phase. In this case, we obtain several orders of magnitude saving in terms of solving time. For example, it takes about 50 seconds (wall-clock time) to solve A6, but about 5 hours without parallel initialization. Because we run s = 100 sub-problems in the exploration phase, another comparison would be based on CPU time, which is given by (100·Exploration+5·Aggregation). Under this measurement, our parallel scheme still outperforms the sequential approach on problem instances A2, A3, A5, A6, A8. Even though our CPU time is longer for some instances, our parallel scheme can be easily applied to thousands of cores. As parallel resources are becoming more and more accessible, it is obvious to see the benefit of this scheme. Note that we can also exploit at the same time the built-in parallelism of CPLEX to solve these instances. However, because CPLEX cannot explore the problem structure explicitly, it cannot achieve significant speed-ups on many instances. For example, it takes 12813.15, 259100.25 and 113475.12 seconds to solve the largest A6, A7 and A8 instances using CPLEX on 12-cores.

3 For this instance, 73 out of 100 subproblem instances complete within 2 hours. Thus the aggregation phase is conducted based on these instances. This exploration time here is calculated based on the slowest of the 73 instances.

Although our focus is on improving the run-time of exact solvers, Table 2 also shows the performance of several stateof-the-art incomplete solvers on these synthetic instances. We implemented FastMiner from [Vaidya et al., 2006], ASSO from [Miettinen et al., 2008], and HPe from [Ene et al., 2008], which are among the most widely used incomplete algorithms. FastMiner and ASSO take the size of the basis K as input, and output K basis sets. They are incomplete in the sense that their solution may contain false positives and false negatives, which are defined as follows. c is a false positive element if c ̸∈ Ci, but c is in one basis set that is a contributor to Ci. c is a false negative element, if c ∈ Ci, but c is not covered by any basis sets contributing for Ci. FastMiner does not provide the information about which basis set contributes to an original set. We therefore give the most conservative assignment: B is a contributor to Ci if and only if B ⊆ Ci. This assignment introduces no false positives. Both FastMiner and ASSO have parameters to tune. Our report are based on the best parameters we found. We report the maximum error rate in Table 2, which is defined as maxCi∈C{(fti + ffi)/|Ci|}, where fti and ffi are the number of false positive and false negative elements at Ci, respectively. As seen from the table, neither of these two algorithms can recover the exact solution. ASSO performs better, but it still has 51.06% error rate on the hardest benchmark. We think the reason why FastMiner performs poorly is because it is usually used in situations where certain number of false positives can be tolerated. HPe is a graph based incomplete algorithm. It is guaranteed to find a complete cover, however it might require a number of basis sets K′ larger than the optimal number K. We implemented both the greedy algorithm and the lattice-based post-improvement for HPe, and we used the best heuristic reported by the authors. As we can see from Table 2, HPe often needs five times more basis sets to cover the entire collection.

The authors in [Ene et al., 2008] implemented the only complete solver we are aware of. Unfortunately, we can not obtain their code, so a direct comparison is not possible. However, the parallel scheme we developed does not make any assumption on the specific complete solver used. We expect other complete solvers (in addition to the MIP one we experimented with) will improve from the initialization information provided by solving subproblems.

Discussion We now provide empirical evidence that justifies and explains the surprising empirical effectiveness of our method. For clarity, we call a basis set found by solving subproblems an s-basis set, and a basis set in the global solution a g-basis set. For any set S, we define the hitting rate as: p(S) = maxB∈B |S ∩B|/|B|, and the inverse hitting rate as: ip(S) = maxB∈B |S ∩ B|/|S|, where B is chosen among all g-basis set. Intuitively, p(S) and ip(S) measure the distance between S and the closest g-basis set. Note that p(S) = 1 (respectively ip(S) = 1) implies the basis S is the superset (respectively subset) of at least one g-basis in the global solution. If p(S) and ip(S) are both 1, then S matches exactly to one g-basis set 4 .

First, we study the overlap between a single s-basis set and all g-basis sets. As shown in Figure 2, across the benchmarks we considered the hitting rate is almost always one (with the lowest mean is for A2, which is 0.9983). This means that the s-basis sets are almost always supersets of at least one gbasis set in the global solution.

Next, we study the relationship between the intersection of multiple s-basis sets and g-basis sets. Figure 3 shows the median hitting rate and inverse hitting rate with respect to different number of s-basis sets involved in one intersection. The error bars show the 10-th and 90-th percentile. The result is averaged over all instances A1 through A8, with equal number of samples obtained from each instance. In the top chart, the s-basis sets involved in one intersection are supersets of one common g-basis set. In this case, the hitting rate is always 1. However, by intersecting only a few (2 or 3) s-basis sets, the inverse hitting rate becomes close to 1 as well, which implies the intersection becomes very close to an exact match of one g-basis set. This is in contrast with the result in the bottom chart, where the intersection is among randomly selected s-basis sets. In this case, when we increase the size of the intersection, fewer and fewer elements remain in the intersection. The bottom chart of Figure 3 shows the percentage of elements left, defined as |∩ki=1 Ai|/maxki=1 |Ai|. When intersecting 5 basis sets, in median case less than 10% elements still remain in the intersection.

4 Assuming no g-basis set is the subset of another g-basis set, which is the case in instances A1 to A8.

The top and bottom charts of Figure 3 provide an empirical explanation for the success of our scheme: as we randomly intersect basis sets from the solutions to the subproblems, some intersections become close to the empty set (as in the bottom chart case), but others converge to one of the g-basis sets in the global solution (as in the upper chart case). In the second case, we obtain good solution candidates for the global problem by intersecting solutions to subproblems.

4.2 Phase Identification Problem

We also apply our parallel scheme to speed up solvers for a variant of the set basis problem with extra constraints. We show how this more general formulation can be applied to the so-called phase identification problem in combinatorial materials discovery [Le Bras et al., 2011].

In combinatorial materials discovery, a thin film is obtained by depositing three metals onto a silicon wafer using guns pointing at three locations. As metals are sputtered on the silicon wafer, different locations have different combinations of the metals, due to their distances from the gun points. As a result, various crystal structures are formed across locations. Researchers then analyze the X-ray diffraction patterns (XRD) at a selected set of sample points. The XRD pattern at one sample point reflects the crystal structure of the underlying material, and is a mixture of one or more basis patterns, each of which characterizes one crystal structure. The overall goal of the phase identification problem is to explain all the XRD patterns using a small number of basis patterns.

The phase identification problem can be formulated as an extended version of the set basis problem. We begin by introducing some terminologies. Similar to [Ermon et al., 2012], we use discrete representations of the XRD signals, where we characterize each XRD pattern with the locations of its peaks. In this model, we define a peak q as a set of (sample point, location) pairs: q = {(si, li)|i = 1, . . . , nq}, where {si|i = 1, . . . , nq} is a set of sample points where peak q is present, and li is the location of peak q at sample point si, respectively. We use the term phase to refer to a basis XRD pattern. Precisely, a phase comprises set of peaks that occur in the same set of sample points. We use the term partial phase to refer to a subset of the peaks and/or a subset of the sample points of a phase. We use lower-case letters p, q, r to represent peaks, and use upper-case letters P, Q, R to represent phases. Given these definitions, the Phase Identification Problem is:

Given A set of X-ray diffraction patterns representing different material compositions and a set of detected peaks for each pattern; and K, the expected number of phases.

Find A set of K phases, characterized as a set of peaks and the sample points in which they are involved.

Subject to Physical constraints that govern the underlying crystallographic process. We use all the constraints in [Ermon et al., 2012]. For example, one physical constraint is that a phase must span a continuous region in the silicon wafer.

Figure 4 shows an illustrative example. In this example, there are 4 peaks for phase α, and 3 peaks for phase β. Peaks in phase α exist in all sample points in the green region, and peaks in phase β exist in purple region. They co-exist in several sample points in the mid-right region of the triangle.

There is an analogy between the Phase Identification Problem and the classical Set Basis Problem. In the Set Basis Problem, each original set is the union of some basis sets. In the Phase Identification Problem, the XRD pattern at a given sample point is a mixture of several phases. Here, the phase is analogous to the basis set, and the XRD pattern at a given sample point is analogous to the original set. Because of this relationship, we employ a similar parallel scheme to solve the Phase Identification Problem, which also includes an exploration phase followed by an aggregation phase.

Exploration Phase

In the Exploration Phase, a set of subproblems are solved in parallel. For the Phase Identification Problem, a subproblem is defined as finding the minimal number of phases to explain a contiguous region of sample points on the silicon wafer.

This is analogous to the exploration phase defined for set basis problem - finding basis for a subset of sets. The reason why we emphasize a contiguous region is because of the underlying physical constraint: the phase found must span a contiguous region in the silicon wafer. Figure 5 shows a sample decomposition into subproblems. Here each colored small region represents a subproblem.

Aggregation Phase

The exploration phase produces a set of partial phases from solving subproblems. We call them partial because each of them describes only a subset of sample points.

As in the Set Basis Problem, we find partial phases can be merged together into larger phases. Figure 4 shows an illustrative example. Formally, two phases A and B may be merged into a new phase C, denoted as C = A ◦ B, which contains all the peaks from A and B whose locations match across all the sample points they both present. The peaks in C then span the union of sample points of A and B. The merge operator ◦ plays the same role as the intersection operator of the Set Basis Problem. Similarly, we define S as the closure of (partial) phases S with respect to the merge operator ◦, which generates all possible merging of the phases in S.

Suppose B0 = ∪si=1Bi is the set of all (partial) phases identified by solving subproblems, where Bi is the set of (partial) phases identified when solving subproblem i. As with the Set Basis Problem, the aggregation phase also has a pre-solving step, and a re-solving step. The pre-solving step takes as input the responses B0 from all subproblems, and extracts a subset of K partial phases from the closure B0 as the candidate solution, which explains as many peaks on the silicon wafer as possible. The re-solving step searches in an iterativedeepening way for an exact solution, starting from the phases close to the candidate solution from the pre-solving step.

As in the pre-solving step of the Set Basis Problem, B0 could be a large space and we are unable to enumerate all items in B0 to find an exact solution. Instead, we take an approximate approach which first expands B0 to a larger set B′ ⊆ B0 using a greedy approach. Then we employ a MixedInteger Program (MIP) formulation that selects the best K phases from B′ which covers the largest number of peaks. The greedy algorithm and the MIP encoding are similar in concept to the ones used in solving the Set Basis Problem, but take into account extra physical constraints.

The Re-solving step expands the search from the presolving step in an iteratively deepening way to achieve completeness. Suppose the pre-solving step produces K phases P ∗ 1 , P ∗ 2 , . . . , P ∗K . In the first round of the re-solving step, the complete solver is initialized such that the first phase must contain all the peaks of P ∗ 1 , the second phase must contain all the peaks of P ∗ 2 , etc. If the solver can find a solution with this initialization, then the solver terminates and returns the results. Otherwise, it usually detects a contradiction very quickly. In this case, we remove some peaks from P ∗ 1 , . . . , P ∗K and re-solve the problem. We continue this resolving process, until all the peaks from the Pre-solving step are removed, in which case the solver is free to explore the entire space without any restrictions. Again, this is highly unlikely in practice. In most cases, the solver is able to find solutions in the first one or two iterations.

We augmented the Satisfiability Modulo Theory formulation as described in [Ermon et al., 2012] with our parallel scheme and use the Z3 solver [De Moura and Bjørner, 2008] in the experiments. We use Z3 directly in the exploration phase, and then use it as a component of an iterative deepening search scheme in the aggregation phase. Due to a rather more imbalanced distribution of the running times across different sub-problems, we only wait for 50% of sub-problem solvers to complete before conducting the aggregation phase.

Table 3 displays the experimental results for the phase identification problem. We run on the same benchmark instances used in the work of Ermon et al [Ermon et al., 2012]. We can see from Table 3 that in all cases the solver completes much faster when initialized with information obtained by parallel subproblem solving. This improvement in the runtime allows us to analyze much bigger problems than previously possible in combinatorial materials discovery.

5 Conclusion

We introduced a novel angle for using parallelism to exploit hidden structure of hard combinatorial problems. We demonstrated empirical success in solving the Set Basis Problem, obtaining over an order of magnitude speedups on certain problem instances. We also identified a novel application area of the Set Basis Problem, concerning the discovery of new materials for renewable energy sources. Future directions include applying this approach to other NP-complete problems, and exploring its theoretical foundations.

Acknowledgments

We are thankful to the anonymous reviewers, Richard Beinstein and Ronan Lebras for their constructive feedback. This work was supported by NSF Expeditions grant 0832782, NSF Infrastructure grant 1059284, NSF Eager grant 1258330, NSF Inspire grant 1344201, and ARO grant W911-NF-14-1-0498.