Polynomial Array

Amplify SDK では NumPy ndarray の多項式版とも言える多項式配列クラスが提供されます。多項式配列を活用することで、高速かつ効率的な定式化が可能です。

参考

このページを読む前に Polynomial にて多項式オブジェクトの使用方法を確認してください。

多項式配列クラス

多項式配列クラスは多項式オブジェクトの多次元配列への拡張です。 NumPy ndarray が数値型に対する多次元配列を提供するのと同じように、それぞれの多項式クラスに対する多次元配列を提供します。

多項式クラスと多項式配列クラスの関係は次の通りです。

多項式クラス	多項式配列クラス
`BinaryPoly`	`BinaryPolyArray`
`IsingPoly`	`IsingPolyArray`
`BinaryIntPoly`	`BinaryIntPolyArray`
`IsingIntPoly`	`IsingIntPolyArray`

NumPy ndarray と同様に、任意の次元に対して初期化された多項式配列を構築したり、Python list 型からの変換が可能です。

from amplify import BinaryPolyArray

# 2行3列のバイナリ多項式配列を初期化
poly_array = BinaryPolyArray(2, 3)  # BinaryPolyArray(shape=(2, 3)) としても同様

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

from amplify import BinaryPoly, BinaryPolyArray

# リストからバイナリ多項式配列を構築
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]]

しかし、上記の方法で多項式配列を構築することは通常は稀でしょう。変数配列を用いた構築で説明したように、SymbolGenerator() を用いて変数配列を構築し多項式を組み立てることを想定しています。

多項式配列による定式化

多項式配列に対するインデックス、スライス、配列演算を活用することで、柔軟で高速な定式化が行えます。ここでは、巡回セールスマン問題の定式化に現れる目的関数の数式を例に、多項式配列の持つ機能を説明します。

巡回セールスマン問題の目的関数は次の通りです。

\[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\) は問題サイズ (都市数) を表します。\(q\) は目的関数の決定変数でバイナリ変数とします。また、\(d\) は都市間の距離を表す距離行列とします。

変数 \(q\) は二つのインデックスを持つため、 \(N \times N\) の二次元配列で表せば良さそうです。同様に距離行列 \(d\) も \(N \times N\) の二次元配列として用意しておきます。距離行列 \(d\) は何らかの入力データから与えられるものですが、簡単のため乱数で生成しておきます。以下では 10 都市間の距離行列を乱数から作成し、 \(10 \times 10\) の変数配列 \(q\) を作成します。

import numpy as np
from amplify import BinarySymbolGenerator

# 問題サイズ (都市数)
ncity = 10

# 距離行列を乱数から作成
d = np.triu(np.random.rand(ncity, ncity), k=1)
d += np.triu(d).T  # 対称行列に変換

# 変数配列の作成
gen = BinarySymbolGenerator()
q = gen.array(ncity, ncity)  # ncity x ncity の変数配列を生成

以上で準備は完了です。ここから、目的関数 \(f\) の構築についていくつかの方法を紹介します。後ろの方法ほど難しくなる一方、実行速度の観点から効率的でコードの記述もシンプルになります。

インデックスを用いた要素取得

目的関数 \(f\) は和の記号を用いて書かれています。定式化補助関数のうち、和の記号に相当する sum_poly() 関数を用いることが出来ます。ただし、目的関数 \(f\) は和の記号が三重にネストしています。そのためラムダ式を用いて次のようにして sum_poly() 関数を多重に呼び出すことで目的関数が計算されます。

from amplify import sum_poly

# 目的関数
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}
    ),
)

一つ目の sum_poly(ncity, lambda n: は \(\sum_n\) に、二つ目の sum_poly(ncity, lambda i: は \(\sum_i\) に、三つ目の sum_poly(ncity, lambda j: は \(\sum_j\) に対応します。和の対象となる式 \(d_{ij} q_{n,i} q_{n+1,j}\) は d[i, j] * q[n, i] * q[(n + 1) % ncity, j] として表されます。

最後のラムダ式では、配列要素へのアクセスに q[n, i] と二次元のインデックスを指定しています。多項式配列の要素アクセスは、 NumPy ndarray と同様に、各次元のインデックスをカンマで区切って [] に与えます。

注釈

q[n][i] としても結果は同じになりますが、二段階で要素を抽出するため q[n, i] と比べて非効率です。多次元配列では各次元のインデックスをカンマで区切って指定することを推奨します。

インデックスに (n + 1) % ncity と剰余を用いているのは \(q_{N,j} = q_{0,j}\) とするためです。剰余を用いない書き方として n に対する和の順番を入れ替える方法があります。多項式配列のインデックス \(n\) に負の値が与えられた場合、指定された次元の長さを \(d\) として、\(n + d\) と解釈されます。これを用いると d[i, j] * q[n, i] * q[n - 1, j] としても同様の結果になります。

注釈

for ループを用いて構築する場合は注意が必要です。次のようにインプレースの加法演算を用いて記述すると 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]

一方、sum() 関数を用いた構築は 強く非推奨 です。計算時間やメモリ使用量の観点から非常に非効率的です。計算時間は100倍以上かかってしまう場合があります。

# 非推奨
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)
    ]
)

