Berlekamp Massey算法

Published: by Creative Commons Licence

  • Tags:

Berlekamp-Massey算法

要了解Berlekamp-Massey算法的原理,可以看一下这篇论文

Berlekamp-Massey算法可以用于寻找线性递推数列的最短递推关系。所谓的线性递推数列是指这样的一个无穷数列a0,a1,,而线性递推关系是指一组系数c0,,ck1k为它的长度,满足对任意ik,都有ak=k1i=0ciaki1

实际上,如果我们记xi=ai,则可以发现最短递推关系实际上是一个最小多项式:

Za(x)=xkk1i=0cixki1

BM算法求解递最短线性递推关系

如果已知递推关系长度不超过n,那么我们需要递推数列的前面前面多少项才能计算出正确的递推关系,记这个数为m。BM算法的论文中指出之后如果由于最短递推关系不再满足而发生重构,那么重构后的递推关系不小于m+12,但是我们提过整个序列的最短递推关系不会超出n项,因此我们选择m=2n就足够了。

下面再讨论一种情况,在模p的意义下,要求我们为一个矩阵n×m的矩阵组成的序列A1,,Ak计算最短线性递推关系。我们可以随机一个n维的列向量um维的行向量v,之后找到序列uA1v,uA2v,,uAkv的最短线性递推关系,可以通过Schwartz–Zippel定理定理推出AiuAiv两个序列的最短线性递推关系不同的概率仅为n+mp

BM算法求解递推数列第n

我们要计算的是an,由于an=xn,而Za(x)是它的零化多项式,因此我们可以通过多项式取模得到g(x)=xn(modZa(x)),得到的结果g(x)的阶数小于k,因此可以直接求解。

取模的时间复杂度为O(k2log2n),如果使用FFT则可以优化到O(klog2klog2n)

BM算法求解矩阵A的最小多项式

An×n矩阵,且只有m个元素非0。设poly(A)为其最小多项式。

考虑这样的向量数列A0v,A1v,A2v。可以发现poly(A)也是数列的零化多项式(令xi=Aiv)。从这样的数列A0v,A1v,A2v中利用BM求最小递推关系,时间复杂度为O(n3)

对于稀疏矩阵,我们这边可以使用一个混淆的方式,就是求数列uA0v,uA1v,uA2v的最小递推关系,它的时间复杂度是O(n(n+m))

上面提到的u,v都是随机选择的向量。

BM算法求解det(A)

An×n矩阵,且只有m个元素非0。设poly(A)为其最小多项式。

根据特征多项式的定义f(t)=det(tA),可知det(A)=(1)nf(0)。因此我们只需要求出矩阵A的特征多项式,求其常数项即可。

已知的求特征多项式的方式是将矩阵转换成上海森堡形式,之后通过多项式插值求解,时间复杂度为O(n3),和直接通过高斯消元求行列式时间复杂度相同。我们不能直接求解最小多项式,因为最小多项式不一定是特征多项式,不过二者的常数项均只有在矩阵A奇异的情况下才会为0

对于稀疏矩阵,我们可以生成一个随机的对角线矩阵D,之后用BM算法计算DA的最小多项式。而DA的最小多项式有很大的概率就是DA的特征多项式,于是我们得到poly(DA),从而得到det(DA)。而det(A)=det(DA)det(D),这是很好算的。这种方法的时间复杂度为O(n(n+m))

BM算法求解A1b

对于线性方程组Ax=b,如果A可逆,则可以通过x=A1b求解。(这种情况非常常见,比如求随机游走的概率和期望的时候)。

对于一般的密集矩阵,我们可以通过直接对A求逆计算结果,时间复杂度为O(n3)

对于稀疏矩阵。我们可以先计算出A的最小多项式g(A)=ki=0ciAi。之后进行推导:

ki=0ciAi=0c0=ki=1ciAiI=1c0ki=1ciAiA1b=1c0k1i=0ci+1Aib

上面的过程我们可以O(n(n+m))求解。

参考资料