随机算法 (Spring 2014)/Introduction and Probability Space and 随机算法 (Spring 2014)/Random Variables: Difference between pages

From TCS Wiki
(Difference between pages)
Jump to navigation Jump to search
imported>Etone
 
imported>Etone
 
Line 1: Line 1:
=Introduction=
=Random Variable=
This course will study ''Randomized Algorithms'', the algorithms that use randomness in computation.
{{Theorem|Definition (random variable)|
;Why do we use randomness in computation?
:A random variable <math>X</math> on a sample space <math>\Omega</math> is a real-valued function <math>X:\Omega\rightarrow\mathbb{R}</math>. A random variable X is called a '''discrete''' random variable if its range is finite or countably infinite.
* Randomized algorithms can be simpler than deterministic ones.
}}
:(median selection, load balancing, etc.)
* Randomized algorithms can be faster than the best known deterministic algorithms.
:(min-cut, checking matrix multiplication, primality testing, etc.)
* Randomized algorithms can do things that deterministic algorithms cannot do.
:(routing, volume estimation, communication complexity, data streams, etc.)
* Randomized algorithms may lead us to smart deterministic algorithms.
:(hashing, derandomization, SL=L, Lovász Local Lemma, etc.)
* Randomness is presented in the input.
:(average-case analysis, smoothed analysis, learning, etc.)
* Some deterministic problems are random in nature.
:(counting, inference, etc.)
* ...


;How is randomness used in computation?
For a random variable <math>X</math> and a real value <math>x\in\mathbb{R}</math>, we write "<math>X=x</math>" for the event <math>\{a\in\Omega\mid X(a)=x\}</math>, and denote the probability of the event by
* To hit a witness/certificate.
:<math>\Pr[X=x]=\Pr(\{a\in\Omega\mid X(a)=x\})</math>.
:(identity testing, fingerprinting, primality testing, etc.)
* To avoid worst case or to deal with adversaries.
:(randomized quick sort, perfect hashing, etc.)
* To simulate random samples.
:(random walk, Markov chain Monte Carlo, approximate counting etc.)
* To enumerate/construct solutions.
:(the probabilistic method, min-cut, etc.)
* To break symmetry.
:(mutual exclusion, leader election, parrallel/distributed algorithms, etc.)
* ...


== Principles in probability theory  ==
The independence can also be defined for variables:
The course is organized by the advancedness of the probabilistic tools. We do this for two reasons: First, for randomized algorithms, analysis is usually more difficult and involved than the algorithm itself; and second, getting familiar with these probability principles will help you understand the true reasons for which the smart algorithms are designed.
{{Theorem
* '''Basics of probability theory''': probability space, events, the union bound, independence, conditional probability.
|Definition (Independent variables)|
* '''Moments and deviations''': random variables, expectation, linearity of expectation, Markov's inequality, variance, second moment method.
:Two random variables <math>X</math> and <math>Y</math> are '''independent''' if and only if
* '''The probabilistic method''': averaging principle, threshold phenomena, Lovász Local Lemma.
::<math>
* '''Concentrations''': Chernoff-Hoeffding bound, martingales, Azuma's inequality, bounded difference method.
\Pr[(X=x)\wedge(Y=y)]=\Pr[X=x]\cdot\Pr[Y=y]
* '''Markov chains and random walks''': Markov chians, random walks, hitting/cover time, mixing time, coupling, conductance.
</math>
:for all values <math>x</math> and <math>y</math>. Random variables <math>X_1, X_2, \ldots, X_n</math> are '''mutually independent''' if and only if, for any subset <math>I\subseteq\{1,2,\ldots,n\}</math> and any values <math>x_i</math>, where <math>i\in I</math>,
::<math>\begin{align}
\Pr\left[\bigwedge_{i\in I}(X_i=x_i)\right]
&=
\prod_{i\in I}\Pr[X_i=x_i].
\end{align}</math>
}}


