Wafik Abassi

Fethi Aloui
email: fethi.aloui@univ-valenciennes.fr

Sassi Ben Nasrallah
email: sassi.bennasrallah@enim.rnu.tn

Jack Legrand
email: jack.legrand@univ-nantes.fr

Larbi Labraga
email: larbi.labraga@univ-valenciennes.fr

USE OF THE POD AND THE COHERENT STRUCTURE DETECTION CRITERIA TO STUDY THE FLOW DYNAMICS DOWNSTREAM A CONFINED SQUARE OBSTACLE

Keywords: Lattice Boltzmann Method (LBM), Experimental measurements, PIV, Vortex shedding, Detection criteria, Confined cylinder wake, POD

Studying the concept of the vortex is an essential tool in the comprehension of fluid dynamics. Despite this, it is still very difficult to find a universal definition of a vortex. Several methods of detection and characterization of vortex structures has been developed and performed. Specifically tailored for PIV data, they present an important topic in modern experimental fluid mechanics. In fact, the performance of the method is related to its effectiveness in the velocity data analysis and to successfully detect and locate the vortices as well as calculate the characteristic vortex parameters.

In the present study, we explore the efficiency of different vortex detection algorithms such as the vorticity ,  2 and Q criteria by studying the vortices structures in the wake of a confined square obstacle.

Proper orthogonal decomposition (POD) of the velocity fields was used to extract the energetic contribution of the different instabilities modes.

These methods were firstly applied to an experimentally two-dimensional instantaneous velocity fields obtained by Particle Image Velocimetry (PIV) technique. Then, in a second step, we tested these criterions on a numerical velocity fields determined from Lattice Boltzmann simulations and compared to experimental results.

INTRODUCTION

The confined flow around bluff bodies has received considerable attention in recent years, due to its fundamental relevance. This type of geometry have a variety of engineering applications in aerodynamics, wind engineering, electronics cooling and flow metering devices, etc. Bluff body cross sections that are often employed are circular and rectangular especially, square. The flow details behind these geometries depend mainly on Reynolds number and blockage ratio. Flow past a square cylinder resembles flow past a circular cylinder as far as instabilities are concerned. However, the vortex formation region is significantly broader and longer for a square cylinder compared to the circular. At low Reynolds number, aspect ratio and end conditions play a significant role in determining the flow properties. The study on the effect of aspect ratio in the literature has been largely focused to a circular cylinder (Rehimi et al. [START_REF] Rehimi | Experimental investigation of a confined flow downstream of a circular cylinder centred between two parallel walls[END_REF], Paranthoën et al. [START_REF] Paranthoen | Characteristics of the near wake of a cylinder at low Reynolds numbers[END_REF]). These experimental works have studied the wake of the cylinder: the effect of the confinement and gap parameter on the different flow regimes and the reorgonisation of the cylinder wake. Sushanta et al. [START_REF] Sushanta | Experimental investigation of flow past a square cylinder at an angle of incidence[END_REF] studied the sensitivity of flow patterns in the wake of square cylinder cross section to aspect ratio and orientation with respect to the mean flow (the incidence angle). Numerical simulations are documented in the works of Breuer et al. and Turki et al. [START_REF] Breuer | Accurate computations of the laminar flow past a square cylinder based on two different methods: Lattice-Bolzamann and finite-volume[END_REF][START_REF] Turki | Modeling and numerical simulation of the flows and the transfers of heat around an obstacle[END_REF]. These works, especially based on 2D laminar flows simulations, have studied the effect of the Von Karman vortex in the wake behind the obstacle. In this paper, we have restricted our attention to the range of Reynolds were exist a laminar Von Karman vortex street. The appears of this coherent structure is conditioned by the These coherent structures play an essential role in heat and mass transfer. Their interaction and the mecanism of generation/dissipation constitues a fundamental topic in fluid mechanics research. Many definitions of the concept of coherent structures have been proposed [START_REF] Lugt | Vortex Flow in Nature and Technology[END_REF][START_REF] Chassaing | Turbulence en mécanique des fluides[END_REF] [START_REF] Lugt | Vortex Flow in Nature and Technology[END_REF][START_REF] Chassaing | Turbulence en mécanique des fluides[END_REF]. However, some definitions are still debated, and the most essential question with coherent structure includes: how do we detected the center of vortex and how do we define its limits? This search of the coherent structure definition has led to the development of the detection criterion. The present study reports experimental measurements of flow patterns in the wake downstream of a square cylinder localazed in a rectangular cross-section and centred in a channel, at low Reynold numbers. The main objective of this work is to study the influence of confinement on different regimes and the reorganisation of the flow in the cylinder wake. For this purpose, in the first step, the vorticity , the Q and  2 criteria have been applied to the PIV velocity fields. It is also essential to add that such configuartion is suitable for the validation of numerical simulations codes. Thus, the same criteria are applied, in the second step to a numerical velocity fields obtained by the Lattice Boltzmann Method simulation [START_REF] Abassi | Lattice Boltzmann method used to simulate an unsteady flow around an obstacle in laminar regime[END_REF]. Proper orthogonal decomposition (POD) of the flow was used for filtring and extracting the energetic contribution of different modes. For the experimental investigation, a square obstacle of 10 mm height was placed in a channel with a rectangular cross section 300x30 mm². The confinement ratio was fixed in this study to r=1/3 and the gap to =1.

GENERALIZED LATTICE BOLTZMANN METHOD

Lattice Boltzmann method is a recent numerical method alternative to the classical techniques such as finite-volume, finite-element and finite-difference for solving Navier-Stokes equations.

This method treats the fluid on a statistical level, simulating the movement and interaction of single particles by solving a discrete velocity Boltzmann equation on a particle distribution function noted f(x,c,t).

With the BGK approximation and single relaxation time, this equation can be written as:

(1)

The discretized form is given by:

(2)

The right term of Eq. ( 2) represents the collision process through a relaxation time , while the left one is the streaming process. Many authors studied different approaches to write the collision process and the simplest one in order to obtain a linearized form with a single relaxation time. This form, was introduced by Bhatnagar, Gross and Krook (BGK) [START_REF] Bhatnagar | A model for collision process in gases[END_REF]. The common notation used for Lattice pattern is refered to the dimension of the problem and the number of speed DnQm where n represents the dimension of the problem and m the speed models.

In 2D problems, there are generally two types of Lattice patterns: square lattice and hexagonal lattice suggested in the literature. Generally, for LBM computation, an equidistant orthogonal Lattice is chosen and on each lattice node, the distribution function f is stored and updated on each new time step. As usual, the 9-velocity square Lattice is the most widely used. In fact, numerical studies proved that both D2Q9 and 7-velocity hexagonal lattice D2Q7 have sufficient Lattice symmetry, which is a dominant requirement for recovery of the correct flow equations [START_REF] Frish | Lattice-gas automata for Navier-Stokes equation[END_REF] and recent studies have shown that D2Q9. It usually gives more accurate results than that based on hexagonal Lattice [START_REF] Skordos | Initial and boundary conditions for the Lattice Boltzmann method[END_REF]. Hence, we have used in this study the 9-velocity square lattice D2Q9, on which each particle moves one Lattice unit only along the eight links indicated with 1-8. The 0 indicates the rest particles with zero speed and the discrete velocity vector of particles is defined by:

(3) where is the Lattice speed, eq f is the equilibrium distribution functions (EDFs) that plays an essential role in the Lattice Boltzmann Method and recovering the Navier-Stokes equations [START_REF] Bhatnagar | A model for collision process in gases[END_REF].

This equilibrium distribution function can be considered as a Gaussian quadrature of the Maxwell-Boltzmann continuous distribution function, which is expanded as a Taylor series in macroscopic velocity to its second order to recover incompressible Navier-Stokes equation. According to He et al. [START_REF] He | Theory of lattice Boltzmann method: from the Boltzmann equation to lattice Boltzmann equation[END_REF], it is given by the following form:

  2 2 4 2 1 22 i eq i i j i s s s c .u c .u u² f ( x ,t ) c c c            (4)
