Knapsack branch and bound: forward filter

Question

I've coded branch and bound to solve the knapsack problem, and use a greedy linear relaxation to find an upper bound on a node I'm currently exploring. What I mean by this is that I first sort the items by their respective value density (computed value/weight, for each item) and then fill the bag according to a greedy algorithm, using a fraction of the last item that would put us over the bag's capacity. This solution is unattainable (we can't take half an item) but provides an upper bound on the maximum value we can attain from considering an item j:

def _linear_relaxation(self, j, cp, cw):
    # we consider how well we can do to increase our current profit (cp) 
    # by looking at items beyond the item that we either just took or left behind 

    weight = cw
    value  = cp

    for (v, w), index in self.ordered:
        if index > j:

            if w + weight <= self.capacity:
                value  += v
                weight += w
            else:
                return value + (v * (self.capacity - w)/float(w))

    return value

My code works great to solve small problems using this relaxation, but starts slowing down when the problem size increases. One optimization (or reduction on a node's upper bound) is to avoid taking items of high value density whose weight exceeds the leftover weight we're trying to fill.

Think about it this way: our greedy linear relaxation bound takes items in decreasing order according to their density, value/weight. If the residual or remaining weight in the bag is less than the item we're considering adding to the bound, it is much better to take items further down our density list. Consider, for example, residual = 1, weight = 2. Wouldn't it be great if we found a small value, so that we could add to our estimate value/2, rather than 2 * value? We'd get a much better estimate on the least upper bound than if we considered the last item with the greatest density, with no regard as to its relative effect on our upper bound.

With this in mind, a better bound can be found as follows:

def _linear_relaxation(self, j, cp, cw):

    weight = cw
    value  = cp

    for (v, w), index in self.ordered:

        if index > j:

            if w + weight <= self.capacity:
                value  += v
                weight += w

            else:
                if w > self.capacity - (w + weight):
                    continue

                return value + (v * (self.capacity - w)/float(w))

    return value

This bound is great for a list of no more than 1000 items (around 11 seconds compared to no solution) but fails on a problem with 10,000 items. My question, how can we we do better? Can the bound be improved on? Or can we more smartly consider what nodes are worth adding, based on the best solution discovered so far.

My question/goal of the review:

To provide some context, the code I've written is included below. Please let me know if you have improvements on its structure (for example, is it better practice to make TreeNode a truly recursive object and give it left and right attributes? Do we gain anything by doing this?). What I'm really looking for, though, is how we can improve the current depth first search implementation (please don't provide a solution that uses a heap, for example, unless you can demonstrate it is an all around better algorithm than what I've written). If you'd like to provide code, that would be great, but I'm really just looking for some ideas on how my current code might be optimized. I'm kind of new to optimization and don't really know where to look for improvements or really how to think about the problem as might a person with a speciality in optimization.

While searching, if we include an item, do we need to recalculate the bound on this node (or would it just be the value of its parent--this solution is still attainable)?

One optimization I should make: put the right node on after putting on the left node. The order in which we consider expanding a path matters, here!

class TreeNode:

    def __init__(self, level, V, W, taken):

        self.level = level
        self.V = V
        self.W = W
        self.taken = taken

    def __str__(self):
        return str((self.level, self.V, self.W, self.taken))


class KP:

    def __init__(self, cap, values, weights):

        self.capacity = cap
        self.values   = values
        self.weights  = weights

        # sort the items according to their value per weight 
        self.ordered = sorted([((v, w), i) for i, (v, w) in enumerate(zip(self.values, self.weights))],
                          key = lambda tup: float(tup[0][0])/tup[0][1], reverse = True)

    # our best initial solution will be the greedy solution 
    def _greedy_solution(self):

        weight = 0
        value  = 0
        taken  = []

        for (v, w), index in self.ordered:

            if w + weight <= self.capacity:   # take items of highest value per unit weight, while we still can
                value  += v
                weight += w
                taken.append(index)

            else:
                return TreeNode(index, value, weight, taken)

        return TreeNode(index, value, weight, taken)


    def solve(self):

        best = self._greedy_solution()
        root = TreeNode(0, 0, 0, [])

        stack = [root]

        while len(stack) > 0:

            current = stack.pop()    
            index   = current.level

            if current.V >= best.V:         # consider whether we can immediately improve our best
                best = current

            if index < len(self.values):    # there are items left to consider


                # we have two branches to consider: if we have room, we take the item
                # if we don't have room, we continue our solution to the next iteration 
                # where we consider the next item

                if current.W + self.weights[index] <= self.capacity:
                    taken = list(current.taken)     # update path of items taken
                    taken.append(index)
                                                    # create the right node
                    right = TreeNode( index + 1, current.V + self.values[index],
                                  current.W + self.weights[index], taken )

                    if right.V > best.V:            # update our best, if possible 
                        best = right

                    bound = self._linear_relaxation(index, right.V, right.W)

                                                    # only continue exploring if 
                                                    # we're guaranteed something better
                    if bound >= best.V:
                        stack.append(right)

                                                    # create the left node, continue
                                                    # solution of the parent, without             
                                                    # taking the item 

                left  = TreeNode( index + 1, current.V, current.W, list(current.taken) )
                bound = self._linear_relaxation(index, left.V, left.W)

                                                    # only continue exploring if 
                                                    # we're guaranteed something better
                if bound >= best.V:                 
                    stack.append(left)


        return best 

Answer 1

Your first upper bound is buggy

return value + (v * (self.capacity - w)/float(w))

should be:

return value + (v * (self.capacity - weight)/float(w))

If w_i << W you are overestimating the upper bound by heaps, which causes your branch and bound to never "bound".

Your second relaxation is not an upper bound

I'm assuming your linear relaxation is for the upper bound. As such your first function to calculate the upper bound is good (bar the bug above). The second function is no longer an upper bound because you are essentially taking a greedy solution as the upper bound.

Example that triggers faulty upper bound:

Define the value and weight of item i as: v_i and w_i respectively. Define the ratio of value to weight for item i as: r_i = v_i / w_i. The knapsack can carry at most W units of mass.

We have three items as follows:

r_0 = k
r_1 = r_2 = k - dk
w_0 = W - dw
w_1 = w_2 = W/2

where k is an arbitrary constant, dk is a small value and dw is a value strictly smaller than W/2 i.e: dw < W/2.

Your upper bound algorithm would pick i=0 first and then be unable to add any other items to the knapsack and the achieved value is:

   V_ub = r_0*w_0 = k*(W-dw) = k*W - dw*K

However the highest attainable value is achieved by taking r_1 and r_2 but not r_0:

   V_optimal = r_1*w_1 + r_2*w_2 = 2*r_1*w_1 = 2*(k-dk)*W/2 = k*W - dk*W

It is clear that when dw*K > dk*W then V_ub < V_optimal and your algorithm will cause your branch and bound algorithm to prune sections of the graph that could contain the optimum.

If you are uncertain of the above inequality you can let dw approach W/2 from 0 and let dk approach 0 from k and see that the limit is trivially:

    W/2*K > 0

which proves that this scenario is possible and likely to happen when you get enough items.

In comparison, the first upper bound would achieve:

V_ub = r_0*w_0 + r_1*(W-w0) = k*(W-dw) + (k-dk)*(W-(W-dw)) = k*W - dw*K + dw*k - dk*dw
     = k*W -dk*dw

As dw < W/2 then V_ub > V_opt and it is proved that this is an upper bound.

Edit/Addendum: Improve your lower bound

You initialize a best node with a greedy solution at the start here: Then you update the best solution here (similarly for left node):

                if right.V > best.V:            # update our best, if possible 
                    best = right

Then you calculate the upper bound with the linear relaxation method:

                bound = self._linear_relaxation(index, right.V, right.W)

And then only if a branch has a upper bound that is better than the best solution's partial value do you continue evaluating the branch:

                if bound >= best.V:
                    stack.append(right)

This efficiently make the best solution's partial value your lower bound. The lower and upper bounds should as tightly as possible bound the best achievable value in that branch. Taking the partial value in the branch this far, grossly underestimates the lower bound. Making your search space larger. Remember, your search space grows with your bounds.

You set the lower bound here (similarly for the left case):

                right = TreeNode( index + 1, current.V + self.values[index],
                              current.W + self.weights[index], taken )

Right.V is used as a lower bound which it technically is, albeit a poor one.

Instead of:

                if right.V > best.V:            # update our best, if possible 
                    best = right

You should do:

    # Somewhere earlier,
    best = self._greedy_solution()
    best_lb = best.V; # Here V is a lower bound on the whole search space. Good!

                # Later on...
                right = TreeNode( index + 1, current.V + self.values[index],
                              current.W + self.weights[index], taken )
                lowerBound = greedy_solution_on_remaining_items(right)
                if lowerBound > best_lb:            # update our best, if possible 
                    best = right
                     best_lb = lowerBound 

And similar checks for the left case.

Answer 2

I would add this as a comment, but as I am new in here, I can't. Mathematically speaking, this piece of code if w > self.capacity - (w + weight): may not make sense, since you already know that w + weight <= self.capacity is false. If my math is not failing me, the "if" statement is always true, since self.capacity - (w + weight) is smaller than zero, since w + weight > self.capacity, and w > 0.

Am I right?

And yet it doesn't give me the optimal solution.

Answer 3

Readability leads to efficiency

My first recommendation is readability. More readable code makes it easier (for you or others) to understand the algorithm and thus easier to think of efficiency improvements. You are free to choose whatever line length suits you, but if you want other people to be able to read your code easily 79 characters is a good guide.

Note: code blocks on Stack Exchange will require horizontal scroll bars beyond 91 characters. Some people (including me...) may be less likely to review code that requires scrolling left and right, when there is other code available to review without the inconvenience.

Blank lines can help readability by separating one block of code lines from another. However, using blank lines excessively means they lose their separating purpose and just stretch your code, actually reducing readability. How much is too much and how much is too little is a personal decision, but this guidance may help.

It's easier to understand a function if it fits on one page. If you have to scroll up and down to view a function then ask yourself if it could be broken down into smaller functions. In this particular case your solve method would fit on one page if it had fewer blank lines. You could remove some of the blank lines, or split the function into smaller ones. I recommend you do both.

---

Algorithm fine tuning?

Consider this method:

    def _greedy_solution(self):

        weight = 0
        value  = 0
        taken  = []

        for (v, w), index in self.ordered:

            if w + weight <= self.capacity:
                value  += v
                weight += w
                taken.append(index)

           else:
               return TreeNode(index, value, weight, taken)

        return TreeNode(index, value, weight, taken)

The highlighted else block here causes the method to return at the first occurrence of an item that would exceed capacity. Is this what you require? I'm not familiar with the algorithm but I expected that the greedy algorithm would continue considering the less dense items to add any which do not exceed capacity. This means the remaining space will always be used for the most dense of the remaining items that fit. Since the items are sorted according to density (that is, value/weight) the first item that exceeds capacity is not necessarily the last item worth considering. There could still be an item with lower density but smaller weight which would fit.

If I've understood correctly, removing the highlighted else block will give you a potentially higher value for your greedy solution, so you can start your solve method with a slightly higher lower bound (best). This should allow you to discard more of the potential branches later on, saving you unnecessary calculation.

---

A look at the solve method

Here is your solve method with some of the blank lines and spaces removed so we can see it all at once for discussion:

    def solve(self):
        best = self._greedy_solution()
        root = TreeNode(0, 0, 0, [])
        stack = [root]

        while len(stack) > 0:
            current = stack.pop()
            index   = current.level
            if current.V >= best.V:
                best = current

            if index < len(self.values):
                if current.W + self.weights[index] <= self.capacity:
                    taken = list(current.taken)
                    taken.append(index)
                    right = TreeNode(index + 1, current.V + self.values[index],
                                  current.W + self.weights[index], taken)

                    if right.V > best.V:
                        best = right

                    bound = self._linear_relaxation(index, right.V, right.W)
                    if bound >= best.V:
                        stack.append(right)

                left  = TreeNode(index + 1, current.V, current.W, list(current.taken))
                bound = self._linear_relaxation(index, left.V, left.W)
                if bound >= best.V:
                    stack.append(left)
        return best 

Consider this section:

                    if right.V > best.V:
                        best = right

                    bound = self._linear_relaxation(index, right.V, right.W)
                    if bound >= best.V:
                        stack.append(right)

If you have already assigned best = right then am I right in thinking that bound >= best.V will always be True? If so you can avoid the call to _linear_relaxation each time right is the new best, like this:

                    if right.V > best.V:
                        best = right
                        stack.append(right)
                    else:
                        bound = self._linear_relaxation(index, right.V, right.W)
                        if bound >= best.V:
                            stack.append(right)

Before adding left or right to your tree, you check if bound >= best.V. Since bound is an upper bound on the possible value obtainable through this branch, it is guaranteed that the branch will not find a better solution that the current best if bound = best.V. For those occasions when equality occurs, abandoning the branch will save you unnecessary work. So I would recommend replacing bound >= best.V with bound > best.V, before appending right and before appending left. There is no point in keeping a branch which can at best provide an equal value solution, since your program is not designed to return multiple solutions.