=Probability Space=
Note that in probability theory, the "mutual independence" is <font color="red">not</font> equivalent with "pair-wise independence", which we will learn in the future.
The axiom foundation of probability theory is laid by [http://en.wikipedia.org/wiki/Andrey_Kolmogorov Kolmogorov], one of the greatest mathematician of the 20th century, who advanced various very different fields of mathematics.


{{Theorem|Definition (Probability Space)|
== Expectation ==
A '''probability space''' is a triple <math>(\Omega,\Sigma,\Pr)</math>.
Let <math>X</math> be a discrete '''random variable'''. The expectation of <math>X</math> is defined as follows.
*<math>\Omega</math> is a set, called the '''sample space'''.  
{{Theorem
*<math>\Sigma\subseteq 2^{\Omega}</math> is the set of all '''events''', satisfying:
|Definition (Expectation)|
*:(K1). <math>\Omega\in\Sigma</math> and <math>\empty\in\Sigma</math>. (The ''certain'' event and the ''impossible'' event.)
:The '''expectation''' of a discrete random variable <math>X</math>, denoted by <math>\mathbf{E}[X]</math>, is given by
*:(K2). If <math>A,B\in\Sigma</math>, then <math>A\cap B, A\cup B, A-B\in\Sigma</math>. (Intersection, union, and diference of two events are events).
::<math>\begin{align}
* A '''probability measure''' <math>\Pr:\Sigma\rightarrow\mathbb{R}</math> is a function that maps each event to a nonnegative real number, satisfying
\mathbf{E}[X] &= \sum_{x}x\Pr[X=x],
*:(K3). <math>\Pr(\Omega)=1</math>.
\end{align}</math>
*:(K4). If <math>A\cap B=\emptyset</math> (such events are call ''disjoint'' events), then <math>\Pr(A\cup B)=\Pr(A)+\Pr(B)</math>.
:where the summation is over all values <math>x</math> in the range of <math>X</math>.
*:(K5*). For a decreasing sequence of events <math>A_1\supset A_2\supset \cdots\supset A_n\supset\cdots</math> of events with <math>\bigcap_n A_n=\emptyset</math>, it holds that <math>\lim_{n\rightarrow \infty}\Pr(A_n)=0</math>.
}}
}}


;Remark
== Linearity of Expectation ==
* In general, the set <math>\Omega</math> may be continuous, but we only consider '''discrete''' probability in this lecture, thus we assume that <math>\Omega</math> is either finite or countably infinite.
Perhaps the most useful property of expectation is its '''linearity'''.
* Sometimes it is convenient to assume <math>\Sigma=2^{\Omega}</math>, i.e. the events enumerates all subsets of <math>\Omega</math>. But in general, a probability space is well-defined by any <math>\Sigma</math> satisfying (K1) and (K2). Such <math>\Sigma</math> is called a <math>\sigma</math>-algebra defined on <math>\Omega</math>.
* The last axiom (K5*) is redundant if <math>\Sigma</math> is finite, thus it is only essential when there are infinitely many events. The role of axiom (K5*) in probability theory is like [http://en.wikipedia.org/wiki/Zorn's_lemma Zorn's Lemma] (or equivalently the [http://en.wikipedia.org/wiki/Axiom_of_choice Axiom of Choice]) in axiomatic set theory.


Useful laws for probability can be deduced from the ''axioms'' (K1)-(K5).
{{Theorem
{{Theorem|Proposition|
|Theorem (Linearity of Expectations)|
# Let <math>\bar{A}=\Omega\setminus A</math>. It holds that <math>\Pr(\bar{A})=1-\Pr(A)</math>.
:For any discrete random variables <math>X_1, X_2, \ldots, X_n</math>, and any real constants <math>a_1, a_2, \ldots, a_n</math>,
# If <math>A\subseteq B</math> then <math>\Pr(A)\le\Pr(B)</math>.
::<math>\begin{align}
\mathbf{E}\left[\sum_{i=1}^n a_iX_i\right] &= \sum_{i=1}^n a_i\cdot\mathbf{E}[X_i].
\end{align}</math>
}}
}}
{{Proof|
{{Proof| By the definition of the expectations, it is easy to verify that (try to prove by yourself):
# The events <math>\bar{A}</math> and <math>A</math> are disjoint and <math>\bar{A}\cup A=\Omega</math>. Due to Axiom (K4) and (K3), <math>\Pr(\bar{A})+\Pr(A)=\Pr(\Omega)=1</math>.
for any discrete random variables <math>X</math> and <math>Y</math>, and any real constant <math>c</math>,
# The events <math>A</math> and <math>B\setminus A</math> are disjoint and <math>A\cup(B\setminus A)=B</math> since <math>A\subseteq B</math>. Due to Axiom (K4), <math>\Pr(A)+\Pr(B\setminus A)=\Pr(B)</math>, thus <math>\Pr(A)\le\Pr(B)</math>.
* <math>\mathbf{E}[X+Y]=\mathbf{E}[X]+\mathbf{E}[Y]</math>;
* <math>\mathbf{E}[cX]=c\mathbf{E}[X]</math>.
The theorem follows by induction.
}}
}}
The linearity of expectation gives an easy way to compute the expectation of a random variable if the variable can be written as a sum.
;Example
: Supposed that we have a biased coin that the probability of HEADs is <math>p</math>. Flipping the coin for n times, what is the expectation of number of HEADs?
: It looks straightforward that it must be np, but how can we prove it? Surely we can apply the definition of expectation to compute the expectation with brute force. A more convenient way is by the linearity of expectations: Let <math>X_i</math> indicate whether the <math>i</math>-th flip is HEADs. Then <math>\mathbf{E}[X_i]=1\cdot p+0\cdot(1-p)=p</math>, and the total number of HEADs after n flips is <math>X=\sum_{i=1}^{n}X_i</math>. Applying the linearity of expectation, the expected number of HEADs is:
::<math>\mathbf{E}[X]=\mathbf{E}\left[\sum_{i=1}^{n}X_i\right]=\sum_{i=1}^{n}\mathbf{E}[X_i]=np</math>.


;Notation
The real power of the linearity of expectations is that it does not require the random variables to be independent, thus can be applied to any set of random variables. For example:
An event <math>A\subseteq\Omega</math> can be represented as <math>A=\{a\in\Omega\mid \mathcal{E}(a)\}</math> with a predicate <math>\mathcal{E}</math>.
:<math>\mathbf{E}\left[\alpha X+\beta X^2+\gamma X^3\right] = \alpha\cdot\mathbf{E}[X]+\beta\cdot\mathbf{E}\left[X^2\right]+\gamma\cdot\mathbf{E}\left[X^3\right].</math>


The predicate notation of probability is
However, do not exaggerate this power!
:<math>\Pr[\mathcal{E}]=\Pr(\{a\in\Omega\mid \mathcal{E}(a)\})</math>.
* For an arbitrary function <math>f</math> (not necessarily linear), the equation <math>\mathbf{E}[f(X)]=f(\mathbf{E}[X])</math> does <font color="red">not</font> hold generally.
* For variances, the equation <math>var(X+Y)=var(X)+var(Y)</math> does <font color="red">not</font> hold without further assumption of the independence of <math>X</math> and <math>Y</math>.


During the lecture, we mostly use the predicate notation instead of subset notation.
==Conditional Expectation ==


== Independence ==
Conditional expectation can be accordingly defined:
{{Theorem
{{Theorem
|Definition (Independent events)|
|Definition (conditional expectation)|
:Two events <math>\mathcal{E}_1</math> and <math>\mathcal{E}_2</math> are '''independent''' if and only if
:For random variables <math>X</math> and <math>Y</math>,
::<math>\begin{align}
::<math>
\Pr\left[\mathcal{E}_1 \wedge \mathcal{E}_2\right]
\mathbf{E}[X\mid Y=y]=\sum_{x}x\Pr[X=x\mid Y=y],
&=
</math>
\Pr[\mathcal{E}_1]\cdot\Pr[\mathcal{E}_2].
:where the summation is taken over the range of <math>X</math>.
\end{align}</math>
}}
}}
This definition can be generalized to any number of events:
 
There is also a '''law of total expectation'''.
{{Theorem
{{Theorem
|Definition (Independent events)|
|Theorem (law of total expectation)|
:Events <math>\mathcal{E}_1, \mathcal{E}_2, \ldots, \mathcal{E}_n</math> are '''mutually independent''' if and only if, for any subset <math>I\subseteq\{1,2,\ldots,n\}</math>,
:Let <math>X</math> and <math>Y</math> be two random variables. Then
::<math>\begin{align}
::<math>
\Pr\left[\bigwedge_{i\in I}\mathcal{E}_i\right]
\mathbf{E}[X]=\sum_{y}\mathbf{E}[X\mid Y=y]\cdot\Pr[Y=y].
&=
</math>
\prod_{i\in I}\Pr[\mathcal{E}_i].
\end{align}</math>
}}
}}


Note that in probability theory, the "mutual independence" is <font color="red">not</font> equivalent with "pair-wise independence", which we will learn in the future.
=Random Quicksort=
Given as input a set <math>S</math> of <math>n</math> numbers, we want to sort the numbers in <math>S</math> in increasing order. One of the most famous algorithm for this problem is the  [http://en.wikipedia.org/wiki/Quicksort Quicksort] algorithm.
* if <math>|S|>1</math> do:
** pick an <math>x\in S</math> as the ''pivot'';
** partition <math>S</math> into <math>S_1</math>, <math>\{x\}</math>, and <math>S_2</math>, where all numbers in <math>S_1</math> are smaller than <math>x</math> and all numbers in <math>S_2</math> are  larger than <math>x</math>;
** recursively sort <math>S_1</math> and <math>S_2</math>;


= Model of Computation =
The time complexity of this sorting algorithm is measured by the '''number of comparisons'''.
Our model of computation extends the standard model ([http://en.wikipedia.org/wiki/Turing_machine Turing machine] or [http://en.wikipedia.org/wiki/Random-access_machine random-access machine]) with access to uniform and independent random bits (fair coin flips). On a fixed input, the behavior of the algorithm is random. To be specific, the output or running time of the algorithm may be random.


== Monte Carlo algorithms ==
For the '''deterministic''' quicksort algorithm, the pivot is picked from a fixed position (e.g. the first number in the array). The worst-case time complexity in terms of number of comparisons is <math>\Theta(n^2)</math>.
[http://en.wikipedia.org/wiki/Monte_carlo_algorithm Monte Carlo algorithms] always returns in finite steps but may output the wrong answer. For decision problems (problems with two answers "yes" and "no"), the Monte Carlo algorithms are further divided into those with ''one-sided errors'' and ''two-sided errors''.


;Monte Carlo algorithms with one-sided errors
We consider the following randomized version of the quicksort.
:These algorithms only make errors in one direction, which may be further divided into two cases:
* if <math>|S|>1</math> do:
:'''False positive''': If the true answer is "yes" then the algorithm returns "yes" with probability 1, and if the true answer is "no" then the algorithm returns "no" with probability at least <math>\epsilon</math>, where <math>0<\epsilon< 1</math> is the ''confidence''. The algorithm may return a wrong "yes" while the true answer is "no".
** ''uniformly'' pick a ''random'' <math>x\in S</math> as the pivot;
:The one-sided error can be reduced by independent repetitions. Run the algorithm independently for <math>t</math> times, output "yes" if all running instances return "yes", and output "no" if otherwise. If the true answer is "yes" this new algorithm returns "yes" since all running instances are guaranteed to do so, and if the true answer is "no" the new algorithm returns "yes" only if all running instances return "yes", whose probability is bounded by <math>(1-\epsilon)^t</math>, which can be reduced to any <math>\delta\in(0,1)</math> by setting <math>t=O\left(\frac{1}{\epsilon}\log\frac{1}{\delta}\right)</math>.
** partition <math>S</math> into <math>S_1</math>, <math>\{x\}</math>, and <math>S_2</math>, where all numbers in <math>S_1</math> are smaller than <math>x</math> and all numbers in <math>S_2</math> are larger than <math>x</math>;
** recursively sort <math>S_1</math> and <math>S_2</math>;


:'''False negative''': If the true answer is "yes" then the algorithm returns "yes" with probability at least <math>\epsilon</math>, and if the true answer is "no" then the algorithm returns "no" with probability 1. The algorithm may return a wrong "no" while the true answer is "yes". The error can be reduced in the same way.
== Analysis of Random Quicksort==
Our goal is to analyze the expected number of comparisons during an execution of RandQSort with an arbitrary input <math>S</math>. We achieve this by measuring the chance that each pair of elements are compared, and summing all of them up due to [http://en.wikipedia.org/wiki/Expected_value#Linearity Linearity of Expectation].


Let <math>a_i</math> denote the <math>i</math>th smallest element in <math>S</math>.
Let <math>X_{ij}\in\{0,1\}</math> be the random variable which indicates whether <math>a_i</math> and <math>a_j</math> are compared during the execution of RandQSort. That is:


;Monte Carlo algorithms with two-sided errors
:<math>
:These algorithms make errors in both directions. If the true answer is "yes" then the algorithm returns "yes" with probability at least <math>\frac{1}{2}+\epsilon</math>, and if the true answer is "no" then the algorithm returns "no" with probability at least <math>\frac{1}{2}+\epsilon</math>, where <math>\epsilon\in \left(0,\frac{1}{2}\right)</math> is a bias.
:The error can be reduced by repetitions and majority vote. Run the algorithm independently for <math>t</math> times, output "yes" if over half running instances return "yes", and output "no" if otherwise. The numbers of "yes"s and "no"s in the <math>t</math> trials follow the Binomial distribution. For each <math>0\le i\le t</math>, the probability that there are precisely <math>i</math> correct answers in <math>t</math> trials is given by
::<math>
{t\choose i}\left(\frac{1}{2}+\epsilon\right)^i\left(\frac{1}{2}-\epsilon\right)^{t-i},
</math>
:and the new algorithm returns a wrong answer only if at most <math>\lfloor t/2\rfloor</math> correct answers in <math>t</math> trials, which is given by
::<math>
\begin{align}
\begin{align}
\sum_{i=0}^{\lfloor t/2\rfloor} {t\choose i}\left(\frac{1}{2}+\epsilon\right)^i\left(\frac{1}{2}-\epsilon\right)^{t-i}
X_{ij} &=
&\le
\begin{cases}
\sum_{i=0}^{\lfloor t/2\rfloor} {t\choose i}\left(\frac{1}{2}+\epsilon\right)^{t/2}\left(\frac{1}{2}-\epsilon\right)^{t/2}\\
1 & a_i\mbox{ and }a_j\mbox{ are compared}\\
&=
0 & \mbox{otherwise}
\left(\frac{1}{4}-\epsilon^2\right)^{t/2}\sum_{i=0}^{\lfloor t/2\rfloor}{t\choose i}\\
\end{cases}.
&\le
\left(\frac{1}{4}-\epsilon^2\right)^{t/2}2^t\\
&=(1-4\epsilon^2)^{t/2},
\end{align}
\end{align}
</math>
</math>
:which can be reduced to any <math>\delta\in(0,1)</math> by setting <math>t=O\left(\frac{1}{\epsilon^2}\log\frac{1}{\delta}\right)</math>.


== Las Vegas algorithms ==
Elements <math>a_i</math> and <math>a_j</math> are compared only if one of them is chosen as pivot. After comparison they are separated (thus are never compared again). So we have the following observations:
[http://en.wikipedia.org/wiki/Las_Vegas_algorithm Las Vegas algorithms] always output correct answers but the running time is random. The time complexity of a Las Vegas algorithm is measure by the expected running time. The concept of Las Vegas algorithm is introduced by Babai in 1979 in his seminal work on graph isomorphsm testing.
 
'''Observation 1:  Every pair of <math>a_i</math> and <math>a_j</math> are compared at most once.'''
 
Therefore the sum of <math>X_{ij}</math> for all pair <math>\{i, j\}</math> gives the total number of comparisons. The expected number of comparisons is <math>\mathbf{E}\left[\sum_{i=1}^n\sum_{j>i}X_{ij}\right]</math>. Due to [http://en.wikipedia.org/wiki/Expected_value#Linearity Linearity of Expectation], <math>\mathbf{E}\left[\sum_{i=1}^n\sum_{j>i}X_{ij}\right] = \sum_{i=1}^n\sum_{j>i}\mathbf{E}\left[X_{ij}\right]</math>.  
Our next step is to analyze <math>\mathbf{E}\left[X_{ij}\right]</math> for each <math>\{i, j\}</math>.


A Las Vegas algorithm can be converted to a Monte Carlo algorithm by truncating. The error of the resulting Monte Carlo algorithm can be made one-sided and can be bounded by Markov's inequality. There is no general way known to convert a Monte Carlo algorithm to a Las Vegas algorithm.
By the definition of expectation and <math>X_{ij}</math>,


=  Checking Matrix Multiplication=
:<math>\begin{align}
[[File: matrix_multiplication.png|thumb|360px|right|The evolution of time complexity <math>O(n^{\omega})</math> for matrix multiplication.]]
\mathbf{E}\left[X_{ij}\right]
&= 1\cdot \Pr[a_i\mbox{ and }a_j\mbox{ are compared}] + 0\cdot \Pr[a_i\mbox{ and }a_j\mbox{ are not compared}]\\
&= \Pr[a_i\mbox{ and }a_j\mbox{ are compared}].
\end{align}</math>


Let <math>\mathbb{F}</math> be a feild (you may think of it as the filed <math>\mathbb{Q}</math> of rational numbers, or the finite field <math>\mathbb{Z}_p</math> of integers modulo prime <math>p</math>). We suppose that each field operation (addition, subtraction, multiplication, division) has unit cost. This model is called the '''unit-cost RAM''' model, which is an ideal abstraction of a computer.
We are going to bound this probability.


Consider the following problem:
'''Observation 2: <math>a_i</math> and <math>a_j</math> are compared if and only if one of them is chosen as pivot when they are still in the same subset.'''
* '''Input''': Three <math>n\times n</math> matrices <math>A</math>, <math>B</math>, and <math>C</math> over the field <math>\mathbb{F}</math>.
* '''Output''': "yes" if <math>C=AB</math> and "no" if otherwise.


A naive way to solve this is to multiply <math>A</math> and <math>B</math> and compare the result with <math>C</math>.
This is easy to verify: just check the algorithm. The next one is a bit complicated.
The straightforward algorithm for matrix multiplication takes <math>O(n^3)</math> time, assuming that each arithmetic operation takes unit time.
The [http://en.wikipedia.org/wiki/Strassen_algorithm Strassen's algorithm] discovered in 1969 now implemented by many numerical libraries runs in time <math>O(n^{\log_2 7})\approx O(n^{2.81})</math>. Strassen's algorithm starts the search for fast matrix multiplication algorithms. The [http://en.wikipedia.org/wiki/Coppersmith%E2%80%93Winograd_algorithm Coppersmith–Winograd algorithm] discovered in 1987 runs in time <math>O(n^{2.376})</math> but is only faster than Strassens' algorithm on extremely large matrices due to the very large constant coefficient. This has been the best known for decades, until recently Stothers got an <math>O(n^{2.3737})</math> algorithm in his PhD thesis in 2010, and independently Vassilevska Williams got an <math>O(n^{2.3727})</math> algorithm in 2012. Both these improvements are based on generalization of Coppersmith–Winograd algorithm. It is unknown whether the matrix multiplication can be done in time <math>O(n^{2+o(1)})</math>.


== Freivalds Algorithm ==
'''Observation 3: If <math>a_i</math> and <math>a_j</math> are still in the same subset then all <math>\{a_i, a_{i+1}, \ldots, a_{j-1}, a_{j}\}</math> are in the same subset.'''
The following is a very simple randomized algorithm due to Freivalds, running in <math>O(n^2)</math> time:


{{Theorem|Algorithm (Freivalds, 1979)|
We can verify this by induction. Initially, <math>S</math> itself has the property described above; and partitioning any <math>S</math> with the property into <math>S_1</math> and <math>S_2</math> will preserve the property for both <math>S_1</math> and <math>S_2</math>. Therefore Claim 3 holds.
*pick a vector <math>r \in\{0, 1\}^n</math> uniformly at random;
 
*if <math>A(Br) = Cr</math> then return "yes" else return "no";
Combining Observation 2 and 3, we have:
}}
The product <math>A(Br)</math> is computed by first multiplying <math>Br</math> and then <math>A(Br)</math>.
The running time of Freivalds algorithm is <math>O(n^2)</math> because the algorithm computes 3 matrix-vector multiplications.  


If <math>AB=C</math> then <math>A(Br) = Cr</math> for any <math>r \in\{0, 1\}^n</math>, thus the algorithm will return a "yes" for any positive instance (<math>AB=C</math>).
'''Observation 4: <math>a_i</math> and <math>a_j</math> are compared only if one of <math>\{a_i, a_j\}</math> is chosen from <math>\{a_i, a_{i+1}, \ldots, a_{j-1}, a_{j}\}</math>.'''
But if <math>AB \neq C</math> then the algorithm will make a mistake if it chooses such an <math>r</math> that <math>ABr = Cr</math>. However, the following lemma states that the probability of this event is bounded.


{{Theorem|Lemma|
And,
:If <math>AB\neq C</math> then for a uniformly random <math>r \in\{0, 1\}^n</math>,
::<math>\Pr[ABr = Cr]\le \frac{1}{2}</math>.
}}
{{Proof| Let <math>D=AB-C</math>. The event <math>ABr=Cr</math> is equivalent to that <math>Dr=0</math>. It is then sufficient to show that for a <math>D\neq \boldsymbol{0}</math>, it holds that <math>\Pr[Dr = \boldsymbol{0}]\le \frac{1}{2}</math>.


Since <math>D\neq \boldsymbol{0}</math>, it must have at least one non-zero entry. Suppose that <math>D_{ij}\neq 0</math>.
'''Observation 5: Every one of <math>\{a_i, a_{i+1}, \ldots, a_{j-1}, a_{j}\}</math> is chosen equal-probably.'''


We assume the event that <math>Dr=\boldsymbol{0}</math>. In particular, the <math>i</math>-th entry of <math>Dr</math> is
This is because the Random Quicksort chooses the pivot ''uniformly at random''.
:<math>(Dr)_{i}=\sum_{k=1}^n D_{ik}r_k=0.</math>
The <math>r_j</math> can be calculated by
:<math>r_j=-\frac{1}{D_{ij}}\sum_{k\neq j}^n D_{ik}r_k.</math>
Once all other entries <math>r_k</math> with <math>k\neq j</math> are fixed, there is a unique solution of <math>r_j</math>. Therefore, the number of <math>r\in\{0,1\}^n</math> satisfying <math>Dr=\boldsymbol{0}</math> is at most <math>2^{n-1}</math>. The probability that <math>ABr=Cr</math> is bounded as
:<math>\Pr[ABr=Cr]=\Pr[Dr=\boldsymbol{0}]\le\frac{2^{n-1}}{2^n}=\frac{1}{2}</math>.
}}


When <math>AB=C</math>, Freivalds algorithm always returns "yes"; and when <math>AB\neq C</math>, Freivalds algorithm returns "no" with probability at least 1/2.
Observation 4 and 5 together imply:


To improve its accuracy, we can run Freivalds algorithm for <math>k</math> times, each time with an ''independent'' <math>r\in\{0,1\}^n</math>, and return "yes" if and only if all running instances returns "yes".
:<math>\begin{align}
\Pr[a_i\mbox{ and }a_j\mbox{ are compared}]
&\le \frac{2}{j-i+1}.
\end{align}</math>


{{Theorem|Freivalds' Algorithm (multi-round)|
{|border="1"
*pick <math>k</math> vectors <math>r_1,r_2,\ldots,r_k \in\{0, 1\}^n</math> uniformly and independently at random;
|'''Remark:''' Perhaps you feel confused about the above argument. You may ask: "''The algorithm chooses pivots for many times during the execution. Why in the above argument, it looks like the pivot is chosen only once?''" Good question! Let's see what really happens by looking closely.
*if <math>A(Br_i) = Cr_i</math> for all <math>i=1,\ldots,k</math> then return "yes" else return "no";
}}


If <math>AB=C</math>, then the algorithm returns a "yes" with probability 1. If <math>AB\neq C</math>, then due to the independence, the probability that all <math>r_i</math> have <math>ABr_i=C_i</math> is at most <math>2^{-k}</math>, so the algorithm returns "no" with probability at least <math>1-2^{-k}</math>. For any <math>0<\epsilon<1</math>, choose <math>k=\log_2 \frac{1}{\epsilon}</math>. The algorithm runs in time <math>O(n^2\log_2\frac{1}{\epsilon})</math> and has a one-sided error (false positive) bounded by <math>\epsilon</math>.
For any pair <math>a_i</math> and <math>a_j</math>, initially <math>\{a_i, a_{i+1}, \ldots, a_{j-1}, a_{j}\}</math> are all in the same set <math>S</math> (obviously!). During the execution of the algorithm, the set which containing <math>\{a_i, a_{i+1}, \ldots, a_{j-1}, a_{j}\}</math> are shrinking (due to the pivoting), until one of <math>\{a_i, a_{i+1}, \ldots, a_{j-1}, a_{j}\}</math> is chosen, and the set is partitioned into different subsets. We ask for the probability that the chosen one is among <math>\{a_i, a_j\}</math>. So we really care about "the last" pivoting before <math>\{a_i, a_{i+1}, \ldots, a_{j-1}, a_{j}\}</math> is split.


=Polynomial Identity Testing (PIT) =
Formally, let <math>Y</math> be the random variable denoting the pivot element. We know that for each <math>a_k\in\{a_i, a_{i+1}, \ldots, a_{j-1}, a_{j}\}</math>, <math>Y=a_k</math> with the same probability, and <math>Y\not\in\{a_i, a_{i+1}, \ldots, a_{j-1}, a_{j}\}</math> with an unknown probability (remember that there might be other elements in the same subset with <math>\{a_i, a_{i+1}, \ldots, a_{j-1}, a_{j}\}</math>). The probability we are looking for is actually
Consider the following problem of '''Polynomial Identity Testing (PIT)''':
<math>\Pr[Y\in \{a_i, a_j\}\mid Y\in\{a_i, a_{i+1}, \ldots, a_{j-1}, a_{j}\}]</math>, which is always <math>\frac{2}{j-i+1}</math>, provided that <math>Y</math> is uniform over <math>\{a_i, a_{i+1}, \ldots, a_{j-1}, a_{j}\}</math>.
* '''Input:''' two polynomials <math>P_1, P_2\in\mathbb{F}[x]</math> of degree <math>d</math>.
* '''Output:''' "yes" if two polynomials are identical, i.e. <math>P_1\equiv P_2</math>, and "no" if otherwise.
The <math>\mathbb{F}[x]</math> denote the [http://en.wikipedia.org/wiki/Polynomial_ring ring of polynomials] over field <math>\mathbb{F}</math>.


Alternatively, we can consider the following equivalent problem:
The '''conditional probability''' rules out the ''irrelevant'' events in a probabilistic argument.
* '''Input:''' a polynomial <math>P\in\mathbb{F}[x]</math> of degree <math>d</math>.
|}
* '''Output:''' "yes" if <math>P\equiv 0</math>, and "no" if otherwise.


The probalem is trivial if <math>P</math> is presented in its explicit form <math>P(x)=\sum_{i=0}^d a_ix^i</math>. But we assume that <math>P</math> is given in product form or as black box.
Summing all up:


A straightforward deterministic algorithm that solves PIT is to query <math>d+1</math> points <math>P(1),P(2),\ldots,P(d+1)</math> and check whether thay are all zero. This can determine whether <math>P\equiv 0</math> by interpolation.
:<math>\begin{align}
\mathbf{E}\left[\sum_{i=1}^n\sum_{j>i}X_{ij}\right]
&=
\sum_{i=1}^n\sum_{j>i}\mathbf{E}\left[X_{ij}\right]\\
&\le \sum_{i=1}^n\sum_{j>i}\frac{2}{j-i+1}\\
&= \sum_{i=1}^n\sum_{k=2}^{n-i+1}\frac{2}{k} & & (\mbox{Let }k=j-i+1)\\
&\le \sum_{i=1}^n\sum_{k=1}^{n}\frac{2}{k}\\
&= 2n\sum_{k=1}^{n}\frac{1}{k}\\
&= 2n H(n).
\end{align}</math>


We now introduce a simple randomized algorithm for the problem.
<math>H(n)</math> is the <math>n</math>th [http://en.wikipedia.org/wiki/Harmonic_number Harmonic number]. It holds that


{{Theorem|Algorithm for PIT|
:<math>\begin{align}H(n) = \ln n+O(1)\end{align}</math>.
*pick <math>x\in\{1,2,\ldots,2d\}</math> uniformly at random;
*if <math>P(x) = 0</math> then return “yes” else return “no”;
}}


This algorithm requires only the evaluation of <math>P</math> at a single point. And if <math>P\equiv 0</math> it is always correct.
Therefore, for an arbitrary input <math>S</math> of <math>n</math> numbers, the expected number of comparisons taken by RandQSort to sort <math>S</math> is <math>\mathrm{O}(n\log n)</math>.
And if <math>P\not\equiv 0</math> then the probability that the algorithm wrongly returns "yes" is bounded as follows.


{{Theorem|Theorem|
: Let <math>P\in\mathbb{F}[x]</math> be a polynomial of degree <math>d</math> over the field <math>\mathbb{F}</math>. Let <math>S\subset\mathbb{F}</math> be an arbitrary set and <math>x\in S</math> is chosen uniformly at random from <math>S</math>. If <math>P\not\equiv 0</math> then
::<math>\Pr[P(x)=0]\le\frac{d}{|S|}.</math>
}}
{{Proof|
A non-zero <math>d</math>-degree polynomial <math>P</math> has at most <math>d</math> distinct roots, thus at most <math>d</math> members <math>x</math> of <math>S</math> satisfy that <math>P(x)=0</math>. Therefore, <math>\Pr[P(x)=0]\le\frac{d}{|S|}</math>.
}}


By the theorem, the algorithm can distinguish a non-zero polynomial from 0 with probability at least <math>1/2</math>. This is achieved by evaluation of the polynomial at only one point and <math>1+\log_2 d</math> many random bits.
= Distributions of Coin Flips =
We introduce several important distributions induced by independent coin flips (independent probabilistic experiments), including: Bernoulli trial, geometric distribution, binomial distribution.


== Communication Complexity of Identity ==
==Bernoulli trial (Bernoulli distribution)==
The [http://en.wikipedia.org/wiki/Communication_complexity communication complexity] is introduced by Andrew Chi-Chih Yao as a model of computation which involves multiple participants, each with partial information of the input.
Bernoulli trial describes the probability distribution of a single (biased) coin flip. Suppose that we flip a (biased) coin where the probability of HEADS is <math>p</math>. Let <math>X</math> be the 0-1 random variable which indicates whether the result is HEADS. We say that <math>X</math> follows the Bernoulli distribution with parameter <math>p</math>. Formally,  
:<math>\begin{align}
X
&=
\begin{cases}
1 & \text{with probability }p\\
0 & \text{with probability }1-p
\end{cases}
\end{align}</math>.


Assume that there are two entities, say Alice and Bob. Alice has a private input <math>a</math> and Bob has a private input <math>b</math>. Together they want to compute a function <math>f(a,b)</math> by communicating with each other. The communication follows a predefined '''communication protocol''' (the "algorithm" in this model) whose logics depends only on the problem <math>f</math> but not on the inputs. The complexity of a communication protocol is measured by the number of bits communicated between Alice and Bob in the worst case.
==Geometric distribution==
Suppose we flip the same coin repeatedly until HEADS appears, where each coin flip is independent and follows the Bernoulli distribution with parameter <math>p</math>. Let <math>X</math> be the random variable denoting the total number of coin flips. Then <math>X</math> has the geometric distribution with parameter <math>p</math>. Formally, <math>\Pr[X=k]=(1-p)^{k-1}p</math>.


The problem of checking identity is formally defined by the function EQ as follows: <math>\mathrm{EQ}:\{0,1\}^n\times\{0,1\}^n\rightarrow\{0,1\}</math> and for any <math>a,b\in\{0,1\}^n</math>,
For geometric <math>X</math>, <math>\mathbf{E}[X]=\frac{1}{p}</math>. This can be verified by directly computing <math>\mathbf{E}[X]</math> by the definition of expectations. There is also a smarter way of computing  <math>\mathbf{E}[X]</math>, by using indicators and the linearity of expectations. For <math>k=0, 1, 2, \ldots</math>, let <math>Y_k</math> be the 0-1 random variable such that <math>Y_k=1</math> if and only if none of the first <math>k</math> coin flipings are HEADS, thus <math>\mathbf{E}[Y_k]=\Pr[Y_k=1]=(1-p)^{k}</math>. A key observation is that <math>X=\sum_{k=0}^\infty Y_k</math>. Thus, due to the linearity of expectations,
:<math>
:<math>
\mathrm{EQ}(a,b)=
\begin{align}
\begin{cases}
\mathbf{E}[X]
1& \mbox{if } a=b,\\
=
0& \mbox{otherwise.}
\mathbf{E}\left[\sum_{k=0}^\infty Y_k\right]
\end{cases}</math>
=
\sum_{k=0}^\infty \mathbf{E}[Y_k]
=
\sum_{k=0}^\infty (1-p)^k
=
\frac{1}{1-(1-p)}
=\frac{1}{p}.
\end{align}
</math>


A trivial way to solve EQ is to let Bob send his entire input string <math>b</math> to Alice and let Alice check whether <math>a=b</math>. This costs <math>n</math> bits of communications.
==Binomial distribution==
Suppose we flip the same (biased) coin for <math>n</math> times, where each coin flip is independent and follows the Bernoulli distribution with parameter <math>p</math>. Let <math>X</math> be the number of HEADS. Then <math>X</math> has the binomial distribution with parameters <math>n</math> and <math>p</math>. Formally, <math>\Pr[X=k]={n\choose k}p^k(1-p)^{n-k}</math>.


It is known that for deterministic communication protocols, this is the best we can get for computing EQ.
A binomial random variable <math>X</math> with parameters <math>n</math> and <math>p</math> is usually denoted by <math>B(n,p)</math>.


{{Theorem|Theorem (Yao 1979)|
As we saw above, by applying the linearity of expectations, it is easy to show that <math>\mathbf{E}[X]=np</math> for an <math>X=B(n,p)</math>.
:Any deterministic communication protocol computing EQ on two <math>n</math>-bit strings costs <math>n</math> bits of communication in the worst-case.
}}


This theorem is much more nontrivial to prove than it looks, because Alice and Bob are allowed to interact with each other in arbitrary ways. The proof of this theorem in Yao's 1979 paper initiates the field of communication complexity.
=Balls into Bins=


If the randomness is allowed, we can solve this problem up to a tolerable probabilistic error with significantly less communications. The inputs <math>a,b\in\{0,1\}^{n}</math> are two strings <math>a=a_0a_1\cdots a_{n-1}, b=b_0b_1\cdots b_{n-1}</math> of <math>n</math> bits. Let <math>k=\lceil\log_2 (2n)\rceil</math> and <math>p\in[2^k,2^{k+1}]</math> be an arbitrary prime number. (Such a prime <math>p</math> always exists.) The input strings <math>a,b</math> can be respectively represented as two polynomials <math>f,g\in\mathbb{Z}_p[x]</math> such that <math>f(x)=\sum_{i=0}^{n-1}a_ix^{i}</math> and <math>g(x)=\sum_{i=0}^{n-1}b_ix^{i}</math> of degree <math>n-1</math>, where all additions and multiplications are modulo <math>p</math>. The randomized communication protocol is given as follows:
== Birthday Problem==
{{Theorem|A randomized protocol for EQ|
There are <math>m</math> students in the class. Assume that for each student, his/her birthday is uniformly and independently distributed over the 365 days in a years. We wonder what the probability that no two students share a birthday.
'''Alice does''':
:* pick <math>x\in[p]</math> uniformly at random;
:* send <math>x</math> and <math>f(x)</math> to Bob;
'''Upon receiving''' <math>x</math> and <math>f(x)</math> '''Bob does''':
:* If <math>f(x)= g(x)</math> return "'''yes'''"; else return "'''no'''".
}}
Repeat this protocol for 100 times.
The total number of bits to communicate is bounded by <math>200\log_2p=O(\log n)</math>. Due to the analysis of the randomized algorithm for PIT, if <math>a=b</math> the protocol is always correct and if <math>a\neq b</math> the protocol fails to report a difference with probability less than <math>2^{-100}</math>.


Due to the [http://en.wikipedia.org/wiki/Pigeonhole_principle pigeonhole principle], it is obvious that for <math>m>365</math>, there must be two students with the same birthday. Surprisingly, for any <math>m>57</math> this event occurs with more than 99% probability. This is called the [http://en.wikipedia.org/wiki/Birthday_problem '''birthday paradox''']. Despite the name, the birthday paradox is not a real paradox.


=Polynomial Identity Testing (PIT)=
We can model this problem as a balls-into-bins problem. <math>m</math> different balls (students) are uniformly and independently thrown into 365 bins (days). More generally, let <math>n</math> be the number of bins. We ask for the probability of the following event <math>\mathcal{E}</math>
Consider the following problem of '''Polynomial Identity Testing (PIT)''':
* '''Input:''' two <math>n</math>-variate polynomials <math>f, g\in\mathbb{F}[x_1,x_2,\ldots,x_n]</math> of degree <math>d</math>.
* '''Output:''' "yes" if <math>f\equiv g</math>, and "no" if otherwise.
The <math>\mathbb{F}[x_1,x_2,\ldots,x_n]</math> is the [http://en.wikipedia.org/wiki/Polynomial_ring#The_polynomial_ring_in_several_variables of multi-variate polynomials] over field <math>\mathbb{F}</math>. The most natural way to represent an <math>n</math>-variate polynomial of degree <math>d</math> is to write it as a sum of monomials:
:<math>f(x_1,x_2,\ldots,x_n)=\sum_{i_1,i_2,\ldots,i_n\ge 0\atop i_1+i_2+\cdots+i_n\le d}a_{i_1,i_2,\ldots,i_n}x_{1}^{i_1}x_2^{i_2}\cdots x_{n}^{i_n}</math>.
The '''degree''' or '''total degree''' of a monomial <math>a_{i_1,i_2,\ldots,i_n}x_{1}^{i_1}x_2^{i_2}\cdots x_{n}^{i_n}</math> is given by <math>i_1+i_2+\cdots+i_n</math> and the degree of a polynomial <math>f</math> is the maximum degree of monomials of nonzero coefficients.


Alternatively, we can consider the following equivalent problem:
* <math>\mathcal{E}</math>: there is no bin with more than one balls (i.e. no two students share birthday).
* '''Input:''' a polynomial <math>f\in\mathbb{F}[x_1,x_2,\ldots,x_n]</math> of degree <math>d</math>.
* '''Output:''' "yes" if <math>f\equiv 0</math>, and "no" if otherwise.


If <math>f</math> is written explicitly as a sum of monomials, then the problem is trivial. Again we allow <math>f</math> to be represented in product form.
We first analyze this by counting. There are totally <math>n^m</math> ways of assigning <math>m</math> balls to <math>n</math> bins. The number of assignments that no two balls share a bin is <math>{n\choose m}m!</math>.


{{Theorem|Example|
Thus the probability is given by:
The [http://en.wikipedia.org/wiki/Vandermonde_matrix Vandermonde matrix] <math>M=M(x_1,x_2,\ldots,x_n)</math> is defined as that <math>M_{ij}=x_i^{j-1}</math>, that is
:<math>\begin{align}
:<math>M=\begin{bmatrix}
\Pr[\mathcal{E}]
1 & x_1 & x_1^2 & \dots & x_1^{n-1}\\
=
1 & x_2 & x_2^2 & \dots & x_2^{n-1}\\
\frac{{n\choose m}m!}{n^m}.
1 & x_3 & x_3^2 & \dots & x_3^{n-1}\\
\end{align}
\vdots & \vdots & \vdots & \ddots &\vdots \\
1 & x_n & x_n^2 & \dots & x_n^{n-1}
\end{bmatrix}</math>.
Let <math>f</math> be the polynomial defined as
:<math>
f(x_1,\ldots,x_n)=\det(M)=\prod_{j<i}(x_i-x_j).
</math>
</math>
It is pretty easy to evaluate <math>f(x_1,x_2,\ldots,x_n)</math> on any particular <math>x_1,x_2,\ldots,x_n</math>, however it is prohibitively expensive to symbolically expand <math>f(x_1,\ldots,x_n)</math> to its sum-of-monomial form.
}}


== Schwartz-Zippel Theorem==
Recall that <math>{n\choose m}=\frac{n!}{(n-m)!m!}</math>. Then
Here is a very simple randomized algorithm, due to Schwartz and Zippel.
:<math>\begin{align}
{{Theorem|Randomized algorithm for multi-variate PIT|
\Pr[\mathcal{E}]
* fix an arbitrary set <math>S\subseteq \mathbb{F}</math> whose size to be fixed;
=
* pick <math>r_1,r_2,\ldots,r_n\in S</math> uniformly and independently at random;
\frac{{n\choose m}m!}{n^m}
* if <math>f(\vec{r})=f(r_1,r_2,\ldots,r_n) = 0</math> then return “yes” else return “no”;
=
}}
\frac{n!}{n^m(n-m)!}
=
\frac{n}{n}\cdot\frac{n-1}{n}\cdot\frac{n-2}{n}\cdots\frac{n-(m-1)}{n}
=
\prod_{k=1}^{m-1}\left(1-\frac{k}{n}\right).
\end{align}
</math>


This algorithm requires only the evaluation of <math>f</math> at a single point <math>\vec{r}</math>. And if <math>f\equiv 0</math> it is always correct.
There is also a more "probabilistic" argument for the above equation. To be rigorous, we need the following theorem, which holds generally and is very useful for computing the AND of many events.


In the Theorem below, we’ll see that if <math>f\not\equiv 0</math> then the algorithm is incorrect with probability at most <math>\frac{d}{|S|}</math>, where <math>d</math> is the degree of the polynomial <math>f</math>.
:::{|border="1"
|By the definition of conditional probability, <math>\Pr[A\mid B]=\frac{\Pr[A\wedge B]}{\Pr[B]}</math>. Thus, <math>\Pr[A\wedge B] =\Pr[B]\cdot\Pr[A\mid B]</math>. This hints us that we can compute the probability of the AND of events by conditional probabilities. Formally, we have the following theorem:
'''Theorem:'''
:Let <math>\mathcal{E}_1, \mathcal{E}_2, \ldots, \mathcal{E}_n</math>  be any <math>n</math> events. Then
::<math>\begin{align}
\Pr\left[\bigwedge_{i=1}^n\mathcal{E}_i\right]
&=
\prod_{k=1}^n\Pr\left[\mathcal{E}_k \mid \bigwedge_{i<k}\mathcal{E}_i\right].
\end{align}</math>
'''Proof:'''  It holds that <math>\Pr[A\wedge B] =\Pr[B]\cdot\Pr[A\mid B]</math>. Thus, let <math>A=\mathcal{E}_n</math> and <math>B=\mathcal{E}_1\wedge\mathcal{E}_2\wedge\cdots\wedge\mathcal{E}_{n-1}</math>, then
:<math>\begin{align}
\Pr[\mathcal{E}_1\wedge\mathcal{E}_2\wedge\cdots\wedge\mathcal{E}_n]
&=
\Pr[\mathcal{E}_1\wedge\mathcal{E}_2\wedge\cdots\wedge\mathcal{E}_{n-1}]\cdot\Pr\left[\mathcal{E}_n\mid \bigwedge_{i<n}\mathcal{E}_i\right].
\end{align}
</math>
Recursively applying this equation to <math>\Pr[\mathcal{E}_1\wedge\mathcal{E}_2\wedge\cdots\wedge\mathcal{E}_{n-1}]</math> until there is only <math>\mathcal{E}_1</math> left, the theorem is proved. <math>\square</math>
|}


{{Theorem|Schwartz-Zippel Theorem|
Now we are back to the probabilistic analysis of the birthday problem, with a general setting of <math>m</math> students and <math>n</math> possible birthdays (imagine that we live in a planet where a year has <math>n</math> days).
: Let <math>f\in\mathbb{F}[x_1,x_2,\ldots,x_n]</math> be a multivariate polynomial of degree <math>d</math> over a field <math>\mathbb{F}</math> such that <math>f\not\equiv 0</math>. Fix any finite set <math>S\subset\mathbb{F}</math>, and let <math>r_1,r_2\ldots,r_n</math> be chosen uniformly and independently at random from <math>S</math>. Then
::<math>\Pr[f(r_1,r_2,\ldots,r_n)=0]\le\frac{d}{|S|}.</math>
}}
{{Proof|
We prove by induction on <math>n</math> the number of variables.


For <math>n=1</math>, assuming that <math>f\not\equiv 0</math>, due to the fundamental theorem of algebra, the degree-<math>d</math> polynomial <math>f(x)</math> has at most <math>d</math> roots, thus
The first student has a birthday (of course!). The probability that the second student has a different birthday is <math>\left(1-\frac{1}{n}\right)</math>. Given that the first two students have different birthdays, the probability that the third student has a different birthday from the first two is <math>\left(1-\frac{2}{n}\right)</math>. Continuing this on, assuming that the first <math>k-1</math> students all have different birthdays, the probability that the <math>k</math>th student has a different birthday than the first <math>k-1</math>, is given by <math>\left(1-\frac{k-1}{n}\right)</math>. So the probability that all <math>m</math> students have different birthdays is the product of all these conditional probabilities:
:<math>\Pr[f(r)=0]\le\frac{d}{|S|}.
:<math>\begin{align}
\Pr[\mathcal{E}]=\left(1-\frac{1}{n}\right)\cdot \left(1-\frac{2}{n}\right)\cdots \left(1-\frac{m-1}{n}\right)
&=
\prod_{k=1}^{m-1}\left(1-\frac{k}{n}\right),
\end{align}
</math>
</math>
which is the same as what we got by the counting argument.


Assume the induction hypothesis for a multi-variate polynomial up to <math>n-1</math> variable.
[[File:Birthday.png|border|450px|right]]


An <math>n</math>-variate polynomial <math>f(x_1,x_2,\ldots,x_n)</math> can be represented as
There are several ways of analyzing this formular. Here is a convenient one: Due to [http://en.wikipedia.org/wiki/Taylor_series Taylor's expansion], <math>e^{-k/n}\approx 1-k/n</math>. Then
:<math>f(x_1,x_2,\ldots,x_n)=\sum_{i=0}^kx_n^{i}f_i(x_1,x_2,\ldots,x_{n-1})</math>,
:<math>\begin{align}
where <math>k</math> is the largest power of <math>x_n</math>, which means that the degree of <math>f_k</math> is at most <math>d-k</math> and <math>f_k\not\equiv 0</math>.
\prod_{k=1}^{m-1}\left(1-\frac{k}{n}\right)
&\approx
\prod_{k=1}^{m-1}e^{-\frac{k}{n}}\\
&=
\exp\left(-\sum_{k=1}^{m-1}\frac{k}{n}\right)\\
&=
e^{-m(m-1)/2n}\\
&\approx
e^{-m^2/2n}.
\end{align}</math>
The quality of this approximation is shown in the Figure.


In particular, we write <math>f</math> as a sum of two parts:
Therefore, for <math>m=\sqrt{2n\ln \frac{1}{\epsilon}}</math>, the probability that <math>\Pr[\mathcal{E}]\approx\epsilon</math>.
:<math>f(x_1,x_2,\ldots,x_n)=x_n^k f_k(x_1,x_2,\ldots,x_{n-1})+\bar{f}(x_1,x_2,\ldots,x_n)</math>,
where both <math>f_k</math> and <math>\bar{f}</math> are polynomials, such that  
* as argued above, the degree of <math>f_k</math> is at most <math>d-k</math> and <math>f_k\not\equiv 0</math>;
* <math>\bar{f}(x_1,x_2,\ldots,x_n)=\sum_{i=0}^{k-1}x_n^i f_i(x_1,x_2,\ldots,x_{n-1})</math>, thus <math>\bar{f}(x_1,x_2,\ldots,x_n)</math> has no <math>x_n^{k}</math> factor in any term.


By the law of total probability, it holds that
==Coupon Collector ==
Suppose that a chocolate company releases <math>n</math> different types of coupons. Each box of chocolates contains one coupon with a uniformly random type. Once you have collected all <math>n</math> types of coupons, you will get a prize. So how many boxes of chocolates you are expected to buy to win the prize?
 
The coupon collector problem can be described in the balls-into-bins model as follows. We keep throwing balls one-by-one into <math>n</math> bins (coupons), such that each ball is thrown into a bin uniformly and independently at random. Each ball corresponds to a box of chocolate, and each bin corresponds to a type of coupon. Thus, the number of boxes bought to collect <math>n</math> coupons is just the number of balls thrown until none of the <math>n</math> bins is empty.
 
{{Theorem
|Theorem|
:Let <math>X</math> be the number of balls thrown uniformly and independently to <math>n</math> bins until no bin is empty. Then <math>\mathbf{E}[X]=nH(n)</math>, where <math>H(n)</math> is the <math>n</math>th harmonic number.
}}
{{Proof| Let <math>X_i</math> be the number of balls thrown while there are ''exactly'' <math>i-1</math> nonempty bins, then clearly <math>X=\sum_{i=1}^n X_i</math>.
 
When there are exactly <math>i-1</math> nonempty bins, throwing a ball, the probability that the number of nonempty bins increases  (i.e. the ball is thrown to an empty bin) is
:<math>p_i=1-\frac{i-1}{n}.
</math>
<math>X_i</math> is the number of balls thrown to make the number of nonempty bins increases from <math>i-1</math> to <math>i</math>, i.e. the number of balls thrown until a ball is thrown to a current empty bin. Thus, <math>X_i</math> follows the [http://en.wikipedia.org/wiki/Geometric_distribution geometric distribution], such that  
:<math>\Pr[X_i=k]=(1-p_i)^{k-1}p_i</math>
 
For a geometric random variable, <math>\mathbf{E}[X_i]=\frac{1}{p_i}=\frac{n}{n-i+1}</math>.
 
Applying the linearity of expectations,
:<math>
:<math>
\begin{align}
\begin{align}
&\Pr[f(r_1,r_2,\ldots,r_n)=0]\\
\mathbf{E}[X]
=
&=
&\Pr[f(\vec{r})=0\mid f_k(r_1,r_2,\ldots,r_{n-1})=0]\cdot\Pr[f_k(r_1,r_2,\ldots,r_{n-1})=0]\\
\mathbf{E}\left[\sum_{i=1}^nX_i\right]\\
&+\Pr[f(\vec{r})=0\mid f_k(r_1,r_2,\ldots,r_{n-1})\neq0]\cdot\Pr[f_k(r_1,r_2,\ldots,r_{n-1})\neq0].
&=
\sum_{i=1}^n\mathbf{E}\left[X_i\right]\\
&=
\sum_{i=1}^n\frac{n}{n-i+1}\\
&=
n\sum_{i=1}^n\frac{1}{i}\\
&=
nH(n),
\end{align}
\end{align}
</math>
</math>
Note that <math>f_k(r_1,r_2,\ldots,r_{n-1})</math> is a polynomial on <math>n-1</math> variables of degree <math>d-k</math> such that <math>f_k\not\equiv 0</math>.
where <math>H(n)</math> is the <math>n</math>th Harmonic number, and <math>H(n)=\ln n+O(1)</math>. Thus, for the coupon collectors problem, the expected number of coupons required to obtain all <math>n</math> types of coupons is <math>n\ln n+O(n)</math>.
By the induction hypothesis, we have
}}
:<math>
 
\begin{align}
----
(*)
 
&\qquad
Only knowing the expectation is not good enough. We would like to know how fast the probability decrease as a random variable deviates from its mean value.
&\Pr[f_k(r_1,r_2,\ldots,r_{n-1})=0]\le\frac{d-k}{|S|}.
 
\end{align}
{{Theorem
|Theorem|
:Let <math>X</math> be the number of balls thrown uniformly and independently to <math>n</math> bins until no bin is empty. Then <math>\Pr[X\ge n\ln n+cn]<e^{-c}</math> for any <math>c>0</math>.
}}
{{Proof| For any particular bin <math>i</math>, the probability that bin <math>i</math> is empty after throwing <math>n\ln n+cn</math> balls is
:<math>\left(1-\frac{1}{n}\right)^{n\ln n+cn}
< e^{-(\ln n+c)}
=\frac{1}{ne^c}.
</math>
</math>


For the second case, recall that <math>\bar{f}(x_1,\ldots,x_n)</math> has no <math>x_n^k</math> factor in any term, thus the condition <math>f_k(r_1,r_2,\ldots,r_{n-1})\neq0</math> guarantees that  
By the union bound, the probability that there exists an empty bin after throwing <math>n\ln n+cn</math> balls is
:<math>f(r_1,\ldots,r_{n-1},x_n)=x_n^k f_k(r_1,r_2,\ldots,r_{n-1})+\bar{f}(r_1,r_2,\ldots,r_n)=g_{r_1,\ldots,r_{n-1}}(x_n)</math>
is a single-variate polynomial such that the degree of <math>g_{r_1,\ldots,r_{n-1}}(x_n)</math> is <math>k</math> and <math>g_{r_1,\ldots,r_{n-1}}\not\equiv 0</math>, for which we already known that the probability <math>g_{r_1,\ldots,r_{n-1}}(r_n)=0</math> is at most <math>\frac{k}{|S|}</math>.
Therefore,
:<math>
:<math>
\begin{align}
\Pr[X\ge n\ln n+cn]
(**)
< n\cdot \frac{1}{ne^c}
&\qquad
=e^{-c}.
&\Pr[f(\vec{r})=0\mid f_k(r_1,r_2,\ldots,r_{n-1})\neq0]=\Pr[g_{r_1,\ldots,r_{n-1}}(r_n)=0\mid f_k(r_1,r_2,\ldots,r_{n-1})\neq0]\le\frac{k}{|S|}
\end{align}
</math>.
Substituting both <math>(*)</math> and <math>(**)</math> back in the total probability, we have
:<math>
\Pr[f(r_1,r_2,\ldots,r_n)=0]
\le\frac{d-k}{|S|}+\frac{k}{|S|}=\frac{d}{|S|},
</math>
</math>
which proves the theorem.
}}


------
=== Stable Marriage ===
We now consider the famous [http://en.wikipedia.org/wiki/Stable_marriage_problem '''stable marriage problem'''] or '''stable matching problem''' (SMP). This problem captures two aspects: allocations (matchings) and stability, two central topics in economics.


In above proof, for the second case that <math>f_k(r_1,\ldots,r_{n-1})\neq 0</math>, we use an "probabilistic arguement" to deal with the random choices in the condition. Here we give a more rigorous proof by enumerating all elementary events in applying the law of total probability. You make your own judgement which proof is better.
An instance of stable marriage consists of:
* <math>n</math> men and <math>n</math> women;
* each person associated with a strictly ordered ''preference list'' containing all the members of the opposite sex.
Formally, let <math>M</math> be the set of <math>n</math> men and <math>W</math> be the set of <math>n</math> women. Each man <math>m\in M</math> is associated with a permutation <math>p_m</math> of elemets in <math>W</math> and each woman <math>w\in W</math> is associated with a permutation <math>p_w</math> of elements in <math>M</math>.


By the law of total probability,
A ''matching'' is a one-one correspondence <math>\phi:M\rightarrow W</math>. We said a man <math>m</math> and a woman <math>w</math> are ''partners'' in <math>\phi</math> if <math>w=\phi(m)</math>.
:<math>
{{Theorem|Definition (stable matching)|
\begin{align}
:A pair <math>(m,w)</math> of a man and woman is a '''blocking pair''' in a matching <math>\phi</math> if <math>m</math> and <math>w</math> are not partners in <math>\phi</math> but
&\Pr[f(\vec{r})=0]\\
:* <math>m</math> prefers <math>w</math> to <math>\phi(m)</math>, and
:* <math>w</math> prefers <math>m</math> to <math>\phi(w)</math>.
:A matching <math>\phi</math> is '''stable''' if there is no blocking pair in it.
}}
 
It is unclear from the definition itself whether stable matchings always exist, and how to efficiently find a stable matching. Both questions are answered by the following proposal algorithm due to Gale and Shapley.
{{Theorem|The proposal algorithm (Gale-Shapley 1962)|
: Initially, all person are not married;
: in each step (called a '''proposal'''):
:* an arbitrary unmarried man <math>m</math> proposes to the woman <math>w</math> who is ranked highest in his preference list <math>p_m</math> among all the women who has not yet rejected <math>m</math>;
:* if <math>w</math> is still single then <math>w</math> accepts the proposal and is married to <math>m</math>;
:* if <math>w</math> is married to another man <math>m'</math> who is ranked lower than <math>m</math> in her preference list <math>p_w</math> then <math>w</math> divorces <math>m'</math> (thus <math>m'</math> becomes single again and considers himself as rejected by <math>w</math>) and is married to <math>m</math>;
:* if otherwise <math>w</math> rejects <math>m</math>;
}}
 
The algorithm terminates when the last single woman receives a proposal. Since for every pair <math>(m,w)\in M\times W</math> of man and woman, <math>m</math> proposes to <math>w</math> at most once.
The algorithm terminates in at most <math>n^2</math> proposals in the worst case.
 
It is obvious to see that the algorithm retruns a macthing, and this matching must be stable. To see this, by contradiction suppose that the algorithm resturns a macthing <math>\phi</math>, such that two men <math>A, B</math> are macthed to two women <math>a,b</math> in <math>\phi</math> respectively, but <math>A</math> and <math>b</math> prefers each other to their partners <math>a</math> and <math>B</math> respectively. By definition of the algorithm, <math>A</math> would have proposed to <math>b</math> before proposing to <math>a</math>, by which time <math>b</math> must either be single or be matched to a man ranked lower than <math>A</math> in her list (because her final partner <math>B</math> is ranked lower than <math>A</math>), which means <math>b</math> must have accepted <math>A</math>'s proposal, a contradiction.
 
 
We are interested in the average-case performance of this algorithm, that is, the expected number of proposals if everyone's preference list is a uniformly and independently random permutation.
 
The following '''principle of deferred decisions''' is quite useful in analysing performance of algorithm with random input.
{{Theorem|Principle of deferred decisions|
:The decision of random choice in the random input can be deferred to the running time of the algorithm.
}}
 
Apply the principle of deferred decisions, the deterministic proposal algorithm with random permutations as input is equivalent to the following random process:
* At each step, a man <math>m</math> choose a woman <math>w</math> uniformly and independently at random to propose, among all the women who have not rejected him yet. ('''sample without replacement''')
 
We then compare the above process with the following modified process:
* The man <math>m</math> repeatedly samples a uniform and independent woman to propose among all women, until he successfully samples a woman who has not rejected him and propose to her. ('''sample with replacement''')
 
It is easy to see that the modified process (sample with replacement) is no more efficient than the original process (sample without replacement) because it simulates the original process if at each step we only count the last proposal to the woman who has not rejected the man. Such comparison of two random processes by forcing them to be related in some way is called [http://en.wikipedia.org/wiki/Coupling_(probability) coupling].
 
Note that in the modified process (sample with replacement), each proposal, no matter from which man, is going to a uniformly and independently random women. And we know that the algorithm terminated once the last single woman receives a proposal, i.e. once all <math>n</math> women have received at least one proposal. This is the coupon collector problem with proposals as balls (cookie boxes) and women as bins (coupons).
Due to our analysis of the coupon collector problem, the expected number of proposals is bounded by <math>O(n\ln n)</math>.
 
== Occupancy Problem ==
Now we ask about the loads of bins. Assuming that <math>m</math> balls are uniformly and independently assigned to <math>n</math> bins, for <math>1\le i\le n</math>, let <math>X_i</math> be the '''load''' of the <math>i</math>th bin, i.e. the number of balls in the <math>i</math>th bin.
 
An easy analysis shows that for every bin <math>i</math>, the expected load <math>\mathbf{E}[X_i]</math> is equal to the average load <math>m/n</math>.
 
Because there are totally <math>m</math> balls, it is always true that <math>\sum_{i=1}^n X_i=m</math>.
 
Therefore, due to the linearity of expectations,
:<math>\begin{align}
\sum_{i=1}^n\mathbf{E}[X_i]
&=
\mathbf{E}\left[\sum_{i=1}^n X_i\right]
=
=
&\sum_{x_1,\ldots,x_{n-1}\in S}\Pr[f(\vec{r})=0\mid \forall i<n, r_i=x_i]\cdot\Pr[\forall i<n, r_i=x_i]\\
\mathbf{E}\left[m\right]
=
=m.
&\sum_{x_1,\ldots,x_{n-1}\in S\atop f_k(x_1,\ldots,x_{n-1})=0}\Pr[f(\vec{r})=0\mid \forall i<n, r_i=x_i]\cdot\Pr[\forall i<n, r_i=x_i]\\
\end{align}</math>
&+\sum_{x_1,\ldots,x_{n-1}\in S\atop f_k(x_1,\ldots,x_{n-1})\neq0}\Pr[f(\vec{r})=0\mid \forall i<n, r_i=x_i]\cdot\Pr[\forall i<n, r_i=x_i]\\
 
\le
Because for each ball, the bin to which the ball is assigned is uniformly and independently chosen, the distributions of the loads of bins are identical. Thus <math>\mathbf{E}[X_i]</math> is the same for each <math>i</math>. Combining with the above equation, it holds that for every <math>1\le i\le m</math>, <math>\mathbf{E}[X_i]=\frac{m}{n}</math>. So the average is indeed the average!
&\sum_{x_1,\ldots,x_{n-1}\in S\atop f_k(x_1,\ldots,x_{n-1})=0}\Pr[\forall i<n, r_i=x_i]\\
 
&+\sum_{x_1,\ldots,x_{n-1}\in S\atop f_k(x_1,\ldots,x_{n-1})\neq 0}\Pr[f(x_1,\ldots,x_{n-1},r_n)=0\mid \forall i<n, r_i=x_i]\cdot\Pr[\forall i<n, r_i=x_i]\\
----
=
 
&\Pr[f_k(r_1,\ldots,r_{n-1})=0]+\sum_{x_1,\ldots,x_{n-1}\in S\atop f_k(x_1,\ldots,x_{n-1})\neq 0}\Pr[f(x_1,\ldots,x_{n-1},r_n)=0]\cdot\Pr[\forall i<n, r_i=x_i].
Next we analyze the distribution of the maximum load. We show that when <math>m=n</math>, i.e. <math>n</math> balls are uniformly and independently thrown into <math>n</math> bins, the maximum load is <math>O\left(\frac{\log n}{\log\log n}\right)</math> with high probability.
 
{{Theorem
|Theorem|
:Suppose that <math>n</math> balls are thrown independently and uniformly at random into <math>n</math> bins. For <math>1\le i\le n</math>, let <math>X_i</math> be the random variable denoting the number of balls in the <math>i</math>th bin. Then
::<math>\Pr\left[\max_{1\le i\le n}X_i \ge\frac{3\ln n}{\ln\ln n}\right] <\frac{1}{n}.</math>
}}
 
{{Proof| Let <math>M</math> be an integer. Take bin 1. For any particular <math>M</math> balls, these <math>M</math> balls are all thrown to bin 1 with probability <math>(1/n)^M</math>, and there are totally <math>{n\choose M}</math> distinct sets of <math>M</math> balls. Therefore, applying the union bound,
:<math>\begin{align}\Pr\left[X_1\ge M\right]
&\le
{n\choose M}\left(\frac{1}{n}\right)^M\\
&=
\frac{n!}{M!(n-M)!n^M}\\
&=
\frac{1}{M!}\cdot\frac{n(n-1)(n-2)\cdots(n-M+1)}{n^M}\\
&=
\frac{1}{M!}\cdot \prod_{i=0}^{M-1}\left(1-\frac{i}{n}\right)\\
&\le \frac{1}{M!}.
\end{align}</math>
 
According to [http://en.wikipedia.org/wiki/Stirling's_approximation Stirling's approximation], <math>M!\approx \sqrt{2\pi M}\left(\frac{M}{e}\right)^M</math>, thus
:<math>\frac{1}{M!}\le\left(\frac{e}{M}\right)^M.</math>
 
[[file:Balls2bins.png|frame|Figure 1]]
 
Due to the symmetry. All <math>X_i</math> have the same distribution.
Apply the union bound again,  
:<math>\begin{align}
\Pr\left[\max_{1\le i\le n}X_i\ge M\right]
&=
\Pr\left[(X_1\ge M) \vee (X_2\ge M) \vee\cdots\vee (X_n\ge M)\right]\\
&\le
n\Pr[X_1\ge M]\\
&\le n\left(\frac{e}{M}\right)^M.
\end{align}
\end{align}
</math>
</math>
We have argued that <math>f_k\not\equiv 0</math> and the degree of <math>f_k</math> is <math>d-k</math>. By the induction hypothesis, we have
 
:<math>
When <math>M=3\ln n/\ln\ln n</math>,
\Pr[f_k(r_1,\ldots,r_{n-1})=0]\le\frac{d-k}{|S|}.
:<math>\begin{align}
\left(\frac{e}{M}\right)^M
&=
\left(\frac{e\ln\ln n}{3\ln n}\right)^{3\ln n/\ln\ln n}\\
&<
\left(\frac{\ln\ln n}{\ln n}\right)^{3\ln n/\ln\ln n}\\
&=
e^{3(\ln\ln\ln n-\ln\ln n)\ln n/\ln\ln n}\\
&=
e^{-3\ln n+3\ln\ln\ln n\ln n/\ln\ln n}\\
&\le  
e^{-2\ln n}\\
&=
\frac{1}{n^2}.
\end{align}
</math>
</math>
And for every fixed <math>x_1,\ldots,x_{n-1}\in S</math> such that <math>f_k(x_1,\ldots,x_{n-1})\neq 0</math>, we have argued that <math>f(x_1,\ldots,x_{n-1},x_n)</math> is a polynomial in <math>x_n</math> of degree <math>k</math>, thus
 
:<math>
Therefore,
\Pr[f(x_1,\ldots,x_{n-1},r_n)=0]\le\frac{k}{|S|},
:<math>\begin{align}
</math>
\Pr\left[\max_{1\le i\le n}X_i\ge \frac{3\ln n}{\ln\ln n}\right]
which holds for all <math>x_1,\ldots,x_{n-1}\in S</math> such that <math>f_k(x_1,\ldots,x_{n-1})\neq 0</math>, therefore the weighted average
&\le n\left(\frac{e}{M}\right)^M
:<math>
&< \frac{1}{n}.
\sum_{x_1,\ldots,x_{n-1}\in S\atop f_k(x_1,\ldots,x_{n-1})\neq 0}\Pr[f(x_1,\ldots,x_{n-1},r_n)=0]\cdot\Pr[\forall i<n, r_i=x_i]
\end{align}
\le\frac{k}{|S|}.
</math>
Substituting these inequalities back to the total probability, we have
<math>
\Pr[f(\vec{r})=0]
\le\frac{d-k}{|S|}+\frac{k}{|S|}
=\frac{d}{|S|}.
</math>
</math>
}}
}}
When <math>m>n</math>, Figure 1 illustrates the results of several random experiments, which show that the distribution of the loads of bins becomes more even as the number of balls grows larger than the number of bins.
Formally, it can be proved that for <math>m=\Omega(n\log n)</math>, with high probability, the maximum load is within <math>O\left(\frac{m}{n}\right)</math>, which is asymptotically equal to the average load.

Revision as of 06:41, 3 March 2014

Random Variable

Definition (random variable)
A random variable [math]\displaystyle{ X }[/math] on a sample space [math]\displaystyle{ \Omega }[/math] is a real-valued function [math]\displaystyle{ X:\Omega\rightarrow\mathbb{R} }[/math]. A random variable X is called a discrete random variable if its range is finite or countably infinite.

For a random variable [math]\displaystyle{ X }[/math] and a real value [math]\displaystyle{ x\in\mathbb{R} }[/math], we write "[math]\displaystyle{ X=x }[/math]" for the event [math]\displaystyle{ \{a\in\Omega\mid X(a)=x\} }[/math], and denote the probability of the event by

[math]\displaystyle{ \Pr[X=x]=\Pr(\{a\in\Omega\mid X(a)=x\}) }[/math].

The independence can also be defined for variables:

Definition (Independent variables)
Two random variables [math]\displaystyle{ X }[/math] and [math]\displaystyle{ Y }[/math] are independent if and only if
[math]\displaystyle{ \Pr[(X=x)\wedge(Y=y)]=\Pr[X=x]\cdot\Pr[Y=y] }[/math]
for all values [math]\displaystyle{ x }[/math] and [math]\displaystyle{ y }[/math]. Random variables [math]\displaystyle{ X_1, X_2, \ldots, X_n }[/math] are mutually independent if and only if, for any subset [math]\displaystyle{ I\subseteq\{1,2,\ldots,n\} }[/math] and any values [math]\displaystyle{ x_i }[/math], where [math]\displaystyle{ i\in I }[/math],
[math]\displaystyle{ \begin{align} \Pr\left[\bigwedge_{i\in I}(X_i=x_i)\right] &= \prod_{i\in I}\Pr[X_i=x_i]. \end{align} }[/math]

Note that in probability theory, the "mutual independence" is not equivalent with "pair-wise independence", which we will learn in the future.

Expectation

Let [math]\displaystyle{ X }[/math] be a discrete random variable. The expectation of [math]\displaystyle{ X }[/math] is defined as follows.

Definition (Expectation)
The expectation of a discrete random variable [math]\displaystyle{ X }[/math], denoted by [math]\displaystyle{ \mathbf{E}[X] }[/math], is given by
[math]\displaystyle{ \begin{align} \mathbf{E}[X] &= \sum_{x}x\Pr[X=x], \end{align} }[/math]
where the summation is over all values [math]\displaystyle{ x }[/math] in the range of [math]\displaystyle{ X }[/math].

Linearity of Expectation

Perhaps the most useful property of expectation is its linearity.

Theorem (Linearity of Expectations)
For any discrete random variables [math]\displaystyle{ X_1, X_2, \ldots, X_n }[/math], and any real constants [math]\displaystyle{ a_1, a_2, \ldots, a_n }[/math],
[math]\displaystyle{ \begin{align} \mathbf{E}\left[\sum_{i=1}^n a_iX_i\right] &= \sum_{i=1}^n a_i\cdot\mathbf{E}[X_i]. \end{align} }[/math]
Proof.
By the definition of the expectations, it is easy to verify that (try to prove by yourself):

for any discrete random variables [math]\displaystyle{ X }[/math] and [math]\displaystyle{ Y }[/math], and any real constant [math]\displaystyle{ c }[/math],

  • [math]\displaystyle{ \mathbf{E}[X+Y]=\mathbf{E}[X]+\mathbf{E}[Y] }[/math];
  • [math]\displaystyle{ \mathbf{E}[cX]=c\mathbf{E}[X] }[/math].

The theorem follows by induction.

[math]\displaystyle{ \square }[/math]

The linearity of expectation gives an easy way to compute the expectation of a random variable if the variable can be written as a sum.

Example
Supposed that we have a biased coin that the probability of HEADs is [math]\displaystyle{ p }[/math]. Flipping the coin for n times, what is the expectation of number of HEADs?
It looks straightforward that it must be np, but how can we prove it? Surely we can apply the definition of expectation to compute the expectation with brute force. A more convenient way is by the linearity of expectations: Let [math]\displaystyle{ X_i }[/math] indicate whether the [math]\displaystyle{ i }[/math]-th flip is HEADs. Then [math]\displaystyle{ \mathbf{E}[X_i]=1\cdot p+0\cdot(1-p)=p }[/math], and the total number of HEADs after n flips is [math]\displaystyle{ X=\sum_{i=1}^{n}X_i }[/math]. Applying the linearity of expectation, the expected number of HEADs is:
[math]\displaystyle{ \mathbf{E}[X]=\mathbf{E}\left[\sum_{i=1}^{n}X_i\right]=\sum_{i=1}^{n}\mathbf{E}[X_i]=np }[/math].

The real power of the linearity of expectations is that it does not require the random variables to be independent, thus can be applied to any set of random variables. For example:

[math]\displaystyle{ \mathbf{E}\left[\alpha X+\beta X^2+\gamma X^3\right] = \alpha\cdot\mathbf{E}[X]+\beta\cdot\mathbf{E}\left[X^2\right]+\gamma\cdot\mathbf{E}\left[X^3\right]. }[/math]

However, do not exaggerate this power!

  • For an arbitrary function [math]\displaystyle{ f }[/math] (not necessarily linear), the equation [math]\displaystyle{ \mathbf{E}[f(X)]=f(\mathbf{E}[X]) }[/math] does not hold generally.
  • For variances, the equation [math]\displaystyle{ var(X+Y)=var(X)+var(Y) }[/math] does not hold without further assumption of the independence of [math]\displaystyle{ X }[/math] and [math]\displaystyle{ Y }[/math].

Conditional Expectation

Conditional expectation can be accordingly defined:

Definition (conditional expectation)
For random variables [math]\displaystyle{ X }[/math] and [math]\displaystyle{ Y }[/math],
[math]\displaystyle{ \mathbf{E}[X\mid Y=y]=\sum_{x}x\Pr[X=x\mid Y=y], }[/math]
where the summation is taken over the range of [math]\displaystyle{ X }[/math].

There is also a law of total expectation.

Theorem (law of total expectation)
Let [math]\displaystyle{ X }[/math] and [math]\displaystyle{ Y }[/math] be two random variables. Then
[math]\displaystyle{ \mathbf{E}[X]=\sum_{y}\mathbf{E}[X\mid Y=y]\cdot\Pr[Y=y]. }[/math]

Random Quicksort

Given as input a set [math]\displaystyle{ S }[/math] of [math]\displaystyle{ n }[/math] numbers, we want to sort the numbers in [math]\displaystyle{ S }[/math] in increasing order. One of the most famous algorithm for this problem is the Quicksort algorithm.

  • if [math]\displaystyle{ |S|\gt 1 }[/math] do:
    • pick an [math]\displaystyle{ x\in S }[/math] as the pivot;
    • partition [math]\displaystyle{ S }[/math] into [math]\displaystyle{ S_1 }[/math], [math]\displaystyle{ \{x\} }[/math], and [math]\displaystyle{ S_2 }[/math], where all numbers in [math]\displaystyle{ S_1 }[/math] are smaller than [math]\displaystyle{ x }[/math] and all numbers in [math]\displaystyle{ S_2 }[/math] are larger than [math]\displaystyle{ x }[/math];
    • recursively sort [math]\displaystyle{ S_1 }[/math] and [math]\displaystyle{ S_2 }[/math];

The time complexity of this sorting algorithm is measured by the number of comparisons.

For the deterministic quicksort algorithm, the pivot is picked from a fixed position (e.g. the first number in the array). The worst-case time complexity in terms of number of comparisons is [math]\displaystyle{ \Theta(n^2) }[/math].

We consider the following randomized version of the quicksort.

  • if [math]\displaystyle{ |S|\gt 1 }[/math] do:
    • uniformly pick a random [math]\displaystyle{ x\in S }[/math] as the pivot;
    • partition [math]\displaystyle{ S }[/math] into [math]\displaystyle{ S_1 }[/math], [math]\displaystyle{ \{x\} }[/math], and [math]\displaystyle{ S_2 }[/math], where all numbers in [math]\displaystyle{ S_1 }[/math] are smaller than [math]\displaystyle{ x }[/math] and all numbers in [math]\displaystyle{ S_2 }[/math] are larger than [math]\displaystyle{ x }[/math];
    • recursively sort [math]\displaystyle{ S_1 }[/math] and [math]\displaystyle{ S_2 }[/math];

Analysis of Random Quicksort

Our goal is to analyze the expected number of comparisons during an execution of RandQSort with an arbitrary input [math]\displaystyle{ S }[/math]. We achieve this by measuring the chance that each pair of elements are compared, and summing all of them up due to Linearity of Expectation.

Let [math]\displaystyle{ a_i }[/math] denote the [math]\displaystyle{ i }[/math]th smallest element in [math]\displaystyle{ S }[/math]. Let [math]\displaystyle{ X_{ij}\in\{0,1\} }[/math] be the random variable which indicates whether [math]\displaystyle{ a_i }[/math] and [math]\displaystyle{ a_j }[/math] are compared during the execution of RandQSort. That is:

[math]\displaystyle{ \begin{align} X_{ij} &= \begin{cases} 1 & a_i\mbox{ and }a_j\mbox{ are compared}\\ 0 & \mbox{otherwise} \end{cases}. \end{align} }[/math]

Elements [math]\displaystyle{ a_i }[/math] and [math]\displaystyle{ a_j }[/math] are compared only if one of them is chosen as pivot. After comparison they are separated (thus are never compared again). So we have the following observations:

Observation 1: Every pair of [math]\displaystyle{ a_i }[/math] and [math]\displaystyle{ a_j }[/math] are compared at most once.

Therefore the sum of [math]\displaystyle{ X_{ij} }[/math] for all pair [math]\displaystyle{ \{i, j\} }[/math] gives the total number of comparisons. The expected number of comparisons is [math]\displaystyle{ \mathbf{E}\left[\sum_{i=1}^n\sum_{j\gt i}X_{ij}\right] }[/math]. Due to Linearity of Expectation, [math]\displaystyle{ \mathbf{E}\left[\sum_{i=1}^n\sum_{j\gt i}X_{ij}\right] = \sum_{i=1}^n\sum_{j\gt i}\mathbf{E}\left[X_{ij}\right] }[/math]. Our next step is to analyze [math]\displaystyle{ \mathbf{E}\left[X_{ij}\right] }[/math] for each [math]\displaystyle{ \{i, j\} }[/math].

By the definition of expectation and [math]\displaystyle{ X_{ij} }[/math],

[math]\displaystyle{ \begin{align} \mathbf{E}\left[X_{ij}\right] &= 1\cdot \Pr[a_i\mbox{ and }a_j\mbox{ are compared}] + 0\cdot \Pr[a_i\mbox{ and }a_j\mbox{ are not compared}]\\ &= \Pr[a_i\mbox{ and }a_j\mbox{ are compared}]. \end{align} }[/math]

We are going to bound this probability.

Observation 2: [math]\displaystyle{ a_i }[/math] and [math]\displaystyle{ a_j }[/math] are compared if and only if one of them is chosen as pivot when they are still in the same subset.

This is easy to verify: just check the algorithm. The next one is a bit complicated.

Observation 3: If [math]\displaystyle{ a_i }[/math] and [math]\displaystyle{ a_j }[/math] are still in the same subset then all [math]\displaystyle{ \{a_i, a_{i+1}, \ldots, a_{j-1}, a_{j}\} }[/math] are in the same subset.

We can verify this by induction. Initially, [math]\displaystyle{ S }[/math] itself has the property described above; and partitioning any [math]\displaystyle{ S }[/math] with the property into [math]\displaystyle{ S_1 }[/math] and [math]\displaystyle{ S_2 }[/math] will preserve the property for both [math]\displaystyle{ S_1 }[/math] and [math]\displaystyle{ S_2 }[/math]. Therefore Claim 3 holds.

Combining Observation 2 and 3, we have:

Observation 4: [math]\displaystyle{ a_i }[/math] and [math]\displaystyle{ a_j }[/math] are compared only if one of [math]\displaystyle{ \{a_i, a_j\} }[/math] is chosen from [math]\displaystyle{ \{a_i, a_{i+1}, \ldots, a_{j-1}, a_{j}\} }[/math].

And,

Observation 5: Every one of [math]\displaystyle{ \{a_i, a_{i+1}, \ldots, a_{j-1}, a_{j}\} }[/math] is chosen equal-probably.

This is because the Random Quicksort chooses the pivot uniformly at random.

Observation 4 and 5 together imply:

[math]\displaystyle{ \begin{align} \Pr[a_i\mbox{ and }a_j\mbox{ are compared}] &\le \frac{2}{j-i+1}. \end{align} }[/math]
Remark: Perhaps you feel confused about the above argument. You may ask: "The algorithm chooses pivots for many times during the execution. Why in the above argument, it looks like the pivot is chosen only once?" Good question! Let's see what really happens by looking closely.

For any pair [math]\displaystyle{ a_i }[/math] and [math]\displaystyle{ a_j }[/math], initially [math]\displaystyle{ \{a_i, a_{i+1}, \ldots, a_{j-1}, a_{j}\} }[/math] are all in the same set [math]\displaystyle{ S }[/math] (obviously!). During the execution of the algorithm, the set which containing [math]\displaystyle{ \{a_i, a_{i+1}, \ldots, a_{j-1}, a_{j}\} }[/math] are shrinking (due to the pivoting), until one of [math]\displaystyle{ \{a_i, a_{i+1}, \ldots, a_{j-1}, a_{j}\} }[/math] is chosen, and the set is partitioned into different subsets. We ask for the probability that the chosen one is among [math]\displaystyle{ \{a_i, a_j\} }[/math]. So we really care about "the last" pivoting before [math]\displaystyle{ \{a_i, a_{i+1}, \ldots, a_{j-1}, a_{j}\} }[/math] is split.

Formally, let [math]\displaystyle{ Y }[/math] be the random variable denoting the pivot element. We know that for each [math]\displaystyle{ a_k\in\{a_i, a_{i+1}, \ldots, a_{j-1}, a_{j}\} }[/math], [math]\displaystyle{ Y=a_k }[/math] with the same probability, and [math]\displaystyle{ Y\not\in\{a_i, a_{i+1}, \ldots, a_{j-1}, a_{j}\} }[/math] with an unknown probability (remember that there might be other elements in the same subset with [math]\displaystyle{ \{a_i, a_{i+1}, \ldots, a_{j-1}, a_{j}\} }[/math]). The probability we are looking for is actually [math]\displaystyle{ \Pr[Y\in \{a_i, a_j\}\mid Y\in\{a_i, a_{i+1}, \ldots, a_{j-1}, a_{j}\}] }[/math], which is always [math]\displaystyle{ \frac{2}{j-i+1} }[/math], provided that [math]\displaystyle{ Y }[/math] is uniform over [math]\displaystyle{ \{a_i, a_{i+1}, \ldots, a_{j-1}, a_{j}\} }[/math].

The conditional probability rules out the irrelevant events in a probabilistic argument.

Summing all up:

[math]\displaystyle{ \begin{align} \mathbf{E}\left[\sum_{i=1}^n\sum_{j\gt i}X_{ij}\right] &= \sum_{i=1}^n\sum_{j\gt i}\mathbf{E}\left[X_{ij}\right]\\ &\le \sum_{i=1}^n\sum_{j\gt i}\frac{2}{j-i+1}\\ &= \sum_{i=1}^n\sum_{k=2}^{n-i+1}\frac{2}{k} & & (\mbox{Let }k=j-i+1)\\ &\le \sum_{i=1}^n\sum_{k=1}^{n}\frac{2}{k}\\ &= 2n\sum_{k=1}^{n}\frac{1}{k}\\ &= 2n H(n). \end{align} }[/math]

[math]\displaystyle{ H(n) }[/math] is the [math]\displaystyle{ n }[/math]th Harmonic number. It holds that

[math]\displaystyle{ \begin{align}H(n) = \ln n+O(1)\end{align} }[/math].

Therefore, for an arbitrary input [math]\displaystyle{ S }[/math] of [math]\displaystyle{ n }[/math] numbers, the expected number of comparisons taken by RandQSort to sort [math]\displaystyle{ S }[/math] is [math]\displaystyle{ \mathrm{O}(n\log n) }[/math].


Distributions of Coin Flips

We introduce several important distributions induced by independent coin flips (independent probabilistic experiments), including: Bernoulli trial, geometric distribution, binomial distribution.

Bernoulli trial (Bernoulli distribution)

Bernoulli trial describes the probability distribution of a single (biased) coin flip. Suppose that we flip a (biased) coin where the probability of HEADS is [math]\displaystyle{ p }[/math]. Let [math]\displaystyle{ X }[/math] be the 0-1 random variable which indicates whether the result is HEADS. We say that [math]\displaystyle{ X }[/math] follows the Bernoulli distribution with parameter [math]\displaystyle{ p }[/math]. Formally,

[math]\displaystyle{ \begin{align} X &= \begin{cases} 1 & \text{with probability }p\\ 0 & \text{with probability }1-p \end{cases} \end{align} }[/math].

Geometric distribution

Suppose we flip the same coin repeatedly until HEADS appears, where each coin flip is independent and follows the Bernoulli distribution with parameter [math]\displaystyle{ p }[/math]. Let [math]\displaystyle{ X }[/math] be the random variable denoting the total number of coin flips. Then [math]\displaystyle{ X }[/math] has the geometric distribution with parameter [math]\displaystyle{ p }[/math]. Formally, [math]\displaystyle{ \Pr[X=k]=(1-p)^{k-1}p }[/math].

For geometric [math]\displaystyle{ X }[/math], [math]\displaystyle{ \mathbf{E}[X]=\frac{1}{p} }[/math]. This can be verified by directly computing [math]\displaystyle{ \mathbf{E}[X] }[/math] by the definition of expectations. There is also a smarter way of computing [math]\displaystyle{ \mathbf{E}[X] }[/math], by using indicators and the linearity of expectations. For [math]\displaystyle{ k=0, 1, 2, \ldots }[/math], let [math]\displaystyle{ Y_k }[/math] be the 0-1 random variable such that [math]\displaystyle{ Y_k=1 }[/math] if and only if none of the first [math]\displaystyle{ k }[/math] coin flipings are HEADS, thus [math]\displaystyle{ \mathbf{E}[Y_k]=\Pr[Y_k=1]=(1-p)^{k} }[/math]. A key observation is that [math]\displaystyle{ X=\sum_{k=0}^\infty Y_k }[/math]. Thus, due to the linearity of expectations,

[math]\displaystyle{ \begin{align} \mathbf{E}[X] = \mathbf{E}\left[\sum_{k=0}^\infty Y_k\right] = \sum_{k=0}^\infty \mathbf{E}[Y_k] = \sum_{k=0}^\infty (1-p)^k = \frac{1}{1-(1-p)} =\frac{1}{p}. \end{align} }[/math]

Binomial distribution

Suppose we flip the same (biased) coin for [math]\displaystyle{ n }[/math] times, where each coin flip is independent and follows the Bernoulli distribution with parameter [math]\displaystyle{ p }[/math]. Let [math]\displaystyle{ X }[/math] be the number of HEADS. Then [math]\displaystyle{ X }[/math] has the binomial distribution with parameters [math]\displaystyle{ n }[/math] and [math]\displaystyle{ p }[/math]. Formally, [math]\displaystyle{ \Pr[X=k]={n\choose k}p^k(1-p)^{n-k} }[/math].

A binomial random variable [math]\displaystyle{ X }[/math] with parameters [math]\displaystyle{ n }[/math] and [math]\displaystyle{ p }[/math] is usually denoted by [math]\displaystyle{ B(n,p) }[/math].

As we saw above, by applying the linearity of expectations, it is easy to show that [math]\displaystyle{ \mathbf{E}[X]=np }[/math] for an [math]\displaystyle{ X=B(n,p) }[/math].

Balls into Bins

Birthday Problem

There are [math]\displaystyle{ m }[/math] students in the class. Assume that for each student, his/her birthday is uniformly and independently distributed over the 365 days in a years. We wonder what the probability that no two students share a birthday.

Due to the pigeonhole principle, it is obvious that for [math]\displaystyle{ m\gt 365 }[/math], there must be two students with the same birthday. Surprisingly, for any [math]\displaystyle{ m\gt 57 }[/math] this event occurs with more than 99% probability. This is called the birthday paradox. Despite the name, the birthday paradox is not a real paradox.

We can model this problem as a balls-into-bins problem. [math]\displaystyle{ m }[/math] different balls (students) are uniformly and independently thrown into 365 bins (days). More generally, let [math]\displaystyle{ n }[/math] be the number of bins. We ask for the probability of the following event [math]\displaystyle{ \mathcal{E} }[/math]

  • [math]\displaystyle{ \mathcal{E} }[/math]: there is no bin with more than one balls (i.e. no two students share birthday).

We first analyze this by counting. There are totally [math]\displaystyle{ n^m }[/math] ways of assigning [math]\displaystyle{ m }[/math] balls to [math]\displaystyle{ n }[/math] bins. The number of assignments that no two balls share a bin is [math]\displaystyle{ {n\choose m}m! }[/math].

Thus the probability is given by:

[math]\displaystyle{ \begin{align} \Pr[\mathcal{E}] = \frac{{n\choose m}m!}{n^m}. \end{align} }[/math]

Recall that [math]\displaystyle{ {n\choose m}=\frac{n!}{(n-m)!m!} }[/math]. Then

[math]\displaystyle{ \begin{align} \Pr[\mathcal{E}] = \frac{{n\choose m}m!}{n^m} = \frac{n!}{n^m(n-m)!} = \frac{n}{n}\cdot\frac{n-1}{n}\cdot\frac{n-2}{n}\cdots\frac{n-(m-1)}{n} = \prod_{k=1}^{m-1}\left(1-\frac{k}{n}\right). \end{align} }[/math]

There is also a more "probabilistic" argument for the above equation. To be rigorous, we need the following theorem, which holds generally and is very useful for computing the AND of many events.

By the definition of conditional probability, [math]\displaystyle{ \Pr[A\mid B]=\frac{\Pr[A\wedge B]}{\Pr[B]} }[/math]. Thus, [math]\displaystyle{ \Pr[A\wedge B] =\Pr[B]\cdot\Pr[A\mid B] }[/math]. This hints us that we can compute the probability of the AND of events by conditional probabilities. Formally, we have the following theorem:

Theorem:

Let [math]\displaystyle{ \mathcal{E}_1, \mathcal{E}_2, \ldots, \mathcal{E}_n }[/math] be any [math]\displaystyle{ n }[/math] events. Then
[math]\displaystyle{ \begin{align} \Pr\left[\bigwedge_{i=1}^n\mathcal{E}_i\right] &= \prod_{k=1}^n\Pr\left[\mathcal{E}_k \mid \bigwedge_{i\lt k}\mathcal{E}_i\right]. \end{align} }[/math]

Proof: It holds that [math]\displaystyle{ \Pr[A\wedge B] =\Pr[B]\cdot\Pr[A\mid B] }[/math]. Thus, let [math]\displaystyle{ A=\mathcal{E}_n }[/math] and [math]\displaystyle{ B=\mathcal{E}_1\wedge\mathcal{E}_2\wedge\cdots\wedge\mathcal{E}_{n-1} }[/math], then

[math]\displaystyle{ \begin{align} \Pr[\mathcal{E}_1\wedge\mathcal{E}_2\wedge\cdots\wedge\mathcal{E}_n] &= \Pr[\mathcal{E}_1\wedge\mathcal{E}_2\wedge\cdots\wedge\mathcal{E}_{n-1}]\cdot\Pr\left[\mathcal{E}_n\mid \bigwedge_{i\lt n}\mathcal{E}_i\right]. \end{align} }[/math]

Recursively applying this equation to [math]\displaystyle{ \Pr[\mathcal{E}_1\wedge\mathcal{E}_2\wedge\cdots\wedge\mathcal{E}_{n-1}] }[/math] until there is only [math]\displaystyle{ \mathcal{E}_1 }[/math] left, the theorem is proved. [math]\displaystyle{ \square }[/math]

Now we are back to the probabilistic analysis of the birthday problem, with a general setting of [math]\displaystyle{ m }[/math] students and [math]\displaystyle{ n }[/math] possible birthdays (imagine that we live in a planet where a year has [math]\displaystyle{ n }[/math] days).

The first student has a birthday (of course!). The probability that the second student has a different birthday is [math]\displaystyle{ \left(1-\frac{1}{n}\right) }[/math]. Given that the first two students have different birthdays, the probability that the third student has a different birthday from the first two is [math]\displaystyle{ \left(1-\frac{2}{n}\right) }[/math]. Continuing this on, assuming that the first [math]\displaystyle{ k-1 }[/math] students all have different birthdays, the probability that the [math]\displaystyle{ k }[/math]th student has a different birthday than the first [math]\displaystyle{ k-1 }[/math], is given by [math]\displaystyle{ \left(1-\frac{k-1}{n}\right) }[/math]. So the probability that all [math]\displaystyle{ m }[/math] students have different birthdays is the product of all these conditional probabilities:

[math]\displaystyle{ \begin{align} \Pr[\mathcal{E}]=\left(1-\frac{1}{n}\right)\cdot \left(1-\frac{2}{n}\right)\cdots \left(1-\frac{m-1}{n}\right) &= \prod_{k=1}^{m-1}\left(1-\frac{k}{n}\right), \end{align} }[/math]

which is the same as what we got by the counting argument.

There are several ways of analyzing this formular. Here is a convenient one: Due to Taylor's expansion, [math]\displaystyle{ e^{-k/n}\approx 1-k/n }[/math]. Then

[math]\displaystyle{ \begin{align} \prod_{k=1}^{m-1}\left(1-\frac{k}{n}\right) &\approx \prod_{k=1}^{m-1}e^{-\frac{k}{n}}\\ &= \exp\left(-\sum_{k=1}^{m-1}\frac{k}{n}\right)\\ &= e^{-m(m-1)/2n}\\ &\approx e^{-m^2/2n}. \end{align} }[/math]

The quality of this approximation is shown in the Figure.

Therefore, for [math]\displaystyle{ m=\sqrt{2n\ln \frac{1}{\epsilon}} }[/math], the probability that [math]\displaystyle{ \Pr[\mathcal{E}]\approx\epsilon }[/math].

Coupon Collector

Suppose that a chocolate company releases [math]\displaystyle{ n }[/math] different types of coupons. Each box of chocolates contains one coupon with a uniformly random type. Once you have collected all [math]\displaystyle{ n }[/math] types of coupons, you will get a prize. So how many boxes of chocolates you are expected to buy to win the prize?

The coupon collector problem can be described in the balls-into-bins model as follows. We keep throwing balls one-by-one into [math]\displaystyle{ n }[/math] bins (coupons), such that each ball is thrown into a bin uniformly and independently at random. Each ball corresponds to a box of chocolate, and each bin corresponds to a type of coupon. Thus, the number of boxes bought to collect [math]\displaystyle{ n }[/math] coupons is just the number of balls thrown until none of the [math]\displaystyle{ n }[/math] bins is empty.

Theorem
Let [math]\displaystyle{ X }[/math] be the number of balls thrown uniformly and independently to [math]\displaystyle{ n }[/math] bins until no bin is empty. Then [math]\displaystyle{ \mathbf{E}[X]=nH(n) }[/math], where [math]\displaystyle{ H(n) }[/math] is the [math]\displaystyle{ n }[/math]th harmonic number.
Proof.
Let [math]\displaystyle{ X_i }[/math] be the number of balls thrown while there are exactly [math]\displaystyle{ i-1 }[/math] nonempty bins, then clearly [math]\displaystyle{ X=\sum_{i=1}^n X_i }[/math].

When there are exactly [math]\displaystyle{ i-1 }[/math] nonempty bins, throwing a ball, the probability that the number of nonempty bins increases (i.e. the ball is thrown to an empty bin) is

[math]\displaystyle{ p_i=1-\frac{i-1}{n}. }[/math]

[math]\displaystyle{ X_i }[/math] is the number of balls thrown to make the number of nonempty bins increases from [math]\displaystyle{ i-1 }[/math] to [math]\displaystyle{ i }[/math], i.e. the number of balls thrown until a ball is thrown to a current empty bin. Thus, [math]\displaystyle{ X_i }[/math] follows the geometric distribution, such that

[math]\displaystyle{ \Pr[X_i=k]=(1-p_i)^{k-1}p_i }[/math]

For a geometric random variable, [math]\displaystyle{ \mathbf{E}[X_i]=\frac{1}{p_i}=\frac{n}{n-i+1} }[/math].

Applying the linearity of expectations,

[math]\displaystyle{ \begin{align} \mathbf{E}[X] &= \mathbf{E}\left[\sum_{i=1}^nX_i\right]\\ &= \sum_{i=1}^n\mathbf{E}\left[X_i\right]\\ &= \sum_{i=1}^n\frac{n}{n-i+1}\\ &= n\sum_{i=1}^n\frac{1}{i}\\ &= nH(n), \end{align} }[/math]

where [math]\displaystyle{ H(n) }[/math] is the [math]\displaystyle{ n }[/math]th Harmonic number, and [math]\displaystyle{ H(n)=\ln n+O(1) }[/math]. Thus, for the coupon collectors problem, the expected number of coupons required to obtain all [math]\displaystyle{ n }[/math] types of coupons is [math]\displaystyle{ n\ln n+O(n) }[/math].

[math]\displaystyle{ \square }[/math]

Only knowing the expectation is not good enough. We would like to know how fast the probability decrease as a random variable deviates from its mean value.

Theorem
Let [math]\displaystyle{ X }[/math] be the number of balls thrown uniformly and independently to [math]\displaystyle{ n }[/math] bins until no bin is empty. Then [math]\displaystyle{ \Pr[X\ge n\ln n+cn]\lt e^{-c} }[/math] for any [math]\displaystyle{ c\gt 0 }[/math].
Proof.
For any particular bin [math]\displaystyle{ i }[/math], the probability that bin [math]\displaystyle{ i }[/math] is empty after throwing [math]\displaystyle{ n\ln n+cn }[/math] balls is
[math]\displaystyle{ \left(1-\frac{1}{n}\right)^{n\ln n+cn} \lt e^{-(\ln n+c)} =\frac{1}{ne^c}. }[/math]

By the union bound, the probability that there exists an empty bin after throwing [math]\displaystyle{ n\ln n+cn }[/math] balls is

[math]\displaystyle{ \Pr[X\ge n\ln n+cn] \lt n\cdot \frac{1}{ne^c} =e^{-c}. }[/math]
[math]\displaystyle{ \square }[/math]

Stable Marriage

We now consider the famous stable marriage problem or stable matching problem (SMP). This problem captures two aspects: allocations (matchings) and stability, two central topics in economics.

An instance of stable marriage consists of:

  • [math]\displaystyle{ n }[/math] men and [math]\displaystyle{ n }[/math] women;
  • each person associated with a strictly ordered preference list containing all the members of the opposite sex.

Formally, let [math]\displaystyle{ M }[/math] be the set of [math]\displaystyle{ n }[/math] men and [math]\displaystyle{ W }[/math] be the set of [math]\displaystyle{ n }[/math] women. Each man [math]\displaystyle{ m\in M }[/math] is associated with a permutation [math]\displaystyle{ p_m }[/math] of elemets in [math]\displaystyle{ W }[/math] and each woman [math]\displaystyle{ w\in W }[/math] is associated with a permutation [math]\displaystyle{ p_w }[/math] of elements in [math]\displaystyle{ M }[/math].

A matching is a one-one correspondence [math]\displaystyle{ \phi:M\rightarrow W }[/math]. We said a man [math]\displaystyle{ m }[/math] and a woman [math]\displaystyle{ w }[/math] are partners in [math]\displaystyle{ \phi }[/math] if [math]\displaystyle{ w=\phi(m) }[/math].

Definition (stable matching)
A pair [math]\displaystyle{ (m,w) }[/math] of a man and woman is a blocking pair in a matching [math]\displaystyle{ \phi }[/math] if [math]\displaystyle{ m }[/math] and [math]\displaystyle{ w }[/math] are not partners in [math]\displaystyle{ \phi }[/math] but
  • [math]\displaystyle{ m }[/math] prefers [math]\displaystyle{ w }[/math] to [math]\displaystyle{ \phi(m) }[/math], and
  • [math]\displaystyle{ w }[/math] prefers [math]\displaystyle{ m }[/math] to [math]\displaystyle{ \phi(w) }[/math].
A matching [math]\displaystyle{ \phi }[/math] is stable if there is no blocking pair in it.

It is unclear from the definition itself whether stable matchings always exist, and how to efficiently find a stable matching. Both questions are answered by the following proposal algorithm due to Gale and Shapley.

The proposal algorithm (Gale-Shapley 1962)
Initially, all person are not married;
in each step (called a proposal):
  • an arbitrary unmarried man [math]\displaystyle{ m }[/math] proposes to the woman [math]\displaystyle{ w }[/math] who is ranked highest in his preference list [math]\displaystyle{ p_m }[/math] among all the women who has not yet rejected [math]\displaystyle{ m }[/math];
  • if [math]\displaystyle{ w }[/math] is still single then [math]\displaystyle{ w }[/math] accepts the proposal and is married to [math]\displaystyle{ m }[/math];
  • if [math]\displaystyle{ w }[/math] is married to another man [math]\displaystyle{ m' }[/math] who is ranked lower than [math]\displaystyle{ m }[/math] in her preference list [math]\displaystyle{ p_w }[/math] then [math]\displaystyle{ w }[/math] divorces [math]\displaystyle{ m' }[/math] (thus [math]\displaystyle{ m' }[/math] becomes single again and considers himself as rejected by [math]\displaystyle{ w }[/math]) and is married to [math]\displaystyle{ m }[/math];
  • if otherwise [math]\displaystyle{ w }[/math] rejects [math]\displaystyle{ m }[/math];

The algorithm terminates when the last single woman receives a proposal. Since for every pair [math]\displaystyle{ (m,w)\in M\times W }[/math] of man and woman, [math]\displaystyle{ m }[/math] proposes to [math]\displaystyle{ w }[/math] at most once. The algorithm terminates in at most [math]\displaystyle{ n^2 }[/math] proposals in the worst case.

It is obvious to see that the algorithm retruns a macthing, and this matching must be stable. To see this, by contradiction suppose that the algorithm resturns a macthing [math]\displaystyle{ \phi }[/math], such that two men [math]\displaystyle{ A, B }[/math] are macthed to two women [math]\displaystyle{ a,b }[/math] in [math]\displaystyle{ \phi }[/math] respectively, but [math]\displaystyle{ A }[/math] and [math]\displaystyle{ b }[/math] prefers each other to their partners [math]\displaystyle{ a }[/math] and [math]\displaystyle{ B }[/math] respectively. By definition of the algorithm, [math]\displaystyle{ A }[/math] would have proposed to [math]\displaystyle{ b }[/math] before proposing to [math]\displaystyle{ a }[/math], by which time [math]\displaystyle{ b }[/math] must either be single or be matched to a man ranked lower than [math]\displaystyle{ A }[/math] in her list (because her final partner [math]\displaystyle{ B }[/math] is ranked lower than [math]\displaystyle{ A }[/math]), which means [math]\displaystyle{ b }[/math] must have accepted [math]\displaystyle{ A }[/math]'s proposal, a contradiction.


We are interested in the average-case performance of this algorithm, that is, the expected number of proposals if everyone's preference list is a uniformly and independently random permutation.

The following principle of deferred decisions is quite useful in analysing performance of algorithm with random input.

Principle of deferred decisions
The decision of random choice in the random input can be deferred to the running time of the algorithm.

Apply the principle of deferred decisions, the deterministic proposal algorithm with random permutations as input is equivalent to the following random process:

  • At each step, a man [math]\displaystyle{ m }[/math] choose a woman [math]\displaystyle{ w }[/math] uniformly and independently at random to propose, among all the women who have not rejected him yet. (sample without replacement)

We then compare the above process with the following modified process:

  • The man [math]\displaystyle{ m }[/math] repeatedly samples a uniform and independent woman to propose among all women, until he successfully samples a woman who has not rejected him and propose to her. (sample with replacement)

It is easy to see that the modified process (sample with replacement) is no more efficient than the original process (sample without replacement) because it simulates the original process if at each step we only count the last proposal to the woman who has not rejected the man. Such comparison of two random processes by forcing them to be related in some way is called coupling.

Note that in the modified process (sample with replacement), each proposal, no matter from which man, is going to a uniformly and independently random women. And we know that the algorithm terminated once the last single woman receives a proposal, i.e. once all [math]\displaystyle{ n }[/math] women have received at least one proposal. This is the coupon collector problem with proposals as balls (cookie boxes) and women as bins (coupons). Due to our analysis of the coupon collector problem, the expected number of proposals is bounded by [math]\displaystyle{ O(n\ln n) }[/math].

Occupancy Problem

Now we ask about the loads of bins. Assuming that [math]\displaystyle{ m }[/math] balls are uniformly and independently assigned to [math]\displaystyle{ n }[/math] bins, for [math]\displaystyle{ 1\le i\le n }[/math], let [math]\displaystyle{ X_i }[/math] be the load of the [math]\displaystyle{ i }[/math]th bin, i.e. the number of balls in the [math]\displaystyle{ i }[/math]th bin.

An easy analysis shows that for every bin [math]\displaystyle{ i }[/math], the expected load [math]\displaystyle{ \mathbf{E}[X_i] }[/math] is equal to the average load [math]\displaystyle{ m/n }[/math].

Because there are totally [math]\displaystyle{ m }[/math] balls, it is always true that [math]\displaystyle{ \sum_{i=1}^n X_i=m }[/math].

Therefore, due to the linearity of expectations,

[math]\displaystyle{ \begin{align} \sum_{i=1}^n\mathbf{E}[X_i] &= \mathbf{E}\left[\sum_{i=1}^n X_i\right] = \mathbf{E}\left[m\right] =m. \end{align} }[/math]

Because for each ball, the bin to which the ball is assigned is uniformly and independently chosen, the distributions of the loads of bins are identical. Thus [math]\displaystyle{ \mathbf{E}[X_i] }[/math] is the same for each [math]\displaystyle{ i }[/math]. Combining with the above equation, it holds that for every [math]\displaystyle{ 1\le i\le m }[/math], [math]\displaystyle{ \mathbf{E}[X_i]=\frac{m}{n} }[/math]. So the average is indeed the average!


Next we analyze the distribution of the maximum load. We show that when [math]\displaystyle{ m=n }[/math], i.e. [math]\displaystyle{ n }[/math] balls are uniformly and independently thrown into [math]\displaystyle{ n }[/math] bins, the maximum load is [math]\displaystyle{ O\left(\frac{\log n}{\log\log n}\right) }[/math] with high probability.

Theorem
Suppose that [math]\displaystyle{ n }[/math] balls are thrown independently and uniformly at random into [math]\displaystyle{ n }[/math] bins. For [math]\displaystyle{ 1\le i\le n }[/math], let [math]\displaystyle{ X_i }[/math] be the random variable denoting the number of balls in the [math]\displaystyle{ i }[/math]th bin. Then
[math]\displaystyle{ \Pr\left[\max_{1\le i\le n}X_i \ge\frac{3\ln n}{\ln\ln n}\right] \lt \frac{1}{n}. }[/math]
Proof.
Let [math]\displaystyle{ M }[/math] be an integer. Take bin 1. For any particular [math]\displaystyle{ M }[/math] balls, these [math]\displaystyle{ M }[/math] balls are all thrown to bin 1 with probability [math]\displaystyle{ (1/n)^M }[/math], and there are totally [math]\displaystyle{ {n\choose M} }[/math] distinct sets of [math]\displaystyle{ M }[/math] balls. Therefore, applying the union bound,
[math]\displaystyle{ \begin{align}\Pr\left[X_1\ge M\right] &\le {n\choose M}\left(\frac{1}{n}\right)^M\\ &= \frac{n!}{M!(n-M)!n^M}\\ &= \frac{1}{M!}\cdot\frac{n(n-1)(n-2)\cdots(n-M+1)}{n^M}\\ &= \frac{1}{M!}\cdot \prod_{i=0}^{M-1}\left(1-\frac{i}{n}\right)\\ &\le \frac{1}{M!}. \end{align} }[/math]

According to Stirling's approximation, [math]\displaystyle{ M!\approx \sqrt{2\pi M}\left(\frac{M}{e}\right)^M }[/math], thus

[math]\displaystyle{ \frac{1}{M!}\le\left(\frac{e}{M}\right)^M. }[/math]
Figure 1

Due to the symmetry. All [math]\displaystyle{ X_i }[/math] have the same distribution. Apply the union bound again,

[math]\displaystyle{ \begin{align} \Pr\left[\max_{1\le i\le n}X_i\ge M\right] &= \Pr\left[(X_1\ge M) \vee (X_2\ge M) \vee\cdots\vee (X_n\ge M)\right]\\ &\le n\Pr[X_1\ge M]\\ &\le n\left(\frac{e}{M}\right)^M. \end{align} }[/math]

When [math]\displaystyle{ M=3\ln n/\ln\ln n }[/math],

[math]\displaystyle{ \begin{align} \left(\frac{e}{M}\right)^M &= \left(\frac{e\ln\ln n}{3\ln n}\right)^{3\ln n/\ln\ln n}\\ &\lt \left(\frac{\ln\ln n}{\ln n}\right)^{3\ln n/\ln\ln n}\\ &= e^{3(\ln\ln\ln n-\ln\ln n)\ln n/\ln\ln n}\\ &= e^{-3\ln n+3\ln\ln\ln n\ln n/\ln\ln n}\\ &\le e^{-2\ln n}\\ &= \frac{1}{n^2}. \end{align} }[/math]

Therefore,

[math]\displaystyle{ \begin{align} \Pr\left[\max_{1\le i\le n}X_i\ge \frac{3\ln n}{\ln\ln n}\right] &\le n\left(\frac{e}{M}\right)^M &\lt \frac{1}{n}. \end{align} }[/math]
[math]\displaystyle{ \square }[/math]

When [math]\displaystyle{ m\gt n }[/math], Figure 1 illustrates the results of several random experiments, which show that the distribution of the loads of bins becomes more even as the number of balls grows larger than the number of bins.

Formally, it can be proved that for [math]\displaystyle{ m=\Omega(n\log n) }[/math], with high probability, the maximum load is within [math]\displaystyle{ O\left(\frac{m}{n}\right) }[/math], which is asymptotically equal to the average load.