A discrete version of CMA-ES
Eric Benhamou Jamal Atif Rida Laraki

[
A discrete version of CMA-ES
Eric Benhamou^†, ^, Jamal Atif^, Rida Laraki^
eric.benhamou@aisquareconnect.com or eric.benhamou@dauphine.eu, 
 jamal.atif@dauphine.fr, rida.laraki@dauphine.fr 
A.I. Square Connect^†  Lamsade, Universite Paris Dauphine, PSL^ 
]


Modern machine learning uses more and more advanced optimization techniques to find optimal hyper parameters. Whenever the objective function is non-convex, non continuous and with potentially multiple local minima, standard gradient descent optimization methods fail. A last resource and very different method is to assume that the optimum(s), not necessarily unique, is/are distributed according to a distribution and iteratively to adapt the distribution according to tested points. These strategies originated in the early 1960s, named Evolution Strategy (ES) have culminated with the CMA-ES (Covariance Matrix Adaptation) ES. It relies on a multi variate normal distribution and is supposed to be state of the art for general optimization program. However, it is far from being optimal for discrete variables.
In this paper, we extend the method to multivariate binomial correlated distributions. For such a distribution, we show that it shares similar features to the multi variate normal: independence and correlation is equivalent and correlation is efficiently modeled by interaction between different variables. We discuss this distribution in the framework of the exponential family. We prove that the model can estimate not only pairwise interactions among the two variables but also is capable of modeling higher order interactions. 
This allows creating a version of CMA ES that can accomodate efficiently discrete variables. We provide the corresponding algorithm and conclude.







§ INTRODUCTION

