K S Sesh Kumar

Francis Bach
email: francis.bach@inria.fr

K S Sesh Kumar
email: sesh-kumar.karri@inria.fr

Convex Relaxations for Learning Bounded Treewidth Decomposable Graphs

We consider the problem of learning the structure of undirected graphical models with bounded treewidth, within the maximum likelihood framework. This is an NP-hard problem and most approaches consider local search techniques. In this paper, we pose it as a combinatorial optimization problem, which is then relaxed to a convex optimization problem that involves searching over the forest and hyperforest polytopes with special structures, independently. A supergradient method is used to solve the dual problem, with a run-time complexity of O(k 3 n k+2 log n) for each iteration, where n is the number of variables and k is a bound on the treewidth. We compare our approach to state-of-the-art methods on synthetic datasets and classical benchmarks, showing the gains of the novel convex approach.

Introduction

Graphical models provide a versatile set of tools for probabilistic modeling of large collections of interdependent variables. They are defined by graphs that encode the conditional independences among the random variables, together with potential functions or conditional probability distributions that encode the specific local interactions leading to globally well-defined probability distributions [START_REF] Bishop | Pattern recognition and machine learning[END_REF][START_REF] Wainwright | Graphical models, exponential families, and variational inference. Found. and Trends in[END_REF][START_REF] Koller | Probabilistic graphical models: principles and techniques[END_REF].

In many domains such as computer vision, natural language processing or bioinformatics, the structure of the graph follows naturally from the constraints of the problem at hand. In other situations, it might be desirable to estimate this structure from a set of observations. It allows (a) a statistical fit of rich probability distributions that can be considered for further use, and (b) discovery of structural relationship between different variables. In the former case, distributions with tractable inference are often desirable, i.e., inference with run-time complexity does not scale exponentially in the number of variables in the model. The simplest constraint to ensure tractability is to impose tree-structured graphs [START_REF] Chow | Approximating discrete probability distributions with dependence trees[END_REF]. However, these distributions are not rich enough, and following earlier work [START_REF] Malvestuto | Approximating discrete probability distributions with decomposable models[END_REF][START_REF] Bach | Thin junction trees[END_REF][START_REF] Narasimhan | PAC-learning bounded tree-width graphical models[END_REF][START_REF] Chechetka | Efficient principled learning of thin junction trees[END_REF][START_REF] Gogate | Learning efficient Markov networks[END_REF][START_REF] Szántai | Discovering a junction tree behind a Markov network by a greedy algorithm[END_REF], we consider models with bounded treewidth, not simply by one (i.e., trees), but by a small constant k.

Beyond the possibility of fitting tractable distributions (for which probabilistic inference has linear complexity in the number of variables), learning bounded-treewidth graphical models is a key to 1 design approximate inference algorithms for graphs with higher treewidth. Indeed, as shown by [START_REF] Saul | Exploiting tractable substructures in intractable networks[END_REF][START_REF] Wainwright | Graphical models, exponential families, and variational inference. Found. and Trends in[END_REF][START_REF] Kolmogorov | Generalized sequential tree-reweighted message passing[END_REF], approximating general distributions by tractable distributions is a common tool in variational inference. However, in practice, the complexity of variational distributions is often limited to trees (i.e., k = 1), since these are the only ones with exact polynomial-time structure learning algorithms. The convex relaxation designed in this paper enables us to augment the applicability of variational inference, by allowing a finer trade-off between run-time complexity and approximation quality.

Apart from trees, learning the structure of a directed or undirected graphical model, with or without constraints on the treewidth, remains a hard problem. Two types of algorithms have emerged, based on the two equivalent definitions of graphical models: (a) by testing conditional independence relationships [START_REF] Spirtes | Causation, prediction, and search[END_REF] or (b) by maximizing the log-likelihood of the data using the factorized form of the distribution [START_REF] Friedman | Being Bayesian about network structure. a bayesian approach to structure discovery in bayesian networks[END_REF]. In the specific context of learning bounded-treewidth graphical models, the latter approach has been shown to be NP-hard [START_REF] Srebro | Maximum likelihood bounded tree-width Markov networks[END_REF] and led to various approximate algorithms based on local search techniques [START_REF] Malvestuto | Approximating discrete probability distributions with decomposable models[END_REF][START_REF] Deshpande | Efficient stepwise selection in decomposable models[END_REF][START_REF] Karger | Learning Markov networks: maximum bounded tree-width graphs[END_REF][START_REF] Bach | Thin junction trees[END_REF][START_REF] Shahaf | Learning thin junction trees via graph cuts[END_REF][START_REF] Szántai | Discovering a junction tree behind a Markov network by a greedy algorithm[END_REF] while the former approach led to algorithms based on independence tests [START_REF] Narasimhan | PAC-learning bounded tree-width graphical models[END_REF][START_REF] Chechetka | Efficient principled learning of thin junction trees[END_REF][START_REF] Gogate | Learning efficient Markov networks[END_REF], which have recovery guarantees when the data-generating distribution has low treewidth. Malvestuto [START_REF] Malvestuto | Approximating discrete probability distributions with decomposable models[END_REF] proposed a greedy heuristic of hyperedge selection with least incremental entropy. Deshpande et al. [START_REF] Deshpande | Efficient stepwise selection in decomposable models[END_REF] proposed a simple edge selection technique that maintains decomposability of the graph while minimizing the KL-divergence to the original distribution. Karger et al. [START_REF] Karger | Learning Markov networks: maximum bounded tree-width graphs[END_REF] proposed the first convex optimization approach to learn the maximum weighted k-windmill, a sub-class of the decomposable graph. Bach et al. [START_REF] Bach | Thin junction trees[END_REF] gave an approach which iteratively refines the hyperedge selection based on KL-divergence using iterative scaling. Shahaf et al. [START_REF] Shahaf | Learning thin junction trees via graph cuts[END_REF] proposed another convex optimization approach with Bethe approximation of the likelihood using graph-cuts. Szántai et al. [START_REF] Szántai | Discovering a junction tree behind a Markov network by a greedy algorithm[END_REF] proposed a hyperedge selection criteria based on high mutual information within a hyperedge. Narasimhan et al. [START_REF] Narasimhan | PAC-learning bounded tree-width graphical models[END_REF] performs independence tests by solving submodular optimization problems and derives a decomposable graph using dynamic programming. Chechetka et al. [START_REF] Chechetka | Efficient principled learning of thin junction trees[END_REF] used the weaker notion of conditional mutual information instead of conditional independence to learn approximate junction trees. Gogate et al. [START_REF] Gogate | Learning efficient Markov networks[END_REF] uses low mutual information criteria to recursively split the state space to smaller subsets until no further splits are possible.

