Kartesisches Produkt von x- und y-Array-Punkten in ein einzelnes Array von 2D-Punkten

147

Ich habe zwei Numpy-Arrays, die die x- und y-Achse eines Gitters definieren. Beispielsweise:

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

Ich möchte das kartesische Produkt dieser Arrays generieren, um Folgendes zu generieren:

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

In gewisser Weise ist das nicht besonders ineffizient, da ich dies viele Male in einer Schleife tun muss. Ich gehe davon aus, dass das Konvertieren in eine Python-Liste und das Verwenden itertools.productund Zurück in ein Numpy-Array nicht die effizienteste Form ist.

Reich
quelle
Mir ist aufgefallen, dass der teuerste Schritt im itertools-Ansatz die endgültige Konvertierung von Liste zu Array ist. Ohne diesen letzten Schritt ist es doppelt so schnell wie Kens Beispiel.
Alexey Lebedev

Antworten:

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

Eine allgemeine Lösung zur Berechnung des kartesischen Produkts von N Arrays finden Sie unter Verwenden von numpy zum Erstellen eines Arrays aller Kombinationen von zwei Arrays.

kennytm
quelle
1
Ein Vorteil dieses Ansatzes besteht darin, dass eine konsistente Ausgabe für Arrays derselben Größe erzeugt wird. Der meshgrid+ dstack-Ansatz ist zwar in einigen Fällen schneller, kann jedoch zu Fehlern führen, wenn Sie erwarten, dass das kartesische Produkt für Arrays derselben Größe in derselben Reihenfolge erstellt wird.
tlnagy
3
@tlnagy, ich habe keinen Fall bemerkt, in dem dieser Ansatz andere Ergebnisse liefert als die von meshgrid+ dstack. Könnten Sie ein Beispiel posten?
senderle
148

Ein kanonischer cartesian_product(fast)

Es gibt viele Ansätze für dieses Problem mit unterschiedlichen Eigenschaften. Einige sind schneller als andere, andere sind allgemeiner. Nach vielen Tests und Optimierungen habe ich festgestellt, dass die folgende Funktion, die eine n-Dimension berechnet cartesian_product, für viele Eingaben schneller ist als die meisten anderen. Für ein paar Ansätze, die etwas komplexer sind, aber in vielen Fällen sogar etwas schneller, siehe die Antwort von Paul Panzer .

Angesichts dieser Antwort ist dies nicht mehr die schnellste Implementierung des kartesischen Produkts numpy, die mir bekannt ist. Ich denke jedoch, dass seine Einfachheit es weiterhin zu einem nützlichen Maßstab für zukünftige Verbesserungen machen wird:

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)

Es ist erwähnenswert, dass diese Funktion auf ix_ungewöhnliche Weise verwendet wird. Während die dokumentierte Verwendung von darin ix_besteht, Indizes in einem Array zu generieren , kommt es nur so vor, dass Arrays mit derselben Form für die Broadcast-Zuweisung verwendet werden können. Vielen Dank an mgilson , der mich dazu inspiriert hat, ix_diesen Weg zu versuchen , und an unutbu , der einige äußerst hilfreiche Rückmeldungen zu dieser Antwort gegeben hat, einschließlich des Vorschlags zur Verwendung numpy.result_type.

Bemerkenswerte Alternativen

Es ist manchmal schneller, zusammenhängende Speicherblöcke in Fortran-Reihenfolge zu schreiben. Das ist die Basis dieser Alternative, cartesian_product_transposedie sich auf einigen Hardwarekomponenten als schneller erwiesen hat als cartesian_product(siehe unten). Die Antwort von Paul Panzer, die dasselbe Prinzip verwendet, ist jedoch noch schneller. Trotzdem füge ich dies hier für interessierte Leser ein:

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

Nachdem ich Panzers Ansatz verstanden hatte, schrieb ich eine neue Version, die fast so schnell ist wie seine und fast so einfach wie 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

Dies scheint einen zeitlich konstanten Overhead zu haben, der es bei kleinen Eingaben langsamer laufen lässt als bei Panzer. Bei größeren Eingaben funktioniert es in allen von mir durchgeführten Tests genauso gut wie seine schnellste Implementierung ( cartesian_product_transpose_pp).

