Out of the box
This tutorial
-
is related to Google’s or-tools example clustering_sat.py.
-
looks at how to divide a set of cities into four equal-sized subsets so that the sum of the distances inside each cluster is minimized.
-
asks whether the SAT-solver-based solution proposed by Google still scales when the number of cities grows.
-
shows what to do when it does not.
The code for this tutorial is here: clustering.py
Motivation
Google’s or-tools provide excellent implementations of many optimization algorithms. fcmaes also has a fast C++-based backend and supports multi-threaded optimization. The two libraries complement each other well. It is easy to build objective functions that use or-tools to solve sub-problems.
Still, it is worth being careful when classifying a problem too early. Small changes in problem size, or a few extra constraints, can move a problem into a different practical category.
The "dividing a set of cities into four equal sized subsets, so that the sum of the distances inside a cluster is minimized" problem seems to fit very well into the "sat-solver" category. But Paul Rulkens showed in Why the majority is always wrong that it is sometimes worth looking "outside the box".
Problem with the sat solver
The scaling issue becomes obvious when we increase the problem size from 40 to 200 cities:
...
repeat = 5
distance_matrix = np.array(distance_matrix)
distance_matrix = np.repeat(distance_matrix, repeat, 0)
distance_matrix = np.repeat(distance_matrix, repeat, 1)Running clustering_sat.py with this change first failed completely on a 16 core AMD 5950 machine with 128 GB RAM. The following modification helped:
solver.parameters.num_search_workers = 32Now the solver can use all 16 CPU cores. After about 90 minutes, and with 118 GB of memory in use, we get a result:
Starting CP-SAT solver v9.3.10497
...
CpSolverResponse summary:
status: OPTIMAL
objective: 82393100
best_bound: 82393100
booleans: 19900
conflicts: 0
branches: 39800
propagations: 0
integer_propagations: 39813
restarts: 39800
lp_iterations: 6783
walltime: 4694.31
usertime: 4694.31
deterministic_time: 21487.6
gap_integral: 35098.9
Group 0 : 0 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 31 32 33 34 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179
Group 1 : 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199
Group 2 : 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99
Group 3 : 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149The CP-SAT solver succeeds, but it is clearly at its limit here. Still, we now have a reference solution for testing an alternative method:
objective: 82393100Out of the box solution
Using continuous optimization for this clustering problem may be unusual in an fcmaes tutorial. Still, it is surprising how well such a simple approach works. A continuous derivative-free algorithm seems like a poor fit for the huge number of possible clusterings of 200 cities. For many algorithms, that intuition is correct (Exercise: check this statement).
This example shows a useful fcmaes pattern. We can encode a discrete assignment with continuous variables and let the optimizer search over the induced permutations.
Defining the objective function is simple:
@njit(fastmath=True)
def sum_distance(s):
dist = 0
g = np.empty((num_groups, group_size), dtype=numba.int32)
for i in range(num_groups):
g[i] = s[i*group_size:(i+1)*group_size]
for i in range(num_groups):
for v in range(group_size):
for w in range(group_size):
dist += distance_matrix[g[i,v],g[i,w]]
return dist
def unique(x):
return np.argsort(x)
def fitness(x):
return sum_distance(unique(x))This is much less work than building the CP-SAT solver configuration. We only need to remember one Python rule: "No nested inner loop without Numba". Otherwise the performance penalty is severe.
Calling the optimizer is straightforward too:
def opt():
print('Num nodes =', num_nodes)
res = retry.minimize(wrapper(fitness),
Bounds([0]*num_nodes, [1]*num_nodes),
optimizer=Bite_cpp(500000, stall_criterion = 3),
num_retries=32)For this kind of problem, we should try biteopt first.
retry.minimize automatically uses all available parallel threads (32 on the AMD 5950),
and we now get the solution in about 3.5 to 8 seconds with less than 2 GB of memory
usage. Note that we use the np.argsort trick to get a sequence of unique integers from
a list of continuous variables.
3.55 2239198 630760.0 82393100.0 [0.48553103356888927, 0.4942163924201614, 0.48884715084924707, ...
[[114 138 133 119 140 134 120 148 113 121 132 102 118 122 141 105 127 111
146 117 137 110 124 107 135 136 142 106 115 129 101 109 131 108 100 143
144 112 145 125 104 116 126 149 128 139 123 147 103 130]
[182 188 44 153 42 199 162 160 156 152 183 41 157 161 40 186 155 185
180 35 45 49 151 154 48 38 181 190 194 43 158 46 196 164 195 150
197 159 163 184 36 37 47 187 189 192 193 198 191 39]
[ 32 11 6 178 0 5 1 175 27 169 29 16 173 2 166 12 20 9
23 165 171 179 8 174 172 15 28 31 168 22 4 177 24 3 30 13
10 21 19 14 167 7 25 33 18 34 176 170 17 26]
[ 72 99 95 62 79 71 87 91 64 58 60 81 55 84 54 52 94 85
67 76 73 92 59 51 83 69 70 86 65 61 68 75 96 89 57 53
66 97 50 78 63 74 80 93 98 77 82 88 90 56]] 82393100Conclusion
-
Google’s or-tools provide excellent implementations of many important algorithms and complement fcmaes well.
-
Use the method that fits best, or a combination of methods, but do not classify a problem too early.
-
Small changes in problem size or extra constraints may change which method works best.
-
A good continuous optimizer can solve some discrete problems surprisingly well.
-
Use the
np.argsorttrick if you need a sequence of unique integers. -
Try biteopt first for this kind of problem.