In this paper, we make the following contributions:

• We provide a novel convex relaxation for learning bounded-treewidth decomposable graphical models from data in polynomial time. This is achieved by posing the problem as a combinatorial optimization problem in Section 2, which is relaxed to a convex optimization problem that involves the graphic and hypergraphic matroids, as shown in Section 3.

• We show in Section 4 how a supergradient ascent method may be used to solve the dual optimization problem, using greedy algorithms as inner loops on the two matroids. Each iteration has a run-time complexity of O(k 3 n k+2 log n), where n is the number of variables. We also show how to round the obtained fractional solution.

• We compare our approach to state-of-the-art methods on synthetic datasets and classical benchmarks in Section 5, showing the gains of the novel convex approach.

Maximum Likelihood Decomposable Graphical Models

In this section, we first review the relevant concepts of decomposable graphs and junction trees; for more details, see [START_REF] Bishop | Pattern recognition and machine learning[END_REF][START_REF] Wainwright | Graphical models, exponential families, and variational inference. Found. and Trends in[END_REF][START_REF] Koller | Probabilistic graphical models: principles and techniques[END_REF]. We then cast the problem of learning the maximum likelihood bounded treewidth graph as a combinatorial optimization problem.

Decomposable graphs and junction trees

We assume we are given an undirected graph G defined on the set of vertices V = {1, 2, . . . , n}.

Let C(G) denote the set of maximal cliques of G (which we will refer to as cliques). We consider n random variables X 1 , . . . , X n (referred to as X), associated with each vertex indexed by V . For simplicity, they are assumed to be discrete, but this is not a restriction as maximum likelihood will use only entropies that can be extended to differentiable entropies [START_REF] Cover | Elements of information theory[END_REF].

The distribution p(x) of X is said to factorize in the graph G, if and only if it factorizes as a product of potentials that depend only on the variables within maximal cliques.

A graph is said to be decomposable if it has a junction tree, i.e., a spanning tree whose vertices are maximal cliques of G (i.e., C(G) is the vertex set) such that:

• the junction tree connects only cliques that have a common element (clique tree property),

• for any vertex i ∈ V , the subgraph of cliques containing i is a tree (running intersection property).

Let T (G) denote the edges of the junction tree over the set of cliques C(G). When the graph G is decomposable, the distribution p(x) of X factorizes in G if and only if it may be written as

p G (x) = C∈C(G) p C (x C ) (C,D)∈T (G) p C∩D (x C∩D ) , (1)
where x is an instance in the domain of X, which we denote by X . p C (x C ) denotes the marginal distribution of random variables belonging to C ∈ C(G) and p C∩D (x C∩D ) denotes the marginal distribution of random variables belonging to the separator set C ∩ D, such that (C, D) ∈ T (G). See Figure 1. The treewidth of G is the maximal size of the cliques in G, minus one.

An alternative representation of decomposable graphs may be obtained by considering hypergraphs.

Hypergraphs are defined by a base set V and a set of hyperedges, i.e., subsets of V . A hypergraph is said to be acyclic if and only if the resulting graph obtained by connecting all elements within an hyperedge is decomposable. Unfortunately, the nice properties of acyclic graphs do not transfer to acyclic hypergraphs. Particulary, the matroid property, which allows exact greedy algorithms, does not hold. In Section 3.2, we will use a lesser known more general notion of acyclicity that will lead to exact greedy algorithms.

Maximum likelihood estimation

Given N observations x 1 , . . . , x N of X, we denote the corresponding empirical distribution of X by p

(x) = 1 N N i=1 δ(x = x i ).
Given the structure of a decomposable graph G, the maximum likelihood distribution that factorizes in G may be obtained by combining the marginal empirical distributions on all maximum cliques and their separators.

Let p(x) denote the empirical distribution and pG (x) denotes the projected distribution on a decomposable graph G. Estimating the maximum likelihood decomposable graph which best approximates p is equivalent to finding the graph, G, which minimizes the KL-divergence between the target distribution and the projected distribution, pG , defined by D(p||p G ).

