加速隨機矩陣計算

Question

我正在創建隨機Toeplitz矩陣來估計它們是可逆的概率。我目前的代碼是

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

一種方式是通過緩存其中的值被投入指標從託普利茨（）函數的調用重複節省一些工作。以下代碼比原始代碼快約30％。其餘的表現是排名計算... 而我不知道是否存在一個更快的秩和計算0和1的toeplitz矩陣。

（更新）的代碼是更快實際上〜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

這兩個構建0和1列表的行：

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

有許多無效。他們不必要地構建，擴展和丟棄大量列表，他們在列表上調用random.choice（）來獲取真正的一個隨機位。我加快他們上升了約500％這樣的：

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分解來求解線性方程。

更換matrix_rank與det計算使用LAPACK的dgetrf行列式，其計算僅輸入矩陣（http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.det.html）的LU分解。因此，matrix_rank和det調用的漸近複雜度因此是O（n^3），即LU分解的複雜度。但是，Toepelitz系統可以在O（n^2）中解決（根據維基百科）。所以，如果你想在大型矩陣上運行你的代碼，編寫一個python擴展來調用專門的庫是有意義的。