arXiv:1706.06182v2 [math.PR] 30 Jun 2017

BERNOULLI CORRELATIONS
AND
CUT POLYTOPES
MARK HUBER AND NEVENA MARIĆ

Abstract. Given n symmetric Bernoulli variables, what can be said
about their correlation matrix viewed as a vector? We show that the
set of those vectors R(Bn ) is a polytope and identify its vertices. Those
extreme points correspond to correlation vectors associated to the discrete uniform distributions on diagonals of the cube [0, 1]n . We also
show that the polytope is affinely isomorphic to a well-known cut polytope CUT(n) which is defined as a convex hull of the cut vectors in a
complete graph with vertex set {1, . . . , n}. The isomorphism is obtained
explicitly as R(Bn ) = 1 − 2 CUT(n). As a corollary of this work, it is
straightforward using linear programming to determine if a particular
correlation matrix is realizable or not. Furthermore, a sampling method
for multivariate symmetric Bernoullis with given correlation is obtained.
In some cases the method can also be used for general, not exclusively
Bernoulli, marginals.
MSC 52B12, 60E05, 62H20
Keywords: Bernoulli distribution, Extreme Correlations, Cut Polytopes.

1. Introduction
Consider the question of admissible correlations among n random variables (X1 , . . . , Xn ) for which the marginal distributions are known. This
topic has a long history, dating partly back to the work of de Finetti [5]
where the problem of maximum negative achievable correlation among n
random variables was studied. The general form of the problem was studied
by Fréchet [10] and Hoeffding [14] in a body of work which grew out of the
problem originally posed by Lévy [20].
The big question is: can we completely describe set of correlation matrices
for a given set of marginal distributions? For n = 2 the answer is completely
known in terms of Fréchet-Hoeffding bounds, so the interesting question is
what happens in higher dimensions. In this work we give a partial answer to
this question, showing that if a particular vector calculated from the target
correlations and marginals falls into the CUT(n) polytope, then there does
exist such a joint distribution. This condition is both necessary and sufficient
in the case of symmetric Bernoulli marginals.
It is well known that a correlation matrix is symmetric positive semidefinite and has all diagonal elements equal to 1. The set of all such matrices
1

2

MARK HUBER AND NEVENA MARIĆ

of order n will be denoted by En×n . This convex compact set is sometimes
called the elliptope, a term coined by Laurent and Poljak [18].
For Gaussian marginals, the entirety of En×n can be realized, but this is
the only nontrivial set of marginals for which the question has been settled.
Surprisingly enough, for other common distributions very little is known. A
reason is perhaps that much of correlation theory has been constructed using
normal marginals [16]. In spite of wide usage of the correlation coefficient, its
bounds in a multivariate non-Gaussian setting has been mostly unexplored.
One case that has been partially explored is that of copulas. A probability
measure µ on [0, 1]n is a copula if all its marginals are uniformly distributed
on [0, 1]. In a recent work of Devroye and Letac [6] it has been shown that
every element in En×n is a correlation matrix for some copula, for n ≤ 9,
but the authors believe that the statement does not hold for n ≥ 10.
In this paper we focus on symmetric Bernoulli variables for multiple reasons. First, it is the simplest distribution, with equally likely binary outcomes. We say that a random variable X has the Bernoulli distribution
(write X ∼ Bern(p)) if P(X = 1) = p and P(X = 0) = 1 − p. The symmetric
Bernoulli distribution is the case when p = 1/2.
The second reason comes from Huber and Marić [15] where this distribution was shown to be in a certain sense the most difficult problem: for general
marginals and some correlations it is possible to transform the problem into
symmetric Bernoulli marginals.
To be precise, the problem is to simulate (X1 , . . . , Xn ) where the correlation between each Xi and Xj are specified along with the marginal distribution of each Xi , or to determine that no such random variables exist.
This problem, in different guises, appears in numerous fields: physics [25],
engineering [17], ecology [8], and finance [19], to name just a few. Due to
its applicability in the generation of synthetic optimization problems, it has
also received special attention by the simulation community [13], [12].
An excellent overview of the developments in the field of generating multivariate probability densities with pre-specified margins can be found in
Dall’Aglio et al. [4], Rüschendorf et al. [24], and Conway [3].
Consider a matrix Σ that is potentially a correlation matrix for a particular choice of marginal distributions. Then what the method of [15] does is
build a second correlation matrix ΣB such that if is possible to have a multivariate distribution with symmetric Bernoulli marginals and correlation
ΣB , then it is possible to build a multivariate distribution with the original
marginals and correlation matrix Σ.
The question of existence for general marginals then becomes the question
of which matrices can be realized as correlation matrices for (B1 , . . . , Bn )
which are all marginally Bern(1/2).
It should be noted that the answer for symmetric Bernoulli marginals
will be a strict subset of En×n , even when n is small. As a simple example

BERNOULLI CORRELATIONS

consider

AND

CUT POLYTOPES

