Iloczyn kartezjański punktów tablicy xiy w pojedynczą tablicę punktów 2D

Mam dwie tablice numpy, które definiują osie x i y siatki. Na przykład:

x = numpy.array([1,2,3])
y = numpy.array([4,5])

Chciałbym wygenerować iloczyn kartezjański tych tablic, aby wygenerować:

array([[1,4],[2,4],[3,4],[1,5],[2,5],[3,5]])

W pewnym sensie nie jest to strasznie nieefektywne, ponieważ muszę to robić wiele razy w pętli. Zakładam, że przekonwertowanie ich na listę Pythona i użycie itertools.productiz powrotem do tablicy numpy nie jest najbardziej wydajną formą.

>>> 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]])

Zobacz Używanie numpy do budowania tablicy wszystkich kombinacji dwóch tablic, aby zapoznać się z ogólnym rozwiązaniem obliczania iloczynu kartezjańskiego N tablic.

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.

Możesz po prostu zrobić normalne rozumienie list w Pythonie

x = numpy.array([1,2,3])
y = numpy.array([4,5])
[[x0, y0] for x0 in x for y0 in y]

co powinno ci dać

[[1, 4], [1, 5], [2, 4], [2, 5], [3, 4], [3, 5]]

Interesowało mnie to również i dokonałem małego porównania wydajności, być może nieco wyraźniejszego niż w odpowiedzi @ senderle.

Dla dwóch tablic (przypadek klasyczny):

Dla czterech tablic:

(Zauważ, że długość tablic to tylko kilkadziesiąt wpisów tutaj.)

Kod do odtworzenia działek:

from functools import reduce
import itertools
import numpy
import perfplot


def dstack_product(arrays):
    return numpy.dstack(numpy.meshgrid(*arrays, indexing="ij")).reshape(-1, len(arrays))


# Generalized N-dimensional products
def cartesian_product(arrays):
    la = len(arrays)
    dtype = numpy.find_common_type([a.dtype for a in 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 = reduce(numpy.multiply, broadcasted[0].shape), len(broadcasted)
    dtype = numpy.find_common_type([a.dtype for a in 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)))


perfplot.show(
    setup=lambda n: 2 * (numpy.arange(n, dtype=float),),
    n_range=[2 ** k for k in range(13)],
    # setup=lambda n: 4 * (numpy.arange(n, dtype=float),),
    # n_range=[2 ** k for k in range(6)],
    kernels=[
        dstack_product,
        cartesian_product,
        cartesian_product_transpose,
        cartesian_product_recursive,
        cartesian_product_itertools,
    ],
    logx=True,
    logy=True,
    xlabel="len(a), len(b)",
    equality_check=None,
)

Opierając się na wzorowej pracy naziemnej @ senderle, opracowałem dwie wersje - jedną dla C i jedną dla układów Fortran - które często są nieco szybsze.

• cartesian_product_transpose_ppjest - w przeciwieństwie do @ senderle, cartesian_product_transposektóry używa zupełnie innej strategii - wersja, cartesion_productktóra wykorzystuje bardziej korzystny układ pamięci transpozycji + kilka bardzo drobnych optymalizacji.
• cartesian_product_ppzachowuje oryginalny układ pamięci. To, co sprawia, że ​​jest szybki, to ciągłe kopiowanie. Ciągłe kopie okazują się być o wiele szybsze, niż kopiowanie całego bloku pamięci, mimo że tylko jego część zawiera prawidłowe dane, a nie tylko kopiowanie ważnych bitów.

Niektóre perfplots. Zrobiłem osobne dla układów C i Fortran, bo to są różne zadania IMO.

Nazwy kończące się na „pp” to moje podejście.

1) wiele drobnych czynników (po 2 elementy)

2) wiele małych czynników (po 4 elementy)

3) trzy czynniki o jednakowej długości

4) dwa czynniki o jednakowej długości

Kod (trzeba wykonać oddzielne przebiegi dla każdej działki b / c Nie mogłem wymyślić, jak zresetować; również muszę odpowiednio edytować / komentować w / wy):

import numpy
import numpy as np
from functools import reduce
import itertools
import timeit
import perfplot

def dstack_product(arrays):
    return numpy.dstack(
        numpy.meshgrid(*arrays, indexing='ij')
        ).reshape(-1, len(arrays))

