Randomized Algorithms: QuickSelect, Reservoir Sampling, and Beyond

---

Why Randomize?

Deterministic algorithms must work correctly on every input. An adversary can craft inputs that trigger worst-case behavior (e.g., sorted input for naive quicksort). Randomized algorithms use random choices to make every input equally likely to be fast, eliminating adversarial worst cases.

Two key benefits:

1. Simpler algorithms — randomization often avoids complex bookkeeping
2. Better expected performance — worst case becomes astronomically unlikely

---

QuickSelect: Kth Smallest Element in Expected O(n)

The Problem

Find the kth smallest element in an unsorted array. Sorting solves this in O(n log n), but QuickSelect does it in expected O(n).

The Algorithm

QuickSelect is like QuickSort, but only recurses into one side:

1. Pick a random pivot.
2. Partition the array around the pivot.
3. If the pivot is at position k, return it.
4. If k is to the left, recurse left. If k is to the right, recurse right.

import random

def quickselect(arr, k):
    """
    Find the kth smallest element (0-indexed).
    Expected O(n) time, O(1) extra space (in-place).
    """
    def select(left, right, k_target):
        if left == right:
            return arr[left]
        
        # Random pivot to avoid worst case
        pivot_idx = random.randint(left, right)
        pivot_idx = partition(left, right, pivot_idx)
        
        if k_target == pivot_idx:
            return arr[pivot_idx]
        elif k_target < pivot_idx:
            return select(left, pivot_idx - 1, k_target)
        else:
            return select(pivot_idx + 1, right, k_target)
    
    def partition(left, right, pivot_idx):
        pivot_val = arr[pivot_idx]
        # Move pivot to end
        arr[pivot_idx], arr[right] = arr[right], arr[pivot_idx]
        
        store_idx = left
        for i in range(left, right):
            if arr[i] < pivot_val:
                arr[store_idx], arr[i] = arr[i], arr[store_idx]
                store_idx += 1
        
        # Move pivot to its final position
        arr[store_idx], arr[right] = arr[right], arr[store_idx]
        return store_idx
    
    return select(0, len(arr) - 1, k)


arr = [8, 3, 5, 1, 9, 2, 7, 4, 6]
print(f"1st smallest: {quickselect(arr[:], 0)}")  # 1
print(f"5th smallest: {quickselect(arr[:], 4)}")  # 5
print(f"9th smallest: {quickselect(arr[:], 8)}")  # 9

Time Complexity Analysis

Expected O(n): On average, a random pivot splits the array into roughly equal halves. The work is:

T(n) = T(n/2) + O(n)
     = n + n/2 + n/4 + ...
     = 2n
     = O(n)

Worst case O(n^2): If we always pick the smallest or largest element as pivot, we only reduce the problem size by 1 each time. But with random pivot selection, the probability of this happening for many rounds is astronomically small.

Iterative Version

def quickselect_iterative(arr, k):
    """Iterative QuickSelect -- avoids recursion stack overflow."""
    left, right = 0, len(arr) - 1
    
    while left < right:
        pivot_idx = random.randint(left, right)
        pivot_val = arr[pivot_idx]
        
        # Move pivot to end
        arr[pivot_idx], arr[right] = arr[right], arr[pivot_idx]
        
        store = left
        for i in range(left, right):
            if arr[i] < pivot_val:
                arr[store], arr[i] = arr[i], arr[store]
                store += 1
        
        arr[store], arr[right] = arr[right], arr[store]
        
        if store == k:
            return arr[store]
        elif store < k:
            left = store + 1
        else:
            right = store - 1
    
    return arr[left]

Three-Way Partition (Handles Duplicates)

def quickselect_3way(arr, k):
    """QuickSelect with 3-way partition for arrays with many duplicates."""
    left, right = 0, len(arr) - 1
    
    while left < right:
        pivot = arr[random.randint(left, right)]
        
        # Three-way partition: [< pivot | == pivot | > pivot]
        lt, gt = left, right
        i = left
        
        while i <= gt:
            if arr[i] < pivot:
                arr[lt], arr[i] = arr[i], arr[lt]
                lt += 1
                i += 1
            elif arr[i] > pivot:
                arr[gt], arr[i] = arr[i], arr[gt]
                gt -= 1
            else:
                i += 1
        
        # arr[lt..gt] are all equal to pivot
        if k < lt:
            right = lt - 1
        elif k > gt:
            left = gt + 1
        else:
            return arr[k]
    
    return arr[left]


arr = [3, 3, 3, 1, 1, 2, 2, 2, 4]
print(quickselect_3way(arr[:], 4))  # 2 (5th smallest)

---

Reservoir Sampling

The Problem

Select k items uniformly at random from a stream of unknown length, using only O(k) memory.

This is crucial for:

• Sampling from a file too large to fit in memory
• Sampling from a live data stream
• A/B testing where the population size is unknown

