python中的Levinson算法

Question

我需要确定是否有很多不同的Toeplitz矩阵是奇异的。例如，我想精确地计算出有多少12乘12 0-1的Toeplitz矩阵是奇异的。这里有一些代码可以做到这一点。在

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

然而，scipy.linalg.det是一种非常低效的方法。原则上，{a1}更快，但我不知道如何实现它。有人能让我开始吗（或者有更快更好的方法吗）？在

Answer 1

我们需要加速toeplitz和det调用：

• 在2**k批量中工作
• 首先创建一个toeplitz索引
• 在numpy1.8中，det是一个通用的ufunc，它可以在一个调用中计算may-det。在

代码：

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

以下是原始代码：

^{pr2}$

下面是优化的代码：

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

以下是n=10的singcount的完整代码：

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