Given a set of taxa S and a complete set of quartet topologies Q over S, the problem of determining a phylogeny that satisfies the maximum number of topologies is called the Maximum Quartet Consistency (MQC) problem. The MQC problem is an NP-problem. MQC has been solved both heuristically and exactly. Exact solutions for MQC include Constraint Programming, Answer Set Programming and Pseudo-Boolean Optimization (PBO). This paper extends the range of solutions for MQC, by developing two new PBO models and also by introducing models based on the optimization of Satisfiability Modulo Theories (SMT). The models were experimentally compared with existing exact solutions. The results show that for instances with small percentage of quartet errors, the models based on SMT can be competitive, whereas for instances with higher number of quartet errors the PB-models are more efficient.
Evolutionary trees (or phylogenies) of a given set of taxa (a set of biological species) are often used for studying the evolutionary relationship between the given taxa. Current phylogenetic analyses are carried out using biological data (such as DNA or protein sequences). Since the amount of available data is different for different species, biologists have to balance the number of taxa used and the amount of available data common to all the species. This is know as the Data Disparity problem [1]. In order to overcome the Data Disparity problem, quartet based methods have been suggested.
Quartet based methods divide the task of inferring an evolutionary tree in two subtasks. First, a set of trees for subsets of four taxa is obtained. Then from the trees obtained a global tree is constructed. Considering only four taxa in the first step, maximizes the amount of available data used, and the confidence on the tree produced. These evolutionary trees for subsets of four taxa are called quartet topologies. The quartet topologies produced may still be conflicting or even missing. Given a set of quartet topologies, the problem of computing an evolutionary tree that agrees with the maximum number of quartet topologies is referred to as the Maximum Quartet Consistency (MQC) problem. In recent years, solutions based on Answer-Set Programming (ASP) have been identified as the most efficient [18]. This paper proposes alternative combinatorial optimization solutions for the MQC problem. The first set of models extend recent work on using Pseudo-Boolean b d [a, b|c, d] a c d [a, c|b, d] b a d c [a, d|b, c] b
Optimization (PBO) [10]. Moreover, the paper also proposes models based on Satisfiability Modulo Theories (SMT).
The paper is organized as follows. The next section introduces the MQC problem and surveys representative solutions. Afterwards, Section 3 summarizes recent work on using PBO for solving the MQC problem, and develops new PBO models. The use of SMT for MQC is described in Section 4. Section 5 analyzes experimental results on real and artificially generated instances. Finally, the paper concludes in Section 6. 2
This section defines the problem of Maximum Quartet Consistency (MQC), and reviews existing solutions, either heuristic or exact.
Consider a given set of taxa S = {s1, . . . , sn}. A phylogeny on S is a tree such that each leave is mapped to a taxon in S. Each internal node in the phylogeny represents an extinct or hypothesized ancestor of the taxa in the subtrees that the internal node connects. A phylogeny may be rooted or unrooted. In a rooted phylogeny a link between two nodes represent a parentchild relation, in which the root of the phylogeny is a common ancestor of all the taxa. A phylogeny is called binary or resolved if every internal node has degree three (in a rooted phylogeny the root node has degree two). In the rest of the paper the term phylogeny refers to an unrooted resolved phylogeny, unless explicitly stated.
A subset of four different taxa is called a quartet. A phylogeny for a quartet is called a quartet topology (simply refered as topology). In this paper, for a quartet {a, b, c, d} ⊆ S three different topologies are considered (up to symmetry) denoted by [a, b|c, d], [a, c|b, d] and [a, d|b, c] as represented in Figure 1. Consider a quartet q ⊆ S and a phylogeny T on S. The topology obtained from T by removing all nodes not in the paths connecting the nodes mapped to taxa in q is called the topology of q derived from T . On the left side of Figure 2, a phylogeny is shown where the paths that connect taxa a, b, c and f have been marked with dotted lines. In this example the derived topology for quartet {a, b, c, f } is [a, b|c, f ] (respresented on the right side of Figure 2).
A topology qt is said to be consistent with a phylogeny T or T satisfies qt if the derived topology for the quartet of qt from T is the same as qt. The topology in Figure 2 is consistent with the phylogeny in the same figure.
Given a set of topologies Q whose set of taxa is S, Q is said to be complete g a b f c e d a b c f if Q contains a topology for each of the C4n quartets of S; otherwise Q is said to be incomplete. The set Q is said to be compatible if there exists a phylogeny T on S that satisfies all the topologies in Q, and T is called the associated phylogeny of Q. The MQC problem can now be defined as: Definition 1 (Maximum Quartet Consistency problem (MQC)). Given a set of taxa S = {s1, . . . , sn} and a set of topologies Q on S, the Maximum Quartet Consistency problem is the problem of determining a phylogeny T on S that satisfies the maximum number of topologies of Q.
The MQC problem is a NP-hard problem [2]. If Q is complete, then MQC admits a polynomial-time approximation scheme [8]. If Q is incomplete then MQC is MAX SNP-hard [8]. Existing approaches for the MQC problem can be categorized as either heuristic or exact. Heuristic approaches return an approximation to the optimal phylogeny. Among the many heuristics proposed over the years, a few are overviewed next. The short quartet methods of Erdos et al. [5], the use of semi-definite programming by BenDor et al. [1], the quartet cleaning algorithms of Berry et al. [2]. More recently, in 2008, Wu et al. [15] proposed three algorithms that return a phylogeny with a probability of success.
This paper focuses on exact solutions for the MQC problem, which guarantee that the phylogeny returned satisfies the maximum number of topologies. Existing exact approaches can be categorized in four main classes. Bendor et al. [1] proposed the use of Dynamic Programming to solve the problem of MQC where the topologies given have weights associated to them. Their objective is to compute a phylogeny with a maximal score, that is, a phylogeny whose sum of weights of the satisfied topologies is maximal.
The Fixed-Parameter Algorithm was proposed by Gramm and Niedermeier [7]. The objective of the algorithm is to compute a phylogeny that fails to satisfy no more than a given k topologies. Their algorithm works in time O(4kn + n4), where k is the given number of errors and n = |S|.
Wu et al. [16] proposed the use of Answer Set Programming (ASP) for MQC. The proposed approach searches for an ultrametric matrix that satisfies the maximum number of topologies. From the ultrametric matrix it is possible to obtain a phylogeny in polynomial time.
More recently, in 2005, Wu et al. [17] presented a lookahead branch-andbound algorithm for MQC. The algorithm can be seen as an improvement over the Fixed-Parameter Algorithm [7], and has the same running time of O(4kn + n4). However, in contrast with the Fixed-Parameter Algorithm, the number of quartet errors is not required to be known.
Finally, in 2008, Morgado and Marques-Silva [10] translated the problem of MQC to Pseudo-Boolean Optimization (PBO) . 3
Pseudo-Boolean Optimization (PBO) models for MQC were first proposed in [10]. This section reviews in some detail the PBO models followed by the description of two new PBO models based on the existing models. 3.1
Given a complete set of topologies Q over a set of taxa S, the objective of the PBO models for MQC is to encode a phylogeny T of S that satisfies the maximum number of topologies in Q. In order to obtain T , and similarly to the use of ASP, the PBO models encode an ultrametric matrix M . From the ultrametric matrix M it is possible to obtain a rooted phylogeny T .
The ultrametric matrix M is a square symmetric matrix n × n (where n = |S|) whose entry values are labels that range from 0 to n − 1. Each row and column of the ultrametric matrix M is indexed by the given taxa, and the value of an entry M (i, j) is the label of the internal node in T corresponding to the lowest common ancestor between leaf nodes si and sj. The labels of the internal nodes of T have the property that for any path from root node to a leaf, the corresponding labels are in decreasing order. To guarantee that M is an ultrametric matrix, M must satisfy the following three properties [16]: (a) all values in the matrix must be between 1 and n, except for the elements on the main diagonal which must be 0, that is M (i, i) = 0; (b) M is symmetric, thus M (i, j) = M (j, i); and (c) for each triple (i, j, l) such that 1 ≤ i, j, l ≤ n then
M (i, j) = M (i, l) ∧
M (i, l) > M (j, l) ∨ M (i, j) = M (j, l) ∧
M (j, l) > M (i, l) ∨ M (j, l) = M (i, l) ∧
M (i, l) > M (i, j)
(
The PBO models can be divided in three main parts: encoding of the ultrametric matrix; encoding of the consistency of the topologies; and encoding of the cost function. 3.1.1
We start by detailing the Basic PBO model. To encode the ultrametric
matrix, the Basic PBO model considers the boolean variables Mi,j,k where
1 ≤ i < j ≤ n and 1 ≤ k ≤ ⌈ n2 ⌉. The semantics of variable Mi,j,k is
Mi,j,k = 1 if M (i, j) = k, otherwise Mi,j,k = 0. For a given pair (i, j) one
and only one k is allowed to have Mi,j,k = 1. To ensure this condition, the
Basic PBO model introduces the following constraints:
n
⌈ 2 ⌉
X Mi,j,k = 1
k=1
(
Consider variable c1i,j,l and the corresponding condition (
c1i,j,l ≡ AN D(c1i(,=j,)l, c1i(,>j,)l)
The values of variables c1i(,=j,)l and c1i(,>j,)l are encoded as comparator circuits on the unary representation of M (i, j), M (i, l) and M (j, l) (that is, using variables Mi,j,k, Mi,l,k and Mj,l,k). Variables c2i,j,l and c3i,j,l are encoded in an analogous way.
According to [16], an ultrametric matrix M satisfies a topology [i, j|l, m] if and only if one of the following conditions is satisfied: (M (i, l) > M (i, j) ∧
M (j, m) > M (i, j)) or (M (i, l) > M (l, m) ∧
M (j, m) > M (l, m))
The Basic PBO model considers a new boolean variable qt that encodes
the consistency of the t-th topology. Suppose that the t-th topology in Q is
topology [i, j|l, m]. The model encodes the value of qt by introducing two
new boolean variables q1t, q2t and associates them with conditions (
The Reduced PBO model is based on the Basic PBO model, but where
auxiliary variables are used that allow a simpler encoding of greater than
conditions (M (i, j) > M (l, m)) using fewer variables. In the Reduced model,
constraint (
that P⌊kl=og02(⌈ n2 ⌉)⌋ Mi,j,k ≤ ⌈ n2 ⌉.
The last PBO model proposed in [10] was the Binary PBO model. The rationale of the Binary PBO model is to replace the encoding of M (i, j) by changing the semantics of variables Mi,j,k. In the Binary PBO model, variable Mi,j,k represents the k-th bit of the binary representation of M (i, j). As a consequence the limits of index k are reduced to 0 ≤ k ≤ ⌊log2(⌈ n2 ⌉)⌋, therefore reducing the number of variables. ⌊log2(⌈ n2 ⌉)⌋ 2k × Mi,j,k.
Variables M (i, j) are now encoded as M (i, j) = Pk=0
Moreover, all the constraints of the Basic PBO model are maintained (using
the new limits of index k) except for constraints (
This section proposes two new PBO models, namely the Unary PBO model and the Binary PBO model with a dedicated comparator. 3.2.1
The rationale of the Unary model is to change the semantics of the selection variables Mi,j,k and merge them with the semantics of the auxiliary variables
k 1 2 3 4 5 6 si,j 1 1 1 1 0 0 k (b) sik,j.
The objective is for the new variables to be able to represent the values in the ultrametric matrix and still allow the simplified encoding of M (i, j) < M (l, m).
Consider for example that n = 12, then k ≤ ⌈ n ⌉ = 6. For some position 2 (i, j) in the matrix, suppose that M (i, j) = 4. A valid valuation of the selection variables and of the auxiliary variables is given in Figure 3 (a).
The Unary model uses only the auxiliary variables sik,j, but with a new semantics, which is if M (i, j) = a for same a, 1 ≤ a ≤ ⌈ n2 ⌉, for all 1 ≤ k ≤ a, then si,j has value 1, and for all 1 < k ≤ ⌈ n2 ⌉, then sik,j has value 0.
k
M (i, j) can then be obtained by: M (i, j) = ( P⌈k=n2⌉1−1[(sik,j − sik,+j1) × k] if (P⌈k=n2⌉1−1[(sik,j − sik,+j1) × k]) 6= 0 ⌈ n2 ⌉ otherwise For the previous example, the new valuation is represented in Figure 3 (b). Given the new semantics, condition M (i, j) > M (l, m) is true if there is a k′ such that sik,′j = 1 and slk,′m = 0. The Reduced model considers new boolean variables e(i,j)(l,m) The semantics of variables e(i,j)(l,m) is e(ki,j)(l,m) has value k k 1 if sik,j = 1 and slk,m = 0, otherwise e(i,j)(l,m) has value 0. Variables e(i,j)(l,m) k k are encoded as AND gates, e(i,j)(l,m) ≡ AN D(sik,j, N OT (slk,m)).
k
The greater than condition M (i, j) > M (l, m) is then encoded by a new variable GT(i,j)(l,m) whose value is given as an OR gate of variables e(ki,j)(l,m), that is: k GT(i,j)(l,m) ≡ OR(e(ki,j)(l,m) : 1 ≤ k ≤ ⌈ ⌉) 2
As with the Reduced model, the cost function and all the constraints of ultrametricity of matrix M and of quartet consistency of the PB model are maintained, but making use of variables GT(i,j)(l,m) as appropriate. 3.2.2
Sinz [14] presented a specialized comparator circuit to represent that a binary word (bmbm−1 . . . b0) is smaller or equal to a number n. An optimization of the Binary PBO model corresponds to use the comparator circuit proposed by [14] for encoding that M (i, j) ≤ ⌈ n2 ⌉. Suppose ⌈ n2 ⌉ is given by the binary word (bmbm−1 . . . b0). The recursive equations of [14] (that obtain a set of clauses) can be modified to obtain a set of sums: Cp(b0) = {1 − Mi,j,0}, b0 = 0 ∅ , b0 = 1 , Cp(br) = {1 − Mi,j,r} ∪ Cp(br−1) , br = 0 {1 − Mi,j,r} ⊗ Cp(br−1), br = 1 where A ⊗ B = {x + y|x ∈ A, y ∈ B}. Then the set of conditions to add to the model is the set {s ≥ 1|s ∈ Cp(bm)}.
All the constraints of the Binary PBO model are maintained except for the comparator circuit that encodes the restriction M (i, j) ≤ ⌈ n2 ⌉, which is now encoded according to the above equations. 4
This section outlines the use of Satisfiability Modulo Theories (SMT) for solving the MQC problem, by describing how to optimize a linear cost function subject to a SMT formula. This solution is referred to as the maximum Satisfiability Modulo Quartet Consistency (mSMQC) model. The mSMQC model follows the same principles of the PBO models of Section 3, but using the modelling capabilities of SMT. Thus the mSMQC model encodes an ultrametric matrix M and checks for the consistency of the topologies given, against the values in M .
The mSMQC model can be divided in four parts:(
Domain Restriction. Consider a given set of topologies Q over the set of taxa S = {s1, . . . , sn}. In the PB-models the solution to the MQC problem was given by encoding the resulting phylogeny with an ultrametric matrix M . In the mSMQC model, each position of the matrix is associated with one integer variable containing the value in that position. For each pair of taxa (si, sj) ∈ S, such that 1 ≤ i < j ≤ n, the mSMQC model considers n the integer variable Mi,j. Variable Mi,j has value k with 1 ≤ k ≤ ⌈ 2 ⌉ if and only if position (i, j) of M has value k.
To constrain the domain of variables Mi,j (1 ≤ i < j ≤ n), the mSMQC model introduces the following constraints: n ((Mi,j ≥ 1) ∧ (Mi,j ≤ ⌈ ⌉))
2
(
All the constraints added to the mSMQC model are constraints of Linear Integer Arithmetic (LIA). 4.2
For the PBO models of Section 3, the maximization was achieved through a cost function. SMT is a decision problem and as such, there is no cost function to be optimized. Nevertheless, recent work addressed optimization problems in SMT, including weighted MAX-SMT [11]. In the following we consider how to maximize the number of topologies satisfied by either using bounds or MAX-SMT. 4.2.1
Let nsqt be the number of satisfied topologies. Thus 1 ≤ nsqt ≤ |Q|, that is at least one topology is satisfied and at most all topologies are satisfied.
The motivation for Bounding by Bottom-Up is to search for the maximum number of satisfied topologies iteratively, by creating instances of SMT and giving them to a SMT solver.
In the Bounding by Bottom-Up approach a lower bound lb nsqt for nsqt is maintained. This lower bound has the property that is always possible to satisfy at least lb nsqt topologies. The initial value of lb nsqt is 1. In each iteration, a SMT instance is created as described in Section 4.1 to which is added the constraint:
|Q|
(X qt > lb nsqt)
r=1
(
The way the lower bound is updated depends on the SMT solver used. If the SMT solver is capable of reporting an assignment to the variables every time time it returns SATISFIABLE, then the new lower bound is given by the number of topologies satisfied by the assignment reported. Otherwise, the lower bound can only be increased by one. 4.2.2
The idea of the Bounding by Top-Down is the opposite of the previous Bounding by Bottom-Up. Instead of iterating through satisfiable SMT instances, the Top-Down approach will iterate through a set of unsatisfiable SMT instances till a satisfiable instance is found.
The Bounding by Top-Down maintains an upper bound up nsqt for nsqt.
This upper bound up nsqt has the property that is impossible to satisfy
(up nsqt + 1) topologies. The initial value of up nsqt is the total number of
topologies |Q|. In each iteration the constraint is added to the model:
|Q|
r=1
(X qt ≥ up nsqt)
(
The idea of the Mixed Bounding approach is to search for the maximum number of satisfied topologies using binary search and augment it with the Bounding by Bottom-Up approach.
The Mixed Bounding approach maintains both the lower bound lb nsqt
and the upper bound up nsqt (as in the previous bounding approaches).
Consider max nsqt to be the exact maximum number of topologies satisfied.
At each iteration, both the following conditions are true: (
In each iteration the mSMQC with the Mixed Bounding approach considers the value in the middle of the bounds. Let mdl nsqt = ⌊ up nsqt+lb nsqt ⌋. 2 If mdl nsqt = lb nsqt then max nsqt is lb nsqt and mSMQC returns. Otherwise mSMQC with Mixed Bounding adds the constraints: |Q| r=1 (X qt ≥ mdl nsqt), |Q| r=1 (X qt < up nsqt)
(
In MAX-SMT the overall formula is an unsatisfied conjunction of SMT subformulas. Each sub-formula is associated with a cost of unsatisfying the sub-formula. The objective is to determine an assignment that minimizes the sum of the costs of the unsatisfied sub-formulas. For example, consider the formula ϕ ≡ (x < 10) ∧ (y > 20) ∧ (x − y > 40) ∧ (y < 30), where subformulas (x − y > 40) and (y < 30) have costs 3 and 10, respectively. Note that not all sub-formulas need to have a cost. The sub-formulas without a cost have to be satisfied, and are called hard sub-formulas. The sub-formulas with a cost (which may or may not be satisfied) are called soft sub-formulas. One MAX-SMT solution for ϕ is x = 5, y = 25 with cost 3.
In the mSMQC model, all the Domain Restriction constraints,
Ultrametric constraints and Quartet Consistency constraints have to be satisfied
(hard formulas). What may or may not be satisfied is that the variable
that encodes the consistency of the topology has value 1. Thus, the soft
sub-formula for the t-th topology is:
(qt = 1)
(
This section analyzes the results obtained from running the existing ASP solution [18], the PBO models described in Section 3, and the SMT-based solutions proposed in Section 4. First, the experimental setup is described, followed by the analysis of the experimental results for two sets of instances.
The first set of instances were obtained from [18] and correspond to a set of artificially generated instances. Each instance has 10 taxa and is associated with a percentage of modified topologies (to account for errors). The percentages considered were 1%, 5%, 10%, 15%, 20%, 25%, and 30%.
The second set of instances was generated from a real example with 52 taxa1. This example corresponds to a complete set of topologies inferred 1This example was kindly provided by Dr. Guohui Lin. 1% 0.37 3.43 5.20 14.04 6.89 2.00 0.78 61.77 4.35 1.04 1.00 2.66 2.30 4.98 4.91 1.96 2.14 34.48 1.51 10.77 0.06 from a set of CCV distances. In order to compare the different models, the instance was divided in 10 sets of 10 taxa each. For a considered set of taxa, all the topologies on that set of taxa were gathered to form a new instance. To introduce a higher number of errors, the same idea of modifying a number of random topologies was used. The same percentages were used.
The encoders/solvers used in the experiments can be divided in 3 categories, phy+ASP-solver, PBO-encoder+PBO-solver and SMT-based solver.
phy+ASP-solver. phy is an encoder of the MQC problem into ASP, which was obtained from [18]. The approach of using phy with an ASPsolver is currently one of the best approaches in the literature. The ASP encoder used in the experiments is derived from the encoder phy [18]. phy receives as arguments, the number of taxa, the maximum number of allowed quartet errors and the instance with the topologies. In the experiments the maximum number of allowed errors was set to the total number of topologies.
Two ASP-solvers were considered for the experiments, SMODELS [13] and CLASP [6]. Other ASP-solvers were considered, but these two were the only that allowed maximizing with phy. The maximization was done by setting the ASP-solvers to enumerate all the stable models.
acterized by first using one of the encoders described in Section 3 to obtain a pbo file. The pbo file is then given to a pbo-solver. Three different PBOsolvers were used in the experiments, the bsolo solver [9], the minisat+ solver [4] and the pueblo solver [12].
SMT-based solver. The models based in SMT described in Section 4 require interaction with the SMT-solver. The SMT-solver used in the exper0% 1.24 1.99 7.32 9.12 3.88 1.70 1.48 93.21 2.81 1.20 2.05 3.09 2.17 3.12 4.78 2.89 2.57 34.20 0.46 0.46 0.04 iments was the Yices solver [3]. Yices offers a C API that allows the interaction with the solver, where it is possible to ask the solver for a satisfying assignment to the variables when the solver declares that the smt-instance is satisfiable. Yices is also able to solve max-smt instances.
The results were obtained in a set of nine Intel Xeon 5160, 3GHz dual core servers with 4GB of RAM. The timeout was set to 3600 seconds. Table 1 presents the results for the first set of instances. The percentages of modified topologies are shown in the top of the columns, and the encoders and solvers used are indicated on the left column. In terms of timeouts, the only timeouts were for instances with 20% or 30%. For 20%, the pbobinary model has timed-out with bsolo and pueblo in one instance, and the pbo-bin-sc model with pueblo has timed-out in one instance. For 30%, the pbo-basic has timed-out seven times with bsolo and once with pueblo. The pbo-reduce with bsolo timeout nine times and pbo-unary has timed out five times with bsolo and once with pueblo.
FromTable 1 and omparing the results of phy with both the ASP-solvers, we can conclude that smodels is on average four times faster than clasp for these instances. Comparing the results among the PBO-models, minisat+ stands out as being the fastest in these instances. The difference among the PB-encoders is not as clear. For instances with percentage of 5% or smaller, the pbo-reduced is the fastest model. On the other hand on instances with percentage of 10% or higher, then the binary models are faster. For the SMT-based approaches, the max-smt approach is the most efficient.
By analyzing all the models and solvers, we can conclude that for instances with percentage of 1% smt max-smt is faster than any other solver. For instances with percentage of 5%, 10% or 15% the phy + smodels approach has better performance. For instances with 20% or higher the pbobinary and pbo-binary-sc with minisat+ present better results.
Table 2 shows for the second set of instances the average running time of the instances per encoder/solver and percentage of modified topologies. The timeouts for these instances occurred for 1%, 25% and 30%. The pbo-binary model with bsolo timeout once for 1%. The pbo-basic with bsolo timed-out once for 25% and nine times for 30%. The pbo-unary with bsolo timeout once for 25% and seven times for 30%. The pbo-reduced with bsolo timeout nine times for 30%. For instances that have 0%, 1% or 5%, smt max-smt is the solver with best performance results, whereas if the instance has 10% of modified topologies, the pbo-reduced encoder with minisat+ is the fastest. For instances with 15% or more, the pbo-binary-sc encoder with minisat+ presents the best performance results. 6
This paper proposes new approaches for solving the Maximum Quartet Consistency (MQC) problem. Building on recent work on using Pseudo-Boolean Optimization (PBO) models for MQC [10], a number of new alternative PBO models is proposed. In addition, the paper shows how maximization variants of Satisfiability Modulo Theories (SMT) can be used for solving the MQC problem. To our best knowledge, this is the first concrete example of using SMT for solving computational problems in bioinformatics.
The experimental results suggest that the SMT models are the most adequate for instances with a small percentage of modified quartet topologies, whereas for instances with a larger percentage of modified quartet topologies the most effective solution is currently based on PBO using the Minisat+ solver. The results in the paper provide evidence that there are a number of effective alternative approaches for solving the MQC problem, including ASP, PBO and SMT.
Future research directions will include better integration of SMT solvers in addressing the MQC problem.
This work was supported in part by Microsoft Research through its PhD Scholarship Programme.