numpy：对多维数组的操作

Question

抱歉标题有点模糊。我有两个相关的问题。

首先，假设我有一个叫“hessian”的函数，它接受两个参数（x, y），并返回一个矩阵。现在我想在一个二维空间中计算这个矩阵，具体来说，我想做类似这样的事情：

x = linspace(1, 4, 100).reshape(-1,1)
y = linspace(1, 4, 100).reshape(1,-1)
H = vectorize(hessian)(x, y)

我希望得到的结果H的形状是(100,100,2,2)。但是这样做不行（出现了ValueError: setting an array element with a sequence）。我想到的唯一方法是

H = array([ hessian(xx, yy) for xx in x.flat for yy in y.flat ]).reshape(100,100,2,2)

有没有更好、更直接的方法呢？

第二个问题是，现在H的形状是(100,100,2,2)，而dominant_eigenvector(X)的功能正如你所想的那样。

U, V = hsplit(array(map(dominant_eigenvector, H.reshape(10000,2,2))), 2)

我又需要使用列表推导来进行迭代，并手动指定形状来重新打包结果。有没有更直接的方法可以实现同样的结果呢？

谢谢！

补充：根据Paul和JoshAdel的建议，我实现了一个可以处理数组的hessian版本，下面是它的代码

def hessian(w1, w2):
    w1 = atleast_1d(w1)[...,newaxis,newaxis]
    w2 = atleast_1d(w2)[...,newaxis,newaxis]
    o1, o2 = ix_(*map(xrange, Z.shape))
    W = Z * pow(w1, o1) * pow(w2, o2)
    A = (W).sum()
    B = (W * o1).sum()
    C = (W * o2).sum()
    D = (W * o1 * o1).sum()
    E = (W * o1 * o2).sum()
    F = (W * o2 * o2).sum()
    return array([[ D/A - B/A*B/A, E/A - B/A*C/A ],
                  [ E/A - B/A*C/A, F/A - C/A*C/A ]])

Z可以看作是一个大约250x150的全局数组。 o1和o2是Z的两个维度的索引，用来计算 像$\sum_{i,j} Z_{ij} * i * j$这样的东西。

这个版本的问题在于中间数组太大了。如果w1和w2是像w1_k和w2_l这样的数组， 那么W就变成了W_{k,l,i,j}，这时numpy会报错ValueError: too big。

Answer 1

你可以试试使用meshgrid，可能需要把xn和yn变成一维的。

x = linspace(1, 4, 100)
y = linspace(1, 4, 100)
xn,yn=meshgrid(x,y)
H = vectorize(hessian)(xn, yn)