Majdi Elhaji

Abdelkrim Zitouni

Rached Tourki

Pierre Boulet

Jean-Luc Dekeyser

Samy Meftali

Modeling Networks-on-Chip at System Level with the MARTE UML profile

The study of Networks on Chips (NoCs) is a research field that primarily addresses the global communication in Systems-on-Chip (SoCs). The selected topology and the routing algorithm play a prime role in the performance of NoC architectures. In order to handle the design complexity and meet the tight time-to-market constraints, it is important to automate most of these NoC design phases. The extension of the UML language called UML profile for MARTE (Modeling and Analysis of Real-Time and Embedded systems) specifies some concepts for model-based design and analysis of real time and embedded systems. This paper presents a MARTE based methodology for modeling concepts of NoC based architectures. It aims at improving the effectiveness of the MARTE standard by clarifying some notations and extending some definitions in the standard, in order to be able to model complex architectures like NoCs.

I. INTRODUCTION

System-on-chip (SoC) designs supply integrated solutions to challenging design problems of embedded systems and consumer electronic domains. Much of the progress in these fields allows the designer to conceive more and more complex electronic systems and reduce time to market. Future systemson-chips will contain tens to hundreds of IP cores according to the International Technology Roadmap for Semiconductors [START_REF]International Technology Roadmap for Semiconductors[END_REF]. The on-chip communication paradigm is introduced as Network-on-Chip (NoC) by Benini and De Micheli [START_REF] Benini | Networks on chips: a new SoC paradigm[END_REF], Dally and Towles [START_REF] Dally | Route packets, not wires: on-chip interconnection networks[END_REF] and Sgroi et al [START_REF] Sgroi | Addressing the system-on-a-chip interconnect woes through communication-based design[END_REF]. Many propositions of NoC architectures for SoCs design have been presented in literature, such as, SPIN [START_REF] Ogras | Application-specific network-on-chip architecture customization via long-range link insertion[END_REF], HERMES [START_REF] Moraes | HERMES: an infrastructure for low area overhead packet-switching networks on chip[END_REF] XPIPES [START_REF] Dall'osso | Xpipes: a latency insensitive parameterized network-on-chip architecture for multiprocessor SoCs[END_REF], OCTAGON [START_REF] Karim | On-chip communication architecture for OC-768 network processors[END_REF], Mesh [START_REF] Holsmark | Design issues and performance evaluation of mesh NoC with regions[END_REF], honeycomb [START_REF] Hemani | Network on chip: An architecture for billion transistor era[END_REF].

Topology and routing algorithm are the two most important aspects that distinguish the various NoC architectures. The choice of topology is very important to provide a high level of performance. On the other hand routing on NoCs has the same principle as routing on any network. Consequently, the system becomes more and more complex, from transistors to layout. The system designer needs with a level of abstraction that focuses more on system functionality rather than low level This work has been supported by the French state/region project contract Campus Intelligence Ambiante design details and there is a demand for the development of high level modeling environments for computer-aided design.

The modeling of highly repetitive structures such as network on chip topologies in graphical form poses a particular challenge. Model driven engineering is a software development methodology where the complete system is modeled at a high abstraction level using a modeling language as UML. The field of real-time embedded software systems is a domain for which extensions to UML are required to provide more precise notations of domain specific phenomena. The UML profile for Modeling and Analysis of Real-Time Embedded systems, MARTE [START_REF]UML profile for MARTE: modeling and analysis of Real-Time embedded systems, version 1.0[END_REF] is the current standard for the SoC domain. We study in this paper how well it is suited to model complex NoC architectures.

We note that there is a few works interested in modeling concepts of NoC. Only the approaches defined in [START_REF] Boulet | Repetitive allocation modeling with MARTE[END_REF] [13] come close to our contribution by using the MDE concepts for modeling repetitive structure such as multi-processor system on chip (MPSoC) but not NoC concepts. In [START_REF] Boulet | Repetitive allocation modeling with MARTE[END_REF] the authors show the need for modeling the distribution of a parallel application onto parallel hardware architecture and describe the notation (MARTE) for modeling only regular architecture. This notation allows distributing computations to processing elements, data to shared or distributed memories, with the aim of clarifying its usage through MPSoC examples and the comparisons to other distribution notations such as in High Performance Fortran. We note that in literature there is a little work interested in modeling NoC with MARTE, defining methodology or providing the effectiveness for notation to model all concepts of NoC such as topology, routing algorithm, switching mode and communication protocol. The aim of this paper consists in proposing a novel methodology for modeling all family of NoC topologies and show that the advanced MARTE can support the modeling of routing algorithm with the use of state machines at the end to propose an extension in the HwCommunication package.

II. NETWORK ON CHIP CHARACTERISTICS

The NoC is a replacement for global interconnect and single bus architectures. The main immediate benefit of a NoC based approach is clearly due to the possibility to reuse the communication network throughout different products. Furthermore, as the complexity of integrated systems keeps growing, a NoC provides enhanced performance such as throughput and scalability in comparison with the bus communication architectures. Thus it is useful to propose a methodology for modeling NoC topologies with the aims of being in line with the evolution of semiconductor technology and reduction of the time to market. Also, from a predictability perspective, the regularity of NoC layout provides well characterized electrical and physical properties. The topology and routing algorithm are the most important aspects which distinguish the diverse NoC architectures.

A. Topologies

The topology refers to the physical structure of the network graph, i.e., how network nodes are physically connected. It defines the connectivity or the routing possibility between nodes, thus having a fundamental impact on the network performance as well as the router structure (number of ports and port width). The trade-off between generality and customization is an important issue when determining a network topology. Each topology can be characterized by a few properties. The degree of a router is the number of links connecting that node to its neighbor vertices. A topology is considered as regular when all routers have the same degree, if not it is considered as irregular. In [START_REF] Srinivasan | Linear-programming-based techniques for synthesis of network-on-chip architectures[END_REF], the number of ports in routers can be synthesized according to the requirement of connectivity. However, the area and power consumption of an irregular network topology may not scale predictably with the topology size. Besides, many network topologies have been proposed where most of them are proposed for minimizing the number of nodes and node degrees. Also there are many research activities for designing NoC topologies. Murali in [START_REF] Murali | SUNMAP: a tool for automatic topology selection and generation for NoCs[END_REF] have developed a tool for automatically selecting an application-specific topology for minimizing average communication delay, area, and power dissipation. Other topologies between regular and irregular are also proposed for NoCs. For example, an interesting NoC topology is the Spirdegon of STMicroelectronics, GEXspedirgon NoC [START_REF] Zid | Nouvelles architectures génériques de NoC= news generics architectures for NoC[END_REF] and Honeycomb [START_REF] Hemani | Network on chip: An architecture for billion transistor era[END_REF]. Spidergon, is one of the NoCs researched at STMicroelectronics and is being proposed as an evolution of the STNoC [START_REF] Coppola | Spidergon: a novel on-chip communication network[END_REF]. It is a packet-switched NoC, inspired from the Octagon NoC [START_REF] Karim | On-chip communication architecture for OC-768 network processors[END_REF], deterministic routing, wormhole switching and TDMA quality of service QoS. In [START_REF] Hemani | Network on chip: An architecture for billion transistor era[END_REF], Hemani and all present a honeycomb NoC topology and its associated methodology as solution to the design productivity problem. Besides they propose a platform to handle the complexity of emerging design of chip architecture and still allow companies to meet the time to market and make profit, but implementations and results are not presented. Another study was made in order to define a generic NoC [START_REF] Schmaltz | A generic network on chip model[END_REF] named GeNoC that is intended to serve as a reference for the design and the validation of high level specifications of communicating virtual modules, but this study is limited to topology and routing concepts. The 2D-mesh is the most used topology due to its simplicity. It consists of horizontal and vertical lines with nodes placed at their intersections. This specific structure is often used because inter-node delays can be predicted at a high level.

Figures 1 and2 show some of topologies that have been proposed in the literature.

B. Routing algorithm

Routing algorithms define the path followed by each message or packet routing on an NoC. A routing algorithm defines how the data are transmitted from sender to receiver. The choice of a routing algorithm depends on several metrics such as minimizing power, minimizing logic and routing tables to achieve a lower area, increasing performance by reducing delay and maximizing traffic utilization of the network. A lot of works aim at improving the hardware block for the routing algorithm in the NoC.

In this paper, routing algorithms are divided into two groups, deterministic and adaptive algorithms. Static routing algorithms ignore the network path diversity and are not sensitive to the network state. They are also simple and inexpensive to implement. Besides, it is often a simple way to provide the ordering of packets. Static routing also permits packets to be split among multiple paths between a source and destination, in a predetermined way. If only a single path is used, static routing usually guarantees in-order delivery of data packets. This eliminates the need for adding bits to packets at the NI (Network Interface), in order to correctly identify and reorder them at the destination.

