Velocizza calcolo a matrice casuale

Question

Sto creando matrici di Toeplitz casuali per stimare la probabilità che siano invertibili. Il mio codice attuale è

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 

Può essere accelerato?

Vorrei aumentare il valore 50000 per ottenere maggiore precisione ma è troppo lento per farlo al momento.

Profiling utilizzando solo for n in xrange(10,14) mostra

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

Un modo è quello di risparmiare un po 'di lavoro da ripetute chiamata della funzione Toeplitz() mettendo in cache gli indici in cui vengono messi i valori. Il seguente codice è ~ 30% più veloce del codice originale. Il resto della performance è nel calcolo del rango ... E non so se esiste un calcolo del rango più veloce per le matrici di toeplitz con 0 e 1 s.

(aggiornamento) Il codice è in realtà ~ 4 volte più veloce se si sostituisce matrix_rank da scipy.linalg.det() == 0 (determinante è più veloce quindi di calcolo rango per le piccole matrici)

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

Questi due linee che creano gli elenchi di 0 e 1:

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

hanno un numero di inefficienze. Costruiscono, espandono e scaricano molte liste inutilmente, e chiamano random.choice() in una lista per ottenere quello che è veramente un solo bit casuale. Li ho accelerato di circa il 500% in questo modo:

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

Sembra che il codice originale sta chiamando la routine di LAPACK dgesdd per risolvere un sistema lineare calcolando prima di decomposizione LU della matrice di input.

Sostituzione matrix_rank con det calcola il determinante usando LAPACK di dgetrf, che calcola solo la decomposizione LU della matrice di ingresso (http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.det.html).

La complessità asintotica di entrambe le matrix_rank e det chiamate sono quindi O (n^3), la complessità di decomposizione LU.

I sistemi di Toepelitz, tuttavia, possono essere risolti in O (n^2) (secondo Wikipedia). Quindi, se si desidera eseguire il codice su matrici di grandi dimensioni, sarebbe opportuno scrivere un'estensione python per chiamare una libreria specializzata.