When facing an optimization problem where there is no access to the objective function's gradient or the objective function's gradient is not very smooth, the state of the art techniques rely on stochastic and derivative free algorithm that change radically the point of view of the optimization program. Instead of a deterministic gradient descent, we take a Bayesian point of view and assumes that the optimum is distributed according to a prior statistical distribution and uses particles or random draws to gradually update our statistical distribution. 
Among these method, the covariance matrix adaptation evolution strategy (CMA-ES; e.g., <cit.>) has emerged as the leading stochastic and derivative-free algorithm for solving continuous optimization problems, i.e., for finding the optimum denoted by 𝐱^* of a real-valued objective function f, defined on a subset of a multi dimensional space of dimension d: ℝ^d.This method generates candidate points {𝐱_i}, i ∈{1, 2, …, λ}, from a multivariate Gaussian distribution. It evaluates their objective function (also called fitness) values {f(𝐱_i)}. As the distribution is characterized by its two first moments, it updates the mean vector and covariance matrix by using the sampled points and their fitness values, {(𝐱_i, f(𝐱_i))}. The algorithm keeps repeating the sampling-evaluation-update procedure (which can be seen like an exploration exploitation method until the distribution contracts to a single point or reaches the maximum of iterations. Convergence is measured either by a very small covariance matrix. The different variations around the original method  to improve the convergence investigate various heuristic method to update the distribution parameters. This strongly determines the behavior and efficiency of the whole algorithm. The theoretical foundation of the CMA-ES are that for continuous variables with given first two moments, the maximum entropy distribution is the normal distribution and that the update is based on a maximum-likelihood estimation, making this method based on a statistical principle.

A natural question is to adapt this method for discrete variables. Surprisingly, this has not been done before as the Gaussian distribution is a continuous time distribution inappropriate to discrete variables. One needs to change the underlying distirbution and also find the way to correlate the marginal distribution which can be tricky. However, we show in this paper that multivariate binomials are the natural discrete counterpart of Gaussian distributions. Hence we are able to change CMA Es to accommodate for discrete variables. This is the subject of this paper. In the section <ref>, we introduce the multivariate binomial distribution. Presenting this distribution in the general setting of exponential family, we can easily derive various properties and connect this distribution to maximum entropy. We also proved that for the assumed correlation structure, independence and correlation are equivalent, which is a also a feature of Gaussian distributions. In section <ref>, we present the algorithm. 






§ PRIMER ON MULTIVARIATE BINOMIALS



 §.§ Intuition

To start building some intuition on multivariate binomials, we start by the simplest case, that is a two dimensional Bernoulli. It is the extension to two dimensions of the univariate Bernoulli distribution. A Bernoulli random variable X, is a discrete variable that takes the value 1 with probability p and 0 otherwise. The usual notation for the probability mass function is


    ℙ(X=x) = p^x (1-p)^1-x , x ∈{0,1}


A natural extension is to consider the random vector X = (X_1, X_2). It takes values in the
Cartesian product space {0, 1}^2 = {0, 1}×{0, 1}. If we denote the joint probabilities
p_ij = P(X_1 = i, X_2 = j) for i,j ∈{0, 1}, then the probability for the bivariate Bernoulli writes:

    ℙ(X  =  x)    =   ℙ(X_1=x_1, X_2=x_2) 
         =     p_11^x_1x_2p_10^x_1(1-x_2)p_01^(1-x_1)x_2p_00^(1-x_1)(1-x_2)

with the side conditions that the joint probabilities are between 0 and 1: for i,j ∈{0, 1}, 0 ≤ p_ij≤ 1 
and they sum to one p_00 + p_10 + p_01 + p_11 = 1

It is however better to write the joint distribution in terms of canonical parameters of the related exponential family. 
Hence, if we define

    θ_1     =   log(p_10/p_00), 
    θ_2    =   log(p_01/p_00), 
    θ_12   =   log(p_11p_00/p_10p_01),


and T(x) the vector of sufficient statistics denoted by

    T(X) = (X_1, X_2, X_1 X_2)^T

we can rewrite the distribution as an exponential family distribution as follows:

    ℙ(X = x)    =   exp( ⟨ T(X) , θ⟩ - A(θ) )

where the log partition function A(θ)  is defined such as the probability normalizes to one. It is very easy to check that A(θ)= -log p_00. We can also relate the initial moment parameters {p_ij} for for i,j ∈{0, 1}, to the canonical parameters as follows

 
The moment parameters can be expressed in terms of the canonical parameters as follows;

    p_00 = 1/1 + exp(θ_1) + exp(θ_2) + exp(θ_1 + θ_2 +
    θ_12), 
    
    p_10 = exp(θ_1)/1 + exp(θ_1) + exp(θ_2) + exp(θ_1 + θ_2 +
    θ_12), 
    
    p_01 = exp(θ_2)/1 + exp(θ_1) + exp(θ_2) + exp(θ_1 + θ_2 +
    θ_12), 
    
    p_11 = exp(θ_1 + θ_2 + θ_12)/1 + exp(θ_1) + exp(θ_2) +
    exp(θ_1 + θ_2 + θ_12).



 
See <ref>


The expression of the distribution in terms of the canonical parameters is particularly useful as it indicates immediately that independence and correlation are equivalent 
as in a Gaussian distribution. We will see that this result generalizes to multivariate binomial in the next subsection <ref>.


The components of the bivariate Bernoulli random vector (X_1, X_2) are independent if and only if θ_12 is zero. 
Like for a normal distribution, independence and correlation are equivalent.


 
See <ref>


The equivalence between correlation and independence was already presented in 
<cit.> where it was referred to as Proposition 2.4.1. The
importance of θ_12  referred to as the cross term or u-terms) is discussed and
called cross-product ratio between X_1 and X_2. In <cit.>  and <cit.>, this cross product ratio
is also identified but called the log odds.

Intuitively, there are similarities between Bernoulli (their sum version that is the Binomial) and the Gaussian. And like for the multivariate Gaussian, we can prove that the  marginal and the conditional Bernoulli are still binomial  as shown by the proposition <ref>, making the analogy between Bernoulli (and soon their independent sum version which is the Binomial) and Gaussian even more striking!


In the bivariate Bernoulli vector whose probability mass function is given by <ref> 


  * the marginal distribution of X_1  is also a univariate Bernoulli whose probability mass function is 

    ℙ(X_1 = x_1) = (p_10 + p_11)^x_1(p_00 +p_01)^(1-x_1).


  * the conditional distribution of X_1 given X_2 is also a univariate Bernoulli whose probability mass function is 

    ℙ(X_1  =  x_1 | X_2 = x_2) =  (p_1 x_2/p_1 x_2 + p_0 x_2)^x_1
          (p_0, x_2/p_1, x_2 + p_0, x_2)^1-x_1




 
See <ref>