3




1
−0.4 −0.4
−0.4
1
−0.4 .
−0.4 −0.4
1
While this matrix is in the elliptope E3×3 , it cannot (see [15]) be the correlation matrix of three random variables with symmetric Bernoulli marginals.
Considering the applications, there is an obvious need to understand what
are the theoretically attainable correlations, so that inference using the correlations drawn from data can be done properly. Chaganty and Joe [2]
claim that there were errors in reports of efficiency calculations for generalized estimating equations, such as Table 1 of Liang and Zeger [22], and with
the analysis of real-life binary data using the current generalized estimating
equations software, caused by the belief that any matrix in En×n is a possible
correlation matrix for a set of binary random variables. In the same paper
they were able to characterize the achievable correlation matrices when the
marginals are Bernoulli. When the dimension is 3 their characterization is
easily checkable (as for the 3 by 3 matrix given above), in higher dimensions
they give a number of inequalities that grows exponentially in the dimension. They also give an approximate method for checking attainability of
the correlation matrix in higher dimensions.
In this paper we give a complete characterization of the correlation matrices for multivariate symmetric Bernoulli distributions by relating them to
the well-known CUT polytope. This approach leads also to a novel sampling
method from the desired marginals and correlations.
Let Bn be a set of all n-variate symmetric Bernoulli distributions, En =
{(i, j) : 1 ≤ i 6= j ≤ n}, and R : Bn → [−1, 1]En the correlation mapping.
Here in place of a correlation matrix we focus on the correlation vector which
contains elements above (or below) the diagonal of a correlation matrix
placed in the same order, row by row. We show that R(Bn ) is a polytope
and identify its vertices.
Let πi for i in {1, . . . , 2n−1 } be the uniform distribution over the end
points of the i-th diagonal of the n-dimensional cube [0, 1]n . Then our main
results is as follows.
Theorem 1. R(Bn ) is a polytope with vertices {R(π1 ), ..., R(π2n−1 )}. That
is, a vector ρ ∈ [−1, 1]En is a correlation vector for some distribution µ ∈ Bn
if and only if it can be written as a convex combination of {R(π1 ), ..., R(π2n−1 )}.
We also uncover a striking relation between this polytope and the cut
polytope CUT(n), which is defined as a convex hull of the cut vectors in
a complete graph with vertex set {1, . . . , n}. The cut polytopes play an
important role in combinatorial optimization, as they can be used to formulate the max-cut problem, which has many applications in various fields [7],
[28]. Vertices of the CUT are all 0-1 vectors, meaning that each coordinate
is either 0 or 1. A relation between the polytopes is given the following
theorem.

4

MARK HUBER AND NEVENA MARIĆ

Theorem 2. For ρ ∈ [−1, 1]En : ρ ∈ R(Bn ) if and only if 1 − 2ρ ∈ CUT(n).
The remainder of the text is organized in the following way. In Section
2 we introduce the problem of Bernoulli correlations via their agreement
probabilities, starting off from a prior work [15] and study distributions πi
as (some) extreme points of Bn . Cut polytopes are introduced in Section
3, where also an important theorem of Avis [1] is analyzed. The theorem
provides a certain probabilistic context of cut polytopes, and we are able
to extend it in connection to Bn . In Section 4 we derive our main results.
The asymmetric Bernoulli case is also discussed here. Section 5 is focused
on applications and how our findings can be used in practice via linear
programming. A sampling method for n-variate symmetric Bernoullis with
given correlation is outlined as well. This section also contains a worked
through example using difficult marginals. Finally in Section 6 we discuss
our results in a larger context.
2. Bernoulli correlations and agreement probabilities
Definition 1. An n-variate symmetric Bernoulli distribution µ is a probability measure on {0, 1}n such that all the marginals are Bern(1/2), that
is
X
1
for k = 1, . . . , n.
µ(x) =
2
n
x∈{0,1} : x(k)=0

Let Bn be a set of all such measures.
Remark 1. Clearly Bn is a convex set: if µ1 , µ2 ∈ Bn then, for any α ∈
(0, 1), αµ1 + (1 − α)µ2 ∈ Bn . Since it is described by a finite set of linear
n
equalities, Bn is closed subset of R2 and is bounded. Moreover Bn is a
polytope [6].
Definition 2. If (B1 , . . . , Bn ) ∼ µ, µ ∈ Bn then for all i, j ∈ {1, . . . , n},
the probability that Bi and Bj have the same value is called the agreement
probability, and denoted as λµ (i, j) = Pµ (Bi = Bj ).
The map that takes µ to λµ is denoted by Λ.
Now consider the correlation between two Bern(1/2) random variables:
ρµ (i, j) = corr(Bi , Bj ) = 4Eµ Bi Bj − 1.
Remark 2. Correlation is related to the agreement probability between the
variables in a linear way:
1
λµ (i, j) = (1 + ρµ (i, j)).
2
In [15] Huber and Marić studied elements of Bn via their agreement probabilities and they were able to provide necessary and sufficient conditions
for an agreement matrix to be attainable up to the dimension 4. It can be
verified that those conditions, given in Theorem 3 of [15], place λ in the
following polytopes:

