跳转至

Berlekamp–Massey 算法

Berlekamp–Massey 算法是一種用於求數列的最短遞推式的算法。給定一個長為 \(n\) 的數列，如果它的最短遞推式的階數為 \(m\)，則 Berlekamp–Massey 算法能夠在 \(O(nm)\) 時間內求出數列的每個前綴的最短遞推式。最壞情況下 \(m = O(n)\)，因此算法的最壞複雜度為 \(O(n^2)\)。

定義

定義一個數列 \(\{a_0 \dots a_{n - 1} \}\) 的遞推式為滿足下式的序列 \(\{r_0\dots r_m\}\)：

\(\sum_{j = 0} ^ m r_j a_{i - j} = 0, \forall i \ge m\)

其中 \(r_0 = 1\)。\(m\) 稱為該遞推式的 階數。

數列 \(\{a_i\}\) 的最短遞推式即為階數最小的遞推式。

做法

與上面定義的稍有不同，這裏定義一個新的遞推係數 \(\{f_0 \dots f_{m - 1}\}\)，滿足：

\(a_i = \sum_{j = 0} ^ {m - 1} f_j a_{i - j - 1}, \forall i \ge m\)

容易看出 \(f_i = -r_{i + 1}\)，並且階數 \(m\) 與之前的定義是相同的。

我們可以增量地求遞推式，按順序考慮 \(\{a_i\}\) 的每一位，並在遞推結果出現錯誤時對遞推係數 \(\{f_i\}\) 進行調整。方便起見，以下將前 \(i\) 位的最短遞推式記為 \(F_i = \{f_{i, j}\}\)。

顯然初始時有 \(F_0 = \{\}\)。假設遞推係數 \(F_{i - 1}\) 對數列 \(\{a_i\}\) 的前 \(i - 1\) 項均成立，這時對第 \(i\) 項就有兩種情況：

  1. 遞推係數對 \(a_i\) 也成立，這時不需要進行任何調整，直接令 \(F_i = F_{i - 1}\) 即可。
  2. 遞推係數對 \(a_i\) 不成立，這時需要對 \(F_{i - 1}\) 進行調整，得到新的 \(F_i\)。

設 \(\Delta_i = a_i - \sum_{j = 0} ^ m f_{i - 1, j} a_{i - j - 1}\)，即 \(a_i\) 與 \(F_{i - 1}\) 的遞推結果的差值。

如果這是第一次對遞推係數進行修改，則説明 \(a_i\) 是序列中的第一個非零項。這時直接令 \(F_i\) 為 \(i\) 個 \(0\) 即可，顯然這是一個合法的最短遞推式。