Before we move on to the generalization, we need to mention a few important facts. 
First, recall that the sum of n independent Bernoulli trials with parameter p is a Binomial with parameter n and p. 
And when we talk about independent sum, it should ring a bell! Independent sum should make you immediately think about the Central Limit Theorem. 
This intuition is absolutely correct and shows the connection between the binomial and Gaussian distribution. 
This is the Moivre Laplace theorem stated below


If the variable B_n follows a binomial distribution with parameters n and p in ]0,1[, then the variable
Z_n = B_n-np/√(np (1-p)) = √(n)B_n / n - p /p (1-p) converges in law to a standard normal law 𝒩 (0,1). 
Another presentation of this result is to say that, for  p ∈ ]0,1[,
as n grows large, for k in the neighborhood of n p we can approximate the binomial distribution by a normal as follows:

    nk p^k (1-p)^n-k≃1/√(2 π n p (1-p))   
    e^-(k - n p)^2/ 2 n p (1-p)



The proof of this theorem is traditionally done with doing a Taylor expansion of the characteristic function. An alternative proof is to use the Sterling formula as well as a Taylor expansion to relate the binomial distribution to the normal one. Historically, de Moivre was the first to establish this theorem in 1733 in the particular case: p =1 / 2. Laplace generalized it in 1812 for any value of p between 0 and 1 and started creating the ground for the central limit theorem that extended this result far beyond. Later on, many more mathematicians generalized and extended this result like Cauchy, Bessel, Poisson but also von Mises, Pólya, Lindeberg, Lévy, Cramér as well as Chebyshev, Markov and Lyapunov.

Second, if we take an infinite sum of Bernoulli, this is a discrete distribution that is the asymptotic limit of the binomial distribution. 
This is also a distribution that is part of the exponential family and is given by the Poisson distribution.


For a large number n of independent Bernoulli trials with probability p such that lim_n→∞ np  = λ, then the corresponding binomial distribution with parameter n and p converges in distribution to the Poisson distribution

 
See <ref>


The two previous results show that binomials, Poisson and Gaussian distributions that are part of the exponential family are closely connected and represent the discrete and continuous version of very similar concepts, namely independent and identically distributed increments. 



 §.§ Maximum entropy

It is also interesting to relate these distributions to maximum Shannon entropy. Let a function : Φ : Ξ→ℝ^d, where Ξ is the space of the random variable X and a vector α∈ℝ^d. It is well known that the maximum entropy distribution whose constraint is given by 𝔼_P [ Φ(X) ] = α is a distribution of the exponential family given by the following theorem


The distribution that maximizes the Shannon entropy : - ∫ p(x) log p(x) dμ(x)
subject to  the constraint 𝔼_P [ ϕ(X) ] = α and the obvious probability constraints ∫ p(x) dμ(x) = 1, p(x) ≥ 0,
is the unique distribution that is part of the exponential family and given by

    p_θ = exp( < θ, Φ(x) > - A(θ) ),

with 

    A(θ ) = log∫exp( < θ, Φ(x) > ) dμ(x)


 
See <ref>



The theorem <ref> works also for discrete distributions. It says that the discrete distribution that maximizes the Shannon entropy  - ∑ p(x) log p(x) subject to the constraint
𝔼_P [ ϕ(X) ] = α and the obvious probability constraints ∑ p(x) = 1, p(x) ≥ 0,
is the unique distribution that is part of the exponential family and given by

    p_θ = exp( < θ, Φ(x) > - A(θ) ),

with 

    A(θ ) = log∫∑ ( < θ, Φ(x) > )



The theorem <ref> implies in particular that the continuous distribution that maximizes entropy with given mean and variance (or equivalently first and second moments) is an exponential family of the form 

    exp( θ_1 x + θ_2 x^2 - A(θ) )

where the log partition function A(θ) is defined to ensure the probability distribution sums to one. This distribution is indeed a normal distribution as it is the exponential of a quadratic form. This is precisely the continuous distribution used in the CMA ES algorithm. Taking a distribution that maximizes the entropy means that we take a distribution that has the less information prior. Or said differently, this is the distribution with the minimal prior structural constraint. If nothing is known, it should therefore be preferred.

Ideally, for our CMA ES discrete adaptation, we would like to find the discrete distribution equivalent of the normal. if we want the discrete distribution with independent increment, we should turn to binomial distributions. Binomials have the other advantage to converge to the normal distribution whenever the discrete parameter converges to a continuous one. Binomials are also distributions that are part of the exponential family.  But we are facing various problems. To keep the discussion simple, let us first look at a single parameter that can take as values all the integer between 0 to n

First of all, we face the issue of controlling the first two moments of our distribution or equivalent to be able to control the mean denoted by μ and the variance denoted by σ^2. Binomial distributions do not have two parameters like normals to be able to adapt to first and second moments constraints as easily as normals.  Indeed for our given parameter n  that is the number of discrete state of our parameter to optimize in the discrete CMA ES, we are only left with a single parameter p for our binomial distribution ℬ(n,p) to accommodate for the constraints. The expectation is given by n p while the variance is given by n p (1-p). If we would like to have a discrete distribution that progressively peaks to the minimum, we would like to be able to force the variance to converge to 0. This will fix the variance to σ^2 = np (1-p). We can easily solve this quadratic equation p^2 - p + σ^2 / n = 0 and use the minimal solution given by 

    p =  1- √( 1 - 4 σ^2 / n)/2

provided that σ^2 ≤ n / 4. As σ will tend to zero, the parameter p will tend to zero. In order to accommodate for the mean constraint, we need a work-around. We see that the discrete parameter is we do not do anything will converge to 0 as p will converge to 0. A solution that is simple is to assume that our discrete parameter is distributed according to 

    (μ + (ℬ(n,p) - n p) )  n

where a  n is a modulo n. It is the remainder of the Euclidean division of a by n. This method will ensure that we sample all possible 0 to n possible value with a mean that is equal to μ and a variance controlled by the parameter p.

Secondly, we would like to use a discrete distribution that maximizes the entropy. This is the case for the continuous version of CMA-ES with the normal distribution. However, for discrete distribution, this maximum entropy condition is not as easy. It is well known that the maximum entropy discrete distribution with a given mean is not the binomial distribution but rather the distribution given by 

    p(X=i) = c   ρ^i

where c = 1 / ∑_i=0^n ρ^i = 1-ρ/1 - ρ^n+1 and where ρ is determined such as ∑_i=0^n c    i   ρ^i = μ which leads to the implicit equation for ρ: 

    (1 + μ) ρ + (μ- (n+1)) ρ^n+1 + (n -μ) ρ^n+2= μ,

using the well known geometric identities: ∑_i=0^n ρ^i =  1-ρ^n+1/1 - ρ and ∑_i=0^n i ρ^i =  ρ1-ρ^n+1/1 - ρ - (n+1) ρ^n+1/1 - ρ___. The distribution is sometimes referred to as the truncated geometric distribution. This is not our desirable binomial distribution. Obviously, we can rule out this truncated geometric distribution as its probability mass function does not make sense for our parameter. The probability mass function is decreasing which is not a desirable feature. Rather, we would like a bell shape, which is the case for our binomial distribution. The tricky question is how to relate our binomial distribution to a maximum entropy principle as this is the case for the normal. 

We can first remark that the binomial distribution is not too far away from a geometric distribution when the number of trials n tends to infinity at least for some terms. Indeed, the probability mass function is given by nk p^k (1-p)^n-k. And using the Sterling formula, we can see that for n large, we can approximate factorial n as follows n ! ∼√(2π n) ( n/e)^n, which leads an asymptotic term similar to the geometric distribution. This gives some hope that there should be a way to relate our binomial distribution to a maximum entropy principle. And the trick here is to reduce the space of possible distributions. It instead of looking at the entire space of distribution, we reduce the space of possible distributions to any Poisson binomial distributions (also referred to in the statistics literature as the generalized binomial distribution), we could find a solution. The latter distribution named after the famous French mathematician Siméon Denis Poisson is the discrete probability distribution of a sum of independent Bernoulli trials that are not necessarily identically distributed. And nicely, restricting the space of possible distribution to any Poisson binomial distributions, theorem <ref> proves that the binomial distribution is the distribution that maximizes the entropy for a given mean.


Among all Poisson binomial distributions with n trials, the distribution that maximizes the Shannon entropy : - ∑ p(x) log p(x)
subject to  the constraint ∑ x p(x) = μ and the obvious probability constraints is the binomial distribution ℬ(n,p)
such that n p = μ

 
See <ref>








 §.§ Multivariate and Correlated Binomials
 
Equipped with the intuition of the first section, we can see the profound connection between multivariate normal and multivariate binomial. 
We will define our multivariate binomial as the sum for n independent trials of multivariate Bernoulli defined as before.

Let X= (X_1, …, X_k) be a k-dimensional random vector of 
possibly correlated binomial random variables that may have different parameters n_i and p_i and
let x = (x_1,… ,x_k) be a realization of X. The joint probability is given naturally by


    ℙ(X_1 = x_1, X_2 = x_2, …, X_K = x_K)  
        = p(0, 0, …,0)^∏_j=1^K (1-x_j)× p(1, 0, …, 0)^x_1 ∏_j=2^K(1-x_j)
         ×  p(0, 1, …, 0)^(1-x_1)x_2 ∏_j=3^K(1-x_j)  ×…×
         × p(1, 1, …, 1) ^∏_j=1^K x_j,



Like for the simple case of section <ref>, we can re-write this joint probability in the exponential form. 
Let us give some notations. 

Let  T(X) be the vector ( X_1, ..., X_k, X_1 X_2, …, X_1 … X_k)^T of size 2^k-1 whose elements represents all the possible 1 to k selection of X_1, …, X_k. 
These  1 to k selections of X_1, …, X_k are all the possible monomial polynomials of X_1, …, X_k of degree 1 to k. By monomial, we mean that we can take only distinct power of  X_1, …, X_k with all of them having an exponent equal to 0 or 1. We also denote by (i_1, …, i_l) an ordered set of 1 ≤ l ≤ k elements of the integers from 1 to k and by Υ_{1,…,k} the set of all the order sets (i_1, …, i_l) with 1 ≤ l ≤ k elements elements.  Υ_{1,…,k} is also the sets of all possible non empty sets with integer elements in {1, …, k}. Similarly,  Υ_{ i_1,…, i_l } is the sets of all possible non empty set with elements in { i_1,…, i_l }. Finally, Υ_{ i_1,…, i_l }^even (respectively  Υ_{ i_1,…, i_l }^odd) is the subset of Υ_{ i_1,…, i_l } for set whose cardinality is even (respectively odd).

We are now able to provide the following proposition that gives the exponential form of the multi variate binomial mass probability function:


The multivariate Bernoulli model has a probability mass function of the exponential form given by

    ℙ(X) = exp( <θ, T(X) > - A(θ))

where the sufficient statistic T(X) is T(X)=(X_1, ..., X_k, X_1 X_2, …, X_1 … X_k), the log partition function A(θ) is A(θ)= -log p(0, 0, …,0) and the coefficients θ are given by:

    θ_i_1,…, i_l   =   log Even / Odd, 
    with  Even    =   ∏_{j_1, .., j_m}∈Υ_{ i_1,.., i_l }^even p( 1  for  all  i_1, …,i_l 
     but   j_1, …, j_m  with 0 
     rest with 0) 
    and   Odd    =   ∏_{j_1, .., j_m}∈Υ_{ i_1,.., i_l }^odd   p( 1  for  all  i_1, …,i_l 
     but  j_1, …, j_m  with 0 
     rest with 0)


Similarly, we can compute the regular probabilities from the canonical parameters as follows:

    p( 1  for  i_1, …, i_l 
     rest with 0) = exp( S^i_1,…, i_l) /D.

where S^i_1,…, i_l is the sum of all the theta parameters indexed by any non empty selection within { i_1,…, i_l }:

    S^i_1,…, i_l =∑_{ i_1, …, i_m}∈Υ_{ i_1,…, i_l }θ_ i_1, …, i_m

with the convention for the empty set, S=1
and D is the normalizing constant such that all the probabilities sum to 1:

    D =∑_ l=0, .. , k 
     1 ≤ i_1  ≤…≤ i_l exp( S^i_1, …, i_l)

with the convention for l=0 that exp( S^i_1, …, i_l)=1.


 
See <ref>



Last but not least, we can extend the result already found for the simple two dimension Bernoulli variable to the general multi dimensional Bernoulli concerning independence and correlation. Recall that one of the important statistical properties for the multivariate Gaussian distribution is the equivalence of independence and no correlation. This is a remarkable properties of the Gaussian (although more could be said about independent and Gaussian as explained for instance in <cit.>).


The independence of a random vector is determined by the separability of coordinates in its probability mass function. If we use the natural (or moment) parameter form of the probability mass function, this is not obvious. However, using the exponential form, the result is almost trivial and is given by the following proposition

The multivariate Bernoulli variable X= (X_1, …, X_k) is independent element-wise if and only if

    θ_i_1, …, i_l = 0     ∀ 1≤ i_1< …< i_l ≤ k, l ≥ 2.



 
See <ref>



The condition of equivalence between independence and no correlation can also be rewritten as

    S^i_1, …, i_l = ∑_k=1^l θ_i_k    ∀ l ≥2.




A general multi variate binomial model implies 2^n-1 parameters, which is way to many when n is large. A simpler model is to impose that only the probabilities involving one state X_i or two states X_i, X_j are non zero. This is in fact the Ising model.








§ ALGORITHM



  §.§.§ CMA-ES estimation

Another radically difference approach is to minimize some cost function depending on the Kalman filter parameters. As opposed to the maximum likelihood approach that tries to find the best suitable distribution that fits the data, this approach can somehow factor in some noise and directly target a cost function that is our final result. Because our model is an approximation of the reality, this noise introduction may leads to a better overall cost function but a worse distribution in terms of fit to the data. 

Let us first introduce the CMA-ES algorithm. Its name stands for covariance matrix adaptation evolution strategy. As it points out, it is an evolution strategy optimization method, meaning that it is a derivative free method that can accomodate non convex optimization problem. The terminology covariance matrix alludes to the fact that the exploration of new points is based on a multinomial distribution whose covariance matrix is progressively determined at each iteration. Hence the covariance matrix adapts in a sense to the sampling space, contracts in dimension that are useless and expands in dimension where natural gradient is steep. This algorithm has led to a large number of papers and articles and we refer to <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.> to cite a few of the numerous articles around CMA-ES. We also refer the reader to the excellent wikipedia page <cit.>.

In order to adapt CMA ES to discrete variables, we change in the algorithm the generation so Gaussian variables into the ones of multi variate binomials as follows:

    x_i ∼ m + ℬ( σ^2 C ) dim

The corresponding algorithm is given in <ref>.





§ CONCLUSION


In this paper, we showed that using multi-variate correlated binomial distribution, we can derive an efficient adaptation of CMA-ES for discrete variable optimization problem using correlated binomials. We have proved that correlated binomials share some similarities with normal distribution in terms of independence and correlation equivalence as well as rich information for correlation structure. In order to avoid too many parameters, we impose that only single state and bi-state probabilities are not null. In the future, we hope to develop additional variations around this CMA-ES version for the combination of discrete and continuous variables mixing potentially multivariate binomial and normal distributions.
plain











§ PROOFS





 §.§ Proof of proposition <ref>

 We can trivially infer all the moment parameters from equations <ref>, <ref> and <ref>.





 §.§ Proof of proposition <ref>

 
The exponential family formulation of the bivariate Bernoulli distribution shows that a necessary and sufficient condition for the distribution to seperable 
into two components with each only depending on x_1 and x_2 respectively is that θ_12=0. This proves the first assertion of proposition <ref>.

Proving equivalence between correlation and independence is the same as proving equivalence between covariance and independence. 
The covariance between X_1 and X_2 is easy to calculate and given by 

    Cov(X_1, X_2)    =    𝔼[ X_1 X_2 ] - 𝔼[ X_1 ] 𝔼[ X_2 ]  
        =    K e^θ_0 + θ_1 + θ_12 - K  (e^θ_0 + e^θ_0 + θ_1 + θ_12) )   K (  e^θ_1 + e^θ_0 + θ_1 + θ_12 ) 
       =     K  e^θ_0 + θ_1 + θ_12 (1 - K e^θ_0 - K e^θ_1 - K e^θ_0 + θ_1 + θ_12) - K^2 e^θ_0 + θ_1
        =     K ^2  e^θ_0 + θ_1 (e^θ_12 - 1 )

