क्या numpy के साथ matrices की एक सरणी कुशलतापूर्वक उलटा करने का कोई तरीका है?

Question

आम तौर पर मैं नीचे दिए गए उदाहरण में for लूप में 3x3 matrices की सरणी को घुमा दूंगा। दुर्भाग्य से for लूप धीमे हैं। क्या ऐसा करने के लिए एक तेज, अधिक कुशल तरीका है?

import numpy as np 
A = np.random.rand(3,3,100) 
Ainv = np.zeros_like(A) 
for i in range(100): 
    Ainv[:,:,i] = np.linalg.inv(A[:,:,i]) 

Answer 1

यह पता चला है कि आप numpy.linalg कोड में दो स्तरों को जला रहे हैं। यदि आप numpy.linalg.inv को देखते हैं, तो आप देख सकते हैं कि यह सिर्फ numpy.linalg.solve (ए, inv (एशिप [0]) पर कॉल है। इसका आपके प्रत्येक पुनरावृत्ति में पहचान मैट्रिक्स को पुनर्निर्माण का प्रभाव है लूप के लिए। चूंकि आपके सभी सरणी एक ही आकार के हैं, यह समय बर्बाद है। पहचान मैट्रिक्स को पूर्व-आवंटित करके इस चरण को छोड़कर समय से 20% (fast_inverse) को रोकता है। मेरा परीक्षण बताता है कि सरणी आवंटित करना या आवंटित करना परिणामों की एक सूची से यह बहुत अंतर नहीं करता है।

एक स्तर को गहरा देखो और आपको लैपैक दिनचर्या पर कॉल मिलती है, लेकिन यह कई सैनिटी चेक में लपेटा गया है। यदि आप इन सब को बाहर निकाल देते हैं और बस लैपैक को कॉल करते हैं लूप के लिए (क्योंकि आप पहले से ही अपने मैट्रिक्स के आयामों को जानते हैं और शायद यह जानते हैं कि यह वास्तविक है, जटिल नहीं है), चीजें बहुत तेजी से चलती हैं (ध्यान दें कि मैंने अपना सरणी बड़ा कर दिया है):

import numpy as np 
A = np.random.rand(1000,3,3) 
def slow_inverse(A): 
    Ainv = np.zeros_like(A) 

    for i in range(A.shape[0]): 
     Ainv[i] = np.linalg.inv(A[i]) 
    return Ainv 

def fast_inverse(A): 
    identity = np.identity(A.shape[2], dtype=A.dtype) 
    Ainv = np.zeros_like(A) 

    for i in range(A.shape[0]): 
     Ainv[i] = np.linalg.solve(A[i], identity) 
    return Ainv 

def fast_inverse2(A): 
    identity = np.identity(A.shape[2], dtype=A.dtype) 

    return array([np.linalg.solve(x, identity) for x in A]) 

from numpy.linalg import lapack_lite 
lapack_routine = lapack_lite.dgesv 
# Looking one step deeper, we see that solve performs many sanity checks. 
# Stripping these, we have: 
def faster_inverse(A): 
    b = np.identity(A.shape[2], dtype=A.dtype) 

    n_eq = A.shape[1] 
    n_rhs = A.shape[2] 
    pivots = zeros(n_eq, np.intc) 
    identity = np.eye(n_eq) 
    def lapack_inverse(a): 
     b = np.copy(identity) 
     pivots = zeros(n_eq, np.intc) 
     results = lapack_lite.dgesv(n_eq, n_rhs, a, n_eq, pivots, b, n_eq, 0) 
     if results['info'] > 0: 
      raise LinAlgError('Singular matrix') 
     return b 

    return array([lapack_inverse(a) for a in A]) 


%timeit -n 20 aI11 = slow_inverse(A) 
%timeit -n 20 aI12 = fast_inverse(A) 
%timeit -n 20 aI13 = fast_inverse2(A) 
%timeit -n 20 aI14 = faster_inverse(A) 

परिणाम प्रभावशाली हैं:

20 loops, best of 3: 45.1 ms per loop 
20 loops, best of 3: 38.1 ms per loop 
20 loops, best of 3: 38.9 ms per loop 
20 loops, best of 3: 13.8 ms per loop 

संपादित करें: मैं क्या में हल लौटे हो जाता है पर बारीकी से पर्याप्त मुड़कर नहीं देखा। यह पता चला है कि 'बी' मैट्रिक्स ओवरराइट किया गया है और अंत में परिणाम शामिल है। यह कोड अब लगातार परिणाम देता है।

Answer 2

लूप वास्तव में विकल्पों की तुलना में बहुत धीमे नहीं हैं और इस मामले में, यह आपकी बहुत मदद नहीं करेगा।

import numpy as np 
A = np.random.rand(100,3,3) #this is to makes it 
          #possible to index 
          #the matrices as A[i] 
Ainv = np.array(map(np.linalg.inv, A)) 

समय इस समाधान अपने समाधान बनाम एक छोटी लेकिन उल्लेखनीय अंतर पैदावार:

# The for loop: 
100 loops, best of 3: 6.38 ms per loop 
# The map: 
100 loops, best of 3: 5.81 ms per loop 

मैं एक बनाने की उम्मीद के साथ numpy दिनचर्या 'vectorize' का उपयोग करने की कोशिश की लेकिन यहाँ एक सुझाव है यहां तक ​​कि क्लीनर समाधान भी है, लेकिन मुझे इसमें दूसरा नजर रखना होगा। सरणी ए में ऑर्डर करने का परिवर्तन शायद सबसे महत्वपूर्ण परिवर्तन है, क्योंकि यह इस तथ्य का उपयोग करता है कि numpy arrays को स्तंभ-वार का आदेश दिया जाता है और डेटा के रैखिक रीडआउट के लिए इस तरह से थोड़ा तेज़ तरीका होता है।