Geek asks and answers

Levinson's algorithm in python

I need to determine if many different greenhouse matrices are singular. I would like to be able to accurately calculate how many 12 by 12 0-1 Teplitz matrices are singularly, for example. Here is the 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, a very inefficient way to do this. In principle, Levinson Recursion is faster, but I don’t see how to implement it. Can anyone get me started (or is there an even faster and better way)?

+3
source share
1 answer

We need to speed up the calls toeplitzand det:

  • works in 2**klot size
  • create index first toeplitz
  • NumPy 1.8, det - ufunc, det .

:

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

:

%%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))

:

%%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())

:

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

:

np.allclose(r1, r2)

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())

:

Wall time: 494 ms

singcount 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

43892, 2,15 .

+3

More articles:


All Articles