BERNOULLI CORRELATIONS

AND

CUT POLYTOPES

5

• n = 2: Interval [0, 1].
• n = 3: Tetrahedron with vertices (1, 1, 1), (1, 0, 0), (0, 1, 0), (0, 0, 1).
• n = 4: (λ12 , λ13 , λ14 , λ23 , λ24 , λ34 ) belongs to a 6-polytope with 8
vertices:(1, 1, 1, 1, 1, 1), (0, 1, 1, 0, 0, 1), (1, 1, 0, 1, 0, 0), (1, 0, 0, 0, 0, 1),
(0, 1, 0, 0, 1, 0), (0, 0, 1, 1, 0, 0), (0, 0, 0, 1, 1, 1), (1, 0, 1, 0, 1, 0). Here
each appropriate 3-dimensional projection (that corresponds to agreement probabilities in any given subset of 3 variables; for example
along first, second, and fourth coordinate) bears a tetrahedron from
the case n = 3, as it should be.
These results stem from [15] but also follow directly from Theorem 1.
2.1. Diagonals. Let x = (x1 , . . . , xn ) ∈ {0, 1}n . We can think of x as a
vertex of the n-dimensional cube [0, 1]n . Let 1 = (1, . . . , 1). The diagonal
associated to the vertex x (and its “opposite” 1 − x) is set
Dx := {(x1 , . . . , xn ), (1 − x1 , . . . , 1 − xn )}.
Definition 3. For every x ∈ {0, 1}n , the (discrete) uniform distribution
over the diagonal Dx is denoted by πx . That is πx (x) = πx (1 − x) = 1/2.
Remark 3. Note that πx ∈ Bn and also πx = π1−x .
The number of vertices x is finite, so they can be ordered. Every vertex
of the n-dimensional cube [0, 1]n can be mapped to a positive integer via
b : {0, 1}n → N:
n
X
xj 2n−j .
b(x1 , x2 , . . . , xn ) = 1 +
j=1

Each set Dx contains one vector y whose first coordinate is 0, and for this
vector b(y) ∈ {1, . . . , 2n−1 }. To simplify the notation, for k ∈ {1, . . . , 2n−1 }
let
πk := the uniform distribution over Dx such that ∃y ∈ Dx with b(y) = k .
n−1

It should be noted that {πx }x∈{0,1}n and {πk }2k=1 refer essentially to the
same set of distributions, with a difference being that in the ordered notation
there are no equal elements, e.g. π1 = π0 = π1 . However, both notations
are useful, and we will be using them concurrently throughout the text.
Before we state the first result highlighting the importance of distributions
π, a few additional definitions are needed. As in Grünbaum [11], v ∈ P is an
extreme point of the polytope P if (∀a, b ∈ P )(∃β ∈ (0, 1))(βa + (1 − βb) =
v ⇒ a = b = v). Let ext(P ) denote the extreme points of the polytope P .
Proposition 1. For every k = 1, . . . , 2n−1 , a measure πk is an extreme
point of Bn . That is {π1 , . . . , π2n−1 } ⊆ ext(Bn ).
Proof. Let x, y ∈ {0, 1}n and πx be a convex linear combination of µ and
ν in Bn . Then for y 6= x, 1 − x, πx (y) = 0 and therefore µ(y) = ν(y) = 0.
Suppose µ(x) = p and µ(1 − x) = 1 − p, for 0 < p < 1. As µ has Bern(1/2)

6

MARK HUBER AND NEVENA MARIĆ

marginals, it follows that p = 1/2 and therefore µ = πx . Analogously, ν is
also equal to πx and hence, πx is an extreme point of Bn .

3. CUT polytopes
Let G = (V, E) be a graph with vertex set V and edge set E. For S ⊆ V
a cut of the graph is a partition (S, S C ) of the vertices. The cut-set consists
of all edges that connect a node in S to a node not in S.
Let Vn = [n] = {1, . . . , n}, En = {(i, j); 1 ≤ i 6= j ≤ n}, and Kn =
(Vn , En ) be a complete graph with the vertex set [n].
Definition 4. For every S ⊆ [n] a vector δ(S) ∈ {0, 1}En , defined as
δ(S)ij =



1,
0,

if |S ∩ {i, j}| = 1
otherwise

