2011-03-15 67 views
2

对不起,标题含糊不清。我有两个相关的问题。首先让我们假设我有一个函数“hessian”,给出两个参数(x,y)返回一个矩阵。我现在想要计算在二维空间上运行的(x,y)矩阵。我想做类似的事情:numpy:在多维阵列上运行

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

与形状(100,100,2,2)的结果H.上述不起作用(ValueError:用序列设置数组元素)。我唯一想出来的是

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) 

我再次需要使用列表解析做迭代和在阵列中手动指定形状重新包装的结果。有没有更直接的方法来达到同样的效果?

谢谢!

编辑:保罗和JoshAdel的建议,我实现了一个版本,粗麻布的,与数组的作品,这里是

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这样的数组,则w变成W_ {k,l,i,j},其中numpy给出ValueError: too big

+3

是否有一个原因'hessian'不能进行修改,以接受阵列? – Paul 2011-03-15 03:44:55

+0

uhm,'hessian'对已经有三个指标的东西进行计算。是的,原则上可以重写它接受阵列我认为,只是...前面的大麻烦 – andreabedini 2011-03-15 04:01:00

+0

我敢打赌,如果你发布'hessian'正在做什么,我们有几个人可以帮你重写它以一种紧凑的方式直接处理阵列。 – JoshAdel 2011-03-15 18:13:13

回答

0

你可以尝试使用meshgrid,也许你有扁平化XN,YN:

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

问题在于返回类型,hessian返回一个2x2矩阵。我得到 “ValueError:使用序列设置数组元素。” – andreabedini 2011-03-24 00:38:10