where in equation <ref>, we have used that the four probabilities sum to one.
Hence, the correlation or the covariance is null for non trivial probabilities if and only if θ_12 = 0, which is equivalent to the independence.





 §.§ Proof of proposition <ref>

 
For the coordinate X_1, it is trivial to see that 

    ℙ(X_1 = 1)     =    P(X_1 = 1, X_2 = 0) + P(X_1 = 1, X_2 = 1) = p_10 + p_11, 
    ℙ(X_1 = 0)    =     p_00 + p_01, 
    ℙ(X_1 = 1)    +    P(X_1 = 0) = 1.

which shows that X_1 follows the univariate Bernoulli distribution with density given by equation (<ref>).
Likewise, it is trivial to see that

    ℙ(X_1 = 0 | X_2 = 0)    =   P(X_1 = 0, X_2 = 0)/P(X_2 = 0) =  p_00/p_00 + p_10, 
    ℙ(X_1 = 1 | X_2 = 0)    =   p_10/p_00 + p_10, 
    ℙ(X_1 = 1 | X_2 = 0)    +    P(X_1 = 0 | X_2 = 0) = 1.

Similar results apply for the condition X_2 = 1, which shows the second result and concludes the proof.





 §.§ Proof of proposition <ref>

 
Let us write the limit for the binomial distribution when number of trials n→∞, and probability of success in trial p → 0 but np →λ remains finite. We have for a given k


    lim_n→∞nk(λn)^k(1-λn)^n-k
          
    =λ^kk!·lim_n→∞[1(1-1n)(1-2n)…(1-k-1n)]·lim_n→∞(1-λn)^n·1lim_n→∞(1-λn)^k
          
    =lim_n→∞(1-λn)^n·λ^kk!·lim_n→∞[1(1-1n)(1-2n)…(1-k-1n)]·1lim_n→∞(1-λn)^k
          
    =e^-λλ^kk!.