Algorithm (k=1)

For each item at position i (0-indexed):

• Keep it with probability 1/(i+1)
• If kept, replace the current sample

def reservoir_sample_one(stream):
    """Select one element uniformly at random from a stream."""
    result = None
    
    for i, item in enumerate(stream):
        # Keep item with probability 1/(i+1)
        if random.randint(0, i) == 0:
            result = item
    
    return result


# Verify uniform distribution
from collections import Counter

stream = [10, 20, 30, 40, 50]
counts = Counter()
for _ in range(100000):
    counts[reservoir_sample_one(stream)] += 1

for item in sorted(counts):
    print(f"  {item}: {counts[item] / 100000:.3f}")
# Each should be ~0.200

Algorithm (General k)

def reservoir_sampling(stream, k):
    """
    Select k elements uniformly at random from a stream.
    O(n) time, O(k) space.
    """
    reservoir = []
    
    for i, item in enumerate(stream):
        if i < k:
            # Fill the reservoir first
            reservoir.append(item)
        else:
            # Replace element at random position with probability k/(i+1)
            j = random.randint(0, i)
            if j < k:
                reservoir[j] = item
    
    return reservoir


# Sample 3 items from a stream of 100
stream = range(100)
for _ in range(5):
    sample = reservoir_sampling(iter(stream), 3)
    print(f"  Sample: {sample}")

Proof of Correctness

For any element at position i (0-indexed, i > k-1):

Probability of being in final reservoir = Probability of being selected * Probability of not being replaced later

P(selected) = k / (i+1)
P(not replaced by next item) = 1 - 1/(i+2)  for each subsequent item
...

The math works out to exactly k/n for each element, proving uniformity.

Weighted Reservoir Sampling

When items have different weights:

import math

def weighted_reservoir(stream, k):
    """
    Weighted reservoir sampling (Efraimidis-Spirakis algorithm).
    Each item has (value, weight).
    """
    import heapq
    
    heap = []  # Min-heap of (key, value)
    
    for value, weight in stream:
        # key = random^(1/weight) -- higher weight = higher key
        key = random.random() ** (1.0 / weight) if weight > 0 else 0
        
        if len(heap) < k:
            heapq.heappush(heap, (key, value))
        elif key > heap[0][0]:
            heapq.heapreplace(heap, (key, value))
    
    return [val for _, val in heap]

---

Randomized QuickSort

Why Random Pivot?

Deterministic QuickSort with “first element as pivot” is O(n^2) on sorted input. Randomized QuickSort picks a random pivot each time, making the expected time O(n log n) regardless of input.

def randomized_quicksort(arr):
    """QuickSort with random pivot selection."""
    def sort(lo, hi):
        if lo >= hi:
            return
        
        # Random pivot
        pivot_idx = random.randint(lo, hi)
        arr[pivot_idx], arr[hi] = arr[hi], arr[pivot_idx]
        pivot = arr[hi]
        
        # Partition
        i = lo
        for j in range(lo, hi):
            if arr[j] <= pivot:
                arr[i], arr[j] = arr[j], arr[i]
                i += 1
        arr[i], arr[hi] = arr[hi], arr[i]
        
        sort(lo, i - 1)
        sort(i + 1, hi)
    
    sort(0, len(arr) - 1)
    return arr


arr = list(range(20, 0, -1))  # Worst case for naive quicksort
print(randomized_quicksort(arr))

Expected Comparisons Analysis

The expected number of comparisons for randomized quicksort is:

E[comparisons] = 2n * H_n = 2n * ln(n) + O(n) ≈ 1.39 * n * log2(n)

This is about 39% more comparisons than the information-theoretic lower bound, but the constant is small and the algorithm is very cache-friendly.

---

Skip Lists

A skip list is a randomized data structure that provides O(log n) expected search, insertion, and deletion — similar to a balanced BST but much simpler to implement.

Structure

Multiple layers of sorted linked lists. Each element appears in layer i+1 with probability p (typically 0.5). The bottom layer contains all elements.

import random

class SkipNode:
    def __init__(self, val, level):
        self.val = val
        self.forward = [None] * (level + 1)