for (1 ≤ i < j ≤ n), is called a cut vector of Kn .
The cut polytope CUT(n) is the convex hull of all cut vectors of Kn :
CUT(n) = conv{δ(S) : S ⊆ [n]}.
Remark 4. Since every cut vector is a vertex of CUT(n), there are 2n−1
vertices of this polytope [28].
Each δ(·) is a 0/1-vector (every coordinate value is either 0 or 1). The
convex hulls of finite sets of 0/1-vectors are called 0/1-polytopes, out of
which cut polytopes are a sub-class. An excellent lecture on 0/1-polytopes,
including CUT, is given in Ziegler [28]. More thorough treatment of cut
polytopes can be found in Deza and Laurent [7]. The cut polytopes play
an important role in combinatorial optimization, as they can be used to
formulate the max-cut problem, which has many applications in various
fields, like statistical physics, in relation to spin glasses [7].
All the symmetries of the cube [0, 1]n transform 0/1-polytopes into 0/1polytopes [28]. In particular, a symmetry is obtained by replacing some
coordinates xi by 1−xi , which is called switching. Two 0/1-polytopes P and
P’ are 0/1- equivalent if a sequence of switching and coordinate permuting
operations can transform P into P’. We define now a polytope that is 0/1equivalent to CUT(n) and also useful for our further analysis.
Definition 5. CUT∗ (n) is a polytope obtained by applying 0/1 switching
operation to all vertices of CUT(n).
A relation of cut polytopes with probability spaces is given in the following
theorem. We cite here a version of the theorem given in the book by Deza
and Laurent [7], but the authorship dates back to Avis [1].

BERNOULLI CORRELATIONS

AND

CUT POLYTOPES

7

Theorem 3 (Avis [1]). For a vector d ∈ REn the following statements are
equivalent
1) d ∈ CUT(n)
2) There exists a probability space (Ω, A, ν) and events A1 , . . . , An ∈ A
such that dij = ν(Ai △Aj ) for all 1 ≤ i < j ≤ n.
Here A△B = (A \ B) ∪ (B \ A) denotes the symmetric difference between
sets A and B.
This theorem is useful for us for the following reason: to every event
in a probability space one can assign a Bernoulli random variable, as its
indicator. In order to obtain a symmetric Bern, the probability of the event
has to be equal to 1/2. As we are going to show, that is exactly what
happens in the setting of the above theorem.
3.1. On Theorem 3. We look into the proof of Theorem 3, given in [7]
Proposition 4.2.1., of the part 1) ⇒ 2) and will show that events Ai , as
defined in that proof, have toP
be of probability 1/2.
P
Let d ∈ CUT(n), then d = S⊆[n] βS δ(S) for some βS ≥ 0 and βS = 1.
Note that δ(S) = δ(S c ) and that this vector ultimately participates only
once in the convex representation of d, so, without loss of generality, one
can assume that βS = βS c . Now, the probability space (Ω, A, ν) is defined
as follows. Let Ω = {S : S ⊆ [n]}, A the family of subsets of Ω and ν a
probability measure, for A ∈ A defined by
X
ν(A) =
βS .
S∈A

Setting Ai = {S ∈ Ω : i ∈ S} one obtains events with the desired property,
namely for which dij = ν(Ai △Aj ). In order to calculate ν(Ai ) we observe
that
X
X
(1)
βS =
βS 1(i ∈ S)
ν(Ai ) =
S∈Ai

S⊆[n]

and for its complement Aci = {S ∈ Ω : i ∈
/ S}
X
X
X
βS 1(i ∈
/ S) =
βS 1(i ∈ S c )
βS =
ν(Aci ) =
S∈Aci

=

X

S⊆[n]

S⊆[n]

βS c 1(i ∈ S c ) = ν(Ai ).

S⊆[n]

Therefore ν(Ai ) = 1/2 for i = 1, . . . , n and we can formulate the following
corollary.

8

MARK HUBER AND NEVENA MARIĆ

Corollary 1. For a vector d ∈ REn the following statements are equivalent
1) d ∈ CUT(n)
2) There exists a probability space (Ω, A, ν) and events A1 , . . . , An ∈ A
1
such that ν(Ai ) = for all 1 ≤ i ≤ n and dij = ν(Ai △Aj ) for all 1 ≤ i < j ≤ n.
2
Remark 5. In the 2) ⇒ 1) part of the proof of the Theorem 3, the probabilities of the events Ai do not play any role, so that direction of the statement
remains true in the above Corollary.
4. Bernoulli agreements and CUT*
To every event A in a probability space can be associated its indicator
1(A), a random variable that takes binary values: 1 if the event occurs and
0 otherwise. The indicator has Bernoulli distribution with parameter that
equals the probability of the event A.
Suppose now that the condition 2) from Corollary 1 is satisfied and define
Di = 1(Ai ) for i = 1, . . . , n. Clearly Di ∼ Bern(1/2) and for all i, j ∈ En
(2)

dij = ν(Ai △Aj ) = Pν (Di 6= Dj ) = 1 − λν (i, j).

From (2) and Corollary 1 then it follows
(3)

λν (i, j) ∈ CU T ∗ (n).

Proposition 2. CU T ∗ (n) is a polytope with vertex set {λπx : x ∈ {0, 1}n }.
Proof. Let x = (x1 , . . . , xn ) ∈ {0, 1}n . Fix i, j ∈ [n] and consider λπx (i, j)
X
πx (y).
λπx (i, j) =
y∈{0,1}n :yi =yj