否則設上一次對遞推係數進行修改時，已考慮的 \(\{a_i\}\) 的項數為 \(k\)。如果存在一個序列 \(G = \{g_0 \dots g_{m' - 1}\}\)，滿足：

\(\sum_{j = 0} ^ {m' - 1} g_j a_{i' - j - 1} = 0, \forall i' \in [m', i)\)

並且 \(\sum_{j = 0} ^ {m' - 1} g_j a_{i - j - 1} = \Delta_i\)，那麼不難發現將 \(F_k\) 與 \(G\) 按位分別相加之後即可得到一個合法的遞推係數 \(F_i\)。

考慮如何構造 \(G\)。一種可行的構造方案是令

\(G = \{0, 0, \dots, 0, \frac{\Delta_i}{\Delta_k}, -\frac{\Delta_i}{\Delta_k}F_k\}\)

其中前面一共有 \(i - k - 1\) 個 \(0\)，且最後的 \(-\frac{\Delta_i}{\Delta_k} F_k\) 表示將 \(F_k\) 每項乘以 \(-\frac{\Delta_i}{\Delta_k}\) 後接在序列後面。

不難驗證此時 \(\sum_{j = 0} ^ {m' - 1} g_j a_{i - j - 1} = \Delta_k \frac{\Delta_i}{\Delta_k} = \Delta_i\)，因此這樣構造出的是一個合法的 \(G\)。將 \(F_i\) 賦值為 \(F_k\) 與 \(G\) 逐項相加後的結果即可。

如果要求的是符合最開始定義的遞推式 \(\{r_i\}\)，則將 \(\{f_j\}\) 全部取相反數後在最開始插入 \(r_0 = 1\) 即可。

從上述算法流程中可以看出，如果數列的最短遞推式的階數為 \(m\)，則算法的複雜度為 \(O(nm)\)。最壞情況下 \(m = O(n)\)，因此算法的最壞複雜度為 \(O(n^2)\)。

在實現算法時，由於每次調整遞推係數時都只需要用到上次調整時的遞推係數 \(F_k\)，因此如果只需要求整個數列的最短遞推式，可以只存儲當前遞推係數和上次調整時的遞推係數，空間複雜度為 \(O(n)\)。

參考實現 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
vector<int> berlekamp_massey(const vector<int> &a) {
  vector<int> v, last;  // v is the answer, 0-based, p is the module
  int k = -1, delta = 0;

  for (int i = 0; i < (int)a.size(); i++) {
    int tmp = 0;
    for (int j = 0; j < (int)v.size(); j++)
      tmp = (tmp + (long long)a[i - j - 1] * v[j]) % p;

    if (a[i] == tmp) continue;

    if (k < 0) {
      k = i;
      delta = (a[i] - tmp + p) % p;
      v = vector<int>(i + 1);

      continue;
    }

    vector<int> u = v;
    int val = (long long)(a[i] - tmp + p) * power(delta, p - 2) % p;

    if (v.size() < last.size() + i - k) v.resize(last.size() + i - k);

    (v[i - k - 1] += val) %= p;

    for (int j = 0; j < (int)last.size(); j++) {
      v[i - k + j] = (v[i - k + j] - (long long)val * last[j]) % p;
      if (v[i - k + j] < 0) v[i - k + j] += p;
    }

    if ((int)u.size() - i < (int)last.size() - k) {
      last = u;
      k = i;
      delta = a[i] - tmp;
      if (delta < 0) delta += p;
    }
  }

  for (auto &x : v) x = (p - x) % p;
  v.insert(v.begin(), 1);

  return v;  // $\forall i, \sum_{j = 0} ^ m a_{i - j} v_j = 0$
}

樸素的 Berlekamp–Massey 算法求解的是有限項數列的最短遞推式。如果待求遞推式的序列有無限項，但已知最短遞推式的階數上界，則只需取出序列的前 \(2m\) 項即可求出整個序列的最短遞推式。（證明略）

應用

由於 Berlekamp–Massey 算法的數值穩定性比較差，在處理實數問題時一般很少使用。為了敍述方便，以下均假定在某個質數 \(p\) 的剩餘系下進行運算。

求向量列或矩陣列的最短遞推式

如果要求向量列 \(\boldsymbol{v}_i\) 的最短遞推式，設向量的維數為 \(n\)，我們可以隨機一個 \(n\) 維行向量 \(\mathbf u^T\)，並計算標量序列 \(\{\boldsymbol{u}^T\boldsymbol{v}_i\}\) 的最短遞推式。由 Schwartz–Zippel 引理，二者的最短遞推式有至少 \(1 - \frac n p\) 的概率相同。

求矩陣列 \(\{A_i\}\) 的最短遞推式也是類似的，設矩陣的大小為 \(n \times m\)，則只需隨機一個 \(1 \times n\) 的行向量 \(\mathbf u^T\) 和一個 \(m \times 1\) 的列向量 \(\boldsymbol{v}\)，並計算標量序列 \(\{\boldsymbol{u}^T A_i \boldsymbol{v}\}\) 的最短遞推式即可。由 Schwartz–Zippel 引理可以類似地得到二者相同的概率至少為 \(1 - \frac{n + m} p\)。

優化矩陣快速冪

設 \(\boldsymbol{f}_i\) 是一個 \(n\) 維列向量，並且轉移滿足 \(\boldsymbol{f}_i = A \boldsymbol{f}_{i - 1}\)，則可以發現 \(\{\boldsymbol{f}_i\}\) 是一個不超過 \(n\) 階的線性遞推向量列。（證明略）

我們可以直接暴力求出 \(\boldsymbol{f}_0 \dots \boldsymbol{f}_{2n - 1}\)，然後用前面提到的做法求出 \(\{\boldsymbol{f}_i\}\) 的最短遞推式，再調用 常係數齊次線性遞推 即可。

如果要求的向量是 \(\boldsymbol{f}_m\)，則算法的複雜度是 \(O(n^3 + n\log n \log m)\)。如果 \(A\) 是一個只有 \(k\) 個非零項的稀疏矩陣，則複雜度可以降為 \(O(nk + n\log n \log m)\)。但由於算法至少需要 \(O(nk)\) 的時間預處理，因此在壓力不大的情況下也可以使用 \(O(n^2 \log m)\) 的線性遞推算法，複雜度同樣是可以接受的。

求矩陣的最小多項式

方陣 \(A\) 的最小多項式是次數最小的並且滿足 \(f(A) = 0\) 的多項式 \(f\)。

實際上最小多項式就是 \(\{A^i\}\) 的最小遞推式，所以直接調用 Berlekamp–Massey 算法就可以了。如果 \(A\) 是一個 \(n\) 階方陣，則顯然最小多項式的次數不超過 \(n\)。

瓶頸在於求出 \(A^i\)，因為如果直接每次做矩陣乘法的話複雜度會達到 \(O(n^4)\)。但考慮到求矩陣列的最短遞推式時實際上求的是 \(\{\boldsymbol{u}^T A^i \boldsymbol{v}\}\) 的最短遞推式，因此我們只要求出 \(A^i \boldsymbol{v}\) 就行了。

假設 \(A\) 有 \(k\) 個非零項，則複雜度為 \(O(kn + n^2)\)。

求稀疏矩陣行列式

如果能求出方陣 \(A\) 的特徵多項式，則常數項乘上 \((-1)^n\) 就是行列式。但是最小多項式不一定就是特徵多項式。

實際上如果把 \(A\) 乘上一個隨機對角陣 \(B\)，則 \(AB\) 的最小多項式有至少 \(1 - \frac {2n^2 - n} p\) 的概率就是特徵多項式。最後再除掉 \(\text{det}\;B\) 就行了。

設 \(A\) 為 \(n\) 階方陣，且有 \(k\) 個非零項，則複雜度為 \(O(kn + n ^ 2)\)。

求稀疏矩陣的秩

設 \(A\) 是一個 \(n\times m\) 的矩陣，首先隨機一個 \(n\times n\) 的對角陣 \(P\) 和一個 \(m\times m\) 的對角陣 \(Q\), 然後計算 \(Q A P A^T Q\) 的最小多項式即可。

實際上不用調用矩陣乘法，因為求最小多項式時要用 \(Q A P A^T Q\) 乘一個向量，所以我們依次把這幾個矩陣乘到向量裏就行了。答案就是最小多項式除掉所有 \(x\) 因子後剩下的次數。

設 \(A\) 有 \(k\) 個非零項，且 \(n \le m\)，則複雜度為 \(O(kn + n ^ 2)\)。

解稀疏方程組

問題：已知 \(A \mathbf x = \mathbf b\), 其中 \(A\) 是一個 \(n \times n\) 的 滿秩 稀疏矩陣，\(\mathbf b\) 和 \(\mathbf x\) 是 \(1\times n\) 的列向量。\(A, \mathbf b\) 已知，需要在低於 \(n^\omega\) 的複雜度內解出 \(x\)。

做法：顯然 \(\mathbf x = A^{-1} \mathbf b\)。如果我們能求出 \(\{A^i \mathbf b\}\)(\(i \ge 0\)) 的最小遞推式 \(\{r_0 \dots r_{m - 1}\}\)(\(m \le n\)), 那麼就有結論

\(A^{-1} \mathbf b = -\frac 1 {r_{m - 1}} \sum_{i = 0} ^ {m - 2} A^i \mathbf b r_{m - 2 - i}\)

（證明略）

因為 \(A\) 是稀疏矩陣，直接按定義遞推出 \(\mathbf b \dots A^{2n - 1} \mathbf b\) 即可。

同樣地，設 \(A\) 中有 \(k\) 個非零項，則複雜度為 \(O(kn + n^2)\)。

參考實現 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
vector<int> solve_sparse_equations(const vector<tuple<int, int, int> > &A,
                                   const vector<int> &b) {
  int n = (int)b.size();  // 0-based

  vector<vector<int> > f({b});

  for (int i = 1; i < 2 * n; i++) {
    vector<int> v(n);
    auto &u = f.back();

    for (auto [x, y, z] : A)  // [x, y, value]
      v[x] = (v[x] + (long long)u[y] * z) % p;

    f.push_back(v);
  }

  vector<int> w(n);
  mt19937 gen;
  for (auto &x : w) x = uniform_int_distribution<int>(1, p - 1)(gen);

  vector<int> a(2 * n);
  for (int i = 0; i < 2 * n; i++)
    for (int j = 0; j < n; j++) a[i] = (a[i] + (long long)f[i][j] * w[j]) % p;

  auto c = berlekamp_massey(a);
  int m = (int)c.size();

  vector<int> ans(n);

  for (int i = 0; i < m - 1; i++)
    for (int j = 0; j < n; j++)
      ans[j] = (ans[j] + (long long)c[m - 2 - i] * f[i][j]) % p;

  int inv = power(p - c[m - 1], p - 2);

  for (int i = 0; i < n; i++) ans[i] = (long long)ans[i] * inv % p;

  return ans;
}

例題

  1. LibreOJ #163. 高斯消元 2
  2. ICPC2021 台北 Gym103443E. Composition with Large Red Plane, Yellow, Black, Gray, and Blue
本页面最近更新：，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：AntiLeaf
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用