मानक विचलन

Question

के लिए NumPy फ़ंक्शन की मेमोरी खपत मैं वर्तमान में काफी बड़े रास्टर डेटा सेट (> 4 जीबी) पर काम करने के लिए जीडीएएल की पायथन बाइंडिंग का उपयोग कर रहा हूं। चूंकि उन्हें एक बार में स्मृति में लोड करना मेरे लिए कोई व्यवहार्य समाधान नहीं है, इसलिए मैं उन्हें छोटे ब्लॉक में पढ़ता हूं और टुकड़े टुकड़े टुकड़े टुकड़े करता हूं। प्रत्येक ब्लॉक पढ़ने के लिए एक नए आवंटन से बचने के लिए मैं buf_obj तर्क (here) का उपयोग कर रहा हूं ताकि मूल्यों को एक पूर्ववर्ती NumPy सरणी में पढ़ा जा सके। एक बिंदु पर मुझे पूरे रास्टर के औसत और मानक विचलन की गणना करना है। स्वाभाविक रूप से मैंने गणना के लिए np.std का उपयोग किया है। हालांकि मेरे कार्यक्रम की स्मृति खपत को प्रोफाइल करके मुझे एहसास हुआ कि np.std के प्रत्येक आमंत्रण के साथ अतिरिक्त स्मृति आवंटित और जारी की गई है।

एक न्यूनतम काम कर उदाहरण जो इस व्यवहार को दर्शाता है:

In [1] import numpy as np 
In [2] a = np.random.rand(20e6) # Approx. 150 MiB of memory 
In [3] %memit np.mean(a) 
peak memory: 187.30 MiB, increment: 0.48 MiB 
In [4] %memit np.std(a) 
peak memory: 340.24 MiB, increment: 152.91 MiB 

NumPy GitHub पर के स्रोत वृक्ष के भीतर तलाशी में पता चला है कि np.std समारोह आंतरिक _methods.py (here) से _var समारोह का आह्वान। एक बिंदु पर _var माध्य से विचलन की गणना करता है और उन्हें बताता है। इसलिए इनपुट सरणी की एक अस्थायी प्रति बनाई गई है। इस प्रकार समारोह अनिवार्य रूप से मानक विचलन की गणना करता है:

mu = sum(arr)/len(arr) 
tmp = arr - mu 
tmp = tmp * tmp 
sd = np.sum(tmp)/len(arr) 

हालांकि यह दृष्टिकोण छोटे इनपुट सरणियों के लिए ठीक है यह निश्चित रूप से बड़े लोगों के लिए जाने के लिए कोई रास्ता नहीं है। चूंकि मैं इस अतिरिक्त प्रतिलिपि से पहले बताई गई स्मृति के छोटे ब्लॉक का उपयोग कर रहा हूं, मेरे प्रोग्राम में स्मृति बिंदु से गेम ब्रेकिंग समस्या नहीं है। हालांकि मुझे क्या बग है कि प्रत्येक ब्लॉक के लिए अगला ब्लॉक पढ़ने से पहले एक नया आवंटन किया जाता है और जारी किया जाता है।

क्या न्यूमपी या साइपीपी के भीतर कोई अन्य फ़ंक्शन है जो वेल्फोर्ड एल्गोरिदम (Wikipedia) जैसे निरंतर स्मृति खपत के साथ औसत और मानक विचलन की एक पास गणना के लिए अपर्याप्तता का उपयोग करता है?

जाने का एक और तरीका _var फ़ंक्शन का एक कस्टम संस्करण लागू करने के लिए एक वैकल्पिक out एक पूर्ववर्ती बफर (जैसे NumPy ufuncs) के लिए तर्क होगा। इस दृष्टिकोण के साथ अतिरिक्त प्रति को समाप्त नहीं किया जाएगा, लेकिन कम से कम स्मृति खपत स्थिर रहेगी और प्रत्येक ब्लॉक में आवंटन के लिए रनटाइम बचाया जाएगा।

संपादित करें: वेज़ोस द्वारा सुझाए गए वेलफोर्ड एल्गोरिदम के साइथन कार्यान्वयन का परीक्षण किया गया।

Cython कार्यान्वयन (kezzos से संशोधित):

cimport cython 
cimport numpy as np 
from libc.math cimport sqrt 

@cython.boundscheck(False) 
def iterative_approach(np.ndarray[np.float32_t, ndim=1] a): 
    cdef long n = 0 
    cdef float mean = 0 
    cdef float M2 = 0 
    cdef long i 
    cdef float delta 
    cdef float a_min = 10000000 # Must be set to Inf and -Inf for real cases 
    cdef float a_max = -10000000 
    for i in range(len(a)): 
     n += 1 
     delta = a[i] - mean 
     mean += delta/n 
     M2 += delta * (a[i] - mean) 
     if a[i] < a_min: 
      a_min = a[i] 
     if a[i] > a_max: 
      a_max = a[i] 
    return a_min, a_max, mean, sqrt(M2/(n - 1)) 

NumPy कार्यान्वयन (मतलब और एसटीडी संभवतः एक समारोह में गणना की जा सकती):

def vector_approach(a): 
    return np.min(a), np.max(a), np.mean(a), np.std(a, ddof=1) 

टेस्ट एक यादृच्छिक डेटा सेट का उपयोग कर परिणाम (बार मिलीसेकंड में, 25 में से सर्वश्रेष्ठ):

---------------------------------- 
| Size | Iterative |  Vector | 
---------------------------------- 
| 1e2 | 0.00529 | 0.17149 | 
| 1e3 | 0.02027 | 0.16856 | 
| 1e4 | 0.17850 | 0.23069 | 
| 1e5 | 1.93980 | 0.77727 | 
| 1e6 | 18.78207 | 8.83245 | 
| 1e7 | 180.04069 | 101.14722 | 
| 1e8 | 1789.60228 | 1086.66737 | 
---------------------------------- 

