यादृच्छिक मैट्रिक्स गणना

Question

मैं गतिशील Toeplitz matrices बना रहा हूं ताकि वे अनुमान लगा सकें कि वे अचूक हैं। मेरा वर्तमान कोड

import random 
from scipy.linalg import toeplitz 
import numpy as np 
for n in xrange(1,25): 
    rankzero = 0 
    for repeats in xrange(50000): 
     column = [random.choice([0,1]) for x in xrange(n)] 
     row = [column[0]]+[random.choice([0,1]) for x in xrange(n-1)] 
     matrix = toeplitz(column, row) 
     if (np.linalg.matrix_rank(matrix) < n): 
      rankzero += 1 
    print n, (rankzero*1.0)/50000 

क्या यह बढ़ाया जा सकता है?

मैं अधिक सटीकता प्राप्त करने के लिए 50000 मूल्य बढ़ाना चाहता हूं लेकिन वर्तमान में ऐसा करना बहुत धीमा है।

का उपयोग कर रूपरेखा केवल for n in xrange(10,14) से पता चलता

400000 9.482 0.000 9.482 0.000 {numpy.linalg.lapack_lite.dgesdd} 
    4400000 7.591 0.000 11.089 0.000 random.py:272(choice) 
    200000 6.836 0.000 10.903 0.000 index_tricks.py:144(__getitem__) 
     1 5.473 5.473 62.668 62.668 toeplitz.py:3(<module>) 
    800065 4.333 0.000 4.333 0.000 {numpy.core.multiarray.array} 
    200000 3.513 0.000 19.949 0.000 special_matrices.py:128(toeplitz) 
    200000 3.484 0.000 20.250 0.000 linalg.py:1194(svd) 
6401273/64.421 0.000 2.421 0.000 {len} 
    200000 2.252 0.000 26.047 0.000 linalg.py:1417(matrix_rank) 
    4400000 1.863 0.000 1.863 0.000 {method 'random' of '_random.Random' objects} 
    2201015 1.240 0.000 1.240 0.000 {isinstance} 
[...] 

Answer 1

एक तरह से अनुक्रमित जहां मूल्यों डाल जा रहा है कैशिंग द्वारा Toeplitz() फ़ंक्शन के बार-बार फोन करने से कुछ काम को बचाने के लिए है। निम्नलिखित कोड मूल कोड से ~ 30% तेज है। शेष प्रदर्शन रैंक गणना में है ... और मुझे नहीं पता कि 0s और 1s के साथ toeplitz matrices के लिए एक तेज रैंक गणना मौजूद है या नहीं।

(अपडेट) कोड वास्तव में ~ 4 गुना तेजी से अगर आप scipy.linalg.det द्वारा matrix_rank की जगह है() == 0 (निर्धारक तेजी से छोटे मैट्रिक्स के लिए पद गणना तो है)

import random 
from scipy.linalg import toeplitz, det 
import numpy as np,numpy.random 

class si: 
    #cache of info for toeplitz matrix construction 
    indx = None 
    l = None 

def xtoeplitz(c,r): 
    vals = np.concatenate((r[-1:0:-1], c)) 
    if si.indx is None or si.l != len(c): 
     a, b = np.ogrid[0:len(c), len(r) - 1:-1:-1] 
     si.indx = a + b 
     si.l = len(c) 
    # `indx` is a 2D array of indices into the 1D array `vals`, arranged so 
    # that `vals[indx]` is the Toeplitz matrix. 
    return vals[si.indx] 

def doit(): 
    for n in xrange(1,25): 
     rankzero = 0 
     si.indx=None 

     for repeats in xrange(5000): 

      column = np.random.randint(0,2,n) 
      #column=[random.choice([0,1]) for x in xrange(n)] # original code 

      row = np.r_[column[0], np.random.randint(0,2,n-1)] 
      #row=[column[0]]+[random.choice([0,1]) for x in xrange(n-1)] #origi 

      matrix = xtoeplitz(column, row) 
      #matrix=toeplitz(column,row) # original code 

      #if (np.linalg.matrix_rank(matrix) < n): # original code 
      if np.abs(det(matrix))<1e-4: # should be faster for small matrices 
       rankzero += 1 
     print n, (rankzero*1.0)/50000 

Answer 2

इन दो लाइनें जो 0s और 1s की सूचियां बनाती हैं:

column = [random.choice([0,1]) for x in xrange(n)] 
row = [column[0]]+[random.choice([0,1]) for x in xrange(n-1)] 

में कई अक्षमताएं हैं। वे अनावश्यक रूप से बहुत सी सूचियों का निर्माण, विस्तार और त्याग करते हैं, और वे वास्तव में केवल एक यादृच्छिक बिट प्राप्त करने के लिए यादृच्छिक.चॉइस() को एक सूची में कॉल करते हैं।

column = [0 for i in xrange(n)] 
row = [0 for i in xrange(n)] 

# NOTE: n must be less than 32 here, or remove int() and lose some speed 
cbits = int(random.getrandbits(n)) 
rbits = int(random.getrandbits(n)) 

for i in xrange(n): 
    column[i] = cbits & 1 
    cbits >>= 1 
    row[i] = rbits & 1 
    rbits >>= 1 

row[0] = column[0] 

Answer 3

यह अपने मूल कोड LAPACK दिनचर्या dgesdd बुला रहा है पहला इनपुट मैट्रिक्स की एक LU अपघटन कंप्यूटिंग द्वारा एक रेखीय प्रणाली को हल करने की तरह दिखता है: मैं इस तरह के बारे में 500% से उन्हें तेज।

det साथ matrix_rank की जगह जो इनपुट मैट्रिक्स (http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.det.html) का केवल LU अपघटन गणना करता निर्धारक LAPACK के dgetrf का उपयोग कर, गणना करता है।

matrix_rank और det कॉल दोनों की एसिम्प्टोटिक जटिलता इसलिए ओ (एन^3), ल्यू अपघटन की जटिलता है।

टोपेलिट्ज सिस्टम, हालांकि, ओ (एन^2) (विकिपीडिया के अनुसार) में हल किया जा सकता है। इसलिए, यदि आप अपने कोड को बड़े मैट्रिस पर चलाना चाहते हैं, तो एक विशेष लाइब्रेरी को कॉल करने के लिए एक पायथन एक्सटेंशन लिखना समझदारी होगी।