Speed up random matrix computation

math, numpy, performance, python, scipy

Solution

One way is to save some work from repeated calling of toeplitz() function by caching the indexes where the values are being put. The following code is ~ 30% faster than the original code. The rest of the performance is in the rank calculation... And I don't know whether there exists a faster rank calculation for toeplitz matrices with 0s and 1s.

(update) The code is actually ~ 4 times faster if you replace matrix_rank by scipy.linalg.det() == 0 (determinant is faster then rank calculation for small matrices)

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

Problem

Original source