If xi = xj , then also 1 − xi = 1 − xj , and the above sum equals to πx (x) +
πx (1 − x) = 1.
Similarly, if xi 6= xj , then 1 − xi 6= 1 − xj and every summand above
equals 0. Hence
(4)

λπx (i, j) = 1(xi = xj ).

To establish a relation between the agreement vector and a cut vector
introduce Sx = {i ∈ [n] : xi = 1}, the set of all coordinates whose value in
x is 1. If exactly one of i, j ∈ [n] is in Sx then xi 6= xj , otherwise xi = xj .
Therefore, from (4) it follows
λπx (i, j) = 1 − 1(xi 6= xj ) = 1 − 1(|Sx ∩ {i, j}| = 1)
Applying Definition 4 we get the desired relationship
(5)

λπx = 1 − δ(Sx ).

∀i, j ∈ [n].

BERNOULLI CORRELATIONS

AND

CUT POLYTOPES

9

Note also that a map from {0, 1}n to P([n]) (the power set of [n]) where
x 7→ Sx is a bijection and that every S ⊆ [n] can be identified as Sx for
some x ∈ {0, 1}n . From the equation (5) and the definition of cut polytope
then we have
CUT(n) = conv{δ(S) : S ⊆ [n]}
= conv{δ(Sx ) : x ∈ {0, 1}n } = conv{1 − λπx : x ∈ {0, 1}n }.
Finally, since CUT∗ (n) is a 0/1 switched image of CUT(n), it follows that
CUT∗ (n) = conv{λπx : x ∈ {0, 1}n }.
As every cut vector is a vertex of CUT(n), then it follows that {λπx : x ∈
{0, 1}n } must be the vertex set of the polytope CUT∗ (n).

The next theorem says the set of all agreement vectors (for n-dim symmetric Bernoulli distributions) is exactly the known polytope CUT∗ (n). Proposition 2 gives a probabilistic interpretation of its extreme points, and consequently, we are able to completely describe every λµ .
Theorem 4. Λ(Bn ) = CUT∗ (n).
Proof. (⇒:) Let µ be a measure in Bn . Then, by (3), λµ ∈ CUT∗ (n) and
by Proposition 2 follows that λµ ∈ conv{λπx : x ∈ {0, 1}n }. Hence
Λ(Bn ) ⊆ CUT∗ (n).
(⇐:) Conversely, suppose that λ ∈ CUT∗ (n) and then we want to show
that there is µ ∈ Bn such that λµ = λ. Recall that for each diagonal
distribution πx there is a unique label k ∈ {1, . . . , 2n−1 } such that πx = πk .
Then write λ as a convex linear combination of the λπk .
λ=

n−1
2X

πk

αk λ ,

k=1

n−1
2X

αk = 1

k=1

Then the distribution is
µ=

n−1
2X

αk πk .

k=1

Being a convex combination of (π1 , . . . , π2n−1 ) ∈ Bn , µ is also in Bn . Fix
i, j ∈ {1, . . . , n} and consider the probability that for (B1 , . . . , Bn ) ∼ µ, that
Bi = Bj . Since µ is a convex mixture of other measures,
µ

λ (i, j) =

n−1
2X

αk λπk (i, j) = λ(i, j).

k=1


The original motivation was to understand the attainable correlation matrices. Recall the correlation mapping R: µ 7→ ρµ , and recall the relation

10

MARK HUBER AND NEVENA MARIĆ

between agreement probability and the correlation, which is obviously a
linear bijection:
R(µ) = 2λµ − 1.
Our original goal was to describe the set R(Bn ) = 2Λ(Bn ) − 1. Since Λ(Bn )
has extreme points corresponding to λπi , so does R(Bn ), proving Theorem 1.
Remark 6. From Proposition 1 we know that for every x ∈ {0, 1}n , πx is
an extreme point of Bn . Since R is an affine transformation we also know
that ext(R(Bn )) ⊆ R(ext(Bn )). Due to Theorem 4, even we do not have a
complete description of Bn (i.e we do not know all its extreme points) we
are still able to completely describe R(Bn ).
Note also that for δ ∈ [0, 1]En , δ ∈ CUT(n) ⇔ 1 − δ ∈ CU T ∗ (n). Now
Theorem 2 follows directly from Theorem 4.
4.1. Relation to the correlation polytope. The polytope R(Bn ) does
describe the set of allowable correlation matrices, but it is not directly a
correlation polytope, which is typically defined (see [9]) as
COR(n) = conv{bbT ∈ Rn×n |b ∈ {0, 1}n }.
However, since there is an affine isomorphism between CUT(n) and COR(n−
1) (see [7]) they can be related in a similar way as we did here.
4.2. Asymmetric n-variate Bernoulli. It should be noted that discovered relation between CUT(n) and Bn does not extend to asymmetric multivariate Bernoulli distributions. It is enough to analyze the bivariate case
with equal marginals.
The correlation between two Bern(p) random variables belongs to the interval [ρmin , 1]. Maximum correlation in case of equal marginals, always
equals to 1 and the minimum correlation ρmin can be calculated using
Fréchet-Hoeffding bounds [10, 14]

