Improved Boltzmann sampling for the Hadamard product of distributions Andrea Sportiello CNRS and LIPN, Universit´ e Paris Nord, Villetaneuse GASCom 2016, La Marana, Corsica June 4th 2016 Andrea Sportiello Sampling the Hadamard product of two distributions
The problem We have two probability distributions, p ( x ) and q ( x ) Andrea Sportiello Sampling the Hadamard product of two distributions
The problem We have two probability distributions, p ( x ) and q ( x ) p ( x ) q ( x ) Define f ( x ) = � y p ( y ) q ( y ) (the Hadamard product of p and q ) Andrea Sportiello Sampling the Hadamard product of two distributions
The problem We have two probability distributions, p ( x ) and q ( x ) We have two black-box algorithms that sample from p and from q p ( x ) q ( x ) Define f ( x ) = � y p ( y ) q ( y ) (the Hadamard product of p and q ) Andrea Sportiello Sampling the Hadamard product of two distributions
The problem We have two probability distributions, p ( x ) and q ( x ) We have two black-box algorithms that sample from p and from q p ( x ) q ( x ) Define f ( x ) = We want to sample from f � y p ( y ) q ( y ) (the Hadamard product of p and q ) Andrea Sportiello Sampling the Hadamard product of two distributions
The problem We have two probability distributions, p ( x ) and q ( x ) We have two black-box algorithms that sample from p and from q p ( x ) q ( x ) Define f ( x ) = We want to sample from f � y p ( y ) q ( y ) Fast (the Hadamard product of p and q ) Andrea Sportiello Sampling the Hadamard product of two distributions
The obvious rejection algorithm The rules of the game are clear. We have no analytic information whatsoever on p and q . We only have the black boxes The obvious rejection algorithm seems to be the only candidate: Algorithm 1: Obvious rejection begin repeat x p ; y q ; until x = y ; return x end Call ( p , q ) = � x p ( x ) q ( x ) The cycle is repeated on average 1 / ( p , q ) times If we have a size parameter n , this may be large, and we want to make better Andrea Sportiello Sampling the Hadamard product of two distributions
Scaling of the size parameter Assume that ( p , p ), ( p , q ), ( q , q ) are all Θ( n − 2 α ) (by Cauchy–Schwarz, that’s ( p , p ) , ( q , q ) = Θ( n − 2 α ) and p , q not asymptotically orthogonal) x ( p ( x ) 3 + q ( x ) 3 ) ≪ ( p , q ) and that � (that’s normally the case, anyhow we have variants of the algorithm x ( p ( x ) k + q ( x ) k ) ≪ ( p , q ), for any finite k ) that work if � The complexity of the obvious rejection algorithm is Θ( n 2 α ) We want to go down to Θ( n α ln n ) or Θ( n α ) Θ( n 2 α ) ✛ ✲ Andrea Sportiello Sampling the Hadamard product of two distributions
A complexity paradigm In algorithm complexity you normally just count operations. Nonetheless, when you have black boxes, it is wise to count separately operations and black-box queries, which are generally much more expensive that a single operation Here we have two black boxes, for simplicity we will assume they have similar complexity � for black-box queries Nice notation: for operations Then, e.g., Θ( n x + n y ) simplifies to Θ( n x ) if x ≥ y , but stays as is if x < y , as we do not specify the / ratio, except for it being > 1 The complexity of the obvious rejection algorithm is Θ( n 2 α ) + n α ln( n ) We want to go down to Θ( n α ) Andrea Sportiello Sampling the Hadamard product of two distributions
The na¨ ıve ‘birthday paradox’ algorithm Let us consider the following algorithm Algorithm 2: Birthday paradox, first try begin repeat ( x 1 , . . . , x k ) p ; ( y 1 , . . . , y k ) q ; until ∃ ! ( i , j ) | x i = y j ; return x i end Best hope: at each run there are ∼ Poiss k 2 ( p , q ) pairs ( i , j ) such that x i = y j , � so if we tune k ∼ 1 / ( p , q ) the cycle costs Θ( k ), and is repeated on average Θ(1) times The complexity drops down to Θ( n α ) Andrea Sportiello Sampling the Hadamard product of two distributions
An easy win? An easy win? Andrea Sportiello Sampling the Hadamard product of two distributions
An easy win? An easy win? . . . No! Andrea Sportiello Sampling the Hadamard product of two distributions
An easy win? An easy win? . . . No! besides a bunch of solvable minor issues (is the number of good pairs really distributed as a Poissonian?) ıvely searching for a good pair takes time k 2 ) (na¨ there is the one big problem: Andrea Sportiello Sampling the Hadamard product of two distributions
An easy win? An easy win? . . . No! besides a bunch of solvable minor issues (is the number of good pairs really distributed as a Poissonian?) ıvely searching for a good pair takes time k 2 ) (na¨ there is the one big problem: The resulting probability distribution is wrong! Andrea Sportiello Sampling the Hadamard product of two distributions
An easy win? An easy win? . . . No! besides a bunch of solvable minor issues (is the number of good pairs really distributed as a Poissonian?) ıvely searching for a good pair takes time k 2 ) (na¨ there is the one big problem: The resulting probability distribution is wrong! The average number of good pairs ( i , j ) with x i = y j = z is in fact proportional to f ( z ), and boosted by a factor k 2 , which is good. . . . . . but knowing that you have a unique good pair gives a bias! However the whole idea reamins valuable, because, if the average number of good pairs is small, having no further pairs is ‘normal’, and this bias shall also be small. Andrea Sportiello Sampling the Hadamard product of two distributions
Poissonisation of the na¨ ıve algorithm Analysing the bias of the previous algorithm is complicated. It gets easier if we ‘Poissonise’ k , i.e. we rather consider: Algorithm 3: Birthday paradox, Poissonised begin repeat sample k p and k q with Poiss k ; ( x 1 , . . . , x k p ) p ; ( y 1 , . . . , y k q ) q ; until ∃ ! ( i , j ) | x i = y j ; return x i end Call ν z = ( a z , b z ), with a z = # { i | x i = z } and b z = # { j | y j = z } The good fact of Poissonisation is that the ν z ’s are independent random variables so that we can easily analyse what is the probability of returning z Andrea Sportiello Sampling the Hadamard product of two distributions
A formula for the bias Fact: prob ( a z , b z ) = Poiss kp ( z ) ( a z ) Poiss kq ( z ) ( b z ) thus the probability that there is a single pair for z is Poiss kp ( z ) (1) Poiss kq ( z ) (1) = e − k ( p ( z )+ q ( z )) k 2 p ( z ) q ( z ) and the probability that there are no pairs for values z ′ � = z is � �� � 1 − 1 − Poiss kp ( z ′ ) (0) 1 − Poiss kq ( z ′ ) (0) = e kp ( z ′ ) + e kq ( z ′ ) − 1 e − k ( p ( z ′ )+ q ( z ′ )) � � as a result, the probability of returning z is proportional to † f ( z ) e kp ( z ) + e kq ( z ) − 1 and we would be fine if we could devise a final rejection procedure for the spurious factor Bias ( z ) := 1 / ( e kp ( z ) + e kq ( z ) − 1) † I.e., up to a factor, like � z ′ φ ( z ′ ), which is the same for all z , although possibly complicated to calculate Andrea Sportiello Sampling the Hadamard product of two distributions
A philosophical digression So, we want to add a final procedure for correcting the spurious factor Bias ( z ) := 1 / ( e kp ( z ) + e kq ( z ) − 1). In other words, for some C > max z ( e kp ( z ) + e kq ( z ) − 1), z f ( z ) e kp ( z ) + e kq ( z ) − 1 but still with � = Θ(1), C I want to sample a Bernoulli random variable of parameter e kp ( z ) + e kq ( z ) − 1 C However, recall that we have no analytic information on p and q (we only have the black boxes!) Can we ever sample Bern ξ without evaluating ξ ? This is not impossible a priori, just remember the famous algorithm for Bern π/ 4 that makes no use of π . . . Andrea Sportiello Sampling the Hadamard product of two distributions
A first solution with static lists Algorithm 4: Hadamard product, with static lists repeat sample k p and k q with Poiss k ; ( x 1 , . . . , x k ) p ; ( y 1 , . . . , y k ) q ; sort the lists above, produce the list of { ( z , a z , b z ) } ; if N • = 0 and N • = 1 then ν ← a z + b z − 1 for the only z ∈ Ω • ; win ← Bern ν/ 3 ; return z ; 4 until win = true ; 3 1 2 2 3 1 2 1 1 3 3 0 0 1 2 3 4 Andrea Sportiello Sampling the Hadamard product of two distributions
Why this algorithm is unbiased Let us first see why this algorithm is correct 4 Call for short φ ( z ) = e − k ( p ( z )+ q ( z )) 3 1 2 2 The algorithm may return z only if 3 no z ′ � = z arrives out of the Ω • region, 1 2 1 1 3 3 each z ′ makes a factor 0 � 1 + kp ( z ′ ) + ( kp ( z ′ )) 2 + kq ( z ′ ) + ( kq ( z ′ )) 2 � φ ( z ′ ) 0 1 2 3 4 2 2 Then, it shall also happen that z falls within the Ω • region this makes a factor � ( kp ( z )) 2 ( kq ( z )) 2 � kp ( z ) kq ( z ) k 2 p ( z ) q ( z ) φ ( z ) 1 + + + + 2 6 2 6 We would have done if the two factors in parenthesis were proportional They are not. For this reason we put a final Bernoulli rejection, that corrects for the factorials Andrea Sportiello Sampling the Hadamard product of two distributions
Recommend
More recommend
