Декартово произведение точек массива x и y на один массив точек 2D

У меня есть два массива, которые определяют оси X и Y сетки. Например:

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.transpose([numpy.tile(x, len(y)), numpy.repeat(y, len(x))])
array([[1, 4],
       [2, 4],
       [3, 4],
       [1, 5],
       [2, 5],
       [3, 5]])

См. Использование numpy для создания массива всех комбинаций из двух массивов для общего решения для вычисления декартового произведения N массивов.

Канонический cartesian_product(почти)

Есть много подходов к этой проблеме с различными свойствами. Некоторые из них быстрее, чем другие, а некоторые более общего назначения. После долгих испытаний и настроек я обнаружил, что следующая функция, которая вычисляет n-мерную cartesian_product, быстрее для большинства входных данных быстрее, чем большинство других. Для пары подходов, которые немного сложнее, но во многих случаях даже немного быстрее, см. Ответ Пола Панцера .

Учитывая этот ответ, это больше не самая быстрая реализация декартового произведения, 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 , который вдохновил меня попробовать использовать ix_этот способ, и unutbu , который предоставил несколько чрезвычайно полезных отзывов об этом ответе, включая предложение использовать numpy.result_type.

Известные альтернативы

Иногда быстрее записать непрерывные блоки памяти в порядке Фортрана. Это основа этой альтернативы, cartesian_product_transposeкоторая на некоторых аппаратных средствах оказалась быстрее, чем cartesian_product(см. Ниже). Однако ответ Пола Панцера, использующий тот же принцип, еще быстрее. Тем не менее, я включаю это здесь для заинтересованных читателей:

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

Поняв подход Panzer, я написал новую версию, которая почти такая же быстрая, как и его, и такая же простая, как 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, а также в ответе Нико Шлёмера .

Тесты против альтернатив

Вот набор тестов, которые показывают повышение производительности, которое предоставляют некоторые из этих функций относительно ряда альтернатив. Все тесты, показанные здесь, были выполнены на четырехъядерном компьютере под управлением Mac OS 10.12.5, Python 3.6.1 и numpy1.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остается конкурентоспособным, пока число входных массивов не превысит (примерно) четыре. После этого, 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для трансляции двух массивов вместе. Но 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изменить форму на этом этапе, мы могли бы передать выходные данные to и изменить форму впоследствии, что экономит некоторую работу:

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

Относительная эффективность этих двух подходов со временем изменилась. В более ранней версии 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) на более новой машине. Два подхода сейчас практически идентичны.

Старый тест

>>> 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Функция , определенная в другом ответе используется для выполнения довольно хорошо для больших входов. (Это то же самое, что и функция называется cartesian_product_recursiveвыше.) Для того , чтобы сравнить cartesianс dstack_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определено в самом начале этого ответа, что само по себе немного медленнее, чем самые быстрые реализации среди ответов на этот вопрос.

Вы можете просто сделать обычное понимание списка в 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]]

Я также заинтересовался этим и провел небольшое сравнение производительности, возможно, несколько понятнее, чем в ответе @ senderle.

Для двух массивов (классический случай):

Для четырех массивов:

(Обратите внимание, что длина массивов здесь составляет всего несколько десятков записей.)

Код для воспроизведения сюжетов:

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

Основываясь на образцовой фундаментальной работе @ senderle, я придумал две версии - одну для C и одну для разметки Fortran - которые часто бывают немного быстрее.

• cartesian_product_transpose_ppэто - в отличие от @ senderle's, cartesian_product_transposeкоторый использует совсем другую стратегию - версияcartesion_product которой используется более выгодная схема транспонирования памяти + некоторые очень незначительные оптимизации.
• cartesian_product_ppпридерживается оригинального расположения памяти. Что делает его быстрым, так это использование непрерывного копирования. Непрерывные копии оказываются намного быстрее, чем копирование полного блока памяти, хотя только часть его содержит действительные данные, предпочтительнее, чем копирование действительных битов.

Несколько перфлотов. Я сделал отдельные для макетов C и Fortran, потому что это разные задачи IMO.

Имена, оканчивающиеся на 'pp' - мои подходы.

1) много крошечных факторов (2 элемента каждый)

2) много мелких факторов (4 элемента каждый)

3) три фактора равной длины

4) два фактора равной длины

Код (нужно делать отдельные прогоны для каждого сюжета, потому что я не могу понять, как выполнить сброс; также нужно правильно редактировать / комментировать / выводить):

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
    )

По состоянию на октябрь 2017 года numpy теперь имеет обобщенную np.stackфункцию, которая принимает параметр оси. Используя его, мы можем получить «обобщенное декартово произведение», используя технику «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)

Обратите внимание на axis=-1параметр. Это последняя (самая внутренняя) ось в результате. Это эквивалентно использованию axis=ndim.

Еще один комментарий: поскольку декартовы произведения взрываются очень быстро, если нам не нужно по какой-то причине реализовать массив в памяти, если продукт очень большой, мы можем захотеть использовать itertoolsи использовать значения на лету.

Некоторое время я использовал ответ @kennytm , но при попытке сделать то же самое в TensorFlow, но обнаружил, что TensorFlow не имеет эквивалентаnumpy.repeat() . После небольшого эксперимента, я думаю, я нашел более общее решение для произвольных векторов точек.

Для 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)

и для 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)

Пакет Scikit-learn имеет быструю реализацию именно этого:

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

Обратите внимание, что соглашение этой реализации отличается от того, что вы хотите, если вы заботитесь о порядке вывода. Для вашего точного заказа вы можете сделать

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

В более общем случае, если у вас есть два двумерных числовых массива a и b, и вы хотите объединить каждую строку a с каждой строкой b (декартово произведение строк, вроде соединения в базе данных), вы можете использовать этот метод :

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

Самое быстрое, что вы можете получить, это либо объединить выражение генератора с функцией 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()))

Выходы (фактически весь итоговый список печатается):

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

или с помощью выражения двойного генератора:

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

Выходы (весь список напечатан):

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

Учтите, что большая часть времени вычислений уходит на команду печати. Расчеты генератора в остальном прилично эффективны. Без печати время расчета составляет:

execution time: 0.079208 s

для выражения генератора + функция карты и:

execution time: 0.007093 s

для выражения двойного генератора.

Если на самом деле вы хотите вычислить фактическое произведение каждой из пар координат, самое быстрое - это решить его как произведение матрицы на пустые места:

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

Выходы:

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

и без печати (в данном случае это не сильно экономит, так как на самом деле распечатывается только крошечный кусочек матрицы):

execution time: 0.003083 s

Это также легко сделать с помощью метода 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')

Результат: массив ([
[1, 4],
[1, 5],
[2, 4],
[2, 5],
[3, 4],
[3, 5]], dtype = int32)

Время выполнения: 0,000155 с

В конкретном случае, когда вам нужно выполнить простые операции, такие как сложение для каждой пары, вы можете ввести дополнительное измерение и позволить вещанию делать свою работу:

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

Я не уверен, есть ли подобный способ на самом деле получить сами пары.