class SkipList:
    """Skip list with expected O(log n) search, insert, delete."""
    
    MAX_LEVEL = 16
    P = 0.5
    
    def __init__(self):
        self.header = SkipNode(-float('inf'), self.MAX_LEVEL)
        self.level = 0
    
    def _random_level(self):
        lvl = 0
        while random.random() < self.P and lvl < self.MAX_LEVEL:
            lvl += 1
        return lvl
    
    def search(self, target):
        """Search for target. Expected O(log n)."""
        current = self.header
        
        for i in range(self.level, -1, -1):
            while current.forward[i] and current.forward[i].val < target:
                current = current.forward[i]
        
        current = current.forward[0]
        return current is not None and current.val == target
    
    def insert(self, val):
        """Insert a value. Expected O(log n)."""
        update = [None] * (self.MAX_LEVEL + 1)
        current = self.header
        
        for i in range(self.level, -1, -1):
            while current.forward[i] and current.forward[i].val < val:
                current = current.forward[i]
            update[i] = current
        
        new_level = self._random_level()
        
        if new_level > self.level:
            for i in range(self.level + 1, new_level + 1):
                update[i] = self.header
            self.level = new_level
        
        new_node = SkipNode(val, new_level)
        for i in range(new_level + 1):
            new_node.forward[i] = update[i].forward[i]
            update[i].forward[i] = new_node
    
    def delete(self, val):
        """Delete a value. Expected O(log n)."""
        update = [None] * (self.MAX_LEVEL + 1)
        current = self.header
        
        for i in range(self.level, -1, -1):
            while current.forward[i] and current.forward[i].val < val:
                current = current.forward[i]
            update[i] = current
        
        target = current.forward[0]
        if target and target.val == val:
            for i in range(self.level + 1):
                if update[i].forward[i] != target:
                    break
                update[i].forward[i] = target.forward[i]
            
            while self.level > 0 and self.header.forward[self.level] is None:
                self.level -= 1
            return True
        return False
    
    def display(self):
        """Print the skip list structure."""
        for i in range(self.level, -1, -1):
            current = self.header.forward[i]
            vals = []
            while current:
                vals.append(str(current.val))
                current = current.forward[i]
            print(f"Level {i}: {' -> '.join(vals)}")


sl = SkipList()
for x in [3, 6, 7, 9, 12, 19, 17, 26, 21, 25]:
    sl.insert(x)

sl.display()
print(f"Search 19: {sl.search(19)}")  # True
print(f"Search 15: {sl.search(15)}")  # False

Skip lists are used in Redis (sorted sets), LevelDB, and MemSQL.

---

Bloom Filters

A Bloom filter is a space-efficient probabilistic data structure for set membership testing. It can say:

• “Definitely not in the set” (no false negatives)
• “Probably in the set” (possible false positives)

import hashlib

class BloomFilter:
    """Simple Bloom filter with k hash functions."""
    
    def __init__(self, size, num_hashes):
        self.size = size
        self.num_hashes = num_hashes
        self.bit_array = [False] * size
    
    def _hashes(self, item):
        """Generate k hash values for the item."""
        hashes = []
        for i in range(self.num_hashes):
            h = hashlib.md5(f"{item}_{i}".encode()).hexdigest()
            hashes.append(int(h, 16) % self.size)
        return hashes
    
    def add(self, item):
        """Add an item to the filter."""
        for h in self._hashes(item):
            self.bit_array[h] = True
    
    def might_contain(self, item):
        """Check if item might be in the set (may return false positive)."""
        return all(self.bit_array[h] for h in self._hashes(item))


# Usage
bf = BloomFilter(size=1000, num_hashes=5)

# Add some items
for word in ["apple", "banana", "cherry", "date"]:
    bf.add(word)

# Test membership
print(bf.might_contain("apple"))    # True (correct)
print(bf.might_contain("banana"))   # True (correct)
print(bf.might_contain("grape"))    # False (correct)
print(bf.might_contain("fig"))      # False (likely correct)

False Positive Probability

With m bits, k hash functions, and n inserted elements:

P(false positive) = (1 - e^(-kn/m))^k

Optimal number of hash functions: k = (m/n) * ln(2)

import math

def bloom_filter_params(n, desired_fp_rate):
    """Calculate optimal Bloom filter parameters."""
    # Optimal bits
    m = int(-n * math.log(desired_fp_rate) / (math.log(2) ** 2))
    # Optimal hash functions
    k = int((m / n) * math.log(2))
    # Actual false positive rate
    fp_rate = (1 - math.exp(-k * n / m)) ** k
    
    return {
        'bits': m,
        'hash_functions': k,
        'fp_rate': fp_rate,
        'bytes': m // 8,
    }


# For 1 million items with 1% false positive rate
params = bloom_filter_params(1_000_000, 0.01)
print(f"Need {params['bits']} bits ({params['bytes']} bytes)")
print(f"Using {params['hash_functions']} hash functions")
print(f"Actual FP rate: {params['fp_rate']:.4f}")

Bloom filters are used in databases (avoid disk reads for non-existent keys), web browsers (safe browsing), and distributed systems (reduce network calls).

---

Monte Carlo vs Las Vegas Algorithms

Property	Monte Carlo	Las Vegas
Correctness	May be wrong	Always correct
Running time	Always bounded	May vary
Example	Miller-Rabin primality	Randomized QuickSort
Guarantee	P(correct) > 1 - epsilon	E[time] = O(f(n))