D(p||p

G ) = x∈X p(x) log p(x) pG (x) ∝ x∈X -p(x) log pG (x) as p(x) is independent of G = x∈X -p(x) log C∈C(G) pC (x C ) (C,D)∈T (G) pC∩D (x C∩D ) from Eq. (1) = x∈X -p(x) log C∈C(G) pC (x C ) - x∈X -p(x) log (C,D)∈T (G) pC∩D (x C∩D ) = C∈C(G) x∈X -p(x) log pC (x C ) - (C,D)∈T (G) x∈X -p(x) log pC∩D (x C∩D ) = C∈C(G) xC∈XC -p C (x C ) log pC (x C ) - (C,D)∈T (G) xC∩D∈XC∩D -p C∩D (x C∩D ) log pC∩D (x C∩D ) = C∈C(G) H(C) - (C,D)∈T (G) H(C ∩ D) , (2)
where H(S) is the empirical entropy of the random variables indexed by the set S ⊆ V , defined by H(S) = xS {-p S (x S ) log pS (x S )}, and where the sum is taken over all possible values of x S .

Note that in this paper, we will not be using a traditional model selection term [START_REF] Friedman | Being Bayesian about network structure. a bayesian approach to structure discovery in bayesian networks[END_REF] as we will only consider models of low tree-width (with a bounded number of parameters).

Combinatorial optimization problem

We now consider the problem of learning a decomposable graph of treewidth less than k. We assume that we are given all entropies H(S) for subsets S of V of cardinality less than k + 1.

Since we do not add any model selection term, without loss of generality [START_REF] Szántai | Hypergraphs as a mean of discovering the dependence structure of a discrete multivariate probability distribution[END_REF], we restrict the search space to the space of maximal junction trees, i.e., junction trees with n -k maximal cliques of size k +1 and n-k -1 separator sets of size k between two cliques of size k +1. Our natural search spaces are thus characterized by D, the set of all subsets of size k + 1 of V , of cardinality n k+1 , and E, the set of all potential edges in a junction tree, i.e.,

E = {(C, D) ∈ D × D, C ∩ D = ∅, |C ∩ D| = k}. The cardinality of E is n k+2 . k+2 2
(number of subsets of size k + 2 times the number of possibility of excluding two elements to obtain a separator).

A decomposable graph will be represented by a clique selection function τ : D → {0, 1} and an edge selection function ρ : E → {0, 1} so that τ (C) = 1 if C is a maximal clique of the graph and ρ(C, D) = 1 if (C, D) is an edge in the junction tree. Both ρ and τ will be referred to as incidence functions or incidence vectors, when seen as elements of {0, 1} D and {0, 1} E . Thus, minimizing the problem defined in Eq. ( 2) is equivalent to minimizing,

P(τ, ρ) = C∈D H(C)τ (C) - (C,D)∈E H(C ∩ D)ρ(C, D), (3)
with the constraint that (τ, ρ) forms a decomposable graph.

At this time, we have merely reparameterized the problem with the clique and edge selection functions. We now consider a set of necessary and sufficient conditions for the pair to form a decomposable graph. Some are convex in (τ, ρ), while some are not. The latter ones will be relaxed in Section 3. From now on, we denote by 1 i∈C the indicator function for i ∈ C (i.e., it is equal to 1 if i ∈ C and zero otherwise).

• Covering V : Each vertex in V must be covered by atleast one of the selected cliques,

∀i ∈ V, C∈D 1 i∈C τ (C) ≥ 1. (4)
• Number of edges: Exactly n -k -1 edges from E must be selected,

(C,D)∈E ρ(C, D) = n -k -1. (5)
• Number of cliques: Exactly n -k cliques from D must be selected,

C∈D τ (C) = n -k. (6)
• Running intersection property: Every vertex, i ∈ V must induce a tree, i.e., the number of selected edges containing the vertex, i, must be equal to the number of selected cliques containing the vertex, i , minus one.

∀i ∈ V,

(C,D)∈E 1 i∈(C∩D) ρ(C, D) - C∈D 1 i∈C τ (C) + 1 = 0. (7)
• Edges between selected cliques: An edge in E is selected by ρ only if the cliques it is incident on is selected by τ .

∀C ∈ D, τ (C) = max D∈D, (C,D)∈E ρ(C, D). (8)
• Acyclicity of ρ: ρ selects edges in E such that they do not have loops, e.g., the blue lines in Figure 1

• Acyclicity of τ : τ selects the hyperedges of V in D such that they are acyclic, i.e., τ represents an acyclic hypergraph of (V, D).

The above constraints encode the classical definition of junction trees. Thus our combinatorial problem is exactly equivalent to minimizing P (τ, ρ) defined in Eq. ( 3), subject to the constraints in in Eq. ( 4), Eq. ( 5), Eq. ( 6), Eq. ( 7), Eq. ( 8), Eq. ( 9) and Eq. [START_REF] Frank | On decomposing a hypergraph into k connected subhypergraphs[END_REF]. Note that the constraint Eq. ( 10) that τ represents an acyclic hypergraph is implied by the other constraints.

Figure 2 shows clique and edge selections in blue which satisfy all these constraints and hence represent a decomposable graph. The clique and edge selections in red violates at least one of these constraints.

Convex Relaxation

We now provide a convex relaxation of the combinatorial problem defined in Section 2.3. The covering constraint in Eq. ( 4), the number of edges and the number of cliques constraints in Eq. ( 5) and Eq. ( 6) respectively, and the running intersection property in Eq. ( 7) are already convex in (τ, ρ).

The constraint in Eq. ( 8) that ∀C ∈ D, τ (C) = max D∈D, (C,D)∈E ρ(C, D) may be relaxed into:

• Edge constraint: selection of edges only if the both the incident cliques are selected, i.e.,

∀C ∈ D, ∀(C, D) ∈ E, ρ(C, D) ≤ τ (C). (11)
• Clique constraint: selection of a clique if at least an edge incident on it is selected, i.e.,

∀C ∈ D, τ (C) ≤ (C,D)∈E ρ(C, D). (12)
We now consider the two acyclicity constraints in Eq. ( 9) and Eq. (10).

