An introduction to Randomized Algorithms
To follow along, you should know:
- Basic discrete probability (independence, inclusion-exclusion, conditional probability, and the law of total probability).
- What a matrix is and how to multiply matrices.
- What a graph is. The basic definition is sufficient.
- Big-O notation is helpful, but optional.
AI Disclosure: I conceived and directed all the visuals in this explainer, which were implemented in JavaScript with the help of a coding agent. The source code for the entire project can be found in the public GitHub repository for reference.
Annotation note: Any underlined text in this explainer is an annotation. Hover your cursor over this example to see additional information. You can also click it.
1 Introduction
Imagine you’re a cracked engineer living in the 80s. You have a Compaq Deskpro 386/20, a state of the art personal computer of the time, that boasts an impressive performance 160,000 Floating point operations (FLOPs) per second.
You’re working on a computational math problem that requires you to compute matrix multiplications of matrices with massive dimensions.
In particular, say you have two 3000 \times 3000 dimensional matrices and you want to compute AB.
To do this, you decide to implement the pen and paper matrix multiplication algorithm in code, and then run the code on your machine. How many Floating point operations will this take?
Let P \in \mathbb{R}^{a \times b}, Q \in \mathbb{R}^{b \times c}, C = PQ \in \mathbb{R}^{a \times c}. C has a \times c entries. To compute the (i,j)^{\text{th}} entry, denoted by C_{i,j}, we do a dot product between the ith row of P and the jth column of Q, that is:
\begin{aligned} C_{i,j} &= \sum_{k=1}^b P_{i,k}Q_{k,j} \\ &= P_{i,1}Q_{1, j} + P_{i, 2}Q_{2, j} + \dots + P_{i, b}Q_{b, j}. \end{aligned}
There are b multiplications, and b-1 additions, giving us a total FLOP count of 2b-1 to compute one entry of C. Thus the total FLOP count to compute C is given as
a \times c \times (2b-1). \tag{1}
In your case you’re multiplying two square d \times d matrices. The total FLOPs in that case:
\begin{aligned} \text{Total FLOPs} &= d \times d \times (2d-1) \\ &= 2d^3-d^2 \\ &\in O(d^3). \end{aligned}
In your case, you want to compute AB.
\begin{aligned} d &= 3000 \\ \text{Total FLOPs} &= 2\times 3000^3 - 3000^2 \\ &= \boxed{53{,}991{,}000{,}000\text{ FLOPs}}. \end{aligned}
The time your computer would take to do this is given by:
\begin{aligned} \text{Time taken by Compaq 386/20} &= \frac{\text{Total FLOPs}}{\text{FLOPs per second}} \\ &= \frac{53{,}991{,}000{,}000}{160{,}000} \\ &= 337{,}443.75\,\text{seconds} \\ &\approx 94\,\text{hours} \\ &\approx 3\,\text{days }22\,\text{hours}. \end{aligned}
Multiplying just one pair of 3000 \times 3000 matrices would take close to four days!
Realizing the impracticality, you go to a supercomputer provider and pay them to compute the answer for you. They give you the answer, the matrix, C.
Now you’re a bit paranoid. You did not write the code for the supercomputer. Maybe their code has a subtle bug that got exposed on your massive matrix multiplication load. How can you be sure AB is indeed C?
This is what we call the matrix multiplication verification problem (MMV).
Here it is stated more formally:
Matrix Multiplication Verification (MMV)
Given three real-valued matrices A, B, C \in \mathbb{R}^{d \times d}, verify whether AB = C.
How can you solve this problem? The only way to be absolutely sure is to compute AB yourself and compare it with C. But the impracticality of doing so was the whole point of using the supercomputer provider. What do you do then? How can you be more confident?
Here’s an algorithm idea -
Freivalds’ Algorithm
- Sample a vector uniformly at random.
- Compute A(Bv) and Cv.
- If A(Bv)=Cv, we say that AB=C, else we say AB\neq C.
Freivalds’ algorithm is an example of what we call a randomized algorithm. The first step of the algorithm involves a random choice where we randomly choose a vector from the set \{0,1\}^d. Let’s state the definition precisely:
Randomized Algorithm
A randomized algorithm is an algorithm that uses one or more random choices as part of its execution.
When I first encountered Freivalds’ in my studies, I was in disbelief. How can multiplying the matrices on a random vector give us any sort of confidence on whether the matrices in question are equal? This couldn’t possibly work. It’s too simple for a problem that seems more complicated. But it somehow does work.
This is a common theme in randomized algorithms, seemingly complicated problems have absurdly simple solutions. This feeling of disbelief and wonder is precisely why randomized algorithms are my favourite class of algorithms.
Let’s unpack this.
For starters, if AB = C, then for any v \in \{0,1\}^d ABv = Cv. So in this case the algorithm always works correctly.
The number of FLOPs taken to carry out one iteration of the algorithm is also better than computing AB directly:
\begin{aligned} \text{Total FLOPs} &= \text{FLOPs for computing } Bv \\ &\quad + \text{FLOPs for computing } A(Bv) \\ &\quad + \text{FLOPs for computing } Cv. \end{aligned}
Since all the products are between d \times d matrices (A, B, C) and d \times 1 matrices (Bv, v), they take the same FLOPs given by:
d \times 1 \times (2d - 1) = 2d^2 - d.
So, total FLOPs are given as:
\begin{aligned} \text{Total FLOPs} &= 3(2d^2-d) \\ &= 6d^2-3d \\ &\in O(d^2). \end{aligned}
Asymptotically this algorithm is faster than computing AB directly by a factor of d (O(d^3) vs O(d^2)). More exactly,
6d^2 - 3d < 2d^3 - d^2 \qquad \text{for all integers } d \ge 4.
In your case,
\begin{aligned} d &= 3000 \\ \text{Total FLOPs} &= 6 \times 3000^2 - 3 \times 3000 \\ &= \boxed{53{,}991{,}000\text{ FLOPs}}. \end{aligned}
The time your computer would take to carry out one iteration is given by:
\begin{aligned} \text{Time taken by Compaq 386/20} &= \frac{\text{Total FLOPs}}{\text{FLOPs per second}} \\ &= \frac{53{,}991{,}000}{160{,}000} \\ &= 337.44375\,\text{seconds} \\ &\approx 5.62\,\text{minutes}. \end{aligned}
Compared to computing AB yourself, which would take 4 days, you can now run this algorithm to verify AB = C in a measly 5.62 minutes! This is amazing! It’s almost too good to be true. Unfortunately, it is indeed too good to be true.
If AB \neq C, we can come up with a vector v such that ABv = Cv. Here’s an example of a Freivalds’ Algorithm iteration that fails:
The illustrative Freivalds iteration is ready to play.
\begin{aligned} A &= \begin{pmatrix} 1 & 1 \\ 0 & 1 \end{pmatrix}, \\ B &= \begin{pmatrix} 1 & 0 \\ 1 & 1 \end{pmatrix}, \\ C &= \begin{pmatrix} 1 & 2 \\ 1 & 1 \end{pmatrix}. \end{aligned}
\begin{aligned} AB &= \begin{pmatrix} 2 & 1 \\ 1 & 1 \end{pmatrix} \neq C, \\ v^{(1)} &= \begin{pmatrix} 1 \\ 1 \end{pmatrix}. \end{aligned}
\begin{aligned} Bv^{(1)} &= \begin{pmatrix} 1 \\ 2 \end{pmatrix}, \\ A\bigl(Bv^{(1)}\bigr) &= \begin{pmatrix} 3 \\ 2 \end{pmatrix}, \\ Cv^{(1)} &= \begin{pmatrix} 3 \\ 2 \end{pmatrix}. \end{aligned}
ABv^{(1)}=Cv^{(1)}.
So what do we do now?
The key is to remember that we want to be more confident about whether AB = C. Being absolutely sure would be the best case, but since it is impractical to compute AB ourselves, we can settle for the next best thing: be extremely confident, accepting a small chance that we’re wrong.
This is our cue to think in probabilities. What is the probability that this algorithm gives us the wrong result? If we can somehow compute this number and it is low enough, we can be happy. Let’s analyse this:
Let D = AB-C. Thus D=0 \Leftrightarrow AB=C.
The algorithm can fail only when AB \neq C, equivalently D \neq 0. We fix the inputs A, B, C to the algorithm such that AB \neq C. All probabilities below then are over the randomly sampled vector v. The algorithm fails if it claims that AB=C. Let F = \{v \in \{0,1\}^d | Dv = 0\}. The probability of failure then is exactly P[F] = P[Dv = 0].
Let r_k be the k^{\text{th}} entry of the d \times 1 matrix Dv.
\begin{aligned} P[F] &= P[Dv=0] \\ &= P[r_1=0 \cap r_2=0 \cap \cdots \cap r_d=0] \\ &\leq P[r_k=0], \quad \text{for any } k \in \{1,2,\ldots,d\} \end{aligned} \tag{2}
Since D \neq 0, there is at least one non-zero entry in D. Say the (i,j)^{\text{th}} entry, denoted by D_{i,j}, is non-zero.
By definition of matrix multiplication:
r_i = \sum_{k=1}^{d} D_{i,k}v_k = D_{i,j}v_j+Y, \quad \text{where } Y := \sum_{\substack{k=1 \\ k\neq j}}^{d} D_{i,k}v_k \text{ is a random variable}
Let’s compute the probability that r_i=0.
For every sampled vector v, exactly one of Y=0 and Y\neq0 occurs. Therefore, the law of total probability gives:
\begin{aligned} P[r_i=0] &= P[(Y \neq 0 \cap r_i=0) \cup (Y=0 \cap r_i=0)] \\ &= P[r_i=0 \mid Y \neq 0]P[Y \neq 0] + P[r_i=0 \mid Y=0]P[Y=0] \end{aligned} \tag{3}
P[r_i=0 \mid Y \neq 0] = P\left[v_j=-\frac{Y}{D_{i,j}} \mid Y \neq 0\right]
If -\frac{Y}{D_{i,j}} \in \{0,1\} then P[v_j = -\frac{Y}{D_{i,j}} \mid Y \neq 0] = \frac{1}{2}. If -\frac{Y}{D_{i,j}} \in \mathbb{R} \setminus \{0,1\}, P[v_j = -\frac{Y}{D_{i,j}} \mid Y \neq 0] = 0. Thus we get: P[r_i=0 \mid Y \neq 0] = P\left[v_j=-\frac{Y}{D_{i,j}} \mid Y \neq 0\right] \leq \frac{1}{2} \tag{4}
P[r_i=0 \mid Y=0] = P[v_j=0] = \frac{1}{2} \tag{5}
Substituting Equation 4 and Equation 5 into Equation 3, we get:
\begin{aligned} P[r_i=0] &\leq \frac{1}{2}P[Y \neq 0] + \frac{1}{2}P[Y=0] \\ &= \frac{1}{2}\left[P[Y \neq 0] + P[Y=0]\right] \\ &= \frac{1}{2} \end{aligned}
From Equation 2:
P[F] \leq P[r_i=0] \leq \frac{1}{2}
Thus
\boxed{P[F] \leq \frac{1}{2}}
So the probability that one iteration of Freivalds’ gives us the wrong answer is at most \frac{1}{2}. That is not that good.
Hence we come to the sledgehammer of Randomized Algorithms: amplification. Sure, one iteration of Freivalds’ is not that reliable. But we can repeat the algorithm as many times as we like. Let’s make this idea formal:
Amplification
Amplification is a technique for reducing the error probability of a randomized algorithm by running it independently multiple times.
If AB \neq C then all it takes for us to know that is for one of the iterations of Freivalds’ to sample a vector v such that A(Bv) \neq Cv. Each iteration is independent of any other iteration: what vector v was sampled in iteration i has no effect on what vector is sampled in any other iteration j.
Let F_i be the event that the i^{\text{th}} iteration fails, and let F^{(k)} be the event that the complete k-iteration algorithm fails. If we repeat the algorithm k times, we fail only if every iteration fails. The probability of failure becomes:
\begin{aligned} P[F^{(k)}] &= P[F_1 \cap F_2 \cap \dots \cap F_k] \\ &= P[F_1] \times P[F_2] \times \dots \times P[F_k] \\ &\leq \left(\frac{1}{2}\right)^k \\ &= \frac{1}{2^k} \end{aligned}
We can thus get the probability of failure arbitrarily small. If you want the probability of failure to be at most \delta, where \frac{1}{2} \geq \delta > 0, we can repeat the algorithm k = \left\lceil -\log_2(\delta) \right\rceil times.
If you wanted to get the probability of failure to be at most \frac{1}{2^{10}} \approx \frac{1}{1000} = 0.001, you need to repeat the algorithm k = 10 times, after which you can be much more confident about your verification, having at most a 0.1\% chance of being wrong!
That is only about 5.62\text{ minutes} \times 10 = 56.2\text{ minutes} of work on your Compaq Deskpro 386/20, still a measly amount compared to the 4 days that it’d take you to compute AB!
This completes our first randomized algorithm. Isn’t it exciting? This problem defines well what a particular class of randomized algorithms do:
- Find a good random choice to make, and make an algorithm around it.
- Analyse the randomized algorithm and hope that the probability bounds of failure are good enough.
- Amplify: Do several independent iterations of the randomized algorithm to get the failure probability to an arbitrarily small non-zero value.
What do I mean by ‘good enough’? We’ll talk about that soon.
Now we move to another application of randomized algorithms - finding the minimum cut of a graph. This is my favourite graph algorithm, so lock in.
2 Finding the Minimum Cut
You’re a computer network administrator and you’ve been tasked by your company to check its computer network for risks that may make it unavailable to users.
Here’s a visual of the computer network you’re managing:
The lines between each computer represent wires connecting them. In particular, any two computers A and B in the network are connected and reachable to each other if there is a path of wires that leads us from A to B.
One metric you care about is the resistance of the network’s connectivity to damaged wires. In particular, what’s the smallest number of wires that, if damaged, would make the network unconnected (meaning that there are at least two computers that cannot reach each other)?
If we think of the computer network as a graph, with the computers being vertices and the wires being edges, this metric amounts to finding the minimum cut of the graph.
Let’s define what a cut is:
Cut
A cut C = (S, T) of a graph G=(V, E) is a partition of the vertices into two disjoint non-empty subsets S \subset V and T = V \setminus S.
Size of a Cut
The size of a cut C=(S,T) is the number of edges with one endpoint in S and the other in T. We say that any such edge belongs to the cut C, denoted by e\in C.
Let’s see different cuts and their cut sizes on the computer network graph.
- Vertex in S
- Vertex in T
- Edge in cut
- Cut boundary
Now that we understand what a cut is, we can formally state the minimum-cut problem:
Minimum Cut
Given an undirected connected graph G, find a minimum cut of G.
As you may observe, the minimum cut on the computer network graph has size 2, and is defined by the sets S = \{1, 2, 3, 4\} and T = \{5\}. Removing the edges of this cut make the graph disconnected, and form the smallest possible set of edges whose removal disconnects the graph. This is precisely the metric you wanted to calculate for the computer network!
This graph has a unique minimum cut, but a graph can have multiple distinct cuts that attain the minimum-cut size.
There are several deterministic algorithms that can find the minimum cut of a graph. If you’ve ever read them (for example the ones based on finding the max-flow), you’d know those algorithms tend to be complicated. Randomization, on the other hand, gives us a really simple and elegant algorithm.
Before I share this algorithm, we need to define what it means to contract an edge. Here’s the formal definition:
Edge Contraction
Contracting an edge e=\{u,v\} in a graph G means merging u and v into a single merged vertex w. Every edge that previously had exactly one of u or v as an endpoint now has w as that endpoint. All edges between u and v disappear.
Edge 1–3 is selected.
A common confusion is what happens if e = \{u, v\} is one of several parallel edges between u and v. As per the definition, any edge that has both u and v as an endpoint will vanish. Thus all the parallel edges disappear.
Now we have all the machinery needed to state our randomized algorithm.
Karger’s Algorithm
- Repeat until G has exactly two vertices:
- Select an edge in G uniformly at random.
- Contract the selected edge.
- Return the unique cut that remains in G.
Here are two different runs of Karger’s algorithm on the computer network graph:
The two Karger executions are ready to play.
2.1 Bounding the Failure Probability
Just like how we upper bounded the failure probability for Freivalds’ algorithm, let’s do the same for Karger. Before we do that, let’s establish some notation.
Let |V| = n and |E| = m.
Let G^{(i)} denote the graph after i-1 contractions. Thus, G^{(1)}=G.
Let
n_i := |V(G^{(i)})|, \qquad m_i := |E(G^{(i)})|
denote the number of vertices and edges in G^{(i)}, respectively. In particular, n_1=n and m_1=m.
It is important to observe three things:
Every edge contraction decreases the number of vertices in the graph by one. Since each iteration of Karger’s algorithm performs exactly one contraction, we have:
n_i=n-i+1.
Moreover, since the algorithm runs until only two vertices remain, it performs exactly n-2 edge contractions.
Any cut in G^{(i)} corresponds to a cut of the same size in G. Thus, the minimum-cut size of G^{(i)} is greater than or equal to the minimum-cut size of G, for all i.
What do I mean? Let’s make a cut (S^{(i)},T^{(i)}) of G^{(i)}. Every vertex v \in V(G^{(i)}) is a nonempty subset of the original vertices. This includes singleton subsets for unmerged vertices, such as \{5\}, and subsets of the form v=\{a,b,c,\dots\} for merged vertices.
Let C = (S, T) be the cut in G that corresponds to the cut (S^{(i)}, T^{(i)}) in graph G^{(i)}. S is constructed from S^{(i)} by taking the union over every vertex in S^{(i)}.
T is constructed from T^{(i)} in the same way.
The edges belonging to the cut (S,T) are precisely the same as the edges belonging to the cut (S^{(i)},T^{(i)}). Hence, the two cuts have the same size.
As an example, let’s look at graph G^{(4)} in the successful Karger run shown above. The cut it defines is S^{(4)}=\{\{1,2,3,4\}\}, T^{(4)}=\{\{5\}\}.
In the original graph, this is the cut S=\{1,2,3,4\}, T=\{5\}.
There were two edges between \{1,2,3,4\} and \{5\}, which correspond exactly to the edges \{3,5\} and \{4,5\} in G that belong to the cut (S,T).
You can see how this cut looks like in both the graphs in the successful run version of the visual above.
For any minimum cut C in graph G, Karger’s algorithm fails to find C if it contracts any edge e \in C. For example, in the failed Karger run shown in the visual, the algorithm contracted the edge \{3,5\} in graph G^{(2)}, which was part of the minimum cut C.
Let C_1, C_2, \dots, C_\ell be all the minimum cuts of graph G=(V,E), each having size \lambda. Karger’s algorithm fails to find minimum cut C_j precisely when it contracts an edge of C_j. Thus, Karger’s algorithm fails to find any minimum cut when it contracts at least one edge belonging to each minimum cut.
Let E_j be the event that Karger’s algorithm fails to find cut C_j. Let F be the event that Karger’s algorithm fails to find any minimum cut. We have:
P[F] = P[E_1 \cap E_2 \cap \dots \cap E_\ell] \leq P[E_j] \quad \text{for any } j \in \{1,2,\dots,\ell\}.
Let’s fix our focus on the j^{\text{th}} minimum cut and find P[E_j].
To do so, we use the following lemma:
Lemma 1 (Minimum-Cut Edge Bound) If an undirected graph G=(V,E) has n vertices, m edges, and a minimum cut of size \lambda, then
m \geq \frac{\lambda n}{2}.
Proof. We proceed by contradiction. Suppose that
m < \frac{\lambda n}{2}.
By the handshake lemma,
2m=\sum_{v\in V}\deg(v).
Therefore,
\sum_{v\in V}\deg(v)<\lambda n.
This implies that there is at least one vertex v\in V such that \deg(v)<\lambda. Consider the cut S=\{v\} and T=V\setminus\{v\}. This cut has size \deg(v)<\lambda, contradicting the fact that G has a minimum cut of size \lambda.
\blacksquare
We find P[E_j] by computing its complement, P[E_j^c]. That is, we compute the probability that Karger’s algorithm finds the cut C_j.
To find the minimum cut, Karger’s algorithm must not contract any edge of C_j. Let A_i be the event that Karger’s algorithm does not contract an edge of C_j on its i^{\text{th}} iteration.
Thus,
\begin{aligned} P[E_j^c] &= P[A_1 \cap A_2 \cap \dots \cap A_{n-2}] \\ &= P[A_1] \times P[A_2 \mid A_1] \times P[A_3 \mid A_2 \cap A_1] \times \dots \times P[A_{n-2} \mid A_{n-3} \cap A_{n-4} \cap \dots \cap A_1] \end{aligned}
Let’s compute P[A_i \mid A_{i-1} \cap A_{i-2} \cap \dots \cap A_1], the probability that Karger’s algorithm does not contract an edge of C_j on the i^{\text{th}} iteration, given that no edge of C_j was contracted in the previous i-1 iterations.
\begin{aligned} P[A_i \mid A_{i-1} \cap A_{i-2} \cap \dots \cap A_1] &= 1 - P[\text{an edge of } C_j \text{ is contracted on iteration } i \mid A_{i-1} \cap A_{i-2} \cap \dots \cap A_1] \\ &= 1 - \frac{\lambda}{m_i} \end{aligned}
Since all the edges of C_j have survived so far, C_j remains a cut of size \lambda in G^{(i)}. By observation 2, the minimum-cut size of G^{(i)} is at least \lambda. Since C_j is still a cut of size \lambda in G^{(i)}, its minimum-cut size is also at most \lambda. Therefore, the minimum-cut size of G^{(i)} is exactly \lambda.
Applying the lemma to G^{(i)} gives:
m_i \geq \frac{\lambda n_i}{2} = \frac{\lambda(n-i+1)}{2}.
Thus,
\begin{aligned} P[A_i \mid A_{i-1} \cap A_{i-2} \cap \dots \cap A_1] &\geq 1 - \frac{\lambda}{\frac{\lambda(n-i+1)}{2}} \\ &= 1 - \frac{2}{n-i+1} \end{aligned}
Thus we can now compute P[E_j^c]:
\begin{aligned} P[E_j^c] & \geq \prod_{i=1}^{n-2} (1 - \frac{2}{n-i+1}) \\ & = \frac{2}{n(n-1)} \end{aligned}
Thus the probability that Karger’s algorithm fails to find minimum cut C_j is:
\begin{aligned} P[E_j] &= 1 - P[E_j^c] \\ &\leq 1 - \frac{2}{n(n-1)} \end{aligned}
We use the fact that 1 + x \leq e^x for all x \in \mathbb{R}, and get:
\begin{aligned} P[E_j] &\leq 1 - \frac{2}{n(n-1)} \\ &\leq e^{-2/(n(n-1))} \end{aligned}
The probability that Karger’s algorithm fails is then bounded as:
\begin{aligned} P[F] &\leq P[E_j] \\ &\leq e^{-2/(n(n-1))} \end{aligned}
What does this exponential upper bound mean? For our graph, n=5, so the probability of failure is at most e^{-2/(5\cdot4)}=e^{-0.1}\approx 0.905.
That absolutely sucks!
What can we do to get this to a more respectable number? Amplify. We run Karger’s algorithm independently k times and return the smallest cut found. The amplified algorithm fails only if none of the k runs finds a minimum cut.
The probability that all k runs fail is:
\begin{aligned} P[\text{all }k\text{ runs fail}] &= P[F]^k \\ &\leq \left(e^{-2/(n(n-1))}\right)^k \\ &=e^{-2k/(n(n-1))}. \end{aligned}
To make this failure probability at most \delta, we require:
e^{-2k/(n(n-1))}\leq\delta.
Taking logs on both sides and solving for k, we get:
k\geq\frac{n(n-1)}{2}\ln\frac{1}{\delta}.
Thus, it suffices to choose:
k=\left\lceil\frac{n(n-1)}{2}\ln\frac{1}{\delta}\right\rceil.
To get to a 1\% failure bound, given by \delta=0.01, for our graph with n=5, we will need:
k=\left\lceil\frac{5(5-1)}{2}\ln\frac{1}{0.01}\right\rceil=47 \text{ independent runs}.
2.2 Running Time
A single run of Karger’s algorithm performs exactly n-2 edge contractions. Unfortunately, an edge contraction is not a trivial constant time operation. With a suitable graph representation, contraction can be implemented in O(n) time. Therefore, one run of Karger’s algorithm takes O(n^2) time.
Running the algorithm independently k times therefore takes O(kn^2) time.
So far, both algorithms that we’ve discussed have a probability of failing. Can there be randomized algorithms that always run correctly? Yes! The next problem describes one such algorithm.
3 Fast sorting
Sorting a set of numbers is among the most fundamental problems in computer science.
In sorting you have a list of numbers and you want to arrange them in increasing or decreasing order.
Assume you get an input list of numbers A[1 \dots n] of size n.
Here’s an algorithm -
- If A contains at most one element, return A.
- Choose an element of A as the pivot p.
- Partition the remaining elements of A into two lists: L, containing every element less than or equal to p, and R, containing every element greater than p.
- Recursively sort L and R using the same steps 1-3.
- Return the concatenation of the sorted list L, the pivot p, and the sorted list R.
As you can see, if we could somehow find the median as the pivot, the algorithm works well. In particular, if we can find the pivot in O(n) time, where n is the size of the list, the time T(n) taken to sort a list of size n is:
T(n) = 2T(\frac{n}{2}) + cn
Using the Master theorem (if you do not know what that is, just trust me on this):
T(n) = O(n \log n)
No comparison-based sorting algorithm can do asymptotically better than n \log n in general, so this running time is asymptotically optimal.
The only problem is, how are we planning on always finding the median element as the pivot?
Here’s the kicker - you don’t have to worry about that to still get the optimal O(n \log n) time bound on average.
Consider this algorithm:
Randomized Quick Sort
- If A contains at most one element, return A.
- Choose an element of A uniformly at random as the pivot p.
- Partition the remaining elements of A into two lists: L, containing every element less than or equal to p, and R, containing every element greater than p.
- Recursively sort L and R using the same steps 1-3.
- Return the concatenation of the sorted list L, the pivot p, and the sorted list R.
Hence, we have a randomized sorting algorithm, which we call Randomized Quick Sort. Here’s the theorem:
Theorem 1 If Randomized Quick Sort is run repeatedly on the same input of n distinct elements, using fresh random pivot choices each time, its average running time is O(n \log n).
I’m going to omit the proof of this theorem here. Quicksort is widely studied, and proofs are readily available; see Wikipedia’s average-case analysis of Quicksort or Chapter 7 of CLRS.
It’s important to note that this algorithm, unlike the previous two, never fails. Given any array, it’ll always give us the sorted array.
This gives us a nice separation of randomized algorithms into two types: Las Vegas and Monte Carlo.
Las Vegas Algorithm
A Las Vegas algorithm always returns a correct answer, but its running time may depend on the random choices it makes.
Monte Carlo Algorithm
A Monte Carlo algorithm may return an incorrect answer, but the probability that it does so is bounded.
Randomized Quick Sort is a Las Vegas algorithm. Freivalds’ algorithm and Karger’s algorithm are Monte Carlo algorithms.
4 Challenges
Before I conclude this explainer, I must warn you. As I hope these algorithms have demonstrated, randomized algorithms can be remarkably simple. The algorithms we’ve discussed have minimal steps, and don’t even do any complex operations over the random choice. The random experiment we do to make our random choice (like sampling from a uniform distribution) is also a simple experiment.
You must not confuse this simplicity with easiness. If you’ve learnt enough math, you’ll know that often the simplest ideas are the hardest things to come up with.
When I mentioned the general three-step process to come up with a randomized algorithm, I used the words “good enough.” Let’s expound on what that means.
For a randomized algorithm of the kind we have discussed to be useful, we need two things:
For the random choice we decide to make, we must be able to come up with an expression that bounds the probability of failure:
This is hard. Even though the algorithms we discussed were simple, the way we went about deriving the failure bounds for the algorithms was not obvious. If you’re inspired after reading the explainer and decide to solve a few problems on constructing a randomized algorithm, you’ll immediately run into this issue.
In general, if one run takes time T and we amplify by running the algorithm k times, the total running time of the algorithm is O(kT). This tells us that, for the algorithm to work practically, one run should be quick enough and the function that bounds the failure probability from above should vanish toward 0 quickly enough on repeated runs of the algorithm. If either of these conditions is false, the algorithm might be impractical because it’d simply take too much time to run for a good enough probability bound!
Let’s talk about the second point in more detail. We go back to the Matrix Multiplication Verification problem. This time, we consider a different randomized algorithm to solve it:
Random Column Verification
- Sample a column index j \in \{1,2,\ldots,d\} uniformly at random.
- Compute the j^{\text{th}} column of AB by computing AB_{:,j}=A(B_{:,j}).
- If AB_{:,j}=C_{:,j}, we say that AB=C; otherwise, we say that AB\neq C.
Let F be the event that the algorithm fails.
Let D = AB - C. We notice that, like Freivalds’ algorithm, this algorithm fails only when AB \neq C \iff D \neq 0. Thus to compute the probability of failure we fix A, B, C such that AB \neq C
If we sample j such that D_{:,j} = AB_{:,j} - C_{:,j} \neq 0 then our algorithm correctly reports that AB \neq C.
It’s only when D_{:,j} = 0 that the algorithm fails by incorrectly saying that AB = C.
Thus F = \left\{j \in \{1, \dots, d\} \mid D_{:,j} = 0\right\}.
P[F] = P[D_{:,j} = 0]
Let t be the number of nonzero columns in D. Since D \neq 0, it contains at least one nonzero column.
1 \leq t \leq d
The probability of failure becomes
\begin{aligned} P[F] &= \frac{d-t}{d} \\ &= 1-\frac{t}{d} \\ & \leq 1 - \frac{1}{d} \\ & = \frac{d-1}{d} \end{aligned}
Compared to the Freivalds’ flat bound of \frac{1}{2}, this absolutely sucks. Under amplification, this bound also decreases extremely slowly, and it will take many runs to reduce the failure probability to an acceptable amount.
For a concrete example, let’s go back to multiplying two 3000 \times 3000 matrices on a Compaq Deskpro 386/20 having 160,000 FLOPs per second.
To compute the j^{\text{th}} column of AB, we multiply the d \times d matrix A by the d \times 1 matrix B_{:,j}. From our previous derivation of the FLOPs for multiplying two compatible matrices, multiplying an a \times b matrix by a b \times c matrix takes:
a \times c \times (2b-1) \text{ FLOPs}.
In this case, a=d, b=d, and c=1. Thus the total FLOPs in one run are:
\begin{aligned} \text{Total FLOPs} &= d \times 1 \times (2d-1) \\ &= 2d^2-d \end{aligned}
In our case,
\begin{aligned} d &= 3000 \\ \text{Total FLOPs} &= 3000 \times 1 \times (2\times 3000-1) \\ &= \boxed{17{,}997{,}000\text{ FLOPs}}. \end{aligned}
If you compare it with the total FLOPs for one run of Freivalds’ algorithm, it is significantly less (\approx 18 million vs \approx 54 million).
This suggests that one run of this algorithm is significantly faster than one run of Freivalds’ algorithm.
However, its probability bounds suck. At d = 3000, the probability we fail is at most \frac{2999}{3000} \approx 0.99967! That’s as bad as it gets.
Let’s see how amplification works for this algorithm.
Let F_i be the event that the i^{\text{th}} run fails, and let F^{(k)} be the event that the amplified k-run algorithm fails. If we run the algorithm k times, we fail only if every run fails. The probability of failure becomes:
\begin{aligned} P[F^{(k)}] &= P[F_1 \cap F_2 \cap \dots \cap F_k] \\ &= P[F_1] \times P[F_2] \times \dots \times P[F_k] \\ &\leq \left(\frac{d-1}{d}\right)^k. \end{aligned}
If we want the probability of failure to be at most \delta, where 1>\delta>0, we need:
\left(\frac{d-1}{d}\right)^k \leq \delta.
Taking the natural logarithm of both sides, we get:
k\ln\left(\frac{d-1}{d}\right) \leq \ln(\delta).
Since \ln\left(\frac{d-1}{d}\right)<0, dividing both sides by it reverses the inequality. Thus:
k \geq \frac{\ln(\delta)}{\ln\left(\frac{d-1}{d}\right)}.
Therefore, to get the probability of failure to be at most \delta, we need to repeat the algorithm
k=\left\lceil\frac{\ln(\delta)}{\ln\left(\frac{d-1}{d}\right)}\right\rceil
times.
For our d=3000 matrices, suppose we want the probability of failure to be at most 10\%, given by \delta=0.1. We will need:
\begin{aligned} k &= \left\lceil\frac{\ln(0.1)}{\ln\left(\frac{2999}{3000}\right)}\right\rceil \\ &= \boxed{6907\text{ independent runs}}. \end{aligned}
Since one run takes 17{,}997{,}000 FLOPs, the total FLOPs are:
\begin{aligned} \text{Total FLOPs} &= 6907 \times 17{,}997{,}000 \\ &= \boxed{124{,}305{,}279{,}000\text{ FLOPs}}. \end{aligned}
The time your computer would take to carry out all 6907 runs is given by:
\begin{aligned} \text{Time taken by Compaq 386/20} &= \frac{\text{Total FLOPs}}{\text{FLOPs per second}} \\ &= \frac{124{,}305{,}279{,}000}{160{,}000} \\ &= 776{,}907.99375\,\text{seconds} \\ &\approx 215.81\,\text{hours} \\ &\approx \boxed{8.99\,\text{days}}. \end{aligned}
That’s even more impractical than the 4 days it’d take you to compute AB directly!
What went wrong with this algorithm? The amplified probability bound \left(\frac{d-1}{d}\right)^k simply decreases too slowly as k grows.
This algorithm succumbs to the second challenge of being practical, because it fails at having a probability bound that vanishes quickly enough.
I hope these two challenges - the second one in particular - have given you a sense of how hard designing a practical, working randomized algorithm is.
5 Conclusion
If you’ve read this explainer till here, thank you for sticking around! I remember the feeling of disbelief and wonder when I first studied Randomized Algorithms, and I hope I’ve been able to instill that feeling through this explainer.
If you enjoyed reading this, I suggest reading Probability and Computing by Michael Mitzenmacher and Eli Upfal. I could not cover a lot of algorithms that I think are really cool. I’ll list them here for you to check: (1) Johnson–Lindenstrauss dimensionality reduction (2) Morris counter
This is an annotation. Annotations provide additional context without interrupting the main explanation.
\{0,1\}^d is the set of all length-d vectors whose entries are either 0 or 1. For example, \{0,1\}^2 = \{(0,0),(0,1),(1,0),(1,1)\}. For d=3, examples include (0,0,0), (1,0,1), and (1,1,1); there are 8 vectors in total.
Dongarra’s LINPACK table reports 0.16 Mflop/s for a Compaq 386/20 with an Intel 80387 math coprocessor, using Microsoft Fortran 4.1 on the n=100 benchmark. Since 1 Mflop/s is 10^6 floating-point operations per second, this is 160{,}000 FLOP/s. Source: Jack J. Dongarra, Performance of Various Computers Using Standard Linear Equations Software, University of Tennessee Computer Science Technical Report CS-89-85; also published in ACM SIGARCH Computer Architecture News.
FLOP stands for floating-point operation, and FLOPs is its plural. A floating-point operation is an arithmetic operation on a floating-point number, such as +, -, \times, or \div. A floating-point number is a computer’s finite representation of a real number. For example, to compute 2 + 2 \times 5, where 2, 2, and 5 are floating-point numbers, a computer performs two floating-point operations: one addition and one multiplication. A rate of 160{,}000 FLOPs per second means 160{,}000 such operations per second.
\mathbb{R}^{a \times b} is the set of all a \times b matrices whose entries are real numbers. In this case, A and B are real-valued matrices with 3000 rows and 3000 columns.
Sampling uniformly at random from a finite set S means that every element of S is selected with the same probability, namely 1/|S|.