which proves that the binomial converges to the Poisson distribution.





 §.§ Proof of theorem <ref>


We follow the proof of Theorem 11.1.1 of <cit.>. If we write the Lagrangian ℒ(p,θ, θ_0, λ)  for the problem 

    maximize   - ∫ p(x) log p(x) dμ(x) 
       subject to  𝔼_P [ ϕ(X) ] = α  and   ∫ p(x) dμ(x) = 1    and    p(x) ≥ 0

where the Lagrange multipliers (θ, θ_0, λ) are for the three constraints, we have

    ℒ(p,θ, θ_0, λ) = ∫  p(x) log p(x) dμ(x)  + ∑_i=1^d θ_i ( α_i -  ∫  p(x) ϕ_i(x) dμ(x)  ) 
    
    + θ_0 ( ∫ p(x) dμ(x) - 1 ) - ∫λ(x) p(x)    dμ(x)

and noticing that the function to optimize is convex and satisfies the Slater’s constraint, we can use Lagrange duality to characterize the solution as the solution of the critical point given by

    1 +log p(x) - < θ, ϕ(x)> + θ_0 - λ(x) = 0

or equivalently, 

    p(x) = exp(  < θ, ϕ(x)> -1 - θ_0 + λ(x) )

As this solution always satisfies the condition p(x) > 0, we have necessarily that the Lagragian multiplier related to the constraint  p(x) > 0 should be null: λ(x) = 0.  The solution should be a probability distribution, which implies that 

    ∫ p(x)  dμ(x) = 1

