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 aims to support the oldest version of NumPy available as a prebuilt binary package, which varies depending on the platform and Python version.

Using NumPy in 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 = np.array([p.addVariable()                  for i in S], dtype=xp.npvar)
y = np.array([p.addVariable(vartype=xp.binary) for i in S], dtype=xp.npvar)
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 NumPy arrays must have the attribute dtype equal to xpress.npvar (abbreviated to xp.npvar here) in order to use the matricial/vectorial form of the comparison (<=, =, >=), arithmetic (+, -, *, /, **), and logic (&, |) operators.

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 a NumPy array 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 = np.array([p.addVariable(name='v({0})'.format(i)) for i in range(20)],
      dtype=xp.npvar).reshape(5,4)
y = np.array([p.addVariable(vartype=xp.binary) for i in range(27)],
      dtype=xp.npvar).reshape(3,3,3)
z = np.array([p.addVariable(lb=-1, ub=1) for i in range(1000)], dtype=xp.npvar)

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 defining the output 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:

x = p.addVariables(1000)
A = scipy.sparse.csr_array(scipy.sparse.random(1000, 1000, density=0.1))
b = np.random.random(1000)
p.addConstraint(xp.Dot(A, x) <= b)

xpress.Dot can compute the product of a 1-dimensional NumPy array of variables or expressions with a sparse matrix of type int32, int64, float32 or float64 in CSR or CSC format. Computing the product of an N-dimensional array with a sparse matrix is not currently supported when N > 1.

See https://docs.scipy.org/doc/scipy/reference/sparse.html for more information about SciPy sparse matrices.