numpy मैट्रिक्स प्रवंचना - उलटा बार मैट्रिक की राशि

Question

मैं निम्नलिखित करते हैं, और अभिसरण जब तक दोहराने की कोशिश कर रहा हूँ:

जहां प्रत्येक एक्स _मैंn x p है, और उनमें से r हैं r x n x p सरणी में samples कहा जाता है। Un x n है, Vp x p है। (मुझे matrix normal distribution का एमएलई मिल रहा है।) आकार सभी संभावित रूप से बड़े-आश हैं; मैं कम से कम r = 200, n = 1000, p = 1000 के आदेश पर चीजों की अपेक्षा कर रहा हूं।

मेरे वर्तमान कोड

V = np.einsum('aji,jk,akl->il', samples, np.linalg.inv(U)/(r*n), samples) 
U = np.einsum('aij,jk,alk->il', samples, np.linalg.inv(V)/(r*p), samples) 

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

मैं की तरह

V = sum(np.dot(x.T, scipy.linalg.solve(A, x)) for x in samples) 

(या इसी तरह की है कि PSD सत्ता शोषण कुछ) कुछ कर रही द्वारा उलटा से बचने के सकता है, लेकिन फिर वहाँ एक अजगर पाश है, और है कि numpy परियों रोना आता है।

मैं भी इस तरह से है कि मैं एक अजगर पाश करने के लिए बिना हर x के लिए solve का उपयोग कर A^-1 x की एक सरणी मिल सकता है में samples देगी कल्पना कर सकता है, लेकिन यह एक बड़ा सहायक सरणी स्मृति की बर्बादी है कि बनाता है।

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

(जैसा कि यह पता चला है, मुझे वास्तव में अंत में matrices U और V की आवश्यकता नहीं है - केवल उनके निर्धारक, या वास्तव में केवल उनके क्रोनकर उत्पाद के निर्धारक। इसलिए यदि किसी के पास कम काम करने के लिए एक चालाक विचार है और अभी भी निर्धारकों बाहर निकलना, कि बहुत सराहना की जाएगी।)

Answer 1

तक किसी को एक अधिक प्रेरित जवाब के साथ आता है, अगर मैं तुम थे, मैं करता हूँ परियों रोना ...

r, n, p = 200, 400, 400 

X = np.random.rand(r, n, p) 
U = np.random.rand(n, n) 

In [2]: %timeit np.sum(np.dot(x.T, np.linalg.solve(U, x)) for x in X) 
1 loops, best of 3: 9.43 s per loop 

In [3]: %timeit np.dot(X[0].T, np.linalg.solve(U, X[0])) 
10 loops, best of 3: 45.2 ms per loop 

तो एक पाइथन लूप होने और सभी परिणामों को एक साथ जोड़ना, 200 9 से अधिक बार 3 9 0 एमएस ले रहा है जो हल करने के लिए 200 सिस्टमों में से प्रत्येक को हल करने के लिए होता है। लूपिंग और संक्षेप मुक्त होने पर आपको 5% से भी कम सुधार मिलेगा। कुछ लोग एक पायथन फ़ंक्शन ओवरहेड भी कॉल कर सकते हैं, लेकिन संभवतः यह समीकरणों को हल करने के वास्तविक समय के खिलाफ नगण्य होने जा रहा है, इससे कोई फर्क नहीं पड़ता कि आप किस भाषा में कोड करते हैं।