CHAPTER  3

ALGORITHM DESIGN

3.0 Recursion
-------------

There is really only one general method of algorithm design - break the problem
down into smaller subproblems. Recursion is just a version of this break-down
technique which breaks a problem into smaller subproblems of the SAME KIND. It
does not mean that the resulting program need have any recursive procedures in
it.

For example, let's use recursion to analyse the sorting problem. The idea is to
reduce the problem of sorting n items to the problem of sorting n-1 items.
There are two approaches. One is to sort the first n-1 items then insert the
n'th item in its correct place, giving the familiar insertion sort:

sort a[1]...a[n]
----------------
if n=1 return
sort a[1]...a[n-1]
insert a[n] in the sorted list a[1]...a[n-1]

Alternatively, we could search for the largest item, put it last, then sort the
first n-1 items, giving a selection sort:

Sort a[1]...a[n]
----------------
if n=1 return
find the largest item and swap it with a[n].
sort a[1]...a[n-1]

Both these algorithms take time f(n) satisfying

f(n) = f(n-1) + O(n)      =>      f(n) = O(n^2)

Insertion sort is simpler and faster than any of the other famous O(n^2)
algorithms (especially bubble sort). It is fast on small problems, and fast
(linear, i.e. O(n)) on arrays which are almost sorted already. This makes it
suitable for combining with more sophisticated algorithms.

Remember, a recursive design should always begin with a simple case which can
be solved easily.

3.1 Divide And Conquer
----------------------

This is another important technique in algorithm design. It is just a more
sophisticated use of recursion in which we try to break the problem into equal
sized pieces. In other words, there are two features: a breakdown into
(1) SIMILAR PROBLEMS of (2) SIMILAR SIZES.