In den folgenden Abschnitten füge ich einige Tests anderer Alternativen hinzu. Diese sind jetzt etwas veraltet, aber anstatt doppelte Anstrengungen zu unternehmen, habe ich beschlossen, sie aus historischem Interesse hier zu lassen. Aktuelle Tests finden Sie in Panzers Antwort sowie in der von Nico Schlömer .

Tests gegen Alternativen

Hier finden Sie eine Reihe von Tests, die die Leistungssteigerung zeigen, die einige dieser Funktionen im Vergleich zu einer Reihe von Alternativen bieten. Alle hier gezeigten Tests wurden auf einem Quad-Core-Computer unter Mac OS 10.12.5, Python 3.6.1 und numpy1.12.1 durchgeführt. Es ist bekannt, dass Variationen von Hardware und Software zu unterschiedlichen Ergebnissen führen, so YMMV. Führen Sie diese Tests selbst durch, um sicherzugehen!

Definitionen:

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

Testergebnisse:

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)

In allen Fällen ist cartesian_productdie am Anfang dieser Antwort definierte Antwort am schnellsten.

Für Funktionen, die eine beliebige Anzahl von Eingabearrays akzeptieren, lohnt es sich, die Leistung auch dann zu überprüfen len(arrays) > 2. (Bis ich feststellen kann, warum cartesian_product_recursivein diesem Fall ein Fehler auftritt, habe ich ihn aus diesen Tests entfernt.)

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)

Wie diese Tests zeigen, cartesian_productbleibt es wettbewerbsfähig, bis die Anzahl der Eingabearrays über (ungefähr) vier steigt. Danach cartesian_product_transposehat eine leichte Kante.

Es sollte wiederholt werden, dass Benutzer mit anderer Hardware und Betriebssystemen möglicherweise andere Ergebnisse sehen. Unutbu-Berichte zeigen beispielsweise die folgenden Ergebnisse für diese Tests mit Ubuntu 14.04, Python 3.4.3 und 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

Im Folgenden gehe ich auf einige Details zu früheren Tests ein, die ich in diesem Sinne durchgeführt habe. Die relative Leistung dieser Ansätze hat sich im Laufe der Zeit für verschiedene Hardware und verschiedene Versionen von Python und geändert numpy. Es ist zwar nicht sofort nützlich für Benutzer, die aktuelle Versionen von verwenden numpy, zeigt jedoch, wie sich die Dinge seit der ersten Version dieser Antwort geändert haben.

Eine einfache Alternative: meshgrid+dstack

Die aktuell akzeptierte Antwort verwendet tileund repeatsendet zwei Arrays zusammen. Aber die meshgridFunktion macht praktisch das Gleiche. Hier ist die Ausgabe von tileundrepeat bevor sie zur Transponierung übergeben wird:

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

Und hier ist die Ausgabe von meshgrid :

In [4]: numpy.meshgrid(x, y)
Out[4]: 
[array([[1, 2, 3],
        [1, 2, 3]]), array([[4, 4, 4],
        [5, 5, 5]])]

Wie Sie sehen können, ist es fast identisch. Wir müssen nur das Ergebnis umformen, um genau das gleiche Ergebnis zu erzielen.

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

Anstatt an dieser Stelle eine Umformung vorzunehmen, könnten wir die Ausgabe von meshgridan übergeben dstackund anschließend umformen, was einige Arbeit spart:

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

Entgegen der Behauptung in diesem Kommentar habe ich keine Beweise dafür gesehen, dass unterschiedliche Eingaben unterschiedlich geformte Ausgaben erzeugen, und wie oben gezeigt, tun sie sehr ähnliche Dinge, so dass es ziemlich seltsam wäre, wenn sie dies tun würden. Bitte lassen Sie mich wissen, wenn Sie ein Gegenbeispiel finden.

Testen meshgrid+dstack gegen repeat+transpose

Die relative Leistung dieser beiden Ansätze hat sich im Laufe der Zeit geändert. In einer früheren Version von Python (2.7) war das Ergebnis mit meshgrid+ dstackbei kleinen Eingaben deutlich schneller. (Beachten Sie, dass diese Tests aus einer alten Version dieser Antwort stammen.) Definitionen:

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