Forest polytope

Given the graph (D, E), the forest polytope is the convex hull of all incidence vectors ρ of subforests of (D, E). Thus, it is exactly the convex hull of all ρ : E → {0, 1} such that ρ satisfies the constraint in Eq. ( 9). We may thus relax it into:

-Tree constraint:

ρ is in the forest polytope of (D, E).

While the new constraint in Eq. ( 13) forms a convex constraint, it is crucial that it may be dealt with empirically in polynomial time. This is made possible by the fact that one may maximize any linear function over that polytope. Indeed, for a weight function w : E × E → R, maximizing

(C,D)∈E w(C, D)ρ(C, D
) is exactly a maximum weight spanning forest problem, and its solution may be obtained by Kruskal's algorithm, i.e., (a) order all (potentially negative) weights w(C, D) and (b) greedily select edges (C, D), i.e., set ρ(C, D) = 1, with higher weights first, as long as they form a forest and as long as the weights are positive. When we add the restriction that the number of edges is fixed (in our case n -k -1), then the algorithm is stopped when exactly the desired number of edges is selected (whether the corresponding weights are positive or not). See, e.g., [START_REF] Schrijver | Combinatorial optimization: Polyhedra and efficiency[END_REF].

The polytope defined above may also be defined as the independence polytope of the graphic matroid, which is the traditional reason why the greedy algorithm is exact [START_REF] Schrijver | Combinatorial optimization: Polyhedra and efficiency[END_REF]. In the next section, we show how this can be extended to hypergraphs.

Hypergraphic matroid

Given the set of potential cliques D over V , we consider functions τ : D → {0, 1} that are equal to one when a clique is selected, and zero otherwise. Ideally, we would like to treat the acyclicity of the associated hypergraph in a similar way than for regular graphs. However, the set of acyclic subgraphs of the hypergraph defined from D does not form a matroid, and thus the polytope defined as the convex hull of all incidence vectors/functions of acyclic hypergraphs may be defined, but the greedy algorithm is not applicable. In order to define what is referred to as the hypergraphic matroid, one needs to relax the notion of acyclicity.

We now follow [START_REF] Lorea | Hypergraphes et matroides[END_REF][START_REF] Frank | On decomposing a hypergraph into k connected subhypergraphs[END_REF][START_REF] Fukunaga | Computing minimum multiway cuts in hypergraphs from hypertree packings[END_REF] and define a different notion of acyclicity for hypergraphs. An hypergraph (V, F ) is an hyperforest if and only if for all A ⊂ V , the number of hyperedges in F contained in A is less than |A| -1. A non-trivially equivalent definition is that we can select two elements in each hyperedege so that the graph with vertex set V and with edge set composed of these pairs is a forest.

Given an hypergraph with hyperedge set D, the set of sub-hypergraphs which are hyperforests forms a matroid. This implies that given a weight function on D, one may find the maximum weight hyperforest with a greedy algorithm that ranks all hyperedges and select them as long as they don't violate acyclicity (with the notion of acyclicity just defined and for which we exhibit a test below).

Testing acyclicity. Checking acyclicity of an hypergraph (V, F ) (which is needed for the greedy algorithm above) may be done by minimizing with respect to

A ⊂ V |A| - G∈F 1 G⊂A .
The hypergraph is an hyperforest if and only if the minimum is greater or equal to one. The minimization of this function may be cast a min-cut/max-flow problem as follows [START_REF] Fukunaga | Computing minimum multiway cuts in hypergraphs from hypertree packings[END_REF]:

-single source, single sink, one node per hyperedge in F , one node per vertex in V , -the source points towards each hyperedge with unit capacity, -each hyperedge points towards the vertices it contains, with infinite capacity, -each vertex points towards the sink, with unit capacity.

Link with decomposability. The hypergraph obtained from the maximal cliques of a decomposable graph can easily be seen to be an hyperforest. But the converse is not true.

Hyperforest polytope. We can now naturally define the hyperforest polytope as the convex hull of all incidence vectors of hyperforests. Thus the constraint in Eq. ( 10) may be relaxed into:

-Hyperforest constraint: τ is in the hyperforest polytope of (V, D).

(14)

Relaxed optimization problem

We can now formulate our combinatorial problem from the constraints in Eq. ( 4), Eq. ( 5), Eq. ( 6), Eq. ( 7), Eq. ( 11), Eq. ( 12), Eq. ( 13) and Eq. ( 14) as follows min P(τ, ρ) subject to

                               τ ∈ {0, 1} D , ρ ∈ {0, 1} E , ∀i ∈ V, C∈D 1 i∈C τ (C) ≥ 1, (C,D)∈E ρ(C, D) = n -k -1, C∈D τ (C) = n -k, ∀i ∈ V, (C,D)∈E 1 i∈(C∩D) ρ(C, D) -C∈D 1 i∈C τ (C) + 1 = 0, ∀C ∈ D, ∀(C, D) ∈ E, ρ(C, D) ≤ τ (C), ∀C ∈ D, τ (C) ≤ (C,D)∈E ρ(C, D),
ρ is in the forest polytope of (D, E), τ is in the hyperforest polytope of (V, D).

(15)

All constraints except the integrality constraints are convex. Let τ -relaxed primal be the partially relaxed primal optimization problem formed by relaxing only the integral constraint on τ in Eq. ( 15), i.e., replacing τ ∈ {0, 1} D by τ ∈ [0, 1] D . Note that this is not a convex problem due to the remaining integral constraint on ρ, but it remains equivalent to the original problem as the following proposition shows.

Proposition 1

