2016-04-13 16 views
2

Sayısal einsum kullanıyorum, sütun vektörleri pts (Şekil 3, N) dizisinin nokta ürününü, bir matris dotps ile sonuçlanacak şekilde hesaplamak için kullanıyorum. N, N), tüm nokta ürünleri ile. Kullandığım kod:Üst üçgen elemanlar sadece NumPy einsum ile işleniyor

dotps = np.einsum('ij,ik->jk', pts, pts) 

Bu çalışır, ancak yalnızca ana köşegen üzerindeki değerlere ihtiyacım var. yani. sonucun çapraz üçgen olmayan üst üçgen kısmı. Sadece bu değerleri einsum ile hesaplamak mümkün mü? veya tüm matrisi hesaplamak için einsum kullanmaktan daha hızlı olan herhangi bir şekilde?

Pts dizimim oldukça büyük olabilir, bu yüzden yalnızca ihtiyacım olan değerleri hesaplayabilirsem, hesaplama hızımı iki katına çıkarırdım.

cevap

3

Alakalı sütunları dilim ve ardından kullanabilirsiniz np.einsum -

R,C = np.triu_indices(N,1) 
out = np.einsum('ij,ij->j',pts[:,R],pts[:,C]) 

Numune koşmak -

edelim zaman yaklaşımımızı ve herhangi olmadığını görmek -

In [109]: N = 5 
    ...: pts = np.random.rand(3,N) 
    ...: dotps = np.einsum('ij,ik->jk', pts, pts) 
    ...: 

In [110]: dotps 
Out[110]: 
array([[ 0.26529103, 0.30626052, 0.18373867, 0.13602931, 0.51162729], 
     [ 0.30626052, 0.56132272, 0.5938057 , 0.28750708, 0.9876753 ], 
     [ 0.18373867, 0.5938057 , 0.84699103, 0.35788749, 1.04483158], 
     [ 0.13602931, 0.28750708, 0.35788749, 0.18274288, 0.4612556 ], 
     [ 0.51162729, 0.9876753 , 1.04483158, 0.4612556 , 1.82723949]]) 

In [111]: R,C = np.triu_indices(N,1) 
    ...: out = np.einsum('ij,ij->j',pts[:,R],pts[:,C]) 
    ...: 

In [112]: out 
Out[112]: 
array([ 0.30626052, 0.18373867, 0.13602931, 0.51162729, 0.5938057 , 
     0.28750708, 0.9876753 , 0.35788749, 1.04483158, 0.4612556 ]) 

ayrıca optimize Performansı geliştirmek için kapsam. Biz np.einsum optimize etme konusunda çok yapabileceği gibi hafıza kısıtlamaları içinde kalmak

In [126]: N = 5000 

In [127]: pts = np.random.rand(3,N) 

In [128]: %timeit np.triu_indices(N,1) 
1 loops, best of 3: 413 ms per loop 

In [129]: R,C = np.triu_indices(N,1) 

In [130]: %timeit np.einsum('ij,ij->j',pts[:,R],pts[:,C]) 
1 loops, best of 3: 1.47 s per loop 

, bu görünmüyor. Yani, odağı np.triu_indices olarak değiştirelim.

N = 4 için elimizde:

In [131]: N = 4 

In [132]: np.triu_indices(N,1) 
Out[132]: (array([0, 0, 0, 1, 1, 2]), array([1, 2, 3, 2, 3, 3])) 

Bu normal bir deseni yaratmak tür olsa değişen bir gibi görünüyor. Bu, 3 ve 5 konumlarında vardiyaları olan kümülatif bir toplamla yazılabilir. jenerik düşünerek, biz böyle o şey kodlama sona ereceğini - Çeşitli N's için

def triu_indices_cumsum(N): 

    # Length of R and C index arrays 
    L = (N*(N-1))/2 

    # Positions along the R and C arrays that indicate 
    # shifting to the next row of the full array 
    shifts_idx = np.arange(2,N)[::-1].cumsum() 

    # Initialize "shift" arrays for finally leading to R and C 
    shifts1_arr = np.zeros(L,dtype=int) 
    shifts2_arr = np.ones(L,dtype=int) 

    # At shift positions along the shifts array set appropriate values, 
    # such that when cumulative summed would lead to desired R and C arrays. 
    shifts1_arr[shifts_idx] = 1 
    shifts2_arr[shifts_idx] = -np.arange(N-2)[::-1] 

    # Finall cumsum to give R, C 
    R_arr = shifts1_arr.cumsum() 
    C_arr = shifts2_arr.cumsum() 
    return R_arr, C_arr 

edelim sefer!

In [133]: N = 100 

In [134]: %timeit np.triu_indices(N,1) 
10000 loops, best of 3: 122 µs per loop 

In [135]: %timeit triu_indices_cumsum(N) 
10000 loops, best of 3: 61.7 µs per loop 

In [136]: N = 1000 

In [137]: %timeit np.triu_indices(N,1) 
100 loops, best of 3: 17 ms per loop 

In [138]: %timeit triu_indices_cumsum(N) 
100 loops, best of 3: 16.3 ms per loop 

Böylece, iyi N's için benziyor, triu_indices tabanlı özelleştirilmiş cumsum bakmaya değer olabilir!

+0

Teşekkürler. Sadece sorumu iyi ifade etmediğimi anladım. Lütfen düzenlemeye bakın. Temel olarak, ihtiyacım olan şey ana diyagonalin üstündeki değerler. yani. sonucun çapraz üçgen olmayan üst üçgen kısmı. – martinako

+0

@martinako Lütfen düzenlemelere göz atın. – Divakar

+0

İşte bu kadar! çok teşekkürler! – martinako