−(1 − p)/p, for p ≥ 1/2
ρmin =
−p/(1 − p), for p ≤ 1/2.
It is clear now that only for p = 1/2, ρmin = −1 and possible correlations
equal to the entire interval [−1, 1], while for any other value of p it is a strict
subinterval of [−1, 1]. From the linear relationship between the agreement
probability λ and ρ, it follows that, again only for p = 1/2 it is true 0 ≤ λ ≤
1, while for other values of p these bounds are not sharp. For example, for
p = 3/4, −1/3 ≤ ρ ≤ 1, and 1/3 ≤ λ ≤ 1.
In two dimensional case the cut polytope is known to be CUT(2) = [0, 1] =
CUT∗ (2) so it corresponds to Λ(Bn ) only in the symmetric case.

BERNOULLI CORRELATIONS

AND

CUT POLYTOPES

11

5. Applications
5.1. Determining feasibility. Once the vertices of the polytope R(Bn )
have been determined, it is straightforward using linear programming to
determine if a particular correlation vector ρ : En → [−1, 1] is realizable.
For k ∈ 1, . . . , 2n−1 , let xk (i) be the ith bit in the binary representation
of the number k − 1. Then let
R(πk )(i, j) = 21(xk (i) = xk (j)) − 1.
The goal is to find nonnegative α1 , . . . , α2n−1 such that, for all i < j
α1 R(π1 )(i, j) + · · · + α2n−1 R(π2n−1 )(i, j) = ρ(i, j)
and α1 + · · · + α2n−1 = 1.

Let ρaug = [ρ 1]. These n2 + 1 equations can be written as the system
M α = ρaug subject to α ≥ 0 for appropriate choice of M and ρaug . It can
be determined if this polytope is nonempty by solving the Phase I linear
program
min z1 + · · · + z2n−1
M α + z = ρaug
α, z ≥ 0.
This linear program has a feasible solution of α = 0 and z = ρaug . Then
the original polytope is nonempty if and only if this linear program has a
solution with objective function value 0.
The free lpsolve package within R was used to test how possible this
would be. On a problem with n = 15 it took roughly eight seconds to
solve on an Intel i5@3.30 Ghz. Of course the number of variables is growing
exponentially with n, so this method will grow rapidly as n grows.
5.1.1. Generating multivariate Bern(1/2) with given correlation vector ρ. The vector (α1 , . . . , α2n−1 ) obtained above, using linear programming, can be used to sample (B1 , . . . , Bn ) with desired correlation in the
following way. Let us label the vertices of n-dimensional cube [0, 1]n whose
first coordinate is 0:
vk := {y ∈ {0, 1}n : πk (y) = 1/2 and y(1) = 0}.
Then for k = 1, . . . , 2n−1
(B1 , . . . , Bn ) = vk or 1 − vk with probability αk /2.
5.2. Example. To illustrate how the above algorithm can be used for general margins, consider the following problem with n = 4. Suppose X1 is
uniform over [0, 1], X2 is exponential with mean 2, X3 equals 1 with probability 0.3 and 4 with probability 0.7, and X4 is a standard normal random

12

MARK HUBER AND NEVENA MARIĆ

variable. The goal is to generate (X1 , X2 , X3 , X4 ) from a multivariate distribution such that they have correlation matrix


1
0.2 −0.1 −0.4
 0.3
1 −0.4 0.3 

Σ=
−0.2 0
1
−0.2
−0.4 0.3 −0.2
1