The combinatorial problem in Eq. ( 15) and the τ -relaxed primal problem are equivalent.

Proof Let us assume (τ * , ρ * ) be a feasible solution for the relaxed primal with 0 < τ * (C) < 1 for some C ∈ D. The edge constraint in Eq. [START_REF] Friedman | Being Bayesian about network structure. a bayesian approach to structure discovery in bayesian networks[END_REF] ensures that there are no incident edges on C selected by ρ * (as ρ * is integral). This violates the clique constraint in Eq. [START_REF] Fukunaga | Computing minimum multiway cuts in hypergraphs from hypertree packings[END_REF]. Therefore, the feasible solutions of relaxed primal are integral. Hence the optimal solutions of the primal and the relaxed primal are identical.

The convex relaxation for the primal optimization problem formed by relaxing the integral constraint on both τ and ρ can now be defined as min P(τ, ρ) subject to

                               τ ∈ [0, 1] D , ρ ∈ [0, 1] E , ∀i ∈ V, C∈D 1 i∈C τ (C) ≥ 1, (C,D)∈E ρ(C, D) = n -k -1, C∈D τ (C) = n -k, ∀i ∈ V, (C,D)∈E 1 i∈(C∩D) ρ(C, D) -C∈D 1 i∈C τ (C) + 1 = 0, ∀C ∈ D, ∀(C, D) ∈ E, ρ(C, D) ≤ τ (C), ∀C ∈ D, τ (C) ≤ (C,D)∈E ρ(C, D),
ρ is in the forest polytope of (D, E), τ is in the hyperforest polytope of (V, D).

(16)

Solving the dual problem

We now show how the convex problem may be minimized in polynomial time. Among the constraints of our convex problem in Eq. ( 15), some are simple linear constraints, some are complex constraints depending on the forest and hyperforest polytopes defined in Section 3. We will define a dual optimization problem by introducing the least possible number of Lagrange multipliers (a.k.a. dual variables) [START_REF] Bertsekas | Nonlinear Programming[END_REF] so that the dual function (and a supergradient) may be computed and maximized efficiently. We introduce the following dual variables:

-Set cover constraints in Eq. ( 4): γ ∈ R V + . -Running intersection property in Eq. ( 7): µ ∈ R V . -Edge constraints in Eq. ( 11): λ ∈ R 2E

+ . -Clique constraints in Eq. ( 12):

η ∈ R D + .
Therefore, the dual variables are (γ, µ, λ, η). Let L(τ, ρ, γ, µ, λ, η) be the Lagrangian relating the primal and dual variables. It is derived from the primal cost function defined in Eq. ( 3) along with the covering constraint, running intersection property, the edge and the clique constraints defined in Eq. ( 4), Eq. ( 7), Eq. ( 11) and Eq. ( 12) respectively. The Lagrangian can be computed from the dual variables (γ, µ, λ, η) as follows:

L(τ, ρ, γ, µ, λ, η) = C∈D H(C)τ (C) - (C,D)∈E H(C ∩ D)ρ(C, D) + i∈V γ i 1 - C∈D 1 i∈C τ (C) + i∈V µ i (C,D)∈E 1 i∈(C∩D) ρ(C, D) - C∈D 1 i∈C τ (C) + 1 + C∈D (C,D)∈E λ CD ρ(C, D) -τ (C) + C∈D η C τ (C) - (C,D)∈E ρ(C, D) = C∈D H(C) - i∈C (µ i + γ i ) - (C,D)∈E λ CD + η C τ (C) - (C,D)∈E H(C ∩ D) - i∈(C∩D) µ i -λ CD -λ DC + η C + η D ρ(C, D) + i∈V (µ i + γ i ), (17)
with the following dual constraints on the Lagrange multipliers

∀i ∈ V, γ i ≥ 0, ∀C ∈ D, ∀(C, D) ∈ E, λ CD ≥ 0, ∀C ∈ D, η C ≥ 0. ( 18
)
We can now derive a dual optimization problem with Q(γ, µ, λ, η) represent the dual cost function, which can be derived from the Lagrangian in Eq. [START_REF] Kolmogorov | Generalized sequential tree-reweighted message passing[END_REF]. We use the the number of edges constraint, the number of cliques constraint, tree constraint and hyperforest constraint given by Eq. ( 5), Eq. ( 6), Eq. ( 13) and Eq. ( 14) respectively in deriving the dual as follows:

Q(γ, µ, λ, η) = inf τ (C)∈[0,1] D C∈D τ (C)=n-k τ ∈ hyperforest polytope of (V,D) H(C) - i∈C (µ i + γ i ) - (C,D)∈E λ CD + η C τ (C) - sup ρ∈[0,1] E (C,D)∈E ρ(C,D)=n-k-1 ρ∈forest polytope of (D,E) (C,D)∈E H(C ∩ D) - i∈(C∩D) µ i -λ CD -λ DC + η C + η D ρ(C, D) + i∈V (µ i + γ i ). ( 19
)
It is decomposed in three parts defined in Eq. ( 21), Eq. ( 22) and Eq. ( 23) respectively :

Q(γ, µ, λ, η) = q 1 (γ, µ, λ, η) + q 2 (γ, µ, λ, η) + q 3 (γ, µ, λ, η), (20)
where 21)

