numpy - 使用 numpy 一维数组的置换副本构建二维数组的最快方法
numpy - fastest way to build 2d array with permuted copies of numpy 1d array
>>> import numpy as np
>>> a = np.arange(5)
>>> b = desired_function(a, 4)
array([[0, 3, 4, 1],
... [1, 2, 1, 3],
... [2, 4, 2, 4],
... [3, 1, 3, 0],
... [4, 0, 0, 2]])
到目前为止我尝试了什么
def repeat_and_shuffle(a, ncols):
nrows, = a.shape
m = np.tile(a.reshape(nrows, 1), (1, ncols))
return m
不知何故，我必须按列高效地随机播放 m[:,1:ncols]。
使用原始数组的随机排列构建新数组。
>>> a = np.arange(5)
>>> n = 4
>>> z = np.array([a]+[np.random.permutation(a) for _ in xrange(n-1)])
>>> z.T
array([[0, 0, 4, 3],
[1, 1, 3, 0],
[2, 3, 2, 4],
[3, 2, 0, 2],
[4, 4, 1, 1]])
>>>
由于随机性，可能出现重复的列。
这是创建此类数组的一种方法：
>>> a = np.arange(5)
>>> perms = np.argsort(np.random.rand(a.shape[0], 3), axis=0) # 3 columns
>>> np.hstack((a[:,np.newaxis], a[perms]))
array([[0, 3, 1, 4],
[1, 2, 3, 0],
[2, 1, 4, 1],
[3, 4, 0, 3],
[4, 0, 2, 2]])
这将创建所需形状的随机值数组，然后按对应值对每列中的索引进行排序。这个索引数组然后用于索引 a.
（使用 np.argsort 创建置换索引列数组的想法来自@jme 的回答 here。）
这是 Ashwini Chaudhary 解决方案的一个版本：
>>> a = numpy.array(['a', 'b', 'c', 'd', 'e'])
>>> a = numpy.tile(a[:,None], 5)
>>> a[:,1:] = numpy.apply_along_axis(numpy.random.permutation, 0, a[:,1:])
>>> a
array([['a', 'c', 'a', 'd', 'c'],
['b', 'd', 'b', 'e', 'a'],
['c', 'e', 'd', 'a', 'e'],
['d', 'a', 'e', 'b', 'd'],
['e', 'b', 'c', 'c', 'b']],
dtype='|S1')
我认为它构思周密并且在教学上很有用（我希望他取消删除它）。但有些令人惊讶的是，它始终是我执行过的测试中最慢的一个。定义：
>>> def column_perms_along(a, cols):
... a = numpy.tile(a[:,None], cols)
... a[:,1:] = numpy.apply_along_axis(numpy.random.permutation, 0, a[:,1:])
... return a
...
>>> def column_perms_argsort(a, cols):
... perms = np.argsort(np.random.rand(a.shape[0], cols - 1), axis=0)
... return np.hstack((a[:,None], a[perms]))
...
>>> def column_perms_lc(a, cols):
... z = np.array([a] + [np.random.permutation(a) for _ in xrange(cols - 1)])
... return z.T
...
对于小数组和几列：
>>> %timeit column_perms_along(a, 5)
1000 loops, best of 3: 272 µs per loop
>>> %timeit column_perms_argsort(a, 5)
10000 loops, best of 3: 23.7 µs per loop
>>> %timeit column_perms_lc(a, 5)
1000 loops, best of 3: 165 µs per loop
对于小数组和许多列：
>>> %timeit column_perms_along(a, 500)
100 loops, best of 3: 29.8 ms per loop
>>> %timeit column_perms_argsort(a, 500)
10000 loops, best of 3: 185 µs per loop
>>> %timeit column_perms_lc(a, 500)
100 loops, best of 3: 11.7 ms per loop
对于大数组和少数列：
>>> A = numpy.arange(1000)
>>> %timeit column_perms_along(A, 5)
1000 loops, best of 3: 2.97 ms per loop
>>> %timeit column_perms_argsort(A, 5)
1000 loops, best of 3: 447 µs per loop
>>> %timeit column_perms_lc(A, 5)
100 loops, best of 3: 2.27 ms per loop
对于大数组和许多列：
>>> %timeit column_perms_along(A, 500)
1 loops, best of 3: 281 ms per loop
>>> %timeit column_perms_argsort(A, 500)
10 loops, best of 3: 71.5 ms per loop
>>> %timeit column_perms_lc(A, 500)
1 loops, best of 3: 269 ms per loop
这个故事的寓意：总是测试！我想对于 非常 大数组，n log n 解决方案（如排序）的缺点在这里可能会变得很明显。但根据我的经验，numpy 的排序实现非常好。我敢打赌，您可以在注意到效果之前上升几个数量级。
假设您最终打算遍历多个一维输入数组，您可以缓存排列索引，然后在使用时仅 take 而不是 permute。即使一维数组的长度不同，这也可以工作：您只需要丢弃太大的排列索引。
实现的粗略（部分测试）代码：
def permute_multi(X, k, _cache={}):
"""For 1D input `X` of len `n`, it generates an `(k,n)` array
giving `k` permutations of `X`."""
n = len(X)
cached_inds = _cache.get('inds',np.array([[]]))
# make sure that cached_inds has shape >= (k,n)
if cached_inds.shape[1] < n:
_cache['inds'] = cached_inds = np.empty(shape=(k,n),dtype=int)
for i in xrange(k):
cached_inds[i,:] = np.random.permutation(n)
elif cached_inds.shape[0] < k:
pass # TODO: need to generate more rows
inds = cached_inds[:k,:] # dispose of excess rows
if n < cached_inds.shape[1]:
# dispose of high indices
inds = inds.compress(inds.ravel()<n).reshape((k,n))
return X[inds]
根据您的使用情况，您可能希望提供一些清除缓存的方法，或者至少提供一些启发式方法，以便在缓存的 n 和 k 变得比大多数缓存大得多时发现共同的输入。请注意，上面的函数给出的是 (k,n) 而不是 (n,k)，这是因为 numpy 默认行是连续的，我们希望 n 维度是连续的——如果你可以强制使用 Fortran 样式希望，或者只是转置输出（它翻转数组内的标志而不是真正移动数据）。
关于这个缓存概念是否在统计上有效，我相信在大多数情况下它可能没问题，因为它大致相当于将函数开始时的种子重置为一个固定常量......但是如果您正在对返回的数组做任何特别花哨的事情，您可能需要在使用这种方法之前仔细考虑。
快速基准测试表明，n=1000 和 k=1000 的（预热后）大约需要 2.2 ms，而完整 k 的 150 ms ]-循环 np.random.permutation。这大约快了 70 倍……但这是我们不调用 compress 的最简单的情况。对于 n=999 和 k=1000，用 n=1000 预热后，它需要额外几毫秒，给出 8ms 总时间，这仍然比 [= 快大约 19 倍14=]-循环。
>>> import numpy as np
>>> a = np.arange(5)
>>> b = desired_function(a, 4)
array([[0, 3, 4, 1],
... [1, 2, 1, 3],
... [2, 4, 2, 4],
... [3, 1, 3, 0],
... [4, 0, 0, 2]])
到目前为止我尝试了什么
def repeat_and_shuffle(a, ncols):
nrows, = a.shape
m = np.tile(a.reshape(nrows, 1), (1, ncols))
return m
不知何故，我必须按列高效地随机播放 m[:,1:ncols]。
使用原始数组的随机排列构建新数组。
>>> a = np.arange(5)
>>> n = 4
>>> z = np.array([a]+[np.random.permutation(a) for _ in xrange(n-1)])
>>> z.T
array([[0, 0, 4, 3],
[1, 1, 3, 0],
[2, 3, 2, 4],
[3, 2, 0, 2],
[4, 4, 1, 1]])
>>>
由于随机性，可能出现重复的列。
这是创建此类数组的一种方法：
>>> a = np.arange(5)
>>> perms = np.argsort(np.random.rand(a.shape[0], 3), axis=0) # 3 columns
>>> np.hstack((a[:,np.newaxis], a[perms]))
array([[0, 3, 1, 4],
[1, 2, 3, 0],
[2, 1, 4, 1],
[3, 4, 0, 3],
[4, 0, 2, 2]])
这将创建所需形状的随机值数组，然后按对应值对每列中的索引进行排序。这个索引数组然后用于索引 a.
（使用 np.argsort 创建置换索引列数组的想法来自@jme 的回答 here。）
这是 Ashwini Chaudhary 解决方案的一个版本：
>>> a = numpy.array(['a', 'b', 'c', 'd', 'e'])
>>> a = numpy.tile(a[:,None], 5)
>>> a[:,1:] = numpy.apply_along_axis(numpy.random.permutation, 0, a[:,1:])
>>> a
array([['a', 'c', 'a', 'd', 'c'],
['b', 'd', 'b', 'e', 'a'],
['c', 'e', 'd', 'a', 'e'],
['d', 'a', 'e', 'b', 'd'],
['e', 'b', 'c', 'c', 'b']],
dtype='|S1')
我认为它构思周密并且在教学上很有用（我希望他取消删除它）。但有些令人惊讶的是，它始终是我执行过的测试中最慢的一个。定义：
>>> def column_perms_along(a, cols):
... a = numpy.tile(a[:,None], cols)
... a[:,1:] = numpy.apply_along_axis(numpy.random.permutation, 0, a[:,1:])
... return a
...
>>> def column_perms_argsort(a, cols):
... perms = np.argsort(np.random.rand(a.shape[0], cols - 1), axis=0)
... return np.hstack((a[:,None], a[perms]))
...
>>> def column_perms_lc(a, cols):
... z = np.array([a] + [np.random.permutation(a) for _ in xrange(cols - 1)])
... return z.T
...
对于小数组和几列：
>>> %timeit column_perms_along(a, 5)
1000 loops, best of 3: 272 µs per loop
>>> %timeit column_perms_argsort(a, 5)
10000 loops, best of 3: 23.7 µs per loop
>>> %timeit column_perms_lc(a, 5)
1000 loops, best of 3: 165 µs per loop
对于小数组和许多列：
>>> %timeit column_perms_along(a, 500)
100 loops, best of 3: 29.8 ms per loop
>>> %timeit column_perms_argsort(a, 500)
10000 loops, best of 3: 185 µs per loop
>>> %timeit column_perms_lc(a, 500)
100 loops, best of 3: 11.7 ms per loop
对于大数组和少数列：
>>> A = numpy.arange(1000)
>>> %timeit column_perms_along(A, 5)
1000 loops, best of 3: 2.97 ms per loop
>>> %timeit column_perms_argsort(A, 5)
1000 loops, best of 3: 447 µs per loop
>>> %timeit column_perms_lc(A, 5)
100 loops, best of 3: 2.27 ms per loop
对于大数组和许多列：
>>> %timeit column_perms_along(A, 500)
1 loops, best of 3: 281 ms per loop
>>> %timeit column_perms_argsort(A, 500)
10 loops, best of 3: 71.5 ms per loop
>>> %timeit column_perms_lc(A, 500)
1 loops, best of 3: 269 ms per loop
这个故事的寓意：总是测试！我想对于 非常 大数组，n log n 解决方案（如排序）的缺点在这里可能会变得很明显。但根据我的经验，numpy 的排序实现非常好。我敢打赌，您可以在注意到效果之前上升几个数量级。
假设您最终打算遍历多个一维输入数组，您可以缓存排列索引，然后在使用时仅 take 而不是 permute。即使一维数组的长度不同，这也可以工作：您只需要丢弃太大的排列索引。
实现的粗略（部分测试）代码：
def permute_multi(X, k, _cache={}):
"""For 1D input `X` of len `n`, it generates an `(k,n)` array
giving `k` permutations of `X`."""
n = len(X)
cached_inds = _cache.get('inds',np.array([[]]))
# make sure that cached_inds has shape >= (k,n)
if cached_inds.shape[1] < n:
_cache['inds'] = cached_inds = np.empty(shape=(k,n),dtype=int)
for i in xrange(k):
cached_inds[i,:] = np.random.permutation(n)
elif cached_inds.shape[0] < k:
pass # TODO: need to generate more rows
inds = cached_inds[:k,:] # dispose of excess rows
if n < cached_inds.shape[1]:
# dispose of high indices
inds = inds.compress(inds.ravel()<n).reshape((k,n))
return X[inds]
根据您的使用情况，您可能希望提供一些清除缓存的方法，或者至少提供一些启发式方法，以便在缓存的 n 和 k 变得比大多数缓存大得多时发现共同的输入。请注意，上面的函数给出的是 (k,n) 而不是 (n,k)，这是因为 numpy 默认行是连续的，我们希望 n 维度是连续的——如果你可以强制使用 Fortran 样式希望，或者只是转置输出（它翻转数组内的标志而不是真正移动数据）。
关于这个缓存概念是否在统计上有效，我相信在大多数情况下它可能没问题，因为它大致相当于将函数开始时的种子重置为一个固定常量......但是如果您正在对返回的数组做任何特别花哨的事情，您可能需要在使用这种方法之前仔细考虑。
快速基准测试表明，n=1000 和 k=1000 的（预热后）大约需要 2.2 ms，而完整 k 的 150 ms ]-循环 np.random.permutation。这大约快了 70 倍……但这是我们不调用 compress 的最简单的情况。对于 n=999 和 k=1000，用 n=1000 预热后，它需要额外几毫秒，给出 8ms 总时间，这仍然比 [= 快大约 19 倍14=]-循环。