Polynomial Array

The Amplify SDK provides a polynomial array class that is a polynomial version of NumPy ndarray. Fast and efficient formulation is possible by utilizing polynomial arrays.

Polynomial Array Class

The polynomial array class is an extension of polynomial objects to multidimensional arrays. It provides a multidimensional array for each polynomial class, just as NumPy ndarray provides a multidimensional array for numeric types.

The relation between the polynomial class and the polynomial array class is as follows:

Polynomial class	Polynomial array class
`BinaryPoly`	`BinaryPolyArray`
`IsingPoly`	`IsingPolyArray`
`BinaryIntPoly`	`BinaryIntPolyArray`
`IsingIntPoly`	`IsingIntPolyArray`

Similar to NumPy ndarray, it is possible to construct a polynomial array initialized for any dimension and to convert from Python list type.

from amplify import BinaryPolyArray

# Initialize 2 rows by 3 columns binary polynomial array
poly_array = BinaryPolyArray(2, 3)  # It also can be BinaryPolyArray(shape=(2, 3))

>>> poly_array
[[0, 0, 0],
 [0, 0, 0]]

from amplify import BinaryPoly, BinaryPolyArray

# Generate a binary polynomial array from a list
poly_list = [
    [BinaryPoly({0: 1}), BinaryPoly({1: 1})],
    [BinaryPoly({2: 1}), BinaryPoly({3: 1})],
    [BinaryPoly({4: 1}), BinaryPoly({5: 1})],
]
poly_array = BinaryPolyArray(poly_list)

>>> poly_array
[[q_0, q_1],
 [q_2, q_3],
 [q_4, q_5]]

However, it is usually rare to generate a polynomial array in the manner described above. As described in Construction using variable arrays, we assume that SymbolGenerator() is used to generate variable arrays and assemble polynomial expressions.

Formulation by Polynomial Array

By utilizing index, slice, and array operations on polynomial arrays, flexible and fast formulations can be performed. Here we will explain the functions of polynomial arrays using as an example of the formulation of the objective function that appears in Formulation of the Traveling Salesman Problem.

The objective function of the Traveling Salesman Problem is as follows:

\[f = \sum_{n=0}^{N-1}{\sum_{i=0}^{N-1}{\sum_{j=0}^{N-1}{ d_{ij} q_{n, i} q_{n+1, j} }}} \quad \left( q_{N,j} = q_{0,j} \right)\]

\(N\) is the problem size (number of cities) in this case. \(q\) denotes the decision variables of the objective function which are \(N \times N\) binary variables, and \(d\) is the distance matrix representing the distance between cities.

Since the variable \(q\) has two indices, it should be represented by a two-dimensional array. Similarly, the distance matrix \(d\) should be prepared as a two-dimensional array of \(N \times N\). The distance matrix \(d\) should be determined by some input data for actual problems, but for simplicity, here we generate it by random numbers. In the following, we create a matrix of distances between 10 cities from random numbers and further a \(10 \times 10\) array of variables \(q\).

import numpy as np
from amplify import BinarySymbolGenerator

# Problem size (Number of cities)
ncity = 10

# Generate distance matrix from random numbers
d = np.triu(np.random.rand(ncity, ncity), k=1)
d += np.triu(d).T  # Convert to symmetric matrix

# Generate variable arrays
gen = BinarySymbolGenerator()
q = gen.array(ncity, ncity)  # Generate ncity x ncity variable array

The preparation is now completed. Next, several methods for constructing the objective function \(f\) will be shown. While the later methods are harder to understand, they are more efficient in terms of execution speed and simpler in terms of code expressions.

Obtain Element Using Index

Since the objective function \(f\) is written using the summation notation, among the formulation auxiliary functions, the sum_poly() function can be used. The objective function \(f\) has triply nested summation notations, therefore it is computed by multiple calls of the sum_poly() function using lambda expressions as follows:

from amplify import sum_poly

# Objective function
f = sum_poly(  # sum_n
    ncity,
    lambda n: sum_poly(  # sum_i
        ncity,
        lambda i: sum_poly(  # sum_j
            ncity, lambda j: d[i, j] * q[n, i] * q[(n + 1) % ncity, j]),  # d_{ij} q_{n, i} q_{n+1, j}
    ),
)