def cartesian_product_transpose_pp(arrays):
    la = len(arrays)
    dtype = numpy.result_type(*arrays)
    arr = numpy.empty((la, *map(len, arrays)), dtype=dtype)
    idx = slice(None), *itertools.repeat(None, la)
    for i, a in enumerate(arrays):
        arr[i, ...] = a[idx[:la-i]]
    return arr.reshape(la, -1).T

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 itertools import accumulate, repeat, chain

def cartesian_product_pp(arrays, out=None):
    la = len(arrays)
    L = *map(len, arrays), la
    dtype = numpy.result_type(*arrays)
    arr = numpy.empty(L, dtype=dtype)
    arrs = *accumulate(chain((arr,), repeat(0, la-1)), np.ndarray.__getitem__),
    idx = slice(None), *itertools.repeat(None, la-1)
    for i in range(la-1, 0, -1):
        arrs[i][..., i] = arrays[i][idx[:la-i]]
        arrs[i-1][1:] = arrs[i]
    arr[..., 0] = arrays[0][idx]
    return arr.reshape(-1, la)

def cartesian_product_itertools(arrays):
    return numpy.array(list(itertools.product(*arrays)))


# 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

### Test code ###
if False:
  perfplot.save('cp_4el_high.png',
    setup=lambda n: n*(numpy.arange(4, dtype=float),),
                n_range=list(range(6, 11)),
    kernels=[
        dstack_product,
        cartesian_product_recursive,
        cartesian_product,
#        cartesian_product_transpose,
        cartesian_product_pp,
#        cartesian_product_transpose_pp,
        ],
    logx=False,
    logy=True,
    xlabel='#factors',
    equality_check=None
    )
else:
  perfplot.save('cp_2f_T.png',
    setup=lambda n: 2*(numpy.arange(n, dtype=float),),
    n_range=[2**k for k in range(5, 11)],
    kernels=[
#        dstack_product,
#        cartesian_product_recursive,
#        cartesian_product,
        cartesian_product_transpose,
#        cartesian_product_pp,
        cartesian_product_transpose_pp,
        ],
    logx=True,
    logy=True,
    xlabel='length of each factor',
    equality_check=None
    )

Od października 2017 r. Numpy ma teraz funkcję ogólną, np.stackktóra przyjmuje parametr osi. Używając go, możemy mieć „uogólniony produkt kartezjański” przy użyciu techniki „dstack and meshgrid”:

import numpy as np
def cartesian_product(*arrays):
    ndim = len(arrays)
    return np.stack(np.meshgrid(*arrays), axis=-1).reshape(-1, ndim)

Uwaga dotycząca axis=-1parametru. Jest to ostatnia (najbardziej wewnętrzna) oś w wyniku. Jest to równoważne użyciu axis=ndim.

Jeszcze jeden komentarz, ponieważ iloczyn kartezjański szybko się wysadza, chyba że z jakiegoś powodu musimy zrealizować tablicę w pamięci, jeśli iloczyn jest bardzo duży, możemy chcieć itertoolswykorzystać wartości w locie.

Używałem odpowiedzi @kennytm przez jakiś czas, ale próbując zrobić to samo w TensorFlow, okazało się, że TensorFlow nie ma odpowiednika numpy.repeat(). Po krótkich eksperymentach, myślę, że znalazłem bardziej ogólne rozwiązanie dla dowolnych wektorów punktów.

Dla numpy:

import numpy as np

def cartesian_product(*args: np.ndarray) -> np.ndarray:
    """
    Produce the cartesian product of arbitrary length vectors.

    Parameters
    ----------
    np.ndarray args
        vector of points of interest in each dimension

    Returns
    -------
    np.ndarray
        the cartesian product of size [m x n] wherein:
            m = prod([len(a) for a in args])
            n = len(args)
    """
    for i, a in enumerate(args):
        assert a.ndim == 1, "arg {:d} is not rank 1".format(i)
    return np.concatenate([np.reshape(xi, [-1, 1]) for xi in np.meshgrid(*args)], axis=1)

i dla TensorFlow:

import tensorflow as tf

def cartesian_product(*args: tf.Tensor) -> tf.Tensor:
    """
    Produce the cartesian product of arbitrary length vectors.

    Parameters
    ----------
    tf.Tensor args
        vector of points of interest in each dimension

    Returns
    -------
    tf.Tensor
        the cartesian product of size [m x n] wherein:
            m = prod([len(a) for a in args])
            n = len(args)
    """
    for i, a in enumerate(args):
        tf.assert_rank(a, 1, message="arg {:d} is not rank 1".format(i))
    return tf.concat([tf.reshape(xi, [-1, 1]) for xi in tf.meshgrid(*args)], axis=1)

