6.046J Complete Lecture Notes

Contents 1 Introduction & Median Finding 1.1 The Course 1.2 Algorithms 1.3 Order Statistics (“Median Finding”) 2 Recap & Interval Scheduling 2.1 Recap of Median Finding 2.2 Interval Scheduling 3 Minimum Spanning Trees I 3.1 Greedy Algorithms 3.2 Graphs 3.3 Minimum Spanning Trees 4 Minimum Spanning Trees II 4.1 Implementing Kruskal’s Algorithm 4.2 Prim’s Algorithm 4.3 Greedy Strategies 5 Fast Fourier Transform 5.1 Multiplication 5.2 Convolution 6 All-Pairs Shortest Paths I 6.1 Dynamic Programming 6.2 Shortest-Path Problems 6.3 The Floyd–Warshall Algorithm 7 All-Pairs Shortest Paths II 7.1 Johnson’s Algorithm 7.2 Linear Programming 8 Randomized Algorithms I 8.1 Randomized Median Finding 8.2 Another Example: Verifying Polynomial Multiplication 8.3 The Markov Bound 9 Randomized Algorithms II 9.1 The Central Limit Theorem and the Chernoff Bound 9.2 Analysis of Q UICKSORT 9.3 Monte Carlo Sampling 9.4 Amplification 10 Hashing and Amortization 10.1 Arrays and Hashing 10.2 Hash Families 10.3 Amortization 11 Amortized Analysis 11.1 Aggregate Analysis 11.2 Accounting Method 11.3 Potential Method 11.4 Example: A Dynamically Resized Table 12 Competitive Analysis 12.1 Online and Offline Algorithms 12.2 Example: A Self-Organizing List 13 Network Flow 13.1 The Ford–Fulkerson Algorithm 13.2 The Max Flow–Min Cut Equivalence 13.3 Generalizations 14 Interlude: Problem Solving 14.1 What to Bring to the Table 14.2 How to Attack a Problem 14.3 Recommended Reading 15 van Emde Boas Data Structure 15.1 Analogy: The Two-Coconut Problem 15.2 Implementation: A Recursive Data Structure 15.3 Solving the Recurrence 16 Disjoint-Set Data Structures 16.1 Linked-List Implementation 16.2 Forest-of-Trees Implementation 17 Complexity and NP-completeness 17.1 Examples 17.2 Complexity 17.3 Example: Algorithm Search 18 Polynomial-Time Approximations 18.1 Vertex Cover 18.2 Set Cover 18.3 Partition 19 Compression and Huffman Coding 19.1 Compression 19.2 The Huffman Algorithm 20 Sublinear-Time Algorithms 20.1 Estimating the Number of Connected Components 20.2 Estimating the Size of a Minimum Spanning Tree 21 Clustering 21.1 Hierarchical Agglomerative Clustering 21.2 Minimum-Radius Clustering 22 Derandomization 22.1 Using Conditional Expectation 22.2 Using Pairwise Independence 23 Computational geometry 23.1 Intersection Problem 23.2 Finding the Closest Pair of Points Index Lecture 1 Introduction & Median Finding Supplemental reading in CLRS: Section 9.3; Chapter 1; Sections 4.3 and 4.5 1.1 The Course Hello, and welcome to 6.046 Design and Analysis of Algorithms. The prerequisites for this course are 1. 6.006 Introduction to Algorithms. This course is designed to build on the material of 6.006. The algorithms in 6.046 are generally more advanced. Unlike 6.006, this course will not require students to implement the algorithms; instead, we make frequent use of pseudocode. 2. Either 6.042/18.062J Mathematics for Computer Science or 18.310 Principles of Applied Mathematics. Students must be familiar with beginning undergraduate mathematics and able to write mathematical proofs. After taking this course, students should be able not just to describe the algorithms they have learned, but also to prove their correctness (where applicable) and rigorously establish their asymptotic running times. The course topics are as follows: • • • • • • • Divide and Conquer Dynamic Programming Greedy Algorithms Graph Algorithms Randomized Algorithms Data Structures Approximation Algorithms. The course textbook is Introduction to Algorithms 3e., by Cormen, Leiserson, Rivest and Stein. We will usually refer to the textbook (or the authors) as “CLRS.” In addition, Profs. Bruce Tidor and Dana Moshkovitz, who taught this course in spring of 2012, prepared hand-written notes for each lecture. What you are reading is an electronic, fleshed-out version of those notes, developed by Ben Zinberg with oversight by Prof. Moshkovitz under the commission of MIT OpenCourseWare. 1.2 Algorithms Before our first example, let’s take a moment to ask: What is the study of algorithms, and why is it important? CLRS offer the following answer: An algorithm is “any well-defined computational procedure that takes some [. . . ] input and produces some [. . . ] output.” They point out that algorithms are an essential component not just of inherently electronic domains such as electronic commerce or the infrastructure of the Internet, but also of science and industry, as exemplified by the Human Genome Project or an airline that wishes to minimize costs by optimizing crew assignments. Studying algorithms allows us both to understand (and put to use) the capabilities of computers, and to communicate our ideas about computation to others. Algorithms are a basis for design: they serve as building blocks for new technology, and provide a language in which to discuss aspects of that technology. The question of whether a certain algorithm is effective depends on the context in which it is to be used. Often relevant are the following two concerns: • Correctness. An algorithm is said to be correct if, for every possible input, the algorithm halts with the desired output. Usually we will want our algorithms to be correct, though there are occasions in which incorrect algorithms are useful, if their rate of error can be controlled. • Efficiency. The best algorithms are ones which not just accomplish the desired task, but use minimal resources while doing so. The resources of interest may include time, money, space in memory, or any number of other “costs.” In this course we will mostly be concerned with time costs. In particular, we will try to give asymptotic bounds on the “number of steps” required to perform a computation as a function of the size of the input. The notion of a “step” is an imprecise (until it has been given a precise definition) abstraction of what on current computers could be processor cycles, memory operations, disk operations, etc. Of course, these different kinds of steps take different amounts of time, and even two different instructions completed by the same processor can take different amounts of time. In this course, we do not have the time to study these subtleties in much detail. More advanced courses, such as 6.854 Advanced Algorithms, take up the subject in further detail. Fortunately, the course material of 6.046 is extremely useful despite its limitations. 1.3 Order Statistics (“Median Finding”) Median finding is important in many situations, including database tasks. A precise statement of the problem is as follows: Problem 1.1. Given an array A = A[1, . . . , n] of n numbers and an index i (1 ≤ i ≤ n), find the ithsmallest element of A. For example, if i = 1 then this amounts to finding the minimum element, and if i = n then we are 1 looking for the maximum element. If n is odd, then setting i = n+ gives the median of A. (If n is 2 ¥ ¦ § 1¨ 1 even, then the median will probably be some average of the cases i = n+ and i = n+ 2 2 , depending on your convention.) Intuitive Approach. We could simply sort the entire array A—the ith element of the resulting array would be our answer. If we use M ERGE -S ORT to sort A in place and we assume that jumping to the ith element of an array takes O(1) time, then the running time of our algorithm is Lec 1 – pg. 2 of 7 (M ERGE -S ORT) (jump to the ith element) O(n lg n) O(1) O(n lg n) + (worst case). Can we expect to do better than this initial approach? If sorting all of A were a necessary step in computing any individual median, then the answer would be no. But as it stands, sorting all of A does much more work than we are asking for—it finds all n medians, whereas we only wanted one. So, even though the best sorting algorithms take Θ(n lg n) time, it does not necessarily follow that median finding cannot be done faster. It turns out that median finding can be done in O(n) time. In practice, this is usually accomplished by a randomized algorithm with linear expected running time, but there also exists a deterministic algorithm with linear worst-case time. In a 1973 paper, Blum, Floyd, Pratt, Rivest and Tarjan proposed the so-called “median-of-medians” algorithm, which we present below. For this algorithm, we will assume that the elements of A are distinct. This is not a very strong assumption, and it is easy to generalize to the case where A is allowed to have duplicate elements.1 Algorithm: S ELECT(A, i) 1. Divide the n items into groups of 5 (plus any remainder). 2. Find the median of each group of 5 (by rote). (If the remainder group has an even number of elements, then break ties arbitrarily, for example by choosing the lower median.) 3. Use S ELECT recursively to find the median (call it x) of these d n/5e medians. 4. Partition around x.∗ Let k = rank(x).† | 5. Ai < x {z k − 1 elements x } ↑ Ak | Ai > x {z n − k elements } • If i = k, then return x. • Else, if i < k, use S ELECT recursively by calling S ELECT(A[1, . . . , k − 1], i).‡ • Else, if i > k, use S ELECT recursively by calling S ELECT(A[k + 1, . . . , i], i − k). ∗ By “partition around x,” we mean reorder the elements of A in place so that all elements prior to x are less than x (and consequently all elements after x are ≥ x—in our case actually > x since there are no duplicate elements). Partitioning can be done in linear time. We won’t show the details here, but they can be found in §7.1 of CLRS. † For a set S of distinct numbers, we define the rank of an element x ∈ S to be the number k such that x is the kth-smallest element of S. Thus, in the set {5, 8, 2, 3}, the rank of 5 is 3. ‡ The array A[1, . . . , k − 1] should be passed by reference—there is certainly no need to waste time and memory by making a new copy. 1 One example of such a way is to equip each element of A with satellite data. Replace each element A[ j] (1 ≤ j ≤ n) with the pair 〈 A[ j], j 〉, thus turning A into an array of ordered pairs. Then we can define an order on A by saying 〈 x1 , j 1 〉 < 〈 x2 , j 2 〉 iff either x1 < x2 , or x1 = x2 and j 1 < j 2 . (This is the so-called “lexicographic order.”) Since this order is strict, it is safe to feed into our algorithm which assumes no duplicate elements. On the other hand, once A is sorted in the lexicographic order, it will also be sorted in the regular order. Lec 1 – pg. 3 of 7 Note a few things about this algorithm. First of all, it is extremely non-obvious. Second, it is recursive, i.e., it calls itself. Third, it appears to do “less work” than the intuitive approach, since it doesn’t sort all of A. Thus, we might suspect that it is more efficient than the intuitive approach. 1.3.1 Proof of Correctness The strategy to prove correctness of an algorithm is a mix of knowing standard techniques and inventing new ones when the standard techniques don’t apply. To prove correctness of our median-ofmedians algorithm, we will use a very common technique called a loop invariant. A loop invariant is a set of conditions that are true when the loop is initialized, maintained in each pass through the loop, and true at termination. Loop invariant arguments are analogous to mathematical induction: they use a base case along with an inductive step to show that the desired proposition is true in all cases. One feasible loop invariant for this algorithm would be the following: At each iteration, the “active subarray” (call it A 0 ) consists of all elements of A which are between min(A 0 ) and max(A 0 ), and the current index (call it i 0 ) has the property that the i 0 th-smallest element of A 0 is the ith-smallest element of A. In order to prove that this is a loop invariant, three things must be shown: 1. True at initialization. Initially, the active subarray is A and the index is i, so the proposition is obviously true. 2. Maintained in each pass through the loop. If the proposition is true at the beginning of an iteration (say with active subarray A 0 = A 0 [1, . . . , n0 ] and active index i 0 ), then the (k − 1)- and (n0 − k)-element subarrays depicted in step 4 clearly also have the contiguity property that we desire. Then, if step 5 calls S ELECT(A 00 , i 00 ), it is easy to see by casework that the i 00 thsmallest element of A 00 equals the i 0 th-smallest element of A 0 , which by hypothesis equals the ith-smallest element of A. 3. True at termination. Since there is nothing special about the termination step (it is just whichever step happens to be interrupted by the returned value), the proof of maintenance in each pass though the loop is sufficient to show maintenance at termination. Now that we have established the loop invariant, it is easy to see that our algorithm is correct. If the final recursive iteration of S ELECT has arguments A 0 and i 0 , then step 5 returns the i 0 th-smallest element of A 0 , which by the loop invariant equals the ith-smallest element of A. Thus we have proven that, if our algorithm terminates at all, then it terminates with the correct output. Moreover, the algorithm must terminate eventually, since the size of the active subarray shrinks by at least 1 with each recursive call and since S ELECT(A 0 , i 0 ) terminates immediately if A 0 has one element. 1.3.2 Running Time There is no single magic bullet for writing proofs. Instead, you should use nonrigorous heuristics to guide your thinking until it becomes apparent how to write down a rigorous proof. In order to prove that our algorithm is efficient, it would help to know that the active subarray in each successive recursive call is much smaller than the previous subarray. That way, we are guaranteed that only relatively few recursive calls are necessary. Therefore, let’s determine an upper bound on the size of A 0 , the new active subarray, given that the old subarray, A, had size n. Lec 1 – pg. 4 of 7 larger elements within each group x smaller elements within each group larger medians smaller medians Figure 1.1. Imagine laying out the groups side by side, ordered with respect to their medians. Also imagine that each group is sorted from greatest to least, top to bottom. Then each of the elements inside the red curve is guaranteed to be greater than x. As in step 3, let x be the median of medians. Then either A 0 doesn’t contain any items greater than x, or A 0 doesn’t contain any items less than x. However, as we will see, there are lots of elements greater than x and there are lots of elements less than x. In particular, each (non-remainder) group whose median is less than x contains at least three elements less than x, and each (non-remainder) group whose median is greater than x contains at least three elements greater than x. Moreover, ¥ ¦ there are at least d n/5e − 1 non-remainder groups, of which at least 21 d n/5e − 2 have median less ¥1 ¦ than x and at least 2 d n/5e − 2 have mean greater than x. Finally, the group with x as its median contains two elements greater than x and two elements less than x (unless n < 5, but you can check the following inequalities in that case separately if you like). Thus ¡ ¦¢ ¡ ¢ ¢ ¡¥ # elts of A less than x ≥ 3 12 d n/5e − 2 + 2 ≥ 3 12 d n/5e − 2 − 12 + 2 > 3 10 n − 6. Similarly, ¡ ¢ # elts of A greater than x ≥ 3 10 n − 6. 3 Therefore A 0 must exclude at least 10 n − 6 elements of A, which means that A 0 contains at most ¡3 ¢ 7 n − 10 n − 6 = 10 n + 6 elements. See Figure 1.1 for an illustration. Let’s now put together all the information we know about the running time of S ELECT. Let T(n) denote the worst-case running time of S ELECT on an array of size n. Then: • Step 1 clearly takes O(n) time. • Step 2 takes O(n) time since it takes O(1) time to find the median of each 5-element group. • Step 3 takes T (d n/5e) time since there are d n/5e submedians. • Step 4 takes O(n) time, as explained in §7.1 of CLRS. ¢ ¡7 • Step 5 takes at most T 10 n + 6 since the new subarray has at most 7 10 n + 6 elements.2 2 Technically, we haven’t proved that T(n) is an increasing function of n. However, if we let T 0 (n) = max {T(m) : m ≤ n}, then T 0 (n) is an increasing function of n, and any asymptotic upper bound on T 0 (n) will also be an asymptotic upper bound on T(n). Lec 1 – pg. 5 of 7 Thus T(n) ≤ T (d n/5e) + T 1.3.3 ¡ 7 10 n + 6 ¢ + O(n). (1.1) Solving the Recurrence How does one go about solving a recurrence like (1.1)? Possible methods include: • Substitution Method: make a guess and prove it by induction. • Recursion Tree Method: branching analysis. • Master Method: read off the answer from a laundry list of already-solved recurrences.3 We will use the substitution method this time. Our guess will be that T(n) is O(n), i.e., that there exists a positive constant c such that T(n) ≤ cn for all sufficiently large n. We won’t choose c just yet; later it will become apparent which value of c we should choose. In order to deal with the base case of an inductive argument, it is often useful to separate out the small values of n. (Later in the argument we will see how this helps.) Since S ELECT always terminates, we can write ( O(1), n < 140 T(n) ≤ ¡7 ¢ T (d n/5e) + T 10 n + 6 + O(n), for all n. In other words, there exist positive constants a, b such that T(n) ≤ b for n < 140 (for example, b = ¡7 ¢ max {T(m) : m < 140}) and T(n) ≤ T (d n/5e) + T 10 n + 6 + an for all n. Once this is done, we can start our inductive argument with the base case n = 140: Base case. T(140) ≤ T (d140/5e) + T ¡ 7 10 140 + 6 ¢ + an ≤ b + b + an = 2b + an = 2b + 140a. b Thus, the base case holds if and only if 2b + 140 ≤ c(140), or equivalently c ≥ a + 70 . Inductive step. Suppose n ≥ 140. Then the inductive hypothesis gives ¡7 ¢ n + 6 + an T(n) ≤ c (d n/5e) + c 10 ¢ ¡7 ¢ ¡ n + 6 + an ≤ c n5 + 1 + c 10 = 9 10 cn + 7c + an = cn + ³ − cn ´ + 7c + an . | 10 {z } Thus all we need is for the part with the brace to be nonpositive. This will happen if and only if n c ≥ n10an −70 . On the other hand, since n ≥ 140, we know n−70 ≤ 2. (This is why we wanted a base case of n = 140. Though, any base case bigger than 70 would have worked.) Thus the part with the brace will be nonpositive whenever c ≥ 20a. 3 The master method is a special case of a general technique known as literature search, which is of fundamental importance in not just computer science but also physics, biology, and even sociology. Lec 1 – pg. 6 of 7 Putting these two together, we find that the base case and the inductive n step will o both be satisfied b b as long as c is greater than both a + 70 and 20a. So if we set c = max 20a, 70 + 1, then we have T(n) ≤ cn for all n ≥ 140. Thus T(n) is O(n). Lec 1 – pg. 7 of 7 Lecture 2 Recap & Interval Scheduling Supplemental reading in CLRS: Section 16.1; Section 4.4 2.1 Recap of Median Finding Like M ERGE -S ORT, the median-of-medians algorithm S ELECT calls itself recursively, with the argument to each recursive call being smaller than the original array A. However, while M ERGE -S ORT runs in O(n lg n) time, S ELECT is able to run in linear time. Why is this? There are two key observations to make: Observation 1. Given an element x of A, we can partition around x in linear time. Thus, all the steps in S ELECT other than the recursive calls take a total of only O(n) time. Observation 2. Finding the median of n elements can be reduced to finding the median of n/5 7 elements and then finding the median of at most 10 n elements. (The true figures in our case are 7 d n/5e and 10 n + 6, but that’s not important for the analysis that follows. For a refresher on where this observation comes from, see the proof of running time and Figure 1.1.) To see why the second observation is important, we’ll do a branching analysis. This will also show how we might have figured out the running time of our algorithm without guessing. Observation 2 is just a translation into English of the recurrence relation ¡ ¢ ¡7 ¢ T(n) ≤ T 15 n + T 10 n + O(n). (2.1) ¡ ¢ ¡7 ¢ Pick a constant c such that T(n) ≤ T 51 n + T 10 n + cn for all sufficiently large n. (Actually, for simplicity, let’s assume the relation holds for all n.1 ) Then we can apply (2.1) to each of the terms ¡ ¢ ¡7 ¢ T 15 n and T 10 n , obtaining ¡ ¢ ¡7 ¢ T(n) ≤ T 51 n + T 10 n + cn h ¡ ¢ ¡7 1 ¢ ¡ ¢i h ¡ ¢ ¡7 7 ¢ ¡ 7 ¢i 7 ≤ T 15 · 51 n + T 10 · 5 n + c 51 n + T 51 · 10 n + T 10 · 10 n + c 10 n + cn | {z } | {z } ≤ ··· ≤ cn + c ¡1 5 ¢ ¡ ¢ 7 7 7 1 7 7 n + c 51 · 15 + 15 · 10 n +··· + 10 + 10 · 5 + 10 · 10 1 Make sure you see why we are allowed to do this. n 1 5n 1 1 5 5n ¡ ¢ 7 10 n 7 1 10 5 n ¡ ¢ 1 7 5 10 n ¡ ¢ 7 7 10 10 n ¡ ¢ Figure 2.1. A recursion tree for S ELECT. Running S ELECT on an array of n elements gives rise to two daughter processes, 7 n. Each of these daughters produces two of its own daughters, with arguments of sizes with arguments of sizes 51 n and 10 ³ ´ ³ ´ 1 1 7 1 n , n , etc. 5 5 10 5 ¢ ¡ ¢ ¡ ¢ ¡ 7 7 2 7 3 n + c 51 + 10 n + c 51 + 10 n +··· = cn + c 15 + 10 " # ∞ ¡ X ¢ 1 7 j =c n (a geometric sum) 5 + 10 j =0 = 10cn. 7 This geometric sum makes it clear that it is very important for 15 + 10 to be less than 1—otherwise, 1 the geometric sum would diverge. If the number 5 in our recurrence had been replaced by 25 , the resulting recurrence would have non-linear solutions and the running time would be T(n) = Θ(n2 ). ¡ ¢ 7 9 On an intuitive level, this makes sense: reducing an input of n entries to an input of 51 + 10 n n = 10 ¡2 7 ¢ 11 entries is a good time-saver, but “reducing” an input of n entries to an input of 5 + 10 n = 10 n entries is not. Finally, let’s make an observation that one is not usually able to make after writing an algorithm: “S ELECT achieves the best possible asymptotic running time for an algorithm that solves the median finding problem.” This statement must be interpreted carefully. What we mean is that any correct algorithm which solves the median finding problem must be Ω(n), since such an algorithm must at least look at each of the n entries of A. However, there are several other considerations that prevent us from being able to conclude that S ELECT is the “best” median finding algorithm: • It is only the asymptotic running time of S ELECT that is optimal. The overhead could be very large. There may exist algorithms with a larger asymptotic running time that are still much faster than S ELECT for all practical purposes because the asymptotic behavior only becomes significant when n is inhumanly large. • There may also exist algorithms with the same asymptotic running time as S ELECT but whose implementations run faster than those of S ELECT because they have smaller overhead (i.e., the constants implicit in the big-O notation are smaller). • There may exist randomized algorithms which run in sublinear time and return the correct answer with high probability, or deterministic algorithms which run in sublinear time and Lec 2 – pg. 2 of 7 return approximately the correct answer.2 2.2 Interval Scheduling Does your planner look like this? seminar group meeting lecture lecture 8am 9 10 Noon 1pm sports game work group lunch 11 party language class lecture 2 3 4 5 6 7 8 9 10 11 Midnight The first step to getting the most out of your day is getting enough sleep and eating right. After that, you’ll want to make a schedule that is in some sense “maximal.” There are many possible precise formulations of this problem; one of the simpler ones is this: ­ ® Input: A list of pairs L = (s 1 , f 1 ), . . . , (s n , f n ) , representing n events with start time s i and end time f i (s i < f i ). ­ ® Output: A list of pairs S = (s i 1 , f i 1 ), . . . , (s i k , f i k ) representing a schedule, such that s i1 < f i1 ≤ s i2 < f i2 ≤ · · · ≤ s i k < f i k (i.e., no scheduled events overlap), where k is maximal. This formulation corresponds to wanting to attend as many events as possible (not allowing late arrivals or early departures, without any regard for how long the events are, without any preference of one event over another, and assuming it takes no time to travel between events). While all of us face problems like this in everyday life, probably very few of us have tried to mathematically rigorize the heuristics we use to decide what to do each day. Often we employ greedy strategies. A greedy strategy can be used in situations where it is easy to tell what is “locally optimal.” The hope is that by making locally optimal decisions at each point, we will arrive at a globally optimal solution.3 A couple of reasonable-sounding greedy strategies are: • Pick the activities that take the least time (minimize f i − s i ). 2 While such a sublinear algorithm is probably out of the question in our case (at least if we have no prior information about A), there are lots of problems in which this does occur. For example, if you believe P 6= NP, we will see later in the course that there exist polynomial-time approximations to some NP-hard problems. 3 An analogous situation is this: Imagine that you live in the Himalayan village of Nagarkot and are trying to find the highest point in the Himalayan mountain range. A greedy strategy would be to hike straight up. Since this strategy is locally optimal, you will eventually arrive at a peak from which it is impossible to go up without first going down. However, your peak will probably not be the global optimum, since Nagarkot is in central Nepal and Mt. Everest is on the Chinese border. Lec 2 – pg. 3 of 7 This strategy fails because attending one short activity may prevent you from attending two other activities. • Each time you find yourself doing nothing, go to the activity that starts the soonest (locally minimize s i ). This strategy fails because it is possible for the earliest activity to preclude lots of other activities. Somewhat surprisingly, there is a greedy strategy that is optimal: • Each time you find yourself doing nothing, go to the activity that ends the soonest (locally minimize f i ). In pseudocode, while there are still activities to consider do Pick (s i , f i ) with the smallest f i 3 Remove all activities that intersect (s i , f i ) ­ ® Claim. The above algorithm outputs a list of pairs (s i 1 , f i 1 ), . . . , (s i k , f i k ) such that 1 2 s i1 < f i1 ≤ s i2 < f i2 ≤ · · · ≤ s i k < f i k . Proof. By definition, s i < f i . Now, if s i 2 < f i 1 , then (s i 2 , f i 2 ) overlaps (s i 1 , f i 1 ). Thus (s i 2 , f i 2 ) must have been removed from consideration after we chose (s i 1 , f i 1 ), a contradiction. The same reasoning with 1 and 2 replaced by j and j + 1 for arbitrary j completes the proof. Claim. Suppose given a list of activities L. If an optimal schedule has k∗ activities, then the above algorithm outputs a schedule with k∗ activities. Proof. By induction on k∗ . If the optimal schedule has only one activity, then obviously the claim holds. Next, suppose the claim holds for k∗ and we are given a set of events whose optimal schedule ® ­ has k∗ + 1 activities. Let S ∗ = S ∗ [1, . . . , k∗ + 1] = (s `1 , f `1 ), . . . , (s `k∗ +1 , f `k∗ +1 ) be an optimal schedule, ­ ® and let S = S[1, . . . , k] = (s i 1 , f i 1 ), . . . , (s i k , f i k ) be the schedule that our algorithm gives. By construc­ ® tion, f i 1 ≤ f `1 . Thus the schedule S ∗∗ = (s i 1 , f i 1 ), (s `2 , f `2 ), . . . , (s `k∗ +1 , f `k∗ +1 ) has no overlap, and is also optimal since it has k∗ + 1 activities. Let L0 be the set of activities with s i ≥ f i 1 . The fact that S ∗∗ is optimal for L implies that S ∗∗ [2, . . . , k∗ + 1] is optimal for L0 . Thus an optimal schedule for L0 has k∗ activities. So by the inductive hypothesis, running our algorithm on L0 produces an optimal schedule. But by construction, running our algorithm on L0 just gives S[2, . . . , k]. Since any two optimal schedules have the same length, we have k = k∗ + 1. So S is an optimal schedule, completing the induction. This algorithm can be implemented with Θ(n lg n) running time as follows: Lec 2 – pg. 4 of 7 Algorithm: U NWEIGHTED -S CHEDULE(L) 1 2 3 4 5 6 7 8 Sort L according to finish time using M ERGE -S ORT S ← 〈〉 curf ← −∞ B the finish time of our current schedule for (s, f ) in L do if s ≥ curf then S.append[(s, f )] curf ← f return S This implementation makes one pass through the list of activities, ignoring those which overlap with the current schedule. Notice that the only part of this algorithm that requires Θ(n lg n) time is the sorting—once L is sorted, the remainder of the algorithm takes Θ(n) time. 2.2.1 Weighted Interval Scheduling Let’s try a more sophisticated formulation of the scheduling problem. Suppose that rather than maximize the total number of events we attend, we want to maximize the total amount of time we are busy. More generally still, suppose each event in L has a weight w > 0 representing how much we prioritize that event. ­ ® Input: A list of triples L = (s 1 , f 1 , w1 ), . . . , (s n , f n , wn ) , representing n events with start time s i , end time f i (s i < f i ), and weight w i > 0. ­ ® Output: A list of triples S = (s i 1 , f i 1 , w i 1 ), . . . , (s i k , f i k , w i k ) representing a schedule, such that s i1 < f i1 ≤ s i2 < f i2 ≤ · · · ≤ s i k < f i k (i.e., no scheduled events overlap), where Pk j =1 w i j is maximal. A greedy algorithm of the sort above will not work for this problem. 1 1 10 1 1 Which is better: attending a short event with low weight (so that you have more free time for the heavy events), or attending a slightly longer event with slightly higher weight? It depends on whether there are exciting opportunities in the future. 1 1 2 2 Here is a sketch of an efficient solution to the weighted scheduling problem: Lec 2 – pg. 5 of 7 100 t1 t4 3 1 t0 B: 1 t2 t3 t5 3 3 1 1 0 0 t0 t1 t2 t3 t4 t5 max{1 + 1, 3} Figure 2.2. An illustration of W EIGHTED -S CHEDULE. Algorithm: W EIGHTED -S CHEDULE(L) Put all the start and finish times into an array A = 〈 s 1 , f 1 , . . . , s n , f n 〉. Sort A, keeping track of which activity each number belongs to (that is, the index i). Make a new array B with one entry for each element of A. The idea is that the value of B[t] will be the maximal weight of a schedule that starts at point t, along with a list of event indices that gives such a schedule. Thus B[t 0 ] (assuming your day starts at t = t 0 ) will give you an optimal schedule for your day. To construct B, first initialize all entries to zero. Then traverse backwards through A. The first time you hit a start point s i first , set B[s i first ] ← w i first . From then on, keep track of the most recent start point you have seen in a new variable i cur . Each time you hit an endpoint f i , set B[ f i ] ← B[s i cur ]. Each time you hit a new start point s i , set © ª B[s i ] ← max B[s i cur ], w i + B[ f i ] . Once you have finished, B[t 0 ] (assuming your day starts at t = t 0 ) will be the weight of an optimal schedule. For ease of exposition, I didn’t say explicitly how to keep track of the actual schedule you are constructing (which you will usually want to do—what use is the maximal weight if you don’t know what schedule gives that weight?). That is an easy implementation detail that you can figure out for yourself. A more subtle issue that I didn’t address is the possibility of duplicate elements of A. One way to deal with this is to, for each new start point s, keep a cache of possible weights of schedules starting at s (using a new variable to keep track of the heaviest schedule in cache). The cache should only be cleared when you hit a start point whose value is strictly less than s (or run out of elements of A), at which point you will know what value to assign to B[s]. Again, it is easy to see that this algorithm runs in linear time after sorting, for a total running time of O(n lg n). This example introduces the algorithmic paradigm of dynamic programming. “Dynamic programming” refers to a strategy in which the solution to a problem is computed by combining the solutions to subproblems which are essentially smaller versions of the same problem. In our example, we started with an optimal schedule for the empty set of events and then added events one at a time (in a strategic order), arriving at an optimal schedule for the full set L of events. Each time we wanted to solve a subproblem, we used the answers to the previous subproblems. Lec 2 – pg. 6 of 7 2.2.2 Conclusions This lecture showed that the way in which a problem is formulated is important. A given real-life problem can be formulated as a precise algorithmic problem in many different, non-equivalent ways, with different solutions. There are several standard approaches that you should have in your toolbox. So far we have seen greedy algorithms and dynamic programming; more will come up later in the course. Lec 2 – pg. 7 of 7 Lecture 3 Minimum Spanning Trees I Supplemental reading in CLRS: Chapter 4; Appendix B.4, B.5; Section 16.2 3.1 Greedy Algorithms As we said above, a greedy algorithm is an algorithm which attempts to solve an optimization problem in multiple stages by making a “locally optimal” decision at each stage. Example. We wish to make 99¢ in change using the minimal number of coins. Most people instinctively use a greedy algorithm: 99¢ −75¢ 24¢ −20¢ 4¢ −4¢ 0¢ = (25¢) × 3+ = (10¢) × 2 + = (1¢) × 4 (no nickels) 3 quarters + 2 dimes + 4 pennies = 9 coins. This greedy algorithm is correct: starting with the coin of largest value, take as many as possible without allowing your total to exceed 99¢. However, plenty of greedy algorithms are not correct. For example, suppose we started with the smallest coin rather than the largest coin: 99¢ = (1¢) × 99 =⇒ 99 pennies. Or, imagine trying to make 15¢ of change if the dime were replaced by an 11¢ piece: Greedy: Optimal: 15¢ = (11¢) × 1 + (5¢) × 0 + (1¢) × 4 15¢ = (11¢) × 0 + (5¢) × 3 + (1¢) × 0 =⇒ =⇒ 5 coins 3 coins. Remark. Greedy algorithms sometimes give a correct solution (are globally optimal). But even when they’re not correct, greedy algorithms can be useful because they provide “good solutions” efficiently. For example, if you worked as a cashier in a country with the 11¢ piece, it would be perfectly reasonable to use a greedy algorithm to choose which coins to use. Even though you won’t always use the smallest possible number of coins, you will still always make the correct amount of change, and the number of coins used will always be close to optimal. 3 2 1 1 3 Figure 3.1. From left to right: A directed graph, an undirected graph, and a weighted undirected graph. The above example shows that the performance of a greedy algorithm (in terms of both correctness and efficiency) depends on both • the structure of the algorithm (starting with big coins vs. small coins) • the structure of the problem (coin values). 3.2 Graphs We assume that you have seen some graph theory before. We will give a quick review here; if you are shaky on graph theory then you should review §B.4 of CLRS. A graph is a mathematical object which consists of vertices and edges. Thus, we may say that a graph is a pair G = (V , E), where V is a set of vertices and E is a set of edges. Or, if the vertex set V has been fixed ahead of time, we may simply say that “E is a graph on V .” In an undirected graph, the edges are unordered pairs e = { u, v} with u, v ∈ V ; whereas in a directed graph, the edges are ordered pairs e = (u, v). A directed edge e = (u, v) is said to emanate “from u to v,” and we may use the notation e = (u → v).1 Our convention (which is by no means standard convention, since there is no standard) will be that all graphs are assumed finite, and are allowed to have loops (edges u → u) but may not have multiple edges between the same pair of vertices. Thus E is just a set, not a multiset (though in a directed graph, u → v and v → u count as different edges). It is often useful to consider graphs with satellite data attached to the edges or vertices (or both). One common instance of this is a weighted undirected graph, which is an undirected graph equipped with a weight function w : E → R. We say that the edge e = { u, v} has weight w(e), or sometimes wuv . For example, let V be the set of cities in the United States, and let E be a complete graph on V , meaning that E contains every pair of vertices { u, v} with u 6= v (but no loops). Let w(u, v) be the distance in miles between u and v. The graph G = (V , E, w) is obviously of great importance to airline booking companies, shipping companies, and presidential candidates in the last few months of compaigning. Consider a fixed set V of vertices. A tree T on V is a graph which is acyclic and connected.2 (See Figure 3.2.) Note that acyclic is a smallness condition, while connectivity is a largeness condition. Thus, trees are “just the right size.” According to the following proposition, that size is |V | − 1. 1 In fact, we will probably use several other notations, depending on our whim; in particular, we will sometimes use the notation of directed graphs in discussions of undirected graphs. It is the author’s job to keep clear the meaning of the notation—you should speak up if you are confused by any notation, since it is more likely the author’s fault than yours. 2 A path in a graph (V , E) is a sequence of vertices 〈v , v , . . . , v 〉, where (v → v n 0 1 j j +1 ) ∈ E for j = 0, . . . , n − 1. The number n is called the length of the path. If v0 = vn , then the path is called a cycle. A graph is called acyclic if it contains no cycles, and is said to be connected if, for every pair u, v of vertices, there exists a path from u to v. Lec 3 – pg. 2 of 7 Figure 3.2. Left: An undirected graph G. Center: Two different spanning trees of G. Right: A subgraph of G which is a tree but not a spanning tree. The gray vertices and edges are just visual aids; they are not part of the graphs. Proposition 3.1 (Theorem B.2 of CLRS, slightly modified). Let V be a finite set of vertices, and let G = (V , T) be an undirected graph on V . The following conditions are equivalent: (i) (ii) (iii) (iv) (v) (vi) G is acyclic and connected. G is acyclic and has at least |V | − 1 edges. G is connected and has at most |V | − 1 edges. G is acyclic, but if any edge is added to T, the resulting graph contains a cycle. G is connected, but if any edge is removed from T, the resulting graph is not connected. For any two vertices u, v ∈ V , there exists a unique simple path from u to v.3 Any one of these six conditions could be taken as the definition of a tree. If G is a tree, then G has exactly |V | − 1 edges. We will not prove the proposition here, but you would do yourself a service to write the proof yourself (or look it up in CLRS), especially if you have not seen proofs in graph theory before. There are certain types of arguments that occur again and again in graph theory, and you will definitely want to know them for the remainder of the course. Definition. Let G = (V , E) be an undirected graph. A spanning tree of G is a subset T ⊆ E such that T is a tree on V . The reason for the word “spanning” is that T must be a tree on all of V . There are plenty of subgraphs of G that are trees but not spanning trees: they are graphs of the form G 0 = (V 0 , T 0 ) where V 0 á V and T 0 is a tree on V 0 (but not on V because it does not touch all the vertices). Proposition 3.2. Let G = (V , E) be an undirected graph. Then G has a spanning tree if and only if G is connected. Thus, every connected graph on V has at least |V | − 1 edges. Proof. If G has a spanning tree T, then for any pair of vertices u, v there is a path from u to v in T. That path also lies in E since T ⊆ E. Thus G is connected. Conversely, suppose G is connected but is not a tree. Then, by Proposition 3.1(v), we can remove an edge from E, obtaining a connected subgraph E 1 . If E 1 is not a tree, then we can repeat the process again. Eventually we will come upon a connected subgraph (V , E k ) of G such that it is impossible to remove an edge from E k , either because we have run out of edges, or because the resulting graph will be disconnected. In the former case, V must have only one vertex, so the empty set of edges is a spanning tree. In the latter case, E k is a spanning tree by Proposition 3.1(v). It follows from Proposition 3.1 that G has at least |V |− 1 edges, since any spanning tree has exactly |V |− 1 edges. 3 A path is said to be simple if all vertices in the path are distinct. A cycle is said to be simple if all vertices are distinct except for the common start/end vertex. Thus, strictly speaking, a simple cycle is not a simple path; extra clarifications will be made in any situation where there might be confusion. Lec 3 – pg. 3 of 7 Figure 3.3. Clustering using an MST. 3.3 Minimum Spanning Trees Given a weighted undirected graph G = (V , E, w), one often wants to find a minimum spanning P tree (MST) of G: a spanning tree T for which the total weight w(T) = (u,v)∈T w(u, v) is minimal. Input: A connected, undirected weighted graph G = (V , E, w) Output: A spanning tree T such that the total weight X w(T) = w(u, v) (u,v)∈T is minimal. For example, if V is the set of buildings in a city and E is the complete graph on V , and if w(u, v) is the distance between u and v, then a minimum spanning tree would be very useful in building a minimum-length fiber-optic network, pipeline network, or other infrastructure for the city. There are many other, less obvious applications of minimum spanning trees. Given a set of data points on which we have defined a metric (i.e., some way of quantifying how close together or similar a pair of vertices are), we can use an MST to cluster the data by starting with an MST for the distance graph and then deleting the heaviest edges. (See Figure 3.3 and Lecture 21.) If V is a set of data fields and the distance is mutual information, then this clustering can be used to find higher-order correlative relationships in a large data set. How would we go about finding an MST for this graph? 6 12 5 9 14 7 8 15 3 10 There are several sensible heuristics: • Avoid large weights. We would like to avoid the 14 and 15 edges if at all possible. Lec 3 – pg. 4 of 7 U u u0 v v0 Figure 3.4. Illustration of the MST property. • Include small weights. If the 3 and 5 edges provide any benefit to us, we will probably want to include them. • Some edges are inevitable. The 9 and 15 edges must be included in order for the graph to be connected. • Avoid cycles, since an MST is not allowed to have cycles. Some more thoughts: • Would a greedy algorithm be likely to work here? Is there something about the structure of the problem which allows locally optimal decisions to be informed enough about the global solution? • Should we start with all the edges present and then remove some of them, or should we start with no edges present and then add them? • Is the MST unique? Or are there multiple solutions? A key observation to make is the following: If we add an edge to an MST, the resulting graph will have a cycle, by Proposition 3.1(iv). If we then remove one of the edges in this cycle, the resulting graph will be connected and therefore (by Proposition 3.1(iii)) will be a spanning tree. If we have some way of knowing that the edge we removed is at least as heavy as the edge we added, then it follows that we will have a minimum spanning tree. This is the idea underlying the so-called “MST property.” Theorem 3.3 (MST Property). Let G = (V , E, w) be a connected, weighted, undirected graph. Let U be a proper nonempty subset of V .4 Let S be the set of edges (x, y) with x ∈ U and y ∈ V \ U. Suppose (u, v) is the lightest edge (or one of the lightest edges) in S. Then there exists an MST containing (u, v). This partition of V into U and V \ U is called a “cut.” An edge is said to “cross” the cut if one endpoint is in U and the other endpoint is in V \ U. Otherwise the edge is said to “respect” the cut. So S is the set of edges which cross the cut and E \ S is the set of edges that respect the cut. An edge in S of minimal weight is called a “light edge” for the cut. 4 That is, ; á U á V . Lec 3 – pg. 5 of 7 Proof. Let T ⊆ E be an MST for G, and suppose (u, v) 6∈ T. If we add the edge (u, v) to T, the resulting graph T 0 has a cycle. This cycle must cross the cut (U, V \ U) in an even number of places, so in addition to (u, v), there must be some other edge (u0 , v0 ) in the cycle, such that (u0 , v0 ) crosses the cut. Remove (u0 , v0 ) from T 0 and call the resulting graph T 00 . We showed above that T 00 is a spanning tree. Also, since (u, v) is a light edge, we have w(u, v) ≤ w(u0 , v0 ). Thus w(T 00 ) = w(T) + w(u, v) − w(u0 , v0 ) ≤ w(T). Since T is a minimum spanning tree and T 00 is a spanning tree, we also have w(T) ≤ w(T 00 ). Therefore w(T) = w(T 00 ) and T 00 is a minimum spanning tree which contains (u, v). Corollary 3.4. Preserve the setup of the MST property. Let T be any MST. Then there exists an edge ¡ ¢ (u0 , v0 ) ∈ T such that u0 ∈ U and v0 ∈ V \ U and T \ {(u0 , v0 )} ∪ {(u, v)} is an MST. ¡ ¢ Proof. In the proof of the MST property, T 00 = T \ {(u0 , v0 )} ∪ {(u, v)}. Corollary 3.5. Let G = (V , E, w) be a connected, weighted, undirected graph. Let T be any MST and let (U, V \ U) be any cut. Then T contains a light edge for the cut. If the edge weights of G are distinct, then G has a unique MST. Proof. If T does not contain a light edge for the cut, then the graph T 00 constructed in the proof of the MST property weighs strictly less than T, which is impossible. Suppose the edge weights in G are distinct. Let M be an MST. For each edge (u, v) ∈ M, consider the graph (V , M \ {(u, v)}). It has two connected components5 ; let U be one of them. The only edge in M that crosses the cut (U, V \ U) is (u, v). Since M must contain a light edge, it follows that (u, v) is a light edge for this cut. Since G has distinct edge weights, (u, v) is the unique light edge for the cut, and every MST must contain (u, v). Letting (u, v) vary over all edges of M, we find that every MST must contain M. Thus, every MST must equal M. Now do you imagine that a greedy algorithm might work? The MST property is the crucial idea that allows us to use local information to decide which edges belong in an MST. Below is Kruskal’s algorithm, which solves the MST problem. Algorithm: K RUSKAL -MST(V , E, w) Initially, let T ← ; be the empty graph on V . Examine the edges in E in increasing order of weight (break ties arbitrarily). • If an edge connects two unconnected components of T, then add the edge to T.∗ • Else, discard the edge and continue. Terminate when there is only one connected component. (Or, you can continue through all the edges.) ∗ By “two unconnected components,” we mean two distinct connected components. Please forgive me for this egregious abuse of language. 5 Given a graph G, we say a vertex v is reachable from u if there exists a path from u to v. For an undirected graph, reachability is an equivalence relation on the vertices. The restrictions of G to each equivalence class are called the connected components of G. (The restriction of a graph G = (V , E) to a subset V 0 ⊆ V is the subgraph (V 0 , E 0 ), where E 0 is the set of edges whose endpoints are both in V 0 .) Thus each connected component is a connected subgraph, and a connected undirected graph has only one connected component. Lec 3 – pg. 6 of 7 Try performing Kruskal’s algorithm by hand on our example graph. Proof of correctness for Kruskal’s algorithm. We will use the following loop invariant: Prior to each iteration, every edge in T is a subset of an MST. • Initialization. T has no edges, so of course it is a subset of an MST. • Maintenance. Suppose T ⊆ T ∗ where T ∗ is an MST, and suppose the edge (u, v) gets added to T. Then, as we will show in the next paragraph, (u, v) is a light edge for the cut (U, V \U), where U is one of the connected components of T. Therefore, by Corollary 3.4, there exist vertices u0 ∈ U ¡ ¢ and v0 ∈ V \ U with (u0 , v0 ) ∈ T ∗ such that T 0 is an MST, where T 0 = T ∗ \ {(u0 , v0 )} ∪ {(u, v)}. But (u0 , v0 ) 6∈ T, so T ∪ {(u, v)} is a subset of T 0 . As promised, we show that (u, v) is a light edge for the cut (U, V \ U), where U is one of the connected components of T. The reasoning is this: If (u, v) is added, then it joins two unconnected components of T. Let U be one of those components. Now suppose that (x, y) ∈ E were some lighter edge crossing the cut (U, V \ U), say x ∈ U and y ∈ V \ U. Since x and y lie in different connected components of T, the edge (x, y) is not in T. Thus, when we examined (x, y), we decided not to add it to T. This means that, when we examined the edge (x, y), x and y lay in the same connected component of Tprev , where Tprev is what T used to be at that point in time. But Tprev ⊆ T, which means that x and y lie in the same connected component of T. This contradicts the fact that x ∈ U and y ∈ V \ U. • Termination. At termination, T must be connected: If T were disconnected, then there would exist an edge of E which joins two components of T. But that edge would have been added to T when we examined it, a contradiction. Now, since T is connected and is a subset of an MST, it follows that T is an MST. To implement Kruskal’s algorithm and analyze its running time, we will first need to think about the data structures that will be used and the running time of operations on them. Exercise 3.1. What if we wanted a maximum spanning tree rather than a minimum spanning tree? Lec 3 – pg. 7 of 7 Lecture 4 Minimum Spanning Trees II Supplemental reading in CLRS: Chapter 23; Section 16.2; Chapters 6 and 21 4.1 Implementing Kruskal’s Algorithm In the previous lecture, we outlined Kruskal’s algorithm for finding an MST in a connected, weighted undirected graph G = (V , E, w): Initially, let T ← ; be the empty graph on V . Examine the edges in E in increasing order of weight (break ties arbitrarily). • If an edge connects two unconnected components of T, then add the edge to T. • Else, discard the edge and continue. Terminate when there is only one connected component. (Or, you can continue through all the edges.) Before we can write a pseudocode implementation of the algorithm, we will need to think about the data structures involved. When building up the subgraph T, we need to somehow keep track of the connected components of T. For our purposes it suffices to know which vertices are in each connected component, so the relevant information is a partition of V . Each time a new edge is added to T, two of the connected components merge. What we need is a disjoint-set data structure. 6 3 12 8 9 2 5 6 14 4 9 8 5 10 7 1 3 7 15 10 Figure 4.1. Illustration of Kruskal’s algorithm. 4.1.1 Disjoint-Set Data Structure A disjoint-set data structure maintains a dynamic collection of pairwise disjoint sets S = {S 1 , . . . , S r } in which each set S i has one representative element, rep[S i ]. Its supported operations are • M AKE -S ET(u): Create new set containing the single element u. – u must not belong to any already existing set – of course, u will be the representative element initially • F IND -S ET(u): Return the representative rep[S u ] of the set S u containing u. • U NION(u, v): Replace S u and S v with S u ∪ S v in S. Update the representative element. 4.1.2 Implementation of Kruskal’s Algorithm Equipped with a disjoint set data structure, we can implement Kruskal’s algorithm as follows: Algorithm: K RUSKAL -MST(V , E, w) 1 B Initialization and setup 2 T ←; for each vertex v ∈ V do M AKE -S ET(v) Sort the edges in E into non-decreasing order of weight B Main loop for each edge (u, v) ∈ E in non-decreasing order of weight do if F IND -S ET(u) 6= F IND -S ET(v) then T ← T ∪ {(u, v)} U NION(u, v) return T 3 4 5 6 7 8 9 10 11 The running time of this algorithm depends on the implementation of the disjoint set data structure we use. If the disjoint set operations have running times T M AKE -S ET , T U NION and T F IND -S ET , and if we use a good O(n lg n) sorting algorithm to sort E, then the running time is O(1) + V T M AKE -S ET + O(E lg E) + 2ET F IND -S ET + O(E) + ET U NION .1 4.1.3 Implementations of Disjoint-Set Data Structure The two most common implementations of the disjoint-set data structure are (1) a collection of doubly linked lists and (2) a forest of balanced trees. In what follows, n denotes the total number of elements, i.e., n = |S 1 | + · · · + |S r |. Solution 1: Doubly-linked lists. Represent each set S i as a doubly-linked list, where each element is equipped with a pointer to its two neighbors, except for the leftmost element which has a “stop” marker on the left and the rightmost element which has a “stop” marker on the right. We’ll take the leftmost element of S i as its representative. 1 Actually, we can do better. In line 10, since we have already computed rep[S ] and rep[S ], we do not need to call u v U NION; we need only call W EAK -U NION, an operation which merges two sets assuming that it has been given the correct representative of each set. So we can replace the ET U NION term with ET W EAK -U NION . Lec 4 – pg. 2 of 7 Si : rep[S i ] u1 u2 u3 ··· uk M AKE -S ET(u) – initialize as a lone node Θ(1) F IND -S ET(u) – walk left from u until you reach the head Θ(n) worst-case U NION(u, v) Θ(n) worst-case – walk right from u and left from v. Reassign pointers so that the tail of S u and the head of S v become neighbors. The representative is updated automatically. These can be improved upon—there exist better doubly-linked list implementations of the disjoint set data structure. Solution 2: Forest of balanced trees.2 u1 Si : rep[S i ] u4 u2 u5 u3 u6 M AKE -S ET(u) – initialize new tree with root node u Θ(1) F IND -S ET(u) – walk up tree from u to root Θ(height) = Θ(lg n) best-case U NION(u, v) O(1) + 2T F IND -S ET – change rep[S v ]’s parent to rep[S u ] The forest of balanced trees will be our implementation of choice. With a couple of clever tricks3 , the running times of the operations can be greatly improved: In the worst case, the improved structure has an amortized (average) running time of Θ(α(n)) per operation4 , where α(n) is the inverse Ackermann function, which is technically unbounded but for all practical purposes should be considered bounded.5 So in essence, each disjoint set operation takes constant time, on average. 2 A rooted tree is a tree with one distinguished vertex u, called the root. By Proposition 3.1(vi), for each vertex v there exists a unique simple path from u to v. The length of that path is called the depth of v. It is common to draw rooted trees with the vertices arranged in rows, with the root on top, all vertices of depth 1 on the row below that, etc. 3 The tricks are called union-by-rank and path compression. For more information, see Lecture 16. 4 In a 1989 paper, Fredman and Saks proved that Θ(α(n)) is the optimal amortized running time. 5 For example, if n is small enough that it could be written down by a collective effort of the entire human population before the sun became a red giant star and swallowed the earth, then α(n) ≤ 4. Lec 4 – pg. 3 of 7 In the analysis that follows, we will not use these optimizations. Instead, we will assume that F IND -S ET and U NION both run in Θ(lg n) time. The asymptotic running time of K RUSKAL -MST is not affected. As we saw above, the running time of K RUSKAL -MST is O(1) Initialize: z }| { O(1) + V T M AKE -S ET +O(E lg E) Loop: 2E T F IND -S ET +O(E) + E T U NION | {z } | {z } O(lg V ) O(lg V ) O(E lg E) + 2O(E lg V ). Since there can only be at most V 2 edges, we have lg E ≤ 2 lg V . Thus the running time of Kruskal’s algorithm is O(E lg V ), the same amount of time it would take just to sort the edges. 4.1.4 Safe Choices Let’s philosophize about Kruskal’s algorithm a bit. When adding edges to T, we do not worry about whether T is connected until the end. Instead, we worry about making “safe choices.” A safe choice is a greedy choice which, in addition to being locally optimal, is also part of some globally optimal solution. In our case, we took great care to make sure that at every stage, there existed some MST T ∗ such that T ⊆ T ∗ . If T is safe and T ∪{(u, v)} is also safe, then we call (u, v) a “safe edge” for T. We have already done the heavy lifting with regard to safe edge choices; the following theorem serves as a partial recap. Proposition 4.1 (CLRS Theorem 23.1). Let G = (V , E, w) be a connected, weighted, undirected graph. Suppose A is a subset of some MST T. Suppose (U, V \ U) is a cut of G that is respected by A, and that (u, v) is a light edge for this cut. Then (u, v) is a safe edge for A. Proof. In the notation of Corollary 3.4, the edge (u0 , v0 ) does not lie in A because A respects the cut ¡ ¢ (U, V \ U). Therefore A ∪ {(u, v)} is a subset of the MST T \ {(u0 , v0 )} ∪ {(u, v)}. 4.2 Prim’s Algorithm We now present a second MST algorithm: Prim’s algorithm. Like Kruskal’s algorithm, Prim’s algorithm depends on a method of determining which greedy choices are safe. The method is to continually enlarge a single connected component by adjoining edges emanating from isolated vertices.6 Algorithm: P RIM -MST(V , E, w) 1 2 3 4 5 6 7 8 Choose an arbitrary start vertex s C ← { s} T ←; while C is not the only connected component of T do Select a light edge (u, v) connecting C to an isolated vertex v T ← T ∪ {(u, v)} C ← C ∪ { v} return T 6 An isolated vertex is a vertex which is not connected to any other vertices. Thus, an isolated vertex is the only vertex in its connected component. Lec 4 – pg. 4 of 7 Proof of correctness for Prim’s algorithm. Again, we use a loop invariant: Prior to each iteration, T is a subset of an MST. • Initialization. T has no edges, so trivially it is a subset of an MST. • Maintenance. Suppose T ⊆ T ∗ where T ∗ is an MST, and suppose the edge (u, v) gets added to T, where u ∈ C and v is an isolated vertex. Since (u, v) is a light edge for the cut (C, V \ C) which is respected by T, it follows by Proposition 4.1 that (u, v) is a safe edge for T. Thus T ∪ {(u, v)} is a subset of an MST. • Termination. At termination, C is the only connected component of T, so by Proposition 3.1(v), T has at least |V | − 1 edges. Since T is also a subset of an MST, it follows that T has exactly |V | − 1 edges and is an MST. The tricky part of Kruskal’s algorithm was keeping track of the connected components of T. In Prim’s algorithm this is easy: except for the special component C, all components are isolated vertices. The tricky part of Prim’s algorithm is efficiently keeping track of which edge is lightest among those which join C to a new isolated vertex. This task is typically accomplished with a data structure called a min-priority queue. 4.2.1 Min-Priority Queue A min-priority queue is a data structure representing a collection of elements which supports the following operations: I NSERT(Q, x) – Inserts element x into the set of elements Q M INIMUM(Q) – Returns the element of Q with the smallest key E XTRACT-M IN(Q) – Removes and returns the element with the smallest key D ECREASE -K EY(Q, x, k) – Decreases the value of x’s key to new value k. With this data structure, Prim’s algorithm could be implemented as follows: Lec 4 – pg. 5 of 7 Algorithm: P RIM -MST(G, s) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 T ←; for each u ∈ G.V do u.key ← ∞ B initialize all edges to “very heavy” B The component C will be a tree rooted at s. Once a vertex u gets added to C, u.π will be a pointer to its parent in the tree. u.π ← N IL s.key ← 0 B this ensures that s will be the first vertex we pick Let Q ← G.V be a min-priority queue while Q 6= ; do u ← E XTRACT-M IN(Q) if u.π 6= N IL then T ← T ∪ {(u, u.π)} B We assume that G is presented in the adjacency-list format∗ for each v ∈ G.adj[u] do if v ∈ Q and w(u, v) < v.key then v.π ← u v.key ← w(u, v) B using D ECREASE -K EY return T ∗ For more information about ways to represent graphs on computers, see §22.1 of CLRS. 4.2.2 Running Time of Prim’s Algorithm Lines 1 through 6 clearly take O(V ) time. Line 7 takes T B UILD -Q UEUE (V ), where T B UILD -Q UEUE (n) is the amount of time required to build a min-priority queue from an array of n elements. Within the “while” loop of line 8, E XTRACT-M IN gets called |V | times, and the instructions in lines 14–16 are run a total of O(E) times. Thus the running time of Prim’s algorithm is O(V ) + T B UILD -Q UEUE (V ) + V T E XTRACT-M IN + O(E)T D ECREASE -K EY . Exercise 4.1. How can we structure our implementation so that line 14 runs in O(1) time? Let’s take a look at the performance of P RIM -MST under some implementations of the minpriority queue, in increasing order of efficiency: Q Array Binary Heap Fibonacci Heap T B UILD -Q UEUE (n) O(n) O(n) O(n) T E XTRACT-M IN O(n) O(lg n) O(lg n) | {z } amortized T D ECREASE -K EY O(1) O(lg n) O(1) | {z } Runtime of P RIM -MST O(V 2 ) O(E lg V ) O(E + V lg V ) amortized Thus, for a dense graph (a graph with Θ(V 2 ) edges), Prim’s algorithm with a Fibonacci heap outperforms Kruskal’s algorithm. The best known MST algorithm to date is a randomized algorithm with Θ(V + E) expected running time, introduced by Karger, Klein and Tarjan in 1995. Lec 4 – pg. 6 of 7 4.3 Greedy Strategies General approach: 1. Structure the problem so that we make a choice and are left with one subproblem to solve. 2. Make a greedy choice and then prove that there exists an optimal solution to the original problem which makes the same greedy choice (“safe choice”). 3. Demonstrate optimal substructure. • After making the greedy choice, combine with the optimal solution of the remaining subproblem, giving an optimal solution to the original problem. Note that this sounds a lot like dynamic programming. Let’s now examine the key properties of greedy algorithms in comparison to those of dynamic programming. 1) Greedy Choice Property – Locally optimal solution leads to globally optimal solution. Each greedy local choice is made independently of the solution to the subproblem. (E.g.: Kruskal’s and Prim’s algorithms can choose a safe edge without having examined the full problem.) In dynamic programming, the local choice depends on the solution to the subproblem—it is a bottom-up solution. 2) Optimal Substructure – The optimal solution to a problem contains optimal solutions to subproblems. Both greedy strategies and dynamic programming exploit (or rely on) optimal substructures. Prim’s algorithm produces (optimal) MSTs on subsets of V on the way to finding the full MST. Lec 4 – pg. 7 of 7 Lecture 5 Fast Fourier Transform Supplemental reading in CLRS: Chapter 30 The algorithm in this lecture, known since the time of Gauss but popularized mainly by Cooley and Tukey in the 1960s, is an example of the divide-and-conquer paradigm. Actually, the main uses of the fast Fourier transform are much more ingenious than an ordinary divide-and-conquer strategy— there is genuinely novel mathematics happening in the background. Ultimately, the FFT will allow us to do n computations, each of which would take Ω(n) time individually, in a total of Θ(n lg n) time. 5.1 Multiplication To motivate the fast Fourier transform, let’s start with a very basic question: How can we efficiently multiply two large numbers or polynomials? As you probably learned in high school, one can use essentially the same method for both: 385 ×426 2310 770 +1540 164010 (3x2 + 8x + 5) × (4x2 + 2x + 6) 18x2 + 48x + 30 6x3 + 16x2 + 10x 4 12x + 32x3 + 20x2 12x4 + 38x3 + 54x2 + 58x + 30 Of these two, polynomials are actually the easier ones to work with. One reason for this is that multiplication of integers requires carrying, while multiplication of polynomials does not. To make a full analysis of this extra cost, we would have to study the details of how large integers can be stored and manipulated in memory, which is somewhat complicated.1 So instead, let’s just consider multiplication of polynomials. Suppose we are trying to multiply two polynomials p, q of degree at most n with complex coefficients. In the high-school multiplication algorithm (see Figure 5.1), each row of the diagram is 1 More ideas are required to implement efficient multiplication of n-bit integers. In a 1971 paper, Schönhage and ³ ´ ∗ Strassen exhibited integer multiplication in O(n lg n · lg lg n) time; in 2007, Fürer exhibited O n lg n · 2O(lg n) time. The iterated logarithm lg∗ n is the smallest k such that lg lg · · · lg n (k times) ≤ 1. It is not known whether integer multiplication can be achieved in O(n lg n) time, whereas by the end of this lecture we will achieve multiplication of degree-n polynomials in O(n lg n) time. n+1 ◦ ◦ ◦ ◦ ◦ ◦ ··· ◦ × ◦ ◦ ··· ◦ ◦ ◦ ◦ ◦ ··· ··· ··· ◦ ◦ n+1 ◦ ··· ◦ ◦ 2n + 1 Figure 5.1. The high-school polynomial multiplication algorithm. obtained by multiplying p by a monomial. Thus any single row would take Ω(n) time to compute individually. There are n + 1 rows, so the computation takes Ω(n2 ) time. The key to improving this time is to consider alternative ways of representing p and q. Fix some N > n. We’ll choose N later, but for now all that matters is that N = O(n). A polynomial of degree at most N − 1 is uniquely determined by its values at N points. So, instead of storing the coefficients of p and q, we could represent p and q as the lists ­ ® ­ ® p z = p(z0 ), . . . , p(z N −1 ) and q z p(z0 ), . . . , p(z N −1 ) for any distinct complex numbers z0 , . . . , z N −1 . In this representation, computing pq is very cheap— the list ­ ® (pq) z = p(z0 )q(z0 ), . . . , p(z N −1 )q(z N −1 ) can be computed in O(n) time. If N > deg pq, then pq is the unique polynomial whose FFT is (pq) z . Algorithm: M ULTIPLY(p, q) 1. Fix some N = O(deg p + deg q) such that N > deg p + deg q. 2. Compute the sample vectors p z = 〈 p(z0 ), . . . , p(z N −1 )〉 and q z = 〈 q(z0 ), . . . , q(z N −1 )〉. 3. Compute the unique polynomial r such that r z = p z · q z , where · denotes entrywise multiplication.∗ 4. Return r. ∗ For example, 〈1, 2, 3〉 · 〈4, 5, 6〉 = 〈4, 10, 18〉. Thus, the cost of finding the coefficients of pq is equal to the cost of converting back and forth between the coefficients representation and the sample-values representation (see Figure 5.2). We are free to choose values of z0 , . . . , z N −1 that make this conversion as efficient as possible. 5.1.1 Efficiently Computing the FFT −1 In the fast Fourier transform, we choose z0 , . . . , z N −1 to be the Nth roots of unity, 1, ω N , ω2N , . . . , ω N , N 2π i/N where ω N = e . We make this choice because roots of unity enjoy the useful property that ω2N = ω N/2 . Lec 5 – pg. 2 of 6 PN × PN FFT / CN × CN Entrywise mult. (fast) Polynomial mult. (slow)  PN o FFT−1  CN Figure 5.2. Commutative diagram showing the cost of multiplication on either side of a fast Fourier transform. As we will see, the fastest way to get from the top-left to the bottom-left is through the FFT. (We can choose N to be a power of 2.) This allows us to divide the task of computing an FFT of size N into computing two FFTs of size N/2, and then combining them in a certain way.2 ¡ ¢ Algorithm: FFT N, p(x) = a N −1 x N −1 + a N −2 x N −2 + · · · + a 1 x + a 0 1 B As part of the specification, we assume N is a power of 2 if N = 1 then return a 0 4 Let 2 3 peven (y) ← a 0 + a 2 y + · · · + a N −2 y N/2−1 , podd (y) ← a 1 + a 3 y + · · · + a N −1 y N/2−1 be the even and odd parts of p 5 B Thus ¡ ¢ ¡ ¢ p(x) = peven x2 + xpodd x2 (5.1) even peven ) ω N/2 ← FFT ¡(N/2, p ¢ odd odd 7 p ωN/2 ← FFT N/2, p 8 B That is, 6 D ³ ´E even even even N/2−1 peven = p (1), p ω , . . . , p ω ( ) N/2 ω N/2 N/2 D ³ ´E odd odd odd odd N/2−1 p ωN/2 = p (1), p (ω N/2 ) , . . . , p ω N/2 . 9 −1 B Because ω N/2 = ω2N , we can calculate the vector p ωN = p(1), p (ω N ) , . . . , p ω N N ­ ¡ ¢® very quickly using (5.1): ­ ® −1 10 ω ← 1, ω N , . . . , ω N B the left side is a bold omega N odd 11 return peven + ω · p , where · denotes entrywise multiplication ω N/2 ω N/2 2 I suppose the “conquer” stage is when we recursively compute the smaller FFTs (but of course, each of these smaller FFTs begins with its own “divide” stage, and so on). After the “conquer” stage, the answers to the smaller problems are combined into a solution to the original problem. Lec 5 – pg. 3 of 6 even Above, we compute p ωN by computing podd ω N/2 and p ω N/2 and combining them in Θ(N) time. Thus, the running time of an FFT of size N satisfies the recurrence T(N) = 2T(N/2) + Θ(N). This recurrence is solved in CLRS as part of the Master Theorem in §4.5. The solution is T(N) = Θ(N lg N). 5.1.2 Computing the Inverse FFT Somewhat surprisingly, the inverse FFT can be computed in almost exactly the same way as the FFT. In this section we will see the relation between the two transforms. If you don’t have a background in linear algebra, you can take the math on faith. Let P N be the vector space of polynomials of degree at most N − 1 with complex coefficients. Then the FFT is a bijective linear map P N → C N . If we use the ordered basis 1, x, . . . , x N −1 for P N and the standard basis for C N , then the matrix of the FFT is   1 1 1 ··· 1 1 ω −1  ³ ´ ω2N · · · ωN N   N ij .  A = ωN = .. .. .. .. ..  0≤ i, j ≤ N −1  . . . .  . 2 (N −1) 2(N −1) N −1 · · · ωN 1 ωN ωN j The jth column of A −1 contains the coefficients of the polynomial which, when evaluated on ω N , i gives 1 and, when evaluated on ω N (i 6= j), gives zero. This polynomial is Q Q i 6= j ³ i 6= j ¡ i x − ωN j ¢ i ωN − ωN ´. We will not show the details here, but after about a page of calculation you can find that A −1 = 1 N  1 1 1  . A= N   .. 1 1 ω− N .. . (N −1) 1 ω− N 1 2 ω− N .. . ··· ··· .. . 2(N −1) ω− N ···  1 (N −1)  ω−  N , ..  .  2 (N −1) ω− N where the bar indicates entrywise complex conjugation. You can certainly check that the i, j entry of 1 N A A is ( NX −1 NX −1 1, i = j `(i − j) i` 1 − j` 1 ωN · N ωN = N ωN = 0, i 6= j, `=0 `=0 and similarly for 1 N A A. Thus FFT−1 (v) = 1 N FFT ¡ ¢ v. (5.2) Again, this is surprising. A priori, we would consider P N and C N to be “different worlds”—like two different implementations of the same abstract data structure, the coefficients representation and Lec 5 – pg. 4 of 6 the sample-values representation of polynomials serve essentially the same role but have different procedures for performing operations. Yet, (5.2) tells us that inverting an FFT involves creating a new polynomial whose coefficients are the sample values of our original polynomial. This of course has everything to do with the fact that roots of unity are special; it would not have worked if we had −1 as our sample points. not chosen 1, ω N , . . . , ω N N ¡ ¢ Algorithm: I NVERSE -FFT N, v = 〈v0 , . . . , v N −1 〉 ­ ® 1 v ← v0 , . . . , v N −1 2 Let p v (x) be the polynomial v N −1 x N −1 + · · · + v1 x + v0 3 return the value of ¡ ¢ 1 N FFT N, p v With all the details in place, the M ULTIPLY algorithm looks like this: Algorithm: FFT-M ULTIPLY(p, q) 1. Let N be the smallest power of 2 such that N − 1 ≥ deg p + deg q. 2. Compute v = FFT(N, p) and w = FFT(N, q), and let u = v · w, where · denotes entrywise multiplication. 3. Compute and return I NVERSE -FFT(N, u). ¡ ¢ It runs in Θ(N lg N) = Θ (deg p + deg q) lg(deg p + deg q) time. 5.2 Convolution Let p(x) = a n−1 x n−1 + · · · + a 1 x + a 0 and q(x) = b n−1 x n−1 + · · · + b 1 x + b 0 . The coefficient of x k in pq is X a ` b k−` , `∈Z2n−1 where the subscripts are interpreted modulo 2n − 1. This sort of operation is known as convolution, and is written as ∗. Convolution can be applied quite generally: if f and g are any functions Z2n−1 → C, then one common definition of the convolution operation is ( f ∗ g)(k) = X 1 f (`)g(k − `). 2n − 1 `∈Z2n−1 For each i, computing ( f ∗ g)(i) requires linear time, but the fast Fourier transform allows us to compute ( f ∗ g)(i) for all i = 0, . . . , 2n − 2 in a total of Θ(n lg n) time. Convolution appears frequently, which is part of the reason that the FFT is useful. The above notion of convolution can easily be generalized to allow f and g to be functions from any group G to any ring in which |G | is a unit. There is also a continuous version of convolution which involves an integral. This continuous version is useful in signal compression. The idea is to let f : R → R be an arbitrary signal and let g be a smooth bump function which spikes up near the origin and is zero elsewhere. (See Figure 5.3.) Then the convolution of f and g is a new function f ∗ g defined by Z ∞ ( f ∗ g)(t) = f (τ)g(t − τ) dτ. −∞ Lec 5 – pg. 5 of 6 Noisy signal f Bump function g Convolution f ∗ g Figure 5.3. A complicated, noisy signal f , a bump function g, and the convolution f ∗ g. The value ( f ∗ g)(t) is essentially an average of the values of f at points close to t. Thus, the funtion f ∗ g inherits smoothness from g while still carrying most of the information from f . As a bonus, if f comes with some unwanted random noise, then g will have much less noise. For example, f may be the raw output of a recording device, and computing f ∗ g may be the first step towards encoding f into a compressed audio format such as MP3. The point is that smooth functions are much easier to describe concisely, so creating a smooth version of f is useful in (lossy) data compression. Lec 5 – pg. 6 of 6 Lecture 6 All-Pairs Shortest Paths I Supplemental reading in CLRS: Chapter 25 intro; Section 25.2 6.1 Dynamic Programming Like the greedy and divide-and-conquer paradigms, dynamic programming is an algorithmic paradigm in which one solves a problem by combining the solutions to smaller subproblems. In dynamic programming, the subproblems overlap, and the solutions to “inner” problems are stored in memory (see Figure 6.1). This avoids the work of repeatedly solving the innermost problem. Dynamic programming is often used in optimization problems (e.g., finding the maximum or minimum solution to a problem). • The key feature that a problem must have in order to be amenable to dynamic programming is that of optimal substructure: the optimal solution to the problem must contain optimal solutions to subproblems. • Greedy algorithms are similar to dynamic programming, except that in greedy algorithms, the solution to an inner subproblem does not affect the way in which that solution is augmented to the solution of the full problem. In dynamic programming the combining process is more sophisticated, and depends on the solution to the inner problem. • In divide-and-conquer, the subproblems are disjoint. Dynamic Programming Divide-and-Conquer problem subproblem 1 subproblem subsubproblem subproblem 2 subproblem 3 problem Figure 6.1. Schematic diagrams for dynammic programming and divide-and-conquer. 6.2 Shortest-Path Problems Given a directed graph G = (V , E, w) with real-valued edge weights, it is very natural to ask what the shortest (i.e., lightest) path between two vertices is. Input: A weighted, directed graph G = (V , E, w) with real-valued edge weights, a start vertex u, and an end vertex v Output: The minimal weight of a path from u to v. That is, the value of n o ( p min w(p) : paths u v if such a path exists δ(u, v) = ∞ otherwise, where the weight of a path p = 〈v0 , v1 , . . . , vk 〉 is w(p) = k X w(v i−1 , v i ). i =1 (Often we will also want an example of a path which achieves this minimal weight.) Shortest paths exhibit an optimal-substructure property: Proposition 6.1 (Theorem 24.1 of CLRS). Let G = (V , E, w) be a weighted, directed graph with realvalued edge weights. Let p = 〈v0 , v1 , . . . , vk 〉 be a shortest path from v0 to vk . Let i, j be indices with ­ ® 0 ≤ i ≤ j ≤ k, and let p i j be the subpath v i , v i+1 , . . . , v j . Then p i j is a shortest path from v i to v j . Proof. We are given that p : v0 p 0i / vi pi j / vj p jk / vk is a shortest path. If there were a shorter path from v i to v j , say v i obtain a shorter path from v0 to vk : 0 p : v0 p 0i / vi p0i j / vj p jk p0i j v j , then we could patch it in to / vk , which is impossible. 6.2.1 All Pairs In recitation, we saw Dijkstra’s algorithm (a greedy algorithm) for finding all shortest paths from a single source, for graphs with nonnegative edge weights. In this lecture, we will solve the problem of finding the shortest paths between all pairs of vertices. This information is useful in many contexts, such as routing tables for courier services, airlines, navigation software, Internet traffic, etc. The simplest way to solve the all-pairs shortest path problem is to run Dijkstra’s algorithm |V | times, once with each vertex as the source. This would take time |V |· T D IJKSTRA , which, depending on the implementation of the min-priority queue data structure, would be ¡ ¢ ¡ ¢ Linear array: O V3 +VE = O V3 ¡ ¢ Binary min-heap: O V E lg V ¡ ¢ Fibonacci heap: O V 2 lg V + V E . Lec 6 – pg. 2 of 5 However, Dijkstra’s algorithm only works if the edge weights are nonnegative. If we wanted to allow negative edge weights, we could instead use the slower Bellman–Ford algorithm once per vertex: ¡ ¢ dense graph ¡ ¢ O V 2 E −−−−−−−−→ O V 4 . We had better be able to beat this! 6.2.2 Formulating the All-Pairs Problem When we solved the single-source shortest paths problem, the shortest paths were represented on the actual graph, which was possible because subpaths of shortest paths are shortest paths and because there was only one source. This time, since we are considering all possible source vertices at once, it is difficult to conceive of a graphical representation. Instead, we will use matrices. Let n = |V |, and arbitrarily label the vertices 1, 2, . . . , n. We will format the output of the all-pairs ¡ ¢ algorithm as an n × n distance matrix D = d i j , where ( d i j = δ(i, j) = weight of a shortest path from i to j if a path exists ∞ if no path exists, ¡ ¢ together with an n × n predecessor matrix Π = π i j , where ( πi j = if i = j or there is no path from i to j N IL, predecessor of j in our shortest path i j otherwise. Thus, the ith row of Π represents the sinlge-source shortest paths starting from vertex i, and the ith ¡ ¢ row of D gives the weights of these paths. We can also define the weight matrix W = w i j , where ( wi j = weight of edge (i, j) if edge exists ∞ otherwise. Because G is a directed graph, the matrices D, Π and W are not necessarily symmetric. The existence of a path from i to j does not tell us anything about the existence of a path from j to i. 6.3 The Floyd–Warshall Algorithm The Floyd–Warshall algorithm solves the all-pairs shortest path problem in Θ(V 3 ) time. It allows negative edge weights, but assumes that there are no cycles with negative total weight.1 The Floyd– Warshall algorithm uses dynamic programming based on the following subproblem: What are D and Π if we require all paths to have intermediate vertices2 taken from the © ª © ª set 1, . . . , k , rather than the full 1, . . . , n ? 1 If there were a negative-weight cycle, then there would exist paths of arbitrarily low weight, so some vertex pairs would have a distance of −∞. In particular, the distance from a vertex to itself would not necessarily be zero, so our algorithms would need some emendations. 2 An intermediate vertex in a path p = 〈v , . . . , v 〉 is one of the vertices v , . . . , v n 0 1 n−1 , i.e., not an endpoint. Lec 6 – pg. 3 of 5 1 1 2 2 i i 3 3 j j k k Figure 6.2. Either p passes through vertex k, or p doesn’t pass through vertex k. As k increases from 0 to n, we build up the solution to the all-pairs shortest path problem. The base case k = 0 (where no intermediate vertices are allowed) is easy: ³ ´ ³ ´ (0) D (0) = d i j and Π(0) = π(0) ij , where ( d (0) ij = 0 if i = j wi j otherwise ) ( π(0) ij and = N IL if i = j or w i j = ∞ i if i 6= j and there is an edge from i to j ) . In general, let D (k) and Π(k) be the corresponding matrices when paths are required to have all their intermediate vertices taken from the set {1, . . . , k}. (Thus D = D (n) and Π = Π(n) .) What we need is a way of computing D (k) and Π(k) from D (k−1) and Π(k−1) . The following observation provides just that. Observation. Let p be a shortest path from i to j when we require all intermediate vertices to come from {1, . . . , k}. Then either k is an intermediate vertex or k is not an intermediate vertex. Thus, one of the following: 1. k is not an intermediate vertex of p. Thus, p is also a shortest path from i to j when we require all intermediate vertices to come from {1, . . . , k − 1}. 2. k is an intermediate vertex for p. Thus, we can decompose p as p: i p ik /k pk j / j, where p ik and p k j are shortest paths when we require all intermediate vertices to come from {1, . . . , k − 1}, by (a slightly modified version of) Proposition 6.1. This tells us precisely that ( d (k) ij = wi j min if k = 0 ³ −1) d (k , ij −1) −1) d (k + d (k ik kj Thus we may compute D by the following procedure: Lec 6 – pg. 4 of 5 ´ if k ≥ 1. Algorithm: F LOYD –WARSHALL(W) ¯ ¯ 1 n ← ¯W.rows¯ 2 3 4 5 6 7 8 9 D (0) ← W Π(0) ← an n × n matrix with all entries initially N IL for k ← 1 to n do ³ ´ ³ ´ (k) (k) Let D (k) = d i j and Π(k) = π i j be new n × n matrices for i ← 1 to n do for j ← 1 to n do ³ ´ (k−1) (k−1) (k−1) d (k) ← min d , d + d ij ij ik kj −1) −1) if d (k) ≤ d (k + d (k then ij ik kj −1) π(k) ← π(k ij ij 10 11 12 else −1) π(k) ← π(k ij kj ¡ ¢ Being in a triply nested loop, lines 8–12 are executed n3 times. Thus, the algorithm runs in Θ n3 = ¡ ¢ ¡ ¢ Θ V 3 time. Note that while the algorithm as written requires Θ V 3 space, reasonable implementa¡ 2¢ tions will only require Θ V space: we only ever need to memorize four matrices at a time, namely D (k) , D (k−1) , Π(k) , Π(k−1) . In fact, we can get away with memorizing even less—see Exercise 25.2-4 of CLRS. Lec 6 – pg. 5 of 5 Lecture 7 All-Pairs Shortest Paths II Supplemental reading in CLRS: Section 25.3 7.1 Johnson’s Algorithm ¡ ¢ The Floyd–Warshall algorithm runs in Θ V 3 time. Recall that, if all edge weights are nonnegative, then repeated application of Dijkstra’s algorithm using a Fibonacci heap gives the all-pairs shortest ¡ ¢ paths in Θ V 2 lg V + V E time. If G is not dense (i.e., if |E | is not on the order of |V |2 ), then Dijkstra’s algorithm asymptotically outperforms the Floyd–Warshall algorithm. So our goal for the first half of this lecture will be to reduce the problem of finding all-pairs shortest paths on an arbitrary graph to that of finding all-pairs shortest paths on a graph with nonnegative edge weights. 7.1.1 Reweighting How might we “eliminate” negative edge weights? The naïve way would be to add some constant N to all the edge weights, where N is chosen large enough to make all the edge weights nonnegative. Unfortunately, this transformation does not preserve shortest paths. Given a path p : P 〈 u = v0 , v1 , . . . , vk = v〉 from u to v with total weight w(p) = ki=1 w(v i−1 , v i ), the weight of p in the new graph would be w(p) + kN. Thus, a very light path in the original graph might not be so attractive in the new graph, if it contained lots of edges (see Figure 7.1). We will need to do something smarter; the change in weight must vary from edge to edge. Our goal is to assign to each edge (u, v) a new weight w(u, b v), such that the new weight function w b has the following two properties: 1. For each pair of vertices u, v ∈ V and each path p : u v, the old weight w(p) is minimal if and only if the new weight w(p) b is minimal (among paths from u to v). 2. For all edges (u, v), the new weight w(u, b v) is nonnegative. −4 −4 −5 1 1 0 Figure 7.1. Adding a constant to the weight of every edge does not preserve shortest paths. A clever, simple solution to this problem is as follows. Let h : V → R be a scalar-valued function on the vertices. Assign w(u, b v) = w(u, v) + h(v) − h(u). This weighting has the property that, for a path p : 〈 u = v0 , v1 , . . . , vk = v〉, the new weight is w(p) b = = k X i =1 k X w b (v i−1 , v i ) w (v i−1 , v i ) + h (v i ) − h (v i−1 ) i =1 = w(p) + k X h (v i ) − h (v i−1 ) i =1 = w(p) + h (vk ) − h (v0 ) = w(p) + h(v) − h(u), by telescoping. Thus, the amount by which p changes depends only on the endpoints u and v, not on the intermediate vertices. Consequently, a path p : u v is minimal with respect to w b if and only if it is minimal with respect to w. Also, if u = v, then w(p) b = w(p), so the weight of cycles is unchanged. Thus, negative-weight cycles are preserved, so the distance from u to v with respect to w b is −∞ if and only if the distance from u to v with respect to w is −∞. Our task now is to choose an appropriate function h so that the edge wieghts w b become nonnegative. There is a magical answer to this question: Let s be a new vertex, and construct an augmented graph G 0 = (V 0 , E 0 ) with V 0 = V ∪{ s}. Take E 0 to consist of all the edges in E, plus a directed edge from s to each vertex in V (see Figure 7.2). Then let h(v) = −δ(s, v), where δ(s, v) is the weight of the shortest path from s to v. Note that Proposition 6.1 implies the “triangle inequality” δ(s, v) ≤ δ(s, u) + δ(u, v) (with equality if one of the shortest paths from s to v passes through u). Since δ(u, v) ≤ w(u, v), we have δ(s, v) ≤ δ(s, u) + w(u, v). Hence w(u, b v) = w(u, v) + h(v) − h(u) = w(u, v) − δ(s, v) + δ(s, u) ¡ ¢ ≥ w(u, v) − δ(s, u) + w(u, v) + δ(s, u) = 0. Applying the weights w b to G, we obtain the following algorithm, due to Johnson: Lec 7 – pg. 2 of 7 0 0 3 4 8 0 7 1 s 2 −4 −5 0 6 0 Figure 7.2. The augmented graph G 0 in Johnson’s algorithm, consisting of G together with a special vertex s from which there emanates an edge of weight 0 to each other vertex. Algorithm: J OHNSON(G) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 Construct the augmented graph G 0 , where G 0 .V = G.V ∪ { s} and G 0 .E = G.E ∪ © ª (s, v) with G 0 .w(s, v) = 0, for each v ∈ G.V if B ELLMAN –F ORD(G 0 , s) = FALSE then B if B ELLMAN –F ORD complains about a negative-weight cycle error “The input graph contians a negative-weight cycle.” else for each vertex v ∈ G.V do Set h(v) ← −δ(s, v), where δ(s, v) is as computed by B ELLMAN –F ORD for each edge (u, v) ∈ G.E do w(u, b v) ← G.w(u, v) + h(v) − h(u) ¡ ¢ ¡ ¢ Let D = d uv and Π = πuv be new n × n matrices for each vertex u ∈ G.V do Run D IJKSTRA(G, w, b u) to compute δb(u, v) for all v ∈ G.V for each vertex v ∈ G.V do B correct the shortest path weights ¡ ¢ d uv ← δb(u, v) − h(v) − h(u) πuv ← the predecessor of v as computed by D IJKSTRA(G, w, b u) above return D The only steps in this algorithm that take more than O(E) time are the call to B ELLMAN –F ORD on line 2 (which takes Θ(V E) time) and the |V | calls to D IJKSTRA on lines 11–16 (each of which ¡ ¢ takes Θ V lg V + E time using a Fibonacci heap). Thus, the total running time of J OHNSON is ¡ 2 ¢ Θ V lg V + V E . Lec 7 – pg. 3 of 7 7.2 Linear Programming Linear programming is a very general class of problem that arises frequently in diverse fields. The goal is to optimize a linear objective function subject to linear constraints (equalities and inequalities). That is, we want to maximize or minimize the function n X f (x1 , x2 , . . . , xn ) = a 1 x1 + a 2 x2 + · · · + a n xn = a jxj j =1 subject to the constraints (i) (i) g(i) (x1 , . . . , xn ) = b(i) 1 x1 + b 2 x2 + · · · + b n x n = and (i) (i) h(i) (x1 , . . . , xn ) = c(i) 1 x1 + c 2 x2 + · · · + c n x n = n X j =1 n X j =1 b(i) x j = p(i) j c(i) x j ≤ q(i) . j In vector form, we want to maximize or minimize aT x subject to the constraints Bx = p and Cx ≤ q (meaning that the ith entry of Cx is less than or equal to the ith entry of q for every i). Here B and C are matrices with n columns (and any number of rows). 7.2.1 Application: Bacterial Growth One of the many applications of linear programming is the problem of maximizing bacterial growth subject to biochemical constraints. Consider the transport pathways of a bacterial cell, which we model as a graph embedded in three-dimensional space (see Figure 7.3). Let x i be the flux through edge i. • At steady-state, the net flux through each vertex of the graph is zero, as wastes are not building up along the cell pathways. This amounts to a linear constraint at each vertex, e.g., x1 − x2 − x3 = 0 in the picture. • Each edge has a minimum (namely zero) and a maximum rate at which molecules can pass (i) through. This amounts to a linear constraint 0 ≤ x i ≤ vmax at each edge. • The objective function is the sum of fluxes contributing to biomass production (e.g., synthesis of amino acids and proteins). The availability of efficient solutions to linear programming problems has been valuable in the engineering of new bacterial strains with improved capabilities, e.g., to produce biofuels. Lec 7 – pg. 4 of 7 Food source Zoom in on a vertex: x1 x2 x3 Flow constraint: x1 − x2 − x3 = 0 Figure 7.3. Bacterial cell, with a graph representing its transport pathways. At steady-state, the net flux through each vertex is zero. 7.2.2 Difference Constraints as Graphs It turns out that shortest path problems are related to some forms of linear programming problems. Consider a system of linear constraints, each having the form x j − xi ≤ b k , in other words, a system of the form Ax ≤ b, where A is an m × n matrix, x is an n-dimensional vector and b is an m-dimensional vector, and where each row of A contains one 1, one −1, and all other entries are zero. As a concrete example, let’s consider the system     1 −1 0 0   0 x1  1  0 0 −1   −1      x2    1 0 −1   ≤  1  . 0    x3    5 −1 0 1 0 x4 −3 0 0 −1 1 This is equalivalent to the inequalities   x1 − x2 ≤ 0             x − x ≤ − 1 1 4     x2 − x4 ≤ 1 .        x3 − x1 ≤ 5          x4 − x3 ≤ −3 Lec 7 – pg. 5 of 7 0 v1 v2 0 0 v0 −1 5 1 0 v4 −3 v3 0 Figure 7.4. Constraint graph for the system of difference constraints Ax ≤ b, where A and b are as in our example. Two plausible solutions to this system are  −5 −3   x =  , 0 −4    0 2   x =  . 5 1 In general, it’s not hard to see that adding any multiple of (1, 1, 1, 1)T to a solution produces another solution. Let’s now construct a graph G, called the constraint graph for this system (see Figure 7.4). The © ª vertex set of G is V = v0 , . . . , vn . Draw a directed edge with weight 0 from v0 to each other vertex. Then, for each difference constraint x j − x i ≤ b, draw an edge from v i to v j with weight b. Theorem 7.1 (Theorem 24.9 of CLRS). Given a system Ax ≤ b of difference constraints, let G = (V , E) be the corresponding constraint graph. If G contians no negative-weight cycles, then   δ (v0 , v1 )    δ (v0 , v2 )    x= ..  .   δ (v0 , vn ) is a feasible solution. If G does contain negative-weight cycles, then there is no feasible solution. Proof. Assume there are no negative-weight cycles. Since shortest paths satisfy the triangle inequality, we have ¡ ¢ ¡ ¢ δ v0 , v j ≤ δ (v0 , v i ) + w v i , v j . | {z } | {z } | {z } xj xi b Rearranging, we obtain x j − x i ≤ b. Next, consider the case of a negative-weight cycle. Without loss of generality, assume the cycle is ­ ® of the form c = v1 , v2 , . . . , vk , v1 . Then any solution x to the system of difference constraints must satisfy Lec 7 – pg. 6 of 7 x2 − x1 x3 − x2 + xk − xk−1 x1 − xk 0 ≤ ≤ .. . ≤ ≤ ≤ w (v1 , v2 ) w (v2 , v3 ) w (vk−1 , vk ) w (vk , v1 ) w(c) < 0, a contradiction. Therefore there cannot be any solutions to the system of difference constraints. Lec 7 – pg. 7 of 7 Lecture 8 Randomized Algorithms I Supplemental reading in CLRS: Chapter 5; Section 9.2 Should we be allowed to write an algorithm whose behavior depends on the outcome of a coin flip? It turns out that allowing random choices can yield a tremendous improvement in algorithm performance. For example, • One can find a minimum spanning tree of a graph G = (V , E, w) in linear time Θ(V + E) using a randomized algorithm. • Given two polynomials p, q of degree n − 1 and a polynomial r of degree n − 2, one can check whether r = pq in linear time Θ(n) using a randomized algorithm. No known deterministic algorithms can match these running times. Randomized algorithms are generally useful when there are many possible choices, “most” of which are good. Surprisingly, even when most choices are good, it is not necessarily easy to find a good choice deterministically. 8.1 Randomized Median Finding Let’s now see how randomization can improve our median-finding algorithm from Lecture 1. Recall that the main challenge in devising the deterministic median-finding algorithm was this: Problem 8.1. Given an unsorted array A = A[1, . . . , n] of n numbers, find an element whose rank is sufficiently close to n/2. While the deterministic solution of Blum, Floyd, Pratt, Rivest and Tarjan is ingenious, a simple randomized algorithm will typically outperform it in practice. Input: An array A = A[1, . . . , n] of n numbers Output: An element x ∈ A such that 1 9 10 n ≤ rank(x) ≤ 10 n. Algorithm: F IND -A PPROXIMATE -M IDDLE(A) Pick an element from A uniformly at random. Note that the above algorithm is not correct. (Why not?) However, it does return a correct answer 8 with probability 10 . good A: Figure 8.1. Wouldn’t it be nice if we had some way of choosing an element that probably lies in the middle? 8.1.1 Randomized Median Finding Let’s put the elegant and simple algorithm F IND -A PPROXIMATE -M IDDLE to use in a randomized version of the median-finding algorithm from Lecture 1: Algorithm: R ANDOMIZED -S ELECT(A, i) 1. Pick an element x ∈ A uniformly at random. 2. Partition around x. Let k = rank(x). 3. • If i = k, then return x. • If i < k, then recursively call R ANDOMIZED -S ELECT(A[1, . . . , k − 1], i). • If i > k, then recursively call R ANDOMIZED -S ELECT(A[k + 1, . . . , i], i − k). Note two things about this algorithm: • Unlike F IND -A PPROXIMATE -M IDDLE, R ANDOMIZED -S ELECT makes several random choices: one in each recursive call. • Unlike F IND -A PPROXIMATE -M IDDLE, R ANDOMIZED -S ELECT is correct. The effect of bad random choices is to prolong the running time, not to generate an incorrect answer. For example, if we wanted to find the middle element (i = n/2), and if our random element x happened to always be the smallest element of A, then R ANDOMIZED -S ELECT would take Θ(n) · Tpartition = Θ(n2 ) time. (To prove this rigorously, a more detailed calculation is needed.) The latter point brings to light the fact that, for randomized algorithms, the notion of “worst-case” running time is more subtle than it is for deterministic algorithms. We cannot expect randomized algorithms to work well for every possible input and every possible sequence of random choices (in that case, why even use randomization?). Instead, we do one of two things: • Construct algorithms that always run quickly, and return the correct answer with high probability. These are called Monte Carlo algorithms. With some small probability, they may return an incorrect answer or give up on computing the answer. • Construct algorithms that always return the correct answer, and have low expected running time. These are called Las Vegas algorithms. Note that the expected running time is an average over all possible sequences of random choices, but not over all possible inputs. An algorithm that runs in expected Θ(n2 ) time on difficult inputs and expected Θ(n) time on easy inputs has worst-case expected running time Θ(n2 ), even if most inputs are easy. F IND -A PPROXIMATE -M IDDLE is a Monte Carlo algorithm, and R ANDOMIZED -S ELECT is a Las Vegas algorithm. Lec 8 – pg. 2 of 5 8.1.2 Running Time To begin our analysis of R ANDOMIZED -S ELECT, let’s first suppose all the random choices happen to be good. For example, suppose that every recursive call of R ANDOMIZED -S ELECT(A, i) returns 1 9 | A | and 10 | A |. Then, the argument to each recursive call of an element whose rank is between 10 R ANDOMIZED -S ELECT is at most nine-tenths the size of the argument to the previous call. Hence, if there are a total of K recursive calls to R ANDOMIZED -S ELECT, and if n is the size of the argument to the original call, then ¡ 9 ¢K n ≥ 1 =⇒ K ≤ log10/9 n. 10 (Why?) Of course, in general we won’t always get this lucky; things might go a little worse. However, it is quite unlikely that things will go much worse. For example, there is a 96% chance that (at least) one of our first 15 choices will be bad—quite high. But there is only a 13% chance that three of our first 15 choices will be bad, and only a 0.06% chance that nine of our first 15 choices will be bad. Now let’s begin the comprehensive analysis. Suppose the size of the original input is n. Then, ¡ 9 ¢r let T r be the number of recursive calls that occur after the size of A has dropped below 10 n, but ¡ 9 ¢r+1 n. For example, if it takes five recursive calls to shrink A down before it has dropped below 10 ¡ 9 ¢2 9 to size 10 n and it takes eight recursive calls to shrink A down to size 10 n, then T0 = 5 and T1 = 3. Of course, each T i is a random variable, so we can’t say anything for sure about what values it will take. However, we do know that the total running time of R ANDOMIZED -S ELECT(A = A[1, . . . , n], i) is T= logX 10/9 n Θ ³¡ r =0 ¢ 9 r 10 n ´ · Tr . Again, T is a random variable, just as each T r is. Moreover, each T r has low expected value: Since each recursive call to R ANDOMIZED -S ELECT has a 4/5 chance of reducing the size of A by a factor of 9 10 , it follows that (for any r) £ ¤ ¡ ¢s Pr T r > s ≤ 51 . Thus, the expected value of T r satisfies E [T r ] ≤ ∞ ¡ ¢ X 5 s s 15 = . 16 s=0 (The actual value 5/16 is not important; just make sure you are able to see that the series converges quickly.) Thus, by the linearity of expectation, " # logX ³¡ ¢ ´ 10/9 n 9 r E [T] = E Θ 10 n · T r r =0 = logX 10/9 n Θ ³¡ r =0 ¢ 9 r 10 n ´ · E [T r ] ³¡ ¢ ´ 10/9 n 5 logX 9 r O 10 n 16 r=0 à ! logX 10/9 n ¡ ¢ 9 r =O n 10 ≤ r =0 Lec 8 – pg. 3 of 5 à ≤O n ∞ ¡ X r =0 ¢ 9 r 10 ! = O(n), since the geometric sum 8.2 ¢ 9 r r =0 10 P∞ ¡ converges (to 10). Another Example: Verifying Polynomial Multiplication Suppose we are given two polynomials p, q of degree n − 1 and a polynomial r of degree 2n − 2, and we wish to check whether r = pq. We could of course compute pq in Θ(n lg n) time using the fast Fourier transform, but for some purposes the following Monte Carlo algorithm would be more practical: Algorithm: V ERIFY-M ULTIPLICATION(p, q, r) 1. Choose x uniformly at random from {0, 1, . . . , 100n − 1}. 2. Return whether p(x) · q(x) = r(x). This algorithm is not correct. It will never report no if the answer is yes, but it might report yes when the answer is no. This would happen if x were a root of the polynomial pq − r, which has degree at most 2n − 2. Since pq − r can only have at most 2n − 2 roots, it follows that £ ¤ 2n − 2 Pr wrong answer ≤ < 0.02. 100n Often, the performance improvement offered by this randomized algorithm (which runs in Θ(n) time) is worth the small risk of (one-sided) error. 8.3 The Markov Bound The Markov bound states that it is unlikely for a nonnegative random variable X to exceed its expected value by very much. Theorem 8.2 (Markov bound). Let X be a nonnegative random variable with positive expected value.1 Then, for any constant c > 0, we have h i 1 Pr X ≥ c · E [X ] ≤ . c Proof. We will assume X is a continuous random variable with probability density function f X ; the discrete case is the same except that the integral is replaced by a sum. Since X is nonnegative, so is E [X ]. Thus we have Z ∞ E [X ] = x f X (x) dx 0 Z ≥ ∞ c E [X ] x f X (x) dx 1 i.e., Pr X > 0 > 0. £ ¤ Lec 8 – pg. 4 of 5 ≥ c E [X ] Z ∞ c E [X ] f X (x) dx h i = c E [X ] · Pr X ≥ c E [X ] . Dividing through by c E [X ], we obtain h i 1 ≥ Pr X ≥ c E [X ] . c Corollary 8.3. Given any constant c > 0 and any Las Vegas algorithm with expected running time T, we can create a Monte Carlo algorithm which always runs in time cT and has probability of error at most 1/c. Proof. Run the Las Vegas algorithm for at most cT steps. By the Markov bound, the probability of not reaching termination is at most 1/c. If termination is not reached, give up. Lec 8 – pg. 5 of 5 Lecture 9 Randomized Algorithms II Supplemental reading in CLRS: Appendix C; Section 7.3 After Lecture 8, several students asked whether it was fair to compare randomized algorithms to deterministic algorithms. Determistic algorithm: • always outputs the right answer • always runs efficiently Randomized algorithm: • may sometimes not output the right answer • may sometimes not run efficiently Is this “fair”? Of course not. We demand less from randomized algorithms. But in exchange, we expect randomized algorithms to do more—to be more efficient, or to be simple and elegant. Separately, we might ask whether it is useful to employ randomized algorithms. The answer here is that it depends: • on the situation • on the available alternatives • on the probability of error or inefficiency. Lots of real-life situations call for randomized algorithms. For example, Google Search and IBM’s Jeopardy!-playing computer Watson both employ randomized algorithms. 9.1 The Central Limit Theorem and the Chernoff Bound Ideally, a randomized algorithm will have the property that £ ¤ n→∞ Pr bad things happening −−−−→ 0. Before further analyzing randomized algorithms, it will be useful to establish two basic facts about probability theory: the central limit theorem and the Chernoff bound. The central limit theorem states that the mean of a large number of independent copies of a random variable is approximately normal, as long as the variable has finite variance. The normal distribution (or Gaussian distribution) of mean 0 and variance 1, denoted N(0, 1), is defined by the probability density function ¡ ± ¢ 1 f N(0,1) (x) = p exp − x2 2 . 2π Figure 9.1. A tail of the normal distribution. Taking an affine transform of this standard Gaussian, we obtain the normal distribution with arbi¡ ¢ trary mean µ and arbitrary variance σ2 , denoted N µ, σ2 . Its probability density function is µ ¶ (x − µ)2 1 exp − . f N (µ,σ2 ) (x) = p 2σ2 2πσ2 Due in part to the following theorem, the Gaussian distribution is extremely important in all of probability theory. Theorem 9.1 (Central Limit Theorem). Suppose X 1 , X 2 , . . . are i.i.d.1 random variables, each with ¢ p ¡ +X n finite mean µ and finite variance σ2 . Let S n = X 1 +... , and let Yn = n S n − µ . Then the variables n Y1 , Y2 , . . . converge to a normal distribution: ¡ ¢ d Yn − → N 0, σ2 . The precise meaning of this theorem is somewhat technical. Again, in effect it means that the average of a large number of independent copies of a random variable X is approximately normally distributed. Thus, the probability that this average exceeds its expected value by r is approximately the probability that a normal random variable with variance σ2 exceeds its expected value by r. While this heuristic reasoning does not rigorously prove any specific bound, there are a number of bounds that have been worked out for certain common distributions of X . One such bound is the Chernoff bound, one version of which is as follows: Theorem 9.2 (Chernoff bound). Let Y ∼ B(n, p) be a random variable representing the total number of heads in a series of n independent coin flips, where each flip has probability p of coming up heads. Then, for all r > 0, we have h i ¡ ± ¢ Pr Y ≥ E [Y ] + r ≤ exp −2r 2 n . The distribution in Theorem 9.2 is called the binomial distribution with parameters (n, p). ¡n¢ m The probability of m heads is m p (1 − p)n−m . There are several different versions of the Chernoff bound—some with better bounds for specific hypotheses, some more general. Section C.5 of CLRS gives a classic proof of the Chernoff bound, ¡ ¢ applying the Markov inequality to the random variable exp α(Y − E [Y ]) for a suitable constant α. We give a different proof here2 , due to Impagliazzo and Kabanets, 2010. 1 independent and identically distributed 2 Of a different theorem, actually. The version we prove here is (stronger than) Exercise C.5-6 of CLRS. Lec 9 – pg. 2 of 9 0 1 2 3 4 5 6 Figure 9.2. Galton’s board is a toy in which, at each peg, the ball has probability 12 of moving left and probability 21 of moving right. The number of the bin into which it falls is a binomial random variable with parameters (n, 12 ). If n is large and you drop a large number of balls into the board, the collection of balls in the bins underneath will start to look like a Gaussian probability density curve. Lec 9 – pg. 3 of 9 Proof of Theorem 9.2. Let Y = X1 + · · · + X n, where X 1 , . . . , X n are independent {0, 1}-valued random variables, each having probability p of being 1. Thus Y is a binomial random variable with parameters (n, p). Let S be a random subset of {1, . . . , n} such that each element of {1, . . . , n} independently has probability q of being included in S, where q is a parameter that will be determined later. The probability that X i = 1 for all i ∈ S is " # " # ¯ h i ^ ^ ¯ Pr X i = 1 = Pr X i = 1 ¯ Y ≥ E [Y ] + r Pr Y ≥ E [Y ] + r i ∈S i ∈S " # h i ^ ¯¯ + Pr ¯ Y < E [Y ] + r Pr Y < E [Y ] + r i ∈S " # ¯ h i ¯ X i = 1 ¯ Y ≥ E [Y ] + r Pr Y ≥ E [Y ] + r . ^ ≥ Pr i ∈S (9.1) Meanwhile, " Pr ^ i ∈S # ¯ ¯ X i = 1 ¯ Y ≥ E [Y ] + r ≥ (1 − q)n−(E [Y ]+r) , (9.2) V since the condition i∈S X i = 1 is equivalent to the condition that S doesn’t include the indices of any zeros, and we are conditioning on the hypothesis that there are at most n − (E [Y ] + r) zeros. (This is the “key idea” of the proof. Notice that we have switched from conditioning on S followed by the X i ’s to conditioning on the X i ’s followed by S.) Combining (9.1) and (9.2), we obtain " # h i ^ Pr X i = 1 ≥ (1 − q)n−(E [Y ]+r) Pr Y ≥ E [Y ] + r . (9.3) i ∈S £V ¤ We could also compute Pr i∈S X i = 1 by conditioning on S in a different way: " # ´ ´ ³ n ³¡ ¢ X ^ n k k n− k · p . Pr Xi = 1 = q (1 − q) k | {z } {z } i ∈S k=0 | ¯ i h h i Pr |S |= k Pr V i ∈S ¯ X i =1 ¯ |S |= k Recall the binomial formula, which states (a + b)n = n ¡ ¢ X n k=0 k a k b n− k . In the present case, we have " Pr # ¡ ¢n Xi = 1 = qp +1− q . ^ (9.4) i ∈S Combining (9.3) and (9.4), we obtain h i ¡ ¢n (1 − q)n−(E [Y ]+r) Pr Y ≥ E [Y ] + r ≤ q p + 1 − q . Rearranging and using the fact that E [Y ] = np, we have equivalently ¡ ¢n h i qp +1− q . Pr Y ≥ E [Y ] + r ≤ (1 − q)n−(np+r) Lec 9 – pg. 4 of 9 (9.5) Equation (9.5) holds for arbitrary q, so we might as well assign a value to q for which the right side is minimal. Let r/n ¢ q= ¡ ; r p + n (1 − p) after some expansion, the right side of (9.5) becomes à !1−( p+ r ) n ! p+ r à n n ¡ ¡ ¢¢ 1 − p p  = exp − n f r ,  ¡ ¢ r r n p+ n 1− p+ n where3 ´ ¡ ´ ¢ ³ ¡ ¢ ³ f (x) = p + x ln 1 + px + 1 − (p + x) ln 1 − 1−x p . Below we will show that f (x) ≥ 2x2 for all x ∈ R such that f (x) ∈ R, so that (9.5) becomes h i ¡ ¡ ¢¢ Pr Y ≥ E [Y ] + r ≤ exp − n f nr ³ ¡ ¢2 ´ ≤ exp −2n nr ³ ´ 2 = exp − 2rn . This will complete the proof. All that remains is to show that f (x) ≥ 2x2 for all x such that f (x) ∈ R—namely, for − p < x < 1 − p, though we will only need the range 0 ≤ x < 1 − p.4 This can be shown by calculus: the first and second derivatives of f are ³ ´ ³ ´ f 0 (x) = ln 1 + px − ln 1 − 1−x p , f 00 (x) = 1 . (1 − (p + x))(p + x) In particular, f 00 (x) is minimized when p + x = 12 , at which point f 00 (x) = 4. Thus, the fundamental theorem of calculus gives Z x Z x f 0 (x) = f 0 (0) + f 00 (t) dt ≥ 4 dt = 4x. | {z } t=0 t=0 0 Applying the fundamental theorem of calculus again, we have Z x Z x f (x) = f (0) + f 0 (t) dt ≥ 4t dt = 2x2 . |{z} t=0 t=0 0 3 The definition of f may seem unmotivated here, but in fact f (x) is properly interpreted as the relative entropy between two Bernoulli random variables with parameters p and p + x, respectively. Relative entropy is a measure of the “distance” between two probability distributions. 4 Note that x = r will always fall within this range, as we must have E [Y ] + r ≤ n. (Or at least, the theorem becomes n trivial when nr > 1 − p.) Lec 9 – pg. 5 of 9 ai pivot Figure 9.3. In a given iteration of Q UICKSORT, there will be several working subarrays, each in the process of being sorted by its own recursive call to Q UICKSORT. 9.2 Analysis of Q UICKSORT Recall the procedure Q UICKSORT, which sorts an array A = 〈a 1 , . . . , a n 〉 in place (see Figure 9.3): Algorithm: Q UICKSORT(A) 1. Choose an element x ∈ A uniformly at random. We will call x the “pivot.” 2. Partition A around x. 3. Sort the elements less than x recursively. 4. Sort the elements greater than x recursively. We saw in recitation that the expected running time of Q UICKSORT is Θ(n lg n). In this lecture we’ll use the Chernoff bound to show that, for some constant c ≥ 1, the probability that Q UICKSORT takes more than cn lg n time on any given input is at most n1 . In fact, what we’ll show is that after time cn lg n, the probability that the working subarray containing any given element a ∈ A consists of more than one element is at most n12 . Once this is established, the union bound shows that the probability that there exists any working subarray with more than one element is at most n1 . That is to say, let A = 〈a 1 , . . . , a n 〉 and let E i be the event that after time cn lg n, the working subarray containing a i has more than one element. Then the union bound states that " # · ¸ µ ¶ n n X [ £ ¤ 1 1 ∃ working subarray with Pr E i ≤ n 2 = . Pr = Pr Ei ≤ more than one element n n i =1 i =1 (In general, the union bound says that the probability of a finite or countably infinite union of events E i is at most the sum of their individual probabilities. The intuition is that equality is achieved when the E i ’s are pairwise disjoint, and any overlap only diminishes the size of the union.) All that remains is to find c such that, for any array A, £ ¤ 1 Pr E i ≤ 2 n for all i = 1, . . . , n. As in Lecture 8, we define a “good” pivot x in an array A to be a pivot such that 1 9 | | | | 10 A ≤ rank(x) ≤ 10 A . 9 as big. Let X i,k be the Thus, a good pivot splits the current array into two subarrays each at most 10 indicator random variable ( 0 if, in the kth iteration, a good pivot was chosen for a i ’s subarray X i,k = 1 otherwise. Lec 9 – pg. 6 of 9 (By “the kth iteration,” we mean all calls to Q UICKSORT at recursive depth k. Thus, the kth iteration will consist of 2k−1 recursive calls to Q UICKSORT, assuming none of the subarrays have shrunk to size 1 yet.) Note two things: £ ¤ • E X i,k = 0.2 for all i and all k. (To reduce the amount of bookkeeping we have to do, we will assume that pivots continue to be chosen even after an element’s subarray is reduced to size 1, and that the probability of choosing a good pivot continues to be 0.8.) • For a fixed k, the variables X 1,k , X 2,k , . . . , X n,k are dependent on each other. However, for a fixed i, the variables X i,1 , X i,2 , . . . are independent of each other. For a fixed i and a fixed K, we have " E K X # X i,k = 0.2K. k=1 Thus, by the Chernoff bound (with r = 0.1K), we have " # ¶ µ K X 2(0.1K)2 X i,k > 0.2K + 0.1K ≤ exp − Pr = exp (−0.02K) . K k=1 Now let K = 100 ln n. We obtain " Pr K X # X i,k > 0.3K ≤ k=1 1 . n2 PK On the other hand, if k=1 X i,k ≤ 0.3K, then at least 0.7K of the first K pivots were good, and the size of the working subarray containing a i after K steps is at most µ ¶0.7K µ ¶70 ln n 9 9 n· = n· = n70 ln(9/10)+1 ≈ n−6.3 < 1 10 10 (i.e., the size is at most 1; the reason for the fractional number is that we are ignoring issues of rounding). Thus, after K = 100 ln n iterations, the probability that the working subarray containing a i has more than one element is at most n12 . So by the union bound, the probability that Q UICKSORT requires more than K iterations is at most n1 . While the sizes of the inputs to each iteration and the total number of iterations required are random, the running time of the non-recursive part of Q UICKSORT is deterministic: It takes Θ(A) time to pick a random element and partition around that element. Let’s say it takes at most τ A time to do this, where τ is some appropriately chosen constant. Now, if the sizes of the inputs to the kth P iteration are m 1 , . . . , m ` , then the quantities m 1 , . . . , m ` are random, but we know `λ=1 m λ = n, so the P` total amount of time required by the kth iteration is at most λ=1 τ m λ = τ n. Thus, the probability that Q UICKSORT takes more than K τ n = (100τ)n ln n time is at most 1/n. This is what we set out to show in the beginning of this section, with c = 100τ. To recap: How can we argue about randomized algorithms? • Identify the random choices in the algorithm. • After fixing the random choices, we have a deterministic algorithm. In this example, the ¯random choices were the pivots. Fixing the sequence of choices of pivots, we ¯ ¯ ¯ were able to analyze A i,k , the size of the subarray containing a i after k iterations, for each i and k (see Figure 9.4). Lec 9 – pg. 7 of 9 Random choices at depth 1 Random choices at depth 2 .. . ¯ ¯ ¯ A 1,1 ¯ = n ¯ ¯ ¯ A 1,2 ¯ = 0.3n .. . ¯ ¯ ¯ A 2,1 ¯ = n ¯ ¯ ¯ A 2,2 ¯ = 0.8n .. . ··· ··· .. . Figure 9.4. The values 0.3n and 0.8n are just examples, and of course depend on the sequence of random choices. Once the sequence of random choices is fixed, we would like to know how many rows down we must go before all entries are 1. 9.3 Monte Carlo Sampling Suppose that out of a large population U, there occurs a certain phenomenon in a subset S ⊆ U. S| Given an element x, we are able to test whether x ∈ S. How can we efficiently estimate ||U | ? For example, how can we efficiently estimate what fraction of the world’s population is left-handed, or what percentage of voters vote Republican? There is an elegant simple solution to this problem: Choose a k-element¯ sample A = { x1 , . . . , xk } ⊆ © ª¯ S| ¯ b | A ∩S | 1 ¯ U uniformly at random. We then estimate ||U | by the quantity S = | A | = k i : x i ∈ S , figuring that the chance that an element a ∈ A belongs to S ought to be the same as the chance that a general element x ∈ U belongs to S. Thus, we are relying heavily on our ability to pick a truly uniform 5 b sample. ³ ´ Assuming this reliance is safe, k S = | A ∩ S | is a binomial random variable with parameters |S | S| k, |U | . In particular, the expected value of Sb is the exact correct answer ||U | . Thus, the Chernoff bound states that for any r > 0, · ¸ h i £ ¤ ¡ ¢ |S | Pr Sb ≥ + r = Pr k Sb ≥ E k Sb + kr ≤ exp −2kr 2 . |U | i h S| Exercise 9.1. How can we use the Chernoff bound to give an upper bound on Pr Sb ≤ ||U − r ? (Hint: | replace S by its complement U \ S.) 9.4 Amplification Suppose we are given a Monte Carlo algorithm which always runs in time T and returns the correct answer with probability 2/3. (Assume there is a unique correct answer.) Then, for any ² > 0, we can ¡ ¢ create a new randomized algorithm which always runs in time O T lg 1² and outputs an incorrect answer with probability at most ². How? ¡ ¢ The plan is to run the original algorithm k times, where k = O lg 1² is a quantity which we will determine later. The algorithm should then output whichever answer occurred most frequently out of the k computed answers. To bound the probability of error, let I j be the indicator random variable P which equals 1 if the jth run returns an incorrect answer, and let I = kj=1 I j . Then I is a binomial ¡ 1¢ random variable with parameters k, 3 . Thus, by the Chernoff bound, Pr · ¸ h i £ ¤ ¡ 1 ¢ we ultimately return ≤ Pr I ≥ 12 k = Pr I ≥ E [I] + 61 k ≤ exp − 18 k . an incorrect answer 5 In practice, it could be quite hard to pick a uniform sample. For example, what if you wanted figure out what proportion of New Yorkers speak Chinese? How much of your random sampling should be done in Chinatown? Without the help of extensive census data, it can be hard to make unbiased choices. Lec 9 – pg. 8 of 9 Thus, the probability of error will be at most ² if we let k = 18 ln 1² . Exercise 9.2. Suppose that instead of returning the correct answer with probability 2/3, our Monte Carlo algorithm returned the correct answer with probability p. What conditions on p allow the above strategy to work? In terms of p and ², how many times must we run the Monte Carlo algorithm? Lec 9 – pg. 9 of 9 Lecture 10 Hashing and Amortization Supplemental reading in CLRS: Chapter 11; Chapter 17 intro; Section 17.1 10.1 Arrays and Hashing Arrays are very useful. The items in an array are statically addressed, so that inserting, deleting, and looking up an element each take O(1) time. Thus, arrays are a terrific way to encode functions © ª 1, . . . , n → T, where T is some range of values and n is known ahead of time. For example, taking T = {0, 1}, we find that an array A of n bits is a great way to store a subset of {1, . . . , n}: we set A[i] = 1 if and only if i is in the set (see Figure 10.1). Or, interpreting the bits as binary digits, we can use an n-bit array to store an integer between 0 and 2n − 1. In this way, we will often identify the set {0, 1}n with the set {0, . . . , 2n − 1}. What if we wanted to encode subsets of an arbitrary domain U, rather than just {1, . . . , n}? Or to put things differently, what if we wanted a keyed (or associative) array, where the keys could be arbitrary strings? While the workings of such data structures (such as dictionaries in Python) are abstracted away in many programming languages, there is usually an array-based solution working behind the scenes. Implementing associative arrays amounts to finding a way to turn a key into an array index. Thus, we are looking for a suitable function U → {1, . . . , n}, called a hash function. Equipped with this function, we can perform key lookup: ª array lookup hash function © U −−−−−−−−−→ 1, . . . , n −−−−−−−−−→ T (see Figure 10.2). This particular implementation of associative arrays is called a hash table. There is a problem, however. Typically, the domain U is much larger than {1, . . . , n}. For any hash function h : U → {1, . . . , n}, there is some i such that at least |Un | elements are mapped to i. The set A: 0 1 0 1 1 0 0 1 0 0 1 0 1 2 3 4 5 6 7 8 9 10 11 12 © ª © ª Figure 10.1. This 12-bit array encodes the set 2, 4, 5, 8, 11 ⊆ 1, . . . 12 . 1 2 ­ key3 , val3 ® h key1 = 5 ­ key1 , val1 ® ¡ ¢ h key2 = 6 ­ key2 , val2 ® ¡ ¢ h key3 = 3 4 ¡ ¢ 7 Figure 10.2. An associative array with keys in U and values in T can be implemented as a (U × T)-valued array equipped with a hash function h : U → {1, . . . , n}. h−1 (i) of all elements mapped to i is called the load on i, and when this load contains more than one of the keys we are trying to store in our hash table we say there is a collision at i. Collisions are problem for us—if two keys map to the same index, then what should we store at that index? We have to store both values somehow. For now let’s say we do this in the simplest way possible: storing at each index i in the array a linked list (or more abstractly, some sort of bucket-like object) ¡ ¢ consisting of all values whose keys are mapped to i. Thus, lookup takes O h−1 (i) time, which may be poor if there are collisions at i. Rather than thinking about efficient ways to handle collisions,1 let’s try to reason about the probability of having collisions if we choose our hash functions well. 10.2 Hash Families Without any prior information about which elements of U will occur as keys, the best we can do is to choose our hash function h at random from a suitable hash family. A hash family on U is a set H of functions U → {1, . . . , n}. Technically speaking, H should come equipped with a probability distribution, but usually we just take the uniform distribution on H , so that each hash function is equally likely to be chosen. If we want to avoid collisions, it is reasonable to hope that, for any fixed x1 , x2 ∈ U (x1 6= x2 ), the values h(x1 ) and h(x2 ) are completely uncorrelated as h ranges through the sample space H . This leads to the following definition: Definition. A hash family H on U is said to be universal if, for any x1 , x2 ∈ U (x1 6= x2 ), we have £ ¤ Pr h(x1 ) = h(x2 ) ≤ n1 . h∈H 1 If you are expecting lots of collisions, a more efficient way to handle things is to create a two-layered hash table, where each element of A is itself a hash table with its own, different hash function. In order to have collisions in a two-layer hash table, the same pair of keys must collide under two different hash functions. If the hash functions are chosen well (e.g., if the hash functions are chosen randomly), then this is extremely unlikely. Of course, if you want to be even more sure that collisions won’t occur, you can make a three-layer hash table, and so on. There is a trade-off, though: introducing unnecessary layers of hashing comes with a time and space overhead which, while it may not show up in the big-O analysis, makes a difference in practical applications. Lec 10 – pg. 2 of 7 Similarly, H is said to be ²-universal if for any x1 6= x2 we have £ ¤ Pr h(x1 ) = h(x2 ) ≤ ². h∈H The consequences of the above hypotheses with regard to collisions are as follows: Proposition 10.1. Let H be a universal hash family on U. Fix some subset S ⊆ U and some element x ∈ U. Pick h ∈ H at random. The expected number of elements of S that map to h(x) is at most 1 + |Sn | . In symbols, h¯ ¡ ¢¯i |S | E ¯ h−1 h(x) ¯ ≤ 1 + . h∈H n If H is ²-universal rather than universal, then the same holds when 1 + |Sn | is replaced by 1 + ² |S |. Proof. For a proposition ϕ with random parameters, let I ϕ be the indicator random variable which equals 1 if ϕ is true and equals 0 otherwise. The fact that H is universal means that for each x0 ∈ U \ { x} we have £ ¤ 1 E I h(x)=h(x0 ) ≤ . h∈H n Thus by the linearity of expectation, we have   h¯ i X ¯ ¡ ¢   E ¯ h−1 h(x) ∩ S ¯ = I x∈S + E  I h(x)=h(x0 )  h∈H = I x∈ S + X x0 ∈ S x0 6= x x0 ∈ S x0 6= x E h∈H £ I h(x)=h(x0 ) ¤ ≤ 1 + |S | · n1 . The reasoning is almost identical when H is ²-universal rather than universal. Corollary 10.2. For a hash table in which the hash function from a universal family, ³ is chosen ´ |S | insertion, deletion, and lookup have expected running time O 1 + n , where S ⊆ U is the set of keys which actually occur. If instead the hash family is ²-universal, then the operations have expected ¡ ¢ running time O 1 + ² |S | . Corollary 10.3. Consider a hash table of size n with keys in U, whose hash function is chosen from a universal hash family. Let S ⊆ U be the set of keys which actually occur. If |S | = O(n), then insertion, deletion, and lookup have expected running time O(1). Let H be a universal hash family on U. If |S | = O(n), then the expected load on each index is O(1). Does this mean that a typical hash table has O(1) load at each index? Surprisingly, the answer is no, even when the hash function is chosen well. We’ll see this below when we look at examples of universal hash families. Examples 10.4. © ª 1. The set of all functions h : U → 1, . . . , n is certainly universal. In fact, we could not hope to get any more balanced than this: Lec 10 – pg. 3 of 7 • For any x ∈ U, the random variable h(x) (where h is chosen at random) is uniformly dis© ª tributed on the set 1, . . . , n . • For any pair x1 6= x2 , the random variables h(x1 ), h(x2 ) are independent. In fact, for ¡ ¢ any finite subset { x1 , . . . , xk } ⊆ U, the tuple h(x1 ), . . . , h(xk ) is uniformly distributed on © ªk 1, . . . , n . ¢ ¡ The load on each index i is a binomial random variable with parameters |S | , n1 . Fact. When p is small and N is large enough that N p is moderately sized, the binomial distribution with parameters (N, p) is approximated by the Poisson distribution with parameter N p. That is, if X is a binomial random variable with parameters (N, p), then £ ¤ (N p)k − N p Pr X = k ≈ e k! (k ≥ 0). In our case, N = |S | and p = n1 . Thus, if L i is the load on index i, then ³ £ ¤ Pr L i = k ≈ ´ |S | k n −|S |/n k! e . For example, if |S | = n, then £ ¤ Pr L i = 0 ≈ e−1 ≈ 0.3679, £ ¤ Pr L i = 1 ≈ e−1 ≈ 0.3679, £ ¤ Pr L i = 2 ≈ 21 e−1 ≈ 0.1839, .. . Further calculation shows that, when |S | = n, we have · ¸ µ ¶ lg n E max L i = Θ . 1≤ i ≤ n lg lg n ³ ´ lg n Moreover, with high probability, max L i does not exceed O lg lg n . Thus, a typical hash table with |S | = n and h chosen uniformly from the set of all functions looks like Figure 10.3: about 37% of the buckets empty, about 37% of the buckets having one element, ³ and´ about 26% of the lg n buckets having more than one element, incuding some buckets with Θ lg lg n elements. 2. In Problem Set 4 we considered the hash family © ª H = h p : p ≤ k and p is prime , © ª © ª where h p : 0, . . . , 2m − 1 → 0, . . . , k − 1 is the function h p (x) = x mod p. In Problem 4(a) you proved that, for each x 6= y, we have £ ¤ m ln k Pr h p (x) = h p (y) ≤ . p k Lec 10 – pg. 4 of 7 1 2 3 4 . . . Maximum load = Θ ³ lg n lg lg n ´ . . . n Figure 10.3. A typical hash table with |S | = n and h chosen uniformly from the family of all functions U → {1, . . . , n}. 3. In Problem Set 5, we fixed a prime p and considered the hash family o n a ∈ Zm H = h~a : ~ p , where h~a : Zm p → Z p is the dot product h~a (~x) = ~x · ~ a= X xi a i (mod p). 4. In Problem Set 6, we fixed a prime p and positive integers m and k and considered the hash family n o H = h A : A ∈ Zkp×m , k where h A : Zm p → Z p is the function h A (~x) = A~x. 5. If H 1 is an ²1 -universal hash family of functions {0, 1}m → {0, 1}k and H 2 is an ²2 -universal hash family of functions {0, 1}k → {0, 1}` , then2 © ª H = H 2 ◦ H 1 = h2 ◦ h1 : h1 ∈ H 1 , h2 ∈ H 2 is an (²1 + ²2 )-universal hash family of functions {0, 1}m → {0, 1}` . To see this, note that for any x 6= x0 , the union bound gives £ ¤ Pr h 2 ◦ h 1 (x) = h 2 ◦ h 1 (x0 ) h 1 ∈H 1 h 2 ∈H 2 2 To fully specify H , we have to give not just a set but also a probability distribution. The hash families H and H 1 2 come with probability distributions, so there is an induced distribution on H 1 ×H 2 . We then equip H with the distribution induced by the map H 1 × H 2 → H , (h 1 , h 2 ) 7→ h 2 ◦ h 1 . You could consider this a mathematical technicality if you wish: if H 1 and H 2 are given uniform distributions (as they typically are), then the distribution on H 1 × H 2 is also uniform. The distribution on H need not be uniform, however: an element of H is more likely to be chosen if it can be expressed in multiple ways as the composition of an element of H 2 with an element of H 1 . Lec 10 – pg. 5 of 7 h ³ ´i = Pr h 1 (x) = h 1 (x0 ) or h 1 (x) 6= h 1 (x0 ) and h 2 ◦ h 1 (x) = h 2 ◦ h 1 (x0 ) h i h i ≤ Pr h 1 (x) = h 1 (x0 ) + Pr h 1 (x) 6= h 1 (x0 ) and h 2 ◦ h 1 (x) = h 2 ◦ h 1 (x0 ) ≤ ²1 + ²2 . In choosing the parameters to build a hash table, there is a tradeoff. Making n larger decreases the likelihood of collisions, and thus decreases the expected running time of operations on the table, but also requires the allocation of more memory, much of which is not even used to store data. In situations where avoiding collisions is worth the memory cost (or in applications other than hash tables, when the corresponding tradeoff is worth it), we can make n much larger than S. Proposition 10.5. Let H be a universal hash family U → {1, . . . , n}. Let S ⊆ U be the the set of keys ¡ ¢ that occur. Then the expected number of collisions is at most |S2 | · n1 . In symbols, " E h∈H # X x6= x0 ∈U I h(x)=h(x0 ) à ! |S | 1 ≤ · . 2 n ¡ ¢ Proof. There are |S2 | pairs of distinct elements in S, and each pair has probability at most causing a collision. The result follows from linearity of expectation. 1 n of Corollary 10.6. If n ≥ 100 |S |2 , then the expected number of collisions is less than 1/200, and the probability that a collision exists is less than 1/200. Proof. Apply the Markov bound. Thus, if n is sufficiently large compared to S, a typical hash table consists mostly of empty buckets, and with high probability, there is at most one element in each bucket. As we mentioned above, choosing a large n for a hash table is expensive in terms of space. While the competing goals of fast table operations and low storage cost are a fact of life if nothing is known about S in advance, we will see in recitation that, if S is known in advance, it is feasible to construct a perfect hash table, i.e., a hash table in which there are no collisions. Of course, the smallest value of n for which this is possible is n = |S |. As we will see in recitation, there are reasonably efficient algorithms to construct a perfect hash table with n = O (|S |). 10.3 Amortization What if the size of S is not known in advance? In order to allocate the array for a hash table, we must choose the size n at creation time, and may not change it later. If |S | turns out to be significantly greater than n, then there will always be lots of collisions, no matter which hash function we choose. Luckily, there is a simple and elegant solution to this problem: table doubling. The idea is to start with some particular table size n = O(1). If the table gets filled, simply create a new table of size 2n and migrate all the old elements to it. While this migration operation is costly, it happens infrequently enough that, on the whole, the strategy of table doubling is efficient. Let’s take a closer look. To simplify matters, let’s assume that only insertions and lookups occur, with no deletions. What is the worst-case cost of a single operation on the hash table? Lec 10 – pg. 6 of 7 • Lookup: O(1), as usual. • Insertion: O(n), if we have to double the table. Thus, the worst-case total running time of k operations (k = |S |) on the hash table is ¡ ¢ O (1 + · · · + k) = O k2 . The crucial observation is that this bound is not tight. Table doubling only happens after the second, fourth, eighth, etc., insertions. Thus, the total cost of k insertions is à ! lg Xk j k · O(1) + O 2 = O (k) + O (2k) = O (k) . j =0 Thus, in any sequence of insertion and lookup operations on a dynamically doubled hash table, the average, or amortized, cost per operation is O(1). This sort of analysis, in which we consider the total cost of a sequence of operations rather than the cost of a single step, is called amortized analysis. In the next lecture we will introduce methods of analyzing amortized running time. Lec 10 – pg. 7 of 7 Lecture 11 Amortized Analysis Supplemental reading in CLRS: Chapter 17 Data structures typically support several different types of operations, each with its own cost (e.g., time cost or space cost). The idea behind amortized analysis is that, even when expensive operations must be performed, it is often possible to get away with performing them rarely, so that the average cost per operation is not so high. It is important to realize that these “average costs” are not expected values—there needn’t be any random events.1 Instead, we are considering the worst-case average cost per operation in a sequence of many operations. In other words, we are interested in the asymptotic behavior of the function C(n) = ¢ 1 ¡ · worst-case total cost of a sequence of n operations n (possibly with some condition on how many times each type of operation may occur). “Worst-case” means that no adversary could choose a sequence of n operations that gives a worse running time. In this lecture we discuss three methods of amortized analysis: aggregate analysis, the accounting method, and the potential method. 11.1 Aggregate Analysis In aggregate analysis, one assumes that there is no need to distinguish between the different operations on the data structure. One simply asks, what is the cost of performing a sequence of n operations, of any (possibly mixed) types? Example. Imagine a stack2 S with three operations: • P USH(S, x) – Θ(1) – pushes object x onto the stack • P OP(S) – Θ(1) – pops and returns the top object of S 1 There could be random events, though. It makes sense to talk about the worst-case amortized expected running time of a randomized algorithm, for example: one considers the average expected cost per operation in the string of operations which has the longest expected running time. 2 A stack is a data structure S with two primitive operations: P USH(S, x), which stores x, and P OP(S), which removes and returns the most recently pushed object. Thus, a stack is “last in, first out,” whereas a queue is “first in, first out.” • M ULTIPOP(S, k) – Θ (min {|S | , k}) – pops the top k items from the stack (or until empty) and returns the last item: while S is not empty and k > 0 do x ← P OP(S) k ← k−1 return x Suppose we wish to analyze the the running time for a sequence of n P USH, P OP, and M ULTIPOP operations, starting with an empty stack. Considering individual operations without amortization, we would say that a M ULTIPOP operation could take Θ (|S |) time, and |S | could be as large as n − 1. So in the hypothetical worst case, a single operation could take Θ(n) time, and n such operations ¡ ¢ strung together would take Θ n2 time. However, a little more thought reveals that such a string of operations is not possible. While a single P OP could take Θ(n) time, it would have to be preceded by Θ(n) P USH operations, which are cheap. Taken together, the Θ(n) P USH operations and the one M ULTIPOP operation take Θ(n) · Θ(1) + 1 · Θ(n) = Θ(n) time; thus, each operation in the sequence takes Θ(1) time on average. In general, if there occur r M ULTIPOP operations with arguments k 1 , . . . , k r , then there must also occur at least k 1 +· · ·+ k r P USH operations, so that there are enough items in the stack to pop. (To simplify notation, we assume that k is never chosen larger than |S |.) Thus, in a string of n operations, the total cost of all non-O(1) operations is bounded above by O(n), so the total cost of all operations is O(n) · O(1) + O(n) = O(n). Thus, the amortized running time per operation is O(1). 11.2 Accounting Method Unlike aggregated analysis, the accounting method assigns a different cost to each type of operation. The accounting method is much like managing your personal finances: you can estimate the costs of your operations however you like, as long as, at the end of the day, the amount of money you have set aside is enough to pay the bills. The estimated cost of an operation may be greater or less than its actual cost; correspondingly, the surplus of one operation can be used to pay the debt of other operations. In symbols, given an operation whose actual cost is c, we assign an amortized (estimated) cost cb. The amortized costs must satisfy the condition that, for any sequence of n operations with actual costs c 1 , . . . , c n and amortized costs cb1 , . . . , cbn , we have n X i =1 cbi ≥ n X ci. i =1 As long as this condition is met, we know that the amortized cost provides an upper bound on the actual cost of any sequence of operations. The difference between the above sums is the total surplus or “credit” stored in the data structure, and must at all times be nonnegative. In this way, the accounting model is like a debit account. Example. Perhaps you have bought pre-stamped envelopes at the post office before. In doing so, you pay up-front for both the envelopes and the postage. Then, when it comes time to send a letter, no additional charge is required. This accounting can be seen as an amortization of the cost of sending a letter: Lec 11 – pg. 2 of 8 Operation Buy an envelope Mail a letter Actual cost c i 5¢ 44¢ Amortized cost cbi 49¢ 0¢ Obviously, for any valid sequence of operations, the amortized cost is at least as high as the actual cost. However, the amortized cost is easier to keep track of—it’s one fewer item on your balance sheet. Example. For the stack in Section 11.1, we could assign amortized costs as follows: Operation P USH P OP M ULTIPOP Actual cost c i 1 1 min {|S | , k} Amortized cost cbi 2 0 0 When an object is pushed to the stack, it comes endowed with enough credit to pay not only for the operation of pushing it onto the stack, but also for whatever operation will eventually remove it from the stack, be that a P OP, a M ULTIPOP, or no operation at all. 11.3 Potential Method The potential method is similar in spirit to the accounting method. Rather than assigning a credit to each element of the data structure, the potential method assigns a credit to the entire data structure as a whole. This makes sense when the total credit stored in the data structure is a function only of its state and not of the sequence of operations used to arrive at that state. We think of the credit as a “potential energy” (or just “potential”) for the data structure. The strategy is to define a potential function Φ which maps a state D to a scalar-valued potential Φ(D). Given a sequence of n operations with actual costs c 1 , . . . , c n , which transform the data structure from its initial state D 0 through states D 1 , . . . , D n , we define heuristic costs cbi = c i + Φ (D i ) − Φ (D i−1 ) . This rule can be seen as the conservation of energy: it says that the surplus (or deficit) cost cbi − c i caused by a given operation must be equal to the change in potential of the data structure caused by that operation. An operation which increases the potential of the system is assigned a positive heuristic cost, whereas an operation which decreases the potential of the system is assigned a negative heuristic cost. Summing the heuristic costs of all n operations, we find n X i =1 cbi = n X à = c i + Φ (D i ) − Φ (D i−1 ) (a telescoping sum) i =1 n X ! c i + Φ (D n ) − Φ (D 0 ) . i =1 Thus, the total credit stored in the data structure is Φ (D n ) − Φ (D 0 ). This quantity must remain nonnegative at all times in order to ensure that the amortized cost provides an upper bound on the actual cost of any sequence of operations. (Any function Φ which satisfies this property can be used as the potential function.) One often chooses the potential function Φ so that Φ (D 0 ) = 0; then one must check that Φ remains nonnegative at all times. Lec 11 – pg. 3 of 8 Example. Continuing the stack example from Sections 11.1 and 11.2, we define the potential of a stack S to be Φ(S) = |S |, the number of elements in the stack. An empty stack has zero potential, and clearly Φ is always nonnegative, so Φ is an admissible potential function. Then the heuristic costs of the stack operations are as follows: • A P USH operation increases the size of S by 1, and has actual cost c P USH = 1. Thus, the amortized cost of a P USH operation is cbP USH = c P USH + Φ (D new ) − Φ (D old ) = 1 + (|S old | + 1) − |S old | = 2. • A P OP operation decreases the size of S by 1, and has actual cost c P OP = 1. Thus, the amortized cost of a P OP operation is cbP OP = c P OP + Φ (D new ) − Φ (D old ) = 1 + (|S old | − 1) − |S old | = 0. • The operation M ULTIPOP(S, k) decreases the size of S by min {|S | , k}, and has actual cost c M ULTIPOP = min {|S | , k}. Thus, the amortized cost of a M ULTIPOP operation is cbM ULTIPOP = c M ULTIPOP + Φ (D new ) − Φ (D old ) = min {|S | , k} + (|S new | − |S old |) = 0. Thus, the amortized costs for this application of the potential method are the same as those we came up with using the accounting method in Section 11.2: Operation P USH P OP M ULTIPOP 11.4 Actual cost c i 1 1 min {|S | , k} Amortized cost cbi 2 0 0 Example: A Dynamically Resized Table We now extend the example we started in Lecture 10: a dynamically doubled table T. This time we will forget about keys, inserting each new value into the next available block in the array T.arr. If T.arr is full, then we create a new array of size 2 |T.arr|, copy over the old elements of T, perform the insertion, and reassign T.arr to point to that array. Likewise, in §11.4.5 we ask you to devise an extension of our data structure which supports deletions (specifically, pops) as well as insertions, and halves the array T.arr when appropriate in addition to doubling it when appropriate. In this way, T never takes up much more space than it has to. For now, assume that T only supports insertions and lookups. We will store in T two things: the underlying array T.arr, and an integer T.num which keeps track of the number of entries in T.arr which have yet been populated. Thus, the array is full and needs to be doubled upon the next insertion if and only if |T.arr| = T.num. So the table lookup will be a simple array lookup operation, which takes O(1) time, and the insertion operation will be defined as follows: Lec 11 – pg. 4 of 8 Algorithm: T ABLE -I NSERT(T, x) 1 2 3 4 5 6 7 8 9 10 11.4.1 if |T.arr| = 0 then T.arr ← new array of size 1 B If table is full, double it if T.num = |T.arr| then arr0 ← new array of size 2 |T.arr| Migrate all items from T.arr to arr0 T.arr ← arr0 B Insert x into the array T.arr[T.num] ← x T.num ← T.num + 1 Actual cost per operation A lookup operation on T is simply an array lookup, which takes O(1) time. The running time of T ABLE -I NSERT, on the other hand, depends on whether the table needs to be doubled, which depends on the number of insertions which have happened prior to the current one. Specifically, starting from an initial table of size zero, the cost c i of the ith insertion is   1 if the table has space        i − 1 to allocate the new table      i − 1 to migrate the old entries ci =    Θ (i) if the table has to be expanded:      1 to make the insertion       1 to reassign T.arr The “Chicken Little” analysis would say that the worst possible string of n operations would consist entirely of operations which have the worst-case single-operation running time Θ(i), so that the sequence has running time ¡ ¢ ¡ ¢ Θ 1 + 2 + · · · + n = Θ n2 . This bound is not tight, however, because it is not possible for every operation to be worst-case. 11.4.2 Aggregate analysis Aggregate analysis of our table’s operations amounts to keeping track of precisely which operations will be worst-case. The table has to be expanded if and only if the table is full, so ( ci = Θ(i) if i − 1 is zero or an exact power of 2 1 otherwise. Thus, a string of n operations, m of which are insertion operations, would take time ¡ ¢ ¡ ¢ T n − m lookups + T m insertions  ³¡ ´ X ¢  = n − m · Θ(1) +  Θ(1)  + k − 1 is not an exact power of 2 Lec 11 – pg. 5 of 8 X k − 1 is an exact power of 2  Θ(k) = ³¡ à ´ ¡ ¡ ¥ ¦¢¢ n − m · Θ(1) + m − 2 + lg(m − 1) · Θ(1) ¢ + blg(m X−1)c ³ ´ Θ 2j ! j =0 ¢ = Θ(n − m) + Θ(m) − Θ(lg m) + Θ(2m) ¡ = Θ(n) + Θ(m) = Θ(n) (since m ≤ n). Thus, a string of n operations takes Θ(n) time, regardless of how many of the operations are insertions and how many are lookups. Thus, aggregate analysis gives an amortized running time of O(1) per operation. 11.4.3 Accounting analysis We can set up a balance sheet to pay for the operations on our table as follows: • Lookups cost $1 • Insertions cost $3 – $1 to insert our element into the array itself – $2 to save up for the next time the array has to be expanded * One dollar will pay for migrating the element itself * One dollar will go to charity. Because the current incarnation of T.arr started with 12 |T.arr| elements pre-loaded, it follows that, when T.arr is full, only 21 |T.arr| elements will have enough money in the bank to pay for their migration. However, there will also be $ 21 |T.arr| of charity money saved up from the 21 |T.arr| most recent insertions. This charity is exactly sufficient to pay for migrating the elements that don’t have their own money. Thus, our balance sheet checks out; our assigned costs of $1 and $3 work. In particular, each operation’s cost is O(1). 11.4.4 Potential analysis In our potential analysis, the potential energy will play the role of money in the bank: when the array needs to be doubled, the potential energy will exactly pay for migrating the old elements into the new array. Thus, we want to find a function which equals zero immediately after doubling the table and grows to size |T.arr| as it is filled up. (To simplify notation, assume that copying an array of size k takes exactly k steps, rather than always using the more precise notation Θ(k).) The simplest function which satisfies these properties, and the one that will be most convenient to work with, is Φ(T) = 2 · T.num − |T.arr| . Right after expansion, we have T.num = 21 |T.arr|, so Φ(T) = 0; right before expasnion, we have T.num = |T.arr|, so Φ(T) = 2 · |T.arr| − |T.arr| = |T.arr|. The actual cost of a lookup is, as we said, c lookup = 1. (Again, to simplify notation, let’s just write 1 instead of Θ(1).) Using the potential Φ, the amortized cost of a lookup is cblookup = c lookup + Φ (Tnew ) − Φ (Told ) = c lookup = 1, Lec 11 – pg. 6 of 8 since Tnew = Told . To find the amortized cost of an insertion, we must split into cases. In the first case, suppose the table does not need to be doubled as a result of the insertion. Then the amortized cost is cb = c + Φ (Tnew ) − Φ (Told ) ¡ ¢ ¡ ¢ = 1 + 2 · Tnew .num − |Tnew .arr| − 2 · Told .num − |Told .arr| ¡ ¢ ¡ ¢ = 1 + 2 Tnew .num − Tnew .num − |Tnew .arr| − |Told .arr| = 1 + 2(1) + 0 = 3. In the second case, suppose the table does need to be doubled as a result of the insertion. Then the actual cost of insertion is |Tnew .num|, so the amortized cost of insertion is cb = c + Φ (Tnew ) − Φ (Told ) ¡ ¢ ¡ ¢ = Tnew .num + 2 · Tnew .num − |Tnew .arr| − 2 · Told .num − |Told .arr| ¡ ¢ ¡ ¢ ¡ ¢ = |Told .arr| + 1 + 2 Tnew .num − Told .num − |Tnew .arr| − |Told .arr| ¡ ¢ ¡ ¢ = |Told .arr| + 1 + 2(1) − 2 |Told .arr| − |Told .arr| = 3. In both cases, the amortized running time of insertion is 3. (If we had picked a less convenient potential function Φ, the amortized running time would probably be different in the two cases. The potential analysis would still be valid if done correctly, but the resulting amortized running times probably wouldn’t be as easy to work with.) 11.4.5 Exercise: Allowing deletions Suppose we wanted our dynamically resized table to support pops (deleting the most recently inserted element) as well as insertions. It would be reasonable to want the table to halve itself when fewer than half the slots in the array are occupied, so that the table only ever takes up at most twice the space it needs to. Thus, a first attempt at a pop functionality might look like this (with the T ABLE -I NSERT procedure unchanged): Algorithm: N AÏVE -P OP(T) 1 2 3 4 5 6 7 8 9 10 if T.num = 0 then error “Tried to pop from an empty table” r ← T.arr[T.num − 1] Unset T.arr[T.num − 1] T.num ← T.num − 1 if T.num ≤ 21 |T.arr| then ¥ ¦ arr0 ← new array of size 12 |T.arr| Migrate all items from T.arr to arr0 T.arr ← arr0 return r Lec 11 – pg. 7 of 8 Unfortunately, the worst-case running time of N AÏVE -P OP is not so good. Exercise 11.1. Given n sufficiently large, produce a sequence of n T ABLE -I NSERT, lookup, and/or ¡ ¢ N AÏVE -P OP operations which have total running time Θ n2 . Thus, in any amortized analysis, at least one of the operations must have amortized running time at least Ω (n). Exercise 11.2. Devise a way of supporting pops which still retains O(1) worst-case amortized running time for each operation, and do an amortized analysis to prove that this is the case. You may redefine T ABLE -I NSERT or add more data fields to T if you wish (and of course you will have to define T ABLE -P OP), but it must still remain the case that, for every possible sequence of operations, T only takes up O (T.num) space in memory at any given point in time. [Hint: An adversary can choose a series of operations for which T ABLE -I NSERT and N AÏVE -P OP run slowly by making sure T.num repeatedly crosses the critical thresholds for table doubling and table halving. (You must have used this fact somehow in your solution to Exercise 11.1.) So your first step will be to figure out how to place decisions about doubling and halving outside the control of your adversary.] Lec 11 – pg. 8 of 8 Lecture 12 Competitive Analysis Supplemental reading in CLRS: None 12.1 Online and Offline Algorithms For some algorithmic problems, the optimal response to a sequence of calls to a procedure depends not just on each input separately, but on the sequence of inputs as a whole. This is analogous to amortized analysis, in which the best running time analysis required us to consider the running time of an entire sequence of operations. When an algorithm expects to be given the entire sequence of inputs ahead of time, it is called an offline algorithm; when the algorithm is expected to give an answer after each individual input without knowledge of future inputs, it is called an online algorithm. Example. Traders are faced with the following algorithmic problem: given today’s stock prices (and the history of stock prices in the past), decide which stocks to buy and which stocks to sell. Obviously, this is an online problem for everybody. An attempt to treat the trading problem as an offline algorithmic problem is called insider trading. Example. Tetris is an online game. The offline version of Tetris would be akin to a jigsaw puzzle, and would not require quick thinking. Of course, we would expect players of the offline game to produce much better packings than players of the online game. Let’s lay out a general setup for analyzing online and offline algorithms. Assume we are given an algorithm A, along with a notion of “cost” C A . That is, if S is a sequence of inputs, then C A (S) is the cost of trusting A’s response to the inputs in S. For example, A might be a stock broker and C A (S) might be the net profit after a sequence of days on the stock market. In that case, clients will probably try to learn as much as the can about the function C A (for example, by observing its past outputs) when deciding whether to hire the stock broker. In general, no online algorithm will perform optimally on all sequences of inputs—in other words, the so-called “God’s algorithm” (which performs optimally on every sequence of inputs) will usually be impossible to describe in an online way. This is obvious in the stock market example: if you chose to buy today, what if prices drop tomorrow? If you chose to sell today, what if prices go up tomorrow? Although an online algorithm may not ever be able to match the performance of God’s algorithm, it may be possible for an online algorithm to perform almost as well as God’s algorithm on every input. Definition. An online algorithm A is said to be α-competitive1 (where α is a positive constant) if there exists a constant k such that, for every sequence S of inputs, we have C A (S) ≤ α · C O PT (S) + k, where O PT is the optimal offline algorithm, A . K . A . God’s algorithm. 12.2 Example: A Self-Organizing List Problem 12.1. Design a data structure L representing a list of n key–value pairs (with distinct keys), satisfying the following constraints: • The key–value pairs are stored in a linked list • There is only one supported operation: A CCESS(x), which returns the element with key x. (It is assumed that the input x is always one of the keys which occurs in L.) The implementation of A CCESS(x) has two phases: 1. Walk through the list until you find the element with key x. The cost of this phase is rankL (x), the rank of the element with key x. 2. Reorder the list by making a sequence of adjacent transpositions2 (in this lecture, the word “transposition” will always refer to adjacent transpositions). The cost of each transposition is 1; the sequence of transpositions performed is up to you, and may be empty. Try to choose a sequence of transpositions which optimizes the performance of A CCESS, i.e., for which A CCESS is α-competitive where α is minimal among all possible online algorithms. 12.2.1 Worst-case inputs No matter how we define our algorithm, an adversary (who knows what algorithm we are using) can always ask for the last element in our list, so that phase 1 always costs n, and the cost of a sequence S of such costly operations is at least ¡ ¢ C A (S) = Ω |S | · n , even disregarding the cost of any reordering we do. 12.2.2 Heuristic: Random inputs As a heuristic (which will seem particularly well-chosen in retrospect), let’s suppose that our inputs are random rather than worst-case. Namely, suppose the inputs are independent random variables, where the key x has probability p(x) of being chosen. It should be believable (and you can prove if 1 This is not the only notion of “competitive” that people might care about. We could also consider additive competitive- ness, which is a refinement of the notion considered here: We say an algorithm A is k-additively competitive if there exists a constant k such that C A (S) ≤ C O PT (S) + k for every sequence S of inputs. Thus, every k-additively competitive algorithm is 1-competitive in the sense defined above. 2 An adjacent transposition is the act of switching two adjacent entries in the list. For example, transposing the elements B and C in the list 〈 A, B, C, D 〉 results in the list 〈 A, C, B, D 〉. Lec 12 – pg. 2 of 5 you wish) that, when |S | is much larger than n, the optimal expected cost is achieved by immediately sorting L in decreasing order of p(x), and keeping that order for the remainder of the operations.3 In that case, the expected cost of a sequence S of A CCESS operations is E [C A (S)] = S X ¡ ¢ p(x) · rankP p(x) · |S | , x∈ L where P = { p(y) : y ∈ L}. In practice, we won’t be given the distribution p ahead of time, but we can still estimate what p would be if it existed. We can keep track of the number of times each key appears as an input, and maintain the list in decreasing order of frequency. This counting heuristic has its merits, but we will be interested in another heuristic: the move-to-front (MTF) heuristic, which works as follows: • After accessing an item x, move x to the head of the list. The total cost of this operation is 2 · rankL (x) − 1: – rankL (x) to find x in the list – rankL (x) − 1 to perform the transpositions.4 To show that the move-to-front heuristic is effective, we will perform a competitive analysis. Proposition 12.2. The MTF heuristic is 4-competitive for self-organizing lists. As to the original problem, we are not making any assertion that 4 is minimal; for all we know, there could exist α-competitive heuristics for α < 4.5 Proof. • • • • Let L i be MTF’s list after the ith access. Let L∗i be O PT’s list after the ith access. Let c i be MTF’s cost for the ith operation. Thus c i = 2 rankL i−1 (x) − 1. Let c∗i be O PT’s cost for the ith operation. Thus c∗i = rankL∗i−1 (x) + t i , where t i is the number of transpositions performed by O PT during the ith operation. We will do an amortized analysis. Define the potential function6,7 ¡ ¢ Φ (L i ) = 2 · # of inversions between L i and L∗i 3 The intuition is that the most likely inputs should be made as easy as possible to handle, perhaps at the cost of making some less likely inputs more expensive. 4 The smallest number of transpositions needed to move x to the head of the list is rank (x) − 1. These rank (x) − 1 L L transpositions can be done in exactly one way. Once they are finished, the relative ordering of the elements other than x is unchanged. (If you are skeptical of these claims, try proving them.) For example, to move D to the front of the list 〈 A, B, C, D, E 〉, we can perform the transpositions D ↔ C, D ↔ B, and D ↔ A, ultimately resulting in the list 〈D, A, B, C, E 〉. 5 The question of whether α = 4 is minimal in this example is not of much practical importance, since the costs used in this example are not very realistic. For example, in an actual linked list, it would be easy to move the nth entry to the head in O(1) time (assuming you already know the location of both entries), whereas the present analysis gives it cost Θ(n). Nevertheless, it is instructive to study this example as a first example of competitive analysis. 6 The notation x < L i y means rankL i (x) < rankL i (y). In other words, we are using L i to define an order x for every real number x). Depending on your programming language (and on whether you cast the edge capacities as integers or floating-point numbers), the symbol ∞ may or may not be supported natively. If it is not supported natively, you will have to either implement it or add extra code to the implementation of F ORD –F ULKERSON so that operations on edge capacities support Lec 13 – pg. 9 of 11 s1 t1 ∞ s2 ∞ ∞ ∞ s t2 ∞ ∞ t ∞ s3 t3 s4 Figure 13.4. We can allow multiple sources and sinks by adding a “master” source that connects to each other source via an edge of capacity ∞, and similarly for the sinks. s1 t1 = t2 s2 t3 s3 Figure 13.5. An instance of the multi-commodity flow problem with three commodities. Lec 13 – pg. 10 of 11 It turns out that the multi-commodity flow problem is NP-complete, even when k = 2 and all the edge capacities are 1. Thus, most computer scientists currently believe that there is no way to solve the multi-commodity flow problem in polynomial time, though it has not been definitively proved yet. both numbers and the symbol ∞. Lec 13 – pg. 11 of 11 Lecture 14 Interlude: Problem Solving Supplemental reading in CLRS: None This lecture was originally given as a pep talk before the take-home exam. In note form, this chapter will be light reading, a break in which we look back at the course material as veterans. 14.1 What to Bring to the Table In most technical undergraduate courses, questions are posed in such a way that all solutions are equally valid; if a student is able to come up with any solution to a given problem, she is considered to be in good shape. In 6.046, this is not the case. For most algorithmic problems, there is an obvious but inefficient solution; your task is to find a better solution. It is worth appreciating that this way of thinking may be new to you, and taking a moment to reflect metacognitively on strategies to help you solve algorithmic problems. There is no recipe for solving problems, but there are many things you can do which often, eventually, lead to good solutions. • Knowledge. In academic research, one will usually pick up the relevant background material over the course of working on a given problem. In an exam setting, you should come to the table with as complete an understanding of the background material as possible; you will not have time to learn many new concepts. In either case, you should strive to understand the relevant background as deeply as possible. You should be able to implement a given algorithm correctly in your computer language of choice. You should be able to answer CLRS-style exercises about the pseudocode: What would be the effect on performance if we changed this line? Would the algorithm still be correct? What would be the effect on performance if we used a different data structure? • Belief. Be optimistic. Remind yourself that you might solve this problem, even if it doesn’t look tractable at first. • Motivation. There is no substitute for a strong motivation to solve the problem at hand. Given that you’re trying to learn algorithms from a set of lecture notes, I probably don’t even need to tell you this. 14.2 How to Attack a Problem We’ll now discuss several tips that might help you to crack a problem open. We’ll keep the following concrete example in the back of our mind: Problem 14.1 (Bipartite Matching). In a group of n heterosexual people, each woman has a list of the men she is interested in marrying, and each man has a list of the women he is interested in marrying. Find an algorithm that creates the largest possible number of weddings. (Each person can marry at most one other person. In a couple, each member must be on the other member’s list.) 14.2.1 Understand the Problem • Note everything that is given. • Find some upper bounds and lower bounds. For example: – If the algorithm has to read the entire input, then the running time is at least Ω(n). (This is not always the case, however! For example, binary search runs in O(lg n) and does not read the entire input.) – If the output has to output k things, then the running time is at least Ω(k). – If the algorithm can be used to perform a comparison-based sort, then the running time is at least Ω(n lg n).1 – There are usually a finite number of possibilities for the output. In the case of our current problem, arrange the men in a line (by alphabetical order of last name, say). For each man who has a wife, place his wife behind him. Place the unmarried women at the front of the line (by alphabetical order of last name, say). In this way, each possible matching gives a different line. Thus, the total number of possible matchings is at most n!, the number of ways to arrange all n people in a line. – Usually there is an obvious algorithm that runs in exponential time. Sometimes it runs faster. – Sometimes it is easy to reduce the problem to another problem whose solution is wellknown. (Be careful, though—your special case may be amenable to a more efficient solution than the general problem!) – Whatever your solution is, its running time must of course fall within the lower and upper bounds you find. • Is there a useful diagram? 1 A comparison-based sorting algorithm is a sorting algorithm which makes no assumptions about the items being sorted except that there is a well-defined notion of “less than” and that, for objects a and b, we can check whether a < b in constant time. The fact we are using here is that any comparison-based sorting algorithm takes at least Ω(n lg n) time. The proof is as follows. Suppose we have a comparison-based sorting algorithm A which takes as input a list L = 〈a 1 , . . . , a n 〉 and outputs a list of indices I = 〈 i 1 , . . . , i n 〉 such that a i 1 ≤ · · · ≤ a i n . Let k be the number of comparison operations performed by A on a worst-case input of length n (where k depends on n). Because each comparison is a yes-or-no question, there are at most 2k possible sets of answers to the at most k comparison operations. Since the output of a comparison-based sorting algorithm must depend only on these answers, there are at most 2k possible outputs I. Of course, if A is correct, then all n! permutations of the indices 〈1, . . . , n〉 must be possible, so 2k ≥ n!. Thus k ≥ lg(n!). Finally, Stirling’s approximation tells us that lg(n!) = Θ(n lg n). Thus the running time of A is Ω(k) = Ω(n lg n). Lec 14 – pg. 2 of 5 • Is there a graph? • Try examples! 14.2.2 Try to Solve the Problem • Can we solve the problem under a set of simplifying assumptions? (What if every vertex has degree 2? What if the edge weights are all 1?) • Are there useful subproblems? • Will a greedy algorithm work? Divide-and-conquer? Dynamic programming? • Can this problem be reduced to a known problem? • Is it similar to a known problem? Can we modify the known algorithm to suit it to this version? • Can the problem be solved with convolution? If so, the fast Fourier transform will improve efficiency. • Is there any further way to take advantage of the given information? • Can randomization help? It is often better to find a fast algorithm with a small probability of error than a slower, correct algorithm. Remember that a hard problem usually cannot be solved in one sitting. Taking breaks and changing perspective help. In the original setting of a pre–take-home exam pep talk, it was important to remind students that “incubation time” is completely necessary and should be part of the exam schedule. Thus, students should start the test early—the away-time when they are not actively trying to solve test questions is crucial. 14.2.3 Work Out the Details • What is the algorithm? (Write the algorithm down in as much detail as possible—it is easy to end up with an incorrect algorithm or an incorrect running time analysis because you forgot about a few of what seemed like unimportant implementation details.) • Why is it correct? • What is its asymptotic running time? 14.2.4 Communicate Your Findings • Don’t be shy about writing up partial solutions, solutions to simplified versions of the problem, or observations. This is the most common sort of writing that occurs in academia: often, the complete solution to a problem is just the finishing touch on a large body of prior work by other scientists. And even more often, no complete solution is found—all we have are observations and solutions to certain tractable special cases. (In an exam setting, this sort of writing could win you significant partial credit.) Lec 14 – pg. 3 of 5 • To be especially easy on your reader, you can format your algorithm description as an essay, complete with diagrams and examples. The essay should begin with a high-level “executive summary” of how the algorithm works. Now that you have been thinking about the bipartite matching problem for a while, here is a solution: Algorithm: Let M be the set of men and let W be the set of women; let s and t be objects representing a source and a sink. Create a bipartite directed graph∗ G with vertex set M ∪ W ∪ { s, t}. Draw an edge from s to each man and an edge from each woman to t. Next, for each man m and each woman w, draw an edge from m to w if m is on w’s list and w is on m’s list. Give all edges capacity 1. The graph G is now a flow network with source s and sink t (see Figure 14.1). Note that there is a bijection between valid matchings and integer flows in G. In one direction, given a matching, we can saturate each edge between a man and a woman who have been matched. The flows out of s and into t are then uniquely determined. In the other direction, given an integer flow f in G, the fact that each edge out of s or into t has capacity 1 means that there is at most one edge out of each man and at most one edge into each woman. Thus, we obtain a valid matching by wedding a man m to a woman w if and only if f (m, w) = 1. The magnitude of f is equal to the number of couples, so the number of couples in any matching is bounded above by the maximum flow in G. This bound is tight if and only if there exists a maximum flow f ∗ in G which is an integer flow. Such a flow does exist, and is found by the Ford–Fulkerson algorithm, by Proposition 13.3. Since the number of couples is at most O(n), we have | f ∗ | = O(n), and the ¡ ¢ ¡ ¢ running time of the Ford–Fulkerson algorithm is O (E · | f ∗ |) = O n2 · n = O n3 . Actually, assuming each person has a relatively small list of acceptable spouses, we ought to be more granular about this bound: we have à ¶! n µ X size of the ith , E =O n+ person’s list i =1 and the running time is at most à à ¶!! n µ X size of the ith T =O n· n+ . person’s list i =1 ∗ A graph G = (V , E) is called bipartite if the vertex set V can be written as the union of two disjoint sets V = V1 t V2 such that every edge in E has one endpoint in V1 and one endpoint in V2 . 14.2.5 Reflect and Improve • Can we achieve a better running time?2 2 For the bipartite matching problem, there exist better algorithms than the one given above. Still, I love this solution. It shows how useful flow networks are, even in contexts that seem to have nothing to do with flow. There are many more examples in which this phenomenon occurs. Lec 14 – pg. 4 of 5 1 1 1 s 1 1 t 1 1 1 1 Figure 14.1. In this flow network, the nodes on the left (other than s) represent women and the nodes on the right (other than t) represent men. All edges have capacity 1. A maximum flow is highlighted in black; it shows that the maximal number of couples is 3. • Can randomization help? • Can amortization help? • Can we use less space? 14.3 Recommended Reading There are many great references and exercise books for algorithmic problem solving and mathematical problem solving in general. Two particular selections are: • George Pólya, How to Solve It • The Princeton Companion to Mathematics, VIII.1: “The Art of Problem Solving.” Another good source of algorithms practice is web sites designed to help people prepare for technical interviews. It is common for tech companies to ask algorithm questions to job applications during interviews; students who plan to apply for software jobs may especially want to practice. Even for those who aren’t looking for jobs, web sites like TechInterview.org can be a useful source of challenging algorthmic exercises. Lec 14 – pg. 5 of 5 Lecture 15 van Emde Boas Data Structure Supplemental reading in CLRS: Chapter 20 Given a fixed integer n, what is the best way to represent a subset S of {0, . . . , n − 1} in memory, assuming that space is not a concern? The simplest way is to use an n-bit array A, setting A[i] = 1 if and only if i ∈ S, as we saw in Lecture 10. This solution gives O(1) running time for insertions, deletions, and lookups (i.e., testing whether a given number x is in S).1 What if we want our data structure to support more operations, though? Perhaps we want to be able not just to insert, delete and lookup, but also to find the minimum and maximum elements of S. Or, given some element x ∈ S, we may want to find the successor and predecessor of x, which are the smallest element of S that is greater than x and the largest element of S that is less than x, respectively. On a bit array, these operations would all take time Θ(|S |) in the worst case, as we might need to examine all the elements of S. The van Emde Boas (vEB) data structure is a clever alternative solution which outperforms a bit array for this purpose: Operation I NSERT D ELETE L OOKUP M AXIMUM, M INIMUM S UCCESSOR, P REDECESSOR Bit array O(1) O(1) O(1) Θ(n) Θ(n) van Emde Boas Θ (lg lg n) Θ (lg lg n) Θ(lg lg n) O(1) Θ(lg lg n) We will not discuss the implementation of S UCCESSOR and P REDECESSOR; those will be left to recitation. 15.1 Analogy: The Two-Coconut Problem The following riddle nicely illustrates the idea of the van Emde Boas data structure. I am somewhat embarrassed to give the riddle because it shows a complete misunderstanding of materials science and botany, but it is a standard example and I can’t think of a better one. Problem 15.1. There exists some unknown integer k between 1 and 100 (or in general, between 1 and some large integer n) such that, whenever a coconut is dropped from a height of k inches or 1 This performance really is unbeatable, even for small n. Each operation on the bit array requires only a single memory access. n h h k h h h Figure 15.1. One strategy for the two-coconut problem is to divide the integers {1, . . . , n} into blocks of size h, and then use the first coconut to figure out which block k is in. Once the first coconut breaks, it takes at most h − 1 more drops to find the value of k. more, the coconut will crack, and whenever a coconut is dropped from a height of less than k inches, the coconut will be completely undamaged and it will be as if we had not dropped the coconut at all. (Thus, we could drop the coconut from a height of k − 1 inches a million times and nothing would happen.) Our goal is to find k with as few drops as possible, given a certain number of test coconuts which cannot be reused after they crack. If we have one coconut, then clearly we must first try 1 inch, then 2 inches, and so on. The riddle asks, what is the best way to proceed if we have two coconuts? An approximate answer to the riddle is given by the following strategy: Divide the n-inch range into b blocks of height h each (we will choose b and h later; see Figure 15.1). Drop the first coconut from height h, then from height 2h, and so on, until it cracks. Say the first coconut cracks at height b 0 h. Then drop the second coconut from heights (b 0 − 1)h + 1, (b 0 − 1)h + 2, and so on, until it cracks. p This method requires at most b + h − 1 = hn + h − 1 drops, which is minimized when b = h = n. Notice that, once the first coconut cracks, our problem becomes identical to the one-coconut version (except that instead of looking for a number between 1 and n, we are looking for a number between (b 0 − 1)h + 1 and b 0 h − 1). Similarly, if we started with three coconuts instead of two, then it would be a good idea to divide the range {1, . . . , n} into equally-sized blocks and execute the solution to the two-cononut problem once the first coconut cracked. Exercise 15.1. If we use the above strategy to solve the three-coconut problem, what size should we p choose for the blocks? (Hint: it is not n.) 15.2 Implementation: A Recursive Data Structure Just as in the two-coconut problem, the van Emde Boas data structure divides the range {0, . . . , n − 1} p p into blocks of size n, which we call clusters. Each cluster is itself a vEB structure of size n. In addition, there is a “summary” structure that keeps track of which clusters are nonempty (see Figure Lec 15 – pg. 2 of 6 vEB: size = 16 vEB of size 4 representing {0, 1, 3} summary = max = 13 min = 2 vEB of size 4 representing {2, 3} clusters = vEB of size 4 representing {1, 3} N IL vEB of size 4 representing {1} Figure 15.2. A vEB structure of size 16 representing the set {2, 3, 5, 7, 13} ⊆ {0, . . . , 15}. 15.2). The summary structure is analogous to the first coconut, which told us what block k had to lie in. 15.2.1 Crucial implementation detail The following implementation detail, which may seem unimportant at first, is actually crucial to the performance of the van Emde Boas data structure: Do not store the minimum and maximum elements in clusters. Instead, store them as separate data fields. Thus, the data fields of a vEB structure V of size n are as follows: • V .size – the size of V , namely n • V .max – the maximum element of V , or N IL if V is empty • V .min – the minimum element of V , or N IL if V is empty p • V .clusters – an array of size n which stores the clusters. For performance reasons, the value stored in each entry of V .clusters will initially be N IL; we will wait to build each cluster until we have to insert something into it. p • V .summary – a van Emde Boas structure of size n that keeps track of which clusters are nonempty. As with the entries of V .clusters, we initially set V .summary ← N IL; we do not build the recursive van Emde Boas structure referenced by V .summary until we have to insert something into it (i.e., until the first time we create a cluster for V ). 15.2.2 Insertions To simplify the exposition, in this lecture we use a model of computation in which it takes constant time to initialize any array (setting all entries equal to N IL), no matter how big the array is.2 Thus, 2 Of course, this is cheating; real computers need an initialization time that depends on the size of the array. Still, there are use cases for which this assumption is warranted. We can preload our vEB structure by creating all possible vEB Lec 15 – pg. 3 of 6 it takes constant time to create an empty van Emde Boas structure, and it also takes constant time to insert the first element into a van Emde Boas structure: Algorithm: V EB-F IRST-I NSERTION(V , x) V .min ← x 2 V .max ← x 1 Using this fact, we will show that the procedure V .I NSERT(x) has only one non–constant-time step. Say V .clusters[i] is the cluster corresponding to x (for example, if n = 100 and x = 64, then i = 6). Then: • If V .clusters[i] is N IL, then it takes constant time to create a vEB structure containing only the element corresponding to x. (For example, if n = 100 and x = 64, then it takes constant time to create a vEB structure of size 10 containing only 4.) We update V .clusters[i] to point to this new vEB structure. Thus, the only non–constant-time operation we have to perform is to update V .summary to reflect the fact that V .clusters[i] is now nonempty. • If V .clusters[i] is empty3 (we can check this by checking whether V .clusters[i].min = N IL), then it takes constant time to insert the appropriate entry into V .clusters[i]. So again, the only non–constant-time operation we have to perform is to update V .summary to reflect the fact that V .clusters[i] is now nonempty. • If V .clusters[i] is nonempty, then we have to make the recursive call V .clusters[i].I NSERT(x). However, we do not need to make any changes to V .summary. In each case, we find that the running time T of I NSERT satisfies the recurrence ¡p ¢ T(n) = T n + O(1). (15.1) As we will see in §15.3 below, the solution to this recurrence is T I NSERT = Θ (lg lg n). 15.2.3 Deletions Similarly, the procedure V .D ELETE(x) requires only one recursive call. To see this, let i be as above. Then: • First suppose V has no nonempty clusters (we can check this by checking whether V .summary is either N IL or empty; the latter happens when V .summary.min = N IL). – If x is the only element of V , then we simply set V .min ← V .max ← N IL. Thus, deleting the only element of a single-element vEB structure takes constant time; we will use this fact later. – Otherwise, V contains only two elements including x. Let y be the other element of V . If x = V .min, then y = V .max and we set V .min ← y. If x = V .max, then y = V .min and we set V .max ← y. substructures up front rather than using N IL for the empty ones. This way, after an initial O(n) preloading time, array initialization never becomes an issue. Another possibility is to use a dynamically resized hash table for V .clusters rather than an array. Each of these possible improvements has its share of extra details that must be addressed in the running time analysis; we have chosen to ignore initialization time for the sake of brevity and readability. 3 This would happen if V .clusters[i] was once nonempty, but became empty due to deletions. Lec 15 – pg. 4 of 6 • Next suppose V has some nonempty clusters. – If x = V .min, then we will need to update V .min. The new minimum takes only constant time to calculate, though: it is y, where y ← V .clusters[V .summary.min].min; in other words, it’s the smallest element in the lowest nonempty cluster of V . After making this update, we need to make a recursive call to V .clusters[V .summary.min].D ELETE(y) to remove the new value of V .min from its cluster. At this point, there are two possibilities: * y is not the only element in its cluster. Then V .summary does not need to be updated. * y is the only element in its cluster. Then we must make a recursive call to V .summary.D ELETE(V .summary.min) to reflect the fact that y’s cluster is now empty. However, since y was the only element in its cluster, it took constant time to remove y from its cluster: as we said above, it takes only constant time to remove the last element from a one-element vEB structure. Either way, there is only one non–constant-time step in V .D ELETE: a recursive call to p D ELETE on a vEB structure of size n. – By entirely the same argument (perhaps in mirror-image), we find that the case x = V .max has identical running time to the case x = V .min. – If x is neither V .min nor V .max, then we must delete x from its cluster and, if this causes x’s cluster to become empty, make a recursive call to V .summary.D ELETE to reflect this update. As above, the recursive call to V .summary.D ELETE will only happen when the deletion of x from its cluster took constant time. In each case, the only non–constant-time step in D ELETE is a single recursive call to D ELETE on a p vEB structure of size n. Thus, the running time T for D ELETE satisfies (15.1), which we repeat here for convenience: ¡p ¢ T(n) = T n + O(1). (copy of 15.1) Again, the solution is T D ELETE = Θ(lg lg n). 15.2.4 Lookups Finally, we consider the operation V .L OOKUP(x), which returns T RUE or FALSE according as x is or is not in V . The implementation is easy: First we check whether x = V .min, then whether x = V .max. If neither of these is true, then we recursively call L OOKUP on the cluster corresponding to x. Thus, the running time T of L OOKUP satisfies (15.1), and we have T L OOKUP = Θ(lg lg n). Exercise 15.2. Go through this section again, circling each step or claim that relies on the decision not to store V .min and V .max in clusters. Be careful—there may be more things to circle than you think! Lec 15 – pg. 5 of 6 15.3 Solving the Recurrence As promised, we now solve the recurrence (15.1), which we repeat here for convenience: T(n) = T 15.3.1 ¡p ¢ n + O(1). (copy of 15.1) Base case Before we begin, it’s important to note that we have to lay down a base case at which recursive structures stop occurring. Mathematically this is necessary because one often uses induction to prove solutions to recurrences. From an implementation standpoint, the need for a base case is obvious: how could a vEB structure of size 2 make good use of smaller vEB substructures? So we will lay down n = 2 as our base case, in which we simply take V to be an array of two bits. 15.3.2 Solving by descent The equation (15.1) means that we start pp on a structure of size n, then pass to a structure of size p 1/2 n = n , then to a structure of size n = n1/4 , and so on, spending a constant amount of time at each level of recursion. So the total running time should be proportional to the number of levels of ` recursion before arriving at our base case, which is the number ` such that n1/2 = 2. Solving for `, we find ` = lg lg n. Thus T(n) = Θ(lg lg n). 15.3.3 Solving by substitution Another way to solve the recurrence is to make a substitution which reduces it to a recurrence that we already know how to solve. Let ¡ ¢ T 0 (m) = T 2m . Taking m = lg n, (15.1) can be rewritten as T 0 (m) = T 0 (m/2) + O(1), which we know to have solution T 0 (m) = Θ(lg m). Substituting back n = 2m , we get T(n) = Θ(lg lg n). Lec 15 – pg. 6 of 6 Lecture 16 Disjoint-Set Data Structures Supplemental reading in CLRS: Chapter 21 (§21.4 is optional) When implementing Kruskal’s algorithm in Lecture 4, we built up a minimum spanning tree T by adding in one edge at a time. Along the way, we needed to keep track of the connected components of T; this was achieved using a disjoint-set data structure. In this lecture we explore disjoint-set data structures in more detail. Recall from Lecture 4 that a disjoint-set data structure is a data structure representing a dynamic collection of sets S = {S 1 , . . . , S r }. Given an element u, we denote by S u the set containing u. We will equip each set S i with a representative element rep[S i ].1 This way, checking whether two elements u and v are in the same set amounts to checking whether rep[S u ] = rep[S v ]. The disjoint-set data structure supports the following operations: • M AKE -S ET(u): Creates a new set containing the single element u. – u must not belong to any already existing set – of course, u will be the representative element initially • F IND -S ET(u): Returns the representative rep[S u ]. • U NION(u, v): Replaces S u and S v with S u ∪ S v in S. Updates representative elements as appropriate. In Lecture 4, we looked at two different implementations of disjoint sets: doubly-linked lists and trees. In this lecture we’ll improve each of these two implementations, ultimately obtaining a very efficient tree-based solution. 16.1 Linked-List Implementation Recall the simple linked-list implementation of disjoint sets that we saw in Lecture 4: 1 This is not the only way to do things, but it’s just as good as any other way. In any case, we’ll need to have some symbol representing each set S i (to serve as the return value of F IND -S ET); the choice of what kind of symbol to use is essentially a choice of notation. We have chosen to use a “representative element” of each set rather than create a new symbol. Si : rep[S i ] u1 u2 u3 uk ··· M AKE -S ET(u) – initialize as a lone node Θ(1) F IND -S ET(u) – walk left from u until you reach the head of S u Θ(n) worst-case U NION(u, v) Θ(n) worst-case – walk right (towards the tail) from u and left (towards the head) from v. Reassign pointers so that the tail of S u and the head of S v become neighbors. The representative is updated automatically. Do we really need to do all that walking? We could save ourselves some time by augmenting each element u with a pointer to rep[S u ], and augmenting rep[S u ] itself with a pointer to the tail. That way, the running time of F IND -S ET would be O(1) and all the walking we formerly had to do in U NION would become unnecessary. However, to maintain the new data fields, U NION(u, v) would have to walk through S v and update each element’s “head” field to point to rep[S u ]. Thus, U NION would still take O(n) time in the worst case. Perhaps amortization can give us a tighter bound on the running time of U NION? At first glance, it doesn’t seem to help much. As an example, start with the sets {1} , {2} , . . . , { n}. Perform U NION(2, 1), followed by U NION(3, 1), and so on, until finally U NION(n, 1), so that S 1 = {1, . . . , n}. In each call to U NION, we have to walk to the tail of S 1 , which is continually growing. The ith call to U NION has to walk through i − 1 elements, so that the total running time of the n − 1 U NION operations is ¡ 2¢ Pn−1 i =1 (i − 1) = Θ n . Thus, the amortized cost per call to U NION is Θ(n). However, you may have noticed that we could have performed essentially the same operations more efficiently by instead calling U NION(1, 2) followed by U NION(1, 3), and so on, until finally U NION(1, n). This way, we never have to perform the costly operation of walking to the tail of S 1 ; we only ever have to walk to the tails of one-element sets. Thus the running time for this smarter sequence of n − 1 U NION operations is Θ(n), and the amortized cost per operation is O(1). The lesson learned from this analysis is that, when performing a U NION operation, it is best to always merge the smaller set into the larger set, i.e., the representative element of the combined set should be chosen equal to the representative of the larger constituent—that way, the least possible amount of walking has to occur. To do this efficiently, we ought to augment each S i with a “size” field, which we’ll call S i .weight (see Figure 16.1). It turns out that the “smaller into larger” strategy gives a significant improvement in the amortized worst-case running time of the U NION operation. We’ll show that the total running time of any sequence of U NION operations on a disjoint-set data structure with n elements (i.e., in which M AKE -S ET is called n times) is O(n lg n). Thus, the running time of m operations, n of which are M AKE -S ET operations, is ¡ ¢ O m + n lg n . To start, focus on a single element u. We’ll show that the total amount of time spent updating u’s “head” pointer is O(lg n); thus, the total time spent on all U NION operations is O(n lg n). When u is added to the structure via a call to M AKE -S ET, we have S u .weight = 1. Then, every time S u merges with another set S v , one of the following happens: Lec 16 – pg. 2 of 7 Si : head u1 u2 tail rep[S i ] u3 weight = 3 Figure 16.1. A linked list augmented with data fields for the head, the tail, and the size (weight). • S u .weight > S v .weight. Then no update to u’s “head” pointer is needed. • S v .weight ≥ S u .weight. Then, we update u’s “head” pointer. Also, in this case, the value of S u .weight at least doubles. Because S u .weight at least doubles every time we update u’s “head” pointer, and because S u .weight can only be at most n, it follows that the total number of times we update u’s “head” pointer is at most lg n. Thus, as above, the total cost of all U NION operations is O(n lg n) and the total cost of any sequence of m operations is O (m + n lg n). Exercise 16.1. With this new augmented structure, do we still need the list to be doubly linked? Which pointers can we safely discard? 16.2 Forest-of-Trees Implementation In addition to the linked-list implementation, we also saw in Lecture 4 an implementation of the disjoint-set data structure based on trees: u1 Si : rep[S i ] u4 u2 u5 u3 u6 M AKE -S ET(u) – initialize new tree with root node u Θ(1) F IND -S ET(u) – walk up tree from u to root Θ(height) = Θ(lg n) best-case U NION(u, v) O(1) + 2T F IND -S ET – change rep[S v ]’s parent to rep[S u ] The efficiency of the basic implementation hinges completely on the height of the tree: the shorter the tree, the more efficient the operations. As the implementation currently stands, the trees could Lec 16 – pg. 3 of 7 u v v.rank u.rank Figure 16.2. Union by rank attempts to always merge the shorter tree into the taller tree, using rank as an estimate (always an overestimate) of height. be unbalanced and F IND -S ET could take as long as Θ(n) in the worst case. However, this behavior can be dramatically improved, as we will see below. 16.2.1 Union by rank When we call U NION(u, v), rep[S v ] becomes a child of rep[S u ]. Merging S v into S u in this way results in a tree of height © ª max height[S u ], height[S v ] + 1 (16.1) (why?). Thus, the way to keep our trees short is to always merge the shorter tree into the taller tree. (This is analogous to the “smaller into larger” strategy in the linked-list implementation.) To help us do this, we will introduce a new data field called rank. If u is the root of a tree (i.e., if u = rep[S u ]), then u.rank will be an upper bound on the height of S u .2 In light of (16.1), the pseudocode for U NION will be as follows (see Figure 16.2): Algorithm: U NION( u e, ve) 1 2 3 4 5 6 7 8 9 u ← F IND -S ET( u e) v ← F IND -S ET(ve) if u.rank = v.rank then u.rank ← u.rank + 1 v.parent ← u else if u.rank > v.rank then v.parent ← u else u.parent ← v The following lemma shows that U NION preserves the fact that rank is an upper bound on height. Lemma 16.1. Suppose initially the following hold: • S u 6= S v 2 If union by rank is the only improvement we use, then u.rank will actually be the exact height of S . But in general, u we wish to allow other improvements (such as path compression) to decrease the height of S u without having to worry about updating ranks. In such cases, the upper bound provided by u.rank may not be tight. Lec 16 – pg. 4 of 7 • S u has height h 1 and S v has height h 2 • rep[S u ].rank = r 1 and rep[S v ].rank = r 2 • h 1 ≤ r 1 and h 2 ≤ r 2 . Suppose we then call U NION(u, v), producing a new set S = S u ∪ S v . Let h be the height of S and let r = rep[S].rank. Then h ≤ r. Proof. First, suppose r 1 > r 2 . Then S v has been merged into S u and r = r 1 . By (16.1), we have © ª h = max h 1 , h 2 + 1 © ª ≤ max r 1 , r 2 + 1 = r1 = r. A similar argument shows that h ≤ r in the case that r 2 > r 1 . Finally, suppose r 1 = r 2 . Then S v has been merged into S u and r = r 1 + 1 = r 2 + 1, so © ª h = max h 1 , h 2 + 1 © ª ≤ max r 1 , r 2 + 1 = r2 + 1 = r. It turns out that the rank of a tree with k elements is always at most lg k. Thus, the worst-case performance of a disjoint-set forest with union by rank having n elements is M AKE -S ET F IND -S ET U NION O(1) Θ(lg n) Θ(lg n). Exercise 16.2. Amortization does not help this analysis. Given sufficiently large n and given m which is sufficiently large compared to n, produce a sequence of m operations, n of which are M AKE -S ET operations (so the structure ultimately contains n elements), whose running time is Θ(m lg n). Exercise 16.3. Above we claimed that the rank of any tree with k elements is at most lg k. Use induction to prove this claim. (You may assume that U NION is the only procedure that modifies ranks. However, you should not assume anything about the height of a tree except that it is less than the rank.) What is the base case? 16.2.2 Path compression The easiest kind of tree to walk up is a flat tree, where all non-root nodes are direct children of the root (see Figure 16.3). The idea of path compression is that, every time we invoke F IND -S ET and walk up the tree, we should reassign parent pointers to make each node we pass a direct child of the root (see Figure 16.4). This locally flattens the tree. With path compression, the pseudocode for F IND -S ET is as follows: Lec 16 – pg. 5 of 7 u1 u4 u2 u5 u6 u3 Figure 16.3. In a flat tree, each F IND -S ET operation requires us to traverse only one edge. u1 u1 u2 u4 u3 u5 u6 u2 u3 u7 u6 u9 u8 u7 u4 u9 u5 u8 Figure 16.4. With path compression, calling F IND -S ET (u 8 ) will have the side-effect of making u 8 and all of its ancestors direct children of the root. Algorithm: F IND -S ET(u) 1 2 3 4 5 6 7 8 A←; n←u while n is not the root do A ← A ∪ { n} n ← n.parent for each x ∈ A do x.parent ← n return n What data structure should we use for A? In an ideal world, where n can truly be arbitrarily large, we would probably want A to be a dynamically doubled array of the kind discussed in Lecture 10. In real life, however, some assumptions can be made. For example, if you have less than a petabyte (1024 TB, or 253 bits) of available memory, then the rank (and therefore the height) of any tree is at ¡ ¢ most lg 253 = 53, and it would be slightly more efficient to maintain A as a static array of size 53 (with an end-marking sentinel value, perhaps). It can be shown that, with path compression (but not union by rank), the running time of any sequence of n M AKE -S ET operations, f F IND -S ET operations, and up to n − 1 U NION operations is ¡ ¡ ¢¢ Θ n + f 1 + log2+ f /n n . 16.2.3 Both improvements together The punch-line of this lecture is that, taken together, union by rank and path compression produce a spectacularly efficient implementation of the disjoint-set data structure. Lec 16 – pg. 6 of 7 Theorem 16.2. On a disjoint-set forest with union by rank and path compression, any sequence of m operations, n of which are M AKE -S ET operations, has worst-case running time ¡ ¢ Θ mα(n) , where α is the inverse Ackermann function. Thus, the amortized worst-case running time of each operation is Θ (α(n)). If one makes the approximation α(n) = O(1), which is valid for literally all conceivable purposes, then the operations on a disjoint-set forest have O(1) amortized running time. The proof of this theorem is in §21.4 of CLRS. You can read it if you like; it is not essential. You might also be interested to know that, in a 1989 paper, Fredman and Saks proved that Θ (α(n)) is the fastest possible amortized running time per operation for any implementation of the disjoint-set data structure. The inverse Ackermann function α is defined by © ª α(n) = min k : A k (1) ≥ n , where (k, j) 7→ A k ( j) is the Ackermann function. Because the Ackermann function is an extremely rapidly growing function, the inverse Ackermann function α is an extremely slowly growing function (though it is true that limn→∞ α(n) = ∞). The Ackermann function A (at least, one version of it) is defined by ( A k ( j) = j+1 ( j +1) A k−1 ( j) if k = 0 (that is, A k−1 iterated j + 1 times) for k ≥ 1. Some sample values of the Ackermann function are A 1 (1) = 3 A 1 ( j) = 2 j + 1 A 2 (1) = 7 A 2 ( j) = 2 j+1 ( j + 1) A 3 (1) = 2047 A 4 (1) À 1080 . By current estimates, 1080 is roughly the same order of magnitude as the number of particles in the observable universe. Thus, even if you are a theoretical computer scientist or mathematician, you will still most likely never end up considering a number n so large that α(n) > 4. Exercise 16.4. Write down a sequence of operations on the disjoint-set forest with union by rank and path compression (including M AKE -S ET operations) which cause a taller tree to be merged into a shorter tree. Why do we allow this to happen? Lec 16 – pg. 7 of 7 Lecture 17 Complexity and NP-completeness Supplemental reading in CLRS: Chapter 34 As an engineer or computer scientist, it is important not only to be able to solve problems, but also to know which problems one can expect to solve efficiently. In this lecture we will explore the complexity of various problems, which is a measure of how efficiently they can be solved. 17.1 Examples To begin, we’ll review three problems that we already know how to solve efficiently. For each of the three problems, we will propose a variation which might not be so easy to solve. • Flow. Given a flow network G with integer capacities, we can efficiently find an integer flow that is optimal using the Edmonds–Karp algorithm. • Multi-Commodity Flow. In §13.3, we considered a varation of network flow in which there are k commodities which need to simultaneously flow through our network, where we need to send at least d i units of commodity i from source s i to sink t i . • Minimum Cut. Given an undirected weighted graph G, we can efficiently find a cut (S, V \ S) of minimum weight. Exercise 17.1. Design an efficient algorithm to solve the minimum cut problem. (Hint in this footnote.1 ) • Maximum Cut. What if we want to find a cut of maximum weight? • Minimum Spanning Tree. Given an undirected weighted graph G = (V , E), we know how to efficiently find a spanning tree of minimum weight. A spanning tree of G is a subgraph G 0 = (V 0 , E 0 ) ⊆ G such that – G 0 is connected and contains no cycles – V0 = V. 1 Hint: Turn G into a flow network and find a maximum flow. In Theorem 13.7, we saw how to turn a maximum flow into a minimum cut. • Steiner Tree. What if, instead of requiring V 0 = V , we require V 0 = S for some given subset S ⊆ V? Each of these three variations is NP-hard. The typical attitude towards NP-hard problems is Don’t expect an efficient algorithm for this problem. However, I prefer the following more optimistic interpretation: If you find an efficient algorithm for this problem, you will get $1 million.2 Unlike the pessimistic interpretation, the above statement is 100% factual. The so-called P vs. NP problem is one of seven important open research questions for which Clay Mathematics Institute is offering a $1 million prize.3 17.2 Complexity So far in the course, we have been ignoring the low-level details of the mathematical framework underlying our analyses. We have relied on the intuitive notion that our ideal computer is “like a real-world computer, but with infinite memory”; we have not worried about explicitly defining what a “step” is. The fact of the matter is that there are many reasonable models of computation which make these notions explicit. Historically the first one was the Turing machine, invented by Alan Turing, considered by many to be the founding father of computer science. The model we have been using throughout the course is similar to a Turing machine; the main difference is that a “step” on our ideal computer closely resembles a processor cycle on a modern computer, whereas the “steps” of a Turing machine involve sequentially reading and writing to the so-called tape which represents its memory. For this lecture, you won’t have to go to the trouble of working with a particular model of computation in full detail; but it is worth noting that such details are important in theoretical computer science and should not be regarded as a triviality. In what follows, we will define the complexity classes P and NP. Before doing so, we will need a couple more definitions: Definition. A decision problem is a computation problem to which the answer is either “yes” or “no.” In mathematical language, we can think of a decision problem as a function whose domain is the set of possible input strings4 and whose range is {0, 1} (with 0 meaning “no” and 1 meaning “yes”). Definition. A complexity class is simply a set of decision problems. Most of the problems we have considered so far in the course are not decision problems but rather search problems—they ask not just whether a solution exists, but also what the solution is. Given a search problem, we can derive decision problems which ask yes-or-no questions about the solution; for example, we might ask: 2 I should warn you though, most computer scientists believe that it is not possible to find one. (In other words, most computer scientists believe that P 6= NP.) 3 One of the questions has already been solved, so currently there are six prizes remaining. 4 Strings over what alphabet? A typical choice is {0, 1} (i.e., binary); another possible choice is the ASCII alphabet. The main reason the choice of alphabet matters is that it determines what “an input of size n” is. The number 255 has size 8 in binary, size 1 in ASCII, and size 255 in unary. An algorithm whose running time is linear with respect to a unary input would be exponential with respect to a binary input. Lec 17 – pg. 2 of 7 Problem 17.1. Given a graph G and an integer k, is there a spanning tree of size less than k? For most real-world applications, search problems are much more important than decision problems. So why do we restrict our attention to decision problems when defining complexity classes? Here are a few reasons: • The answer to a decision problem is simple. • The answer to a decision problem is unique. (A search problem might have multiple correct answers, e.g., a given graph might have multiple minimum spanning trees.) • A decision problem which asks about the answer to a search problem is at most as difficult as the search problem itself. For example, if we can find a minimum spanning tree efficiently, then we can certainly also solve Problem 17.1 efficiently. 17.2.1 P and NP The existence of many different models of computation is part of the reason for the following definition: Definition. “Efficient” means “polynomial-time.” An algorithm is polynomial-time if there exists a constant r such that the running time on an input of size n is O (n r ). The set of all decision problems which have polynomial-time solutions is called P. Polynomial time is the shortest class of running times that is invariant across the vast majority of reasonable, mainstream models of computation. To see that shorter running times need not be invariant, consider the following program: Read the first bit of memory 2 Read the nth bit of memory 1 In our model of computation, which has random access to memory, this would take constant time. However, in a model of computation with only serial access to memory (such as a Turing machine), this would take linear time. It is true, though, that any polynomial-time program in our model is also polynomial-time on a Turing machine, and vice versa. Search problems have the property that, once a solution is found, it can be verified quickly. This verifiability is the motivation for the complexity class NP. Definition. A decision problem P is in NP if there exists a polynomial-time algorithm A(x, y) such that, for every input x to the problem P, P(x) = 1 ⇐⇒ there exists some y such that A(x, y) = 1. The string y is called a witness or certificate; the algorithm A is called a verifier or a nondeterministic algorithm.5 For example, in Problem 17.1, the witness y could be the spanning tree itself—we can certainly verify in polynomial time that a given object y is a spanning tree of size less than k. 5 The abbrevation NP stands for “nondeterministic polynomial-time.” The reason for this name is as follows. Imagine receiving x (the input to the problem P) but leaving the choice of y unspecified. The result is a set of possible running times of A, one for each choice of y. The problem P is in NP if and only if at least one of these possible running times is bounded by a polynomial p (| x|) in the size of x. (The choice of y can depend on x, but p cannot depend on x.) Lec 17 – pg. 3 of 7 A B Figure 17.1. A Cook reduction of A to B is a program that would run in polynomial time on an oracle machine—that is, a Turing machine equipped with an oracle for B. If some day we find an efficient algorithm for B, then we can create an efficient algorithm for A by replacing the oracle with a subroutine. Proposition 17.2. P ⊆ NP. Proof. Given a decision problem P, view P as a function whose domain is the set of strings and whose range is {0, 1}. If P can be computed in polynomial time, then we can just take A(x, y) = P(x). In this case, the verifier just re-solves the entire problem. The converse to the above proposition is a famous open problem: Problem 17.3 (P vs. NP). Is it true that P = NP? The vast majority of computer scientists believe that P 6= NP, and so the P vs. NP problem is sometimes called the P 6= NP problem. If it were true that P = NP, then lots of problems that seem hard would actually be easy: one such example is the algorithm search problem described in §17.3. 17.2.2 Polynomial-time reductions It is possible to know that one problem is “at least as hard” as another without knowing exactly how hard each problem is. If problem A can be polynomial-time reduced to problem B, then it stands to reason B is at least as hard as A. There are two different notions of polynomial-time reduction, which we now lay out. Definition. We say that the decision problem A is Karp-reducible to the decision problem B if there exists a polynomial-time computable function f such that, for every input x, A(x) = 1 ⇐⇒ B ( f (x)) = 1. Definition. Problem A is Cook-reducible to problem B if there exists an algorithm which, given an oracle6 for B, solves A in polynomial time. (See Figure 17.1.) Note that a Karp reduction is also a Cook reduction, but not vice versa. Historically, Karp reductions and Cook reductions correspond to different traditions. Definition. A problem is NP-hard if every problem in NP can be Cook-reduced to it. 6 An oracle for B is a magical black-box which solves any instance of problem B in one step. A Cook reduction is a program which is allowed to do any of the normal things a computer program does, and is also allowed to query the oracle. Lec 17 – pg. 4 of 7 Problems with efficient solution (or equivalently, NP-complete problems) Every decision problem is NP-hard Problems with efficient solution P = NP NP-complete NP-hard P 6= NP Figure 17.2. Left: If P = NP, then every decision problem is NP-hard (why?). Right: If P 6= NP. Definition. A problem is NP-complete if it is both NP-hard and in NP. Using the notion of NP-completeness, we can make an analogy between NP-hardness and big-O notation: O( f (n)) on the order of at most f (n) Θ( f (n)) tightly on the order of f (n) Ω( f (n)) on the order of at least f (n) in NP at most as hard as an NP-complete problem NP-complete exactly as hard as any other NP-complete problem NP-hard at least as hard as an NP-complete problem Showing that a given problem is in NP is relatively straightforward (or at least, it is clear what the proof should look like): one must give a polynomial-time verifier. By contrast, it is much less clear how one might show that a given problem is NP-hard. One strategy is to reduce another NP-hard problem to it. But this strategy only works if one already knows certain problems to be NP-hard; it could not have been used as the first ever proof that a problem was NP-hard. That first proof was accomplished by Cook in 1971: Theorem 17.4 (Cook’s Theorem). 3SAT is NP-complete. Problem 17.5 (3SAT). Given a set of atomic statements x1 , . . . , xn , a literal is either an atom x i or its negation ¬ x i . A clause is the disjunction (“or”) of a finite set of literals. The 3SAT problem asks, given a propositional formula ϕ(x1 , . . . , xn ) which is the “and” of finitely many clauses of length 3, does there exist an assignment of either T RUE or FALSE to each x i which makes ϕ(x1 , . . . , xn ) evaluate to T RUE? For example, one instance of 3SAT asks whether there exists an assignment of x1 , . . . , x4 which makes the proposition (x ∨ x2 ∨ ¬ x3 ) ∧ (¬ x1 ∨ x3 ∨ x4 ) ∧ (¬ x1 ∨ ¬ x2 ∨ ¬ x3 ) ∧ (x1 ∨ x3 ∨ ¬ x4 ) . | 1 {z } | {z } | {z } | {z } clause clause clause clause evaluate to T RUE. The answer to this instance happens to be “yes,” as shown by the assignment x2 7→ T RUE, x1 , x3 , x4 7→ FALSE. Lec 17 – pg. 5 of 7 17.3 Example: Algorithm Search In this section, we solve an algorithmic problem whose output is itself an algorithm: Problem 17.6. Given an algorithmic problem P and a function T(n), find an algorithm which runs in time at most T(n), if such an algorithm exists. Output not just a description of the algorithm, but also a proof of correctness and running time analysis. Proposition 17.7. There exists a “meta-algorithm” which solves Problem 17.6 (but runs forever if no algorithm exists). If P = NP, then the running time of this meta-algorithm is polynomial in the size of the shortest possible output. The formats of both the input and the output of this algorithm deserve some explanation. By “an algorithmic problem,” we mean a mathematical description of the relationship between the input and the output. For example, we could express the notion of “a flow with magnitude at least k” in symbols as   ∀ u, v ∈ V 0 ≤ f (u, v) ≤ c(u, v)           X X    ∀ u ∈ V \ { s, t} f (u, v) = f (v, u)  . v∈V v∈V       X X      f (s, v) − f (v, s) ≥ k   v∈V v∈V In the case of a decision problem, this would be more along the lines of ( 1 if there exists f satisfying the above 0 otherwise. Regarding the output of our meta-algorithm, we need a format for “a description of the algorithm.” One possibility is a text file containing the source code of an implementation in a particular programming language. Next, we need a format for “a proof of correctness and running time analysis.” For this we appeal to a machine-checkable proof language of the sort used by theorem-checking software. One such software suite is Coq.7 If you’ve never seen Coq before, I suggest you check it out! The key to the proof of Proposition 17.7 is to consider the following problem: Problem 17.8. Given an algorithmic problem P, a function T(n) and an integer k, does there exist a solution to Problem 17.6 in which the output has length at most k? Problem 17.8 is in NP, so if P = NP then there exists a polynomial-time algorithm for it. Proof of Proposition 17.7. Consider the following search problem: Problem 17.9. Given an algorithmic problem P, a function T(n), an integer k, and a prefix string s, does there exist a solution to Problem 17.8 in which the output starts with s? 7 http://coq.inria.fr/ Lec 17 – pg. 6 of 7 Problem 17.9 is in NP, so if P = NP then there exists a polynomial-time algorithm for it. Thus, if we let |P | denote the length of the description of P, |T | the length of the definition of the function T, and | s| the length of s, then there exist constants a, b, c, d such that Problem 17.9 has a solution A(P, T, k, s) which runs in time ³¯ ¯ ¯ ¯ ¯ ¯d ´ a b O ¯ P ¯ ¯ T ¯ k c ¯ s¯ . We can use this algorithm to solve Problem 17.8 by probing longer and longer prefixes s. For example, supposing we use the 26-letter alphabet A, . . . , Z for inputs, we would proceed as follows: • Run A(P, T, k, A). If the answer is 0, then run A(P, T, k, B). Proceed in this way until you find a prefix which returns 1. If all 26 letters return 0, then the answer to Problem 17.8 is “no.” • Otherwise, let’s say for the sake of a concrete example that A(P, T, k, F) = 1. Then, run A(P, T, k, FA), A(P, T, k, FB), etc., until you find a two-letter prefix that returns 1. • Proceeding in this way, the prefix s will eventually become the answer to Problem 17.8. The above procedure solves Problem 17.8. Because the length of s ranges from 1 to at most k, the running time is ´ ³ ¯ ¯ ¯ ¯ ¯ ¯a ¯ ¯ b ¯ ¯a ¯ ¯ b a b O 26 ¯P ¯ ¯T ¯ k c 1d + 26 ¯P ¯ ¯T ¯ k c 2d + · · · + 26 ¯P ¯ ¯T ¯ k c k d ³ ´´ ³¯ ¯ ¯ ¯ a b =O ¯P ¯ ¯T ¯ k c 1d + 2d + · · · + k d ³¯ ¯ ¯ ¯ ´ a b =O ¯P ¯ ¯T ¯ k c k d +1 ³¯ ¯ ¯ ¯ ´ a b =O ¯P ¯ ¯T ¯ k c+d +1 . Thus, we have done more than just show that Problem 17.8 is in P. We have shown that, for some ¡ ¢ constants α, β, γ, there exists a solution to Problem 17.8 which runs in time O |P |α |T |β kγ and also returns the full algorithm-and-proofs, not just 0 or 1. To conclude the proof, our meta-algorithm is to run the above procedure for k = 1, 2, 3, . . .. If the shortest possible algorithm-and-proofs has length `, then the running time of the meta-algorithm is ³¯ ¯ ¯ ¯ ³¯ ¯ ¯ ¯ ´ ¯ ¯α ¯ ¯β ¯ ¯α ¯ ¯β ´ α β α β O ¯P ¯ ¯T ¯ 1γ + ¯P ¯ ¯T ¯ 2γ + · · · + ¯P ¯ ¯T ¯ `γ = O ¯P ¯ ¯T ¯ `γ+1 . Lec 17 – pg. 7 of 7 Lecture 18 Polynomial-Time Approximations Supplemental reading in CLRS: Chapter 35 except §35.4 If you try to design algorithms in the real world, it is inevitable that you will come across NP-hard problems. So what should you do? 1. Maybe the problem you want to solve is actually less general than the NP-hard problem. • Perhaps the input satisfies certain properties (bounded degree, planarity or other geometric structure, density/sparsity. . . ) • Remember, if you are able to reduce your problem to an NP-hard problem, that doesn’t mean your problem is NP-hard—it’s the other way around! 2. Maybe you can settle for an approximation. 3. Maybe an exponential algorithm is not so bad. • There might be an O (c n ) algorithm with c ≈ 1. In this lecture we will focus on item 2: approximation algorithms. As the name suggests, an approximation algorithm is supposed to return an answer that is close to the correct answer. There are several types of approximations one might consider, each based on a different notion of closeness. Definition. Suppose A is an algorithm for the optimization problem Π. (So the answer to Π is a number which we are trying to maximize or minimize, e.g., the weight of a spanning tree.) Given input x, we denote the output of A by A(x), and the optimal answer by Π(x). The following are senses in which A can be an “approximation algorithm” to the problem Π: • The most common type of approximation is multiplicative. For a positive number α ∈ R (which may depend on the input), we say that A is a multiplicative α-approximation (or simply an α-approximation) to Π if, for every input x, ½ ¾ αΠ(x) ≤ A(x) ≤ Π(x) if Π is a maximization problem . Π(x) ≤ A(x) ≤ αΠ(x) if Π is a minimization problem Of course, we must have 0 < α ≤ 1 for a maximization problem and α ≥ 1 for a minimization problem. If someone talks about a “2-approximation” to a maximization problem, they are probably referring to what here would be called a 21 -approximation. • For a positive number β ∈ R, we say that A is a additive β-approximation to Π if, for every input x, Π(x) − β ≤ A(x) ≤ Π(x) + β. • There are lots of other possibilities.1 There is no need to worry about memorizing different definitions of approximation. When someone talks about an “approximation,” it is their job to define precisely in what way it approximates the answer. In this lecture we will see approximation algorithms for three NP-hard problems: Problem Vertex Cover Set Cover Partition Approximation factor 2 ln n + 1 (where n is the total number of elements) 1 + ² for any given parameter ² > 0 Note that the (1 + ²)-approximation algorithm to Partition is a so-called polynomial-time approximation scheme (PTAS). Given a parameter ² > 0, the PTAS generates a polynomial-time algorithm which approximates the Partition problem to a factor of 1 + ². Naturally, the polynomials will get larger as we decrease ², since we are asking for a better approximation. The approximations to Vertex Cover and Set Cover in this lecture are conjectured to be optimal, in the sense that giving a better approximation would be NP-hard.2,3 There exist better approximations to Partition, though.4 18.1 Vertex Cover NP-hard problem (Vertex Cover). Given an undirected graph G = (V , E), a vertex cover of G is a subset V 0 ⊆ V such that, for every edge (u, v) ∈ E, we have either u ∈ V 0 or v ∈ V 0 . The Vertex Cover problem asks, given a graph G, what is the smallest possible vertex cover? Approximation Algorithm: 1 2 3 4 5 6 V0 ← ; while E is nonempty do Pick any (u, v) ∈ E V 0 ← V 0 ∪ { u, v} Remove from E all edges touching u or v Remove u and v from V 1 As a final example, we might combine multiplicative and additive approximations to define an “α, β–affine-linear approximation,” in which α−1 Π(x) − β ≤ A(x) ≤ αΠ(x) + β for all inputs x, where α ≥ 1. 2 If the so-called “Unique Games Conjecture” holds, then it is NP-hard to α-approximate Vertex Cover for any constant α < 2. It was recently proved that it is NP-hard to (c ln n)-approximate Set Cover for any constant c < 1. 3 Proofs of such optimality typically use the probabilistically checkable proofs (PCP) theorem. The PCP theorem states that a problem is in NP if and only if solutions to the problem can be verified (with high probability) by a certain kind of randomized algorithm. 4 There exists a PTAS for Partition in which the running time depends only polynomially on 1/²; this sort of PTAS is known as a fully polynomial-time ¡ ¢approximation scheme (FPTAS). As a consequence, for any given integer r there exists a polynomial-time 1 + 1/n r -approximation to Partition. Moreover, there exists a pseudo-polynomial–time algorithm which solves the Partition problem exactly. A pseudo-polynomial–time algorithm is one whose running time is polynomial in the size of the input and the numeric value of the output (in this case, the “cost” of the optimal partition). Lec 18 – pg. 2 of 6 Figure 18.1. Example run of the approximation algorithm for Vertex Cover. The edges chosen by line 3 are highlighted in orange and the chosen vertices are circled. The algorithm returns a cover of size 4, whereas the optimal covering has size 3. Figure 18.2. The Set Cover problem gives us a covering of a universe set U and asks us to find a small subcovering. How many of the sets in this picture are required to cover the array of dots? It is clear that this algorithm always returns a vertex cover. Moreover, if G is given in the adjacency-list format, then the algorithm runs in linear time. Claim. The cover returned by this algorithm is at most twice the size of an optimal cover. Proof. Let E 0 denote the set of edges that get chosen by line 3. Notice that no two edges in E 0 share an endpoint. Thus, every vertex cover must contain at least one endpoint of each edge in E 0 . Meanwhile, the output of our algorithm consists precisely of the two endpoints of each edge in E 0 . 18.2 Set Cover NP-hard problem (Set Cover). Given a set U and subsets S 1 , . . . , S m ⊆ U such that indices I ⊆ {1, . . . , m}, with | I | minimal, such that [ S i = U. i∈ I Approximation Algorithm: while U is not empty do 2 Pick the largest subset S i 3 Remove all elements of S i from U and from the other subsets 4 return a list of the sets we chose 1 Lec 18 – pg. 3 of 6 Sm i =1 S i = U, find The running time for a good implementation of this algorithm is à ! m ¯ ¯ X ¯S i ¯ O i =1 (see Exercise 35.3-3 of CLRS). Claim. The above algorithm gives a (ln |U | + 1)-approximation. Proof. Assume the optimal cover has size k. Let U i denote the value of U (our universe set) after i iterations of the loop. Clearly, for all i, the set U i can be covered by k sets. So one of these k sets | | | | must contain at least Uki elements, and therefore the set chosen on line 2 has size at least Uki . Thus, we have ¯ ¯ ¡ ¢¯ ¯ ¯U i+1 ¯ ≤ 1 − 1 ¯U i ¯ k for all i. This implies that ¯ ¯ ¡ ¢ ¯ ¯ ¯U i ¯ ≤ 1 − 1 i ¯U0 ¯ k for all i. Since 1 − 1k ≤ e−1/k , it follows that ¯ ¯ ¯ ¯ ¯U i ¯ ≤ e− i/k ¯U0 ¯ . In particular, letting n = |U0 |, we have ¯ ¯ ¯Uk(ln n+1) ¯ < 1 ¯ ¯ ¯Uk(ln n+1) ¯ = 0. =⇒ Thus, the loop exits after at most k(ln n + 1) iterations. 18.3 Partition NP-hard problem (Partition). Given a sorted list of numbers s 1 ≥ s 2 ≥ · · · ≥ s n , partition the indices {1, . . . , n} into two sets {1, . . . , n} = A t B such that the “cost” ( ) X X max si, si i ∈B i∈ A is minimized. Example. For the partition 12 10 9 7 4 3 2 it is optimal to take A = {1, 2, 7} and B = {3, 4, 5, 6}, so that X i∈ A s i = 24 and X s i = 23. i ∈B Lec 18 – pg. 4 of 6 Approximation Algorithm: 1 2 3 4 5 6 7 8 9 10 11 B ² is a parameter: for a given ², the running time is polynomial in n when we view ² as a constant m ← b1/²c By brute force, find an optimal partition {1, . . . , m} = A 0 t B0 for s 1 , . . . , s m A ← A0 B ← B0 for i ← m + 1 to n do P P if j∈ A s j ≤ j∈B s j then A ← A ∪ { i} else B ← B ∪ { i} return 〈 A, B〉 Note that this algorithm always returns a partition. The running time of a reasonable implementation is5 ¢ ¡ Θ 2m + n = Θ(n). Claim. The above algorithm gives a (1 + ²)-approximation. Proof. Without loss of generality,6 assume X si ≥ X si. i ∈B 0 i∈ A 0 Let H= n 1X si. 2 i=1 Notice that solving the partition problem amounts to finding a set A ⊆ {1, . . . , n} such that P P as close to H as possible. Moreover, since i∈ A s i + i∈B s i = 2H, we have ( ) X X D max si, si = H + , 2 i ∈B i∈ A where P i∈ A s i is ¯ ¯ ¯X X ¯¯ ¯ D=¯ s − s ¯. ¯ i ∈ A i i ∈B i ¯ • Case 1: X s i > H. i∈ A 0 P P In that case, the condition on line 7 is always false, seeing as i∈ A 0 s i > i∉ A 0 s i . So A = A 0 and B = {1, . . . , n} \ A 0 . In fact, approximation aside, I claim that this must be the optimal P P partition. To see this, first note that i∈B s i < H < i∈ A 0 s i . If there were a better partition 5 Of course, the dependence on ² is horrible. In a practical context, we would have to choose ² small enough to get a useful approximation but big enough that our algorithm finishes in a reasonable amount of time. This would be helped if we replaced line 3 with a more efficient subroutine. 6 We can have our algorithm check whether this equation holds, and switch the roles of A 0 and B0 if it doesn’t. Lec 18 – pg. 5 of 6 {1, . . . , n} = A ∗ t B∗ , then taking A 00 = A ∗ ∩{1, . . . , m} and B00 = B ∩{1, . . . , m} would give a partition of {1, . . . , m} in which ( ) X X X max si, si < si, i ∈B00 i ∈ A 00 i∈ A 0 contradicting the brute-force solution on line 3. • Case 2: X s i ≤ H. i∈ A 0 Note that, if A ever gets enlarged (i.e., if the condition on line 7 is ever true), then we have P D ≤ s i (where i is as in line 6). And if A is never enlarged, then it must be the case that i∈B s i P never exceeded i∈ A s i until the very last iteration of lines 6–10, in which s n is added to B. In that case we have D ≤ s n (why?). Either way, we must have D ≤ s m+1 and consequently ( max ) X i∈ A si, X si H+ = i ∈B D 2 ≤ H + 21 s m+1 . Thus, cost of the algorithm’s output cost of the algorithm’s output ≤ optimal cost H ≤ H + 12 s m+1 = 1+ = 1+ H 1 2 s m+1 H 1 2 s m+1 P 1 m+1 2 i =1 s i s m+1 = 1 + Pm+1 i =1 s i ≤ 1+ s m+1 (m + 1) s m+1 = 1+ 1 m+1 < 1 + ². Lec 18 – pg. 6 of 6 Lecture 19 Compression and Huffman Coding Supplemental reading in CLRS: Section 16.3 19.1 Compression As you probably know at this point in your career, compression is a tool used to facilitate storing large data sets. There are two different sorts of goals one might hope to achieve with compression: • Maximize ease of access, manipulation and processing • Minimize size—especially important when storage or transmission is expensive. Naturally, these two objectives are often at odds with each other. In this lecture we will focus on the second objective. In general, data cannot be compressed. For example, we cannot losslessly represent all m-bit strings using (m − 1)-bit strings, since there are 2m possible m-bit strings and only 2m−1 possible (m − 1)-bit strings. So when is compression possible? • If only a relatively small number of the possible m-bit strings appear, compression is possible. • If the same “long” substring appears repeatedly, we could represent it by a “short” string. • If we relax the requirement that every string have a unique representation, then compression might work but make “similar” strings identical. 19.1.1 Lossless and lossy Initial Compressed Reconstituted Compress Decompress Document −−−−−−−−−−−→ Document −−−−−−−−−−−−→ Document D C D0 In lossless compression, we require that D = D 0 . This means that the original document can always be recovered exactly from the compressed document. Examples include: • Huffman coding • Lempel–Ziv (used in gif images) In lossy compression, D 0 is close enough but not necessarily identical to D. Examples include: • mp3 (audio) • jpg (images) • mpg (videos) 19.1.2 Adaptive and non-adaptive Compression algorithms can be either adaptive or non-adaptive. • Non-adaptive – assumes prior knowledge of the data (e.g., character frequncies). • Adaptive – assumes no knowledge of the data, but builds such knowledge. 19.1.3 Framework For the remainder of this lecture, we consider the following problem: Input: Known alphabet (e.g., English letters a,b,c,. . . ) Sequence of characters from the known alphabet (e.g., “helloworld”) We are looking for a binary code—a way to represent each character as a binary string (each such binary string is called a codeword). Output: Concatenated string of codewords representing the given string of characters. 19.1.4 Fixed-length code In a fixed-length code, all codewords have the same length. For example, if our alphabet is © ª a, b, c, d, e, f, g, h , then we can represent each of the 8 characters as a 3-bit string, such as   b c d e f g h   a ↓ ↓ ↓ ↓ ↓ ↓ ↓ ↓ .   000 001 010 011 100 101 110 111 Such a code is easy to encode and decode: ¯ ¯ ¯ baba = 001¯000¯001¯000 Just as in DNA encoding and decoding, it is important to keep track of register: the deletion or insertion of a single bit into the binary sequence will cause a frame shift, corrupting all later characters in the reconstituted document. 19.1.5 Run-length encoding Imagine you are given the string of bits 111 |{z} 000 |{z} 11 00000 000000 | {z } . | {z } |{z} 6 3 3 2 5 Rather than create codewords, we could simply store this string of bits as the sequence 〈6, 3, 3, 2, 5〉. This strategy is used in fax-machine transmission, and also in jpeg. Lec 19 – pg. 2 of 7 19.1.6 Variable-length code If we allow our codewords to have different lengths, then there is an obvious strategy: Use shorter codewords for more frequent characters; use longer codewords for rarer characters. For example, consider a six-letter alphabet with the following character frequencies: Character Frequency Codeword a 45% 0 (1 bit) b 13% 101 c d 12% 16% 100 111 (3 bits) e f 9% 5% 1101 1100 (4 bits) Using this code, the average number of bits used to encode 100 characters is (45)1 + (13)3 + (12)3 + (16)3 + (9)4 + (5)4 = 224. Compare this to a 3-bit fixed-length code, in which it would take 300 bits to encode 100 characters. Notice that I didn’t use 1, 01, or 10 as codewords, even though they would have made the encoding shorter. The reason for this is that we want our code to be uniquely readable: it should never be ambiguous as to how to decode a compressed document.1 One standard way to make uniquely readable codes is to use prefix coding: no codeword occurs as a prefix (initial substring) of another codeword. This allows unambiguous, linear-time decoding: 101 111 1100 0 100 1101 d b a f c e Prefix coding means that we can draw our code as a binary tree, with the leaves representing codewords (see Figure 19.1). 19.2 The Huffman Algorithm The problem of §19.1.3 amounts to the following. We are given an alphabet {a i } with frequencies © ª © ª f (a i ) . We wish to find a set of binary codewords C = c(a 1 ), . . . , c(a n ) such that the average number of bits used to represent the data is minimized: B(C) = n X ¯ ¯ f (a i ) ¯ c(a i )¯ . i =1 Equivalently, if we represent our code as a tree T with leaf nodes a 1 , . . . , a n , then we want to minimize B(T) = n X f (a i ) d(a i ), i =1 where d(a i ) is the depth of a i , which is also equal to the number of bits in the codeword for a i . The following algorithm, due to Huffman, creates an optimal prefix tree for a given set of characters C = {a i }. Actually, the Huffman code is optimal among all uniquely readable codes, though we don’t show it here. 1 For an example of non-unique readibility, suppose we had assigned to “d” the codeword 01 rather than 111. Then the string 0101 could be decoded as either “ab” or “dd.” Lec 19 – pg. 3 of 7 100% 0 1 45% 55% a=0 0 1 25% 0 30% 1 12% 0 13% c = 100 1 14% 16% b = 101 d = 111 0 1 5% 9% f = 1100 e = 1101 Figure 19.1. Representation of a binary code as a binary tree. Algorithm: H UFFMAN -T REE(C) 1 2 3 4 5 6 7 8 9 10 19.2.1 n ← |C | Q ← C B a min–priority queue keyed by frequency for i ← 1 to n − 1 do Allocate new node z z.left ← x ← E XTRACT-M IN(Q) z.right ← y ← E XTRACT-M IN(Q) z.freq ← x.freq + y.freq Q.I NSERT(z) B Return the root of the tree return E XTRACT-M IN(Q) Running time Initially, we build a minqueue Q with n elements. Next, the main loop runs n − 1 times, with each iteration consisting of two E XTRACT-M IN operations and one I NSERT operation. Finally, we call E XTRACT-M IN one last time; at this point Q has only one element left. Thus, the running time for Q a binary minheap would be Θ (n) | {z } B UILD -Q UEUE + Θ (n lg n) + | {z } Loop O(1) | {z } = Θ (n lg n) . E XTRACT-M IN If instead we use a van Emde Boas structure for Q, we achieve the running time Θ (n lg lg n) . Lec 19 – pg. 4 of 7 f: 5% e: 9% c: 12% b: 13% c: 12% b: 13% 14% c: 12% 25% c: 12% d: 16% a: 45% a: 45% 25% d: 16% e: 9% f: 5% a: 45% e: 9% f: 5% 14% d: 16% b: 13% 30% a: 45% 14% b: 13% f: 5% d: 16% e: 9% a: 45% 55% 25% 30% c: 12% 14% b: 13% f: 5% d: 16% e: 9% 100% a: 45% 55% 25% c: 12% 30% 14% b: 13% f: 5% d: 16% e: 9% Figure 19.2. Example run of the Huffman algorithm. The six rows represent the state of the graph each of the six times line 3 is executed. Lec 19 – pg. 5 of 7 x y a y a b T a b x b T0 y x T 00 Figure 19.3. Illustration of the proof of Lemma 19.1. 19.2.2 Correctness The Huffman algorithm is a greedy algorithm: at each stage, we merge together the two nodes of lowest frequency. Lemma 19.1 (CLRS Lemma 16.2). Suppose x, y are the two most infrequent characters of C (with ties broken arbitrarily). Then there exists an optimal prefix code for C with codewords for x and y of the same length and differing only in the last bit. (In other words, there exists an optimal tree in which x and y are siblings.) Proof. Let T be an optimal tree. Let a be a leaf of maximal depth. If a has no sibling, then deleting a from the tree (and using a’s parent to represent the character formerly represented by a) produces a better code, contradicting optimality. So a has a sibling b, and since a has maximum depth, b is a leaf. So a and b are nodes of maximal depth, and without loss of generality we can say a.freq ≤ b.freq. Also without loss of generality, say x.freq ≤ y.freq. Switch the leaves a and x; call the resulting tree T 0 (see Figure 19.3). Then B(T 0 ) ≤ B(T), seeing as x.freq ≤ a.freq. So it must be the case that x.freq = a.freq and T 0 is optimal as well. Similarly, switch the leaves b and y in T 0 ; call the resulting tree T 00 . It must be the case that y.freq = b.freq and T 00 is optimal. In T 00 , the leaves x and y are siblings. Lemma 19.2 (CLRS Lemma 16.3). Given an alphabet C, let: • • • • x, y be the two most infrequent characters in C (with ties broken arbitrarily) z be a new symbol, with z.freq ← x.freq + y.freq ¡ ¢ C 0 = C \ { x, y} ∪ { z} T 0 be an optimal tree for C 0 . Then T is an optimal tree for C, where T is a copy of T 0 in which the leaf for z has been replaced by an internal node having x and y as children. Proof. Suppose T is not optimal. By Lemma 19.1, let T be an optimal tree for C in which x and y are siblings. Delete x and y from T , and label their parent (now a leaf) with the symbol z. Call the Lec 19 – pg. 6 of 7 resulting tree T 0 . Notice that T 0 is a tree for C 0 , and furthermore, B(T 0 ) = B(T ) − x.freq − y.freq < B(T) − x.freq − y.freq = B(T 0 ). This contradicts the optimality of T 0 . We conclude that T must have been optimal in the first place. Corollary 19.3. The Huffman algorithm is correct. Proof sketch. Let T be the tree produced by the Huffman algorithm. Construct a new tree U as follows. Initially, let U consist of just one vertex. Then, perform a series of transformations on U. In each transformation, append two children to what was previously a leaf node. In this way, we can eventually transform U into T. Moreover, Lemma 19.2 guarantees that, at each stage, U is an optimal tree for the alphabet formed by its leaves. Lec 19 – pg. 7 of 7 Lecture 20 Sublinear-Time Algorithms Supplemental reading in CLRS: None If we settle for approximations, we can sometimes get much more efficient algorithms. In Lecture 18, we saw polynomial-time approximations to a few NP-hard problems. In what follows, we will concern ourselves with sublinear-time (that is, o(n)-time) approximation algorithms to problems whose exact solutions require at least linear time. We will see two examples: 1. Estimating the number of connected components of a graph G 2. Estimating the size of a minimum spanning tree of a graph G, given that all of G’s edge weights are integers between 1 and z for some given constant z.1 20.1 Estimating the Number of Connected Components In this section we exhibit a polynomial-time approximation scheme for counting the connected components of a graph. Given a graph G (in adjacency-list format) with n vertices and parameters ², δ > 0, we exhibit a randomized approximation algorithm which, with probability 1 − δ, provides an additive ² n-approximation. That is, ¸ ¶¯ ¶ µ · ¯µ ¯ ¯ output of correct number of ¯ ¯ > ²n − Pr ¯ connected components ¯ algorithm ≤ δ. (20.1) The running time of the algorithm is Poly (1/², lg(1/δ)).2 Given a vertex v, let m v denote the number of vertices in v’s connected component. Take a moment to convince yourself of the following lemma: Lemma 20.1. The number of connected components in G is X 1 . v∈V m v 1 The archaic Greek letter z (digamma) stood for the sound /w/. Although digamma fell out of use by the time of Classical Greek, it is still occasionally used in Greek numerals (which are used in Greece in approximately the same way that we use Roman numerals) to represent the number 6. 2 This notation means that the running time is polynomial in 1/² and lg(1/δ); i.e., there exist positive constants r , r 1 2 ³³ ´ ³ ´ ´ such that the running time of the algorithm is O r2 1 r1 lg δ1 . ² The idea of the algorithm is to approximate m v by a quantity m e v which can be computed in constant time. Then, for a set K of k randomly chosen vertices, we have µ ¶ X 1 X 1 n X 1 number of connected = ≈ ≈ . components ev k v∈ K m ev v∈V m v v∈V m Algorithm: A PPROX -#CC (G, ², δ) ³ ´ 1 1 1 For some k = Θ 2 lg δ , pick k vertices v1 , . . . , v k at random ² 2 for i ← 1 to k do 3 Set ½ ¾ 2 m e v i ← min m v i , ² B computed using breadth-first search 4 return the value of k nX 1 k i=1 m e vi The key line to examine is line 3. Exercise 20.1. (i) The running time of a breadth-first search depends on the size of the graph. So how can there be a bound on the running time of line 3 that doesn’t depend on n? (ii) What is the running time of line 3? Show that the total running time of the algorithm is ´ ³ ´ ³ O ²12 k = O ²14 lg δ1 . 20.1.1 Correctness We will prove (20.1) in two steps. First we will prove ¯ ¯ ¯X 1 X 1 ¯¯ ²n ¯ − ; ¯ ¯ ≤ ¯v∈V m e v v∈V m v ¯ 2 (20.2) then we will prove ¯ # "¯ ¯n X k X 1 ¯¯ 1 ²n ¯ ≤ δ. Pr ¯ − ¯ ≥ ¯ k i=1 m e vi v∈V m ev¯ 2 Combining these two, we obtain ¯ ¯ ¯n X k X 1 ¯¯ 1 ¯ − ¯ ¯ ¯ k i=1 m e v i v∈V m v ¯ ¯ ¯ ¯ ¯ ¯n X ¯X 1 k X 1 ¯¯ X 1 ¯¯ 1 ¯ ¯ ≤¯ − − ¯ + ¯ ¯, ¯ k i=1 m ¯v∈V m e v i v∈V m ev¯ e v v∈V m v ¯ and with probability 1 − δ, ··· ≤ ²n 2 + ²n 2 = ² n, which is (20.1). Lec 20 – pg. 2 of 6 (20.3) Proof of (20.2). This follows from the fact that 0 ≤ ² 1 1 1 1 − < ≤ ¡2¢ = m e v mv m ev 2 ² for each v ∈ V . Proof of (20.3). We use Hoeffding’s inequality, a relative of the Chernoff bound which is stated as P follows. Given independent real-valued random variables X 1 , . . . , X k , let Y = 1k ki=1 X i . Suppose a, b ∈ R are constants such that always3 a ≤ X i ≤ b for each i. Then Hoeffding’s inequality states that for any η > 0, we have ¶ µ ¯ h¯ i −2kη2 ¯ ¯ Pr ¯Y − E [Y ]¯ ≥ η ≤ 2 exp . (b − a)2 1 1P 1 We take X i = m e v , which gives E [Y ] = n v∈V m e v ; we take a = 0 and b = 1 and η = ²/2. Then Hoeffding’s i equality becomes ¯ # "¯ ¯1 X ³ ³ 2 ´´ k 1 1 X 1 ¯¯ ² ¯ ≤ 2 exp −2k ²4 . − Pr ¯ ¯ ≥ ¯ k i=1 m e vi n v∈V m ev¯ 2 ³ ´ lg(1/δ) δ) Thus, for a suitable k = Θ ²2 (namely, k = 2 ln(2/ ), we have ²2 ¡ ¢ · · · ≤ 2 exp − ln δ2 = δ. This is equivalent to (20.3). 20.2 Estimating the Size of a Minimum Spanning Tree In this section, we exhibit an algorithm which solves the following problem: Input: An undirected weighted graph G = (V , E, w) with n vertices all edge weights are integers in the set {1, . . . , z}, where z ≥ 2 is a given parameter all vertices have degree at most d given in adjacency-list format Parameters ², δ > 0 Output: A number t such that, with probability at least 1 − δ, (1 − ²) w∗ ≤ t ≤ (1 + ²) w∗ , where w∗ is the weight of a minimum spanning tree. ¡ ¢ The running time of the algorithm is Poly 1² , lg δ1 , z, d . 3 Or at least, with probability 1. Lec 20 – pg. 3 of 6 20.2.1 Motivation Imagine running Kruskal’s MST algorithm on a graph G = (V , E, w) whose edge weights are all integers from the set {1, . . . , z}. The procedure would look like this: 1 2 3 4 5 6 T ←; for i ← 1 to z do while there exists an edge of weight i which has its endpoints in different connected components do add the edge to T if |T | = n − 1 then return T The number of connected components in T starts at n and decreases by 1 every time we add an edge. This leads to the following insight, which will be crucial for us: Observation 20.2. Let G (i) = (V , E (i) ), where E (i) = { e ∈ E : w(e) ≤ i }, and let T (i) be the restriction of T to G (i) . Let c(i) be the number of connected components in G (i) . Lines 3–6 (plus induction) guarantee that the number of connected components in T (i) is also c(i) . Moreover, the number of connected components in T (i) is equal to ³ ´ n − # edges in T (i) (this is true for any forest with n vertices). Thus, we have ¶ µ ´ ³ # edges in T with = # edges in T (i) = n − c(i) . weight at most i Observation 20.2 plus some clever algebra lead to the following lemma: Lemma 20.3. With G (i) , T (i) and c(i) as above, we have w(T) = n − z + zX −1 c(i) . i =1 Proof. Let ¶ µ # edges it T with Ai = weight exactly i and ¶ µ # edges in T with . Bi = weight at least i Then µ # edges in T with B i = |T | − weight at most i − 1 ³ ´ = (n − 1) − n − c(i−1) = c(i−1) − 1. Now, the clever algebra trick is to notice that w(T) = z X i =1 i · Ai = z X i =1 Lec 20 – pg. 4 of 6 Bi. ¶ (Make sure you see how this works.) Completing the computation, w(T) = z X Bi i =1 = z ³ X c(i−1) − 1 ´ i =1 = −z + zX −1 c(i) i =0 zX −1 = n−z+ c(i) . i =1 Algorithm: MST-A PPROX(G, ², δ, z, d) for i ← 1 to z − 1 do 2 B Let G (i) denote the subgraph of G consisting of those edges whose weight is at most i ³ ´ δ 3 ĉ(i) ← A PPROX -#CC G (i) , 2²z , z 4 return the value of zX −1 |G.V | − z + ĉ(i) 1 i =1 The tricky part here is that we cannot compute G (i) and store it in memory, as that would take Ω (n + G.E) time. So how does line 3 work? We must modify A PPROX -#CC to so as to ignore any edges of weight greater than i. However, this modification forces us to reconsider the running time of A PPROX -#CC, and is the reason for the dependence on d. Exercise 20.2. (i) Suppose we modify A PPROX -#CC so as to ignore any edges of weight greater than i. Use an example to show that, if we treat ², δ and z as constants but do not allow any dependence on d, then the breadth-first search on line 3 of A PPROX -#CC has worst-case running time Ω(n). (Hint: Every time we ignore an edge, it takes one step.) (ii) Despite part (i), we can put a good bound on the running time of the modified version of A PPROX -#CC if we allow the bound to depend on d. Show that the modified breadth-first δ search in line 3 of A PPROX -#CC (with 2²z and z standing in for what in A PPROX -#CC are ³ ´ dz denoted ² and δ) takes O ² time, and thus that the total running time of A PPROX -MST is O 20.2.2 ³ d z4 ²3 ´ lg z . δ Correctness The return value of A PPROX -MST is n−z+ zX −1 ĉ(i) , i =1 Lec 20 – pg. 5 of 6 while in §20.2.1 we showed that w∗ = n − z + zX −1 c(i) . i =1 Thus, to show correctness, we need to show that ¯ "¯ # ¯ ¯zX zX −1 −1 ¯ ¯ ĉ(i) ¯ > ²w∗ ≤ δ. Pr ¯ c(i) − ¯ ¯ i=1 i =1 By the union bound, it suffices to show that ·¯ ¯ ¯ ¯ Pr ¯ c(i) − ĉ(i) ¯ > ² z−1 w δ ¸ ∗ ≤ z−1 for each i. Now, because each edge has weight at least 1, we know that w∗ ≥ n − 1. So it suffices to show that ·¯ ¸ ¯ δ ² ¯ ¯ Pr ¯ c(i) − ĉ(i) ¯ > (n − 1) ≤ . z−1 z−1 The correctness of A PPROX -#CC guarantees that ·¯ ¯ ¯ ¯ Pr ¯ c(i) − ĉ(i) ¯ > so we are fine as long as ² z−1 ² ¸ 2z (n − 1) ≥ n ² 2z ≤ n, which holds whenever n > 1. That is good enough. Lec 20 – pg. 6 of 6 δ z , Lecture 21 Clustering Supplemental reading in CLRS: None Clustering is the process of grouping objects based on similarity as quantified by a metric. Each object should be similar to the other objects in its cluster, and somewhat different from the objects in other clusters. Clustering is extremely useful; it is of fundamental importance in data analysis. Some applications include • • • • Scientific data from a wide range of fields Medical data, e.g., for patient diagnosis Identifying patterns of consumer behavior Categorizing music, movies, images, genes, etc. Clustering is conceptually related to • Unsupervised learning – the notion that objects which produce similar measurements may share some intrinsic property. • Dimensionality reduction. Depending on the application, we may have to carefully choose which metric to use out of many possible metrics. Or, we may have to apply a transformation to our data before we can get good clustering. But ultimately, the (admittedly vague) problem we are trying to solve is this: Figure 21.1. Left: Some data sets are relatively straightforward to cluster; most people would cluster these objects in basically the same way. Right: It’s not always so easy, though. How would you make two clusters from this data set? Input: Instance data D = {~x1 ,~x2 , . . . ,~xn } Desired number of clusters, k Distance metric1 d(~x i ,~x j ) Output: Assignment of instance data D to clusters C = {C 1 , . . . , C k }. 21.1 Hierarchical Agglomerative Clustering The main idea of hierarchical agglomerative clustering is to build up a graph representing the cluster set as follows: • Initially, each object (represented as a vertex) is in its own cluster. • Each time an edge is added, two clusters are merged together. • Stops when we have k clusters. The following is an implementation of hierarchical agglomerative clustering known as singlelinkage clustering. Algorithm: S INGLE -L INKAGE -C LUSTER(D, d, k) 1. Let H be an undirected graph with one vertex for each object and no edges. © ª 2. Sort the set of unordered pairs { u, v} : u, v ∈ D, u 6= v by distance: ¡ ¡ ¢¢ d (pair 1) ≤ d (pair 2) ≤ · · · ≤ d pair n2 . 3. Loop from i = 1 to ¡ n¢ 2 : • If the two members of pair i are in different connected components: – Merge their clusters. – Join pair i with an edge in H. • Exit the loop and halt when there are only k components left. Although we have couched the description of this algorithm in the language of graph theory, the use of a graph data structure is not essential—all we need is a disjoint-set data structure to store the connected components. The advantage of computing the graph is that the graph gives us more information about the influence each object had on cluster formation. 21.1.1 Running time ¡ ¢ ¡ ¢ In step 1, we initialize the graph and make n calls to M AKE -S ET. Next, there are n2 = Θ n2 ¡ 2 ¢ ¡ ¢ unordered pairs of objects; sorting them in step 2 takes Θ n lg n2 = Θ n2 lg n time. Finally, each iteration of the loop in step 3 makes two calls to F IND -S ET and at most one call to U NION. The loop ¡ ¢ is iterated at most O n2 times. ¡ ¢ Thus, the total number of operations performed on the disjoint-set data structure is n + O n2 = ¡ 2¢ O n . If we use a good implementation of the disjoint-set data structure (such as a disjoint-set forest 1 What we call here a “distance metric” corresponds to the notion of a metric space in mathematics. That is, our distance metric is a symmetric, positive-definite scalar-valued function of two arguments which satisfies the triangle inequality. Lec 21 – pg. 2 of 7 B (C ) Figure 21.2. Illustration of the quantity β (C ), which is a measure of the amount of space between clusters. ¡ ¢ with union by rank and path compression), these operations will take less time than the Θ n2 lg n sorting. Thus, the total running time of S INGLE -L INKAGE -C LUSTER is ¡ ¢ Θ n2 lg n . 21.1.2 Correctness Discussion This is where we would usually argue for the correctness of our algorithm. In the current case, there is no precise notion of correctness because we didn’t state exactly what properties our clusters should have—just that there are k of them and they represent some sort of nearness-based grouping. So instead, we’ll discuss the properties of this clustering and consider other possibilities. • This algorithm is essentially a special case of Kruskal’s MST algorithm. For k = 1, running S INGLE -L INKAGE -C LUSTER is exactly the same as running Kruskal’s algorithm on the complete graph whose edges are weighted by distance. • For k > 1, the graph produced is an MST with the k − 1 heaviest edges removed. • This algorithm maximizes the “spacing” between clusters. More precisely, the minimum distance between a pair of clusters is maximized, in the sense that the quantity © ª β(C ) = min d (~x,~y) : ~x, ~y not in the same cluster (21.1) is maximized (see Figure 21.2), as we show below. Proposition 21.1. Let C = {C 1 , . . . , C k } be the output of S INGLE -L INKAGE -C LUSTER. Then the quantity β(C ), as defined in (21.1), is maximized. Proof. Let d ∗ = β(C ). Then d ∗ is the distance between the first pair of vertices not considered by the algorithm. All edges in the graph H were considered, so their weights are all ≤ d ∗ . © ª Consider a different clustering C 0 = C 10 , . . . , C 0k . There must exist some C r ∈ C which is not a subset of any C 0s ∈ C 0 . (Make sure you can see why.2 ) So we can choose some ~x,~y ∈ C r and some s such that ~x ∈ C 0s and ~y ∉ C 0s . Because C r is a connected component of H, there exists a path ~x ~y in 2 What would it mean if every C were a subset of some C 0 ? This should not sit comfortably with the fact that C 0 6= C . r s Lec 21 – pg. 3 of 7 C 0t C 0s Cr ~z ~x ~ w ~y Figure 21.3. Illustration of the proof of Proposition 21.1. ~ ) which connects a vertex in C 0s to a vertex in H. Somewhere in that path there must be an edge (~z, w 0 some other cluster C t (see Figure 21.3). Thus ¡ ¢ ~ ) ≤ d∗, β(C 0 ) ≤ distance between C 0s and C 0t ≤ d(~ z, w ~ ) is an edge in H. since (~z, w 21.1.3 Other distance criteria The procedure S INGLE -L INKAGE -C LUSTER is designed to always merge together the two “closest” clusters, as determined by the distance criterion ¡ ¢ d min C i , C j = min ~x∈C i ,~y∈C j d (~x,~y) . Other forms of hierarchical agglomerative clustering use different distance criteria, such as3 ¡ ¢ d max C i , C j = max ~x∈C i ,~y∈C j d (~x,~y) X 1 d (~x,~y) |C i | · |C j | ~x∈C i ,~y∈C j à ! ¡ ¢ 1 X 1 X ~x, ~y . d centroid C i , C j = d |C i | ~x∈C i |C j | ~y∈C j ¡ ¢ d mean C i , C j = Note that new algorithms are needed to implement these new criteria. Also, different criteria tend to produce quite different clusters. In general, hierarchical agglomerative clustering tends to perform worse on higher-dimensional data sets. 3 The last of these, d x belong to a vector space (hence, can be added together and centroid , assumes that the points ~ scaled). If you care about this sort of thing, you probably already noticed my choice of notation ~x which suggests that the ambient space is Rm . This is the most important case, but d centroid makes sense in any normed vector space. (And the rest of this lecture makes sense in an arbitrary metric space.) Lec 21 – pg. 4 of 7 r r r Figure 21.4. Given a positive integer k (in this case k = 3) and a desired radius r, we attempt to create k clusters, each having radius at most r. 21.2 Minimum-Radius Clustering Hierarchical agglomerative clustering focuses on making sure that objects in different clusters are far apart. By contrast, minimum-radius clustering focuses on making sure that objects in the same cluster are close together. Thus, rather than maximizing the minimum distance β(C ), we will now try to minimize a maximum (actually, maximin) distance. Namely, we will try to assign a “center” C.center to each cluster C ∈ C so as to minimize the quantity F(C ) = max min d (~x, C.center) . ~x∈D C ∈C In principle, given k, we need only find the best possible set of k centers—then, each vertex will belong to the cluster centered at whichever of these k points is closest to it. Thus, we can view the quantity F(C ) which we are trying to minimize as depending only on the choice of centers. The problem of minimizing F(C ) in this way is quite difficult. We will make things easier on ourselves by assuming that we are given a goal radius r to strive for (rather than trying to figure out what the minimum possible radius is) (see Figure 21.4). ¡ ¢ Input: Instance data D = {~x1 ,~x2 , . . . ,~xn } with a distance metric d ~x i ,~x j Desired number of clusters, k Desired radius, r Output: Ideally, a set of k points Γ = {~ c 1 , . . . ,~ c k } ⊆ D (the cluster centers) such that each ~x ∈ D is within distance r of some element of Γ. 21.2.1 First strategy: Approximate k The following algorithm gives the correct r but approximates k: assuming a clustering with the desired k and r is possible, the number of clusters is at most k(ln n + 1). Lec 21 – pg. 5 of 7 Algorithm: Imagine each ~x ∈ D as the center of a cluster of radius r. The set © ª S i = points in D within distance r of ~x i consists of all points which should be allowed to belong to this cluster. We now face the following problem: we have a collection of sets S 1 , . . . , S n ⊆ D, and we wish to find a size-k subcollection S i 1 , . . . , S i k such that k [ S i j = D. j =1 This is the set-cover problem, which is NP-hard. However, we can obtain an approximate solution using the approximation algorithm of §18.2: • Greedily choose the set with the largest coverage until all points are covered. • Assuming a set cover exists, this algorithm gives a multiplicative α-approximation, where α = ln |D | + 1. The running time of this algorithm is à O n X ! |S i | , i =1 ¡ ¢ which, depending on the geometric configuration of D, could be anywhere from O(n) to O n2 . 21.2.2 Second strategy: Approximate r Alternatively, we could strictly meet the k requirement but approximate r. Assuming that a clustering with the desired k and r is possible, the following algorithm produces a set of k clusters with radius at most 2r. Algorithm: 1 2 3 4 5 6 7 8 i←0 while D is not empty do i ← i+1 Pick an arbitrary ~x ∈ D Define cluster C i to be all points within distance 2r of ~x C i .center ← ~x D ← D \ Ci return {C 1 , . . . , C i } The costliest line is line 5, which takes O(n) time. As we show below in Lemma 21.2, line 5 is executed at most k times (assuming that a clustering with the desired k and r is possible). Thus, the running time of this algorithm is O(kn). Lemma 21.2. Suppose that a clustering with the desired k and r is possible. Then, after k passes through the loop, D will be empty. © ª Proof. Let C ∗ = C 1∗ , . . . , C ∗k be a clustering with the desired k and r. Let ~x be as in line 4, and let C ∗j ∈ C ∗ be such that ~x ∈ C ∗j . Let ~ c = C ∗j .center. We claim that C i (as defined on line 5) is a superset Lec 21 – pg. 6 of 7 of C ∗j ∩ D. (The intersection with D is necessary because D may have changed during prior iterations of line 7.) To prove this, note that for any ~y ∈ C ∗j ∩ D, we have by the triangle inequality d (~x,~y) ≤ d (~x,~ c) + d (~ c,~y) ≤ r + r = 2r. Thus ~y ∈ C i , which means that C ∗j ∩ D ⊆ C i . Thus, each time line 7 is executed, a new element of C ∗ is removed from consideration. By the time line 7 has been executed k times, all k of the clusters in C ∗ will have disappeared, so that D is empty. Lec 21 – pg. 7 of 7 Lecture 22 Derandomization Supplemental reading in CLRS: None Here are three facts about randomness: • Randomness can speed up algorithms (e.g., sublinear-time approximations) • Under reasonable assumptions, we know that the speed-up can’t be too large. Advances in complexity theory in recent decades have provided evidence for the following conjecture, which is generally believed to be true by complexity theorists:1 Conecture. Suppose there exists a randomized algorithm A which, on an input of size n, runs in time T and outputs the correct answer with probability at least 2/3. Then there exists a deterministic algorithm A 0 which runs in time Poly(n, T) and always outputs the correct answer. • In practice, randomization often doesn’t buy any (or much) speed-up. Moreover, the randomized algorithm is often based on novel ideas and techniques which can equally well be used in a deterministic algorithm. As we said in Lecture 9, a randomized algorithm can be viewed as a sequence of random choices along with a deterministic algorithm to handle the choices. Derandomization amounts to deterministically finding a possible sequence of choices which would cause the randomized algorithm to output the correct answer. In this lecture we will explore two basic methods of derandomization: 1. Conditional expectations. As we walk down the decision tree (see Figure 22.1), we may be able to use probabilistic calculations to guide our step. In this way, we can make choices that steer us toward a preponderance of correct answers at the bottom level. 2. Pairwise independence. Another way to find the correct answer is to simply run the randomized algorithm on all possible sequences of random choices and see which answer is reported most often (see Figure 22.2). In general this is not practical, as the number of sequences is exponential in the number of random choices, but sometimes it is sufficient to check only a relatively small collection of sequences (e.g., those given by a universal hash family). In what follows, we will use the big-cut problem to explore each of these two methods of derandomization: 1 See, for example, Impagliazzo and Wigderson, “P = BPP unless E has sub-exponential circuits: Derandomizing the XOR Lemma” 〈http://www.math.ias.edu/~avi/PUBLICATIONS/MYPAPERS/IW97/proc.pdf〉. head head . . . tail tail head etc. tail . . . ··· Possible outputs, with correct outputs more concentrated in some regions than in others. Figure 22.1. Probabilistic computations can often be used as a heuristic to find the regions in which correct answers are most concentrated. Lec 22 – pg. 2 of 5 input · · · (true randomness)· · · 0000110010 · · · algorithm output input 1110110110 1101000010 0001110111 uniform among N possible strings algorithm 1011101010 0100110101 output Figure 22.2. Sometimes the deterministic component of a randomized algorithm can’t distinguish between a truly random sequence of choices and a well-chosen pseudorandom sequence of choices. Or, sometimes we can simulate a random choice from a large sample set (e.g., the set of all binary strings) using a random choice from a small sample set (e.g., a fixed set of N binary strings). If we can do this with N sufficiently small, it may be feasible to find the answer deterministically by running the randomized algorithm on all N of the possible random inputs. Lec 22 – pg. 3 of 5 Problem 22.1 (Big Cut). Given an undirected graph G = (V, E) with no loops, find a cut (S, V \ S) which crosses at least half of the edges in E. The solution to this problem is approximated by a simple randomized algorithm: Algorithm: 1 2 3 4 5 S←; for each v ∈ V do Flip a coin If heads, then S ← S ∪ {v} return (S, V \ S) The running time is Θ(V ). Proof of approximation. For each edge (u, v) ∈ E, we have · ¸ £ ¤ 1 (u ∈ S and v ∉ S) or Pr (u, v) crosses the cut = Pr = . (u ∉ S and v ∈ S) 2 Thus, if I (u,v) denotes the indicator random variable which equals 1 if (u, v) ∈ S, then we have " # X X ¯ ¯ £ ¤ £ ¤ E # edges crossed by S = E I (u,v) = E I (u,v) = 1 ¯E ¯ . (u,v)∈E (u,v)∈E 2 Chernoff’s bound then tells ¯ ¯ us that, when E is large, it is extremely likely that the number of edges crossed is at least 0.499 ¯E ¯. 22.1 Using Conditional Expectation We can deterministically simulate the randomized approximation to Problem 22.1 by choosing to include a vertex v if and only if doing so (and making all future decisions randomly) would result in a higher expected number of crossed edges than choosing to exclude v. In more detail, say V = {v1 , . . . vn }, and suppose we have already decided whether S should contain each of the vertices v1 , . . . , vk . We pretend that all future decisions will be random (even though it’s not true), so that we can use probabilistic reasoning to guide our choice regarding vk+1 . Ideally, our choice should maximize the value of ¯ £ ¤ E # edges crossed ¯ our existing decisions about v1 , . . . , vk+1 . (22.1) future random choices Consider the partition E = E1 t E2 t E3, where © ª E 1 = (v i , v j ) ∈ E : i, j ≤ k © ª E 2 = (v i , vk+1 ) ∈ E : i ≤ k © ª E 3 = (v i , v j ) ∈ E : j > k + 1 . Whether we choose to include or exclude vk+1 from S, the number of crossed edges in E 1 is unaffected. Moreover, since each edge in E 3 has at least one of its vertices yet to be randomly assigned to S or Lec 22 – pg. 4 of 5 V \ S, the expected number of crossed edges in E 3 is 21 |E 3 | regardless of where we decide to put vk+1 . So in order to maximize (22.1), we should choose to put vk+1 in whichever set (either S or V \ S) produces more crossings with edges in E 2 . (It will always be possible to achieve at least 12 |E 2 | crossings.) We can figure out what the right choice is by simply checking each edge in E 2 . £ ¤ Proceeding in this way, the value of E # edges crossed starts at 12 |E | and is nondecreasing with each choice. Once we make the nth choice, there are no more (fake) random choices left, so we have ¡ 22.2 ¢ £ ¤ # edges crossed = E # edges crossed ≥ 1¯ ¯ 2 E . ¯ ¯ Using Pairwise Independence The sequence of random choices in the randomized approximation to Problem 22.1 is equivalent to picking the set S uniformly at random from the collection S = P (V ) of all subsets of V . In hindsight, the only reason we needed randomness was to ensure that £ ¤ 1 Pr S doesn’t cross (u, v) ≤ S ∈S 2 for each (u, v) ∈ E. (22.2) Recall from §10.1 that a subset of V can be represented as a function V → {0, 1}. In this way, S becomes a hash family H : V → {0, 1}, and (22.2) becomes £ ¤ 1 Pr h(u) = h(v) ≤ h∈H 2 for each (u, v) ∈ E. This will be satisfied as long as H is universal. So rather than taking H to be the collection {0, 1}V of all functions V → {0, 1}, we can take H to be a much smaller universal hash family; there exist universal hash families of size O(V ).2 In this way, we have · E h∈H ¸ # edges crossed by the ≥ cut corresponding to h 1¯ ¯ 2 E . ¯ ¯ In particular, this guarantees that there exists some h ∈ H such that ¶ # edges crossed by the ≥ cut corresponding to h µ 1¯ ¯ 2 E . ¯ ¯ We can simply check every h ∈ H until we find it. Exercise 22.1. What are the running times of the two deterministic algorithms in this lecture? 2 One such universal hash family is described as follows. To simplify notation, assume V = {1, . . . , n}. Choose some prime p ≥ n with p = O(n) and let © ª H = h a,b : a, b ∈ Z p , a 6= 0 , where h a,b (x) = (ax + b) mod p. Lec 22 – pg. 5 of 5 Lecture 23 Computational geometry Supplemental reading in CLRS: Chapter 33 except §33.3 There are many important problems in which the relationships we wish to analyze have geometric structure. For example, computational geometry plays an important role in • • • • • Computer-aided design Computer vision Computer animation Molecular modeling Geographic information systems, to name just a few. 23.1 Intersection Problem Given a set of line segments S 1 , . . . , S n in the plane, does there exist an intersection? Humans, with their sophisticated visual systems, are particularly well-equipped to solve this problem. The problem is much less natural (but still important, given the above applications) for a computer. Input: Line segments S = {S 1 , . . . , S n } in the plane, represented by the coordinates of their endpoints Output: Goal I: Detect. Determine whether there exists an intersection. Goal II: Report. Report all pairs of intersecting segments. ¡ ¢ The obvious algorithm is to check each pair of line segments for intersection. This would take Θ n2 ¡ 2¢ time. For Goal II, this can’t be improved upon, as it could take Θ n time just to report all the intersections (see Figure 23.1). We will focus on Goal I and show that there exists a Θ (n lg n) algorithm. For ease of exposition, we will assume that our line segments are in general position, i.e., • No two endpoints have the same x-coordinate. (In particular, there are no perfectly vertical segments.) • There are no instances of three segments intersecting at a single point. Both of these assumptions are unnecessary and can be dropped after a couple of minor (but careful) modifications to the algorithm. ¡ ¢ ¡ n ¢ n2 ¡ 2¢ Figure 23.1. If the segments S 1 , . . . , S n form a square grid, there could be as many as n 2 2 = 4 = Θ n intersections. In other situations, it could happen that every pair of segments intersects. 23.1.1 Sweep line and preorder Suppose you have a sheet of paper with the line segments S 1 , . . . , S n drawn on it. Imagine dragging a vertical ruler across the page from left to right (see Figure 23.2). This is exactly what we have in mind when we say, For each x ∈ R, let L x be the vertical line at distance x from the y-axis.1 We imagine x increasing over time, so that L x represents a sweep line which moves across the plane from left to right. For each x ∈ R, there is a natural preorder2 on the set of non-vertical line segments in the plane: for two non-vertical segments a, b, we say a ≥ x b if the intersection of a (or its extension, if need be) with L x lies above the intersection of b with L x . This relation needn’t be antisymmetric because it could happen that a and b intersect L x at the same point. In that case, we say that both a ≥ x b and b ≥ x a; but of course, this does not imply that a = b. Of crucial importance is the following obvious fact: Given x ∈ R and line segments a and b, there exists an O(1)-time algorithm to check whether a ≥ x b.3 Also important is the following observation: Observation 23.1. Suppose line segments a and b intersect at the point P = (x∗ , y∗ ). Then one of the following two statements holds: (i) For x > x∗ , we have a ≥ x b. For x < x∗ , we have b ≥ x a. (ii) For x > x∗ , we have b ≥ x a. For x < x∗ , we have a ≥ x b. Conversely, if a and b do not intersect, then one of the following two statements holds: (i) For all x such that L x intersects a and b, we have a ≥ x b. (ii) For all x such that L x intersects a and b, we have b ≥ x a. Another way to view the second part of Observation 23.1 is as follows: 1 It took me a while to figure out how to write down a definition of L that didn’t include the equation “x = x.” x 2 A preorder is a binary relation that is reflexive and transitive, but not necessarily symmetric or antisymmetric. 3 It is not so obvious, however, what the best such algorithm is. Most naïve attempts involve calculating a slope, which is numerically unstable. To remedy this, the best algorithms use cross products to compute orientations without using division. Lec 23 – pg. 2 of 8 L 0.1 L 1.2 L 2.1 b L 3.1 L 3.8 I d c a a b a b a c b d c b c d Figure 23.2. A sweep line traversing from left to right. Notice that b ≥ x d ≥ x c ≥ x a for all x ∈ I, but b ≥3.8 c ≥3.8 d. Thought Experiment 23.2. Consider a finite set of line segments; for concreteness let’s say we are looking at the four segments {a, b, c, d }. There is some interval I (perhaps empty) such that, for each x ∈ I, L x intersects all four segments (see Figure 23.2). For each x ∈ I, the relation ≥ x induces a preorder structure on the set {a, b, c, d }. Observation 23.1 implies that, if there are no intersections among {a, b, c, d }, then this preorder structure does not depend on the choice of x. Thought Experiment 23.2 leads to the following strategy, which is quite ingenious. To introduce the strategy, we first define the preorder data structure. The preorder data structure stores a dynamic collection T of objects along with a preorder on those objects, which we denote by ≥. The supported operations are4 • I NSERT(a) – inserts a into the collection • D ELETE(a) – deletes a from the collection • A BOVE(a) – returns the element immediately above a. If no such element exists, returns N IL. • B ELOW(a) – returns the element immediately below a. If no such element exists, returns N IL. The obvious issue with using this data structure in conjunction with our plane sweep is that it only stores a single preorder relation, which we have denoted above by ≥ (and the usual counterparts ≤, >, <), whereas the relation ≥ x has the potential to change as x varies. However, this needn’t be a problem, as long as: 4 For ease of exposition, we have obscured one detail from view. It is possible for T to contain two elements d, e such that d ≥ e and e ≥ d. In order for the operations A BOVE(d) and B ELOW(d) to work correctly, there must be some rule for breaking ties. To this end, T will also internally store an order º on its elements, such that b º a always implies b ≥ a (but not necessarily conversely). Then, “the element immediately above a” is the element b such that b  a and, whenever c  a, we always have c º b. (Similarly for “the element immediately below a.”) The choice of order º does not matter, as long as it is persistent through insertions and deletions. Lec 23 – pg. 3 of 8 Condition 23.3. During the lifetime5 of any particular element a, the ranking of a with respect to the other elements of T does not ever need to be changed. We are now prepared to give the algorithm: Algorithm: D ETECT-I NTERSECTION(S ) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Sort the endpoints of segments in S by x-coordinate T ← new preorder structure B Run the sweepline from left to right for each endpoint (x, y), by order of x-coordinate do if (x, y) is the left endpoint of a segment a then Insert a into T, using the relation ≥ x to decide where in the preorder a should belong∗ Check a for intersection with T.A BOVE(a) Check a for intersection with T.B ELOW(a) if there was an intersection then return T RUE else if (x, y) is the right endpoint of a segment a then Check T.A BOVE(a) for intersection with T.B ELOW(a) if there was an intersection then return T RUE T.D ELETE(a) return FALSE ∗ This step needs justification. See the proof of correctness. Proof of correctness. First, note the need to justify line 6. In implementation form, T will probably store its entries in an ordered table or tree; then, when T.I NSERT (a) is called, it will take advantage ¡ ¯ ¯¢ of T’s internal order by using a binary search to insert a in Θ lg ¯T ¯ time. In line 6, we are asking T to use the relation ≥ x to decide where to place a. This is fine, provided that ≥ x agrees with T’s internal preorder when considered as a relation on the elements that currently belong to T. (Otherwise, a binary search would not return the correct result.6 ) So, over the course of this proof, we will argue that Whenever line 6 is executed, the relation ≥ x agrees with T’s internal preorder when considered as a relation on the elements that currently belong to T. From Thought Experiment 23.2 we see that the above statement holds as long as the sweep line has not yet passed an intersection point. Thus, it suffices to prove the following claim: Claim. D ETECT-I NTERSECTION terminates before the sweep line passes an intersection point. The claim obviously holds if there are no intersections, so assume there is an intersection. Let’s say the leftmost intersection point is P, where segments a and b intersect. Assume without loss of generality that a’s left endpoint lies to the right of b’s left endpoint (so that b gets inserted into T 5 By “the lifetime of a,” I mean the period of time in which a is an element of T. 6 Actually, in that situation, there would be no correct answer; the very notion of “binary search” on an unsorted list does not make sense. Lec 23 – pg. 4 of 8 b c P a Figure 23.3. If the left endpoint of c lies between a and b to the left of P and the right endpoint of c lies to the right of P, then c must intersect either a or b. before a does). If a ever becomes adjacent to b in T, then either lines 7–8 or line 12 will detect that a and b intersect, and the algorithm will halt. So we are free to assume that a and b are never adjacent in T until after the sweep line has passed P. This means that there exists a line segment c such that the left endpoint of c lies between a and b7 and to the left of P, and the right endpoint of c lies to the right of P. Geometrically, this implies that c must intersect either a or b, and that intersection point must lie to the left of P (see Figure 23.3). But this is impossibe, since we assumed that P was the leftmost intersection point. We conclude that the claim holds. Actually, we have done more than prove the claim. We have shown that, if there exists an intersection, then an intersection will be reported. The converse is obvious: if an intersection is reported, then an intersection exists, since the only time we report an intersection is after directly checking for one between a specific pair of segments. 23.1.2 Running time The running time for this algorithm depends on the implementation of the preorder data structure. CLRS chooses to use a red–black tree8 , which has running time O (lg n) per operation. Thus, D ETECT-I NTERSECTION has total running time Θ (n lg n) . 23.2 Finding the Closest Pair of Points Our next problem is simple: Input: A set Q of points in the plane Output: The two points of Q whose (Euclidean) distance from each other is shortest. ¡ ¢ ¡ ¢ The naïve solution is to proceed by brute force, probing all n2 pairs of points and taking Θ n2 time. In what follows, we will exhibit a subtle divide-and-conquer algorithm which runs in Θ (n lg n) time. 7 That is, if the left endpoint of c is (x, y), then either a ≤ c ≤ b or a ≥ c ≥ b. x x x x 8 For more information about red–black trees, see Chapter 13 of CLRS. Lec 23 – pg. 5 of 8 In §23.1, the pseudocode would not have made sense without a few paragraphs of motivation beforehand. By contrast, in this section we will give the pseudocode first; the proof of correctness will elucidate some of the strange-seeming choices that we make in the algorithm. This is more or less how the algorithm is presented in §33.4 of CLRS. The algorithm begins by pre-sorting the points in Q according to their x- and y-coordinates: Algorithm: C LOSEST-PAIR(Q) X ← the points of Q, sorted by x-coordinate 2 Y ← the points of Q, sorted by y-coordinate 3 return C LOSEST-PAIR -H ELPER(X , Y ) 1 Most of the work is done by the helper function C LOSEST-PAIR -H ELPER, which makes recursive calls to itself: Algorithm: C LOSEST-PAIR -H ELPER(X , Y ) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 if | X | ≤ 3 then Solve the problem by brute force and return x∗ ← the median x-coordinate of X Let X L ⊆ X consist of those points with x-coordinate ≤ x∗ Let X L ⊆ X consist of those points with x-coordinate > x∗ Let YL ⊆ Y consist of those points which are in X L Let YR ⊆ Y consist of those points which are in X R B Find the closest two points in the left half 〈 p L , q L 〉 ← C LOSEST-PAIR -H ELPER(X L , YL ) B Find the closest two points in the right half 〈 p R , q R 〉 ← C LOSEST-PAIR -H ELPER(X R , YR ) δR ← distance from p L to q L δR ← distance from p R to q R δ ← min {δL , δR } Y 0 ← those points in Y whose x-coordinate is within δ of x∗ B Recall that Y 0 is already sorted by y-coordinate ²←∞ ¯ ¯ for i ← 1 to ¯Y 0 ¯ − 1 do p ← Y 0 [i] ¯ ¯ª © for q in Y 0 [i + 1, . . . , min i + 7, ¯Y 0 ¯ ] do if ² > distance from p to q then ² ← distance from p to q p∗ ← p q∗ ← q if ² < δ then return 〈 p∗ , q∗ 〉 else if δR < δL then return 〈 p R , q R 〉 else return 〈 p L , q L 〉 Lec 23 – pg. 6 of 8 Right half Left half Figure 23.4. C LOSEST-PAIR -H ELPER divides the set X into a left half and a right half, and recurses on each half. 23.2.1 Running time Within the procedure C LOSEST-PAIR -H ELPER, everything except the recursive calls runs in linear time. Thus the running time of C LOSEST-PAIR -H ELPER satisfies the recurrence T(n) = 2 · T ¡n¢ 2 + O(n), (23.1) which (by the Master Theorem in §4.5 of CLRS) has solution T(n) = O (n lg n) . (23.2) Note that, if we had decided to sort within each recursive call to C LOSEST-PAIR -H ELPER, the O(n) term in (23.1) would have instead been an O (n lg n) term and the solution would have been ¡ ¢ T(n) = O n (lg n)2 . This is the reason for creating a helper procedure to handle the recursive calls: it is important that the lists X and Y be pre-sorted so that recursive calls need only linear-time operations. Note also that, if instead of lines 20–24 we had simply checked the distance between each pair of ¡ ¢ points in Y 0 , the O(n) term in (23.1) would have instead been an O n2 term, and the solution would ¡ ¢ have been T(n) = O n2 . 23.2.2 Correctness C LOSEST-PAIR -H ELPER begins by recursively calling itself to find the closest pairs of points on the left and right halves. Thus, lines 15–24 are ostensibly an attempt to check whether there exists a pair of points 〈 p∗ , q∗ 〉, with one point on the left half and one point on the right half, whose distance is less than that of any two points lying on the same half. What remains to be proved is that lines 15–24 do actually achieve this objective. By the time we reach line 15, the variable δ stores the shortest distance between any two points that lie on the same half. It is easy to see that there will be no problems if δ truly is the shortest possible distance. The case we need to worry about is that in which the closest pair of points—call it 〈 p, q〉—has distance less than δ. In such a case, the x-coordinates of p and q would both have to be within δ of x∗ ; so it suffices to consider only points within a vertical strip V of width 2δ centered at x = x∗ (see Figure 23.5). These are precisely the points stored in the array Y 0 on line 15. Say p = Y 0 [i] and q = Y 0 [ j], and assume without loss of generality that i < j. In light of lines 20–24 we see that, in order to complete the proof, we need only show that j − i ≤ 7. Lec 23 – pg. 7 of 8 V δ δ Figure 23.5. It suffices to consider a vertical strip V of width 2δ centered at x = x∗ . δ SL SR Figure 23.6. Each of S L and S R can hold at most 4 points. (Actually, in completely ideal geometry, S R cannot contain 4 points because its left boundary is excluded. But since the coordinates in a computational geometry problem are typically given as floating point numbers, we are not always guaranteed correct handling of edge cases.) Say p = (p x , p y ). Let S L be the square (including boundaries) of side length δ whose right side lies along the vertical line x = x∗ and whose bottom side lies along the horizontal line y = p y . Let S R be the square (excluding the left boundary) of side length δ whose left side lies along the vertical line x = x∗ and whose bottom side lies along the horizontal line y = p y . It is evident that q must lie within either S L or S R (see Figure 23.6). Moreover, any two points in the region S L are separated by a distance of at least δ; the same is true for any two points in the region S R . Thus, by a geometric argument9 , S L and S R each contain at most four points of Y 0 . In total, then, S L ∪ S R contains at most eight points of Y 0 . Two of these at-most-eight points are p and q. Moreover, since Y 0 consists of all points in V sorted by y-coordinate, it follows that the at-most-eight points of S L ∪ S R occur consecutively in Y 0 . Thus, j − i ≤ 7. 9 One such argument is as follows. Divide S into northeast, northwest, southeast and southwest quadrants. Each L quadrant contains at most one point of Y 0 (why?), so S L contains at most four points of Y 0 . Lec 23 – pg. 8 of 8 Index α-approximation, 18:1 α-competitive, 12:2 ²-universal, 10:3 accounting method, 11:2 Ackermann function, 16:7 acyclic, 3:2 additive approximation, 18:2 aggregate analysis, 11:1 amortized analysis, 10:7 amortized cost, 10:7, see also amortized analysis approximation, see α-approximation approximation algorithm, 18:1 associative array, 10:1 augmented flow, 13:3 augmenting path, 13:3 big-cut problem, 22:1 binary code, 19:2 binomial distribution, 9:2 bipartite graph, 14:4 capacity of a cut in a flow network, 13:7 of an edge in a flow network, 13:1 certificate, 17:3 clustering, 21:1 codeword, 19:2 collision, 10:2 competitive, see α-competitive competitive analysis, 12:3 complete graph, 3:2 complexity class, 17:2 compression, 19:1 connected components, 3:6 connected graph, 3:2 convolution, 5:5 Cook reduction, 17:4 correct algorithm, 1:2 cut of a flow network, 13:7 of a graph, 3:5 cycle, 3:2 decision problem, 17:2 depth of a node in a rooted tree, 4:3 directed graph, 3:2 disjoint-set data structure, 4:2, 16:1 divide-and-conquer, 5:1 dynamic programming, 2:6, 6:1 Edmonds–Karp algorithm, 13:6 “efficient” algorithm, 17:3 flow, 13:1 across a cut in a flow network, see net flow flow network, 13:1 flow value, 13:2 Ford–Fulkerson algorithm, 13:2 FPTAS, see fully polynomial-time approximation scheme fully polynomial-time approximation scheme, 18:2 Gaussian distribution, 9:1 God’s algorithm, 12:1 graph, 3:2 greedy algorithm, 3:1 hash family, 10:2 hash function, 10:1 hash table, 10:1 integer flow, 13:6 inverse Ackermann function, 16:7 isolated vertex, 4:4 Karp reduction, 17:4 keyed array, see associative array Kruskal’s MST algorithm, 3:6 Las Vegas algorithm, 8:2 linear programming, 7:4 load on an index in a hash table, 10:2 loop invariant, 1:4 lossless compression, 19:1 lossy compression, 19:1 maximum flow problem, 13:1 min-priority queue, 4:5 minimum spanning tree, 3:4 Monte Carlo algorithm, 8:2 move-to-front, 12:3 MST, see minimum spanning tree MTF, see move-to-front multi-commodity flow, 13:9 multiplicative approximation, 18:1 net flow, 13:7 nondeterministic algorithm, 17:3 normal distribution, 9:1 NP, 17:3 NP-complete, 17:5 NP-hard, 17:4 offline algorithm, 12:1 online algorithm, 12:1 optimal substructure, 6:1 oracle, 17:4 P, 17:3 P vs. NP, 17:2 path, 3:2 perfect hash table, 10:6 Poisson distribution, 10:4 polynomial-time, 17:3 in multiple variables, 20:1 polynomial-time approximation scheme, 18:2 polynomial-time reduction, 17:4 potential method, 11:3 prefix coding, 19:3 preorder, 23:2, see also preorder data structure preorder data structure, 23:3 Prim’s MST algorithm, 4:4 pseudo-polynomial–time, 18:2 PTAS, see polynomial-time approximation scheme rank, 1:3 reachable, 3:6 reduction, see polynomial-time reduction residual capacity, 13:3 residual network, 13:3 restriction of a graph to a set of vertices, 3:6 rooted tree, 4:3 safe choice, 4:4 search problem, 17:2 simple cycle, 3:3 simple path, 3:3 spanning tree, 3:3 stack, 11:1 sweep line, 23:2 table doubling, 10:6 tree, 3:2 Turing machine, 17:2 undirected graph, 3:2 union bound, 9:6 uniquely readable, 19:3 universal hash family, 10:2 van Emde Boas data structure, 15:1 vEB, see van Emde Boas data structure verifier, 17:3 weighted undirected graph, 3:2 witness, 17:3