Bei mäßig großen Eingaben sah ich eine deutliche Beschleunigung. Ich habe diese Tests jedoch mit neueren Versionen von Python (3.6.1) und numpy(1.12.1) auf einem neueren Computer wiederholt. Die beiden Ansätze sind jetzt fast identisch.

Alter 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

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

Wie immer YMMV, aber dies deutet darauf hin, dass diese in neueren Versionen von Python und Numpy austauschbar sind.

Verallgemeinerte Produktfunktionen

Im Allgemeinen können wir erwarten, dass die Verwendung integrierter Funktionen für kleine Eingaben schneller ist, während für große Eingaben eine speziell entwickelte Funktion möglicherweise schneller ist. Weiterhin für ein verallgemeinertes n-dimensionales Produkt tileundrepeat wird nicht helfen, weil sie keine klaren höherdimensionalen Analoga haben. Es lohnt sich also, auch das Verhalten von speziell entwickelten Funktionen zu untersuchen.

Die meisten relevanten Tests erscheinen am Anfang dieser Antwort, aber hier sind einige der Tests, die mit früheren Versionen von Python und numpyzum Vergleich durchgeführt wurden.

Die cartesianin einer anderen Antwort definierte Funktion wurde für größere Eingaben verwendet. (Es ist das gleiche wie die Funktion aufgerufen cartesian_product_recursiveoben.) Um zu vergleichen , cartesianzu dstack_prodctverwenden wir nur zwei Dimensionen.

Auch hier zeigte der alte Test einen signifikanten Unterschied, während der neue Test fast keinen zeigt.

Alter 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

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

Nach wie vor dstack_productschlägt cartesianin kleineren Maßstäben.

Neuer Test ( redundanter alter Test nicht gezeigt )

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)

Diese Unterscheidungen sind meiner Meinung nach interessant und es lohnt sich, sie aufzuzeichnen. aber sie sind am Ende akademisch. Wie die Tests am Anfang dieser Antwort gezeigt haben, sind alle diese Versionen fast immer langsamer als cartesian_productam Anfang dieser Antwort definiert - was selbst etwas langsamer ist als die schnellsten Implementierungen unter den Antworten auf diese Frage.

senderle
quelle
1
und das Hinzufügen dtype=objectin arr = np.empty( )würde die Verwendung verschiedener Typen im Produkt ermöglichen, z arrays = [np.array([1,2,3]), ['str1', 'str2']].
user3820991
Vielen Dank für Ihre innovativen Lösungen. Ich dachte nur, Sie möchten wissen, dass einige Benutzer möglicherweise cartesian_product_tranposeschneller finden als cartesian_productje nach Betriebssystem, Python oder Numpy-Version. Zum Beispiel auf Ubuntu 14.04, python3.4.3, numpy 1.14.0.dev0 + b7050a9, %timeit cartesian_product_transpose(x500,y500)Erträge , 1000 loops, best of 3: 682 µs per loopwährend die %timeit cartesian_product(x500,y500)Renditen 1000 loops, best of 3: 1.55 ms per loop. Ich finde auch, cartesian_product_transposedass es schneller sein kann, wenn len(arrays) > 2.
Unutbu
Zusätzlich cartesian_productgibt ein Array von Gleitkommazahlen dtype während cartesian_product_transposegibt ein Array von derselben dtype wie der erste (ausgestrahlt) -Array. Die Möglichkeit, dtype bei der Arbeit mit ganzzahligen Arrays beizubehalten, kann für Benutzer ein Grund sein, dies zu bevorzugen cartesian_product_transpose.
Unutbu
@unutbu Nochmals vielen Dank - wie ich hätte wissen müssen, bringt das Klonen des dtype nicht nur Bequemlichkeit; In einigen Fällen wird der Code um weitere 20 bis 30% beschleunigt.
senderle
1
@senderle: Wow, das ist schön! Außerdem kam mir gerade dtype = np.find_common_type([arr.dtype for arr in arrays], [])der Gedanke , dass so etwas verwendet werden könnte, um den gemeinsamen d-Typ aller Arrays zu finden, anstatt den Benutzer zu zwingen, das Array zu platzieren, das den d-Typ zuerst steuert.
Unutbu
44