q 1 (γ, µ, λ, η) = inf τ (C)∈[0,1] D C∈D τ (C)=n-k τ ∈ hyperforest polytope of (V,D) C∈D H(C) - i∈C (µ i + γ i ) - (C,D)∈E λ CD + η C τ (C). (
q 2 (γ, µ, λ, η) = -sup ρ∈[0,1] E (C,D)∈E ρ(C,D)=n-k-1 ρ∈ forest polytope of (D,E) (C,D)∈E H(C ∩ D) - i∈(C∩D) µ i -λ CD -λ DC + η C + η D ρ(C, D). ( 22
)
q 3 (γ, µ, λ, η) = i∈V (µ i + γ i ). ( 23
)
Therefore, the dual optimization problem using the dual cost function defined in Eq. ( 19) and the dual constraints defined in Eq. ( 18) is given by max Q(γ, µ, λ, η) subject to

   ∀i ∈ V, γ i ≥ 0, ∀C ∈ D, ∀(C, D) ∈ E, λ CD ≥ 0, ∀C ∈ D, η C ≥ 0. ( 24
)
The dual functions q 1 (γ, µ, λ, η) and q 2 (γ, µ, λ, η) may be computed using the greedy algorithms defined in Section 3.1 and Section 3.2; q 1 can be evaluated in O(r log(r)), where r is the cardinality of the space of cliques, D, i.e., n k+1 and q 2 can be evaluated in O(m log(m)), where m is the cardinality of feasible edges, E, i.e., n k+2 . k+2 2 . This complexity is due to sorting the edges and hyperedges based on their weights. This leads to an overall complexity of O(k 3 n k+2 log n) per iteration of the projected supergradient method which we now present.

Projected supergradient ascent. The dual optimization problem defined by maximizing Q(γ, µ, λ, η) can be solved using the projected supergradient method. In each iteration t of the algorithm, the dual cost function, Q(γ t , µ t , λ t , η t ), is evaluated through estimation of q 1 and q 2 by solving Eq. ( 21) and

Algorithm 1 Projected Supergradient

Input: clique and edge entropies H, step-size constant a and number of iterations T Output: sequence of clique and edge selections over iterations (τ t , ρ t ) Initialize γ 0 = 0, µ 0 = 0, λ 0 = 0, η 0 = 0 for t = 0 to T do solve Eq. ( 21) and evaluate q 1 (γ t , µ t , λ t , η t ) to obtain τ t solve Eq. ( 22) and evaluate q 2 (γ t , µ t , λ t , η t ) to obtain ρ t update dual variables, (γ t+1 , µ t+1 , λ t+1 , η t+1 ) using supergradients and stepsize:

α t = a √ t γ t+1 i = γ t i + α t 1 -C∈D 1 i∈C τ t (C) + µ t+1 i = µ t i + α t (C,D)∈E 1 i∈(C∩D) ρ t (C, D) -C∈D 1 i∈C τ t (C) + 1 λ t+1 CD = λ t CD + α t τ t (C) -ρ t (C, D) + η t+1 C = η t C + α t (C,D)∈E ρ t (C, D) -τ t (C) + end for
Eq. ( 22) respectively. In the process of solving these equations, the corresponding primal variables (τ t , ρ t ) are also estimated and allows the computations of the supergradients of Q (i.e., opposites of subgradients of -Q) [START_REF] Bertsekas | Nonlinear Programming[END_REF]. As shown in Algorithm 1, a step is made toward the direction of the supergradient and projection onto the positive orthant is performed for dual variables that are constrained to be nonnegative. With step sizes α t proportional to 1/ √ t, this algorithm is known to converge to a dual optimal solution [START_REF] Nedić | Approximate primal solutions and rate analysis for dual subgradient methods[END_REF] at rate 1/ √ t. Moreover, the average of all visited primal variables, i.e., after t steps, (τ t , ρt ) = 1 t t u=0 (τ u , ρ u ) is known to be approximately primal-feasible (i.e., it satisfies all the linear constraints that were dualized up to a small constant that is also going to zero at rate 1/ √ t). The convergence to primal feasibility is illustrated in Figure 4(a), where, on one of the synthetic examples from Section 5, the different constraint violations. Note that these are not the number of each of these constraints violated but the maximum value by which they are violated. It can be observed that the constraint violations reduce to zero over iterations.

Proposition 2 If k = 1, the convex relaxation in Eq. ( 16) is equivalent to Eq. [START_REF] Karger | Learning Markov networks: maximum bounded tree-width graphs[END_REF].

Proof If k = 1, all the cliques in the clique space contain only 2 vertices, i.e., ∀C ∈ D, |C| = 2 and the number of elements in the feasible edges is only 1, i.e., ∀(C, D) ∈ E, |C ∩ D| = 1.

Solving the convex relaxation defined in Eq. ( 16) is equivalent to solving the dual defined in Eq. [START_REF] Saul | Exploiting tractable substructures in intractable networks[END_REF]. On solving the dual variables, the optimal dual solution is given by

∀i ∈ V, µ i = H({i}), ∀i ∈ V, γ i = 0, ∀C ∈ D, ∀(C, D) ∈ E, λ CD = 0, ∀C ∈ D, η C = 0, ( 25
)
where

H({i}) = -p i (x i ) log(p i (x i )).
The optimal solution to the dual problem is given by

Q * (γ, µ, λ, η) = inf τ (C)∈[0,1] D C∈D τ (C)=n-k τ ∈ hyperforest polytope of (V,D) C∈D H(C) - i∈C H({i}) τ (C) + i∈V H({i}) = inf τ (C)∈[0,1] D C∈D τ (C)=n-k τ ∈ hyperforest polytope of (V,D) -I(C).τ (C) + i∈V H({i}), (26)
where ∀C ∈ D, I(C) = i∈C H({i})-H(C), which defines the mutual information of the elements in the clique, i.e., an edge if k = 1. The constraints in Eq. ( 26) define a spanning tree polytope [START_REF] Schrijver | Combinatorial optimization: Polyhedra and efficiency[END_REF] and the optimal solution is a maximal information spanning tree, which is given by Chow-Liu trees [START_REF] Chow | Approximating discrete probability distributions with dependence trees[END_REF].

They also form the optimal solution to the non-convex primal optimization defined in Eq. ( 15). Approximate Greedy Primal Solution. We describe an algorithm to project from the average of a sequence of fractional primary infeasible solutions, estimated during the iterations of projective supergradient, to an integral primary feasible solution. "AddClique" adds all the edges of a clique to the adjacency matrix. "checkGraphDecomposability" checks if the maximal cardinality search is a perfect elimination ordering. For decomposable graphs the maximal cardinality search yields a perfect elimination ordering [START_REF] Golumbic | Algorithmic Graph Theory and Perfect Graphs[END_REF]. We refer to this as decomposability test in this paper. "getNumberConnectedComponents" gives the number of connected components in the graph using breadth-first search. Note that the projection only uses the average clique selection function, τm , to obtain the primary feasible solutions, τ m . The corresponding edge selection, ρ m , can be estimated from clique selection, τ m , by selecting the edges between consecutive cliques of the perfect sequence of selected cliques [START_REF] Lauritzen | Graphical Models[END_REF]. The time complexity of the projection algorithm is O(n k+2 ). This is due to decomposability test with run time complexity O(n k+1 ), that is performed on adding O(n) cliques.

Note that the entropy of any graph, G, is independent of the mean of the normal distribution, hence we consider only the covariance matrix.

We use the graph structures representing a chain junction tree as in Figure 3-(a) and a star junction tree as in Figure 3-(b) to analyze the performance of our algorithm for decomposable covariance matrices generated with different correlations.

Table 1 and Table 2 show the performance of our algorithm on these two graphs. Decomposable covariance matrices are generated as above with different values of the correlation parameter d (all averaged over ten different random covariance matrices). We show the difference between the cost function in Eq. ( 3) and the optimal entropy, i.e., the one of the actual structure represented by the covariance matrices. The differences in the table are multiplied by 10 3 for brevity.

The first column ∆Dual represents the optimal value of our convex relaxation (obtained from the dual function), while the second column ∆Dual r represents the optimal value by replacing the hyperforest constraint by the simply τ ∈ [0, 1] D . We can see from the two tables, that the two values are strictly negative (i.e., we indeed have a relaxation) and that the hyperforest constraint is key to obtaining tighter relaxations. Note that the associated solutions are only fractional.

The third column ∆Primal represents the cost function obtained by projection of the optimal fractional solution of the hyperforest constraint, using Approximate Greedy Primal Solution algorithm. The fourth column ∆Primal r represents the cost function obtained by projecting the optimal fractional solution of the hypercube constraint, i.e., the corresponding primal feasible solution related to ∆Dual r . They are compared to a simple greedy algorithm in the fifth column that sorts all mutual information and keep adding the cliques with largest mutual information as long as decomposability is maintained. Although the relaxation is not tight, our rounding scheme leads empirically to the optimal solution when the correlations are strong enough (i.e., large values of d) and outperform the simple greedy algorithm. Performance Comparison. We compare the quality of the graph structures learned by the proposed algorithm with the ones produced by Ordering Based Search (OBS) [START_REF] Teyssier | Ordering-based search: A simple and effective algorithm for learning bayesian networks[END_REF], the combinatorial optimization algorithm proposed by Karger and Srebro (Karger+Srebro) [START_REF] Karger | Learning Markov networks: maximum bounded tree-width graphs[END_REF], the Chow-Liu trees (Chow-Liu) [START_REF] Chow | Approximating discrete probability distributions with dependence trees[END_REF] and different variations of PAC-learning based algorithms (PAC-JT, PAC-JT+local) [START_REF] Chechetka | Efficient principled learning of thin junction trees[END_REF]. We use a real-world dataset, TRAFFIC [START_REF] Krause | Near-optimal nonmyopic value of information in graphical models[END_REF] and an artificial dataset, ALARM [START_REF] Beinlich | The ALARM Monitoring System: A Case Study with Two Probabilistic Inference Techniques for Belief Networks[END_REF] to compare the performances of these algorithms. This ALARM dataset was sampled from a known Bayesian network [START_REF] Beinlich | The ALARM Monitoring System: A Case Study with Two Probabilistic Inference Techniques for Belief Networks[END_REF] of 37 nodes, which has a treewidth equal to 4. We learn an approximate decomposable graph of treewidth 3. The TRAFFIC dataset is the traffic flow information every 5 minutes for a month at 8000 locations in California [START_REF] Krause | Near-optimal nonmyopic value of information in graphical models[END_REF]. The traffic flow information is collected at 32 locations in San Francisco Bay area and the values are discretized into 4 bins. We learn an approximate decomposable graph of treewidth 3 using our approach. Empirical entropies are computed from the generated samples of each data set and we infer the underlying structure from them using our algorithm. show the log-likelihoods of structures learnt using various algorithms on Traffic and Alarm datasets respectively. Note that the performance is better with higher values as we compare log-likelihoods. These figures illustrate the gains of the convex approach over the earlier non-convex approaches.

Conclusion and Future Work

In this paper, we have provided a convex relaxation to the problem of finding the maximum likelihood decomposable graph with bounded treewidth, with a polynomial-time optimization algorithm, which empirically outperforms previously proposed algorithms. We are currently exploring two avenues for improvements: (a) design sufficient conditions for tightness of our relaxation, following the recent literature on relaxation of variable selection problems [START_REF] Candès | Decoding by linear programming[END_REF], and (b) use heuristics to speed-up the algorithms to allow application to larger graphs.

Figure 1 :

 1 Figure 1: (a) A decomposable graph on the set of vertices V = {1, 2, 3, 4, 5, 6, 7, 8, 9} having treewidth 2.(b) A junction tree embedded on the decomposable graph representing the maximal cliques by blue dots and the separator sets by blue lines. (c) The corresponding junction tree representation of the decomposable graph with ovals representing the maximal cliques and the rectangles representing the corresponding separator set.

  -(b) cannot form loops, ρ represents a subforest of the graph (D, E).

Algorithm 2

 2 Approximate Greedy Primal Solution Input: primal infeasible sequence τ t for Algorithm 1, treewidth k, number of Vertices n, set of cliques D and integer m such that 0 < m ≤ T Output: approximate discrete primal feasible solution τ m after m iterations of Algorithm 1 Initialize Adjacency Matrix Adj = zeros(n, n), τm = 1 m m t=0 τ t and τ m = zeros(length(τ m )) order = Sorted indices in the descending order τm repeat Initialize decomposable = f alse, treewidth = 0, numConnectedComponents = 0, i = 1 update T estAdj = AddClique(Adj, D(order(i))) update [decomposable, treewidth] = checkGraphDecomposability(T estAdj) if decomposable = true and treewidth ≤ k then update Adj = T estAdj update τ m (Order(i)) = 1 end if [numConnectedComponents] = getNumberConnectedComponents(T estAdj) update i = i + 1 until decomposable = true, treewidth = k, numConnectedComponents = 1, i = length(order)

Figure 3 :

 3 Figure 3: Graph representing (a) chain junction tree, (b) star junction tree, with an embedded junction tree in green and its junction tree representation in blue.

Figure 4 :

 4 Figure 4: Left: (a) Upper bound of constraint violations for d=2 and a chain junction tree. Right: Log likelihood of the structures learnt using various algorithms on (b) TRAFFIC and (c) ALARM datasets with k = 3 except Chow-Liu (k = 1).

Figure 4 (

 4 b) and Figure4(c)

Table 1 :

 1 Performance on chain junction trees. See text for details.

	d	∆Dual	∆Dual r	∆Primal ∆Primal r ∆Greedy
	1 -0.7±0.1 -32.7±16.4 0.2±0.1	0.4±0.1	0.2±0.1
	2 -0.4±0.1 -23.4±9.6	0	0.3±0.2	0.5±0.2
	4 -1.1±0.1 -31.2±9.7	0	0.3±0.1	1.9±0.3
	8 -0.6±0.1 -23.9±9.8	0	0.2±0.1	7.9±0.3
	16 -1.9±0.2	-3.4±2.7	0	0	25.6±1.2
	32 -3.9±0.5	-3.2±0.3	0	0	57.3±1.5

Table 2 :

 2 Performance on star junction trees. See text for details.

	d	∆Dual	∆Dual r	∆Primal ∆Primal r ∆Greedy
	1 -0.8±0.1 -31.4±13.4 0.2±0.1	0.5±0.1	0.9±0.1
	2 -0.5±0.2 -26.6±13.3	0	0.4±0.1	0.4±0.3
	4 -0.3±0.0 -16.6±4.1	0	0.2±0.1	1.7±0.2
	8 -0.4±0.0 -16.0±9.6	0	0	6.9±0.3
	16 -1.2±0.5	-3.1±0.3	0	0	26.3±1.5
	32 -6.8±0.4	-8.5±1.2	0	0	58.3±1.9

Acknowledgements We acknowledge support from the European Research Council grant SIERRA (project 239993), and would like to thank Anton Chechetka, Carlos Guestrin, Nathan Srebro and Percy Liang for sharing code and datasets. We would also like to thank other members of the SIERRA and WILLOW project-teams for helpful discussions.

Experiments and Results

In this section, we show the performance of the proposed algorithm on synthetic datasets and classical benchmarks.

Decomposable covariance matrices. In order to easily generate controllable distributions with entropies which are easy to compute, we use several decomposable graphs and we consider a Gaussian vector with covariance matrix Σ, generated as follows:

• sample a matrix Z of dimensions n × d ′ with entries uniform in [0, 1] and consider the matrix

where Z is a random matrix of dimensions n × d ′ , I is the n-dimensional identity matrix and d is a parameter to determine the correlations between the nodes of the graph, which takes values in {0, d ′ }. In our experiments, we choose d ′ to be 128. We have tight correlations between the nodes with higher values of d.

• normalize Σ ′ to unit diagonal, and

• The normalized random positive definite covariance matrix, Σ ′ , is projected onto a decomposable graph G as follows:

where the operator [(Σ ′ X ) -1 ] n gives an n × n matrix whose columns and rows representing the set X ⊆ V are filled by (Σ ′ X ) -1 and the rest of the elements of the matrix are filled with zeroes. The matrix, Σ, thus generated represents the covariance matrix of a multivariate Gaussian on a decomposable graph, G.

The projection ensures the following relationship between the random positive definite matrix, Σ ′ and the projected covariance matrix Σ:

where A is the adjacency matrix of the decomposable graph G onto which Σ ′ was projected.

The entropy of a multivariate Gaussian with a covariance matrix, Σ, is given by 1 2 log(2πe) n |Σ|, where |Σ| denotes the determinant of the covariance matrix. However, for Gaussian distribution that is factored in G ∈ G:

where Σ X is the sub-matrix of the covariance matrix whose rows and columns belong to the set X ⊆ V . Therefore, for any multivariate decomposable Gaussian graphical model, G: