Zeta-Puppis.com

For a prob­lem I’m work­ing on I got stuck onto the clas­si­cal sit­u­a­tion of local max­i­mum. After trying to work around the prob­lem in sev­eral more or less cre­ative ways, I thought of the sim­u­lated anneal­ing algo­rithm. Con­sid­er­ing it’s been a while since I last saw it I tried to search for it on the web and sur­pris­ingly there is not much stuff about it, and the few bits I found are often con­trad­dic­tory. After quite a lot of dig­ging I decided to write about it here. As a warn­ing I should prob­a­bly say that there will be dig­ging into some basic sta­tis­tic and com­plex­ity analy­sis, as well as a quick formal intro­duc­tion to the prob­lem of the knap­sack. You should be able to follow even if you don’t know noth­ing about those topics, but having some foun­da­tions in these areas would be of great help.

Let begin with the knap­sack prob­lem. This is a clas­sic com­bi­na­to­r­ial com­puter sci­ence prob­lem known to be NP-​complete, mean­ing that the exact opti­mal solu­tion cannot be found in poly­no­mial time. This often means that most of the times we are happy of a good solu­tion, assum­ing it’s not so far from the opti­mal one. In the sim­pli­est pos­si­ble terms you are a thief and you’re in a room with a set of objects that are worth some­thing but you have only one knap­sack, and that knap­sack can carry at most a cer­tain weight, so you have to choose care­fully what objects to steal in order to max­i­mize the earn­ings. For exam­ple con­sider the fol­low­ing sit­u­a­tion: you can carry at most 5kg, and there is one laptop and a 4kg safe with pure dia­monds within. You can’t carry both of them so you have to choose what’s better to carry on, the dia­monds or the laptop. A smart thief would choose the dia­monds since their value is con­sid­er­ably higher than the laptop.

There are few vari­a­tions of the same prob­lem but most common one is named “0-1”: you can’t split the weight over two or more car­ri­ers or bags but either you take the whole weight or you leave the object where it is. Math­e­mat­i­cally talk­ing, con­sider a set of objects, each item is worth and weights with . Then the goal is to max­i­mize the fol­low­ing function:

But keep­ing the fol­low­ing con­straint (being the max­i­mum weight we can carry):

The first prob­lem we have to face then is how we should gen­er­ate the items and how the value of a solu­tion should be cal­cu­lated. The fol­low­ing exam­ple, as the other that will follow, it’s writ­ten in Python but it’s quite easy to under­stand so port­ing to another lan­guage wouldn’t be that hard. I chosen to gen­er­ate 50 objects with values that range from 1 to 99$ (both ends included) using a uni­form dis­tri­b­u­tion (if you don’t know much about sta­tis­tic, it means that all the values are equally dis­trib­uted among the objects). The same with the weights except they range from 1 to 20 (the choice of the weight’s unit mea­sure is left to you).

def generate_items(n_items=100):
    "Generate a list of items that could be stealed"
    items = []
    for n in range(0, n_items):
        # use a uniform distribution both for values and for weights
        cost = random.randint(1, 100)
        weight = random.uniform(1, 20)
        items.append((cost, weight))
    return items

So the items are noth­ing other than a list of pairs in the format for every object . Prob­a­bly a better and more real­is­tic dataset would have used a normal dis­tri­b­u­tion for the weights but it’s triv­ial to change the gen­er­a­tion to func­tion to work in that way. For our pur­poses the uni­form dis­tri­b­u­tion does its job.

As we said above the prob­lem is NP-​complete so we usu­ally need to visit the whole search space to get the opti­mal solu­tion which can be quite big when as the number of objects grows. Here comes the sim­u­lated anneal­ing. We won’t visit the search space exten­sively, but we’d rather gen­er­ate solu­tions. Indeed, it is a sto­chas­tic heuris­tic search algo­rithm. An heuris­tic is a func­tion that mea­sures how much some­thing is good or bad, and sto­chas­tic means that we move more or less in a random way into the search space. In prac­ti­cal terms it’s not greedy as that it doesn’t always follow what the heuris­tic says but rather ran­domly search where the heuris­tic func­tion points to. For exam­ple con­sider the needle in the haystack sit­u­a­tion: an exaus­tive search method would take every straw piece, check that it’s not a needle, put it apart and repeat those moves until you don’t find the needle. You don’t want to pro­ceed in this way, you’re more likely to end in less time if you look ran­domly in the haystack and if some­things stings you while you’re hold­ing straw then search into that straw piece, because there may be the needle in there. In the knap­sack prob­lem we do have the heuris­tic, and it’s the func­tion above that, for each solu­tion (admis­si­ble or not), says how much it’s worth. In this case we define an admis­si­ble solu­tion as one that it’s not too heavy.

An useful (and clas­sic as well) exam­ple is the one of the blind hill climb­ing. You’re blind and stuck on a hill and you need to reach the top. A good prin­ci­ple that could lead you to the top is to touch the ter­rain and always follow the rising path. It could, because if you’re on a rock than you surely haven’t reached the top but the prin­ci­ple above doesn’t apply: you reached a local max­i­mum (invert the things and you get the same thing for a local min­i­mum). Sim­u­lated anneal­ing avoids these prob­lems by trying wors­en­ing moves from time to time: even if this may not sound like a good move it helps avoid­ing the prob­lems we described above.

In the func­tion above once in the middle we could choose to take the left max­i­mum (which is a local max­i­mum). Using hill climb­ing we’d be stuck on that because we wouldn’t try other paths.

Sim­u­lated anneal­ing takes its name from the same process that metals go through when cool­ing from a melt­ing point. Indeed, the cool­ing process con­sists of sev­eral par­ti­cles that changes energy states (this state­ment may not be accu­rate or be inex­act at all, but please for­give me as I never stud­ied those things and I all know in this field comes from sim­u­lated anneal­ing algo­rithm itself), in par­tic­u­lar we can cal­cu­late the tran­sic­tion prob­a­bil­ity from one state to another. Con­sid­er­ing two energy states and and a tem­per­a­ture , switch­ing from to has probability:

Where is a con­stant called Boltzmann’s con­stant. But then, how do we apply those state­ments to our prob­lem (or a com­bi­na­to­r­ial search prob­lem, in gen­eral)? The most impor­tant con­cept to grasp is the energy switch­ing one. As a par­ti­cle change state, a solu­tion might change. Indeed for the knap­sack prob­lem there are many admis­si­ble solu­tion, each one with an asso­ci­ated earn­ing. Of course we’d prefer the one with the higher earn­ings (and sim­u­lated anneal­ing will help us find that) but still it’s per­fectly accept­able to go from a solu­tion to another as long as the other solu­tion con­tin­ues to be admissible.

So here it is what the sim­u­lated anneal­ing does: if you find a better item go on an take that path (under this cir­cum­n­stance, behaves just like the hill climb­ing), oth­er­wise change state with prob­a­bil­ity (you may notice that the Boltzmann’s con­stant is miss­ing, indeed that con­stant applies mostly to ther­mo­dy­namic when deal­ing with dif­fer­ent metals). The effect that the tem­per­a­ture scal­ing has is that at higher tem­per­a­tures it’ll try wors­en­ing moves quite often while on lower tem­per­a­tures the prob­a­bil­ity to do wors­en­ing moves is lower so when tem­per­a­ture tends towards 0 it behaves quite like the hill climb­ing algorithm.

def simulated_annealing(solution, items, max_weight):
    "Apply the simulated annealing for solving the knapsack problem"
    best = solution
    best_value = compute_cost(solution, items)[0]
    current_sol = solution
    temperature = 1.0
 
    while True:
        current_value = compute_cost(best, items)[0]
 
        for i in range(0, COOLING_STEPS):
            moves = generate_moves(current_sol, items, max_weight)
            idx = random.randint(0, len(moves) - 1)
            random_move = moves[idx]
 
            delta = compute_cost(random_move, items)[0] - compute_cost(best, items)[0]
 
            if delta > 0:
                best = random_move
                best_value = compute_cost(best, items)[0]
                current_sol = random_move
            else:
                if math.exp(delta / float(temperature)) > random.random():
                    current_sol = random_move
 
        temperature = TEMP_ALPHA * temperature
        if current_value >= best_value or temperature <= 0:
            break

And finally, that is the sim­u­lated anneal­ing. You start from a tem­per­a­ture of 1.0 then you have a cer­tain number of cool­ing steps, in every one of them you extract a random item from the neigh­bours and, if the item is better than the cur­rent best item then it becomes the new best item (and the new local solu­tion). If the new item’s value is worst than the cur­rent best then update the cur­rent local solu­tion with the prob­a­bil­ity expressed above. After the cool­ing steps the tem­per­a­ture is decreased with an expo­nen­tial decay (usu­ally it is ).

In the exam­ple above you don’t wait for the tem­per­a­ture to be 0 but you leave the loop if after all the cool­ing steps there hasn’t been any improve­ment. How big must be the number of cool­ing steps and the value it’s a fine tuning prob­lem. A dif­fer­ent approach that allow to get rid of the cool­ing steps is to make the tem­per­a­ture get cold slowly ( and is a very small value like 0.01). Besides, a common improve­ment is to cache the moves within the cool­ing steps unless you found a new best value or changed state.

You can see that the most crit­i­cal points are the neigh­bour gen­er­a­tion and the cost com­pu­ta­tion. While the neigh­bour gen­er­a­tion could be cached like I said above, the cost com­pu­ta­tion could be replaced with a prob­a­bil­ity esti­mate in order to reduce the time per cool­ing step.

But how do you apply the algo­rithm? If an empty solu­tion is accept­able then you can just start with that and let the neighbour’s gen­er­a­tor to create a solu­tion, but usu­ally you start with a greedy solu­tion (found through the hill climb­ing, for exam­ple) or from a random one.

Here fol­lows the com­plete code. I used 1000 cool­ing steps and a value of 0.8. I start from a random solu­tion whose behav­iour is not that bad con­sid­er­ing how much time it takes to com­pute the solu­tion. Indeed often random algo­rithms per­form really well on com­bi­na­to­r­ial algo­rithms, see sto­chas­tic search for some more informations.

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
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
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
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120

#!/usr/bin/env python
import math
import operator
import pprint
import random
import sys
 
COOLING_STEPS = 1000
TEMP_ALPHA = 0.8
 
random.seed()
 
def generate_items(n_items=100):
    "Generate a list of items that could be stealed"
    items = []
    for n in range(0, n_items):
        # use a uniform distribution both for values and for weights
        cost = random.randint(1, 100)
        weight = random.uniform(1, 20)
        items.append((cost, weight))
    return items
 
def main(args):
    items = generate_items()
    pprint.pprint(items)
 
    start_sol = generate_random_solution(items, max_weight=40)
    print "Random solution: %s" % start_sol
    print "value: (cost: %d, weight: %f)" % compute_cost(start_sol, items)
 
    solution = simulated_annealing(start_sol, items, max_weight=40)
    print "Final solution: %s" % solution
    print "value: (cost: %d, weight: %f)" % compute_cost(solution, items)
 
    return False
 
def generate_random_solution(items, max_weight):
    "Generate a starting random solution"
 
    # generate a random solution by adding a random item
    # until we don't get over the weight
    solution = []
    while compute_cost(solution, items)[1] <= max_weight:
        idx = random.randint(0, len(items) - 1)
        # skip duplicates
        if idx not in solution:
            solution.append(idx)
    # last item makes us get over the weight so simply remove it
    # we'll look for better results after
    solution = solution[:-1]
    return solution
 
def simulated_annealing(solution, items, max_weight):
    "Apply the simulated annealing for solving the knapsack problem"
    best = solution
    best_value = compute_cost(solution, items)[0]
    current_sol = solution
    temperature = 1.0
 
    while True:
        current_value = compute_cost(best, items)[0]
 
        for i in range(0, COOLING_STEPS):
            moves = generate_moves(current_sol, items, max_weight)
            idx = random.randint(0, len(moves) - 1)
            random_move = moves[idx]
 
            delta = compute_cost(random_move, items)[0] - \
                    compute_cost(best, items)[0]
 
            if delta > 0:
                best = random_move
                best_value = compute_cost(best, items)[0]
                current_sol = random_move
            else:
                if math.exp(delta / float(temperature)) > random.random():
                    current_sol = random_move
 
        temperature = TEMP_ALPHA * temperature
        if current_value >= best_value or temperature <= 0:
            break
    return best
 
def generate_moves(solution, items, max_weight):
    """
    Generate all the ammissible moves starting from the input
    solution
    """
    moves = []
    # try to add another item and save as a possible move
    for idx, item in enumerate(items):
        if idx not in solution:
            move = solution[::]
            move.append(idx)
 
            if compute_cost(move, items)[1] <= max_weight:
                moves.append(move)
 
    # try to remove one item
    for idx, item in enumerate(solution):
        move = solution[::]
        del move[idx]
        if move not in moves:
            moves.append(move)
 
    return moves
 
def compute_cost(solution, items):
    """
    Return a tuple in the format (id_item1, id_item2, ...)
    for the input solution
    """
    cost, weight = 0, 0
    for item in solution:
        cost += items[item][0]
        weight += items[item][1]
    return (cost, weight)
 
if __name__ == '__main__':
    sys.exit(main(sys.argv))

The results are supris­ing. For three dif­fer­ent set of 100 items, each one with its own value and weight, here are the results:

$ python sa.py
Random solution: [98, 71, 95]
value: (cost: 44, weight: 27.001685)
Final solution: [71, 95, 67, 9, 41, 33, 27]
value: (cost: 229, weight: 39.791386)

$ python sa.py
Random solution: [38, 16, 62, 31]
value: (cost: 124, weight: 36.863846)
Final solution: [38, 16, 62, 31, 5]
value: (cost: 194, weight: 38.970745)

Random solution: [61, 44, 48, 38]
value: (cost: 293, weight: 30.357135)
Final solution: [61, 44, 48, 38, 37, 5, 2]
value: (cost: 421, weight: 39.331549)

We usu­ally don’t leave much free weight and the qual­ity of the solu­tions found is quite good con­sid­er­ing that the time required to com­pute that is ~0.5s.