CODEWITHSUNDEEP
pythonalgorithmmathnumpy
Simd
Simd

Reputation: 21343

Levinson algorithm in python

I need to determine if a lot of different Toeplitz matrices are singular. I would like to be able to count exactly how many 12 by 12 0-1 Toeplitz matrices are singular for example. Here is some code that does this.

import itertools
from scipy.linalg import toeplitz
import scipy.linalg
import numpy as np

n = 12
singcount = 0
for longtuple in itertools.product([0,1], repeat = 2*n-1):
    A = toeplitz(longtuple[0:n], longtuple[n-1:2*n-1])
    if (np.isclose(scipy.linalg.det(A),0)):
        singcount +=1
print singcount

However scipy.linalg.det is a very inefficient way to do this. In principle Levinson Recursion is faster but I cannot see how to implement it. Can anyone get me started (or is there an even faster and better way)?

Upvotes: 3

Views: 3145

Answers (1)

HYRY
HYRY

Reputation: 97331

We need to speedup toeplitz and det calls:

  • work in 2**k batch size
  • create a toeplitz index first
  • in NumPy 1.8, det is a general ufunc, which can calculate may det in one call.

Code:

import itertools
import numpy as np
from scipy.linalg import toeplitz, det

Here is the original code:

%%time
n = 12
todo = itertools.islice(itertools.product([0,1], repeat = 2*n-1), 0, 2**16)
r1 = []
for longtuple in todo:
    A = toeplitz(longtuple[0:n], longtuple[n-1:2*n-1])
    r1.append(det(A))

Here is the optimized code:

%%time
batch = 2**10
todo = itertools.islice(itertools.product([0,1], repeat = 2*n-1), 0, 2**16)
idx = toeplitz(range(n), range(n-1, 2*n-1))

r2 = []
while True:
    rows = list(itertools.islice(todo, 0, batch))
    if not rows:
        break
    rows_arr = np.array(rows)
    A = rows_arr[:, idx]
    r2.extend(np.linalg.det(A).tolist())

Here is the time result:

original: Wall time: 4.65 s
optimized: Wall time: 646 ms

We the check the result:

np.allclose(r1, r2)

You can increase the speed by unpackbits():

%%time
r3 = []
todo = np.arange(0, 2**16).astype(np.uint32).byteswap().view(np.uint8).reshape(-1, 4)
for i in range(todo.shape[0]//batch):
    B = np.unpackbits(todo[i*batch:(i+1)*batch], axis=-1)
    rows_arr = B[:, -23:]
    A = rows_arr[:, idx]
    r3.extend(np.linalg.det(A).tolist())

the time is:

Wall time: 494 ms

Here is the full code for the singcount for n=10:

%%time
count = 0
batch = 2**10
n = 10
n2 = 10*2-1
idx = toeplitz(range(n), range(n-1, 2*n-1))
todo = np.arange(0, 2**n2).astype(np.uint32).byteswap().view(np.uint8).reshape(-1, 4)
for i in range(todo.shape[0]//batch):
    B = np.unpackbits(todo[i*batch:(i+1)*batch], axis=-1)
    rows_arr = B[:, -n2:]
    A = rows_arr[:, idx]
    det = np.linalg.det(A)
    count += np.sum(np.isclose(det, 0))
print count

The output is 43892, and it took 2.15s on me PC.

Upvotes: 3

Related Questions