The first sum_poly(ncity, lambda n: corresponds to \(\sum_n\), the second sum_poly(ncity, lambda i: to \(\sum_i\), and the third sum_poly(ncity, lambda j: to \(\sum_j\). The expression \(d_{ij} q_{n,i} q_{n+1,j}\), which is the target of the sum, is expressed by d[i, j] * q[n, i] * q[(n + 1) % ncity, j].

Note, in the last lambda expression, two-dimensional indices are used to access the array elements as seen in q[n, i]. Element access to polynomial arrays is similar to NumPy ndarray, where the indices of each dimension are separated by commas and given in [].

Note

The result is the same as q[n][i], but it is not as efficient as q[n, i] because it accesses the element in two steps. For multidimensional arrays, it is recommended to specify the indices of each dimension separated by commas.

The reason for using the modulo expression (n + 1) % ncity for the index is to make the identification \(q_{N,j} = q_{0,j}\). Without using the modulo operator, the same result can be reached by rearranging the order of the summations over n. If the index \(n\) of a polynomial array is a negative value, it is interpreted as \(n + d\), where \(d\) is the length of the specified axes. Using this, the same result is obtained as d[i, j] * q[n, i] * q[n - 1, j].

Note

Please be careful when you use for loops for formulation of a function. The following in-place additive operation will produce almost the same computational efficiency as sum_poly():

f = 0
for n in range(ncity):
    for i in range(ncity):
        for j in range(ncity):
            f += d[i, j] * q[n, i] * q[n - 1, j]

On the other hand, construction using the sum() function is strongly not recommended. It is very inefficient in terms of computation time and memory usage. Computation time can take more than 100 times longer.

# Deprecated
f = sum(
    [
        d[i, j] * q[n, i] * q[n - 1, j]
        for n in range(ncity)
        for i in range(ncity)
        for j in range(ncity)
    ]
)

In the above notation using list comprehensions, replacing sum() with sum_poly() improves computation time. However, memory consumption is still high. Also, input to sum_poly() using generator expressions is currently not supported.

Array Operation Using Slices

In the index-based method, the sum was calculated by accessing each array element one by one. As a more efficient method of calculation, we consider the numerical operations on arrays, by using sub-arrays of multidimensional arrays.

Amplify SDK’s polynomial array has slicing and broadcasting similarly to NumPy ndarray. For example, to retrieve only one of the rows of the \(10 \times 10\) two-dimensional variable array q, we specify the index of the first dimension.

>>> q[0]    # Can be q[0, :] as well
[q_0, q_1, q_2, q_3, q_4, q_5, q_6, q_7, q_8, q_9]

Furthermore, to retrieve only one of the columns, we specify the index of the second dimension of the array.

>>> q[:, 0]
[ q_0, q_10, q_20, q_30, q_40, q_50, q_60, q_70, q_80, q_90]

To convert a row vector to a column vector, we give it a new axis.

>>> q[0, :, np.newaxis]
[[q_0],
 [q_1],
 [q_2],
 [q_3],
 [q_4],
 [q_5],
 [q_6],
 [q_7],
 [q_8],
 [q_9]]

Note

A sub-array extracted by slicing returns a view, not a copy of the original array. Setting a value to the view will also be reflected in the original array. Nested views can also be created, but for efficiency reasons, the number of views is currently limited to 3 levels.

To write down the objective function \(f\) as operations on arrays, we first fix the value of \(n\).

\[f_n = {\sum_{i=0}^{N-1}{\sum_{j=0}^{N-1}{d_{ij} q_{n, i} q_{n+1, j} }}}\]

The above expression corresponds to the operation of computing the direct product of \(q_n\) and \(q_{n+1}\), take the Hadamard product (product of the elements) of the direct product and \(d\), and finally summing up the matrix elements.

\[f_n = \sum_{ \mathrm{all \ elements} } \ {d \circ \left( q_{n} \otimes q_{n+1} \right) }\]

This can be expressed in terms of polynomial array operations. First, \(q_{n} \otimes q_{n+1}\) is the product of column and row vectors, that is, an \(N \times N\) matrix. For example, for the case n = 0 it is computed as follows:

from amplify import newaxis

q1 = q[1, :]    # One-dimensional array of length ncity
q0t = q[0, :, newaxis]  # Two-dimensional array of ncity x 1
f0 = (d * (q0t * q1)).sum()

For the operations between arrays of unmatched shape, the dimension is stretched by broadcasting. First, q0t is a two-dimensional array of ncity x 1 and q1 is a one-dimensional array of length ncity, so q1 is interpreted as a two-dimensional array of 1 x ncity. Next, the dimensions of q0t, q1 are compared, and attempts are made to complement dimensions that do not match in length by iterating in the direction of that dimension. That is, q0t is repeated ncity times in the row direction to bring the number of columns from 1 to ncity, and q1 is repeated ncity times in the column direction to bring the number of rows from 1 to ncity. They are then interpreted as ncity x ncity matrices, and the product of the elements can be computed. Then the product of the elements with the distance matrix d is computed and the sum of all elements is computed with the sum() method.

For convenience of explanation, \(q_{n} \otimes q_{n+1}\) was calculated first, but since the broadcast function also applies to the product of a matrix and a vector, there is no need to distinguish between the direct product and the Hadamard product, and the associativity holds. Therefore, the same result is obtained as follows:

f0 = (d * q0t * q1).sum()

With the above, the objective function \(f\) is formulated as follows by computing the sum over all n.

from amplify import sum_poly, newaxis

f = sum_poly(ncity, lambda n: (d * q[n, :, newaxis] * q[n - 1]).sum())

For simplicity, the order of the sums is rearranged and written without modulo operator. Using array operations in this way simplifies the formulation. It is also several times faster than nested sum_poly() functions or for loops.

Einstein Summation Convention

Lastly, this part explains how to write a formulation entirely in terms of operations on arrays. The Amplify SDK provides a function equivalent to NumPy’s einsum() for polynomial arrays. This function is called Einstein summation convention. By specifying the target arrays and a string representing their subscripts, this function returns a scalar or array of polynomials, resulting of the sum of the products of their elements.

For example, let’s consider the following case of computing a four-dimensional array \(G\) where the array elements \(g_{p,q,r,s}\) are defined from a two-dimensional array \(A,B,C,D\).

\[g_{p,q,r,s} = \sum_{i, j, k, l} A_{ip} B_{jq} C_{kr} D_{ls}\]

To generate a polynomial array, call the einsum() function or sum_poly() (alias) as follows:

from amplify import einsum

# Create a four-dimensional array G with A, B, C, and D as two-dimensional arrays
G = einsum("ip,jq,kr,ls->pqrs", A, B, C, D) # Same if sum_poly is used instead of einsum

In the contraction rule ip,jq,kr,ls->pqrs, the left-hand side of -> describes comma-separated subscripts of the input arrays, and the right-hand side represents the subscript of the output array. The characters appearing only on the left-hand side are summed, and the array with the dimension represented by the character appearing on the right-hand side is output. If no subscripts are left on the right-hand side, a scalar is output instead of an array. Please refer to Formulation Auxiliary Functions for examples.

There is only one problem in computing the objective function \(f\) of the Traveling Salesman Problem using this Einstein summation convention. Since no expression can be given for the subscripts, n + 1 in \(q_{n+1,j}\) should be adjusted. So we create a new variable array qr so that \(q_{n, j}' = q_{n + 1, j}\) and n can be used instead of n + 1 for the subscripts.

\[f = \sum_{n=0}^{N-1}{\sum_{i=0}^{N-1}{\sum_{j=0}^{N-1}{ d_{ij} q_{n, i} q_{n, j}' }}}\]

To create a variable array qr, the roll() method can be used to rotate the elements along a specified axis. In summary, the formulation of the objective function \(f\) is described in the following very simple form:

from amplify import sum_poly

qr = q.roll(-1, axis=0) # Displace q by -1 in the column direction (axis=0)
f = sum_poly("ij,ni,nj->", d, q, qr)    # einsum can be used instead of sum_poly

In terms of calculation speed, it is the fastest of the methods introduced so far. Compared to the first nested sum_poly() function and for loops, a speedup of up to 10 times may be expected.

NumPy Compatible Methods and Functions

The correspondence between the methods and functions of polynomial arrays and NumPy arrays is as follows. Restrictions and specifications for some functions may vary.

Attributes

Polynomial Array	NumPy	Description
`data`	`numpy.ndarray.data`	Read only
`size`	`numpy.ndarray.size`	Read only
`itemsize`	`numpy.ndarray.itemsize`	Read only
`ndim`	`numpy.ndarray.ndim`	Read only
`strides`	`numpy.ndarray.strides`	Read only
`shape`	`numpy.ndarray.shape`	Cannot be set in view
`T`	`numpy.ndarray.T`	Read only, Return copy if view cannot be created

Methods

Polynomial Array	NumPy	Description
`transpose()`	`numpy.ndarray.transpose()`	Return copy if view cannot be created
`view()`	`numpy.ndarray.view()`
`copy()`	`numpy.ndarray.copy()`
`diagonal()`	`numpy.ndarray.diagonal()`	Return copy
`fill()`	`numpy.ndarray.fill()`
`flatten()`	`numpy.ndarray.flatten()`
`ravel()`	`numpy.ndarray.ravel()`	Return copy
`reshape()`	`numpy.ndarray.reshape()`	Return copy
`roll()`	`numpy.roll()`	Provided as a method, not a function Return copy if view cannot be created
`sum()`	`numpy.ndarray.sum()`
`tolist()`	`numpy.ndarray.tolist()`

Special Methods

Polynomial Array (w/o Logic)	NumPy	Description
`__len__()`	`numpy.ndarray.__len__()`
`__getitem__()`	`numpy.ndarray.__getitem__()`	Not supported for Advanced Indexing
`__setitem__()`	`numpy.ndarray.__setitem__()`	Not supported for Advanced Indexing
`__add__()`	`numpy.ndarray.__add__()`
`__radd__()`
`__iadd__()`	`numpy.ndarray.__iadd__()`
`__sub__()`	`numpy.ndarray.__sub__()`
`__rsub__()`
`__isub__()`	`numpy.ndarray.__isub__()`
`__mul__()`	`numpy.ndarray.__mul__()`
`__rmul__()`
`__imul__()`	`numpy.ndarray.__imul__()`
`__truediv__()`	`numpy.ndarray.__truediv__()`	When performing division of `float` type on an integer coefficient polynomial, the coefficient type is not changed from integer.
`__itruediv__()`	`numpy.ndarray.__itruediv__()`	When performing division of `float` type on an integer coefficient polynomial, the coefficient type is not changed from integer.
`__floordiv__()`	`numpy.ndarray.__floordiv__()`	Provided only for integer coefficient polynomials
`__ifloordiv__()`	`numpy.ndarray.__ifloordiv__()`	Provided only for integer coefficient polynomials
`__pow__()`	`numpy.ndarray.__pow__()`	Only for integers greater than or equal to 0
`__ipow__()`	`numpy.ndarray.__ipow__()`	Only for integers greater than or equal to 0
`__and__()`	`numpy.ndarray.__and__()`	Provided only for binary variables
`__rand__()`		Provided only for binary variables
`__iand__()`	`numpy.ndarray.__iand__()`	Provided only for binary variables
`__or__()`	`numpy.ndarray.__or__()`	Provided only for binary variables
`__ror__()`		Provided only for binary variables
`__ior__()`	`numpy.ndarray.__ior__()`	Provided only for binary variables
`__xor__()`	`numpy.ndarray.__xor__()`	Provided only for binary variables
`__rxor__()`		Provided only for binary variables
`__ixor__()`	`numpy.ndarray.__ixor__()`	Provided only for binary variables
`__eq__()`	`numpy.ndarray.__eq__()`
`__ne__()`	`numpy.ndarray.__ne__()`
`__invert__()`	`numpy.ndarray.__invert__()`	Provided only in binary variables
`__pos__()`	`numpy.ndarray.__pos__()`
`__neg__()`	`numpy.ndarray.__neg__()`
`__int__()`	`numpy.ndarray.__int__()`	Only if the number of elements is 1 and polynomials can be converted to numbers
`__float__()`	`numpy.ndarray.__float__()`	Only if the number of elements is 1 and polynomials can be converted to numbers
`__iter__()`
`__repr__()`	`numpy.ndarray.__repr__()`
`__str__()`	`numpy.ndarray.__str__()`

Functions

Polynomial Array	NumPy	Description
`einsum()`	`numpy.einsum()`	ellipsis (`...`) is not yet supported.

Polynomial Array Methods and Functions

Attribute	Description
`isview`	Read only, Return `True` if the array is a view, `False` otherwise

Method

Description

numpy()

Convert to numpy.ndarray

Only if polynomials of all array elements can be converted to numbers

decode()

Return an array with the variables of the polynomial replaced by the given solution [reference]