numpy: operating on multidimensional arrays

Sorry for title's vagueness. I have two related questions.

First, let's say I have a function "hessian" that given two parameters (x, y) returns a matrix. I now want to compute that matrix for (x,y) running over a two dimensional space. I'd like to do something like:

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

with the resulting H of shape (100,100,2,2). The above doesn't work (ValueError: setting an array element with a sequence). The only thing I came up with is

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

is there a better, more direct, way ?

Second, now H has shape (100,100,2,2) and dominant_eigenvector(X) does exactly what you think.

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

I again need to use list comprehension to do the iteration and repack the result in an array specifying manually the shape. Is there a more direct way to achieve the same result ?

Thanks!

edit: as suggested by Paul and JoshAdel, I implemented a version of hessian that works with arrays, here it is

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 can be considered a global array of roughly 250x150. o1 and o2 index the two dimensions of Z to compute things like $\sum_{i,j} Z_{ij} * i * j$.

The problem with this version is that intermediate arrays are just too big. If w1 and w2 are arrays like w1_k w2_l W becomes W_{k,l,i,j} on which numpy gives ValueError: too big.

You could try to use meshgrid, maybe you have to flatten xn, yn:

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