Adaptive routing algorithms, also named dynamic algorithms, use information about network traffic and/or channel status to avoid congested or faulty regions of the network. In dynamic routing, routing decisions are made according to the current state of the network, considering a load on links. As such it is possible that the path between the source and destination changes over time, as traffic conditions and requirements of the application change. This adaptive behavior requires additional hardware resources that control the state of the network and routing paths. Besides, dynamic routing is able to utilize alternate paths when some directions become congested.

Routing algorithms can be implemented in various ways. The most interesting ways consist of either looking at a routing table or executing a finite-state machine in software or hardware. In this paper we focus on the modeling of routing algorithms as a state machine to demonstrate the ability to model the routing on NoCs in MARTE.

III. PROPOSED METHODOLOGY

Designing an efficient NoC architecture, while satisfying the application performance constraints, is a complex process. The design issues require several abstraction levels, ranging from high-level application modeling to physical layout level implementation. Some of the most important phases in designing the NoC include: NoC topology modeling for the application, choosing the routing algorithm, communication protocol and mode switching.

This methodology can be expressed in two phases. In the first phase, we begin by detailed concepts relating to a NoC, such as mathematical studies of topology and routing algorithm. Second, specify the package we will use and make the relation between hardware architecture concepts and notation that also exist in MARTE. We illustrate this below with the two main characteristics of NoCs identified above: topology and routing algorithm.

A. Topology Modeling

Several topologies, such as hypercube, star, ring and tree, can be modeled easily with the MARTE profile using the RSM (Repetitive Structure Modeling) package (see [START_REF]UML profile for MARTE: modeling and analysis of Real-Time embedded systems, version 1.0[END_REF], annex E, pages 517-533). The more complex such as the GEXspidergon and honeycomb are not as easy but we will demonstrate below that this RSM package is powerful enough to model them. In order to model these topologies in a compact way, one should make a mathematical study of the topology to identify some information of the graph such as degree, valence, the number of links, the number of input/output ports, then classify routers by their numbers of port and factorize routers if possible.

The modeling consists in using a class (and part) by router kind, express the number of routers using the « shaped » stereotype and connect them which links stereotyped with a link topology, either « interRepetition » for the simple cases or « reshape » for the more complex kinds.

We give below two examples of complex topology modeling.

1) Modeling of honeycomb Topology: The honeycomb mesh (see Figure 2c), based on hexagonal plan tessellation network. Hexagonal tessellations were used in literature for various applications such as cellular phone. Hexagonal and honeycomb mesh are dual graph one can obtain the honeycomb mesh by joining the center of each triangle in hexagonal mesh with center of neighboring triangle. One hexagon is a honeycomb mesh of size 1, the honeycomb mesh of size 2 is obtained by adding 6 hexagons to the boundary edge of hexagon size 1. A honeycomb mesh size m is obtained from honeycomb size m -1. the number of nodes of the honeycomb mesh of size t is 6m 2 , and the number of links is 9m 2 -3m, the diameter of honeycomb mesh of size t is 4m -1. The Honeycomb Architecture is constructed based on an elementary hexagonal network that can be described, in parallaxis syntax as following:

Configuration Hex [1....2], [1...3]; Connection:

East: Hex[i, j] → Hex[i, j + 1] West: Hex[i, j] → Hex[i, j -1] North: Hex[i, j] → Hex[i -1, j] South: Hex[i, j] → Hex[i, j -1]
The Honeycomb graph is constructed by multiplying this algorithm. This architecture is modeled on the same following algorithm to model GEXspidergon; in a customizable way on the environment GASPARD by using the RSM package, in order to take into account the modeling of link topologies. The concerned topology is composed of the repetition of a single element (Router or port). Each potential instance of this element is connected to other potential instances of the same element. In our case each instance is connected to neighbors located at north, south, east and west. The inter Repetition topology enables to specify the position of every neighbor of every potential instance of a model element with a multidimensional shape. With this architecture it is to the best of our knowledge the first time that a reshape connector is used with the same port as source and destination, allowing to represent regular but not uniform dependencies between repetitions of a part.

Figure 3 shows the modeling of honeycomb topology. 2) Modeling of GEXspidergon Topology: The GEXspidergon (see Figure 2a) can be seen as a hierarchical and globally regular, locally regular topology. It is an academic topology for NoC. This study presents a generic NoC architecture based on a configurable router. This router integrates a sophisticated dynamic arbiter, the wormhole routing technique and can be configured in a manner that allows it to be used in many possible NoC topologies such as Mesh 2-D, Tree and Polygon architectures. This makes it possible to improve the quality of service (QoS) required by the proposed NoC. This study [START_REF] Coppola | Spidergon: a novel on-chip communication network[END_REF] shows that the Spidergon architecture is characterized by the lower latency and saturation. This architecture is constructed based on an elementary polygon network which is a combination of the star and the ring architectures. This elementary network is formed by 4R+1 (R = 1, 2, etc.) routers including a central router that is connected with the 4R peripheral routers via point to point links. The peripheral routers are connected to each other in the form of a ring. The elementary network is thus characterized by its valence (m = 4R) that represents the number of the peripheral routers. These routers necessitate 2m links to be connected to the central router. Each peripheral router is connected to 4 input/output ports and the central router is connected to m+1 input/output ports. This elementary network can be described in Parallaxis [START_REF] Braunl | Parallaxis-III: a structured data-parallel programming language[END_REF] The GEXspidergon graph is constructed by iterating this algorithm in two dimensions as illustrated in Figure 2. Routers are categorized into four groups according to their degree:

• Corner Routers.

• Middle Routers.

• Horizontal Routers, package of two.

• Vertical Routers, package of two. Figure 4 shows a possible model of this topology. To express the link topology between vertical/horizontal routers and middle routers there was some problem since the considered tiles are not sets of regularly spaced point and in MARTE/RSM the points of the tile must be regularly spaced (as they are built from the reference point of the tile by the linear combination of the column vectors of the fitting matrix). We had to use two reshape links to express such a topology. These two reshapes are detailed below.

Reshape1: VR to MR

    0 0 0 2     paving=     1 0 0 1 0 0 0 0     fitting=     0 0 1 0     Target Tiler: origin=   0 0 0   paving=   1 0 0 1 0 0   fitting=   0 0 1   Reshape2: VR to MR Repetition Space={n, m} Pattern Shape={2} Source Tiler: origin=     0 1 0 0     paving=     1 0 0 1 0 0 0 0     fitting=     0 0 1 0     Target Tiler: origin=   0 0 5   paving=   1 0 0 1 0 0   fitting=   0 0 1  
With this case study we have shown that we can model the complex regular topologies with MARTE/RSM in a compact way, exploiting fully the regularity of the architecture of the NoC.

B. Modeling Routing Algorithms

In this paper routing algorithms, are divided into two groups, deterministic and adaptive algorithms. We are interested in modeling an example for each family.

1) Deterministic XY algorithm: The XY routing is considered as a kind of deterministic routing algorithms. It can be used in a 3 × 3 Mesh topology NoC. In this architecture each router is identified by its coordinates (x, y).

This algorithm begins by comparing the current router address (C x , C y ) to the destination router address (D x , D y ) of the data, stored in the header. Packet must attain the port of the router when (C x , C y ) is equal to (D x , D y ) address. If not, D x is compared to Cx address. Packet will be routed to the west port when C x > D x , to east when C x < D x and if C x ∼ = D x the packet is horizontally routed. When the last condition is true, the D y address is compared to the C y . When C y > D y the packet is routed to the North but if C y < D y it will be routed to the South. Packet will be blocked in the Fifo when the chosen port is busy. This explication can be described with the following algorithm. 2) Adaptive routing: In adaptive routing, paths are dynamically changed during the routing process according to network conditions, such as the presence of congestion. If a channel is busy, the other channel has priority over the busy channel. However, if both output channels are not used, a cost function (CF) decides which output channel will be used. This cost function is required when an adaptive routing is implemented. The adaptive routing algorithm increases the complexity of the hardware architecture of the router. That is why it will be significant to model this algorithm at a high level of abstraction. In this case a typical example of adaptive routing algorithm is modeled, the double Y-channel routing algorithm. This algorithm applies to a 2D mesh where router are connected by one bidirectional links in the X dimension and by two bidirectional links in the Y dimension. The network is divided into two sub-networks. The X1 sub-network uses paths in the ascending X dimension. The X2 sub-network uses paths in the descending X dimension.

We consider S and D the source and destination router of a packet. Their coordinates along the X axis are noted S x and D x . If, at some router the destination of a packet is on the right the packet uses the X1 sub-network; if the destination is on the left the packet used the X2 sub-network. When D x = S x a cost function decides the channel to use. In each subnetwork, several minimal paths are possible. A packet crosses the network through a single sub-network. This explication can be described with the following algorithm. For the QoS there are a lot of work that integrate several services in the routing algorithm such as a dynamic arbiter, guaranteed throughput and TDMA. The designer can consider all these concepts in the modeling. On the other side, the switching technique can be an enumerated type in the MARTE profile.

