我有两个 numpy 数组，它们定义了网格的 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 数组并不是最有效的形式。
【问题讨论】：
标签： python numpy cartesian-product
cartesian_product（几乎）对于这个问题，有许多不同属性的方法。有些比其他更快，有些更通用。经过大量测试和调整后，我发现以下计算 n 维 cartesian_product 的函数对于许多输入来说比大多数其他函数更快。对于稍微复杂一些但在许多情况下甚至更快的方法，请参阅Paul Panzer 的答案。
鉴于这个答案，这不再是我所知道的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_ 的记录使用是 generate indices into 一个数组，恰好可以使用具有相同形状的数组进行广播分配。非常感谢mgilson，他启发了我尝试以这种方式使用ix_，并感谢unutbu，他为此答案提供了一些非常有用的反馈，包括使用numpy.result_type 的建议。
有时以 Fortran 顺序写入连续的内存块会更快。这是此替代方案cartesian_product_transpose 的基础，它在某些硬件上已被证明比cartesian_product 更快（见下文）。但是，使用相同原理的 Paul Panzer 的答案更快。尽管如此，我还是在这里为感兴趣的读者提供了这个：
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 的回答以及 Nico Schlömer 的回答。
以下是一组测试，这些测试显示了其中一些功能相对于许多替代方案提供的性能提升。此处显示的所有测试均在四核机器上执行，运行 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，然后再进行整形，这样可以节省一些工作：
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]])
与this comment 中的说法相反，我没有看到任何证据表明不同的输入会产生不同形状的输出，并且如上所示，它们做的事情非常相似，所以如果他们这样做会很奇怪。如果您找到反例，请告诉我。
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 上执行的测试以供比较。
another answer 中定义的 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 慢——这本身比这个问题的答案中最快的实现要慢一些.
【讨论】：
dtype=object 添加到arr = np.empty( ) 将允许在产品中使用不同的类型，例如arrays = [np.array([1,2,3]), ['str1', 'str2']].
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 时可能会更快。
cartesian_product 返回一个浮点 dtype 数组，而 cartesian_product_transpose 返回一个与第一个（广播）数组具有相同 dtype 的数组。使用整数数组时保留 dtype 的能力可能是用户青睐cartesian_product_transpose 的一个原因。
dtype = np.find_common_type([arr.dtype for arr in arrays], []) 这样的东西可以用来找到所有数组的共同dtype，而不是强迫用户首先放置控制dtype的数组。
>>> 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]])
有关计算 N 个数组的笛卡尔积的通用解决方案，请参阅 Using numpy to build an array of all combinations of two arrays。
【讨论】：
meshgrid + dstack 方法虽然在某些情况下更快，但如果您希望以相同顺序为相同大小的数组构造笛卡尔积，则可能会导致错误。
meshgrid + dstack 产生的结果不同的任何情况。可以举个例子吗？
你可以在 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]]
【讨论】：
[(x0, y0) for x0 in x for y0 in y] 做 [[(x0, y0) for x0 in x] for y0 in y]
我也对此感兴趣，并做了一些性能比较，可能比@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 的 cartesian_product_transpose 完全不同的策略不同 - cartesion_product 的一个版本使用更有利的转置内存布局 + 一些非常小的优化。cartesian_product_pp 坚持原来的内存布局。让它快速的原因是它使用连续复制。结果证明，连续副本的速度要快得多，即使仅部分包含有效数据，复制整个内存块也比仅复制有效位更可取。一些性能图。我为 C 和 Fortran 布局制作了单独的布局，因为这些是 IMO 不同的任务。
以“pp”结尾的名称是我的方法。
1) 许多微小的因素（每个 2 个元素）
2) 许多小因素（每个 4 个元素）
3) 三个等长因子
4) 两个等长的因子
代码（需要为每个绘图 b/c 单独运行我不知道如何重置；还需要适当地编辑/注释/注释）：
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
)
【讨论】：
arrays的大小超过一定大小时，会出现MemoryError。在这种情况下，我希望这个函数产生更小的结果块。我已经就此事发布了question。你能解决我的问题吗？谢谢。
截至 2017 年 10 月，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。
另外一个评论，因为笛卡尔积爆炸得很快，除非我们出于某种原因需要在内存中实现数组，如果积很大，我们可能要使用@ 987654325@ 并即时使用这些值。
【讨论】：
我用了一段时间@kennytm answer，但是当尝试在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]
【讨论】：
更一般地说，如果你有两个二维 numpy 数组 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
用于双生成器表达式。
如果你真正想要的是计算每个坐标对的实际乘积，最快的是将其求解为numpy矩阵乘积：
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
【讨论】：
foo = a[:,None]*b更快。在没有print(foo) 的情况下使用您的计时方法，它是 0.001103 秒与 0.002225 秒。使用 timeit，它是 304 μs 对 1.6 ms。众所周知，矩阵比 ndarray 慢，所以我用 np.array 尝试了你的代码，但它仍然比广播慢（1.57 毫秒）。
这也可以通过使用 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 s
【讨论】：
在需要对每一对进行加法等简单操作的特定情况下，可以引入一个额外的维度，让广播来完成这项工作：
>>> 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]])
我不确定是否有任何类似的方法可以实际获取这些对。
【讨论】：
dtype 是float，你可以使用(a[:, None, None] + 1j * b[None, :, None]).view(float)，这非常快。
我参加聚会有点晚了，但我遇到了这个问题的一个棘手变体。 假设我想要几个数组的笛卡尔积，但笛卡尔积最终比计算机的内存大得多（但是，使用该积完成的计算速度很快，或者至少是可并行的）。
显而易见的解决方案是将这个笛卡尔积分成块，并一个接一个地处理这些块（以一种“流式”方式）。您可以使用itertools.product 轻松做到这一点，但速度非常慢。此外，这里提出的解决方案（尽可能快）都没有给我们这种可能性。我提出的解决方案使用 Numba，比这里提到的“规范”cartesian_product 稍快。这很长，因为我试图尽可能地优化它。
import numba as nb
import numpy as np
from typing import List
@nb.njit(nb.types.Tuple((nb.int32[:, :],
nb.int32[:]))(nb.int32[:],
nb.int32[:],
nb.int64, nb.int64))
def cproduct(sizes: np.ndarray, current_tuple: np.ndarray, start_idx: int, end_idx: int):
"""Generates ids tuples from start_id to end_id"""
assert len(sizes) >= 2
assert start_idx < end_idx
tuples = np.zeros((end_idx - start_idx, len(sizes)), dtype=np.int32)
tuple_idx = 0
# stores the current combination
current_tuple = current_tuple.copy()
while tuple_idx < end_idx - start_idx:
tuples[tuple_idx] = current_tuple
current_tuple[0] += 1
# using a condition here instead of including this in the inner loop
# to gain a bit of speed: this is going to be tested each iteration,
# and starting a loop to have it end right away is a bit silly
if current_tuple[0] == sizes[0]:
# the reset to 0 and subsequent increment amount to carrying
# the number to the higher "power"
current_tuple[0] = 0
current_tuple[1] += 1
for i in range(1, len(sizes) - 1):
if current_tuple[i] == sizes[i]:
# same as before, but in a loop, since this is going
# to get called less often
current_tuple[i + 1] += 1
current_tuple[i] = 0
else:
break
tuple_idx += 1
return tuples, current_tuple
def chunked_cartesian_product_ids(sizes: List[int], chunk_size: int):
"""Just generates chunks of the cartesian product of the ids of each
input arrays (thus, we just need their sizes here, not the actual arrays)"""
prod = np.prod(sizes)
# putting the largest number at the front to more efficiently make use
# of the cproduct numba function
sizes = np.array(sizes, dtype=np.int32)
sorted_idx = np.argsort(sizes)[::-1]
sizes = sizes[sorted_idx]
if chunk_size > prod:
chunk_bounds = (np.array([0, prod])).astype(np.int64)
else:
num_chunks = np.maximum(np.ceil(prod / chunk_size), 2).astype(np.int32)
chunk_bounds = (np.arange(num_chunks + 1) * chunk_size).astype(np.int64)
chunk_bounds[-1] = prod
current_tuple = np.zeros(len(sizes), dtype=np.int32)
for start_idx, end_idx in zip(chunk_bounds[:-1], chunk_bounds[1:]):
tuples, current_tuple = cproduct(sizes, current_tuple, start_idx, end_idx)
# re-arrange columns to match the original order of the sizes list
# before yielding
yield tuples[:, np.argsort(sorted_idx)]
def chunked_cartesian_product(*arrays, chunk_size=2 ** 25):
"""Returns chunks of the full cartesian product, with arrays of shape
(chunk_size, n_arrays). The last chunk will obviously have the size of the
remainder"""
array_lengths = [len(array) for array in arrays]
for array_ids_chunk in chunked_cartesian_product_ids(array_lengths, chunk_size):
slices_lists = [arrays[i][array_ids_chunk[:, i]] for i in range(len(arrays))]
yield np.vstack(slices_lists).swapaxes(0,1)
def cartesian_product(*arrays):
"""Actual cartesian product, not chunked, still fast"""
total_prod = np.prod([len(array) for array in arrays])
return next(chunked_cartesian_product(*arrays, total_prod))
a = np.arange(0, 3)
b = np.arange(8, 10)
c = np.arange(13, 16)
for cartesian_tuples in chunked_cartesian_product(*[a, b, c], chunk_size=5):
print(cartesian_tuples)
这将以 5 个 3-uples 的块输出我们的笛卡尔积：
[[ 0 8 13]
[ 0 8 14]
[ 0 8 15]
[ 1 8 13]
[ 1 8 14]]
[[ 1 8 15]
[ 2 8 13]
[ 2 8 14]
[ 2 8 15]
[ 0 9 13]]
[[ 0 9 14]
[ 0 9 15]
[ 1 9 13]
[ 1 9 14]
[ 1 9 15]]
[[ 2 9 13]
[ 2 9 14]
[ 2 9 15]]
如果您愿意了解这里正在做什么，njitted 函数背后的直觉是在一个奇怪的数字基础中枚举每个“数字”，其元素将由输入数组的大小组成（而不是以常规二进制、十进制或十六进制为基数的相同数字）。
显然，这种解决方案对于大型产品很有趣。对于小型的，开销可能有点高。
注意：由于 numba 仍在大量开发中，我使用 numba 0.50 来运行它，使用 python 3.6。
【讨论】：
还有一个：
>>>x1, y1 = np.meshgrid(x, y)
>>>np.c_[x1.ravel(), y1.ravel()]
array([[1, 4],
[2, 4],
[3, 4],
[1, 5],
[2, 5],
[3, 5]])
【讨论】：
.meshgrid，我还使用了.ravel()，它将n维数组展平为一维k向量（k将是前一个数组的大小），然后它就来了.c_ 将其输入沿第二个轴粘合在一起（它可以连接两个以上的输入）。关于差异，我不太确定这是否比其他人更快，更差或其他任何东西，我只是喜欢简洁。
受Ashkan's answer启发，你也可以试试下面的。
>>> x, y = np.meshgrid(x, y)
>>> np.concatenate([x.flatten().reshape(-1,1), y.flatten().reshape(-1,1)], axis=1)
这将为您提供所需的笛卡尔积！
【讨论】：