गैर तुच्छ numpy में temporaries बिना बाहरी उत्पादों की रकम

Question

वास्तविक समस्या मैं हल करना चाहते हैं, एन इकाई वैक्टर का एक सेट और एम वैक्टर का एक और सेट इकाई से प्रत्येक के लिए गणना दिया जाता है की औसत वैक्टर एम वैक्टरों में से प्रत्येक के साथ डॉट उत्पाद का पूर्ण मूल्य। अनिवार्य रूप से यह दो matrices के बाहरी उत्पाद की गणना कर रहा है और बीच में फंसे एक पूर्ण मूल्य के साथ औसत और औसत।

के लिए एन और एम नहीं बहुत बड़ी इस कठिन नहीं है और वहाँ आगे बढ़ने के लिए कई मायनों (नीचे देखें) कर रहे हैं। समस्या यह है कि एन और एम बड़े पैमाने पर बनाए गए अस्थायी हैं और प्रदान किए गए दृष्टिकोण के लिए व्यावहारिक सीमा प्रदान करते हैं। क्या यह गणना अस्थायी बनाने के बिना किया जा सकता है? मेरे पास मुख्य कठिनाई पूर्ण मूल्य की उपस्थिति के कारण है। क्या ऐसी गणनाओं को "थ्रेडिंग" के लिए सामान्य तकनीकें हैं?

एक उदाहरण के रूप में निम्नलिखित कोड

N = 7 
M = 5 

# Create the unit vectors, just so we have some examples, 
# this is not meant to be elegant 
phi = np.random.rand(N)*2*np.pi 
ctheta = np.random.rand(N)*2 - 1 
stheta = np.sqrt(1-ctheta**2) 
nhat = np.array([stheta*np.cos(phi), stheta*np.sin(phi), ctheta]).T 

# Create the other vectors 
m = np.random.rand(M,3) 

# Calculate the quantity we desire, here using broadcasting. 
S = np.average(np.abs(np.sum(nhat*m[:,np.newaxis,:], axis=-1)), axis=0) 

यह बहुत अच्छा है पर विचार करें, एस अब लंबाई एन की एक सरणी है और वांछित परिणाम होता है। दुर्भाग्यवश इस प्रक्रिया में हमने कुछ संभावित विशाल सरणी बनाई हैं।

np.sum(nhat*m[:,np.newaxis,:], axis=-1) 

का परिणाम एक एम एक्स एन सरणी है। अंतिम परिणाम, निश्चित रूप से, आकार एन का आकार है। एन और एम के आकारों को बढ़ाना शुरू करें और हम जल्दी से स्मृति त्रुटि में भाग लें।

जैसा कि ऊपर बताया, अगर निरपेक्ष मान आवश्यकता नहीं थी तो हम, के रूप में इस प्रकार अब einsum()

T = np.einsum('ik,jk,j', nhat, m, np.ones(M))/M 

का उपयोग कर यह काम करता है और काफी बड़ी एन और एम के लिए भी जल्दी से काम करता है आगे बढ़ सके। विशिष्ट समस्या के लिए मुझे abs() शामिल करने की आवश्यकता है, लेकिन एक अधिक सामान्य समाधान (शायद अधिक सामान्य ufunc) भी ब्याज का होगा।

Answer 1

कुछ टिप्पणियों के आधार पर ऐसा लगता है कि साइथन का उपयोग करने का सबसे अच्छा तरीका है। मैंने मूर्खतापूर्वक साइथन का उपयोग करने में कभी नहीं देखा है। यह कामकाजी कोड का उत्पादन करने के लिए अपेक्षाकृत आसान साबित होता है।

कुछ खोजों के बाद मैंने निम्नलिखित साइथन कोड को एक साथ रखा। यह नहीं सबसे सामान्य कोड है, शायद इसे लिखने का सबसे अच्छा तरीका नहीं है, और शायद इसे अधिक कुशल बनाया जा सकता है। फिर भी, मूल प्रश्न में einsum() कोड से केवल 25% धीमी है, इसलिए यह बहुत बुरा नहीं है! इसे मूल प्रश्न में किए गए सरणी के साथ स्पष्ट रूप से काम करने के लिए लिखा गया है (इसलिए इनपुट सरणी के अनुमानित मोड)।
चेतावनी के बावजूद यह मूल समस्या का एक उचित कुशल समाधान प्रदान करता है और इसी तरह की स्थितियों में शुरुआती बिंदु के रूप में कार्य कर सकता है।

import numpy as np 
cimport numpy as np 
import cython 
DTYPE = np.float64 
ctypedef np.float64_t DTYPE_t 
cdef inline double d_abs (double a) : return a if a >= 0 else -a 

