Berlekamp–Massey Algorithm

The Berlekamp–Massey algorithm is used to find the shortest recurrence relation of a sequence. Given a sequence of length \(n\), if the order of its shortest recurrence is \(m\), then the Berlekamp–Massey algorithm can compute the shortest recurrence for each prefix of the sequence in \(O(nm)\) time. In the worst case, \(m = O(n)\), so the worst-case complexity of the algorithm is \(O(n^2)\).

Definition

Define a recurrence relation of a sequence \(\{a_0 \dots a_{n - 1} \}\) as a sequence \(\{r_0\dots r_m\}\) satisfying:

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

where \(r_0 = 1\). \(m\) is called the order of the recurrence relation.

The shortest recurrence relation of the sequence \(\{a_i\}\) is the one with the smallest order.

Algorithm

Slightly different from the definition above, here we define a new recurrence coefficient sequence \(\{f_0 \dots f_{m - 1}\}\) satisfying:

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

It is easy to see that \(f_i = -r_{i + 1}\), and the order \(m\) is the same as in the previous definition.

We can compute the recurrence relation incrementally, considering each term of \(\{a_i\}\) in order, and adjusting the recurrence coefficients \(\{f_i\}\) when the recurrence result is incorrect. For convenience, let \(F_i = \{f_{i, j}\}\) denote the shortest recurrence relation for the first \(i\) terms.

Clearly, initially \(F_0 = \{\}\). Suppose the recurrence coefficients \(F_{i - 1}\) hold for the first \(i - 1\) terms of the sequence \(\{a_i\}\). Then for the \(i\)-th term, there are two cases:

The recurrence coefficients also hold for \(a_i\). In this case, no adjustment is needed; simply set \(F_i = F_{i - 1}\).
The recurrence coefficients do not hold for \(a_i\). Then \(F_{i - 1}\) must be adjusted to obtain the new \(F_i\).

Let \(\Delta_i = a_i - \sum_{j = 0} ^ m f_{i - 1, j} a_{i - j - 1}\), i.e., the difference between \(a_i\) and the recurrence result of \(F_{i - 1}\).

If this is the first time adjusting the recurrence coefficients, it means \(a_i\) is the first non-zero term in the sequence. In this case, simply set \(F_i\) to \(i\) zeros, which is clearly a valid shortest recurrence relation.