The first thing to do is to check if this correlation matrix is even possible
with these marginals. Recall that if a random variable X has cdf FX , and
F −1 (u) = inf{a : FX (a) ≥ u}, then for U ∼ Unif([0, 1]), we have that
FX−1 (U ) ∼ X. This is called the inverse transform method of converting
uniforms into variables with the same distribution as X.
Therefore, if a random vector (U1 , U2 , U3 , U4 ) is drawn where the Ui
are marginally uniform, then (FX−11 (U1 ), FX−12 (U2 ), FX−13 (U3 ), FX−14 (U4 )) has the
correct marginals for the Xi .
Note that if U is uniform over [0, 1], then so is 1 − U . The variables U
and 1 − U are said to be antithetic variates. With that in mind, one way to
draw (U1 , U2 , U3 , U4 ) is to set Ui = U or Ui = 1 − U for each i.
This is where a symmetric Bernoulli random vector (B1 , B2 , B3 , B4 ) enters
the picture. For each i, if Bi = 1 then Ui = U and Xi = FX−1i (U ), and if
Bi = 0 then Ui = 1 − U and Xi = FX−1i (1 − U ). In other words,
Xi = Bi FX−1i (U ) + (1 − Bi )FX−1i (1 − U ).
In this way, the problem of generating (X1 , X2 , X3 , X4 ) with the correct
marginals and correlations can be transformed into a problem of generating
(B1 , B2 , B3 , B4 ) with symmetric Bernoulli marginals and perhaps different
correlations. This idea was explored thoroughly in [15].
The Fréchet-Hoeffding bound [10, 14] then says that the correlation between Xi and Xj is maximized when this method is used with Bi = Bj and
minimized with Bi = 1 − Bj . For X2 , FX−12 (u) = −2 ln(1 − u). For X3 ,
FX−13 (u) = 1(u ≤ 0.3) + 4 · 1(u > 0.3). A straightforward calculation then
gives the maximum and minimum possible correlation between X2 and X3
(see e.g. Whitt [26]).
corr(FX−12 (u), FX−13 (u)) = 0.544828 . . .
corr(FX−12 (u), FX−13 (1 − u)) = −0.78818 . . . .
Since the target correlation of −0.4 lies within [−0.78818, 0.544828], it is
possible to achieve.
Because −0.4 falls in the interval, it can be written as a convex linear
combination of −0.78818 and 0.544828:
−0.4 = −0.78818(1 − 0.291209) + .544828(0.291209).
This means that we want a 0.291209 chance that U2 = U3 , and a 1 −
0.291209 chance that U2 6= U3 , or equivalently, a 0.291209 chance that B2 =
B3 and a 1 − 0.291209 chance that B2 6= B3 .

BERNOULLI CORRELATIONS

AND

CUT POLYTOPES

13

Let B2 and B3 be symmetric Bernoulli random variables with a 0.291209
probability of being equal. Then they have a correlation of 2(0.291209)−1 =
−0.417582. Performing this calculation for all the different pairs of random
variables gives a correlation matrix for the symmetric Bernoulli random
variables (B1 , B2 , B3 , B4 ) of


1
0.230940 −0.125988 −0.409330
 0.230940
1
−0.417582 0.332154 
.
ΣB = 
−0.125988 −0.417582
1
−0.263598
−0.409330 0.332154 −0.263598
1
This correlation matrix (i.e. correlation vector ρ = (0.230940, -0.125988,
-0.409330, -0.417582, 0.332154, -0.263598)) corresponds to a linear program
which was solved using lpsolve in R, and the solution was found to be
(α1 , . . . , α8 ), which to four significant figures is

0.04337 0.1284 0.2450 0.1985 0.006894 0.2582 0 0.1193 .

Consider the sixth component 0.2582. Since 6 − 1 = 5 has binary expansion 0101, this means that there is a 0.2582 chance that U2 = U4 = U
and U1 = U3 = 1 − U . Note that both (B1 , B2 , B3 , B4 ) = (0, 1, 0, 1) and
(B1 , B2 , B3 , B4 ) = (1, 0, 1, 0) give rise to U1 = U3 and U2 = U4 . So in
order to make a symmetric Bernoulli distribution, there is equal probability of either of those two vectors occurring. That is, P((B1 , B2 , B3 , B4 ) =
(0, 1, 0, 1)) = P((B1 , B2 , B3 , B4 ) = (1, 0, 1, 0)) = 0.2582/2.
This algorithm for drawing a multivariate distribution for this example
can then be given as follows.
(1) Draw U uniformly over [0, 1]. Draw Y so that P(Y = j) = αj , j =
1, . . . , 8.
(2) Write Y − 1 in binary notation. That is, find bi ∈ {0, 1} so that
P
Y − 1 = 4i=1 2i−1 bi .
(3) For each i ∈ {1, . . . , 4} if bi = 1 then let Ui = U , otherwise let
Ui = 1 − U .
(4) Let X1 = U1 , X2 = −2 ln(1−U2 ), X3 = 1(U3 ≤ 0.3)+4·1(U3 > 0.3),
X4 = Φ−1 (U4 ).
Here Φ is used in the usual fashion as the cdf of a standard normal.
6. Discussion
The set of n×n correlation matrices, the elliptope En×n is a nonpolyhedral
convex set with a nonsmooth boundary. The extreme points of the elliptope
have not been explicitly determined, but there exist characterization results
on the rank one and two extreme points, done by Ycart [27] (see also Li
and Tam [21] and Parthasarathy [23]). Laurent and Poljak [18] proved
that cut matrices (analogous to cut vectors) are actually vertices (that is,
extreme points of rank one) of the elliptope and that En×n can be seen as
a nonpolyhedral relaxation of the cut polytope. In view of theorems proved

14

MARK HUBER AND NEVENA MARIĆ