@cython.boundscheck(False) 
@cython.wraparound(False) 
def process_vectors (np.ndarray[DTYPE_t, ndim=2, mode="fortran"] nhat not None, 
        np.ndarray[DTYPE_t, ndim=2, mode="c"] m not None) : 
    if nhat.shape[1] != m.shape[1] : 
     raise ValueError ("Arrays must contain vectors of the same dimension") 
    cdef Py_ssize_t imax = nhat.shape[0] 
    cdef Py_ssize_t jmax = m.shape[0] 
    cdef Py_ssize_t kmax = nhat.shape[1] # same as m.shape[1] 
    cdef np.ndarray[DTYPE_t, ndim=1] S = np.zeros(imax, dtype=DTYPE) 
    cdef Py_ssize_t i, j, k 
    cdef DTYPE_t val, tmp 
    for i in range(imax) : 
     val = 0 
     for j in range(jmax) : 
      tmp = 0 
      for k in range(kmax) : 
       tmp += nhat[i,k] * m[j,k] 
      val += d_abs(tmp) 
     S[i] = val/jmax 
    return S 

Answer 2

यह काफी धीमा है, लेकिन बड़े इंटरमीडिएट मैट्रिक्स को नहीं बनाता है।

vals = np.zeros(N) 
for i in xrange(N): 
    u = nhat[i] 
    for v in m: 
     vals[i]+=abs(np.dot(u,v)) 
    vals[i]=vals[i]/M 

संपादित करें: लूप के बाहर एम द्वारा विभाजित ले जाया गया।

संपादित 2: नया विचार, पोस्टरिटी और प्रासंगिक टिप्पणी के लिए पुराना रखना।

m2 = np.average(m,0) 
vals = np.zeros(N) 
for i in xrange(N): 
    u=nhat[i] 
    vals[i]=abs(np.dot(u,m2)) 

यह तेज़ है, लेकिन कभी-कभी अलग-अलग मूल्य देता है, मैं काम कर रहा हूं लेकिन शायद यह औसत समय में कैसे मदद कर सकता है।

संपादित करें 3: आह, यह पूर्ण मूल्य वस्तु है।हम्म

>>> S 
array([ 0.28620962, 0.65337876, 0.37470707, 0.46500913, 0.49579837, 
     0.29348924, 0.27444208, 0.74586928, 0.35789315, 0.3079964 , 
     0.298353 , 0.42571445, 0.32535728, 0.87505053, 0.25547394, 
     0.23964505, 0.44773271, 0.25235646, 0.4722281 , 0.33003338]) 
>>> vals 
array([ 0.2099343 , 0.6532155 , 0.33039334, 0.45366889, 0.48921527, 
     0.20467291, 0.16585856, 0.74586928, 0.31234917, 0.22198642, 
     0.21013519, 0.41422894, 0.26020981, 0.87505053, 0.1199069 , 
     0.06542492, 0.44145805, 0.08455833, 0.46824704, 0.28483342]) 

time to compute S: 0.000342130661011 seconds 
time to compute vals: 7.29560852051e-05 seconds 

संपादित 4: ठीक है, अगर आप अपने इकाई वैक्टर इस तेज चलाना चाहिए, मी में वैक्टर संभालने के लिए ज्यादातर सकारात्मक मूल्यों है हमेशा सकारात्मक है कि वे अपने डमी डेटा में कर रहे हैं।

m2 = np.average(m,0) 
vals = np.zeros(N) 
for i in xrange(N): 
    u=nhat[i] 
    if u[0] >= 0 and u[1] >= 0 and u[2] >= 0: 
     vals[i] = abs(np.dot(u,m2)) 
    else: 
     for j in xrange(M): 
      vals[i]+=abs(np.dot(u,m[j])) 
     vals[i]/=M 

Answer 3

मुझे नहीं लगता कि आपके सटीक संचालन को तेज करने के लिए कोई आसान तरीका (साइथन और इसी तरह के बाहर) है। लेकिन आप यह विचार करना चाहेंगे कि आपको वास्तव में गणना करने की आवश्यकता है कि आप क्या गणना कर रहे हैं। अगर शुद्ध मान का मतलब के बजाय आप root mean square इस्तेमाल कर सकते हैं, आप अब भी किसी भी तरह भीतरी उत्पादों के परिमाण औसत होगा, लेकिन आप के रूप में एक भी गोली में यह हो सकता है: यह रूप में ही है

rms = np.sqrt(np.einsum('ij,il,kj,kl,k->i', nhat, nhat, m, m, np.ones(M)/M)) 

कर:

rms_2 = np.sqrt(np.average(np.einsum('ij,kj->ik', nhat, m)**2, axis=-1)) 

हाँ, यह आप के लिए वास्तव में क्या नहीं कहा है, लेकिन मैं इसे के रूप में करीब है के रूप में आप एक vectorized दृष्टिकोण के साथ मिल जाएगा डर लग रहा है। यदि आप इस सड़क पर जाने का फैसला करते हैं, तो देखें कि np.einsum बड़े N और M के लिए कितना अच्छा प्रदर्शन करता है: इसमें बहुत से पैरामीटर और इंडेक्स पारित होने पर नीचे गिरने की प्रवृत्ति होती है।