Pakiet Scikit-learn ma szybką implementację dokładnie tego:

from sklearn.utils.extmath import cartesian
product = cartesian((x,y))

Zauważ, że konwencja tej implementacji różni się od tego, czego chcesz, jeśli zależy Ci na kolejności danych wyjściowych. Aby dokładnie zamówić, możesz to zrobić

product = cartesian((y,x))[:, ::-1]

Mówiąc bardziej ogólnie, jeśli masz dwie 2d tablice numpy a i b i chcesz połączyć każdy wiersz a z każdym wierszem b (iloczyn kartezjański wierszy, coś w rodzaju złączenia w bazie danych), możesz użyć tej metody :

import numpy
def join_2d(a, b):
    assert a.dtype == b.dtype
    a_part = numpy.tile(a, (len(b), 1))
    b_part = numpy.repeat(b, len(a), axis=0)
    return numpy.hstack((a_part, b_part))

Najszybciej można uzyskać połączenie wyrażenia generatora z funkcją map:

import numpy
import datetime
a = np.arange(1000)
b = np.arange(200)

start = datetime.datetime.now()

foo = (item for sublist in [list(map(lambda x: (x,i),a)) for i in b] for item in sublist)

print (list(foo))

print ('execution time: {} s'.format((datetime.datetime.now() - start).total_seconds()))

Wyniki (właściwie cała wynikowa lista jest drukowana):

[(0, 0), (1, 0), ...,(998, 199), (999, 199)]
execution time: 1.253567 s

lub używając podwójnego wyrażenia generatora:

a = np.arange(1000)
b = np.arange(200)

start = datetime.datetime.now()

foo = ((x,y) for x in a for y in b)

print (list(foo))

print ('execution time: {} s'.format((datetime.datetime.now() - start).total_seconds()))

Wyniki (drukowana cała lista):

[(0, 0), (1, 0), ...,(998, 199), (999, 199)]
execution time: 1.187415 s

Weź pod uwagę, że większość czasu obliczeń przypada na polecenie drukowania. Poza tym obliczenia generatora są przyzwoicie wydajne. Bez drukowania czasy obliczeń to:

execution time: 0.079208 s

dla wyrażenia generatora + funkcja mapy i:

execution time: 0.007093 s

dla wyrażenia podwójnego generatora.

Jeśli tak naprawdę chcesz obliczyć rzeczywisty iloczyn każdej z par współrzędnych, najszybszym rozwiązaniem jest rozwiązanie jako iloczyn liczbowy macierzy:

a = np.arange(1000)
b = np.arange(200)

start = datetime.datetime.now()

foo = np.dot(np.asmatrix([[i,0] for i in a]), np.asmatrix([[i,0] for i in b]).T)

print (foo)

print ('execution time: {} s'.format((datetime.datetime.now() - start).total_seconds()))

Wyjścia:

 [[     0      0      0 ...,      0      0      0]
 [     0      1      2 ...,    197    198    199]
 [     0      2      4 ...,    394    396    398]
 ..., 
 [     0    997   1994 ..., 196409 197406 198403]
 [     0    998   1996 ..., 196606 197604 198602]
 [     0    999   1998 ..., 196803 197802 198801]]
execution time: 0.003869 s

i bez nadruku (w tym przypadku nie oszczędza dużo, ponieważ drukowany jest tylko malutki kawałek matrycy):

execution time: 0.003083 s

Można to również łatwo zrobić za pomocą metody itertools.product

from itertools import product
import numpy as np

x = np.array([1, 2, 3])
y = np.array([4, 5])
cart_prod = np.array(list(product(*[x, y])),dtype='int32')

Wynik: tablica ([
[1, 4],
[1, 5],
[2, 4],
[2, 5],
[3, 4],
[3, 5]], dtype = int32)

Czas wykonania: 0,000155 s

W konkretnym przypadku, gdy musisz wykonać proste operacje, takie jak dodawanie na każdej parze, możesz wprowadzić dodatkowy wymiar i pozwolić, aby nadawanie wykonało zadanie:

>>> a, b = np.array([1,2,3]), np.array([10,20,30])
>>> a[None,:] + b[:,None]
array([[11, 12, 13],
       [21, 22, 23],
       [31, 32, 33]])

Nie jestem pewien, czy istnieje podobny sposób na uzyskanie samych par.