ऐसा लगता है कि यह पुनरावर्तक है 10000+ तत्वों के साथ बड़े डेटा सेट के लिए छोटे डेटा सेट और न्यूमपी वेक्टर (संभवतः सिम त्वरित) दृष्टिकोण के साथ साइथन का उपयोग करके पिप्रोच तेज है। सभी परीक्षण पायथन 2.7.9 और न्यूपी संस्करण 1.9.2 के साथ किए गए थे।

ध्यान दें कि वास्तविक मामले में ऊपरी कार्यों में रास्टर के एक ब्लॉक के आंकड़ों की गणना करने के लिए उपयोग किया जाएगा।सभी ब्लॉक के लिए मानक विचलन और साधन विकिपीडिया (here) में सुझाई गई पद्धति के साथ संयुक्त किए जाने हैं। इसका लाभ यह है कि रास्टर के सभी तत्वों को सारांशित करने की आवश्यकता नहीं है और इस तरह फ्लोट ओवरफ्लो समस्या से बचा जाता है (कम से कम किसी बिंदु पर)।

Answer 1

बचाव के लिए साइथन!

a = np.random.rand(10000).astype(np.float) 
print std_welford(a) 
%timeit -n 10 -r 10 std_welford(a) 

Cython कोड

0.288327455521 
10 loops, best of 10: 59.6 µs per loop 

मूल कोड

0.289605617397 
10 loops, best of 10: 18.5 ms per loop 

: परीक्षण करने के लिए इस का उपयोग करते हुए

%%cython 
cimport cython 
cimport numpy as np 
from libc.math cimport sqrt 

@cython.boundscheck(False) 
def std_welford(np.ndarray[np.float64_t, ndim=1] a): 
    cdef int n = 0 
    cdef float mean = 0 
    cdef float M2 = 0 
    cdef int a_len = len(a) 
    cdef int i 
    cdef float delta 
    cdef float result 
    for i in range(a_len): 
     n += 1 
     delta = a[i] - mean 
     mean += delta/n 
     M2 += delta * (a[i] - mean) 
    if n < 2: 
     result = np.nan 
     return result 
    else: 
     result = sqrt(M2/(n - 1)) 
     return result 

: यह ऊपर एक अच्छा गति को प्राप्त होता है

Numpy

0.289493223504 
10 loops, best of 10: 29.3 µs per loop 

आसपास 300x की रफ्तार वृद्धि std तो। अभी भी numpy संस्करण के रूप में उतना अच्छा नहीं है ..

Answer 2

मुझे संदेह है कि आपको numpy में ऐसा कोई फ़ंक्शन मिलेगा। numpy का राजन डी'एट्रे यह है कि यह vector processor निर्देश सेट का लाभ उठाता है - बड़ी मात्रा में डेटा के समान निर्देश प्रदर्शन करता है। मूल रूप से numpy गति दक्षता के लिए स्मृति दक्षता व्यापार करता है। हालांकि, पाइथन की स्मृति गहन प्रकृति के कारण, numpy पूरी तरह से सरणी के साथ डेटा प्रकार को जोड़कर और प्रत्येक व्यक्तिगत तत्व को जोड़कर कुछ स्मृति क्षमता प्राप्त करने में सक्षम है।

गति में सुधार करने का एक तरीका है, लेकिन फिर भी कुछ मेमोरी ओवरहेड बलिदानों में मानक विचलन की गणना करता है।

import numpy as np 

def std(arr, blocksize=1000000): 
    """Written for py3, change range to xrange for py2. 
    This implementation requires the entire array in memory, but it shows how you can 
    calculate the standard deviation in a piecemeal way. 
    """ 
    num_blocks, remainder = divmod(len(arr), blocksize) 
    mean = arr.mean() 
    tmp = np.empty(blocksize, dtype=float) 
    total_squares = 0 
    for start in range(0, blocksize*num_blocks, blocksize): 
     # get a view of the data we want -- views do not "own" the data they point to 
     # -- they have minimal memory overhead 
     view = arr[start:start+blocksize] 
     # # inplace operations prevent a new array from being created 
     np.subtract(view, mean, out=tmp) 
     tmp *= tmp 
     total_squares += tmp.sum() 
    if remainder: 
     # len(arr) % blocksize != 0 and need process last part of array 
     # create copy of view, with the smallest amount of new memory allocation possible 
     # -- one more array *view* 
     view = arr[-remainder:] 
     tmp = tmp[-remainder:] 
     np.subtract(view, mean, out=tmp) 
     tmp *= tmp 
     total_squares += tmp.sum() 

    var = total_squares/len(arr) 
    sd = var ** 0.5 
    return sd 

a = np.arange(20e6) 
assert np.isclose(np.std(a), std(a)) 

गति को दिखा रहा है --- बड़ा blocksize, बड़ा गति को। और काफी कम मेमोरी ओवरहेड। पूरी तरह से कम स्मृति ओवरहेड 100% सटीक नहीं है।

In [70]: %timeit np.std(a) 
10 loops, best of 3: 105 ms per loop 

In [71]: %timeit std(a, blocksize=4096) 
10 loops, best of 3: 160 ms per loop 

In [72]: %timeit std(a, blocksize=1000000) 
10 loops, best of 3: 105 ms per loop 

In [73]: %memit std(a, blocksize=4096) 
peak memory: 360.11 MiB, increment: 0.00 MiB 

In [74]: %memit std(a, blocksize=1000000) 
peak memory: 360.11 MiB, increment: 0.00 MiB 

In [75]: %memit np.std(a) 
peak memory: 512.70 MiB, increment: 152.59 MiB