Kanoniczny cartesian_product(prawie)
Istnieje wiele podejść do tego problemu o różnych właściwościach. Niektóre są szybsze niż inne, a niektóre mają bardziej ogólne zastosowanie. Po wielu testach i poprawkach odkryłem, że następująca funkcja, która oblicza n-wymiar cartesian_product, jest szybsza niż większość innych dla wielu danych wejściowych. Aby zapoznać się z parą podejść, które są nieco bardziej złożone, ale w wielu przypadkach są nawet nieco szybsze, zobacz odpowiedź Paula Panzera .
Biorąc pod uwagę tę odpowiedź, nie jest to już najszybsza implementacja produktu kartezjańskiego numpy, o której wiem. Myślę jednak, że jego prostota nadal będzie stanowić przydatny punkt odniesienia dla przyszłych ulepszeń:
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)
Warto wspomnieć, że ta funkcja używa ix_w nietypowy sposób; podczas gdy udokumentowane użycie ix_polega na generowaniu indeksów w tablicy, tak się składa, że tablice o tym samym kształcie mogą być używane do rozgłaszania. Wielkie podziękowania dla mgilson , który zainspirował mnie do spróbowania ix_tego sposobu, oraz dla unutbu , który udzielił niezwykle pomocnych opinii na temat tej odpowiedzi, w tym sugestii użycia numpy.result_type.
Godne uwagi alternatywy
Czasami szybsze jest pisanie ciągłych bloków pamięci w kolejności Fortran. To podstawa tej alternatywy, cartesian_product_transposektóra okazała się szybsza na niektórych urządzeniach niż cartesian_product(patrz poniżej). Jednak odpowiedź Paula Panzera, która opiera się na tej samej zasadzie, jest jeszcze szybsza. Mimo to zamieszczam to tutaj dla zainteresowanych czytelników:
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
Po zrozumieniu podejścia Panzera, napisałem nową wersję, która jest prawie tak szybka jak jego i prawie tak prosta, jak 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
Wydaje się, że ma to stałe obciążenie, które sprawia, że działa wolniej niż Panzer przy małych nakładach. Ale w przypadku większych danych wejściowych we wszystkich testach, które przeprowadziłem, działa równie dobrze, jak jego najszybsza implementacja ( cartesian_product_transpose_pp).
W kolejnych sekcjach zamieszczam kilka testów innych alternatyw. Są one teraz nieco nieaktualne, ale zamiast dublować wysiłek, zdecydowałem się je tutaj zostawić z interesu historycznego. Aby zapoznać się z aktualnymi testami, zobacz odpowiedź Panzera, a także odpowiedź Nico Schlömera .
Testy pod kątem alternatyw
Oto zestaw testów, które pokazują wzrost wydajności, jaki zapewniają niektóre z tych funkcji w porównaniu z wieloma alternatywami. Wszystkie pokazane tutaj testy zostały przeprowadzone na czterordzeniowym komputerze z systemem Mac OS 10.12.5, Python 3.6.1 i numpy1.12.1. Wiadomo, że różnice w sprzęcie i oprogramowaniu dają różne wyniki, więc YMMV. Przeprowadź te testy dla siebie, aby mieć pewność!
Definicje:
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)]
Wyniki testu:
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)
We wszystkich przypadkach, cartesian_productjak określono na początku tej odpowiedzi, jest najszybsza.
W przypadku funkcji, które akceptują dowolną liczbę tablic wejściowych, warto również sprawdzić wydajność len(arrays) > 2. (Dopóki nie mogę ustalić, dlaczego cartesian_product_recursivew tym przypadku jest wyświetlany błąd, usunąłem go z tych testów).
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)
Jak pokazują te testy, cartesian_productpozostaje konkurencyjny, dopóki liczba tablic wejściowych nie wzrośnie powyżej (w przybliżeniu) czterech. Potem cartesian_product_transposema niewielką przewagę.
Warto powtórzyć, że użytkownicy z innym sprzętem i systemami operacyjnymi mogą zobaczyć inne wyniki. Na przykład unutbu raportuje następujące wyniki tych testów przy użyciu Ubuntu 14.04, Python 3.4.3 i 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
Poniżej przedstawiam kilka szczegółów dotyczących wcześniejszych testów, które przeprowadziłem w ten sposób. Względna wydajność tych podejść zmieniała się w czasie dla różnych urządzeń i różnych wersji języków Python i numpy. Chociaż nie jest to od razu przydatne dla osób korzystających z aktualnych wersji numpy, ilustruje, jak wiele się zmieniło od czasu pierwszej wersji tej odpowiedzi.
Prosta alternatywa: meshgrid+dstack
Aktualnie zaakceptowana odpowiedź wykorzystuje tilei repeatdo nadawania razem dwóch tablic. Ale meshgridfunkcja robi praktycznie to samo. Oto wynik działania tilei repeatprzed przekazaniem do transpozycji:
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])]
A oto wynik meshgrid:
In [4]: numpy.meshgrid(x, y)
Out[4]:
[array([[1, 2, 3],
[1, 2, 3]]), array([[4, 4, 4],
[5, 5, 5]])]
Jak widać, jest prawie identyczny. Musimy tylko zmienić kształt wyniku, aby uzyskać dokładnie ten sam wynik.
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])]
Zamiast zmieniać kształt w tym momencie, moglibyśmy przekazać wynik meshgriddo, dstacka następnie zmienić kształt, co oszczędza trochę pracy:
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]])
Wbrew twierdzeniom zawartym w tym komentarzu , nie widziałem dowodów na to, że różne dane wejściowe będą dawać różnie ukształtowane wyniki, a jak pokazuje powyższe, robią bardzo podobne rzeczy, więc byłoby dość dziwne, gdyby tak było. Daj mi znać, jeśli znajdziesz kontrprzykład.
Testowanie meshgrid+ dstackvs. repeat+transpose
Względna wydajność tych dwóch podejść zmieniała się w czasie. We wcześniejszej wersji Pythona (2.7) wynik użycia meshgrid+ dstackbył zauważalnie szybszy w przypadku małych danych wejściowych. (Zauważ, że te testy pochodzą ze starej wersji tej odpowiedzi). Definicje:
>>> 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)
...
W przypadku wejścia o średnim rozmiarze zauważyłem znaczne przyspieszenie. Ale powtórzyłem te testy z nowszymi wersjami Pythona (3.6.1) i numpy(1.12.1) na nowszej maszynie. Te dwa podejścia są teraz prawie identyczne.
Stary test
>>> 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
Nowy test
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)
Jak zawsze YMMV, ale sugeruje to, że w najnowszych wersjach Pythona i numpy są one wymienne.
Uogólnione funkcje produktu
Ogólnie rzecz biorąc, możemy oczekiwać, że korzystanie z funkcji wbudowanych będzie szybsze w przypadku małych danych wejściowych, podczas gdy w przypadku dużych danych wejściowych funkcja specjalnie zbudowana może być szybsza. Ponadto dla uogólnionego produktu n-wymiarowego tilei repeatnie pomoże, ponieważ nie mają wyraźnych analogów wyższego wymiaru. Dlatego warto zbadać również zachowanie funkcji specjalnie zaprojektowanych.
Większość odpowiednich testów pojawia się na początku tej odpowiedzi, ale oto kilka testów przeprowadzonych na wcześniejszych wersjach Pythona i numpydla porównania.
cartesianFunkcja zdefiniowana w innym odpowiedź wykorzystywane do wykonywania bardzo dobrze dla większych nakładów. (To samo jak funkcja o nazwie cartesian_product_recursivepowyżej). W celu porównania cartesiando dstack_prodct, używamy tylko dwa wymiary.
Tutaj ponownie stary test wykazał znaczną różnicę, podczas gdy nowy test prawie nie wykazuje żadnej.
Stary test
>>> 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
Nowy test
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)
Jak poprzednio, dstack_productnadal bije cartesianw mniejszych skalach.
Nowy test ( nie pokazano nadmiarowego starego testu )
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)
Myślę, że te rozróżnienia są interesujące i warte odnotowania; ale w końcu są akademickimi. Jak pokazały testy na początku tej odpowiedzi, wszystkie te wersje są prawie zawsze wolniejsze niż cartesian_productzdefiniowane na samym początku tej odpowiedzi - co samo w sobie jest nieco wolniejsze niż najszybsze implementacje wśród odpowiedzi na to pytanie.
