랜덤 행렬 계산 속도 향상

Question

랜덤 토플 리츠 행렬을 만들어 역전 가능성을 추정합니다. 현재 코드는

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을 가진 토플 리츠 행렬에 대한 더 빠른 순위 계산이 있는지 여부를 알 수 없습니다.

(갱신) 당신이 scipy.linalg.det에 의해 matrix_rank을 대체하면 코드가 실제로 ~ 4 배 빠른() == 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)] 

는 inefficiences의 번호를 가지고있다. 많은 목록을 불필요하게 작성, 확장 및 삭제합니다. 목록에서 random.choice()를 호출하여 실제로 하나의 임의 비트를 얻습니다.

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

그것은 원래의 코드를 먼저 입력 행렬의 LU 분해를 계산하여 선형 시스템을 해결하기 위해 LAPACK 루틴 dgesdd를 호출 같습니다 :이 같은 약 500 %까지 그들을주었습니다. det와 matrix_rank 장착

입력 행렬 (http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.det.html)만을 LU 분해를 계산 LAPACK의 dgetrf을 사용하여 결정을 계산한다.

따라서 matrix_rank과 det 호출의 점근 적 복잡성은 O (n^3)이며 이는 LU 분해의 복잡성입니다.

그러나 Toepelitz 시스템은 O (n^2)에서 해결할 수 있습니다 (Wikipedia에 따라). 따라서 코드를 큰 행렬에서 실행하려면 파이썬 확장을 작성하여 특수 라이브러리를 호출하는 것이 좋습니다.