The last thing that we need to model NoCs at a high level of abstraction is to decide which modeling artifact we will use to model a router. Indeed, this element is the basic building block of NoCs and with which the topology can be described and that should be characterized by a behavioral state machine describing the routing algorithm. It should thus be a structured class with an adequate stereotype.

The HwCommunication package in the MARTE profile defines some concepts for modeling the communication infrastructure of embedded systems as « HwEndPoint » which may present the network interface of the NoC, « HwBridge » to make connection between resources and the « HwArbiter » to control the communication. But in this package the router (the basic building of NoCs as identified above) is not found. That is why we propose to add the « HwRouter » stereotype as a specialization of the « HwMedia » stereotype in the HwCommunication package. This addition would enable the modeling of NoCs at an abstraction level useful for specifying and analyzing the communication possibilities of the current (and probably future) networks used in recent systems-on-chip.

IV. CONCLUSION

We have proposed in this paper a methodology for modeling NoCs with the MARTE standard profile. We have identified two main characteristics of NoCs that are useful to design, analyse and understand the structure and behavior of NoCs: the topology of the network and the routing algorithm.

Several topologies have been proposed in the literature and we have shown that the Repetitive Structure Modeling package of MARTE is powerful enough to model as well the simple regular ones as the more complex ones. To achieve this, we have proposed the first use of the « Reshape » connector to connect some ports of the same part, enabling the description of regular but non uniform connections between repeated parts. The proposed modeling methodology enables to take full advantage of the regularity in the topology to factorize the model. Such factorization could then be used to generate factorized analysis or synthesis code and thus tackle large networks (think many-cores with several hundreds or thousands of routers) as easily are small ones.

The second characteristics that we have identified is the routing algorithm. We have proposed to model it with a standard UML behavioral state machine. This state machine should be attached to a component stereotyped « HwRouter », a new stereotype that we propose to identify the basic building block of NoCs. With this minor addition, the MARTE profile is complete enough to model a very large number of NoC proposals (actually, all those that we have found), especially those that exhibit some regularity enabling a scaling of the network. Some work is on-going to synthesize such networks in VHDL from such models.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Examples of standard topologies

Fig. 3 .

 3 Fig. 3. Modeling of a honeycomb topology Mathematical models for Reshape and InterRepetition are described below: Inter Repetition={1,0} Reshape include Tiler source and Tiler target: Pattern shape={0}, Repetition space={0} Source Tiler origin= 0 0 paving= 2 1 0 1 fitting= 0

Fig. 4 .

 4 Fig. 4. Modeling of a Gexspidergon topology

  Destination router: (D x , D y ) Current router: (C x , C y ) If (D x > C x ) // Move East Return E Else if (D x < C x ) // Move West Return w Else if (D x = C x ) // same direction If (D y < C y ) // Move South Return S Else if (D y > C y ) // Move North Return N Else if (D y = C y ) //same router Return C Thus the routing algorithm can be modeled via MARTE profile as a hardware component where the behavior is described with a behavioral state machine MARTE as shown in Figure 5.

Fig. 5 .

 5 Fig. 5. A high level abstraction of XY routing via UML/MARTE profile

  Destination router: (D x , D y ) Source router: (S x , S y ) If (D x > S x ) // Move in X1 Return X1 Else if (D x < S x ) //Move in X2 Return X2 Else if (D x = C x ) // CF Return CF This algorithm can be simply modeled in state machine via UML/MARTE profile as shown in Figure 6.

Fig. 6 .

 6 Fig. 6. A high level abstraction of adaptive routing via UML/MARTE profile


Christian Schäfer

Parametric families on large binary spaces

In the context of adaptive Monte Carlo algorithms, we cannot directly generate independent samples from the distribution of interest but use a proxy which we need to be close to the target.

Generally, such a proxy distribution is a parametric family on the sampling spaces of the target distribution. For continuous sampling problems in high dimensions, we often use the multivariate normal distribution as a proxy for we can easily parametrise it by its moments and quickly sample from it.

Our objective is to construct similarly flexible parametric families on binary sampling spaces too large for exhaustive enumeration. The binary sampling problem seems more difficult than its continuous counterpart since the choice of a suitable proxy distribution is not obvious.

Parametric families and Monte Carlo

Adaptive Monte Carlo

A Monte Carlo algorithm is said to be adaptive it is able to adjust, sequentially and automatically, its sampling distribution to the problem at hand. Precisely, the algorithms are able to incorporate information obtained from past simulations to improve the sampling distribution q in terms of nearness to the target distribution π.