which imposes that

    ∫exp(  -1 - θ_0  ) dμ(x) = ∫exp(  < θ, ϕ(x)> ) dμ(x)


or equivalently, writing in the exponential form θ_0 -1 = A(θ) =  log∫exp(  < θ, ϕ(x)> )dμ(x), 
we have that p satisfies

    p_θ(x) = exp(  < θ, ϕ(x) - A(θ) >)

which shows that the distribution is part of the exponential family.  

To prove its uniqueness, we use the fact that the Shannon entropy is related to the Kullback Leibler divergence D_kl as follows:

    H(P) =  - ∫ p(x) log p(x) dμ(x)  = - D_kl( P  P_θ_0) - H(P_θ_0)

which concludes the proof as the  Kullback Leibler divergence D_kl( P  P_θ_0) > 0 unless P= P_θ_0






 §.§ Proof of theorem <ref>


We will prove a stronger result that the entropy H(p_1, …, p_n) of a generalized binomial distribution with parameters n, p_1, …, p_n is Schur concave (see <cit.> for a definition and some properties). A straight consequence of Schur concavity is that the function is maximum for the constant function as follows:

    H(p_1, …, p_n) ≤ H(p̅, …, p̅)

with p̅ = ∑_i=1^n p_i/n. This will prove that the regular binomial distribution satisfies the maximum entropy principle.

