Beschleunigung zufällige Matrix Berechnung

Question

Ich erstelle zufällige Toeplitz-Matrizen, um die Wahrscheinlichkeit zu schätzen, dass sie invertierbar sind. Mein aktueller Code ist

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 

Kann dies beschleunigt werden?

Ich möchte den Wert 50000 erhöhen, um mehr Genauigkeit zu erhalten, aber es ist derzeit zu langsam, um dies zu tun.

mit Profil nur for n in xrange(10,14) zeigt

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

Eine Möglichkeit ist durch das Zwischenspeichern der Indizes einige Arbeit von wiederholten Aufrufen von toeplitz() Funktion zu speichern, wo die Werte gesetzt werden. Der folgende Code ist ~ 30% schneller als der ursprüngliche Code. Der Rest der Leistung ist in der Rangberechnung ... Und ich weiß nicht, ob es eine schnellere Rangberechnung für Toeplitz-Matrizen mit 0s und 1s gibt.

(Update) Der Code ist eigentlich ~ 4 mal schneller, wenn Sie matrix_rank von scipy.linalg.det() == 0 (Determinante für kleine Matrizen dann Rang Berechnung schneller ist) ersetzen

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

Diese beiden Zeilen, die die Listen von 0s und 1s aufbauen:

haben eine Reihe von Ineffizienzen. Sie erstellen, erweitern und verwerfen viele Listen unnötigerweise und sie rufen random.choice() auf einer Liste auf, um nur ein zufälliges Bit zu erhalten. Ich eilte sie um etwa 500% wie folgt auf:

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

Es sieht aus wie Ihre Original-Code wird die lapack Routine dgesdd Aufruf ein lineares System, indem man zuerst das Berechnen einer LU-Zerlegung der Eingangsmatrix zu lösen.

Ersetzen matrix_rank mit det berechnet die Determinante lapack die Verwendung dgetrf, die nur die LU-Zerlegung der Eingangsmatrix (http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.det.html) berechnet. Die asymptotische Komplexität der beiden Aufrufe matrix_rank und det ist daher O (n^3), die Komplexität der LU-Zerlegung.

Toepelitz-Systeme können jedoch in O (n^2) gelöst werden (nach Wikipedia). Wenn Sie also Ihren Code auf großen Matrizen ausführen möchten, wäre es sinnvoll, eine Python-Erweiterung zu schreiben, um eine spezielle Bibliothek aufzurufen.