Some important classes of adaptive Monte Carlo algorithms are Adaptive Importance Sampling (e.g. [START_REF] Cappé | Adaptive importance sampling in general mixture classes[END_REF], Adaptive Markov chain Monte Carlo (e.g. [START_REF] Andrieu | A tutorial on adaptive MCMC[END_REF], Sequential Monte Carlo [START_REF] Del Moral | Sequential monte carlo samplers[END_REF] and the Cross-Entropy method [START_REF] Rubinstein | The Cross-Entropy Method: A Unified Approach to Combinatorial Optimization[END_REF].

For the sampling distribution, we usually select a suitable parametric family q = q θ and adjust its parameter θ during the course of the adaptive algorithm. In continuous sampling spaces, good results are often achieved using normal distributions N (µ, Σ), for they reproduce the marginals and covariance structure of the target. * CREST and Université Paris Dauphine • christian.schafer@ensae.fr 1.2 Data from the target distribution In the sequel, let d > 0 denote the dimension of the binary space B d = {0, 1} d . Adaptive Monte Carlo algorithms are generally able to produce a, not necessarily independent and possibly weighted, sample w = (w 1 , . . . , w n ) ∈ [0, 1] n , X = (x 1 , . . . , x n ) ⊺ ∈ B n×d from the target distribution π we want to emulate using a binary model. We define the index set D = {1, . . . , d} and denote by xi

de f = ∑ n k=1 w k x k,i , xi, j de f
= ∑ n k=1 w k x k,i x k, j , i, j ∈ D (1) the weighted first and second sample moments. We further define by

r i, j de f = xi, j -xi x j xi (1 -xi ) x j (1 -x j )
, i, j ∈ D.

(2)

the weighted sample correlation.

Suitable parametric families

We first frame some properties making a parametric family suitable as sampling distribution in adaptive Monte Carlo algorithms.

(a) For reasons of parsimony, we want to construct a family of distributions with at most dim(θ ) ≤ d(d + 1)/2 parameters.

(b) Given a sample X = (x 1 , . . . , x n ) ⊺ from the target distribution π, we need to estimate θ * such that the binary model q θ * is close to π.

(c) We need to generate samples Y = (y 1 , . . . , y m ) ⊺ from the model q θ . We need the rows of Y to be independent.

(d) For some algorithms, we need to evaluate the probability q θ (y). For instance, we need q θ (y) to compute importance weights or acceptance ratios in the context of Importance Sampling or Markov chain Monte Carlo, respectively.

(e) Analogously to the multivariate normal, we need our calibrated binary model q θ * to reproduce the marginals and covariance structure of π.

Distributions on binary spaces

Before we embark on the discussion of binary models, we make some observations which hold true for every binary distribution. The notation and results introduced in this section will be used throughout the rest of this work. Here, we denote by π some generic distribution on B d Moments We use the short notation, u I (γ)

de f = ∏ i∈I γ i , I ⊆ D,
for the product of all components index by I with ∏ i∈ / 0 = 1. Since u I (γ) = 1 iff γ i = 1 for all i ∈ I, u I is the indicator function for the unit vector 1 |I| . We can characterize every distribution on B d by 2 d -1 full probabilities

p I de f = P π γ I = 1, γ D\I = 0 , I ⊆ D
or by 2 d -1 cross-moments, that is marginal probabilities,

m I de f = E π (u I (γ)) = P π (γ I = 1) , I ⊆ D.
In the following, we assume that m i ∈ (0, 1) for all i ∈ D, since for m i ∈ {0, 1}, the component γ i = m i is constant and therefore not part of the sampling problem.

For the product of components normalized to have zero mean and unit variance, we write v I (γ)

de f = ∏ k∈I (γ i -m i )/ m i (1 -m i ), I ⊆ D.
Note that E π (v i, j ) is the correlation between γ i and γ j . Therefore, we call

c I de f = E π (v I (γ))
the correlation of order |I|.

Marginals

We use the notation

π I (γ I ) = ∑ ξ ∈B d-|I| π(γ I , ξ ), I ⊆ D.
for the marginal distributions. Note the connection to the cross-moments

π I (1 |I| ) = ∑ ξ ∈B d-|I| π(1 |I| , ξ ) = ∑ γ∈B d u I (γ) π(γ) = m I . ( 3
)
Representations For any function f : (0, 1) → R, we can write

f (π(γ)) = ∑ I⊆D ∏ i∈I γ i ∏ i∈D\I (1 -γ i ) f (p I ).
If f has an inverse f -1 , there are thus coefficients a I such that π(γ) = f -1 (∑ I⊆D a I u I (γ)).

(4)

Constraints

The general constraints on binary data are

(∑ i∈I m i -|I| + 1) ∨ 0 ≤ m I ≤ min {m K | K ⊆ I} , (5)
where the upper bound is the monotonicity of the measure, and the lower bound follows from

|I| -1 = ∑ γ∈B d (|I| -1)π(γ) ≥ ∑ γ∈B d (∑ i∈I γ i -u I (γ)) π(γ) = ∑ i∈I m i -m I .
In fact, m I is a |I|-dimensional copula with respect to the expectations m i for i ∈ I, see Nelsen (2006, p.45), and the inequalities (5) correspond to the Fréchet-Hoeffding bounds.

Sampling For sampling from a binary distribution π, we apply the chain rule factorization

π(γ) = π {1} (γ 1 ) ∏ d i=2 π {1:i} (γ i | γ 1:i-1 ) = π {1} (γ 1 ) ∏ d i=2 π {1:i-1} (γ 1:i-1 )/π {1:i} (γ 1:i ), (6)
which permits to sample a random vector component-wise, conditioning on the entries we already generated. We do not even need to compute the full decomposition ( 6), but only the conditional probabilities π {1:i}

(γ i = 1 | γ 1:i-1 ) defined by π {1:i} (γ 1:i-1 , 1) π {1:i} (γ 1:i-1 , 1) + π {1:i} (γ 1:i-1 , 0) . ( 7
)
The full probability π(γ) is then computed as a by-product of the sampling Procedure 1.

Procedure 1 Sampling via chain rule factorization

y = (0, . . . , 0), p ← 1 for i = 1 . . . , d do r ← π {1:i} (γ i = 1 | γ 1:i-1 ) sample y i ∼ B r p ← p • r if y i = 1 p • (1 -r) if y i = 0 end for return y, p

Product models

The simplest non-trivial distributions on B d are certainly those having independent components.

Definition

For a vector m ∈ (0, 1) d of marginal probabilities, we define the product model

q m (γ) de f = ∏ i∈D m γ i (1 -m i ) 1-γ = ∏ i∈D (1 -m i ) exp (∑ i∈D logit(m i )) . ( 8
)
The second representation using the logit function

logit : (0, 1) → R, logit(p) = log p -log(1 -p) (9)
is useful to identify the product model as special case of more complex models. Later, we rather write B m instead of q m for the product model, since (8) is the generalization of the Bernoulli distribution to d dimensions.

Properties

We check the requirement list from Section 1.3:

(a) The product model is parsimonious with dim(θ ) = d.

(b) The maximum likelihood estimator m * is the sample mean (1).

(c) We easily sample from q m , since (6) holds trivially.

(d) We easily evaluate the probability of a product of independent components.

(e) The model q m does not reproduce any dependencies we might observe in the data X.

The last point is a weakness which makes this simple model impractical when adaptive Monte Carlo algorithms are applied to challenging sampling problems. The product model q m is often not flexible enough to come sufficiently close to the target distribution π. Therefore, the rest of this paper deals with ideas on how to sample binary vectors with a given dependence structure.

Beyond the product model

There are, to our knowledge, two main strategies to produce binary vectors with correlated components.

(1) We can construct a generalized linear model which permits computation of its marginal distributions. We apply the chain rule factorization (6) and write q θ as

q θ (γ) = q θ (γ 1 ) ∏ d i=2 q θ (γ i | γ 1:i-1 ), (10)
which allows us to sample vectors component-wise.

(2) We sample from a multivariate auxiliary distribution h θ and map the samples into B d . We call

q θ (γ) = τ -1 (γ) h θ (v)dv (11)
a copula model, although we refrain from working with explicit uniform marginals [START_REF] Mikosch | Copulas: Tales and facts[END_REF].

In the following, we first study a few generalized linear models and then review a some copula approaches.

Linear models

Taking f the identity mapping in (4), we obtain a full linear representation π(γ) = ∑ I⊆D a I u I (γ).

However, we cannot give a useful interpretation of the coefficients a I . [START_REF] Bahadur | A representation of the joint distribution of responses to n dichotomous items[END_REF] derived the following representation:

Proposition 1. Define the index set

I de f = ∪ k∈D {I ⊆ D | |I| = k} .
Then we can write any binary distribution as

π(γ) = B m (γ) (1 + ∑ I∈I k v I (γ) c I ),
where m = (m 1 , . . . , m d ) are the marginal probabilities.

Proof. For convenience, we give the proof proposed by Bahadur in Appendix 10.1.

This decomposition, first discovered by Lazarsfeld, is a special case of a more general interaction theory [START_REF] Streitberg | Lancaster interactions revisited[END_REF] and allows for a reasonable interpretation of the parameters. Indeed, we have a product model times a correction term 1 + ∑ I∈I k v I (γ) c I where the coefficients are higher order correlations.

Definition

We can try to construct a more parsimonious model by removing higher order interaction terms. For additive approaches, however, we face the problem that a truncated representations do not necessarily define probability distributions since they might not be non-negative.

Still, for a symmetric matrix A, we define the d(d + 1)/2 parameter model

q A,a 0 (γ) = µ(a 0 + γ ⊺ Aγ), ( 12
)
where µ is a normalization constant and we set a 0 = -(min γ∈B d γ ⊺ Aγ ∧ 0). Since a 0 is the solution of an NP hard quadratic unconstrained binary optimisation problem, this definition is of little practical value.

Moments

In virtue of the linear structure, we can derive explicit expressions for the cross-moments and marginal distributions, explicit meaning that the complexity is polynomial in d. The proofs are basic but rather tedious, so we moved them to the appendix section.

Next, we give a general formula yielding all crossmoments, including the normalization constant.

Proposition 2. For a set of indices I ⊆ D, we can write the corresponding cross-moment as

m I = 1 2 |I| + ∑ i∈I 2 ∑ j∈D a i, j + ∑ d j∈I\{i} a i, j 2 |I| (4a 0 + 1 ⊺ A1 + tr [A])
.

For a proof see Appendix 10.2

Corollary 1. The normalization constant is

µ = 2 -d+2 (4a 0 + 1 ⊺ A1 + tr [A]) -1 ,
and the expected value is

E q A,a 0 (γ i ) = 1 2 + ∑ d k=1 a i,k 4a 0 + 1 ⊺ A1 + tr [A]
.

The mean m i is close to 1/2 unless the row a i dominates the matrix. Therefore, if A is non-negative definite, the marginal probabilities m i can hardly take values at the extremes of the unit interval.

Marginals

For the marginal distributions

q (1:k) A,a 0 (γ 1:k ) = ∑ ξ ∈B d-(k+1) q A,a 0 (γ 1:k , ξ )
there are explicit and recursive formulae. Hence, we can compute the chain rule decomposition (6) which in turn allows to sample from the model.

Proposition 3. For the marginal distribution holds q

(1:k)

A,a 0 (γ 1:k ) = µ2 d-k-2 s k (γ 1:k ),
where

s k (γ 1:k ) = 4a 0 + ∑ k i=1 γ i ∑ k j=1 γ j a i, j + ∑ d j=k+1 a i, j + ∑ d i=k+1 ∑ d j=k+1 a i, j + ∑ d i=k+1 a i,i .
For a proof see Appendix 10.3

Recall the connection between marginal distributions and moments we observed in (3). For γ I = 1 we obtain

s I (1 k ) = 4a 0 + 4 ∑ i∈I (∑ j∈I a i, j + ∑ j∈I c a i, j ) + ∑ i∈I c ∑ j∈I c a i, j + ∑ i∈I c a i,i = 4a 0 + ∑ i∈D ∑ j∈D a i, j + ∑ i∈D a i,i + 3 ∑ i∈I ∑ j∈I a i, j + 2 ∑ i∈I ∑ j∈I c a i, j -∑ i∈I a i,i = 4a 0 + 1 ⊺ A1 + tr [A] + ∑ i∈I 2 ∑ j∈D a i, j + ∑ j∈I\{i} a i, j ,
and π I (1 k ) = µ2 d-|I|-2 s I (1 k
) is indeed the expression for the cross-moments in Proof of Proposition 2.

Fitting the model

Given a sample X = (x 1 , . . . , x n ) ⊺ ∼ π from the target distribution, we can determine a 0 and a matrix A such that the model q A,a 0 fits the first and second sampling moments

x{i, j} = n -1 ∑ n k=1 x k,i x k, j , i, j ∈ D
by solving a linear system of dimension d(d + 1)/2 + 1. We first use the bijection

τ : D × D → {1, . . . , d(d + 1)/2}, τ(i, j) = i(i -1)/2 + j to map symmetric matrices into R (d+1) d/2
. Precisely, for the matrices A and X, we define the vectors

âτ(i, j) de f = a i, j , xτ(i, j) de f = xi, j
and the design matrix (l) .

ŝτ(i, j),τ(k,l) de f = 2 1 {i, j} (k)+1 {i, j,k}
Note that |â| = 1 ⊺ A1+tr [A].
We then equate the distribution moments to the sample moments and normalize such that

2 d-2 (I a 0 + 1 /4 Ŝâ) = x, 2 d-2 (4a 0 + |â|) = 1. ( 13
)
The solution of the linear system

â * a * 0 = 2 -d+2 1 /4 Ŝ 1 4 1 ⊺ 1 -1 x 1
is finally transformed back into a symmetric matrix A * . Since the design matrix does not depend on the data, fitting several models to different data on the same space B d is extremely fast.

Properties

We check the requirement list from Section 1.3:

(a) The linear model is sufficiently parsimonious having dimension dim(θ

) = d(d + 1)/2.
(b) We can fit the parameters A and a 0 via method of moments. However, the fitted function q A * ,a * 0 (γ) is usually not a distribution.

(c) We can sample via chain rule factorization.

(d) We can evaluate q A,a 0 (y) via chain rule factorization while sampling.

(e) The model q A,a 0 reproduces the mean and correlations of the data X.

Since in applications, the fitted matrix A * is hardly ever positive definite, we cannot use the linear model in an adaptive Monte Carlo context. As other authors [START_REF] Park | A Simple Method for Generating Correlated Binary Variates[END_REF][START_REF] Emrich | A method for generating high-dimensional multivariate binary variates[END_REF] remark, additive representations like Proposition 1 are instructive but we cannot derive practical models from them.

Log-linear models

If π(γ) > 0 for all γ ∈ B d , we can use f = log in ( 4) and obtain a full log-linear representation

π(γ) = exp ∑ I⊆D a I u I (γ) .
Note that we assume the probability mass function π is assumed to be log-linear in the parameters a I . In the context of contingency tables the term "log-linear model" refers to the assumption that the marginal probabilities m I are loglinear in the higher order marginals.

Remark Contingency table analysis is a well studied approach to modeling discrete data [START_REF] Bishop | Discrete multivariate analysis: Theory and Practice[END_REF][START_REF] Christensen | Log-linear models and logistic regression[END_REF]. For binary data, the underlying sampling distribution is assumed to be multinomial which requires an enumeration of the state space we want to avoid. [START_REF] Gange | Generating Multivariate Categorical Variates Using the Iterative Proportional Fitting Algorithm[END_REF] uses the Iterative Proportional Fitting algorithm [START_REF] Haberman | Algorithm AS 51: Log-linear fit for contingency tables[END_REF] from log-linear interaction theory to construct a binary distribution with given marginal probabilities. The fitting procedures, however, require storage of all configurations π I (γ I ) and the construction of the joint posterior from the fitted marginal probabilities. The method is powerful and exact but computationally infeasible even for moderate dimensions.

Definition

Removing higher order interaction terms, we can construct a d(d + 1)/2 parameter model q µ,A (γ)

de f = µ exp(γ ⊺ Aγ), ( 14
)
where A is a symmetric matrix. We immediately recognize the product model (8) as the special case µ [START_REF] Cox | A note on the quadratic exponential binary distribution[END_REF] refer to this version of the log-linear model as quadratic exponential model.

= ∏ i∈D (1 -m i ) d and A = diag [logit(m)].

Marginals

The moments or marginal distributions of q A are sums of exponentials which, in general, do not simplify to expressions that are polynomial in d. Therefore, we cannot perform a chain rule factorization (6) to sample from the model. [START_REF] Cox | A note on the quadratic exponential binary distribution[END_REF] proposed the following second degree Taylor approximations to the marginal distributions which are again of the form (14).

Proposition 4. We write the parameter A as

A = A ′ b ⊺ b c , ( 15
)
and define the parameters

Ãd-1 = A ′ + 1 + tanh( c 2 ) diag [b] + 1 2 sech 2 ( c 2 )bb ⊺ , μd-1 = µ(1 + exp(c))
Then q Ãd-1 (γ 1:d-1 ) is the second degree Taylor approximation to the marginal distribution q A 1:d-1 (γ 1:d-1 ). For a proof see Appendix 10.4.

If we recursively compute q Ãd-1 , . . . , q Ã1 , we can derive approximate conditional probabilities using (7). Precisely, we have

q Ãi (γ i = 1 | γ 1:i-1 ) = logit -1 ( ci + b⊺ i γ 1:i-1 ), (16)
where logit -1 (x) = (1 + exp(-x)) -1 and ci , bi are parts of the matrix Ãi according to the notation introduced in (15). In particular, ( 16) is a logistic regression. We come back to this class of models in the following Section 6. We can sample from the proxy q Ã(γ )

de f = ∏ i∈D q Ãi (γ i | γ 1:i-1 ) ≈ q A (γ),
which is close to the original log-linear model. The goodness of the approximation might be improved by judicious permutation of the components. The approximation error is hard to control, however, since we repeatedly apply the second degree approximation and propagate initial errors.

Fitting the model

As in section 4.4, we use the bijection

τ : D × D → {1, . . . , d(d + 1)/2}, τ(i, j) = i(i -1)/2 + j
to map symmetric matrices into R (d+1) d/2 . Precisely, for the matrices A and X, we define the vectors âτ(i, j)

de f = a i, j , xτ(i, j) de f = xi, j .
We let y k = log π(x k ) for k = 1, . . . , n and fit the model solving the least square problem min â∈R (d+1) d/2 Xây 2 which yields the parameters

a * i, j = [( X⊺ X) -1 X⊺ y] τ(i, j) .
Note that in most adaptive Monte Carlo algorithms that involve importance sampling or Markov transitions, the probabilities π(x k ) of the target distribution are already computed such that the fitting procedure is rather fast.

Properties

We check the requirement list from Section 1.3: (c) We can sample from an approximation q Ã(γ ) ≈ q A (γ) to the log-linear model. However, we cannot control the approximation error.

(d) We can evaluate q A,a 0 (y) up to the normalisation constant µ which suffices for most adaptive Monte Carlo methods.

(e) The model q A,a 0 reproduces the mean and correlations of the data X.

Logistic models

In the previous section we saw that even for a rather simple non-linear model we cannot derive closed-form expressions for the marginal probabilities. Therefore, instead of computing the marginals for a d-dimensional model q θ (γ), we directly fit univariate models

q b i (γ i = 1 | γ 1:i-1 ), i ∈ D
to the conditional probabilities π(γ i = 1 | γ 1:i-1 ) of the target function. Precisely, we postulate the logistic relation

logit(P π (γ i = 1)) = b i,i + ∑ i-1 j=1 b i, j γ j , i ∈ D
for the marginal probabilities of the target distribution π. We defined the logit function in (9).

Definition

For a d-dimensional lower triangular matrix B, we define the logistic model as

q B (γ) de f = ∏ i∈D B p(b i,i +b ⊺ i,1:i-1 γ 1:i-1 ) (γ i ) (17) = exp ∑ i∈D (γ i -1)(b i,i + b ⊺ i,1:i-1 γ 1:i-1 ) -log(1 + exp(b i,i + b ⊺ i,1:i-1 γ 1:i-1 ))
where B p is the Bernoulli distribution and

p(x) = logit -1 (x) = (1 + exp(-x)) -1
the inverse-logit function. We identify the product model B m as the special case

B = diag [logit(m)].
The logistic model is not a log-linear model. Note that there are d! possible logistic models and we arbitrarily pick one while there should be a permutation σ (D) of the components which is optimal in a sense of nearness to the data. In practice, however, changing the parametrisation does not seem to have a noticeably impact on the quality of the adaptive Monte Carlo algorithm.

Sparse logistic regressions

The major drawback of all multiplicative models is the fact that they do not have closed-form likelihood-maximizers such that the parameter estimation requires costly iterative fitting procedures. Therefore, we construct a sparse version of the logistic regression model which we can estimate faster than the saturated model.

Instead of fitting the saturated model q(γ i | γ 1:i-1 ), we preferably work with a more parsimonious regression model like q(γ i | γ L i ) for some index set L i ⊆ {1, . . . , i -1}, where the number of predictors #L i is typically smaller than i -1.

We solve this nested variable selection problem using some simple, fast to compute criterion. For ε about 1 /100, we define the index set

I de f = {i = 1, . . . , d | xi / ∈ (ε, 1 -ε) }.
which identifies the components which have, according to the data, a marginal probability close to either boundary of the unit interval.

We do not fit a logistic regression for the components i ∈ I. We rather set L i = / 0 and draw them independently, that is we set b i,i = logit( xi ) and b i,-i = 0 which corresponds to logistic model without predictors. The reason is twofold. Firstly, interactions do not really matter if the marginal probability is excessively small or large. Secondly, these components are prone to cause complete separation in the data or might even be constant.

For the conditional distribution of the remaining components I c = D \ I, we construct parsimonious logistic regressions. For δ about 1 /10, we define the predictor sets

L i de f = { j = 1, . . . , i -1 | δ < r i, j }, i ∈ I c ,
which identifies the components with index smaller than i and significant mutual association.

Fitting the model

Given a sample X = (x 1 , . . . , x n ) ⊺ ∼ π from the target distribution we regress y [i] = X i on the columns Z [i] = (X 1:i-1 , 1), where the column Z [i] i yields the intercept to complete the logistic model.

We maximise the log-likelihood function ℓ(b) = ℓ(b | y, Z) of a weighted logistic regression model by solving the first order condition ∂ ℓ/∂ β = 0. We find a numerical solution via Newton-Raphson iterations

- ∂ 2 ℓ(b [r] ) ∂ bb ⊺ (b [r+1] -b [r] ) = ∂ ℓ(b [r] ) ∂ b , r > 0, (18)
starting at some b [0] ; see Procedure 2 for the exact terms.

Other updating formulas like Iteratively Reweighted Least Squares or quasi-Newton iterations should work as well.

Procedure 2 Fitting the weighted logistic regressions

Input: w = (w 1 , . . . , w n ), X = (x 1 , . . . , x n ) ⊺ , B ∈ R d×d for i ∈ I c do Z ← (X L i , 1), y ← X i , b [0] ← B i,L i ∪{i} repeat p k ← logit -1 (Z k b [r-1] ) for all k = 1, . . . , n q k ← p k (1 -p k ) for all k = 1, . . . , n b [r] ← (Z ⊺ diag [w] diag [q] Z + εI n ) -1 × (Z ⊺ diag [w]) diag [q] Z b [r-1] + (y -p) until |b [r] j -b [r-1] j | < 10 -3 for all j B i,L i ∪{i} ← b end for return B
Sometimes, the Newton-Raphson iterations do not converge because the likelihood function is monotone and thus has no finite maximizer. This problem is caused by data with complete or quasi-complete separation in the sample points [START_REF] Albert | On the existence of maximum likelihood estimates in logistic regression models[END_REF]. There are several ways to handle this issue.

(a) We just halt the algorithm after a fixed number of iterations and ignore the lack of convergence. Such proceeding, however, might cause uncontrolled numerical problems.

(b) [START_REF] Firth | Bias reduction of maximum likelihood estimates[END_REF] proposes to use a Jeffrey's prior on b. The penalized log-likelihood does have a finite maximizer but requires computing the derivatives of the Fisher information matrix.

(c) We just add a simple quadratic penalty term εβ ⊺ β to the log-likelihood to ensure the target-function is convex and does not cause numerical problems.

(d) As we notice that some terms of b i are growing beyond a certain threshold, we move the component i from the set of components with associated logistic regression model I c to the set of independent components I.

In practice, we recommend to combine the approaches (c) and (d). In Procedure 2, we did not elaborate how to handle non-convergence, but added a penalty term to the log-likelihood, which causes the extra εI n in the Newton-Raphson update. Since we solve the update equation via Cholesky factorizations, adding a small term on the diagonal ensures that the matrix is indeed numerically decomposable.

Properties

We check the requirement list from Section 1.3: (c) We can sample y ∼ q B via chain rule factorization.

(d) We can exactly evaluate q B (y).

(e) The model q B reproduces the dependency structure of the data X although we cannot explicitly compute the marginal probabilities.

Gaussian copula models

In the preceding sections, we discussed three approaches based on generalized linear models. Now we turn to the second class of models we call copula models. Let h θ be a family of auxiliary distributions on X and τ : X → B d a mapping into the binary state space. We can sample from the copula model

q θ (γ) = τ -1 (γ) h θ (v) dv
by setting y = h(v) for a draw v ∼ h θ from the auxiliary distribution.

Definition

Apparently, non-normal parametric distributions s θ with at most d(d -1)/2 dependence parameters either have a very limited dependence structure or rather unfavourable properties [START_REF] Joe | Families of m-variate distributions with given margins and m (m-1)/2 bivariate dependence parameters[END_REF]. Therefore, the normal distribution

h Σ (v) = (2π) -d/2 |Σ| -1/2 exp(-1 /2 v ⊺ Σ -1 v), with mapping τ : R d → B d τ µ (v) = (1 (∞,µ i ] (v 1 ), . . . , 1 (∞,µ d ] (v d )),
appears to be the natural and almost the only option for h θ . This model has already been discussed repeatedly in the literature [START_REF] Emrich | A method for generating high-dimensional multivariate binary variates[END_REF][START_REF] Leisch | On the generation of correlated artificial binary data[END_REF][START_REF] Cox | On some models for multivariate binary variables parallel in complexity with the multivariate Gaussian distribution[END_REF].

Moments

For I ⊆ D, the cross-moment or marginal probabilities is

m I = ∑ γ∈B d q µ,Σ (1 I , γ D\I ) = ∪ γ∈B d {τ -1 µ (1 I ,γ D\I )} h Σ (v) dv = × i∈I {τ -1 µ i (1)} h Σ (v) dv = × i∈I (-∞,µ i ] h Σ (v) dv,
where we used (3) in the first line. Thus, the first and second moment of q (µ,Σ) are

m i = Φ 1 (µ i ), m i, j = Φ 2 (µ i , µ j ; σ i, j )
where Φ 1 (v i ) and Φ 2 (v i , v j ; σ i, j ) denote the cumulative dis- tribution functions of the univariate and bivariate normal distributions with zero mean, unit variance and correlation coefficient σ i, j ∈ [-1, 1].

Sparse Gaussian copulae

We can speed up the parameter estimation and improve the condition of Σ, if we work with a parsimonious Gaussian copula. We can apply the same criterion we already introduced for the sparse logistic regression model. For ε about 1 /100, we define the index set

I de f = {i = 1, . . . , d | xi / ∈ (ε, 1 -ε) }.
which identifies the components which have a marginal probability close to either boundary of the unit interval.

We do not fit a any correlation parameters for the components in I but set σ i, j = 0 for all j ∈ D\{i}. Firstly, the corre- lation does not really matter if the marginal probability is excessively small or large. Secondly, we fit the parameter Σ by separately adjusting the bivariate correlations σ i, j , and com- ponents with high correlations and extreme marginal probability lower the chance that Σ is positive definite.

For the remaining components I c = D \ I, we construct parsimonious Gaussian copula. For δ about 1 /10, we define the association set

A de f = {i, j} ∈ I c × I c | δ < r i, j , i = j
which identifies the components with significant correlation. For i, j ∈ D × D \ L we also set σ i, j = 0 to accelerate the estimation procedure.

Fitting the model

We fit the model q (µ,Σ) to the data by adjusting µ and Σ to the sample moments. Precisely, we solve the equations

Φ 1 (µ i ) = xi , i ∈ D (19) Φ 2 (µ i , µ j ; σ i, j ) = xi, j , (i, j) ∈ A (20)
with sample mean xi and xi, j as defined in (1). We easily solve ( 19) by setting

µ i = Φ -1 1 ( xi ), i ∈ D.
The difficult task is computing a feasible correlation matrix from (20). Recall the standard result (Johnson et al., 2002, p.255)

∂ Φ 2 (y 1 , y 2 ; σ ) ∂ σ = h σ (y 1 , y 2 ), (21)
where h σ denotes the density of the bivariate normal distribution. We obtain the following Newton-Raphson iteration

α r+1 = α r - Φ 2 (µ i , µ j ; α r ) -xi, j h α r (µ i , µ j ) , (i, j) ∈ A, (22)
starting at some α 0 ∈ (-1, 1). We use a fast series approx- imation [START_REF] Drezner | On the computation of the bivariate normal integral[END_REF][START_REF] Divgi | Computation of univariate and bivariate normal probability functions[END_REF] to evaluate Φ 2 (µ i , µ j ; α). These approximations are critical when α r comes very close to either boundary of [-1, 1]. The Newton iteration might repeatedly fail when restarted at the corresponding boundary r 0 ∈ {-1, 1}. This is yet another reason why it is preferable to work with a sparse Gaussian copula. In any event, Φ 2 (y 1 , y 2 ; σ ) is monotonic in σ since (21), and we can switch to bi-sectional search if necessary.

Procedure 3 Fitting the dependency matrix

Input: xi , xi, j for all i, j ∈ D µ i = Φ -1 ( xi ) for all i ∈ D Σ = I d for (i, j) ∈ A do repeat σ [r+1] i, j ← σ [r] i, j - Φ 2 (µ i , µ j ; σ [r] i, j ) -xi, j h σ [r] i, j (µ i , µ j ) until |σ [r] i, j -σ [r-1] i, j | < 10 -3 end for if not Σ ≻ 0 then Σ ← (Σ + |λ | I d )/(1 + |λ |) return µ, Σ
A rather discouraging shortcoming of the Gaussian copula model is that locally fitted correlation matrices Σ might not be positive definite for d ≥ 3. This is due to the fact that an elliptical copula, like the Gaussian, can only attain the bounds (5) for d < 3, but not for higher dimensions.

We propose two ideas to obtain an approximate, but feasible parameter:

(1) We replace Σ by Σ * = (Σ + |λ | I)/(1 + |λ |), where λ is the smallest eigenvalue of the dependency matrix Σ. This approach evenly lowers the local correlations to a feasible level and is easy to implement on standard software. Alas, we make an effort to estimate d(d -1)/2 dependency parameters, and in the end we might not get more than an product model.

(2) We can compute the correlation matrix Σ * which minimizes the distance Σ * -Σ F , where

A 2 F = tr [AA ⊺ ].
In other words, we construct the projection of Σ into the set of correlation matrices. [START_REF] Higham | Computing the nearest correlation matrix -a problem from finance[END_REF] proposes an Alternating Projections algorithm to solve nearestcorrelation matrix problems. Yet, if Σ is rather far from the set of correlation matrices, computing the projection is expensive and, according to our experience, leads to troublesome distortions in the correlation structure.

Properties

We check the requirement list from Section 1.3: (a) The Gaussian copula model is sufficiently parsimonious with dim(θ ) = d(d + 1)/2.

(b) We can fit the parameters µ and Σ via method of moments. The parameter Σ is not always be positive definite which might require additional effort it feasible.

(c) We can sample y ∼ q (µ,Σ) using y = τ µ (v) with v ∼ h Σ .

(d) We cannot evaluate q B (y) since this requires computing a high-dimensional integral expression.

(e) The model q (µ,Σ) exactly reproduces the mean and correlation structure of the data X.

We cannot use the Gaussian copula model in the context of Importance Sampling or Markov chain Monte Carlo, since evaluation of q (µ,Σ) (y) is not possible. This model can be quite useful, however, in other adaptive Monte Carlo algorithms, for instance the Cross-Entropy method [START_REF] Rubinstein | Optimization of computer simulation models with rare events[END_REF] for combinatorial optimisation.

Poisson reduction models

Let N = {1, . . . , n} denote another index set with n ≫ d. Approaches to generating binary vectors that do not rely on the chain rule factorization ( 6) are usually based on combinations of independent random variables

v = (v 1 , . . . , v n ) ∼ ⊗ k∈N h θ k .
We define index sets M = {S i ∈ N | i ∈ D} and generate the entry y i via τ i :

X |S i | → {0, 1} , τ i (v) = f (v S i ), i ∈ D.
In the context of Gaussian copulae, the auxiliary distributions h θ k = h θ are d independent standard normal variables. [START_REF] Park | A Simple Method for Generating Correlated Binary Variates[END_REF] propose the following model based on sums of independent Poisson variables.

Definition

We define a Poisson model q (S ,λ ) with auxiliary distribution

h λ (v) = ∏ k∈N (λ v k k e -λ k )/v k ! and mapping τ : N n 0 → B d τ S (v) = (1 {0} (∑ k∈S 1 v k ), . . . , 1 {0} (∑ k∈S d v k )).

Moments

For an index set I ∈ D, the cross-moments or marginal probabilities are

m I = P ∀i ∈ I : ∑ k∈S i v k = 0 = exp(-∑ k∈∩ i∈I S i λ k ).
Therefore, fitting via method of moments is possible.

Proposition 5. For γ ∈ B d , define the index sets

D 0 = {i ∈ D | γ i = 0} , D 1 = {i ∈ D | γ i = 1} ,
q (S ,λ ) (γ) = ∑ v∈τ -1 (γ) h λ (v) = m D 0 1 -∑ |D 0 | t=1 (-1) t-1 ∑ I⊆I t exp(-∑ k∈∩ i∈I S i \∪ j∈D 1 S j λ k ) .
For a proof see Appendix 10.5.

Fitting the model

We need to determine the family of index sets M and the Poisson parameters λ = (λ 1 , . . . , λ n ) such that the resulting model q (S ,λ ) is optimal in terms of distance to the mean and correlation. Obviously, we face a rather difficult combinatorial problem. [START_REF] Park | A Simple Method for Generating Correlated Binary Variates[END_REF] describe a greedy algorithm, based on convolutions of Poisson variables, that finds at least some feasible combination of S and λ .

Properties

We check the requirement list from Section 1.3:

(a) The Poisson reduction model is not necessarily parsimonious. The number of parameters dim(θ ) is determined by the fitting algorithm.

(b) We fit the model via method of moments using a fast but non-optimal greedy algorithm.

(c) We sample y ∼ q (S ,λ ) using y = τ S (v) with v ∼ h λ .

(d) We cannot evaluate q (S ,λ ) (y) since it requires summation of 2 d-|y| -1 terms using an inclusion-exclusion principle which is computationally not feasible.

(e) The model q (S ,λ ) reproduces the mean and certain correlation structures of the data X. We cannot sample negative correlations.

Since the model is limited to certain patterns of non-negative correlations, we cannot use it as general-purpose model in adaptive Monte Carlo algorithms. It might be useful, however, if we know that the target distribution π has strictly non-negative correlations.

9 Archimedean copula models [START_REF] Genest | A primer on copulas for count data[END_REF] discuss in detail the potentials and pitfalls of applying copula theory, which is well developed for bivariate, continuous random variables, to multivariate discrete distribution. Yet, there have been earlier attempts to sample binary vectors via copulae: [START_REF] Lee | Generating Random Binary Deviates Having Fixed Marginal Distributions and Specified Degrees of Association[END_REF] describes how to construct an Archimedean copula, more precisely the Frank family, (see e.g. Nelsen (2006, p.119)), for sampling multivariate binary data. Unfortunately, most results in copula theory do not easily extend to high dimensions. Indeed, we need to solve a non-linear equation for each component when generating a random vector from the Frank copula, and [START_REF] Lee | Generating Random Binary Deviates Having Fixed Marginal Distributions and Specified Degrees of Association[END_REF] acknowledges that this is only applicable for d ≤ 3. For low-dimensional problems, however, we can just enumerate the solution space B d and draw from an alias table [START_REF] Walker | An efficient method for generating discrete random variables with general distributions[END_REF], which somewhat renders the Archimedean copula approach an interesting exercise, but without much practical value in Monte Carlo applications.

. Recall that I = ∪ k∈D {{i 1 , . . . , i k } ⊆ D | i 1 < • • • < i k } and v I (γ) = ∏ i∈I [(γ i -m i )/ m i (1 -m i )] with m i > 0 for all i ∈ D. We define an inner product ( f , g) de f = E B m ( f (γ)g(γ)) = ∑ γ∈B d f (γ)g(γ) ∏ i∈D m γ i i (1 -m i ) 1-γ i
on the vector space of real-valued functions on B d . The set S = {v I (γ) | I ∈ I } is orthonormal, since

(v I , r J ) = ∏ i∈I∩J E B m (γ i -m i ) 2 m i (1 -m i ) ∏ i∈(I∪J)\(I∩J) E B m γ i -m i m i (1 -m i ) = 0 for I = J 1 for I = J,
= ( f , 1) + ∑ I∈I v I ( f , v I ). Since (π/B m , v I ) = ∑ γ∈B d (π(γ)/B m)(γ) v I (γ)B m (γ) = E π (v I (γ)) = c I , we obtain π(γ)/B m (γ) = 1 + ∑ I∈I v I (γ) c I for f = π/B m which concludes the proof.

Proposition 2

Proof. We first derive two auxiliary results to structure the proof.

Lemma 1. For a set I ⊆ D of indices it holds that

∑ γ∈B d ∏ k∈I∪{i, j} γ k = 2 d-|I|-2+1 I (i)+1 I∪{i} ( j) .
(i) + 1 D\(I∪{i}) ( j) = |I| + 2 -1 I (i) -1 I∪{i} ( j) which implies Lemma 1. Lemma 2. ∑ i∈D ∑ j∈D 2 1 I (i)+1 I∪{i} ( j) a i, j = 1 ⊺ A1 + tr [A] + ∑ i∈I 2 ∑ j∈D a i, j + ∑ j∈I\{i} a i, j
Straightforward calculations:

2 1 I (i)+1 I∪{i} ( j) = (1 + 1 I (i))(1 + 1 I∪{i} ( j)) = (1 + 1 I (i))(1 + 1 I ( j) + 1 {i} ( j) -1 I∩{i} ( j)) = 1 + 1 I (i) + 1 I ( j) + 1 I (i)1 I ( j) + 1 {i} ( j) + 1 I (i)1 {i} ( j) -1 I∩{i} ( j) -1 I (i)1 I∩{i} ( j) = 1 + 1 {i} ( j) + 1 I (i) + 1 I ( j) + 1 I×I (i, j) -1 I∩{i} ( j),
where we used 1

I (i)1 {i} ( j) = 1 I (i)1 I (i)1 {i} ( j) = 1 I (i)1 I ( j)1 {i} ( j) = 1 I (i)1 I∩{i} ( j)
in the second line. Thus, we have

∑ i∈D ∑ j∈D 2 1 I (i)+1 I∪{i} ( j) a i, j = ∑ i∈D ∑ j∈D 1 + 1 {i} ( j) + 1 I (i) + 1 I ( j) + 1 I×I (i, j) -1 I∩{i} ( j) a i, j = ∑ i∈D ∑ j∈D a i, j + ∑ j∈D a j, j + ∑ i∈I ∑ j∈D a i, j + ∑ i∈D ∑ j∈I a i, j + ∑ i∈I ∑ j∈I a i, j -∑ i∈I a j, j = 1 ⊺ A1 + tr [A] + ∑ k∈I 2 ∑ l∈D a k,l + ∑ l∈I a k,l -a k,k = 1 ⊺ A1 + tr [A] + ∑ k∈I 2 ∑ l∈D a k,l + ∑ l∈I\{k} a k,l
The last line is the assertion of Lemma 2.

Using the two Lemmata above, we find a convenient expression for the cross-moment

m I = ∑ γ∈B d (∏ k∈I γ k ) µ(a 0 + γ ⊺ Aγ) = µ ∑ γ∈B d a 0 + ∑ γ∈B d (∏ k∈I γ k ) ∑ i∈D ∑ j∈D γ i γ j a i, j = µ 2 d-|I| a 0 + ∑ i∈D ∑ j∈D a i, j ∑ γ∈B d (∏ k∈I∪{i, j} γ k ) (Lemma 1)
= µ 2 d-|I| a 0 + ∑ i∈D ∑ j∈D 2 d-|I∪{i, j}| a i, j = µ2 d-|I|-2 4a 0 + ∑ i∈D ∑ j∈D 2 1 I (i)+1 I∪{i} ( j) a i, j (Lemma 2) = µ2 d-|I|-2 4a 0 + 1 ⊺ A1 + tr [A] + ∑ i∈I 2 ∑ j∈D a i, j + ∑ j∈I\{i} a i, j

Since m / 0 = 1 by definition, we the normalization constant is µ = 2 -d+2 (4a 0 + 1 ⊺ A1 + tr [A]) -1 , which allows us to write down the normalized cross-moments

m I = 1 2 |I| + ∑ i∈I 2 ∑ j∈D a i, j + ∑ j∈I\{i} a i, j 2 |I| (4a 0 + 1 ⊺ A1 + tr [A]
) .

The proof is complete.

Proposition 3

Proof. We easily margin out the last component d, having I = {1, . . . , d -t}, q (d-1)

A,a 0 (γ I )µ -1 = q (d)

A,a 0 (γ I , 1) + q = 4 ∑ i∈I ∑ j∈I γ i γ j a i, j + 4 ∑ i∈I ∑ j∈I c γ i a i, j + ∑ i∈I c ∑ j∈I c a i, j + ∑ i∈I c a i,i = 4 ∑ i∈I γ i (∑ j∈I γ j a i, j + ∑ j∈I c a i, j ) + ∑ i∈I c ∑ j∈I c a i, j + ∑ i∈I c a i,i

The proof is complete.

Proposition 4

Proof. For convenience of notation, let γ -= (γ 1 , . . . , γ d-1 ). Note that q A (γ) = µ exp(γ ⊺ -A ′ γ -+ γ d (2b ⊺ γ -+ c)). The marginal distribution is therefore

π(γ -) = µ exp(γ ⊺ -A ′ γ -) 1 + exp(2γ ⊺ -b + c) = µ exp γ ⊺ -A ′ γ -+ γ ⊺ -b + c /2 exp(-γ ⊺ -b -c /2) + exp(γ ⊺ -b + c /2) = µ exp γ ⊺ -A ′ γ -+ γ ⊺ -b + c /2 2 cosh γ ⊺ -b + c /2 .
The marginal log mass function is thus log π(γ -) = log(2µ)

+ c /2 + γ ⊺ -A ′ γ -+ γ ⊺ -b + logcosh γ ⊺ -b + c /2 .
For log cosh we can use a Taylor approximation log cosh(γ ⊺ -b + c /2) ≈ log cosh( c /2) + γ ⊺ -b tanh( c /2) + 1 /2 (γ ⊺ -b) 2 sech 2 ( c /2) to obtain log π(γ -) ≈ log(2µ cosh( c /2))

+ c /2 + γ ⊺ -A ′ γ -+ 1 + tanh( c /2) γ ⊺ -b + 1 /2 sech 2 ( c /2)(γ ⊺ -b) 2
Since γ -is a binary vector, we have γ ⊺ -b = γ ⊺ -diag [b] γ -and can thus rewrite the inner products as

γ ⊺ -A ′ γ -+ γ ⊺ -b + (γ ⊺ -b) 2 = tr A ′ γ -γ ⊺ -+ diag [b] γ -γ ⊺ -+ bb ⊺ γ -γ ⊺ - = γ ⊺ -(A ′ + diag[b] + bb ⊺ )γ -.
We let denote µ * = 2µ cosh( c /2) exp( c /2) = µ(exp(c /2) + exp( c /2)) exp( c /2) = µ (1 + exp(c))

and

A * = A ′ + (1 + tanh( c /2)) diag [b] + 1 /2 sech 2 ( c /2)bb ⊺
to form the approximation π(γ -) ≈ µ * exp(γ ⊺ -A * γ -) which completes the proof.

Proposition 5

Proof. Straightforward calculations using the inclusion-exclusion principle for the union of events:

q (S ,λ ) (γ) = ∑ v∈τ -1 (γ) h λ (v) = P h λ ∩ i∈D 1 {0} ∑ k∈S i v k = γ i = P h λ ∩ i∈D 1 ∩ k∈S i {v k = 0} , ∩ i∈D 0 ∪ k∈S i {v k > 0} = P h λ ∩ i∈D 1 ∩ k∈S i {v k = 0} P h λ ∩ i∈D 0 ∪ k∈S i \∪ j∈D 1 S j {v k > 0} = P q (S ,λ ) (γ D 1 = 1) 1 -P h λ ∪ i∈D 0 ∩ k∈S i \∪ j∈D 1 S j {v k = 0} = m D 0 1 -∑ |D 0 | t=1 (-1) t-1 ∑ I⊆I t P h λ ∩ i∈I ∩ k∈S i \∪ j∈D 1 S j {v k = 0} = m D 0 1 -∑ |D 0 |
t=1 (-1) t-1 ∑ I⊆I t exp(-∑ k∈∩ i∈I S i \∪ j∈D 1 S j λ k ) .

The proof is complete.

  (a) The log-linear model is sufficiently parsimonious with dim(θ ) = d(d + 1)/2. (b) We can fit the parameter A via minimum least squares.

  (a) The logistic regression model is sufficiently parsimonious with dim(θ ) = d(d + 1)/2. (b) We can fit the parameters b i via likelihood maximisation for all i ∈ D. The fitting is computationally intensive but feasible.

  and the families of subsets I t = {I ∈ D 1 | |I| = t}. We can write the mass function of the Poisson model as

For

  an index set M ⊆ D, we have the sum formula ∑ γ∈B d ∏ k∈M γ k = 2 d-|M| . If we have an empty set M = / 0 the sum equals 2 d and each time we add a new index i ∈ D \ M to M half of the addends vanish. The number of elements in M = I ∪ {i, j} is the number of elements in I plus one if i / ∈ I and again plus one if i = j and j / ∈ I. Written using indicator function, we have |I ∪ {i, j}| = |I| + 1 D\I

  0 (γ I , 0) µ -1 = 2a 0 + (γ I , 1) ⊺ A(γ I , 1) + (γ I , 0) ⊺ A(γ I , 0) = 2a 0 + tr [A [(γ I , 1)(γ I , 1) ⊺ + (γ I , 0)(γ I , 0) ⊺ ]] = 2a 0 + tr A 2γ I γ ⊺ ,we obtain for I = {1, . . . , d -t} and I c = D \ I q (d-t)A,a 0 (γI )µ -1 = 2 t a 0 + 2 t-2 tr A t + I t = tr [A [(2 γ I , 1 t )(2 γ I , 1 t ) ⊺ + diag[0 I , 1 t ]]] = [(2 γ I , 1 t ) ⊺ A(2 γ I , 1 t ) + tr [Adiag [0 I , 1 t )]]]

  There are 2 d -1 elements in S and B m (γ) > 0 which implies that S ∪ {1} is an orthonormal basis of the real-valued function on B d . It follows that each function f : B d → R has exactly one representation as linear combination of functions in S ∪ {1} which is f