For example, consider searching for an item in a sorted array a[1]...a[n]. A linear search, looking at each element in turn is O(n). A binary chop is better. The idea is to break the search into two equal sized searches: To search a[i]...[j] for value v -------------------------------- if i > j return "not found" look at the middle element k = (i + j) div 2 if v = a[k] return "k" if v < a[k] then search a[i]...a[k-1] if v > a[k] then search a[k+1]...a[j] WARNING - if you try to "optimize" the 3-way test into 2-way test, you are liable to produce an incorrect algorithm. The time taken is: f(n) = f(n/2) + O(1) => f(n) = O(log(n)) Now let's try to break a sorting problem into two equal sized sorting problems. One approach is to sort the first half, sort the second half and combine the results, giving a merge sort: Sort items a[i]...a[j] ---------------------- if there is only one item, return sort the first half of the items a[1]...a[k] where k = (i+j) div 2 sort the second half of the items a[k+1]...a[j] merge the two sorted sequences Assuming that merging of sequences of length n can be done in time O(n), (see questions) the time f(n) to sort n items using this algorithm is: f(n) = 2f(n/2) + O(n) => f(n) = O(n*log(n)) Alternatively, we could try to split the numbers into low numbers and high numbers first, then sort the low numbers and sort the high numbers. We get the quicksort algorithm: Sort items a[i]...a[j] ---------------------- 1. if there is only one item, return 2. choose a "medium" value v 3. arrange for a[i]...a[k] to contain values <= v and for a[k+1]...a[n] to be > v (whatever k happens to be) 4. sort a[i]...a[k] 5. sort a[k+1]...a[n] The questions remaining in this algorithm are how to choose v at step (2), and how to do the splitting at step 3. Step 3 is done by another recursive design:
Split a[i]...a[j] into those <= v and those > v ----------------------------------------------- if i > j then return if a[i] <= v then split a[i+1]...a[j] and return if a[j] > v then split a[i]...a[j-1] and return now we know a[i] and a[j] are mis-placed swap a[i] and a[j] split a[i+1]...a[j-1] This splitting can be implemented with simple loops rather than a recursive procedure. (Increment i until a[i]>v, decrement j until a[j]<=v, swap a[i] and a[j], repeat). The splitting algorithm clearly examines each item once, so is O(n). More formally, the time taken is f(n) = f(n-1) + O(1) => f(n) = O(n) As for choosing v, we can choose a random item from a[i]...a[j], or if we assume the items are in random order to start with, we can just use v = a[i]. A popular compromise is to use the median of a[i], a[j] and a[(i+j) div 2]. On average, half the items will be less than v and half greater, so the time for quicksort ON AVERAGE is: f(n) = 2f(n/2) + O(n) => f(n) = O(n*log(n)) Thus Divide & Conquer design gives us two O(n*log(n)) algorithms for sorting. Which sorting algorithm should be used in practice? Mergesort is a useful tool in sorting large files. In memory, insertion sort is best for small problems (say up to 10 items) or when speed is unimportant. Otherwise use quicksort. Either way, take the code from (e.g.) Sedgewick. There are various ways of optimising quicksort. One problem is that it spends a lot of its time sorting small pieces of the array, at which insertion sort is better. Thus a good scheme is to change the first line to "if the number of items is less than 10 return". After the quicksort, insertion sort can be used once on the whole array to finish the job - this will be quick since the array is nearly sorted already. Also, one of the two recursive calls in the quicksort algorithm can easily be replaced by a loop. 3.2 Divide & Conquer - Multiplying Polynomials ---------------------------------------------- As another example, consider multiplying polynomials, say with n coefficients. p = 1 + 2*x + 3*x^2 + 4*x^3 + 5*x^4 + 6*x^5 q = 6 + 5*x + 4*x^2 + 3*x^3 + 2*x^4 + 1*x^5 The obvious algorithm involves multiplying each term of one poly by each term of the other, giving an O(n^2) algorithm. The divide and conquer approach suggests splitting each polynomial in two and reducing the problem to that of multiplying polynomials of length n/2:
p = 1+2*x+3*x^2 + x^3*(4+5*x+6*x^2) ----------- ----------- a b q = 6+5*x+4*x^2 + x^3*(3+2*x+1*x^2) ----------- ----------- c d The product p*q is: a*c + x^3*(a*d + b*c) + x^6*(b*d) If we have to perform all 4 multiplications a*c, a*d, b*c, b*d then we have gained nothing. However, a*d + b*c = (a + b) * (c + d) - a*c - b*d and so the product p*q can be calculated using only three half-length multiplications, namely a*c, b*d and (a+b)*(c+d). The additions etc. take O(n) time as they can be done coefficient by coefficient, and so the time taken is: f(n) = 3*f(n/2) + O(n) => f(n) = O(n^log(3)) = O(n^1.59) This is better than O(n^2), but only for large n because of the overheads. There is a O(n*log(n)*log(log(n))) algorithm for multiplying large polynomials based on the Fast Fourier Transform, but it is too complex to be practical except for VERY large problems. For matrix calculations it has been shown that the costs of multiplying two matrices, inverting a matrix, finding the determinant of a matrix etc. are all the same - a fast algorithm for one results in fast algorithms for the others. The normal multiplication method for NXN matrices takes O(N^3) time. There is a divide & conquer method called Strassen's algorithm for multiplying large NXN matrices in O(N^log7) = O(N^2.81...) based on the fact that 2 by 2 matrices can be multiplied using only seven multiplications instead of the usual 8. It becomes noticeably better than the standard method at around N=100. There is a better algorithm still (about O(N^2.795)), but the best possible power of N is not known. Note that the problem size is, strictly speaking, n = N^2 and so the standard method is O(n^(3/2)) = O(n^1.5) and Strassen's is O(n^log7/2) = O(n^1.40...). 3.3 Divide & Conquer - Graphics ------------------------------- As another example, consider Warnock's algorithm for producing a two-dimensional picture from 3-dimensional information (the hidden surface problem). The aim is (for example) to produce a picture of a cube:
----------------------------------- | | | ___________ | | /__________/| | | | || | | | || | (six squares) | | || | | |__________|/ | | | | | | | ----------------------------------- given a description of the position of the cube, and the viewing position and direction. Suppose the cube is represented as six squares of different colours, Front, Back, Top, Bottom, Left, and Right. The viewing rectangle is split into 4 equal pieces, and the squares are also cut into pieces. ------------------|----------------- | (4 pieces) | (4 pieces) | | _______|____ | | /_______|___/| | | | | || | | | | || | ------------------|----------------- | | | || | | |_______|___|/ | | | | | (1 piece) | (4 pieces) | | | | ------------------|----------------- Any piece which is completely covered by another piece can be discarded (see lower left corner). As this process is repeated, a rectangle will eventually have only one piece (or none) in it, and so can be easily be "coloured in" on the view screen. If a rectangle becomes very small without becoming simple enough to colour in a single colour, it can be coloured by some simple heuristic. Algorithms of this type are very common and useful in graphics, particularly as a pre-processing stage, or for use on multiple processors. 3.4 Removing Recursion ---------------------- Recursion can always be removed from an algorithm if desired. This removes the cost of the procedure calls involved. Advice - always remove recursion if you can do it without using your own stack. If you have to use your own stack, you might make the program slower, since the compiler may have access to a special fast hardware stack. Removing recursion only changes the speed of an algorithm by a constant factor, so it does not change the analysis. The easiest kind of recursion to remove is called HEAD- or TAIL-RECURSION, where there is a recursive call that appears very near the beginning or end of the procedure. For example the simplest cases are:
proc (n) proc (n) -------- -------- if n=0 return becomes for i := 1 to n do proc (n-1) etc..... etc..... proc (n) proc (n) -------- -------- if n=0 return becomes for i := n downto 1 do etc..... etc..... proc (n-1) These assume that n is not altered during the procedure. As another example, the last call in quicksort can be removed as follows: Sort items a[i]...a[j] ---------------------- 1. while i+1 < j (number of items > 1) 2. arrange for a[i]...a[k] to be less than a[k+1]...a[j] 3. sort a[i]...a[k] 4. i := k+1 (instead of: sort a[k+1]...a[j]) More complicated recursions can always be removed in the same way that a compiler does, i.e. by using a stack of local variables. This can often be thought of as "things left to be done". For example Quicksort becomes (using a stack of ranges): Sort a stack of ranges ---------------------- while the stack is not empty do take a range i,j from the stack while i+1 < j arrange for a[i]...[k] < a[k+1]...[j] add i,k to the stack i := k+1 3.5 Dynamic programming ----------------------- This is a fancy name for a simple idea, again built on recursion. The idea is to speed up algorithms at the expense of space by keeping a table of previous smaller results. For example, in calculating factorial(n), it helps to keep an array of results, up to the largest previously requested, say r. When factorial(n) is requested, if n <= r, then the value in the array is returned. If n > r, then entries r+1 to n are filled in, using one multiplication each. This speeds the algorithm up from O(n) to an average of O(1). Another example is pagination. Suppose there is a paragraph of text which consists of n words which must be split into lines in the best possible way. Suppose we have an array w[1]...w[n] of words. We are given a function LineCost(i,j) which measures the cost of putting words w[i]...w[j] on one line. The cost is high if the result would be aesthetically un-pleasant, and (effectively) infinite if it would be impossible. A straightforward recursive design gives:
To find min cost of paragraphing words w[1]...w[n] -------------------------------------------------- for i := 1 to n-1 do find min cost of paragraphing w[1]...w[i] add LineCost(i+1,n) return the minimum of the n-1 costs found In other words, for each i we find the best way of splitting the first i words, and putting the remaining words on one more line. If it were implemented directly as a recursive procedure, the complexity would be roughly f(n) = f(1)+f(2)+f(3)+...+f(n-1) = O(1+1+2+4+8+16+...+2^n) = O(2^n) If instead we keep a table of smaller results, so that ParaCost[i] becomes a table lookup instead of a recursive call, we get: To find min cost of paragraphing words w[1]...w[n] -------------------------------------------------- for j := 1 to n do for i := 1 to j-1 do find Paracost[i] + Linecost(i+1,j) Paracost[j] := minimum of the costs found return Paracost[n] Now f(n) = O(n^3), assuming that LineCost takes time O(n). This is the approach which the TeX typesetting package takes. 3.6 Dynamic Programming - Knapsack ---------------------------------- Another example is the knapsack problem. Suppose there are items of n types varying in size and value, and a knapsack of size K, e.g. size each (say cc) value each (say \$) stack of gold bars 12 24 heap of silver coins 7 13 tree of banknotes 2 1 knapsack 28 What is the greatest value of items that can be put in the knapsack? Any number of items of one type can be used. An approximate answer can be obtained using a "greedy" algorithm by starting with the items which are most valuable for their size. In the above example, gold is best (2 \$/cc) then silver (13/7 \$/cc) then notes (1/2 \$/cc). To fill the knapsack, we put in as much gold as possible (2 bars) then silver coins (none) then notes (2) for a total value of 48+2 = 50. Unfortunately, this algorithm does not give the best answer, as it is better to put in 4 silver coins for a total value of 4*13 = 52. For an exact answer, we must turn to recursion. There are n choices of first item, and for each choice we find the best way of filling the remaining space. Then we take the maximum over the n choices.
To find the max value of a knapsack of size K --------------------------------------------- for i := 1 to n do find value[i] plus max value of a knapsack of size (K-size[i]) return the maximum of the values found The time taken f(n) is again at least a sum of n previous values f(i), giving f(n) = O(2^n). Again if we tabulate the answers as knapsack[j], we get: To find max value of knapsack of size K --------------------------------------- for j := 1 to K do for i := 1 to n do if size[i] <= K then find value[i] + knapsack[j-size[i]] knapsack[j] := maximum value found return knapsack[K] which takes O(K*n) time. If we keep track of the best item as well as the max value at each stage, we can construct the optimum contents as well as the max value of the contents. The method can be adapted to the case when there are a limited number of items of each type. Note that this method is proportional to K, and is no good if K is large, say exponential in n. If the numbers in the problem are real instead of integer, this is equivalent to K being large. 3.7 Dynamic Programming - Sequences of Matrices ----------------------------------------------- Another example comes from multiplying matrices of different sizes. Suppose Mab stands for a matrix with 'a' rows and 'b' columns. The cost of multiplying Mab by Mbc is a*b*c (since there are 'a' rows of Mab to multiply by 'c' columns of Mbc, and each row/col has 'b' items in it). The cost of multiplying many matrices depends on the order. If a=2, b=3, c=1, d=4 then (Mab * Mbc) * Mcd costs a*b*c + a*c*d = 14 Mab * (Mbc * Mcd) costs b*c*d + a*b*d = 36 so the first order is better. Suppose we know the size of each of n matrices M1...Mn, and we want to find the minimum cost of multiplying them all together, i.e. of calculating M1*M2*...*Mn, we have a recursive design: To find min cost of multiplying M1*M2*...*Mn -------------------------------------------- for i := 1 to n-1 do find min cost of multiplying M1*...*Mi find min cost of multiplying Mi+1*...*Mn add these to the cost of multiplying (M1*...*Mi) by (Mi+1*...*Mn) return the minimum of these n-1 costs A two-dimensional table cost[i,j] is needed to avoid an exponential algorithm. Cost[i,j] is the cost of multiplying Mi*...*Mj. When calculating cost[i,j], we need to have available previous values of cost[i,k] and cost[k+1,j].
We can ensure this by filling in cost[i,j] in order of increasing j-i. First we fill in the diagonal i=j. Second we fill in the diagonal j-i=1. Third, we fill in the diagonal j-i=2, and so on. j --> ------------ i | 1 2 3 4 5 | 1 2 3 4 | | 1 2 3 v | 1 2 | 1 To find min cost of multiplying M1*...*Mn ----------------------------------------- initialise cost[i,i] to zero for L := 1 to n-1 do for i := 1 to n-L do calculate cost[i,i+L] as follows for j := i to i+L-1 do calculate cost[i,j] + cost[j+1,i+L] add cost of multiplying (Mi*...*Mj) by (Mj+1*...*Mi+L) cost[i,i+L] := min of values found return cost[1,n] This gives us an O(n^3) algorithm. Again, if we record with cost[i,j] the best place to split Mi...Mj, we can reconstruct an optimum order of multiplications. 3.8 Dynamic Programming - Optimal Search Trees ---------------------------------------------- Suppose we have a set of strings s[1]...s[n] in order, with known frequencies of occurrence f[1]...f[n]: s: A B C D E F G f: .20 .10 .05 .15 .30 .15 .05 We know that there will be a large number of searches for these strings, the strings appearing with the frequencies given (e.g. 20% of the searches are for string A). We want to create a binary search tree from them such as D --- 1 comparison / \ / \ A F --- 2 comparisons \ / \ B E G --- 3 comparisons \ C --- 4 comparisons before searching begins. We want the tree with minimum average search time. We measure the cost of a tree by multiplying the number of comparisons needed to find a string (distance from the root + 1) by the frequency for the string and summing. We want the tree with smallest cost. (This is like Huffman trees, except we must preserve the order). All nodes to the left of a node must be earlier, and all nodes to the right must be later.
If we choose s[i] as root, its left subtree must contain s[1]...s[i-1] and its right subtree must contain s[i+1]...s[n]. Thus we try each string s[1]...s[n] as root. Let weight[i,j] be the sum of the frequencies f[i]+...+f[j]. A recursive design gives: Find optimum time for s[i]...s[j] --------------------------------- for k := i to j do find optimum time t1 for s[i]...s[k-1] find optimum time t2 for s[k+1]...s[j] calculate t1 + t2 + weight[i,j] return the minimum of the n times found Note that every string in the tree takes one more comparison than it does in the two subtrees, so the total time is t1+(f(i)+..+f(k-1)) for the left tree, plus f(k) for the root plus t2+(f(k+1)+..+f(j)) for the right tree, giving a total of t1+t2+weight[i,j]. Again, this gives an exponential algorithm, and we need a two-dimensional table cost[i,j] = "optimum time for s[i]...s[j]". Again, we must fill it in order of increasing length j-i. Find optimum time for s[1]...s[n] --------------------------------- for L := 1 to n-1 do for i := 1 to n-L do find cost[i,i+L] as follows for j := i to i+L do t1 := cost[i,j-1] t2 := cost[j+1,i+L] find f[j] + t1 + weight[i,j-1] + t2 + weight[j+1,i+L] cost[i,i+L] := minimum of the costs found return cost[1,n] Again, this gives an O(n^3) algorithm. 3.9 Backtracking ---------------- Some problems are intractable, that is, there is no known efficient algorithm, and we have to resort to an exhaustive search for possible solutions. Such algorithms are generally exponential, and suffer from "combinatorial explosion". That is, there is a problem size, say 20 items, beyond which there is no practical hope of solution. To cope with 21 items would require a computer 20 times as powerful, and to cope with 22 items would require 400 times the power. The aim in designing exhaustive search algorithms is to put off the explosion as long as possible. For example, if we can improve an O(2^n) algorithm to O(2^(n/2)), we may put off the explosion from size 20 to size 40. The most famous problem of this kind is the Travelling Salesperson Problem (which always used to be called the Travelling Salesman Problem and is now sometimes called the Weary Traveller Problem). Given a number of towns and inter-town distances, find the shortest route which visits each town once. For example:
A B C D E .... A 5 6 7 8 B 4 3 3 C 2 4 D 3 E ... The statement of the problem says that the traveller must visit each town only once. There are circumstances where you might think it better to visit a town more than once. Suppose town T is visited twice, the second time being between visits to A and B (...T...ATB...). Then this must be because the sum of distances AT + TB is smaller than distance AB, otherwise it would be better to leave out the second visit. But that means that the shortest distance from A to B (which just HAPPENS to go through T) hasn't been recorded properly. Thus we may assume that the triangle inequality holds (direct routes are always shortest) and that the traveller visits towns only once. The most obvious algorithm for this is to generate all the n! permutations of the towns, calculate the length of each tour, and keep track of the best one so far. This is O(n!), but remember that success is measured by cutoff point and not by asymptotic behaviour. There are obviously gains to made through symmetry, but symmetry is a rather special problem. A more general technique is PRUNING. Suppose we generate the permutations as follows: begin with each town A,B,C,D... in turn, and for each generate all permutations of the n-1 remaining items giving ABCD... ABDC... ACBD... ACDB... ... When all possibilities of one type have been exhausted, we "backtrack" and try the next type. Note this is just another form of recursion. We may as well keep track of the length of the partial path. Suppose the best tour found so far is 20 and the length of ACBD is >= 20. Then we can ignore all permutations beginning ACBD... This is called pruning because it is cutting off a whole subtree of the search tree: . / \ A / \ AC / ACB / ACBD (this subtree is pruned) //\\
Pruning is more successful if an answer close to optimal is found early. Thus it pays to make the order of search one which makes this likely. For example, having decided on ACBD, the next town chosen could be the one closest to D. Note that the majority of nodes in a tree are near the bottom, so pruning can be very effective. It may be worth applying a good heuristic algorithm first to get a good solution to aid the pruning. Using very sophisticated techniques, the explosion point for this problem is around 100-200. 3.10 Heuristic Methods For Intractable Problems ----------------------------------------------- If faced with intractable problems as above, with problem sizes over the explosion point, we have to resort to approximate methods. These are difficult to design and almost impossible to analyse. Perhaps the commonest method is to look for a LOCAL OPTIMUM. In the travelling salesperson case, we could begin with a reasonable solution and make small changes. When no more small changes can improve the solution, we have a local maximum - no solution "close by" is any better. One quite successful method is to split the current solution into (say) 3 paths and investigate all possible ways of reconnecting the paths: ABCDEFGHI JKLMNOPQR STUVWXYZ ABCDEFGHI JKLMNOPQR ZYXWVUTS ABCDEFGHI RQPONMLKJ STUVWXYZ ABCDEFGHI RQPONMLKJ ZYXWVUTS ABCDEFGHI STUVWXYZ JKLMNOPQR ABCDEFGHI STUVWXYZ RQPONMLKJ ABCDEFGHI ZYXWVUTS JKLMNOPQR ABCDEFGHI ZYXWVUTS RQPONMLKJ It is difficult to say how well the method performs as there is no practical method of finding the best solution to compare against.