Monte Carlo Example: Approximate Pi

def estimate_pi(num_samples):
    """Estimate pi using Monte Carlo method."""
    inside_circle = 0
    
    for _ in range(num_samples):
        x = random.random()
        y = random.random()
        if x * x + y * y <= 1:
            inside_circle += 1
    
    return 4 * inside_circle / num_samples


for n in [1000, 10000, 100000, 1000000]:
    pi_est = estimate_pi(n)
    print(f"n={n:>8}: pi ≈ {pi_est:.4f} (error: {abs(pi_est - math.pi):.4f})")

Monte Carlo Example: Randomized Primality Testing

def miller_rabin(n, k=20):
    """
    Miller-Rabin primality test.
    Returns True if n is probably prime (error probability < 4^(-k)).
    """
    if n < 2:
        return False
    if n < 4:
        return True
    if n % 2 == 0:
        return False
    
    # Write n-1 as 2^r * d
    r, d = 0, n - 1
    while d % 2 == 0:
        r += 1
        d //= 2
    
    for _ in range(k):
        a = random.randrange(2, n - 1)
        x = pow(a, d, n)
        
        if x == 1 or x == n - 1:
            continue
        
        for _ in range(r - 1):
            x = pow(x, 2, n)
            if x == n - 1:
                break
        else:
            return False
    
    return True


print(miller_rabin(104729))           # True (prime)
print(miller_rabin(104730))           # False
print(miller_rabin(2**61 - 1))        # True (Mersenne prime)

Las Vegas Example: Randomized QuickSort

Always produces a correctly sorted array, but running time varies. Expected O(n log n), worst case O(n^2) with vanishingly small probability.

---

Randomized Algorithm for Minimum Cut (Karger’s Algorithm)

import random
from collections import defaultdict
import copy

def karger_min_cut(graph_edges, n):
    """
    Karger's randomized algorithm for minimum cut.
    Run O(n^2 log n) times for high probability of finding the actual min cut.
    """
    # Represent graph as adjacency list with parallel edges
    def contract():
        parent = list(range(n))
        
        def find(x):
            while parent[x] != x:
                parent[x] = parent[parent[x]]
                x = parent[x]
            return x
        
        def union(x, y):
            px, py = find(x), find(y)
            parent[px] = py
        
        edges = graph_edges[:]
        random.shuffle(edges)
        
        vertices = n
        edge_idx = 0
        
        while vertices > 2 and edge_idx < len(edges):
            u, v = edges[edge_idx]
            edge_idx += 1
            
            if find(u) != find(v):
                union(u, v)
                vertices -= 1
        
        # Count crossing edges
        cut = 0
        for u, v in graph_edges:
            if find(u) != find(v):
                cut += 1
        
        return cut
    
    # Run multiple times and take the minimum
    min_cut = float('inf')
    trials = n * n  # n^2 trials for good probability
    
    for _ in range(trials):
        cut = contract()
        min_cut = min(min_cut, cut)
    
    return min_cut

---

Summary: When to Use Randomization

Problem	Deterministic	Randomized	Improvement
Kth element	O(n log n) sort	O(n) QuickSelect	Log factor
Sampling k from n	O(n) + shuffle	O(n) reservoir	Works on streams
Sorting	O(n log n)	O(n log n) expected	No adversarial worst case
Balanced BST	AVL/Red-Black (complex)	Skip list (simple)	Simpler code
Set membership	Hash set O(n) space	Bloom filter O(n) bits	Space savings
Primality	AKS O(n^6)	Miller-Rabin O(k log^2 n)	Much faster

---

Practice Problems

1. Kth Largest Element (LeetCode 215) — QuickSelect.
2. Shuffle an Array (LeetCode 384) — Fisher-Yates shuffle.
3. Random Pick with Weight (LeetCode 528) — Weighted sampling.
4. Linked List Random Node (LeetCode 382) — Reservoir sampling (k=1).
5. Random Pick Index (LeetCode 398) — Reservoir sampling for duplicates.
6. Implement Rand10 Using Rand7 (LeetCode 470) — Rejection sampling.
7. Monte Carlo integration — Estimate area under a curve.
8. Skip List Design (LeetCode 1206) — Implement a skip list.

---

Key Takeaways

• QuickSelect finds the kth element in expected O(n) by partitioning around a random pivot and only recursing into one side.
• Reservoir sampling selects k items uniformly from a stream of unknown length in O(n) time and O(k) space.
• Randomized QuickSort with random pivot avoids adversarial worst cases, giving expected O(n log n) on any input.
• Skip lists provide a simple alternative to balanced BSTs with O(log n) expected operations using randomized level assignment.
• Bloom filters trade a small false positive probability for massive space savings in set membership testing.
• Monte Carlo algorithms are fast but may err; Las Vegas algorithms are always correct but may be slow. Both use randomization to great effect.