グリッドの x 軸と y 軸を定義する 2 つの numpy 配列があります。例えば：
x = numpy.array([1,2,3])
y = numpy.array([4,5])
これらの配列のデカルト積を生成して生成したいと思います。
array([[1,4],[2,4],[3,4],[1,5],[2,5],[3,5]])
ループでこれを何度も行う必要があるため、それほど非効率的ではありません。それらをPythonリストに変換しitertools.product、numpy配列を使用して戻すことは、最も効率的な形式ではないと想定しています。
cartesian_product（ほぼ）この問題には、さまざまなプロパティを使用して多くのアプローチがあります。他のものより高速なものもあれば、より汎用的なものもあります。多くのテストと調整を行った結果、n次元を計算する次の関数は、cartesian_product多くの入力で他のほとんどの関数よりも高速であることがわかりました。少し複雑ですが、多くの場合は少し高速な2つのアプローチについては、PaulPanzerによる回答を参照してください。
その答えを考えると、これは私が知っているデカルト積の最速の実装ではありません。numpyただし、その単純さは、将来の改善のための有用なベンチマークになり続けると思います。
def cartesian_product(*arrays):
la = len(arrays)
dtype = numpy.result_type(*arrays)
arr = numpy.empty([len(a) for a in arrays] + [la], dtype=dtype)
for i, a in enumerate(numpy.ix_(*arrays)):
arr[...,i] = a
return arr.reshape(-1, la)
この関数がix_通常とは異なる方法で使用することは言及する価値があります。文書化された使用法は配列にインデックスix_を生成することですが、同じ形状の配列をブロードキャスト割り当てに使用できるのは偶然です。この方法を試してみるように促してくれたmgilsonと、この回答について、使用の提案など、非常に役立つフィードバックを提供してくれたunutbuに感謝します。 ix_numpy.result_type
連続したメモリブロックをFortranの順序で書き込む方が速い場合があります。これがこの代替案の基礎であり、一部のハードウェアでは（以下を参照）cartesian_product_transposeよりも高速であることが証明されています。cartesian_productただし、同じ原理を使用するPaul Panzerの回答は、さらに高速です。それでも、興味のある読者のためにこれをここに含めます：
def cartesian_product_transpose(*arrays):
broadcastable = numpy.ix_(*arrays)
broadcasted = numpy.broadcast_arrays(*broadcastable)
rows, cols = numpy.prod(broadcasted[0].shape), len(broadcasted)
dtype = numpy.result_type(*arrays)
out = numpy.empty(rows * cols, dtype=dtype)
start, end = 0, rows
for a in broadcasted:
out[start:end] = a.reshape(-1)
start, end = end, end + rows
return out.reshape(cols, rows).T
パンツァーのアプローチを理解するようになった後、私は彼とほぼ同じくらい速く、ほぼ同じくらい簡単な新しいバージョンを書きましたcartesian_product：
def cartesian_product_simple_transpose(arrays):
la = len(arrays)
dtype = numpy.result_type(*arrays)
arr = numpy.empty([la] + [len(a) for a in arrays], dtype=dtype)
for i, a in enumerate(numpy.ix_(*arrays)):
arr[i, ...] = a
return arr.reshape(la, -1).T
これには一定時間のオーバーヘッドがあり、小さな入力の場合はPanzerよりも実行が遅くなるようです。しかし、より大きな入力の場合、私が実行したすべてのテストで、彼の最速の実装と同じように実行されます（cartesian_product_transpose_pp）。
次のセクションでは、他の選択肢のいくつかのテストを含めます。これらは今ややや時代遅れですが、重複した努力ではなく、歴史的な興味からここに残すことにしました。最新のテストについては、Panzerの回答とNicoSchlömerの回答を参照してください。
これは、これらの機能のいくつかがいくつかの選択肢と比較して提供するパフォーマンスの向上を示す一連のテストです。numpyここに示すすべてのテストは、Mac OS 10.12.5、Python 3.6.1、および1.12.1を実行しているクアッドコアマシンで実行されました。ハードウェアとソフトウェアのバリエーションは異なる結果を生み出すことが知られているので、YMMV。これらのテストを自分で実行して確認してください。
定義：
import numpy
import itertools
from functools import reduce
### Two-dimensional products ###
def repeat_product(x, y):
return numpy.transpose([numpy.tile(x, len(y)),
numpy.repeat(y, len(x))])
def dstack_product(x, y):
return numpy.dstack(numpy.meshgrid(x, y)).reshape(-1, 2)
### Generalized N-dimensional products ###
def cartesian_product(*arrays):
la = len(arrays)
dtype = numpy.result_type(*arrays)
arr = numpy.empty([len(a) for a in arrays] + [la], dtype=dtype)
for i, a in enumerate(numpy.ix_(*arrays)):
arr[...,i] = a
return arr.reshape(-1, la)
def cartesian_product_transpose(*arrays):
broadcastable = numpy.ix_(*arrays)
broadcasted = numpy.broadcast_arrays(*broadcastable)
rows, cols = numpy.prod(broadcasted[0].shape), len(broadcasted)
dtype = numpy.result_type(*arrays)
out = numpy.empty(rows * cols, dtype=dtype)
start, end = 0, rows
for a in broadcasted:
out[start:end] = a.reshape(-1)
start, end = end, end + rows
return out.reshape(cols, rows).T
# from https://stackoverflow.com/a/1235363/577088
def cartesian_product_recursive(*arrays, out=None):
arrays = [numpy.asarray(x) for x in arrays]
dtype = arrays[0].dtype
n = numpy.prod([x.size for x in arrays])
if out is None:
out = numpy.zeros([n, len(arrays)], dtype=dtype)
m = n // arrays[0].size
out[:,0] = numpy.repeat(arrays[0], m)
if arrays[1:]:
cartesian_product_recursive(arrays[1:], out=out[0:m,1:])
for j in range(1, arrays[0].size):
out[j*m:(j+1)*m,1:] = out[0:m,1:]
return out
def cartesian_product_itertools(*arrays):
return numpy.array(list(itertools.product(*arrays)))
### Test code ###
name_func = [('repeat_product',
repeat_product),
('dstack_product',
dstack_product),
('cartesian_product',
cartesian_product),
('cartesian_product_transpose',
cartesian_product_transpose),
('cartesian_product_recursive',
cartesian_product_recursive),
('cartesian_product_itertools',
cartesian_product_itertools)]
def test(in_arrays, test_funcs):
global func
global arrays
arrays = in_arrays
for name, func in test_funcs:
print('{}:'.format(name))
%timeit func(*arrays)
def test_all(*in_arrays):
test(in_arrays, name_func)
# `cartesian_product_recursive` throws an
# unexpected error when used on more than
# two input arrays, so for now I've removed
# it from these tests.
def test_cartesian(*in_arrays):
test(in_arrays, name_func[2:4] + name_func[-1:])
x10 = [numpy.arange(10)]
x50 = [numpy.arange(50)]
x100 = [numpy.arange(100)]
x500 = [numpy.arange(500)]
x1000 = [numpy.arange(1000)]
試験結果：
In [2]: test_all(*(x100 * 2))
repeat_product:
67.5 µs ± 633 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
dstack_product:
67.7 µs ± 1.09 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
cartesian_product:
33.4 µs ± 558 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
cartesian_product_transpose:
67.7 µs ± 932 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
cartesian_product_recursive:
215 µs ± 6.01 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
cartesian_product_itertools:
3.65 ms ± 38.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
In [3]: test_all(*(x500 * 2))
repeat_product:
1.31 ms ± 9.28 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
dstack_product:
1.27 ms ± 7.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
cartesian_product:
375 µs ± 4.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
cartesian_product_transpose:
488 µs ± 8.88 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
cartesian_product_recursive:
2.21 ms ± 38.4 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
cartesian_product_itertools:
105 ms ± 1.17 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
In [4]: test_all(*(x1000 * 2))
repeat_product:
10.2 ms ± 132 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
dstack_product:
12 ms ± 120 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
cartesian_product:
4.75 ms ± 57.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
cartesian_product_transpose:
7.76 ms ± 52.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
cartesian_product_recursive:
13 ms ± 209 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
cartesian_product_itertools:
422 ms ± 7.77 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
すべての場合において、cartesian_productこの回答の冒頭で定義されているように、最速です。
任意の数の入力配列を受け入れる関数の場合len(arrays) > 2、パフォーマンスもチェックする価値があります。cartesian_product_recursive（この場合にエラーがスローされる理由を特定できるまで、これらのテストからエラーを削除しました。）
In [5]: test_cartesian(*(x100 * 3))
cartesian_product:
8.8 ms ± 138 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
cartesian_product_transpose:
7.87 ms ± 91.5 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
cartesian_product_itertools:
518 ms ± 5.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
In [6]: test_cartesian(*(x50 * 4))
cartesian_product:
169 ms ± 5.1 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
cartesian_product_transpose:
184 ms ± 4.32 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
cartesian_product_itertools:
3.69 s ± 73.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
In [7]: test_cartesian(*(x10 * 6))
cartesian_product:
26.5 ms ± 449 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
cartesian_product_transpose:
16 ms ± 133 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
cartesian_product_itertools:
728 ms ± 16 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
In [8]: test_cartesian(*(x10 * 7))
cartesian_product:
650 ms ± 8.14 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
cartesian_product_transpose:
518 ms ± 7.09 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
cartesian_product_itertools:
8.13 s ± 122 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
これらのテストが示すようにcartesian_product、入力配列の数が（おおよそ）4を超えるまで競争力を維持します。その後cartesian_product_transpose、わずかなエッジがあります。
他のハードウェアおよびオペレーティングシステムを使用しているユーザーには、異なる結果が表示される可能性があることを繰り返す価値があります。たとえば、unutbuは、Ubuntu 14.04、Python 3.4.3、およびnumpy1.14.0.dev0+b7050a9を使用したこれらのテストで次の結果が表示されることを報告しています。
>>> %timeit cartesian_product_transpose(x500, y500)
1000 loops, best of 3: 682 µs per loop
>>> %timeit cartesian_product(x500, y500)
1000 loops, best of 3: 1.55 ms per loop
以下では、これらの方針に沿って実行した以前のテストについていくつか詳しく説明します。これらのアプローチの相対的なパフォーマンスは、さまざまなハードウェアとさまざまなバージョンのPythonおよびで時間の経過とともに変化しましたnumpy。の最新バージョンを使用している人にはすぐには役立ちませんがnumpy、この回答の最初のバージョンから状況がどのように変化したかを示しています。
meshgrid+dstack現在受け入れられている回答はtile、とを使用しrepeatて2つのアレイを一緒にブロードキャストします。しかし、meshgrid関数は実質的に同じことをします。tile転置にrepeat渡される前の出力は次のとおりです。
In [1]: import numpy
In [2]: x = numpy.array([1,2,3])
...: y = numpy.array([4,5])
...:
In [3]: [numpy.tile(x, len(y)), numpy.repeat(y, len(x))]
Out[3]: [array([1, 2, 3, 1, 2, 3]), array([4, 4, 4, 5, 5, 5])]
そして、これが次の出力ですmeshgrid：
In [4]: numpy.meshgrid(x, y)
Out[4]:
[array([[1, 2, 3],
[1, 2, 3]]), array([[4, 4, 4],
[5, 5, 5]])]
ご覧のとおり、ほぼ同じです。まったく同じ結果を得るには、結果の形状を変更するだけで済みます。
In [5]: xt, xr = numpy.meshgrid(x, y)
...: [xt.ravel(), xr.ravel()]
Out[5]: [array([1, 2, 3, 1, 2, 3]), array([4, 4, 4, 5, 5, 5])]
ただし、この時点で形状を変更するのではなく、の出力をに渡して後で形状を変更することで、作業を節約できmeshgridますdstack。
In [6]: numpy.dstack(numpy.meshgrid(x, y)).reshape(-1, 2)
Out[6]:
array([[1, 4],
[2, 4],
[3, 4],
[1, 5],
[2, 5],
[3, 5]])
このコメントの主張に反して、異なる入力が異なる形状の出力を生成するという証拠は見られませんでした。上記のように、それらは非常に類似したことを行うので、それらが行われた場合は非常に奇妙になります。反例を見つけたら教えてください。
meshgrid+dstack対repeat+transposeこれら2つのアプローチの相対的なパフォーマンスは、時間の経過とともに変化しました。以前のバージョンのPython（2.7）では、meshgrid+dstackを使用した結果は、小さな入力に対して著しく高速でした。（これらのテストはこの回答の古いバージョンからのものであることに注意してください。）定義：
>>> def repeat_product(x, y):
... return numpy.transpose([numpy.tile(x, len(y)),
numpy.repeat(y, len(x))])
...
>>> def dstack_product(x, y):
... return numpy.dstack(numpy.meshgrid(x, y)).reshape(-1, 2)
...
適度なサイズの入力の場合、大幅なスピードアップが見られました。しかし、新しいマシンで、Pythonの最新バージョン（3.6.1）とnumpy（1.12.1）を使用してこれらのテストを再試行しました。現在、2つのアプローチはほぼ同じです。
古いテスト
>>> x, y = numpy.arange(500), numpy.arange(500)
>>> %timeit repeat_product(x, y)
10 loops, best of 3: 62 ms per loop
>>> %timeit dstack_product(x, y)
100 loops, best of 3: 12.2 ms per loop
新しいテスト
In [7]: x, y = numpy.arange(500), numpy.arange(500)
In [8]: %timeit repeat_product(x, y)
1.32 ms ± 24.7 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
In [9]: %timeit dstack_product(x, y)
1.26 ms ± 8.47 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
いつものように、YMMVですが、これは、Pythonとnumpyの最近のバージョンでは、これらが交換可能であることを示唆しています。
一般に、組み込み関数を使用すると、入力が小さい場合は高速になると予想されますが、入力が大きい場合は、専用関数の方が高速になる可能性があります。さらに、一般化されたn次元の製品の場合、tile明確repeatな高次元の類似体がないため、役に立ちません。したがって、専用関数の動作も調査する価値があります。
関連するテストのほとんどはこの回答の冒頭に表示されますが、以前のバージョンのPythonで実行されたテストのいくつかnumpyを比較のために示します。
別の回答cartesianで定義された関数は、より大きな入力に対してかなりうまく機能するために使用されました。（上記の関数と同じです。）と比較するために、2次元のみを使用します。cartesian_product_recursivecartesiandstack_prodct
ここでも、古いテストでは大きな違いが見られましたが、新しいテストではほとんど違いが見られませんでした。
古いテスト
>>> x, y = numpy.arange(1000), numpy.arange(1000)
>>> %timeit cartesian([x, y])
10 loops, best of 3: 25.4 ms per loop
>>> %timeit dstack_product(x, y)
10 loops, best of 3: 66.6 ms per loop
新しいテスト
In [10]: x, y = numpy.arange(1000), numpy.arange(1000)
In [11]: %timeit cartesian([x, y])
12.1 ms ± 199 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
In [12]: %timeit dstack_product(x, y)
12.7 ms ± 334 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
以前のように、それでも小さなスケールでdstack_productビートします。cartesian
新しいテスト（冗長な古いテストは表示されていません）
In [13]: x, y = numpy.arange(100), numpy.arange(100)
In [14]: %timeit cartesian([x, y])
215 µs ± 4.75 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
In [15]: %timeit dstack_product(x, y)
65.7 µs ± 1.15 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
これらの違いは、興味深いものであり、記録する価値があると思います。しかし、彼らは最終的には学術的です。この回答の冒頭のテストが示したように、これらのバージョンはすべて、ほとんどの場合cartesian_product、この回答の冒頭で定義されたよりも遅くなります。これは、この質問に対する回答の中で最速の実装よりも少し遅いです。
>>> numpy.transpose([numpy.tile(x, len(y)), numpy.repeat(y, len(x))])
array([[1, 4],
[2, 4],
[3, 4],
[1, 5],
[2, 5],
[3, 5]])
N 配列のデカルト積を計算するための一般的なソリューションについては、numpy を使用して 2 つの配列のすべての組み合わせの配列を作成するを参照してください。
Pythonで通常のリスト内包表記を行うことができます
x = numpy.array([1,2,3])
y = numpy.array([4,5])
[[x0, y0] for x0 in x for y0 in y]
あなたに与えるべきもの
[[1, 4], [1, 5], [2, 4], [2, 5], [3, 4], [3, 5]]