Sie können in Python einfach ein normales Listenverständnis durchführen

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

was sollte dir geben

[[1, 4], [1, 5], [2, 4], [2, 5], [3, 4], [3, 5]]
Ozooxo
quelle
28

Ich war auch daran interessiert und habe einen kleinen Leistungsvergleich durchgeführt, vielleicht etwas klarer als in @ senderles Antwort.

Für zwei Arrays (der klassische Fall):

Geben Sie hier die Bildbeschreibung ein

Für vier Arrays:

Geben Sie hier die Bildbeschreibung ein

(Beachten Sie, dass die Länge der Arrays hier nur ein paar Dutzend Einträge beträgt.)


Code zur Reproduktion der Diagramme:

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,
)
Nico Schlömer
quelle
17

Aufbauend auf @ senderles beispielhafter Basisarbeit habe ich zwei Versionen entwickelt - eine für C und eine für Fortran-Layouts - die oft etwas schneller sind.

  • cartesian_product_transpose_ppist - im Gegensatz zu @ senderle's cartesian_product_transpose, das eine völlig andere Strategie verwendet - eine Version davon cartesion_product, die das günstigere Transponierungsspeicherlayout + einige sehr geringfügige Optimierungen verwendet.
  • cartesian_product_ppbleibt beim ursprünglichen Speicherlayout. Was es schnell macht, ist die Verwendung von zusammenhängendem Kopieren. Aneinandergrenzende Kopien sind so viel schneller, dass das Kopieren eines vollständigen Speicherblocks, obwohl nur ein Teil davon gültige Daten enthält, dem Kopieren nur der gültigen Bits vorzuziehen ist.

Einige Perfplots. Ich habe separate Layouts für C- und Fortran-Layouts erstellt, da dies IMO unterschiedliche Aufgaben sind.

Namen, die mit 'pp' enden, sind meine Ansätze.

1) viele winzige Faktoren (jeweils 2 Elemente)

Geben Sie hier die Bildbeschreibung einGeben Sie hier die Bildbeschreibung ein

2) viele kleine Faktoren (jeweils 4 Elemente)

Geben Sie hier die Bildbeschreibung einGeben Sie hier die Bildbeschreibung ein

3) drei Faktoren gleicher Länge

Geben Sie hier die Bildbeschreibung einGeben Sie hier die Bildbeschreibung ein

4) zwei Faktoren gleicher Länge

Geben Sie hier die Bildbeschreibung einGeben Sie hier die Bildbeschreibung ein

Code (für jeden Plot müssen separate Läufe durchgeführt werden, da ich nicht herausfinden konnte, wie er zurückgesetzt werden soll; außerdem muss er entsprechend bearbeitet / kommentiert werden):

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
    )
Paul Panzer
quelle
Vielen Dank, dass Sie diese hervorragende Antwort geteilt haben. Wenn die Größe von arraysin cartesian_product_transpose_pp (Arrays) eine bestimmte Größe überschreitet, MemoryErrortritt dies auf. In dieser Situation möchte ich, dass diese Funktion kleinere Ergebnisblöcke liefert. Ich habe eine Frage zu diesem Thema gestellt. Können Sie meine Frage beantworten? Vielen Dank.
Sun Bear
13

Seit Oktober 2017 verfügt numpy über eine generische np.stackFunktion, die einen Achsenparameter übernimmt. Mit ihm können wir ein "verallgemeinertes kartesisches Produkt" unter Verwendung der "dstack and meshgrid" -Technik haben:

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

Hinweis zum axis=-1Parameter. Dies ist die letzte (innerste) Achse im Ergebnis. Es ist gleichbedeutend mit der Verwendung axis=ndim.

Ein weiterer Kommentar: Da kartesische Produkte sehr schnell in die Luft jagen, sollten wir die Werte im laufenden Betrieb verwenden und verwenden , es sei denn, wir müssen das Array aus irgendeinem Grund im Speicher realisieren. Wenn das Produkt sehr groß ist itertools.

MrDrFenner
quelle
8

