Przyspieszenie losowych obliczeń macierzy

Question

Tworzę losowe macierze Toeplitz, aby oszacować prawdopodobieństwo, że są odwracalne. Mój obecny kod to:

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 

Czy można to przyspieszyć?

Chciałbym zwiększyć wartość 50000, aby uzyskać większą dokładność, ale obecnie jest to zbyt wolne.

profilowania za pomocą tylko for n in xrange(10,14) pokazuje

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

Jednym ze sposobów jest, aby zaoszczędzić trochę pracy w wyniku wielokrotnego wywoływania funkcji Toeplitz() poprzez buforowanie indeksy gdzie wartości są wprowadzane. Poniższy kod jest ~ 30% szybszy niż oryginalny kod. Reszta wydajności jest w obliczeniach rangi ... I nie wiem, czy istnieje szybsze obliczanie rang dla macierzy toeplitz z 0s i 1s.

(aktualizacja) Kod jest rzeczywiście ~ 4 razy szybciej, jeśli zastąpić matrix_rank przez scipy.linalg.det() == 0 (wyznacznikiem jest szybszy następnie obliczeniu rang dla małych matrycach)

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

Te dwa linie tworzące listy 0 i 1:

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

mają wiele nieefektywności. Budują, rozwijają i odrzucają niepotrzebnie wiele list, a następnie wywołują random.choice() na liście, aby uzyskać naprawdę jeden losowy bit. I pędził je o około 500%, jak to:

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

wygląda jak oryginalny kod jest wywołanie Lapack rutynowe dgesdd rozwiązać układ liniowy najpierw wyliczania rozkładu LU macierzy wejściowej.

Wymiana matrix_rank z det oblicza determinantę pomocą LAPACK na dgetrf, który oblicza tylko rozkład logicznej matrycy wejściowej (http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.det.html).

Asymptotyczna złożoność obu wywołań matrix_rank i det to zatem O (n^3), złożoność dekompozycji LU.

Systemy Toepelitz można jednak rozwiązać w O (n^2) (według Wikipedii). Tak więc, jeśli chcesz uruchomić swój kod na dużych macierzach, byłoby sensownie napisać rozszerzenie python, aby zadzwonić do specjalistycznej biblioteki.