Using Python numerical libraries

Topics covered in this chapter:

The NumPy library allows for creating and using arrays of any order and size for efficiency and compactness purposes. This chapter shows how to take advantage of the features of NumPy in the creation of optimization problems. The Xpress Python interface supports NumPy versions 1.19 and later.

Using NumPy with the Xpress Python interface

NumPy arrays can be used as usual when creating variables, functions (linear and quadratic) of variables, and constraints. All functions described in this manual that take lists or tuples as arguments can take array's, i.e., NumPy array objects, as well, as in the following example:

import numpy as np
import xpress as xp
N = 20
S = range(N)
p = xp.problem()
x = xp.array([p.addVariable()                  for i in S])
y = xp.array([p.addVariable(vartype=xp.binary) for i in S])
constr1 = x <= y
p.addConstraint(constr1)

The above script imports both NumPy and the Xpress Python interface, then declares two arrays of variables and creates the set of constraints x_i ≤ y_i for all i in the set S.

The arrays must be of class xpress.ndarray in order to use the matrix/vector forms of the comparison (<=, ==, >=) operators. This is to avoid the default behaviour of NumPy arrays, where an expression like a <= b will result in a Boolean array, instead of an array of constraint objects. Ways to create instances of this class are described at xpress.ndarray.

NumPy allows for multiarrays with one or more 0-based indices. Given that declaring a NumPy multiarray of variables can result in a long line of code, the problem.addVariables function in its simplest usage returns an xpress.ndarray of variables with one or more indices. Consider the following three array declarations:

import numpy as np
import xpress as xp
p = xp.problem()
x = xp.array([p.addVariable(name='v({0})'.format(i)) for i in range(20)]).reshape(5,4)
y = xp.array([p.addVariable(vartype=xp.binary) for i in range(27)]).reshape(3,3,3)
z = xp.array([p.addVariable(lb=-1, ub=1) for i in range(1000)])

These can be written equivalently in the compact form as

import numpy as np
import xpress as xp
p = xp.problem()
x = p.addVariables(5, 4, name='v')
y = p.addVariables(3, 3, 3, vartype=xp.binary)
z = p.addVariables(1000, lb=-1, ub=1)

The only side effect is that the assigned names change. In order to preserve the naming convention of the Xpress library, one can specify the parameter setting name='' in the call to p.addVariables. This also makes the creation of large arrays of variables much faster. We use this shorter notation in the remainder of this chapter.

The main advantage of using NumPy operations is the ability to replicate them on each element of an array, taking into account all broadcasting features. For example, the following script ``broadcasts'' the right-hand side 1 to all elements of the array, thus creating the set of constraints x_i + y_i ≤ 1 for all i in the set S.

constr2 = x + y <= 1

All these operations can be carried out on arrays of any number of dimensions, and can be aggregated at any level. The following example shows two three-dimensional array of variables involved in two systems of constraints: the first has two variables per each of the 200 constraints, while the second has 10 constraints and 20 variables in each constraint.

z = p.addVariables(4, 5, 10)
t = p.addVariables(4, 5, 10, vartype=xp.binary)
p.addConstraint(z**2 <= 1 + t)
p.addConstraint(xp.Sum(z[i, j, k] for i in range(4) for j in range(5)) <= 4
                 for k in range(10))

Finally, a note on sums of multi-dimensional NumPy arrays: in keeping with the way NumPy arrays are handled, the sum of a multi-dimensional array results in a scalar expression with the xpress.Sum operator. The result of such a sum is exemplified by the following code:

>>> a = np.array([[1, 2, 3], [4, 5, 6]])
>>> a
array([[1, 2, 3],
       [4, 5, 6]])
>>> a.sum()
21

Products of NumPy arrays

The dot product is a useful operator for carrying out aggregate operations on vectors, matrices, and tensors. The dot operator in NumPy allows for reducing, along one axis of a multi-dimensional arrays, data such as floating points or integer values.

The application of the dot product of NumPy of two multi-dimensional arrays of dimensions (i₁, i₂, ..., i_k') and (j₁, j₂, ..., j_k''), respectively, requires that i_k' = j_k''-1, i.e., the size of the last dimension of the first array must match the size of the penultimate dimension of the second vector. For instance, the following dot product is valid:

import numpy as np
a = np.random.random((4,6))
b = np.random.random((6,2))
c = np.dot(a,b)

and the result is a 4x2 matrix. The Xpress Python interface has its own dot product operator, which can be used for all similar operations on variables and expression. The rules for applying the Xpress dot operator are the same as for the native Python dot product, with one extra feature: there is no limit on the number of arguments, hence the following example is correct as per the restrictions on the dimensions, albeit it yields a nonconvex constraint.

coeff_pre = np.random.random((6,3,7))
p = xp.problem()
x = p.addVariables(4, 7, 5)
y = p.addVariables(2, 5, 8)
coeff_post = np.random.random((6, 8, 7))
p.addConstraint(xp.Dot(coeff_pre, x, y, coeff_post) >= 0)

Similar to the NumPy dot product, the Xpress dot product has an out parameter for suppyling an array in which to store the product.

The following script defines two constraints: the first restricts the squared norm ||z|| = z · z of the vector z of variables to be at most one. It does so by applying the dot operator on the vector itself. The second constraint (t-z)'Q(t-z) ≤ 1 restricts the quadratic form on the left-hand side to be at most 1.

p.addConstraint(xp.Dot(z, z) <= 1) # restrict norm of z to 1

Q = np.random.random(N, N)    # create a random 20x20 matrix
p.addConstraint(xp.Dot((t-z), Q, (t-z)) <= 1)

As for the Sum operator, when handling variables or expressions, it is advised to use the Dot operator in the Xpress module rather than the native Python operator, for reasons of efficiency.

Products with SciPy sparse matrices

The xpress.Dot operator supports the most common SciPy sparse matrix formats, allowing arrays of sparse expressions and constraints to be constructed efficiently. For more information, see the xpress.Dot documentation.

Using Pandas with the Xpress Python interface

When building a model, it is convenient to store data in Pandas data frames. It can also be helpful to store the objects representing the decision variables, expressions involving those variables, and constraints in the data frames. The Xpress Python interface provides a special Pandas dtype, 'xpressobj' for use with series containing Xpress objects. This dtype is required to use the vector forms of the comparison (<=, ==, >=) operators and logical (&, |) operators. With the default dtype (i.e. 'object'), an expression like a <= b will result in a Boolean series, instead of a series of constraint objects.

The following code shows how Pandas could be used to store and aggregate the data, decision variables and constraint for a problem with several knapsack constraints. Note the use of the 'xpressobj' dtype, to ensure that the constraints are created correctly.

prob = xp.problem()

# Knapsack items
items = pd.DataFrame({
    'id': ['I1', 'I2', 'I3', 'I4'],
    'category': ['A', 'A', 'B', 'B'],
    'weight': [5, 12, 11, 2],
    'value': [9, 4, 10, 5],
})

# There is a different weight limit for each category of items
categories = pd.Series({
  'category': ['A', 'B'],
  'capacity': [10, 20]
})

# Create decision variables
items['pick'] = pd.Series(prob.addVariables(len(items), name='pick'), dtype='xpressobj')

# Create a constraint for each category, ensuring that the selected items in that category do not exceed its capacity
prob.addConstraint(
  (items.weight * items.pick).groupby(items.category).sum() <= categories.capacity
)

This results in the following constraints:

5 * pick(0) + 12 * pick(1) <= 10
11 * pick(2) + 2 * pick(3) <= 20