This local EDF is computed at each time step for every node from the local flow velocity, the fluid density  , the Lattice geometry weighting factor i w and the sound velocity c s .

In this equation, i w represent the weight coefficients and c s is called the sound velocity in LBM, which is related to c and i w .

For D2Q9 model, we have [START_REF] Abassi | Lattice Boltzmann method used to simulate an unsteady flow around an obstacle in laminar regime[END_REF]: i w =1/9 for the principal directions: i=1,2,3,4 i w =1/36 for diagonals: i=5,6,7,8 and 0 w = 4/9. At each simulation time step, the density  and the velocity field u are respectively defined in terms of the particle distribution function as follows:

(

The kinematic viscosity of the fluid is given by:

2 1 2 s x        (7)
The great progress of the LB Method during the last decades is due to its many advantages for fluid dynamics simulation.

In fact, only a loop of two steps is implemented during programming: the collision and streaming step which can recover the nonlinear macroscopic advection terms. Also, LBM is a linear equation, where nonlinearity is implicitly imbedded in the collision operator. For incompressible flow, the pressure satisfies a simple state equation unlike the common CFD method, solving Poisson equation becomes unnecessary [START_REF] Hunt | Stream and convergence zones in Turbulent flows[END_REF].

The pressure is then defined by the following state equation as:

(8)

EXPERIMENTAL SET UP AND MEASUREMENT TECHNIQUES

The experimental device is a hydraulic channel made up of a transparent Plexiglas (figure. 1). The length of the testsection is L=3 m, the width is l=0.3 m and H=0.03 m in height. A square obstacle (of dimensions 0.01x0.01 m²) is placed in the test-section perpendicularly to the mean flow direction at a distance l o =1.6m from the entry. The gap is γ=Δ/D=1, were Δ is the distance between the square and the nearest wall. The confinement ratio r=D/H is fixed at r=1/3. The mean flow velocity flow in the channel section upstream the obstacle was in the range 0.2 to 1.8 cm/s.

In the exit of the test-section, two Brooks flowmeters (range from 0 to 0.530 m 3 /h and from 0 to 1.17 m 3 /h) are assembled in parallel in order to measure the liquid flowrate for a range of Reynolds numbers included between 20 < Re < 400. The liquid is then recovered in a storage tank of 0.5 m 3 of capacity. A pump allows filling a tank of 0.15 m 3 capacity, where the water level is maintained constant by an overflow.

A PIV system has been used to determine instantaneous velocity fields downstream the square. It was composed of a CCD camera of 1600x1184 pixels resolution (Dantec Dynamics Flow sense), a pulsed ND-Yag laser (New wave solo). The whole system is driven by the "DANTEC" software "Flow Manager". The flow seeding was performed using polyamide particles with 20 m of diameter.

Figure 1. Schematic view of the used hydraulic channel Figure 2. Physical model

A 45° mirror was used to reflect the horizontal laser plane vertically at z = 0. The physical field of view of camera is 81.73x 61 mm². A mask was applied to retain only 30 mm in the y-direction. For the measurements, the interrogation's zones were taken equal to 32x32 pixels. An adaptive cross correlation was used with an overlap of 50%. The sampling frequency was fe=15Hz for a duration of 20 to 40s.

COHERENT STRUCTURES DETECTION

Many criteria used in practice are mainly deduced from the velocity gradient tensor. The objective of these criteria is to characterise the full coherent structure.

The vorticity tensor, as an invariant Galilean, appears as the most widely choice to identify the coherent structures in the flow. Hussain et al. [START_REF] Hussain | Coherent Structures and Turbulence[END_REF] proposed that the coherent structure is a turbulent zone in which, instantaneous vorticity is correlated with the phase of the flow. The vector   is used to compute the rotation rate in the flow. It is defined as

  U       2
1
where   is the rotationel operator.

Strawn et al. [START_REF] Strawn | Computer visualisation of vortex wake systems[END_REF] defined the center of the vortex as a local maximum of the vorticity module written for a 2D flow as follow:

              y U x U x y 2 1  (9)
Jeong and Hussain [START_REF] Jeong | On the Identification of a vortex[END_REF] studied the use of vorticity  and show that the shear and the rotation of the vortex are not distinguished by this criterion. In our case, the presence of shear at the upper and lower walls of the channel deforms the vorticity pattern of the coherent structures. Based an this observation, another criteria is proposed. Known under the name of the Q criterion, it was proposed by Hunt et al. [START_REF] Hunt | Stream and convergence zones in Turbulent flows[END_REF]. It can be interpretd as a balance between the rotation rate and the strain rate and can not be affected by the local shear. For a 2D flow, this Q creterion is defined as follow:

                                                    2 2 . 2 2 1 y U x U y U x U Q x y
x y [START_REF] Strawn | Computer visualisation of vortex wake systems[END_REF] The vortex structures are identified by the positive iso-values of Q, while their centres are identified by the maximum values of this criterion.

Graftieaux et al. [START_REF] Graftieaux | Combining PIV, POD and vortex identification algorithms for the study of unsteady turbulent swirling flows[END_REF] developed a kinematic criterion noted  1 for the vortex center identification considering only the topology of the flow. Although it is not a Galilean invariant, this criterion has presented a simple and robust method for the vortex structures identification. Thereafter this criteria was improved to became a Galilean invariant. This new criterion named  2 has the advantage to consider the local vortex convection. The extremum is used for the vortex identification. This criterion gives accurate results in comparison with the vorticity magnitude or the Q criterion and it is defined as:

    dS P U M U PM e P U M U PM S S M z . ) ( ) ( . . ) ( ) ( 1 2        (11)
where

dS M U S P U S M . ) ( 1 ) ( ~  
, e z is a unit vector perpendicular to the velocity map, P is the point where the criterion is calculated and S is the surface arrounding the point P. S should be as small as possible since for a large S, a small scale is filtred. For our case the size of S is: S=8.x  y.

PROPER ORTHOGONAL DECOMPOSITON OF THE FLOW (POD)

The velocity fields obtained from PIV measurements are generally noised. By applying vorticity  or velocity gradient tensor based on the velocity derivatives, high noise frequencies are amplified.

Using filter is not adequate in many cases. For these reasons, the POD applied to the flow velocity fields seems as a good solution. It takes into consideration the energetic distribution of the flow for the acquisition duration.

Proper Orthogonal Decomposition (POD) brings new opportunities because this post-processing extracts spatial or temporal structures using a rigorous mathematical decomposition basis, which breaks free from all arbitrary steps. It became popular for aerodynamics research in the 2000s, starting with Tang et al. [START_REF] Tang | System identification and proper orthogonal decomposition method applied to unsteady aerodynamics[END_REF] works, although it was first proposed in 1960 by Lumley [START_REF] Lumley | Atmospheric turbulence and radio wave propagation[END_REF] in the context of coherent structure's investigations.

The basic purpose of POD is to extract the most probabilistic structures in an energetic point of view from a statistical series of signals. Its principle is the creation of a mathematical model that decouples the spatial from the temporal variations in the flow.

This modal base needs only few modes to describe all principal aspects of the signal (majority of the large scale behaviour).

Sirovich [START_REF] Sirovich | Turbulence and dynamics of coherent structures part1: coherent structures[END_REF] has introduced the method of snapshots as a way to reduce the computational requirements of the POD calculations, especially when the spatial data size is much larger than the number of images like PIV measurements. For these reasons, we used snapshot in this work.

Considering the measurements obtained by PIV, it yields, generally, N snapshots of a 2D section of the flow field taken at the sequence time t 1 ,…,t N . These snapshots will usually give feature information on the velocity vectors [U(x,y,t

m n n N n j m j n t a t a , 1 ) ( ) ( . ) ( ). (      (14)       m n if m n if m n 0 1 ,



where n  is the eigenvalue associated to the spatial mode ) , (

) ( y x n  .
The modes verify the relationship:

m n D m n m n dy dx y x y x , ) ( ) ( ) ( ) ( . . ) , ( ). , ( ,         (15)
To filter the intantaneous PIV velocity fields, the equation ( 13) is truncated at a reduced order model M << N, which is chosen as a compromise between model simplicity and model accuracy. The filterd velocity is written as followed:

     M n n k n k filt y x t a y x U t z y x U 1 ) ( ) ( . ) , ( ). ( ) , ( ) , , , ( (16)

RESULTS AND DISCUSSIONS

Flow modes behind the obstacle

An example of a velocity profile obtained by 2D PIV measurements upstream of the square obstacle is shown in figure 3. In order to verify that the studied flow is developed and symmetrical in the channel, we compare the obtained profile with the theoretically solution of Lundgren and Sparrow [START_REF] Lundgren | Pressure drop due to the entrance region in ducts of arbitrary cross section[END_REF]. The measurements are in good agreement with the theoretical results.

To characterize the flow and the various coherent structures present in the flow, we used the PIV measurements obtained by varying the Reynolds number. In fact, for a fixed gap and confinment ratio, only the Reynolds number controls the flow regime. In our case, it is defined as

 / . Re D U m 
where m U is the average axial velocity in the inlet channel flow and  the kinematic viscosity of the flow. The regimes downstream of the square obstacle, similar to those of an obstacle placed in an infinite medium appears by varying the Reynolds number. These schemes are also comparable to those obtained numerically by various authors cited in the literature. In the second part, the velocity profiles were compared with numerical ones obtained by the Lattice boltzmann Method simulation. An example of experimental velocity profile, superposed to the numerical one obtained by LBM simulation, is illustrated by the figure 5.

For the Reynolds number Re=80 (above the critical regime), and for a longitudinal position close to the obstacle x/D=2, the axial and radial velocity profiles are in accordance with the numerical results. The velocity profile of U x according to the y/D position, presents two peaks on the higher and the lower part of the obstacle (y=±D), and then converges towards zeros on the two wall sides. These peaks show well the presence of two vortices on the same level of the x value.

Coherent structure detection Criteria

In this section, different criteria were applied at first to the PIV velocity fields, then to the numerical velocity fields obtained by the LBM simulation. For this reason, we fixed a specific Reynolds number Re=90 for a periodic regime, which is characterized by an instable wake and the appearance of the von Karman vortices. The streamlines and the axial velocity magnitude contours (figure 6) show clearly the oscillatory form of the wake. However, the von Karman vortices are not well highlighted as they are embedded in the mean flow. The vorticity fields brings up two families vortex well detectable (figure 7). Using the vorticity  as a creterion for the structure's detection, is not effective in the presence of shearing in the channel near the wall. Hence, we adopted Q and Γ 2 criteria in order to correctly identify the vortices. Figure 8 shows the Q criterion applied to the PIV velocity fields, at first, then to the numerical (LBM) velocity fields. This criterion brings up the structures that we can define their limits by the nearest values around the vortex centre, which is represented by the maximum of Q values. The comparaison between the two results shows a good agreement.

The Γ 2 criterion applied to PIV and numerical velocity fields (synchronized at the same time), shows the presence of two types of vortices (figure 9). A first family (P 1 , P 2 ) mainly composed of two counter-rotating vortices similar to those of von Kármán and alternatively detached behind the obstacle. Containment, created by the channel walls, leads to a significant change in the vortex trajectory as well as in the structure of the wake.

For example, the passage of the detached vortex on the upper side (P 1 ), creates a change in the space properties near the top wall. This latter, opposes it by creating a vortex with inverse direction, (P 1 ') which is ejected by forcing (P 1 ) through the channel along the y-axis towards the bottom wall. The same phenomenon occurs for detached vortex (P 2 ) on the lower side and the vortex (P 2 ') is ejected towards the higher wall. This proves that a vortex tends to roll and not slide on a wall. It should be noted that the distribution of these vortices is different from the infinite medium case. In fact, the vortex (P 1 ) detached from the upper wall does not reach directly the low wall of the channel, but it encounters, first, a vortex (P 2 ), which creates a coalescence between the vortices (at x=6D for this case). It is to be noted also that the confinement reduces the expansion of the vortices in the transverse direction of the flow.

An example of  2 profile at a specific zone of the channel obtained from the PIV measurements, is shown in figure 10. It is in good accordance with result obtained when the  2 criterion is applied to velocity fields obtained by the LBM simulation.

POD DECOMPOSITION OF THE FLOW

The POD was applied to the experimental PIV measurements obtained for a number of snapshots N=300. It was verified that number of snapshots N=300, is largely sufficient for the statistcal convergence of the data for a given Reynolds number. The four first modes are directly related to the von Karman structures and are the most energetic (figure 11). The reconstruction of filtered data on the basis of this modes permits to recuperate more than 90% of the fluctuating flow energy as % 90

4 1                 i N i i i   .
It should be noted that the POD was applied to two positions separetly. For this reason, the issued spatial modes are interpreted independently. The spatial modes

CONCLUSIONS

In this work, we presented a study of the wake behind a square obstacle in a confined environment. Experimental results of velocity fields issued from PIV measurements were used and compared with numerical ones obtained by Lattice Boltzmann simulations. The comparison between results shows a good accordance. The presence of the square obstacle in channel creates, from the Re > Re c the von Kármán vortices interacting with the walls to form secondary vortices advected downstream the obstacle. These results are important to show the fundamental difference with von Kármán vortices generated in unconfined medium. Thus, we applied, firstly, the structure criteria detection (, Q and Γ 2 ) to the instantaneous experimental and numerical velocity fields. The objective is to properly extract and describe the dynamics of vortex structures present in the flow, and which play a significant role in transfer and transport phenomena. The obtained results exhibit the interactions between a vortex and a wall, which are characterized by the creation of a secondary one at the wall with opposite vorticity signs. This latter is advected to the downstream then coalesce with the opposite vortices. The decomposition of the velocity field that has been implemented in this research has shown that the relevant structures of the field can be described with a limited number of modes (four) though capturing more than 90% of the kinetic energy of the complete field (mean and fluctuating fields). Subdomain 1 (near the obstacle) has revealed the presence of spatial modes paired where the von Karaman vortex street is mainly intense. The LBM simulation results are validated with our experimental ones which highlights the accuracy of this LBGK method for simulation of flow Dynamics behind bluff bodies.

  is the height averged velocity, D the obstacle height and  the kinematic viscosity) and by two additional parameters due to the wall effect. The latter are the gap parameter D    ( is the distance between the obstacle and the nearest wall) and the confinment ratio H D r  (which H is the height of the channel).

Figure 3 .

 3 Figure 3. Velocity profiles in the upstream of the obstacle (z=0)

Figure 4 .

 4 Figure 4. Instantaneous Streamlines and axial velocity contours PIV measurements: (a) Re=50; (b) Re=80; (c) Re=120

Figure 5 .

 5 Figure 5. Experimental and numerical radial velocity profiles at x/D=2 for Re=80: (a) U x ; (b) U y

Figure 6 .Figure 7 .

 67 Figure 6. Instantaneous streamlines and velocity magnitude contours obtained from the LBM simulation for Re=90

Figure 8 .

 8 Figure 8. Q criterion for Re=90: (a) PIV velocity field; (b) LBM simulation velocity field

Figure 9 .

 9 Figure 9.  2 criterion for Re=90: (a) PIV velocity field; (b) LBM simulation velocity field

Figure 10 .

 10 Figure 10. Comparison between the  2 experimental and numerical LBM results in the horizontal centreline of the channel for Re=90



  are paired for Re=90 and present a spatial periodicity (figure12). The right columns of figure12show the spatial modes for the measurements taken in the position 2 far from of the obstacle. The spatial modes relative to the zone near the obstacle (position 1) are shown in the left column of figure12. The coherent mode can be clearly observed mainly in the position 1 where the von Karman vortex street is more intense.

Figure 11 .Figure 12 .

 1112 Figure 11. Distribution of total energy starting from the spatial mode 0

ACKNOWLEDGMENTS

This work was supported by the laboratory GEPEA (University of Nantes), the laboratory TEMPO (University of Valenciennces) and the "Agence Nationale de la Recherche (ANR), France" by the grant n° ANR-08-BLAN-0184-02. These supports are gratefully acknowledged.


Fulya Trösser

Simon De Givry
email: simon.de-givry@inrae.fr

George Katsirelos

Improved Acyclicity Reasoning for Bayesian Network Structure Learning with Constraint Programming

Bayesian networks are probabilistic graphical models with a wide range of application areas including gene regulatory networks inference, risk analysis and image processing. Learning the structure of a Bayesian network (BNSL) from discrete data is known to be an NP-hard task with a superexponential search space of directed acyclic graphs. In this work, we propose a new polynomial time algorithm for discovering a subset of all possible cluster cuts, a greedy algorithm for approximately solving the resulting linear program, and a generalised arc consistency algorithm for the acyclicity constraint. We embed these in the constraint programmingbased branch-and-bound solver CPBayes and show that, despite being suboptimal, they improve performance by orders of magnitude. The resulting solver also compares favourably with GOBNILP, a state-of-the-art solver for the BNSL problem which solves an NP-hard problem to discover each cut and solves the linear program exactly.

Introduction

Towards the goal of explainable AI, Bayesian networks offer a rich framework for probabilistic reasoning. Bayesian Network Structure Learning (BNSL) from discrete observations corresponds to finding a compact model which best explains the data. It defines an NP-hard problem with a superexponential search space of Directed Acyclic Graphs (DAG). Several constraint-based (exploiting local conditional independence tests) and score-based (exploiting a global objective formulation) BNSL methods have been developed in the past.

Complete methods for score-based BNSL include dynamic programming [Silander and Myllymäki, 2006], heuristic search [START_REF] Yuan | Learning optimal bayesian networks: A shortest path perspective[END_REF][START_REF]Xiannian Fan and Changhe Yuan. An improved lower bound for bayesian network structure learning[END_REF], maximum satisfiability [Berg et al., 2014], branch-andcut [START_REF] Bartlett | [END_REF] and constraint programming [ [START_REF] Van Beek | Machine learning of bayesian networks using constraint programming[END_REF]. Here, we focus on the latter two.

GOBNILP [START_REF] Bartlett | [END_REF] is a state-of-theart solver for BNSL. It implements branch-and-cut in an in-teger linear programming (ILP) solver. At each node of the branch-and-bound tree, it generates cuts that improve the linear relaxation. A major class of cuts generated by GOBNILP are cluster cuts, which identify sets of parent sets that cannot be used together in an acyclic graph. In order to find cluster cuts, GOBNILP solves an NP-hard subproblem created from the current optimal solution of the linear relaxation.

CPBayes [ [START_REF] Van Beek | Machine learning of bayesian networks using constraint programming[END_REF]] is a constraint programming-based (CP) method for BNSL. It uses a CP model that exploits symmetry and dominance relations present in the problem, subproblem caching, and a pattern database to compute lower bounds, adapted from heuristic search [Fan and [START_REF]Xiannian Fan and Changhe Yuan. An improved lower bound for bayesian network structure learning[END_REF]. van Beek and Hoffmann showed that CPBayes is competitive with GOBNILP in many instances. In contrast to GOBNILP, the inference mechanisms of CPBayes are very lightweight, which allows it to explore many orders of magnitude more nodes per time unit, even accounting for the fact that computing the pattern databases before search can sometimes consume considerable time. On the other hand, the lightweight pattern-based bounding mechanism can take into consideration only limited information about the current state of the search. Specifically, it can take into account the current total ordering implied by the DAG under construction, but no information that has been derived about the potential parent sets of each vertex, i.e., the current domains of parent set variables.

In this work, we derive a lower bound that is computationally cheaper than that computed by GOBNILP. We give in Section 3 a polynomial-time algorithm that discovers a class of cluster cuts that provably improve the linear relaxation. In Section 4, we give a greedy algorithm for solving the linear relaxation, inspired by similar algorithms for MaxSAT and Weighted Constraint Satisfaction Problems (WCSP). Finally, in Section 5 we give an algorithm that enforces generalised arc consistency on the acyclicity constraint, based on previous work by van Beek and Hoffmann, but with improved complexity and practical performance. In Section 6, we show that our implementation of these techniques in CPBayes leads to significantly improved performance, both in the size of the search tree explored and in runtime.

the reader to existing literature [Papadimitriou and Steiglitz, 1998;[START_REF] Rossi | [END_REF] for more.

Constraint Programming

A constraint satisfaction problem (CSP) is a tuple V, D, C , where V is a set of variables, D is a function mapping variables to domains and C is a set of constraints. An assignment A to V ⊆ V is a mapping from each v ∈ V to D(v). A complete assignment is an assignment to V . If an assignment maps v to a, we say it assigns v = a. A constraint is a pair S, P , where S ⊆ V is the scope of the constraint and P is a predicate over V ∈S D(V ) which accepts assignments to S that satisfy the constraint. For an assignment A to S ⊇ S, let A | S be the restriction of A to S. We say that A satisfies c = S, P if A| S satisfies c. A problem is satisfied by A if A satisfies all constraints.

For a constraint c = S, P and for v ∈ S, a ∈ D(v), v = a is generalized arc consistent (GAC) for c if there exists an assignment A that assigns v = a and satisfies c. If for all v ∈ S, a ∈ D(v), v = a is GAC for c, then c is GAC. If all constraints are GAC, the problem is GAC. A constraint is associated with an algorithm f c , called the propagator for c, that removes (or prunes) values from the domains of variables in S that are not GAC.

CSPs are typically solved by backtracking search, using propagators to reduce domains at each node and avoid parts of the search tree that are proved to not contain any solutions. Although CSPs are decision problems, the technology can be used to solve optimization problems like BNSL by, for example, using branch-and-bound and embedding the bounding part in a propagator. This is the approach used by CPBayes.

Integer Linear Programming

A linear program (LP) is the problem of finding

min{c T x | x ∈ R n ∧ Ax ≥ b ∧ x ≥ 0}
where c and b are vectors, A is a matrix, and x is a vector of variables. A feasible solution of this problem is one that satisfies x ∈ R n ∧Ax ≥ b∧x ≥ 0 and an optimal solution is a feasible one that minimizes the objective function c T x. This can be found in polynomial time. A row A i corresponds to an individual linear constraint and a column A T j to a variable. The dual of a linear program P in the above form is another linear program D:

max{b T y | y ∈ R m ∧ A T y ≤ c ∧ y ≥ 0}
where A, b, c are as before and y is the vector of dual variables. Rows of the dual correspond to variables of the primal and vice versa. The objective value of any dual feasible solution is a lower bound on the optimum of P . When P is satisfiable, its dual is also satisfiable and the values of their optima meet. For a given feasible solution x of P , the slack of constraint i is slack

x(i) = A T i x-b i . Given a dual feasible solution ŷ, slack D ŷ (i)
is the reduced cost of primal variable i, rc ŷ (i). The reduced cost rc ŷ (i) is interpreted as a lower bound on the amount that the dual objective would increase over b T ŷ if x i is forced to be non-zero in the primal.

An integer linear program (ILP) is a linear program in which we replace the constraint x ∈ R n by x ∈ Z n and it is an NP-hard optimization problem.

Bayesian Networks

A Bayesian network is a directed graphical model B = G, P where G = V, E is a directed acyclic graph (DAG) called the structure of B and P are its parameters. A BN describes a normalised joint probability distribution. Each vertex of the graph corresponds to a random variable and presence of an edge between two vertices denotes direct conditional dependence. Each vertex v i is also associated with a Conditional Probability Distribution P (v i | parents(v i )). The CPDs are the parameters of B.

The approach which we use here for learning a BN from data is the score-and-search method. Given a set of multivariate discrete data I = {I 1 , . . . , I N }, a scoring function σ(G | I) measures the quality of the BN with underlying structure G. The BNSL problem asks to find a structure G that minimises σ(G | I) for some scoring function σ and it is NP-hard [START_REF] Chickering | [END_REF]. Several scoring functions have been proposed for this purpose, including BDeu [START_REF] Buntine | Theory refinement on bayesian networks[END_REF][START_REF] Heckerman | [END_REF] and BIC [START_REF] Schwarz | Estimating the dimension of a model[END_REF][START_REF] Lam | Using new data to refine a bayesian network[END_REF]. These functions are decomposable and can be expressed as the sum of local scores which only depend on the set of parents (from now on, parent set) of each vertex:

σ F (G | I) = v∈V σ v F (parents(v) | I) for F ∈ {BDeu, BIC}.
In this setting, we first compute local scores and then compute the structure of minimal score. Although there are potentially an exponential number of local scores that have to be computed, the number of parent sets actually considered is often much smaller, for example because we restrict the maximum cardinality of parent sets considered or we exploit dedicated pruning rules [START_REF] De Campos | Properties of bayesian dirichlet scores to learn bayesian network structures[END_REF][START_REF] De Campos | [END_REF]. We denote P S(v) the set of candidate parent sets of v and P S -C (v) those parent sets that do not intersect C. In the following, we assume that local scores are precomputed and given as input, as is common in similar works. We also omit explicitly mentioning I or F , as they are constant for solving any given instance.

Let C be a set of vertices of a graph G. C is a violated cluster if the parent set of each vertex v ∈ C intersects C. Then, we can prove the following property: Property 1. A directed graph G = V, E is acyclic if and only if it contains no violated clusters, i.e., for all C ⊆ V , there exists v ∈ C, such that parents(v) ∩ C = ∅.

The GOBNILP solver [START_REF] Bartlett | [END_REF] formulates the problem as the following 0/1 ILP:

min v∈V,S⊆V \{v} σ v (S)x v,S (1)
s.t. S∈P S(v) x v,S = 1 ∀v ∈ V (2) v∈C,S∈P S -C (v) x v,S ≥ 1 ∀C ⊆ V (3) x v,S ∈ {0, 1} ∀v ∈ V, S ∈ P S(v) (4)
Algorithm 1: Acyclicity Checker

acycChecker (V, D) order ← {} changes ← true while changes do changes ← f alse foreach v ∈ V \ order do if ∃S ∈ D(v) s.t. (S ∩ V ) ⊆ order then 1 order ← order + v changes ← true return order
This ILP has a 0/1 variable x v,S for each candidate parent set S of each vertex v where x v,S = 1 means that S is the parent set of v. The objective (1) directly encodes the decomposition of the scoring function. The constraint (2) asserts that exactly one parent set is selected for each random variable. Finally, the cluster inequalities (3) are violated when C is a violated cluster. We denote the cluster inequality for cluster C as cons(C) and the 0/1 variables involved as varsof (C). As there is an exponential number of these, GOBNILP generates only those that improve the current linear relaxation and they are referred to as cluster cuts. This itself is an NPhard problem [Cussens et al., 2017], which GOBNILP also encodes and solves as an ILP. Interestingly, these inequalities are facets of the BNSL polytope [Cussens et al., 2017], so stand to improve the relaxation significantly.

The CPBayes solver [ [START_REF] Van Beek | Machine learning of bayesian networks using constraint programming[END_REF] models BNSL as a constraint program. The CP model has a parent set variable for each random variable, whose domain is the set of possible parent sets, as well as order variables, which give a total order of the variables that agrees with the partial order implied by the DAG. The objective is the same as (1). It includes channelling constraints between the set of variables and various symmetry breaking and dominance constraints. It computes a lower bound using two separate mechanisms: a component caching scheme and a pattern database that is computed before search and holds the optimal graphs for all orderings of partitions of the variables. Acyclicity is enforced using a global constraint with a bespoke propagator. The main routine of the propagator is acycChecker (Algorithm 1), which returns an order of all variables if the current set of domains of the parent set variables may produce an acyclic graph, or a partially completed order if the constraint is unsatisfiable. This algorithm is based on Property 1.

Briefly, the algorithm takes the domains of the parent set variables as input and greedily constructs an ordering of the variables, such that if variable v is later in the order than v , then v / ∈ parents(v )1 . It does so by trying to pick a parent set S for an as yet unordered vertex such that S is entirely contained in the set of previously ordered vertices2 . If all assignments yield cyclic graphs, it will reach a point where all remaining vertices are in a violated cluster in all possible graphs, and it will return a partially constructed order. If there exists an assignment that gives an acyclic graph, it will be possible by property 1 to select from a variable in V \ order a parent set which does not intersect V \ order, hence is a subset of order. The value S chosen for each variable in line 1 also gives a witness of such an acyclic graph.

An immediate connection between the GOBNILP and CP-Bayes models is that the ILP variables x v,S , ∀S ∈ P S(v) are the direct encoding [Walsh, 2000] of the parent set variables of the CP model. Therefore, we use them interchangeably, i.e., we can refer to the value S in D(v) as x v,S .

Restricted Cluster Detection

One of the issues hampering the performance of CPBayes is that it computes relatively poor lower bounds at deeper levels of the search tree. Intuitively, as the parent set variable domains get reduced by removing values that are inconsistent with the current ordering, the lower bound computation discards more information about the current state of the problem. We address this by adapting the branch-and-cut approach of GOBNILP. However, instead of finding all violated cluster inequalities that may improve the LP lower bound, we only identify a subset of them.

Consider the linear relaxation of the ILP ( 1)-( 4), restricted to a subset C of all valid cluster inequalities, i.e., with equation (4) replaced by 0 ≤ x v,S ≤ 1∀v ∈ V, S ∈ P S(v) and with equation (3) restricted only to clusters in C. We denote this LP C . We exploit the following property of this LP. Theorem 1 gives a class of cluster cuts, which we call RCclusters, for reduced-cost clusters, guaranteed to improve the lower bound. Importantly, this requires only a feasible, perhaps sub-optimal, solution. Example 1 (Running example). Consider a BNSL instance with domains as shown in Table 1 and let C = ∅. Then, ŷ = 0 leaves the reduced cost of every variable to exactly its primal objective coefficient. The corresponding x assigns 1 to variables with reduced cost 0 and 0 to everything else. These are both optimal solutions, with cost 0 and x is integral, so it is also a solution of the corresponding ILP. However, it is not a solution of the BNSL, as it contains several cycles, including C = {0, 2, 3}. The cluster inequality cons(C) is violated in the primal and allows the dual bound to be increased.

Variable Domain Value Cost 0 {2} 0 1 {2, 4} 0 {} 6 2 {1, 3} 0 {} 10 3 {0} 0 {} 5 4 {2, 3} 0 {3} 1 {2} 2 {} 3
Table 1: BNSL instance used as running example.

Algorithm 2: Lower bound computation with RC-

clusters lowerBoundRC (V, D, C) ŷ ← DualSolve(LP C (D)) while True do 2 C ← V \ acycChecker(V, D rc C,ŷ ) 3 if C = ∅ then return cost(ŷ), C C ← minimise(C) C ← C ∪ {C} ŷ ← DualImprove(ŷ, LP C (D), C)
We Proof. (⇒) Let C be such a cluster. Since for all x v,S ∈ varsof (C), none of these are in D rc C,ŷ , so cons(C) is violated and hence there is no acyclic assignment.

(⇐) Consider once again acycChecker, in Algorithm 1. When it fails to find a witness of acyclicity, it has reached a point where order V and for the remaining variables C = V \ order, all allowed parent sets intersect C. So if acycChecker is called with D rc C,ŷ , all values in varsof (C) have reduced cost greater than 0, so C is an RC-cluster.

Theorem 2 shows that detecting unsatisfiability of D rc C,ŷ is enough to find an RC-cluster. Its proof also gives a way to extract such a cluster from acycChecker.

Algorithm 2 shows how theorems 1 and 2 can be used to compute a lower bound. It is given the current set of domains and a set of clusters as input. It first solves the dual of LP C (D), potentially suboptimally. Then, it uses acycChecker iteratively to determine whether there exists an RC-cluster C under the current dual solution ŷ. If that cluster is empty, there are no more RC-clusters, and it terminates and returns a lower bound equal to the cost of ŷ under LP C (D) and an updated pool of clusters. Otherwise, it minimises C (see section 3.1), adds it to the pool of clusters and solves the updated LP. It does this by calling DualImprove, which solves LP C (D) exploiting the fact that only the cluster inequality cons(C) has been added. Example 2. Continuing our example, consider the behavior of acycChecker with domains D rc ∅,ŷ after the initial dual solution ŷ = 0. Since the empty set has non-zero reduced cost for all variables, acycChecker fails with order = {}, hence C = V . We postpone discussion of minimization for now, other than to observe that C can be minimized to C 1 = {1, 2}. We add cons(C 1 ) to the primal LP and set the dual variable of C 1 to 6 in the new dual solution ŷ1 . The reduced costs of x 1,{} and x 2,{} are decreased by 6 and, importantly, rc ŷ1 (x 1,{} ) = 0. In the next iteration of lowerBoundRC, acycChecker is invoked on D rc {C1},ŷ1 and returns the cluster {0, 2, 3, 4}. This is minimized to C 2 = {0, 2, 3}. The parent sets in the domains of these variables that do not intersect C 2 are x 2,{} and x 3,{} , so minrc(C 2 ) = 4, so we add cons(C 2 ) to the primal and we set the dual variable of C 2 to 4 in ŷ2 . This brings the dual objective to 10. The reduced cost of x 2,{} is 0, so in the next iteration acycChecker runs on D rc {C1,C2},ŷ2 and succeeds with the order {2, 0, 3, 4, 1}, so the lower bound cannot be improved further. This also happens to be the cost of the optimal structure. Theorem 3. Algorithm 2 terminates but is not confluent.

Proof. It terminates because there is a finite number of cluster inequalities and each iteration generates one. In the extreme, all cluster inequalities are in C and the test at line 3 succeeds, terminating the algorithm.

To see that it is not confluent, consider an example with 3 clusters C 1 = {v 1 , v 2 }, C 2 = {v 2 , v 3 } and C 3 = {v 3 , v 4 } and assume that the minimum reduced cost for each cluster is unit and comes from x 2,{4} and x 3,{1} , i.e., the former value has minimum reduced cost for C 1 and C 2 and the latter for C 2 and C 3 . Then, if minimisation generates first C 1 , the reduced cost of x 3,{1} is unaffected by DualImprove, so it can then discover C 3 , to get a lower bound of 2. On the other hand, if minimisation generates first C 2 , the reduced costs of both x 2,{4} and x 3,{1} are decreased to 0 by DualImprove, so neither C 1 nor C 3 are RC-clusters under the new dual solution and the algorithm terminates with a lower bound of 1.

Related Work. The idea of performing propagation on the subset of domains that have reduced cost 0 has been used in the VAC algorithm for WCSPs [START_REF] Cooper | [END_REF]. Our method is more light weight, as it only performs propagation on the acyclicity constraint, but may give worse bounds. The bound update mechanism in the proof of theorem 1 is also simpler than VAC and more akin to the "disjoint core phase" in core-guided MaxSAT solvers [START_REF] Morgado | [END_REF].

Cluster Minimisation

It is crucial for the quality of the lower bound produced by Algorithm 2 that the RC-clusters discovered by acycChecker are minimised, as the following example shows. Empirically, omitting minimisation rendered the lower bound ineffective. Example 3. Suppose that we attempt to use lowerBoundRC without cluster minimization. Then, we use the cluster given by acycChecker, C 1 = {0, 1, 2, 3, 4}. We have minrc(C 1 ) = 3, given from the empty parent set value of all variables. This brings the reduced cost of x 4,{} to 0. It then proceeds to find the cluster C 2 = {0, 1, 2, 3} with minrc(C 2 ) = 2 and decrease the reduced cost of x 3,{} to 0, then C 3 = {0, 1, 2} with minrc(C 3 ) = 1, which brings the reduced cost of x 1,{} to 0. At this point, acycChecker succeeds with the order {4, 3, 1, 2, 0} and lowerBoundRC returns a lower bound of 6, compared to 10 with minimization. The order produced by acycChecker also disagrees with the optimum structure.

Therefore, when we get an RC-cluster C at line 2 of algorithm 2, we want to extract a minimal RC-cluster (with respect to set inclusion) from C, i.e., a cluster C ⊆ C, such that for all ∅ ⊂ C ⊂ C , C is not a cluster.

Minimisation problems like this are handled with an appropriate instantiation of QuickXPlain [Junker, 2004]. These algorithms find a minimal subset of constraints, not variables. We can pose this as a constraint set minimisation problem by implicitly treating a variable as the constraint "this variable is assigned a value" and treating acyclicity as a hard constraint.

However, the property of being an RC-cluster is not monotone. For example, consider the variables {v 1 , v 2 , v 3 , v 4 } and ŷ such that the domains restricted to values with 0 reduced cost are {{v 2 }}, {{v 1 }}, {{v 4 }}, {{v 3 }}, respectively. Then {v 1 , v 2 , v 3 , v 4 }, {v 1 , v 2 } and {v 3 , v 4 } are RC-clusters. but {v 1 , v 2 , v 3 } is not because the sole value in the domain of v 3 does not intersect {v 1 , v 2 , v 3 }. We instead minimise the set of variables that does not admit an acyclic solution and hence contains an RC-cluster. A minimal unsatisfiable set that contains a cluster is an RC-cluster, so this allows us to use the variants of QuickXPlain. We focus on RobustXPlain, which is called the deletion-based algorithm in SAT literature for minimising unsatisfiable subsets [START_REF][END_REF]. The main idea of the algorithm is to iteratively pick a variable and categorise it as either appearing in all minimal subsets of C, in which case we mark it as necessary, or not, in which case we discard it. To detect if a variable appears in all minimal unsatisfiable subsets, we only have to test if omitting this variable yields a set with no unsatisfiable subsets, i.e., with no violated clusters. This is given in pseudocode in Algorithm 3. This exploits a subtle feature of acycChecker as described in Algorithm 1: if it is called with a subset of V, it does not try to place the missing variables in the order and allows parent sets to use these missing variables. Omitting variables from the set given to acycChecker acts as omitting the constraint that these variables be assigned a value

Solving the Cluster LP

Solving a linear program is in polynomial time, so in principle DualSolve can be implemented using any of the commercial or free software libraries available for this. However, solving Algorithm 3: Find a minimal RC-cluster subset of C

MinimiseCluster (V, D, C) N = ∅ while C = ∅ do Pick c ∈ C C ← C \ {c} C ← V \ acycChecker(N ∪ C, D) if C = ∅ then N ← N ∪ {c} else C ← C \ N return N
this LP using a general LP solver is too expensive in this setting. As a data point, solving the instance steel BIC with our modified solver took 25,016 search nodes and 45 seconds of search, and generated 5, 869 RC-clusters. Approximately 20% of search time was spent solving the LP using the greedy algorithm that we describe in this section. CPLEX took around 70 seconds to solve LP C with these cluster inequalities once. While this data point is not proof that solving the LP exactly is too expensive, it is a pretty strong indicator. We have also not explored nearly linear time algorithms for solving positive LPs [START_REF] Orecchia ; Zeyuan Allen-Zhu | Nearly-linear time positive LP solver with faster convergence rate[END_REF].

Our greedy algorithm is derived from theorem 1. Observe first that LP C with C = ∅, i.e., only with constraints (2) has optimal dual solution ŷ0 that assigns the dual variable y v of S∈P S(v) x v,S = 1 to min S∈P S(v) σ v (S). That leaves at least one of x v,S , S ∈ P S(v) with reduced cost 0 for each v ∈ V . DualSolve starts with ŷ0 and then iterates over C. Given ŷi-1 and a cluster C, it sets ŷi = ŷi-1 if C is not an RC-cluster. Otherwise, it increases the lower bound by c = minrc(C) and sets ŷi = ŷi-1 ∪ {y C = c}. It remains to specify the order in which we traverse C.

We sort clusters by increasing size |C|, breaking ties by decreasing minimum cost of all original parent set values in varsof (C). This favours finding non-overlapping cluster cuts with high minimum cost. In section 6, we give experimental evidence that this computes better lower bounds.

DualImprove can be implemented by discarding previous information and calling DualSolve(LP C (D)). Instead, it uses the RC-cluster C to update the solution without revisiting previous clusters.

In terms of implementation, we store varsof (C) for each cluster, not cons(C). During DualSolve, we maintain the reduced costs of variables rather than the dual solution, otherwise computing each reduced cost would require iterating over all cluster inequalities that contain a variable. Specifically, we maintain ∆ v,S = σ v (S) -rc ŷ (x v,S ). In order to test whether a cluster C is an RC-cluster, we need to compute minrc(C). To speed this up, we associate with each stored cluster a support pair (v, S) corresponding to the last minimum cost found. If rc ŷ (v, S) = 0, the cluster is not an RC-cluster and is skipped. Moreover, parent set domains are sorted by increasing score σ v (S), so S S ⇐⇒ σ v (S) > σ v (S ). We also maintain the maximum amount of cost transferred to the lower bound,

∆ max v = max S∈D(v) ∆ v,S
Algorithm 4: GAC propagator for acyclicity

Acyclicity-GAC (V, D) O ← acycChecker(V, D) if O V then return Failure foreach v ∈ V do changes ← true i ← O -1 (v) pref ix ← {O 1 , . . . , O i-1 } 4 while changes do changes ← f alse foreach w ∈ O \ (pref ix ∪ {v}) do if ∃S ∈ D(w) s.t. S ⊆ pref ix then pref ix ← pref ix ∪ {w} changes ← true Prune {S | S ∈ D(v) ∧ S pref ix} return Success for every v ∈ V . We stop iterating over D(v) as soon as σ v (S) -∆ max v is greater than or equal to the current mini- mum because ∀S S, σ v (S ) -∆ v,b ≥ σ v (S) -∆ max v .
In practice, on very large instances 97.6% of unproductive clusters are detected by support pairs and 8.6% of the current domains are visited for the rest 3 .

To keep a bounded-memory cluster pool, we discard frequently unproductive clusters. We throw away large clusters with a productive ratio #productive #productive+#unproductive smaller than 1 1,000 . Clusters of size 10 or less are always kept because they are often more productive and their number is bounded.

GAC for the Acyclicity Constraint

Previously, [START_REF] Van Beek | Machine learning of bayesian networks using constraint programming[END_REF] showed that using acycChecker as a subroutine, one can construct a GAC propagator for the acyclicity constraint by probing, i.e., detecting unsatisfiability after assigning each individual value and pruning those values that lead to unsatisfiability. acycChecker is in O(n 2 d), so this gives a GAC propagator in O(n 3 d 2 ). We show here that we can enforce GAC in time O(n 3 d), a significant improvement given that d is usually much larger than n.

Suppose acycChecker finds a witness of acyclicity and returns the order O = {v 1 , . . . , v n }. Every parent set S of a variable v that is a subset of {v | v ≺ O v} is supported by O. We call such values consistent with O. Consider now S ∈ D(v i ) which is inconsistent with O, therefore we have to probe to see if it is supported. We know that during the probe, nothing forces acycChecker to deviate from {v 1 , . . . , v i-1 }. So in a successful probe, acycChecker constructs a new order O which is identical to O in the first i -1 positions and in which it moves v i further down. Then all values consistent with O are supported. This suggests that instead of probing each value, we can probe different orders.

Acyclicity-GAC, shown in Algorithm 4, exploits this insight. It ensures first that acycChecker can produce a valid 3 See the supplementary material for more.

order O. For each variable v, it constructs a new order O from O so that v is as late as possible. It then prunes all parent set values of v that are inconsistent with O . Let O be the order O p followed by Q p , followed by v, followed by O s , keeping only the first occurrence of each variable when there are duplicates. O is a valid order: O p is witnessed by the assignment that witnesses O, Q p by the assignment that witnesses Q, v by S (as in Q) and O s by the assignment that witnesses O. It also supports S, as required.

Complexity is dominated by repeating O(n) times the loop at line 4, which is a version of acycChecker so has complexity O(n 2 d) for a total O(n 3 d).

6 Experimental Results

Benchmark Description and Settings

The datasets come from the UCI Machine Learning Repository4 , the Bayesian Network Repository5 , and the Bayesian Network Learning and Inference Package6 . Local scores were computed from the datasets using B. Malone's code7 . BDeu and BIC scores were used for medium size instances (less than 64 variables) and only BIC score for large instances (above 64 variables). The maximum number of parents was limited to 5 for large instances (except for accidents.test with maximum of 8), a high value that allows even learning complex structures [START_REF] Scanagatta | [END_REF]. For example, jester.test has 100 random variables, a sample size of 4, 116 and 770, 950 parent set values. For medium instances, no restriction was applied except for some BDeu scores (limit sets to 6 or 8 to complete the computation of the local scores within 24 hours of CPU-time [Lee and [START_REF] Beek | An experimental analysis of anytime algorithms for bayesian network structure learning[END_REF]).

We have modified the C++ source of CPBayes v1.1 by adding our lower bound mechanism and GAC propagator. We call the resulting solver ELSA and have made it publicly available. For the evaluation, we compare with GOBNILP v1.6.3 using SCIP v3.2.1 with cplex v12.7.0. All computations were performed on a single core of Intel Xeon E5-2680 v3 at 2.50 GHz and 256 GB of RAM with a 1-hour (resp. 10-hour) CPU time limit for medium (resp. large) size instances. We used default settings for GOBNILP with no approximation in branch-and-cut (limits/gap = 0). We used the same settings in CPBayes and ELSA for their preprocessing phase (partition lower bound sizes l min , l max and local search number of restarts r min , r max ). We used two

Evaluation

In Table 2 we present the runtime to solve each instance to optimality with GOBNILP, CPBayes, and ELSA with default settings, without the GAC algorithm and without sorting the cluster pool (leaving clusters in chronological order, rather than the heuristic ordering presented in Section 4). For the instances with V ≤ 64 (resp. > 64), we had a time limit of 1 hour (resp. 10 hours). We exclude instances that were solved within the time limit by GOBNILP and have a search time of less than 10 seconds for CPBayes and all variants of ELSA. We also exclude 8 instances that were not solved to optimality by any method. This leaves us 17 instances to analyse here out of 69 total. More details are given in the supplemental material, available from the authors' web pages.

Comparison to GOBNILP. CPBayes was already proven to be competitive to GOBNILP [ [START_REF] Van Beek | Machine learning of bayesian networks using constraint programming[END_REF]. Our results in Table 2 confirm this while showing that neither is clearly better. When it comes to our solver ELSA, for all the variants, all instances solved within the time limit by GOBNILP are solved, unlike CPBayes. On top of that, ELSA solves 9 more instances optimally.

Comparison to CPBayes. We have made some low-level performance improvements in preprocessing of CPBayes, so for a more fair comparison, we should compare only the search time, shown in parentheses. ELSA takes several orders of magnitude less search time to optimally solve most instances, the only exception being the bnetflix instances. ELSA also proved optimality for 8 more instances within the time limit. Gain from GAC. The overhead of GAC pays off as the instances get larger. While we do not see either a clear improvement nor a downgrade for the smaller instances, search time for ELSA improves by up to 47% for larger instances compared to ELSA \ GAC.

Gain from Cluster Ordering. We see that the ordering heuristic improves the bounds computed by our greedy dual LP algorithm significantly. Compared to not ordering the clusters, we see improved runtime throughout and 3 more instances solved to optimality.

Conclusion

We have presented a new set of inference techniques for BNSL using constraint programming, centered around the expression of the acyclicity constraint. These new techniques exploit and improve on previous work on linear relaxations of the acyclicity constraint and the associated propagator. The resulting solver explores a different trade-off on the axis of strength of inference versus speed, with GOBNILP on one extreme and CPBayes on the other. We showed experimentally that the trade-off we achieve is a better fit than either extreme, as our solver ELSA outperforms both GOBNILP and CPBayes. The major obstacle towards better scalability to larger instances is the fact that domain sizes grow exponentially with the number of variables. This is to some degree unavoidable, so our future work will focus on exploiting the structure of these domains to improve performance.

  Theorem 1. Let ŷ be a dual feasible solution of LP C with dual objective o. Then, if C is a cluster such that C / ∈ C and the reduced cost rc of all variables varsof (C) is greater than 0, there exists a dual feasible solution ŷ of LP C∪C with dual objective o ≥ o + minrc(C) where minrc(C) = min x∈varsof (C) rc ŷ (x). Proof. The only difference from LP C to LP C∪C is the extra constraint cons(C) in the primal and corresponding dual variable y C . In the dual, y C only appears in the dual constraints of the variables varsof (C) and in the objective, always with coefficient 1. Under the feasible dual solution ŷ∪{y C = 0}, these constraints have slack at least minrc(C), by the definition of reduced cost. Therefore, we can set ŷ = ŷ ∪ {y C = minrc(C)}, which remains feasible and has objective o = o + minrc(C), as required.

  consider the problem of discovering RC-clusters within the CP model of CPBayes. First, we introduce the notation LP C (D) which is LP C with the additional constraint x v,S = 0 for each S / ∈ D(v). Conversely, D rc C,ŷ is the set of domains minus values whose corresponding variable in LP C (D) has non-zero reduced cost under ŷ, i.e., D rc C,ŷ = D where D (v) = {S | S ∈ D(v) ∧ rc ŷ (x v,S ) = 0}. With this notation, for values S / ∈ D(v), x v,S = 1 is infeasible in LP C (D), hence effectively rc ŷ (x v,S ) = ∞. Theorem 2. Given a collection of clusters C, a set of domains D and ŷ, a feasible dual solution of LP C (D), there exists an RC-cluster C / ∈ C if and only if D rc C,ŷ does not admit an acyclic assignment.

  . The complexity of MinimiseCluster is O(n 3 d), where n = |V | and d = max v∈V |D(v)|, a convention we adopt throughout.

Theorem 4 .

 4 Algorithm 4 enforces GAC on the Acyclicity constraint in O(n 3 d). Proof. Let v ∈ V and S ∈ D(v). Let O = {O 1 , . . . , O n } and Q = {Q 1 , . . . , Q n } be two valid orders such that O does not support S whereas Q does. It is enough to show that we can compute from O a new order O that supports S by pushing v towards the end. Let O i = Q j = v and let O p = {O 1 , . . . , O (i-1) }, Q p = {Q 1 , . . . , Q (j-1) } and O s = {O i+1 , . . . , O n }.

  different settings depending on problem size |V |: l min = 20, l max = 26, r min = 50, r max = 500 if |V | ≤ 64, else l min = 20, l max = 20, r min = 15, r max = 30.

Table 2 :

 2 Comparison of ELSA against GOBNILP and CPBayes. Time limit for instances above the line is 1h, for the rest 10h. Instances are sorted by increasing total domain size. For variants of CPBayes we report in parentheses time spent in search, after preprocessing finishes. † indicates a timeout.

	Instance	|V |	|ps(v)| GOBNILP	CPBayes	ELSA	ELSA \ GAC	ELSA chrono
	carpo100 BIC	60	424	0.6	78.5 (29.7)	40.6 (0.0)	40.7 (0.0)	40.6 (0.0)
	alarm1000 BIC 37	1003	1.2	204.2 (172.9)	27.8 (0.7)	28.8 (1.5)	29.9 (2.7)
	flag BDe	29	1325	4.4	19.0 (18.1)	0.9 (0.1)	0.9 (0.1)	1.3 (0.5)
	wdbc BIC	31	14614	99.8	629.8 (576.6)	48.9 (1.6)	49.1 (1.7)	50.3 (3.1)
	kdd.ts	64	43584	327.6	†	1314.5 (158.2)	1405.4 (239.5)	1663.2 (512.4)
	steel BIC	28	93027	† 1270.9 (1218.9)	98.0 (49.2)	99.2 (50.1)	130.0 (81.2)
	kdd.test	64	152873	1521.7	†	1475.3 (120.6)	1515.9 (128.5)	1492.4 (109.5)
	mushroom BDe 23	438186	†	176.4 (56.0)	135.4 (33.7)	137.0 (35.0)	133.7 (31.9)
	bnetflix.ts	100 446406	†	629.0 (431.4)	1065.1 (878.4)	1111.4 (931.0)	1132.4 (936.3)
	plants.test	69	520148	†	† 18981.9 (17224.0) 30791.2 (29073.0)	†
	jester.ts	100 531961	†	†	10166.0 (9697.9) 14915.9 (14470.1) 23877.6 (23325.7)
	accidents.ts	111 568160	1274.0	†	2238.7 (904.5)	2260.3 (986.1)	2221.1 (904.8)
	plants.valid	69	684141	†	†	12347.6 (8509.7) 19853.1 (15963.1)	†
	jester.test	100 770950	†	† 17637.8 (16979.2) 21284.0 (20661.9)	†
	bnetflix.test	100 1103968	† 3525.2 (3283.8)	8197.7 (7975.6)	8057.3 (7841.4)	7915.0 (7686.3)
	bnetflix.valid	100 1325818	† 1456.6 (1097.0)	9282.0 (8950.3)	10220.5 (9898.4)	9619.7 (9257.4)
	accidents.test	111 1425966	4975.6	†	3661.7 (641.5)	4170.1 (1213.6)	3805.2 (687.6)

We treat order as both a sequence and a set, as appropriate.

When propagating the acyclicity constraint it always holds that a ∩ V = a, so this statement is true. In section

3.1, we use the algorithm in a setting where this is not always the case.

http://archive.ics.uci.edu/ml

http://www.bnlearn.com/bnrepository

https://ipg.idsia.ch/software.php?id=132

http://urlearning.org

Acknowledgements

We thank the GenoToul (Toulouse, France) Bioinformatics platform for its support. This work has been partly funded by the "Agence nationale de la Recherche" (ANR-16-CE40-0028 Demograph project and ANR-19-PIA3-0004 ANTI-DIL chair of Thomas Schiex).