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.
See also
Before reading this page, please check Polynomial to see how to use polynomial objects.
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 

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:
\(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 twodimensional array. Similarly, the distance matrix \(d\) should be prepared as a twodimensional 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, twodimensional 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 inplace 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 indexbased 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 subarrays 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\) twodimensional 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 subarray 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\).
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.
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, :] # Onedimensional array of length ncity
q0t = q[0, :, newaxis] # Twodimensional 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 twodimensional array of ncity x 1
and q1
is a onedimensional array of length ncity
, so q1
is interpreted as a twodimensional 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 fourdimensional array \(G\) where the array elements \(g_{p,q,r,s}\) are defined from a twodimensional array \(A,B,C,D\).
To generate a polynomial array, call the einsum()
function or sum_poly()
(alias) as follows:
from amplify import einsum
# Create a fourdimensional array G with A, B, C, and D as twodimensional 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 lefthand side of >
describes commaseparated subscripts of the input arrays, and the righthand side represents the subscript of the output array. The characters appearing only on the lefthand side are summed, and the array with the dimension represented by the character appearing on the righthand side is output. If no subscripts are left on the righthand 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.
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.
See also
For the rest of the formulation of the Traveling Salesman Problem, please refer to Formulation of Constraints.
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 


Read only 


Read only 


Read only 


Read only 


Read only 


Cannot be set in view 


Read only, Return copy if view cannot be created 
Methods
Polynomial Array 
NumPy 
Description 

Return copy if view cannot be created 

Return copy 

Return copy 

Return copy 

Provided as a method, not a function Return copy if view cannot be created 

Special Methods
Polynomial Array (w/o Logic) 
NumPy 
Description 




Not supported for Advanced Indexing 


Not supported for Advanced Indexing 




















When performing division of 


When performing division of 


Provided only for integer coefficient polynomials 


Provided only for integer coefficient polynomials 


Only for integers greater than or equal to 0 


Only for integers greater than or equal to 0 


Provided only for binary variables 


Provided only for binary variables 


Provided only for binary variables 


Provided only for binary variables 


Provided only for binary variables 


Provided only for binary variables 


Provided only for binary variables 


Provided only for binary variables 


Provided only for binary variables 






Provided only in binary variables 






Only if the number of elements is 1 and polynomials can be converted to numbers 


Only if the number of elements is 1 and polynomials can be converted to numbers 






Functions
Polynomial Array 
NumPy 
Description 

ellipsis ( 
Polynomial Array Methods and Functions
Attribute 
Description 


Read only, Return 
Method 
Description 

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

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