here it is clear that the vertices of En×n correspond precisely to symmetric
Bernoulli correlations.
Acknowledgments. We are grateful to Gérard Letac for sharing his ideas
with us and for inspiring discussions.
References
[1] D. Avis. Some Polyhedral Cones Related to Metric Spaces. PhD thesis, Stanford University, 1977.
[2] N. R. Chaganty and H. Joe. Range of correlation matrices for dependent Bernoulli
random variables. Biometrika, 93:197–206, 2006.
[3] D. Conway. Multivariate Distributions with specified marginals. PhD thesis, Stanford
University, 1979.
[4] G. Dall’Aglio, S. Kotz, and G. Salinetti. Advances in probability distributions with
given marginals: beyond the copulas. Università degli studi di Roma-La Sapienza,
Rome, Italy, 1991.
[5] B. de Finetti. A proposito di correlazione. Supplemento Statistico ai Nuovi problemi
di Politica Storia ed Economia, (3), 1937.
[6] L. Devroye and G. Letac. Copulas with prescribed correlation matrix. In In Memoriam Marc Yor-Séminaire de Probabilités XLVII, pages 585–601. Springer, 2015.
[7] M. M. Deza and M. Laurent. Geometry of cuts and metrics. Algorithms Combin., 15,
1997.
[8] C. T. dos Santos Dias, A. Samaranayaka, and B. Manly. On the use of correlated beta
random variables with animal population modelling. Ecological modelling, 215(4):293–
300, 2008.
[9] S. Fiorini, S. Massar, S. Pokutta, H. R. Tiwary, and R. de Wolf. Linear vs. semidefinite
extended formulations: Exponential separation and strong lower bounds. volume 44.
ACM, 2012.
[10] M. Fréchet. Sur les tableaux de corrélation dont les marges sont données. Annales de
l’Université de Lyon, 4, 1951.
[11] B. Grünbaum. Convex Polytopes, volume 221 of Graduate Texts in Mathematics.
Springer, 2003.
[12] S. G. Henderson, B. A. Chiera, and R. M. Cooke. Generating “dependent” quasirandom numbers. In Simulation Conference, 2000. Proceedings. Winter, volume 1,
pages 527–536. IEEE, 2000.
[13] R. R. Hill and C. H. Reilly. Composition for multivariate random variables. In J. Tew,
S. Manivannan, D. Sadowski, and A Seila, editors, Proceedings of the 1994 Winter
Simulation Conference, pages 332–339, 1994.
[14] W. Hoeffding. Masstabinvariante korrelatiostheorie. Schriften des Mathematischen
Instituts und des Instituts für Angewandte Mathematik der Universitat Berlin, 5:179–
233, 1940.
[15] M. Huber and N. Marić. Simulation of multivariate distributions with fixed marginals
and correlations. J. Appl. Probab., 52(2):602–608, 2015. arXiv:1311.2002.
[16] S. Kotz and D. D. Mari. Correlation and dependence. World Scientific, 2001.
[17] D. G. Lampard. A stochastic process whose successive intervals between events form
a first-order markov chain. Journal of Applied Probability, 5:648668, 1968.
[18] M. Laurent and S. Poljak. On a positive semidefinite relaxation of the cut polytope.
Linear Algebra and its Applications, 223:439–461, 1995.
[19] A. J. Lawrance and P. A. W. Lewis. A new autoregressive time series model in
exponential variables (near (1)). Advances in Applied Probability, 13(04):826–845,
1981.

BERNOULLI CORRELATIONS

AND

CUT POLYTOPES

15

[20] P. Lévy. Distance de deux variables aléatoires et distance de deux lois de probabilité.
Traité de calcul des probabilités et de ses applications by Emile Borel, I(III):286–292,
1937.
[21] C-K. Li and B-S. Tam. A note on extreme correlation matrices. SIAM Journal on
Matrix Analysis and Applications, 15(3):903–908, 1994.
[22] K. Y. Liang and S. L. Zeger. Longitudinal data analysis using generalized linear
models. Biometrika, 73(1):13–22, 1986.
[23] K. R. Parthasarathy. On extremal correlations. Journal of statistical planning and
inference, 103(1):173–180, 2002.
[24] L. Ruschendorf, B. Schweizer, and M. Taylor. Editor’s Preface: Distributions with
fixed marginals and related topics. IMS Lecture Notes-Monograph Series, 28, 1996.
[25] O. E. Smith and S. I. Adelfang. Gust model based on the bivariate gamma probability
distribution. Journal of Spacecraft and Rockets, 18:545549, 1981.
[26] W. Whitt. Bivariate distributions with given marginals. Ann. Statist., 4(6):1280–1289,
11 1976.
[27] B. Ycart. Extreme points in convex sets of symmetric matrices. Proceedings of the
American Mathematical Society, 95(4):607–612, 1985.
[28] G. M. Ziegler. Lectures on 0/1 polytopes. In Polytopes - combinatorics and computation, pages 1–41. Birkhäuser Basel, 2000.
Claremont McKenna College
E-mail address: mhuber@cmc.edu
University of Missouri-St. Louis
E-mail address: maric@math.umsl.edu

