An Ising Optimization algorithm for solving the 0-1 Knapsack Problem

CHANGELOG.md

@@ -25,6 +25,7 @@ Added
 -   Algorithm interface and result classes (#849) 
 -   Chemistry FCIDump file driver (#859)
 -   Chemistry stack automatic Z2 symmetry reduction (#870)
+-   Ising Optimization: The 0-1 Knapsack problem (#878)
 
 Changed
 -------

qiskit/optimization/ising/__init__.py

@@ -2,7 +2,7 @@
 
 # This code is part of Qiskit.
 #
-# (C) Copyright IBM 2019.
+# (C) Copyright IBM 2019, 2020.
 #
 # This code is licensed under the Apache License, Version 2.0. You may
 # obtain a copy of this license in the LICENSE.txt file in the root directory
@@ -36,6 +36,7 @@
    tsp
    vehicle_routing
    vertex_cover
+   knapsack
 
 Automatic Ising Model Generator from DoCPLEX Model
 ==================================================

qiskit/optimization/ising/knapsack.py

@@ -0,0 +1,225 @@
+# -*- coding: utf-8 -*-
+
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2020.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""
+Convert knapsack parameters instances into Pauli list
+The parameters are a list of values a list of weights and a maximum weight of the knapsack.
+
+In the Knapsack Problem we are given a list of objects that each has a weight and a value.
+We are also given a maximum weight we can carry. We need to pick a subset of the objects
+so as to maximize the total value without going over the maximum weight.
+
+If we have the weights w[i], the values v[i] and the maximum weight W_max.
+We express the solution as a binary array x[i]
+where we have a 1 for the items we take in the solution set.
+We need to maximize sum(x[i]*v[i]) while respecting W_max >= sum(x[i]*w[i])
+
+"""
+
+import logging
+import math
+import numpy as np
+
+from qiskit.quantum_info import Pauli
+from qiskit.aqua.operators import WeightedPauliOperator
+
+
+logger = logging.getLogger(__name__)
+
+
+def get_operator(values, weights, max_weight):
+    """
+    Generate Hamiltonian for the knapsack problem.
+
+    Notes:
+        To build the cost function for the Hamiltonian we add a term S
+        that will vary with our solution. In order to make it change wit the solution
+        we enhance X with a number of additional bits X' = [x_0,..x_{n-1},y_{n}..y_{n+m-1}].
+        The bytes y[i] will be the binary representation of S.
+        In this way the optimizer will be able to optimize S as well as X.
+
+        The cost function is
+        $$C(X') = M(W_{max} - \\sum_{i=0}^{n-1} x_{i}w_{i} - S)^2 - \\sum_{i}^{n-1} x_{i}v_{i}$$
+        where S = sum(2**j * y[j]), j goes from n to n+log(W_max).
+        M is a number large enough to dominate the sum of values.
+
+        Because S can only be positive, when W_max >= sum(x[i]*w[i])
+        the optimizer can find an S (or better the y[j] that compose S)
+        so that it will take the first term to 0.
+        This way the function is dominated by the sum of values.
+        If W_max < sum(x[i]*w[i]) then the first term can never be 0
+        and, multiplied by a large M, will always dominate the function.
+
+        The minimum value of the function will be that where the constraint is respected
+        and the sum of values is maximized.
+
+    Args:
+        values (list of non-negative integers) : a list of values
+        weights (list of non-negative integers) : a list of weights
+        max_weight (non negative integer) : the maximum weight the knapsack can carry
+
+    Returns:
+        WeightedPauliOperator: operator for the Hamiltonian
+        float: a constant shift for the obj function.
+
+    Raises:
+        ValueError: values and weights have different lengths
+        ValueError: A value or a weight is negative
+        ValueError: All values are zero
+        ValueError: max_weight is negative
+
+    """
+    if len(values) != len(weights):
+        raise ValueError("The values and weights must have the same length")
+
+    if any(v < 0 for v in values) or any(w < 0 for w in weights):
+        raise ValueError("The values and weights cannot be negative")
+
+    if all(v == 0 for v in values):
+        raise ValueError("The values cannot all be 0")
+
+    if max_weight < 0:
+        raise ValueError("max_weight cannot be negative")
+
+    y_size = int(math.log(max_weight, 2)) + 1 if max_weight > 0 else 1
+    n = len(values)
+    num_values = n + y_size
+    pauli_list = []
+    shift = 0
+
+    # pylint: disable=invalid-name
+    M = 10 * np.sum(values)
+
+    # term for sum(x_i*w_i)**2
+    for i in range(n):
+        for j in range(n):
+            coefficient = -1 * 0.25 * weights[i] * weights[j] * M
+            pauli_op = _get_pauli_op(num_values, [j])
+            pauli_list.append([coefficient, pauli_op])
+            shift -= coefficient
+
+            pauli_op = _get_pauli_op(num_values, [i])
+            pauli_list.append([coefficient, pauli_op])
+            shift -= coefficient
+
+            coefficient = -1 * coefficient
+            pauli_op = _get_pauli_op(num_values, [i, j])
+            pauli_list.append([coefficient, pauli_op])
+            shift -= coefficient
+
+    # term for sum(2**j*y_j)**2
+    for i in range(y_size):
+        for j in range(y_size):
+            coefficient = -1 * 0.25 * (2 ** i) * (2 ** j) * M
+
+            pauli_op = _get_pauli_op(num_values, [n + j])
+            pauli_list.append([coefficient, pauli_op])
+            shift -= coefficient
+
+            pauli_op = _get_pauli_op(num_values, [n + i])
+            pauli_list.append([coefficient, pauli_op])
+            shift -= coefficient
+
+            coefficient = -1 * coefficient
+            pauli_op = _get_pauli_op(num_values, [n + i, n + j])
+            pauli_list.append([coefficient, pauli_op])
+            shift -= coefficient
+
+    # term for -2*W_max*sum(x_i*w_i)
+    for i in range(n):
+        coefficient = max_weight * weights[i] * M
+
+        pauli_op = _get_pauli_op(num_values, [i])
+        pauli_list.append([coefficient, pauli_op])
+        shift -= coefficient
+
+    # term for -2*W_max*sum(2**j*y_j)
+    for j in range(y_size):
+        coefficient = max_weight * (2 ** j) * M
+
+        pauli_op = _get_pauli_op(num_values, [n + j])
+        pauli_list.append([coefficient, pauli_op])
+        shift -= coefficient
+
+    for i in range(n):
+        for j in range(y_size):
+            coefficient = -1 * 0.5 * weights[i] * (2 ** j) * M
+
+            pauli_op = _get_pauli_op(num_values, [n + j])
+            pauli_list.append([coefficient, pauli_op])
+            shift -= coefficient
+
+            pauli_op = _get_pauli_op(num_values, [i])
+            pauli_list.append([coefficient, pauli_op])
+            shift -= coefficient
+
+            coefficient = -1 * coefficient
+            pauli_op = _get_pauli_op(num_values, [i, n + j])
+            pauli_list.append([coefficient, pauli_op])
+            shift -= coefficient
+
+    # term for sum(x_i*v_i)
+    for i in range(n):
+        coefficient = 0.5 * values[i]
+
+        pauli_op = _get_pauli_op(num_values, [i])
+        pauli_list.append([coefficient, pauli_op])
+        shift -= coefficient
+
+    return WeightedPauliOperator(paulis=pauli_list), shift
+
+
+def get_solution(x, values):
+    """
+    Get the solution to the knapsack problem
+    from the bitstring that represents
+    to the ground state of the Hamiltonian
+
+    Args:
+        x (numpy.ndarray): the ground state of the Hamiltonian.
+        values (numpy.ndarray): the list of values
+
+    Returns:
+        numpy.ndarray: a bit string that has a '1' at the indexes
+         corresponding to values that have been taken in the knapsack.
+         i.e. if the solution has a '1' at index i then
+         the value values[i] has been taken in the knapsack
+    """
+    return x[:len(values)]
+
+
+def knapsack_value_weight(solution, values, weights):
+    """
+    Get the total wight and value of the items taken in the knapsack.
+
+    Args:
+        solution (numpy.ndarray) : binary string that represents the solution to the problem.
+        values (numpy.ndarray) : the list of values
+        weights (numpy.ndarray) : the list of weights
+
+    Returns:
+        tuple: the total value and weight of the items in the knapsack
+    """
+    value = np.sum(solution * values)
+    weight = np.sum(solution * weights)
+    return value, weight
+
+
+def _get_pauli_op(num_values, indexes):
+    pauli_x = np.zeros(num_values, dtype=np.bool)
+    pauli_z = np.zeros(num_values, dtype=np.bool)
+    for i in indexes:
+        pauli_z[i] = not pauli_z[i]
+
+    return Pauli(pauli_z, pauli_x)

test/optimization/test_knapsack.py

@@ -0,0 +1,75 @@
+# -*- coding: utf-8 -*-
+
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2020.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+""" Test Knapsack Problem """
+
+from test.optimization import QiskitOptimizationTestCase
+import numpy as np
+
+from qiskit.optimization.ising import knapsack
+from qiskit.optimization.ising.common import sample_most_likely
+from qiskit.aqua.algorithms import NumPyMinimumEigensolver
+
+
+class TestTSP(QiskitOptimizationTestCase):
+    """Knapsack Ising tests."""
+
+    @staticmethod
+    def _run_knapsack(values, weights, max_weight):
+        qubit_op, _ = knapsack.get_operator(values, weights, max_weight)
+
+        algo = NumPyMinimumEigensolver(qubit_op)
+        result = algo.run()
+        x = sample_most_likely(result.eigenstate)
+
+        solution = knapsack.get_solution(x, values)
+        value, weight = knapsack.knapsack_value_weight(solution, values, weights)
+
+        return solution, value, weight
+
+    def test_knapsack(self):
+        """ Knapsack test """
+        values = [10, 40, 50, 75]
+        weights = [5, 10, 3, 12]
+        max_weight = 20
+
+        solution, value, weight = self._run_knapsack(values, weights, max_weight)
+
+        np.testing.assert_array_equal(solution, [1, 0, 1, 1])
+        np.testing.assert_equal(value, 135)
+        np.testing.assert_equal(weight, 20)
+
+    def test_knapsack_zero_max_weight(self):
+        """ Knapsack zero max weight test """
+        values = [10, 40, 50, 75]
+        weights = [5, 10, 3, 12]
+        max_weight = 0
+
+        solution, value, weight = self._run_knapsack(values, weights, max_weight)
+
+        np.testing.assert_array_equal(solution, [0, 0, 0, 0])
+        np.testing.assert_equal(value, 0)
+        np.testing.assert_equal(weight, 0)
+
+    def test_knapsack_large_max_weight(self):
+        """ Knapsack large max weight test """
+        values = [10, 40, 50, 75]
+        weights = [5, 10, 3, 12]
+        max_weight = 1000
+
+        solution, value, weight = self._run_knapsack(values, weights, max_weight)
+
+        np.testing.assert_array_equal(solution, [1, 1, 1, 1])
+        np.testing.assert_equal(value, 175)
+        np.testing.assert_equal(weight, 30)