上記のリスト内包表記を用いる書き方において、 sum() を sum_poly() で置き換えると計算時間が改善します。しかしメモリ消費は依然大きいです。また、ジェネレータ式を用いた sum_poly() への入力は現在のところ非対応です。

スライスを用いた配列演算

インデックスを用いた方法では配列要素一つ一つにアクセスして和を計算しました。さらに効率的に計算するための方法として、多次元配列の一部を部分配列として取り出すことで、配列同士の数値演算で計算できないかを考えます。

Amplify SDK の多項式配列は、NumPy ndarray と同様にスライスやブロードキャストが備わっています。例えば、先ほど作成した \(10 \times 10\) の二次元変数配列 q の行のみを取り出すには、一次元目のインデックスを指定します。

>>> q[0]    # q[0, :] としても同様
[q_0, q_1, q_2, q_3, q_4, q_5, q_6, q_7, q_8, q_9]

また、列のみを取り出すには、配列の二次元目のインデックスを指定します。

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

行ベクトルを列ベクトルに変換するには新しい軸を与えます。

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

注釈

スライスで取り出した部分配列は元の配列のコピーではなくビューを返します。ビューに値をセットすると元の配列にも反映されます。ネストしたビューを作成することも可能ですが、効率の問題から現在3段階までのビューの作成に制限しています。

目的関数 \(f\) を配列同士の演算と解釈するために、まず \(n\) については固定して考えます。

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

上記は、\(q_n\) と \(q_{n+1}\) の直積と \(d\) とのアダマール積 (要素同士の積) を計算し、最後に行列要素の総和を取る演算に一致します。

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

これを多項式配列の演算で表現します。まず \(q_{n} \otimes q_{n+1}\) は列ベクトルと行ベクトルの積です。例えば n = 0 の場合について次のようにして計算されます。

from amplify import newaxis

q1 = q[1, :]    # 長さ ncity の一次元配列
q0t = q[0, :, newaxis]  # ncity x 1 の二次元配列
f0 = (d * (q0t * q1)).sum()

形状の一致しない配列同士の演算ではブロードキャストにより次元が伸張されます。まず q0t は ncity x 1 の二次元配列、q1 は長さ ncity の一次元配列なので、q1 は 1 x ncity の二次元配列と解釈されます。次に q0t, q1 の各次元を比べ、長さが一致しない次元についてはその次元の方向への反復により補完を試みます。つまり、q0t については列数を 1 から ncity にするために行方向に ncity 回繰り返され、q1 については行数を 1 から ncity にするために列方向に ncity 回繰り返され、双方を ncity x ncity の行列と解釈し要素同士の積が計算されます。その後、距離行列 d との要素積を計算し、sum() メソッドで全ての要素の和を計算します。

説明の便宜上 \(q_{n} \otimes q_{n+1}\) を先に計算しましたが、ブロードキャスト機能は行列とベクトルの積にも適用されるので、直積とアダマール積は区別の必要が無く、結合法則が成立します。そのため次のように書いても結果は同様です。

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

以上により、全ての n について和を計算することで、次のように目的関数 \(f\) が定式化されます。

from amplify import sum_poly, newaxis

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

ここで簡便のために和の順番を入れ替えて剰余を用いない書き方としました。このように配列演算を用いると、定式化をシンプルに記述することが出来ます。また、ネストした sum_poly() 関数や for ループと比べ計算速度において数倍の高速化が見込めます。

アインシュタインの縮約記法

最後に、完全に配列同士の演算だけで定式化を記述する方法を紹介します。Amplify SDK では多項式配列に対する NumPy の einsum() に相当する関数が提供されます。この関数はアインシュタインの縮約記法と呼ばれ、配列の添え字を表す文字列と要素積の和を計算する配列を与えることで、多項式のスカラーまたは配列が返却されます。

例えば、次のように二次元配列 \(A,B,C,D\) から配列要素 \(g_{p,q,r,s}\) が定義される四次元配列 \(G\) を計算する場合を考えます。

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

これは次のように einsum() 関数または sum_poly() (エイリアス) を呼びだすことで多項式配列が生成されます。

from amplify import einsum

# A, B, C, D を二次元配列として四次元配列 G を作成
G = einsum("ip,jq,kr,ls->pqrs", A, B, C, D) # einsum の代わりに sum_poly を用いても同様

縮約規則 ip,jq,kr,ls->pqrs は -> の左側にカンマ区切りで各入力配列の添え字を表し、右側に出力配列に残す添え字を記述します。左側にのみ現れる文字は和を計算し、右側に現れる文字で表された次元を持つ配列が出力されます。もし右側に添え字を残さなければ、配列では無くスカラーを出力します。実行例については定式化補助関数を参照してください。