Our proof of the Schur concavity uses the same trick as in <cit.>, namely the usage of elementary symmetric functions. 
Let us denote by (X_i)_i=1,…,n the independent Bernoulli variables with parameter p_i and their sum S_n= ∑_i=1^n X_i the variable for the canonical Poisson binomial variable. Its probability mass function writes as

    π_k^n =(S_n=k) = ∑_A∈ F_k∏_i∈ A p_i ∏_j∈ A^c (1-p_j)


where F_k is the set of all subsets of k integers selected from {1,2,3,...,n}. The entropy H(p_1, …, p_n) is permutation symmetric, hence to prove Schur concavity, it suffices to show that the cross term (p_1 - p_2) (∂ H/p_1 - ∂ H/p_2) is negative. Let us compute


    ∂ H/p_1 - ∂ H/p_2 = - ∑_k=0^n (1+ logπ_k^n) (∂π_k^n/p_1 - ∂π_k^n/p_2)


We can notice that for k ≥ 2 and k ≤ n-2

    π_k^n = p_1 p_2 π_k-2^n-2 + (p_1 (1-p_2) + (1-p_1) p_2)  π_k-1^n-2 +  (1-p_1) (1-p_2)   π_k^n-2

where π_j^n-2 = π_j^n-2( p_3, …, p_n). The equation (<ref>) can be extended for k=0, 1 or k=n-1, n with the convention that π_j^n= 0 for j<0 and π_j^n= 0 for j>n. Hence equation (<ref>) is valid for any k. Deriving equation (<ref>) with respect to p_i  leads to


    ∂π_k^n/p_1 - ∂π_k^n/p_2 = - (p_1-p_2) (  π_k-2^n-2 - 2  π_k-1^n-2 +  π_k^n-2 )