Ich habe eine Weile @kennytm answer verwendet , aber als ich versuchte, dasselbe in TensorFlow zu tun, stellte ich fest, dass TensorFlow kein Äquivalent zu hat numpy.repeat(). Nach ein wenig Experimentieren habe ich eine allgemeinere Lösung für beliebige Punktvektoren gefunden.

Für 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)

und für 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)
Sean McVeigh
quelle
6

Das Scikit- Lernpaket bietet eine schnelle Implementierung genau dieser:

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

Beachten Sie, dass sich die Konvention dieser Implementierung von der gewünschten unterscheidet, wenn Sie sich um die Reihenfolge der Ausgabe kümmern. Für Ihre genaue Bestellung können Sie tun

product = cartesian((y,x))[:, ::-1]
jmd_dk
quelle
Ist das schneller als die Funktion von @ senderle?
CS95
@ cᴏʟᴅsᴘᴇᴇᴅ Ich habe nicht getestet. Ich hatte gehofft, dass dies zB in C oder Fortran implementiert wurde und somit ziemlich unschlagbar ist, aber es scheint mit NumPy geschrieben zu sein. Als solche ist diese Funktion praktisch, sollte aber nicht wesentlich schneller sein als das, was man mit NumPy-Konstrukten selbst konstruieren kann.
jmd_dk
4

Im Allgemeinen können Sie diese Methode verwenden, wenn Sie zwei 2d-Numpy-Arrays a und b haben und jede Zeile von a mit jeder Zeile von b verketten möchten (ein kartesisches Produkt aus Zeilen, ähnlich einem Join in einer Datenbank) ::

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))
Jonathan
quelle
3

Am schnellsten können Sie entweder einen Generatorausdruck mit der Kartenfunktion kombinieren:

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

Ausgaben (tatsächlich wird die gesamte resultierende Liste gedruckt):

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

oder mithilfe eines doppelten Generatorausdrucks:

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

Ausgaben (ganze Liste gedruckt):

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

Beachten Sie, dass der größte Teil der Rechenzeit in den Druckbefehl fließt. Die Generatorberechnungen sind ansonsten recht effizient. Ohne Druck sind die Berechnungszeiten:

execution time: 0.079208 s

für Generatorausdruck + Kartenfunktion und:

execution time: 0.007093 s

für den Doppelgeneratorausdruck.

Wenn Sie tatsächlich das tatsächliche Produkt jedes der Koordinatenpaare berechnen möchten, lösen Sie es am schnellsten als Numpy-Matrix-Produkt:

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

Ausgänge:

 [[     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

und ohne Drucken (in diesem Fall spart es nicht viel, da nur ein winziges Stück der Matrix tatsächlich ausgedruckt wird):

execution time: 0.003083 s
Mosegui
quelle
Für die Produktberechnung ist der äußere Produktrundfunk foo = a[:,None]*bschneller. Wenn Sie Ihre Timing-Methode ohne verwenden print(foo), sind es 0,001103 s gegenüber 0,002225 s. Mit timeit sind es 304 μs gegenüber 1,6 ms. Matrix ist bekanntermaßen langsamer als ndarray, daher habe ich Ihren Code mit np.array ausprobiert, aber er ist immer noch langsamer (1,57 ms) als Broadcasting.
Syockit
2

Dies kann auch einfach mit der Methode itertools.product erfolgen

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

Ergebnis: Array ([
[1, 4],
[1, 5],
[2, 4],
[2, 5],
[3, 4],
[3, 5]], dtype = int32)

Ausführungszeit: 0,000155 s

Muhammad Umar Amanat
quelle
1
Sie müssen nicht numpy anrufen. Damit funktionieren auch einfache alte Python-Arrays.
Coddy
0

In dem speziellen Fall, in dem Sie einfache Vorgänge ausführen müssen, z. B. das Hinzufügen für jedes Paar, können Sie eine zusätzliche Dimension einführen und den Rundfunk die Aufgabe erledigen lassen:

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

Ich bin mir nicht sicher, ob es einen ähnlichen Weg gibt, um die Paare tatsächlich selbst zu bekommen.

Caagr98
quelle
Wenn dtypeist floatkann man tun , (a[:, None, None] + 1j * b[None, :, None]).view(float)was überraschend schnell.
Paul Panzer