このアインシュタインの縮約記法を用いて、巡回セールスマン問題の目的関数 \(f\) を計算する上で一つだけ問題があります。添え字に関しては式を与えられないため、\(q_{n+1,j}\) の n + 1 についての工夫が必要です。そこで、新たに \(q_{n, j}' = q_{n + 1, j}\) となるような変数配列 qr を作成し、添え字を n + 1 ではなく n で利用出来るようにします。

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

変数配列 qr を作成するには、指定された軸について要素を回転させる roll() メソッドが利用できます。以上をまとめると、次のような非常にシンプルな形で目的関数 \(f\) の定式化が記述されます。

from amplify import sum_poly

qr = q.roll(-1, axis=0) # q を列方向 (axis=0) に -1 ずらす
f = sum_poly("ij,ni,nj->", d, q, qr)    # sum_poly の代わりに einsum を用いても同様

計算速度についても、これまで紹介した方法の中で最も高速です。最初のネストした sum_poly() 関数や for ループと比べると、10倍程度の高速化が見込める場合があります。

参考

巡回セールスマン問題の定式化の続きは、制約条件の定式化を参照してください。

NumPy 互換メソッドと関数

多項式配列のメソッドや関数と NumPy 配列の対応関係は次の通りです。一部の機能について制限や仕様が異なる場合があります。

アトリビュート

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`	ビューにはセット出来ない
`T`	`numpy.ndarray.T`	Read only, ビューを作成できない場合コピーを返却

メソッド

Polynomial Array	NumPy	Description
`transpose()`	`numpy.ndarray.transpose()`	ビューを作成できない場合コピーを返却
`view()`	`numpy.ndarray.view()`
`copy()`	`numpy.ndarray.copy()`
`diagonal()`	`numpy.ndarray.diagonal()`	コピーを返却
`fill()`	`numpy.ndarray.fill()`
`flatten()`	`numpy.ndarray.flatten()`
`ravel()`	`numpy.ndarray.ravel()`	コピーを返却
`reshape()`	`numpy.ndarray.reshape()`	コピーを返却
`roll()`	`numpy.roll()`	関数ではなくメソッドで提供ビューを作成できない場合コピーを返却
`sum()`	`numpy.ndarray.sum()`
`tolist()`	`numpy.ndarray.tolist()`

特殊メソッド

Polynomial Array (w/o Logic)	NumPy	Description
`__len__()`	`numpy.ndarray.__len__()`
`__getitem__()`	`numpy.ndarray.__getitem__()`	Advanced Indexing には未対応
`__setitem__()`	`numpy.ndarray.__setitem__()`	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__()`	整数係数多項式に対して `float` 型の除算を行った場合、係数型は整数型から変更されない
`__itruediv__()`	`numpy.ndarray.__itruediv__()`	整数係数多項式に対して `float` 型の除算を行った場合、係数型は整数型から変更されない
`__floordiv__()`	`numpy.ndarray.__floordiv__()`	整数係数多項式に対してのみ提供
`__ifloordiv__()`	`numpy.ndarray.__ifloordiv__()`	整数係数多項式に対してのみ提供
`__pow__()`	`numpy.ndarray.__pow__()`	0以上の整数のみ
`__ipow__()`	`numpy.ndarray.__ipow__()`	0以上の整数のみ
`__and__()`	`numpy.ndarray.__and__()`	バリナリ変数でのみ提供
`__rand__()`		バリナリ変数でのみ提供
`__iand__()`	`numpy.ndarray.__iand__()`	バリナリ変数でのみ提供
`__or__()`	`numpy.ndarray.__or__()`	バリナリ変数でのみ提供
`__ror__()`		バリナリ変数でのみ提供
`__ior__()`	`numpy.ndarray.__ior__()`	バリナリ変数でのみ提供
`__xor__()`	`numpy.ndarray.__xor__()`	バリナリ変数でのみ提供
`__rxor__()`		バリナリ変数でのみ提供
`__ixor__()`	`numpy.ndarray.__ixor__()`	バリナリ変数でのみ提供
`__eq__()`	`numpy.ndarray.__eq__()`
`__ne__()`	`numpy.ndarray.__ne__()`
`__invert__()`	`numpy.ndarray.__invert__()`	バリナリ変数でのみ提供
`__pos__()`	`numpy.ndarray.__pos__()`
`__neg__()`	`numpy.ndarray.__neg__()`
`__int__()`	`numpy.ndarray.__int__()`	要素数が1で多項式を数値に変換できる場合のみ
`__float__()`	`numpy.ndarray.__float__()`	要素数が1で多項式を数値に変換できる場合のみ
`__iter__()`
`__repr__()`	`numpy.ndarray.__repr__()`
`__str__()`	`numpy.ndarray.__str__()`

関数

Polynomial Array	NumPy	Description
`einsum()`	`numpy.einsum()`	エリプシス (`...`) には未対応

多項式配列のメソッドと関数

Attribute	Description
`isview`	Read only, 配列がビューなら `True` そうで無いなら `False` を返す

Method

Description

numpy()

numpy.ndarray に変換する

すべての配列要素の多項式が数値に変換できる場合のみ

decode()

多項式の変数を与えられた解で置き換えた配列を返却する [参照]