Numpy: декартово произведение точек массива 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 и возвращение к массиву не самая эффективная форма.

 Alexey Lebedev21 июн. 2012 г., 21:09
Я заметил, что самым дорогим шагом в подходе itertools является окончательное преобразование из списка в массив. Без этого последнего шага он вдвое быстрее, чем пример Кена.
 coldspeed08 дек. 2018 г., 03:47
Этот вопрос был помечен как дубликат новой QnA как часть усилий по объединению всех знаний о слиянии на этом сайте. Информация об API слияния, соединения и конкатата фрагментирована и ее трудно найти из-за множества неисследуемых заголовков. С каноническим постом один родительский пост может быть связан с другими дубликатами и наоборот.

Ответы на вопрос(11)

Решение Вопроса
>>> 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 массивов.

 tlnagy27 июл. 2015 г., 20:35
Преимущество этого подхода в том, что он выдает согласованный вывод для массивов одинакового размера.meshgrid + dstack подход, хотя в некоторых случаях он быстрее, может привести к ошибкам, если вы ожидаете, что декартово произведение будет построено в том же порядке для массивов одинакового размера.
 senderle17 июл. 2017 г., 00:26
@ tlnagy, я не заметил ни одного случая, когда этот подход дает результаты, отличные от результатов, полученныхmeshgrid + dstack. Не могли бы вы опубликовать пример?
Канонический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, а также Нико Шлёмер S.

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

Вот набор тестов, которые показывают повышение производительности, которое предоставляют некоторые из этих функций относительно ряда альтернатив. Все тесты, показанные здесь, были выполнены на четырехъядерном компьютере под управлением Mac OS 10.12.5, Python 3.6.1 иnumpy 1.12.1. Известно, что различия в аппаратном и программном обеспечении дают разные результаты, поэтому YMMV. Запустите эти тесты для себя, чтобы быть уверенным!

Definitions:

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 иnumpy 1.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 и изменить его после, что экономит некоторую работу:

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

Вопреки утверждению в этот комментарий Я не видел доказательств того, что разные входы будут давать выходы разной формы, и, как показано выше, они делают очень похожие вещи, поэтому было бы довольно странно, если бы они это сделали. Пожалуйста, дайте мне знать, если вы найдете контрпример.

Testingmeshgrid + 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 для сравнения

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

 unutbu18 июл. 2017 г., 20:34
@ senderle: Вау, это мило! Кроме того, мне пришло в голову, что-то вродеdtype = np.find_common_type([arr.dtype for arr in arrays], []) можно использовать для нахождения общего dtype всех массивов вместо того, чтобы заставлять пользователя сначала размещать массив, управляющий dtype.
 unutbu17 июл. 2017 г., 19:49
Большое спасибо за ваши инновационные решения. Просто подумал, что вы хотели бы знать, что некоторые пользователи могут найтиcartesian_product_tranpose быстрее, чемcartesian_product в зависимости от операционной системы компьютера, версии Python или NumPy. Например, в Ubuntu 14.04, python3.4.3, numpy 1.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. Я также нахожуcartesian_product_transpose может быть быстрее, когдаlen(arrays) > 2.
 unutbu17 июл. 2017 г., 19:49
Кроме того,cartesian_product возвращает массив dtype с плавающей точкой, аcartesian_product_transpose возвращает массив того же типа d, что и первый (широковещательный) массив. Возможность сохранения dtype при работе с целочисленными массивами может быть причиной, по которой пользователи предпочитаютcartesian_product_transpose.
 user382099104 мар. 2015 г., 15:23
и добавлениеdtype=object вarr = np.empty( ) позволит использовать различные типы в продукте, например,arrays = [np.array([1,2,3]), ['str1', 'str2']].
 senderle17 июл. 2017 г., 20:59
@ unutbu еще раз спасибо - как я должен был знать, клонирование dtype не просто добавляет удобства; в некоторых случаях он ускоряет код еще на 20-30%.

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

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

Result: массив (
[1, 4],
[1, 5],
[2, 4],
[2, 5],
[3, 4],
[3, 5]], dtype = int32)

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

The Scikit учитьсяакет @ имеет быструю реализацию именно этого:

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

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

product = cartesian((y,x))[:, ::-1]
 coldspeed26 мар. 2018 г., 21:52
Это быстрее, чем функция @ senderle?
 jmd_dk26 мар. 2018 г., 22:04
@ cᴏʟᴅsᴘᴇᴇᴅ Я не проверял. Я надеялся, что это было реализовано, например, в С или Фортран и, таким образом, в значительной степени непобедимы, нопохоже на т будет написано с использованием NumPy. Как таковая, эта функция удобна, но не должна быть значительно быстрее, чем та, которую можно построить с помощью конструкций NumPy.

По состоянию на октябрь 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 и используйте значения на лету.

тельности, возможно, несколько яснее, чем в ответе @ 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: 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
    )

отве некоторое время, но при попытке сделать то же самое в 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)

если у вас есть два двумерных числовых массива 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()))

Outputs:

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

я придумал две версии - одну для C и одну для разметки Fortran - которые часто бывают немного быстрее.

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

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

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

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

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

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

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

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

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
    )

Ваш ответ на вопрос