This example uses an algorithm to find a number's prime factors. For other interesting algorithms, see my book Essential Algorithms, Second Edition. The book describes many fascinating algorithms in pseudocode with source code available for download in Python and C#.

Title: Use dynamic programming to solve the traveling salesman problem

In dynamic programming, you save previously calculated values to make solving future values easier. This example builds a collection of carefully designed shortest paths that it can use to build shortest paths that include more points. In the end, it uses the saved paths to find a shortest path that visits every point and returns to the start point, thus solving the traveling salesman problem (TSP).

The Idea

First, let's start the tour at point 0. It doesn't really matter what point you start at because every point will be in the final tour.

Next, suppose S is a subset of the points that includes points 0 and j.

Let costs(S, j) represent the shortest path from point 0 to point j while visiting all of the points in S without visiting any point more than once.

The picture on the right shows the idea. Set S includes points 0 and j. In this example, there are also a few other points in the set and a few not in the set. A path starts at point 0, visits all of the points in set S, and ends at point j.

The algorithm builds the costs values from the bottom up by using values for smaller sets to calculate values for larger sets.

For the base case, the smallest set S that holds point 0 holds only point 0. In other words, S is {0}. Then costs(S, 0) equals 0 because a path from point 0 to point 0 has 0 length.

Now consider larger sets. For a set S holding more than just the point 0, the value costs(S, 0) is infinite. That value represents a path starting at point 0, visiting all of the points in set S, and then ending at point 0 without visiting any point twice. But if the path starts and ends with point 0, it visits that point twice so it's not allowed.

Now suppose we have a different set S that contains point 0 and some other points, one of which is point j. The shortest path from point 0 to point j must have some point, let's call it i, that comes right before i in the path. In other words, the path goes 0 → ... → i → j.

The shortest such path will have length equal to the shortest path from point 0 to point i plus the distance from point i to point j. Mathematically that distance is costs(S - {j}, i) + distance(i, j).

Now to find costs(S, j), we examine every possible value of i in S, calculate its value costs(S - {j}, i) + distance(i, j), and use the smallest of those values. The value of i that gives the smallest value is the last point in the shortest path from 0 to j.

Mathematically here's the equation.

The Algorithm

1. Set costs[(s, 0)] = (0, [0]).
2. For set_size = 2 to num_points:
   1. Loop through all possible sets S containing set_size points, one of which is point 0 and do the following:
      1. Set costs(S, 0) to infinity. (No path can start at point 0, visit other points, and return to point 0 without visiting point 0 twice.)
      2. Loop through points j ∈ S (not counting point 0) and do the following:
         1. Let s_minus_j be S - {j}.
         2. Loop through points i ∈ s_minus_j (again not counting point 0) and find the minimum of costs(s_minus_j, i) + distances(i, j).
         3. Save that minimum value as costs(S, j).
3. After calculating all of the shortest paths, loop through all points i and find the minimum of costs(S, i) + distances(i, 0). This is the path from 0 to i and then back to 0.

The shortest of the paths from 0 to i and back to 0 is an optimal tour.

The Code

import itertools def make_dynamic_tour_generator(num_points, distances, **kwargs): '''Use dynamic programming to return an optimal tour.''' global num_calls num_calls = 0 # costs is a dictionary where costs[(s, j)] returns a tuple giving # the length and best tour that visits all points in set s, # starting at point 0 and ending at point j. # Actually s will be a frozenset so it and therefore (s, j) are immutable # and can thus be used as dictionary keys. costs = {} # c[{0}, 0] = 0 s = frozenset({0}) costs[(s, 0)] = (0, [0]) # Make a list of all point indices *except 0*. indices = list(range(1, num_points)) # Build values for larger sets. for set_size in range(2, num_points + 1): # Loop through all sets of size set_size. for subset in itertools.combinations(indices, set_size - 1): # Make a frozenset holding the subset plus 0. s = frozenset(subset + (0,)) # A tour cannot both start and end at point 0 if set_size > 1. costs[(s, 0)] = (math.inf, None) # Loop through subsets minus j where j != 0. for j in s: # Skip j == 0. if j == 0: continue s_minus_j = s - {j} # Find the shortest tour costs[s_minus_j, i] --> j. best_length = math.inf best_tour = None for i in s_minus_j: test_length = costs[(s_minus_j, i)][0] + distances[i][j] if test_length < best_length: # This is shorter than the best so far. Save it. best_length = test_length best_tour = costs[(s_minus_j, i)][1] + [j] # Save the result. costs[(s, j)] = (best_length, best_tour) # Print the costs matrix. Only do this for small examples! #for value in costs: # print(f'{value}: {costs[value]}') # Find an optimal tour. best_length = math.inf best_tour = None all_set = frozenset(range(num_points)) for i in indices: # Check the path 0 ... i --> 0. test_length, test_tour = costs[(all_set, i)] test_length += distances[i][0] if test_length < best_length: best_length = test_length best_tour = test_tour print(f'costs size: {len(costs)}', end='') yield best_length, best_tour

First, the program makes costs a dictionary. Its keys are tuples holding a frozenset and a value j so costs(S, j) is as described in the algorithm. A dictionary's keys must be hashable and normal sets are not, so the program uses frozensets, which are hashable.

Second, the dictionary's values are tuples containing both a path's length and a list holding the points in the path. For example, costs(S, j) holds the length of the path from point 0 to point j plus a list giving the shortest path. (Most of the examples I've seen online only hold the shortest distance. Many of the examples I've seen are also confusing and/or just plain wrong.)

Next, the code uses itertools.combinations to generate sets of a desired size. For example, itertools.combinations(a_list, 3) generates all of the three-item subsets of the items in a_list. You could write some code to generate the combinations yourself, but this algorithm is confusing enough already.

Finally, to generate a set of size K containing point 0, the program uses itertools.combinations to find subsets of K - 1 points in the range 1 to num_points and then adds point 0 to the result.

Notice the commented code after the comment Print the costs matrix. Only do this for small examples! That code prints the costs values. They're interesting but be warned that the costs dictionary grows extremely quickly as the problem size grows so printing it out can take a long time. The output is still reasonable for 4 points, but I wouldn't print the costs for much larger problems.

Conclusion

All of that means the costs dictionary grows large quickly and it takes a while to build it.

The following table shows the times in some trials plus the size of the costs dictionary.

As you can see, the dynamic programming algorithm is much faster than the other approaches. The difference should grow larger as the problem size increases but at some point the costs dictionary starts to use a significant fraction of the system's total memory and that slows things down.

Keep in mind that the costs keys are not simple numbers but rather frozensets and numbers so they take up quite a bit of room. The values are also lengths and lists of points, so they take up a bunch of room, too. To solve the 20-point problem, the program peaked at around 2.1 GB of memory. If you had special-purpose hardware with tons of memory, you could probably solve bigger problems faster, but the limit is clearly approaching on my computer.

Download the example to experiment and see additional details.