Combining equations (<ref>) and (<ref>), we have

    (p_1 - p_2) (∂ H/p_1 - ∂ H/p_2)    =    (p_1 - p_2)^2  ∑_k=0^n (1+ logπ_k^n) (π_k-2^n-2 - 2  π_k-1^n-2 +  π_k^n-2) 
        =     (p_1 - p_2)^2  ∑_k=0^n (logπ_k^n) (π_k-2^n-2 - 2  π_k-1^n-2 +  π_k^n-2) 
        =     (p_1 - p_2)^2  ∑_k=0^n-2π_k^n-2logπ_k+2^n π_k^n / (π_k+1^n)^2


Recall a result that is allegedly attributed to Newton about elementary symmetric functions. Denote the product

    (x+a_1)(x+a_2)… (x+a_n)    =    x^n+ c_1 x^n-1+ c_2 x^n-2+… +c_n

with c_k the kth elementary function of the a's. We have

    c_k-1 c_k+1 < c_k^2

unless all a are equal (see for instance <cit.> theorem 51 page 52 section 2.22). Let us take the function

    ∏_i=1^n  (1-p_i) (x + p_1/1-p_1)(x + p_2/1-p_2)… (x+ p_n/1-p_n)    =   ∏_i=1^n  (1-p_i) (x^n+ c_1 x^n-1+ c_2 x^n-2+… +c_n) 
       =    x^n + π_1^n x^n-1+ π_2^n x^n-2+ … + π_n^n

we have therefore that 

    π_k^n = ∏_i=1^n  (1-p_i) c_k

Combining equations (<ref>) and (<ref>), we proved that

    π_k+2^n π_k^n < (π_k+1^n)^2

which concludes the proof

	




 §.§ Proof of proposition<ref>


Comparing the equations (<ref>) and (<ref>), and using the provided sufficient statistics given in  proposition<ref>, we see by identification that the parameters θ should be given by the equation (<ref>) with the terms with a plus sign given by equation (<ref>) and the terms with a negative sign given by equation (<ref>)

Similarly, if we take equation (<ref>),  (<ref>) and (<ref>), we can notice that the probabilities given sum to one, that D= 1 / p_0, …, 0 and that from the expression giving the theta's (<ref>), we back out the probabilities. This concludes the proof.




 §.§ Proof of proposition<ref>


The proof of proposition <ref> is immediate using the exponential form as the independence is equivalent to the separability which is equivalent to equation (<ref>)






§ PSEUDO CODE



