Tăng tốc độ tính toán ma trận ngẫu nhiên

Question

Tôi đang tạo các ma trận Toeplitz ngẫu nhiên để ước tính xác suất chúng không thể đảo ngược. Mã hiện tại của tôi là

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 

Điều này có thể được tăng tốc không?

Tôi muốn tăng giá trị 50000 để có được độ chính xác cao hơn nhưng hiện quá chậm.

Profiling chỉ sử dụng for n in xrange(10,14) lãm

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

Một cách là để tiết kiệm một số công việc từ gọi lặp đi lặp lại của Toeplitz() chức năng của bộ nhớ đệm các chỉ số mà các giá trị đang được đặt. Các mã sau đây là ~ 30% nhanh hơn so với mã ban đầu. Phần còn lại của hiệu suất là trong tính toán xếp hạng ... Và tôi không biết liệu có tồn tại một phép tính xếp hạng nhanh hơn cho ma trận toeplitz bằng 0 và 1 hay không.

(cập nhật) Mã này thực sự là ~ 4 lần nhanh hơn nếu bạn thay thế matrix_rank bởi scipy.linalg.det() == 0 (yếu tố quyết định là nhanh hơn sau đó tính toán xếp hạng cho các ma trận nhỏ)

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

Hai các dòng tạo danh sách 0 và 1:

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

có một số yếu tố không hiệu quả. Chúng xây dựng, mở rộng và loại bỏ rất nhiều danh sách một cách không cần thiết, và chúng gọi random.choice() trên một danh sách để có được những gì thực sự chỉ là một bit ngẫu nhiên. Tôi đẩy họ lên khoảng 500% như thế này:

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

Dường như mã ban đầu của bạn được gọi LAPACK thói quen dgesdd để giải quyết một hệ thống tuyến tính bằng cách đầu tiên tính toán một phân hủy LU của ma trận đầu vào.

Thay thế matrix_rank bằng det tính toán yếu tố quyết định bằng cách sử dụng phân đoạn LU của ma trận đầu vào (http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.det.html).

Độ phức tạp tiệm cận của cả hai cuộc gọi matrix_rank và det do đó là O (n^3), độ phức tạp của quá trình phân hủy LU. Tuy nhiên,

Các hệ thống Toepelitz có thể được giải quyết bằng O (n^2) (theo Wikipedia). Vì vậy, nếu bạn muốn chạy mã của bạn trên các ma trận lớn, nó sẽ có ý nghĩa để viết một phần mở rộng python để gọi một thư viện chuyên ngành.