2011-03-15 101 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