加速随机矩阵计算

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扩展来调用专门的库是有意义的。