创建一个元素为索引函数的numpy矩阵

Question

我该如何创建一个numpy矩阵，让它的元素根据它们的位置来决定呢？

比如说，一个乘法表：a[i,j] = i*j，这里的i和j就是矩阵的位置。

如果不使用numpy和python的常规方式，可以先创建一个全是零的数组，然后再用循环来填充。

不过，肯定有更好的方法可以做到这一点，不用循环。

更棒的是，能够直接创建这个矩阵，而不是先创建再填充。

Answer 1

关于乘法

np.multiply.outer(np.arange(5), np.arange(5))  # a_ij = i * j

还有一般情况

np.frompyfunc(
    lambda i, j: f(i, j), 2, 1
).outer(
    np.arange(5),
    np.arange(5),
).astype(np.float64)  # a_ij = f(i, j)

基本上，你是通过 np.frompyfunc 创建一个 np.ufunc，然后用索引进行外部操作。

编辑

不同解决方案之间的速度比较。

小矩阵：

Eyy![1]: %timeit np.multiply.outer(np.arange(5), np.arange(5))
100000 loops, best of 3: 4.97 µs per loop

Eyy![2]: %timeit np.array( [ [ i*j for j in xrange(5)] for i in xrange(5)] )
100000 loops, best of 3: 5.51 µs per loop

Eyy![3]: %timeit indices = np.indices((5, 5)); indices[0] * indices[1]
100000 loops, best of 3: 16.1 µs per loop

大矩阵：

Eyy![4]: %timeit np.multiply.outer(np.arange(4096), np.arange(4096))
10 loops, best of 3: 62.4 ms per loop

Eyy![5]: %timeit indices = np.indices((4096, 4096)); indices[0] * indices[1]
10 loops, best of 3: 165 ms per loop

Eyy![6]: %timeit np.array( [ [ i*j for j in xrange(4096)] for i in xrange(4096)] )
1 loops, best of 3: 1.39 s per loop

Answer 2

这里有一种方法可以做到这一点：

>>> indices = numpy.indices((5, 5))
>>> a = indices[0] * indices[1]
>>> a
array([[ 0,  0,  0,  0,  0],
       [ 0,  1,  2,  3,  4],
       [ 0,  2,  4,  6,  8],
       [ 0,  3,  6,  9, 12],
       [ 0,  4,  8, 12, 16]])

进一步解释一下，numpy.indices((5, 5)) 会生成两个数组，这两个数组分别包含一个 5x5 数组的 x 和 y 坐标，长得像这样：

>>> numpy.indices((5, 5))
array([[[0, 0, 0, 0, 0],
        [1, 1, 1, 1, 1],
        [2, 2, 2, 2, 2],
        [3, 3, 3, 3, 3],
        [4, 4, 4, 4, 4]],

       [[0, 1, 2, 3, 4],
        [0, 1, 2, 3, 4],
        [0, 1, 2, 3, 4],
        [0, 1, 2, 3, 4],
        [0, 1, 2, 3, 4]]])

当你把这两个数组相乘时，numpy 会把这两个数组在每个位置上的值相乘，然后返回结果。

Answer 3

一个通用的解决方案是使用 np.fromfunction()

根据文档：

numpy.fromfunction(function, shape, **kwargs)

这个函数通过在每个坐标上执行一个函数来构建一个数组。因此，结果数组在坐标 (x, y, z) 处的值是 fn(x, y, z)。

下面的代码片段应该能生成所需的矩阵。

import numpy as np

np.fromfunction(lambda i, j: i*j, (5,5))

输出：

array([[  0.,   0.,   0.,   0.,   0.],
       [  0.,   1.,   2.,   3.,   4.],
       [  0.,   2.,   4.,   6.,   8.],
       [  0.,   3.,   6.,   9.,  12.],
       [  0.,   4.,   8.,  12.,  16.]])

这个函数的第一个参数是一个可调用的函数，它会在每个坐标上执行。如果 foo 是你作为第一个参数传入的函数，那么 foo(i,j) 就会是坐标 (i,j) 处的值。这种情况在更高维度中也是适用的。你可以通过 shape 参数来修改坐标数组的形状。

编辑：

根据评论中提到的使用自定义函数，比如 lambda x,y: 2*x if x > y else y/2，以下代码可以正常工作：

import numpy as np

def generic_f(shape, elementwise_f):
    fv = np.vectorize(elementwise_f)
    return np.fromfunction(fv, shape)


def elementwise_f(x , y):
    return 2*x if x > y else y/2

print(generic_f( (5,5), elementwise_f))

输出：

[[0.  0.5 1.  1.5 2. ]
 [2.  0.5 1.  1.5 2. ]
 [4.  4.  1.  1.5 2. ]
 [6.  6.  6.  1.5 2. ]
 [8.  8.  8.  8.  2. ]]

用户需要传入一个标量函数，这个函数定义了逐元素的操作。 np.vectorize 用于将用户定义的标量函数向量化，并传递给 np.fromfunction()。