Combinations
A great many problems in math and computer science have recursive solutions. In this section of the lecture we will study a problem from mathematics that has an elegant recursive solution.
The combination symbol 
or "n choose k" has many applications in mathematics. In probability theory
counts the number of ways that we can get k heads in a sequence of n flips of a fair coin.
is computed via the formula
Two special cases are easy to compute.
One way to compute a combination is to compute all the factorials involved and then divide the denominator into the numerator. This is not practical for even moderately large value of n because the factorials grow very quickly as a function of n. Instead, the most common method for computing a combination makes use of this recursive relationship:
This recursive relationship coupled with the two special cases makes it possible to construct a recursive function definition in Python to compute a combination.
def C(n,k):
if k == 1:
return n
if k == n:
return 1
return C(n-1,k-1)+C(n-1,k)
This works, and does compute the correct value for a combination.
This approach does have one flaw. Consider what happens when we try to compute C(8,5).
C(8,5) calls C(7,4) and C(7,5) C(7,4) calls C(6,3) and C(6,4) C(7,5) calls C(6,4) and C(6,5) C(6,3) calls C(5,2) and C(5,3) C(6,4) calls C(5,3) and C(5,4) C(6,4) calls C(5,3) and C(5,4) C(6,5) calls C(5,4) and C(5,5)
What you start to notice after a while is that the same functions are getting called multiple times. This effect becomes so severe for large n and k that the simple recursive function shown above becomes immensely inefficient.
We can construct a program that illustrates just how bad this problem is. The trick is to declare a global two dimensional array
counts = np.zeros((18,18),dtype=int)
that stores information about how many times we call C with each pair of parameter values n and k.
We then modify the function that computes the combination to record how many times it gets called with each pair of parameters.
def C(n,k):
counts[n][k] += 1;
if k == 1:
return n
if k == n:
return 1
return C(n-1,k-1)+C(n-1,k)
At the end of the program we dump this count data to a file.
def saveCounts():
with open('counts.txt','w') as f:
for row in counts:
for col in row:
f.write('{:5d}'.format(col))
f.write('\n')
The results are shocking. Here is the set of counts that results when we try to compute C(16,10):
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 1287 1287 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 495 1287 792 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 165 495 792 462 0 0 0 0 0 0 0 0 0 0 0 0 0
0 45 165 330 462 252 0 0 0 0 0 0 0 0 0 0 0 0
0 9 45 120 210 252 126 0 0 0 0 0 0 0 0 0 0 0
0 1 9 36 84 126 126 56 0 0 0 0 0 0 0 0 0 0
0 0 1 8 28 56 70 56 21 0 0 0 0 0 0 0 0 0
0 0 0 1 7 21 35 35 21 6 0 0 0 0 0 0 0 0
0 0 0 0 1 6 15 20 15 6 1 0 0 0 0 0 0 0
0 0 0 0 0 1 5 10 10 5 1 0 0 0 0 0 0 0
0 0 0 0 0 0 1 4 6 4 1 0 0 0 0 0 0 0
0 0 0 0 0 0 0 1 3 3 1 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
What we can read off from this table is that C(2,1), C(2,2), and C(3,2) each got called a total of 1287 times in the process of computing C(16,10). This problem of redundant function calls gets exponentially worse as n and k get larger. By the time you get to n = 30 this simple recursive method for computing C(n,k) becomes completely impractical.
Remembering results
The simple fix needed to eliminate redundant calls to a recursive function is called memoization. The idea is to remember any values we compute, so that the next time we are asked to compute a value we can look it up instead of recomputing it. This works best for recursive functions with a moderate number of integer parameters. In that case, we can remember values we have computed earlier by storing the values in a table. Whenever we are asked to compute a value, we start by consulting the table.
In the case of the function to compute combinations, we proceed as follows. We start by making a global two dimensional array to hold the remembered function values.
c = np.zeros((100,100),dtype=int)
We then rewrite the recursive function to use the table. Whenever we are asked to compute C(n,k) for a particular n and k, we start by checking the table to see if has already been computed. If it hasn't already been computed, we compute it and save the value in the table before returning with the result.
def C(n,k):
result = c[n][k]
if result == 0:
if k == 1:
result = n
elif k == n:
result = 1
else:
result = C(n-1,k-1)+C(n-1,k)
c[n][k] = result
return result
Replacing recursion with a table
An interesting side effect of the memoized solution is that it ends up filling up a table c[n][k] with values for C(n,k).
A more direct approach to doing the combinations problem efficiently is to just write code that fills the table with the right values.
We can do this with a function
def initTable(): # Fill the table with correct values # for c[n][k] for all values of n and k
Once the c array is properly filled with values, we can rewrite the function to compute combinations
def C(n,k):
if k <= n:
return c[n][k]
return 0
To write the initTable function we simply have to replicate the logic in the original recursive solution
def C(n,k):
if k == 1:
return n
if k == n:
return 1
return C(n-1,k-1)+C(n-1,k)
in a set of loops. There are two things we need to do. The first is to write a pair of loops that fill in the c[n][k] entries for all n and k in the base cases. The base cases for combinations are k = 1 and k = n:
def initTable():
# Fill in the base cases
for n in range(0,MAX_N):
c[n][k] = 1;
c[n][1] = n;
c[n][0] = 1;
Finally, we construct a loop that fills in c[n][k] for all of the remaining entries. We have to be careful when constructing the loop to make sure that the formula for c[n][k] only uses entries that the loop has already filled in. Here is the correct loop structure to do this.
MAX_N = 100
MAX_K = 100
c = np.zeros(MAX_N,MAX_K)
def initTable():
# Fill in the base cases
for n in range(0,MAX_N):
c[n][k] = 1
c[n][1] = n
c[n][0] = 1
# Fill in the remaining entries
for n in range(2,MAX_N):
for k in range(2,n):
c[n][k] = c[n-1][k-1]+c[n-1][k]
for k in range(n+1,MAX_K)
c[n][k] = 0
Romberg integration
In an earlier lecture we used the trapezoid rule to estimate the area under a curve. The basic idea is to divide the range of integration into a large number of subintervals. Over each of these subintervals we construct a trapezoid whose area approximates the area under the curve over that subinterval.
The area of one such trapezoid is h (f(xi) + f(xi+1))/2. Adding all of these trapezoid areas up gives us the trapezoid area estimate for N subintervals:
The trapezoid rule forms the starting point for a more sophisticated method. The first step in this method is a technique called the Richardson extrapolation.
The Richardson extrapolation begins with a consideration of the error term in the trapezoid rule. When we deploy the trapezoid rule to estimate an integral with N subintervals and a step size of h, the trapezoid rule typically produces results that are accurate to a factor that is proportional to h2.
The notation O(h2) I used here means that the error is of order h2. Another way to write this fact is to write the error term as a combination of various powers of h:
(This is somewhat beyond the scope of the discussion here, but it does turn out that the error term from the trapezoid rule contains only even powers of h.)
The next step is deploy the Richardson extrapolation, which is a strange, but effective algebraic trick. In this trick we write down the error relationship twice, once for N trapezoids as we saw above, and a second time for 2 N trapezoids. Increasing the number of trapezoids by a factor of 2 reduces the step size from h to h/2, so we get these two expressions for the error:
Multiplying the second equation by 4 and subtracting the first equation from it produces
We can rewrite this as
or
By using this simple algebraic trick, we have switched from a method that has an error term with magnitude O(h2) to a method with an error term with magnitude O(h4).
Better yet, this trick can be iterated. I will spare you the details, but iterating this trick is the basis for an integration technique called Romberg integration. Here is a summary of the Romberg technique. The function R(k,j) represents the jth iteration of this technique starting from a trapezoid rule with N = 2k-1 steps.
R(a,b,k,1) = trapezoidArea(a,b,2k-1)
Here is a Python program that uses this recursive formula to compute integral estimates.
import math
# This is the function we will work with for this example
def f(x):
return math.sqrt(4.0 - x**2)
# Compute an approximation to the integral
# of f(x) from a to b by using trapezoids
# with steps of size (b-a)/N
def trapezoidArea(a, b, N):
h = (b - a) / N
x = np.arange(a, b, h)
return h * (f(x).sum()) - h * f(a) / 2 + h * f(b) / 2
# Recursively defined Romberg formula
def romberg(a,b,k,j):
if j == 1:
return trapezoidArea(a,b,2**(k-1))
elif k < j:
return 0
else:
return romberg(a,b,k,j-1) + (romberg(a,b,k,j-1)-romberg(a,b,k-1,j-1))/(4**(j-1)-1)
print(" j estimate error");
for j in range(2,6):
estimate = romberg(0.0,2.0,4*j,j)
error = math.fabs(math.pi - estimate)
print("{:>3d} {:f} {:.16f}".format(j,estimate,error))
It takes a little bit of effort to unpack what this program is doing. The loop at the bottom is computing R(a,b,4 j,j) for j in the range from 2 to 5. For a given value of j, this method starts with a trapezoid rule estimate using 2k-1 = 24 j - 1 steps. When j = 5 this means using 220 - 1 = 524288 steps.
This program has one last flaw in it, and this flaw seriously bogs down the computation starting at around j= 4. The problem here is one of excessive recursion. As you can see from the code above, a typical call to the romberg function triggers three additional recursive calls. If each of those in turn triggers more calls to the romberg function, you can end up making a huge number of calls this function.
To counteract this problem, we can use a simple caching strategy. Each time we compute a value for romberg(a,b,k,j) for a particular combination of k and j, we should make a note of that result. The next time we end up needing that value we can simply look up the result without having to recompute anything.
The final version of the program makes use of this caching strategy. We set up a two dimensional list R that can store values for many combinations of k and j. The romberg function then hands most of the work off to a recursive internal function, r_int. r_int does an optimized calculation by first checking the cache to see if a result has already been recorded for this combination of k and j. If so, it immediately stops and returns the cached value. If not, it falls back to the recursive formula to compute a result, and then caches that result before returning an answer.
import math
# This is the function we will work with for this example
def f(x):
return math.sqrt(4.0 - x**2)
# Compute an approximation to the integral
# of f(x) from a to b by using trapezoids
# with steps of size (b-a)/N
def trapezoidArea(a, b, N):
h = (b - a) / N
x = np.arange(a, b, h)
return h * (f(x).sum()) - h * f(a) / 2 + h * f(b) / 2
# Recursively defined Romberg formula with caching
# r_int will compute R[k,j] but then remember that
# result in the R list to speed up subsequent calculations.
def r_int(a,b,k,j,R):
if R[k][j] != -1: # Have we seen this k,j before?
return R[k][j] # If so, just return the cached value
# If not, compute, cache, and return the result
if j == 1:
R[k][j] = trapezoidArea(a,b,2**(k-1))
elif k < j:
R[k][j] = 0
else:
R[k][j] = r_int(a,b,k,j-1,R) + (r_int(a,b,k,j-1,R)-r_int(a,b,k-1,j-1,R))/(4**(j-1)-1)
return R[k][j]
def romberg(a,b,k,j):
R = [[-1 for c in range(j+1)] for r in range(k+1)]
return r_int(a,b,k,j,R)
print(" N estimate error");
for j in range(4, 27, 2):
estimate = romberg(0.0, 2.0, j, j)
error = np.fabs(np.pi - estimate)
print("{:>8d} {:f} {:.16f}".format(2**(j-1), estimate, error))
Note that since we have a more efficient method now we can safely expand the range of j that we compute with. Here are the results that this program returns.
N estimate error
8 3.124218 0.0173744893594270
32 3.139447 0.0021459042519050
128 3.141325 0.0002678879206544
512 3.141559 0.0000334785532052
2048 3.141588 0.0000041846140313
8192 3.141592 0.0000005230705571
32768 3.141593 0.0000000653836278
131072 3.141593 0.0000000081729490
524288 3.141593 0.0000000010216188
2097152 3.141593 0.0000000001277005
8388608 3.141593 0.0000000000159650
33554432 3.141593 0.0000000000019940
Programming Assignment
In one of the first lectures in this course we saw an example of the polynomial interpolation problem. Given a list (x0 , y0) , (x1 , y1) , … , (xn , yn) of points to interpolate, we want to construct an nth degree polynomial that passes through the points.
Newton solved this problem by looking for a polynomial that takes a special form, called the Newton polynomial:
Pn(x) = a0 + a1(x - x0) + a2(x - x0)(x - x1) + ⋯ + an(x - x0)(x - x1)⋯(x - xn-1)
Newton was able to come up with a recursive relationship that the coefficients of this polynomial satisfied. He introduced a special notation
Using this notation the coefficients of the Newton polynomial are
ak = f[0,k]
Write a program that can read a list of data points from a text file and construct a Newton polynomial that interpolates the points. Each line in the text file contains one x value and one y value. Your program should read those data points and store them in two lists, xs and ys. Your program will then prompt the user to enter a value for x and then compute and print the value of the Newton polynomial at that x.
Your program should contain at least one function, f(i,j), to compute the factors f[i,j]. You may also want to write a second function, term(k,x) to compute the kth term of the interpolating polynomial:
f[0,k] (x- x0)(x - x1)⋯(x - xk-1)
To test your program, here is some test data.
| x | y |
|---|---|
| 1.0 | 0.7651977 |
| 1.3 | 0.6200860 |
| 1.6 | 0.4554022 |
| 1.9 | 0.2818186 |
| 2.2 | 0.1103623 |
The polynomial of degree four that passes through these points has value approximately 0.5118200 at x = 1.5.