Otherwise, let \(k\) be the number of considered terms of \(\{a_i\}\) when the recurrence coefficients were last adjusted. If there exists a sequence \(G = \{g_0 \dots g_{m' - 1}\}\) satisfying:

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

and \(\sum_{j = 0} ^ {m' - 1} g_j a_{i - j - 1} = \Delta_i\), then it is not hard to see that adding \(F_k\) and \(G\) element-wise yields a valid recurrence coefficient sequence \(F_i\).

Consider how to construct \(G\). One feasible construction is:

\(G = \{0, 0, \dots, 0, \frac{\Delta_i}{\Delta_k}, -\frac{\Delta_i}{\Delta_k}F_{k-1}\}\)

where there are exactly \(i - k - 1\) zeros in front, and \(-\frac{\Delta_i}{\Delta_k} F_{k-1}\) means appending \(F_{k-1}\) with each term multiplied by \(-\frac{\Delta_i}{\Delta_k}\).

It is easy to verify that \(\sum_{j = 0} ^ {m' - 1} g_j a_{i - j - 1} = \Delta_k \frac{\Delta_i}{\Delta_k} = \Delta_i\), so this construction yields a valid \(G\). Then set \(F_i\) to the element-wise sum of \(F_k\) and \(G\).

If the recurrence relation \(\{r_i\}\) as defined initially is required, simply negate all \(\{f_j\}\) and prepend \(r_0 = 1\).

From the above algorithm, it can be seen that if the order of the shortest recurrence of the sequence is \(m\), the complexity of the algorithm is \(O(nm)\). In the worst case, \(m = O(n)\), so the worst-case complexity is \(O(n^2)\).

When implementing the algorithm, since only the recurrence coefficients \(F_k\) from the last adjustment are needed each time, if only the shortest recurrence of the entire sequence is required, we can store only the current recurrence coefficients and those from the last adjustment, resulting in \(O(n)\) space complexity.

Reference Implementation

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$
}

The naive Berlekamp–Massey algorithm finds the shortest recurrence for a finite sequence. If the sequence for which the recurrence is sought is infinite, but an upper bound \(m\) on the order of the shortest recurrence is known, then taking the first \(2m\) terms of the sequence is sufficient to compute the shortest recurrence for the entire sequence. (Proof omitted.)

Applications

Due to poor numerical stability, the Berlekamp–Massey algorithm is rarely used for real-number problems. For convenience, the following assumes operations are performed modulo a prime \(p\).

Finding the Shortest Recurrence for a Vector or Matrix Sequence

To find the shortest recurrence for a vector sequence \(\boldsymbol{v}_i\), where the vectors have dimension \(n\), we can randomly choose an \(n\)-dimensional row vector \(\mathbf u^T\) and compute the shortest recurrence for the scalar sequence \(\{\boldsymbol{u}^T\boldsymbol{v}_i\}\). By the Schwartz–Zippel lemma, the shortest recurrence of the two sequences is the same with probability at least \(1 - \frac n p\).

Finding the shortest recurrence for a matrix sequence \(\{A_i\}\) is similar. If the matrices are of size \(n \times m\), randomly choose a \(1 \times n\) row vector \(\mathbf u^T\) and an \(m \times 1\) column vector \(\boldsymbol{v}\), then compute the shortest recurrence for the scalar sequence \(\{\boldsymbol{u}^T A_i \boldsymbol{v}\}\). By the Schwartz–Zippel lemma, the probability that the two recurrences are the same is at least \(1 - \frac{n + m} p\).

Optimizing Matrix Exponentiation

Let \(\boldsymbol{f}_i\) be an \(n\)-dimensional column vector with transition \(\boldsymbol{f}_i = A \boldsymbol{f}_{i - 1}\). Then \(\{\boldsymbol{f}_i\}\) is a linear recurrence vector sequence of order at most \(n\). (Proof omitted.)

We can directly compute \(\boldsymbol{f}_0 \dots \boldsymbol{f}_{2n - 1}\) by brute force, then use the method mentioned earlier to find the shortest recurrence of \(\{\boldsymbol{f}_i\}\), and finally invoke Linear Recurrence with Constant Coefficients.

If the target vector is \(\boldsymbol{f}_m\), the complexity is \(O(n^3 + n\log n \log m)\). If \(A\) is a sparse matrix with \(k\) non-zero entries, the complexity can be reduced to \(O(nk + n\log n \log m)\). However, since the algorithm requires at least \(O(nk)\) time for preprocessing, the \(O(n^2 \log m)\) linear recurrence algorithm is also acceptable when the constraints are not too tight.

Finding the Minimal Polynomial of a Matrix

The minimal polynomial of a square matrix \(A\) is the polynomial \(f\) of minimal degree such that \(f(A) = 0\).

In fact, the minimal polynomial is exactly the shortest recurrence of \(\{A^i\}\), so we can directly call the Berlekamp–Massey algorithm. If \(A\) is an \(n \times n\) matrix, the degree of the minimal polynomial is clearly at most \(n\).

The bottleneck is computing \(A^i\), because directly performing matrix multiplication each time would yield \(O(n^4)\) complexity. However, note that when finding the shortest recurrence of a matrix sequence, we actually find the shortest recurrence of \(\{\boldsymbol{u}^T A^i \boldsymbol{v}\}\), so we only need to compute \(A^i \boldsymbol{v}\).

Assuming \(A\) has \(k\) non-zero entries, the complexity is \(O(kn + n^2)\).

Computing the Determinant of a Sparse Matrix

If the characteristic polynomial of a square matrix \(A\) can be found, then the determinant is the constant term multiplied by \((-1)^n\). However, the minimal polynomial is not necessarily the characteristic polynomial.

In fact, if we multiply \(A\) by a random diagonal matrix \(B\), then the minimal polynomial of \(AB\) is the characteristic polynomial with probability at least \(1 - \frac {2n^2 - n} p\). Finally, divide by \(\text{det}\;B\).

Let \(A\) be an \(n \times n\) matrix with \(k\) non-zero entries. The complexity is \(O(kn + n ^ 2)\).

Computing the Rank of a Sparse Matrix

Let \(A\) be an \(n\times m\) matrix. First, randomly choose an \(n\times n\) diagonal matrix \(P\) and an \(m\times m\) diagonal matrix \(Q\), then compute the minimal polynomial of \(Q A P A^T Q\).

Actually, matrix multiplication is not needed, because when computing the minimal polynomial, we only need to multiply \(Q A P A^T Q\) by a vector, so we can sequentially multiply these matrices with the vector. The answer is the degree of the minimal polynomial after removing all factors of \(x\).

Let \(A\) have \(k\) non-zero entries, and assume \(n \le m\). The complexity is \(O(kn + n ^ 2)\).

Solving Sparse Linear Systems

Problem: Given \(A \mathbf x = \mathbf b\), where \(A\) is an \(n \times n\) full-rank sparse matrix, \(\mathbf b\) and \(\mathbf x\) are \(1\times n\) column vectors. \(A\) and \(\mathbf b\) are known, and we need to solve for \(x\) in complexity lower than \(n^\omega\).

Solution: Clearly \(\mathbf x = A^{-1} \mathbf b\). If we can find the shortest recurrence \(\{r_0 \dots r_{m - 1}\}\) (\(m \le n\)) of \(\{A^i \mathbf b\}\) (\(i \ge 0\)), then we have:

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

(Proof omitted.)

Since \(A\) is sparse, we can directly compute \(\mathbf b \dots A^{2n - 1} \mathbf b\) by definition.

Again, let \(A\) have \(k\) non-zero entries. The complexity is \(O(kn + n^2)\).

Reference Implementation

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;
}

Example Problems

本页面最近更新：，更新历史
发现错误？想一起完善？在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用