Berlekamp-Massey 알고리즘

Berlekamp-Massey 알고리즘은 특정한 DP의 점화식을 찾아주는 알고리즘이다. $10^{18}$ 번째 피보나치 수를 찾기 위해서 행렬 곱셈을 짜고, 타일 채우기 문제를 풀기 위해서 수많은 점화식과 씨름하던 옛 시간은 이젠 안녕. 이제는 백트래킹 짜고 하드코딩해서 넣으면 끝난다.

이 글은 알고리즘의 구현법, 동작 원리나 증명에 대해서 거의 설명하지 않는다. 그 이유는 내가 구현법과 동작 원리, 증명을 모르기 때문이다. 알고리즘 구현은 여기에서 복붙해서 사용하면 된다. 이론적 배경지식이 상당히 깊지만, 그 활용도가 매우 높기 때문에, 일단 이해하지 말고 작동법부터 제대로 깨우친 후, 나중에 다시 돌아와서 방법을 이해하는 것을 추천한다.

1967년 이 알고리즘을 개발한 수학자 Elwyn Berlekamp가 최근 (2019년 4월 9일) 사망했습니다. 고인의 명복을 빕니다. A final game with Elwyn Berlekamp

Introduction

링크에 주어진 코드에는 크게 네 가지 부분이 있다.

  • berlekamp_massey 함수. $n$개의 DP 결과값이 주어졌을 때, Berlekamp-Massey 알고리즘을 사용하여 해당 DP의 점화식을 찾아주는 루틴. $O(n^2)$ 시간 복잡도에 작동한다 (정확히는, 답이 되는 점화식의 크기가 $k$ 일때 $O(nk+n\log mod)$)
  • get_nth 함수. 주어진 DP 점화식을 사용하여 $n$번째 항을 구하는 루틴.
  • get_min_poly 함수. 3번 챕터에서 설명한다. 지금은 몰라도 상관없다.
  • det 함수. 4번 챕터에서 설명한다. 지금은 몰라도 상관없다.

berlekamp_massey 함수에 구현된 Berlekamp-Massey는 입력 수열을 받아서, 위 수열을 만족하는 가장 짧은 점화식 $DP_n = \sum_{i = 1}^{k}{DP_{n - i} \times C_i}$ 을 찾는다. 이 때 DP 점화식의 형태는 정확히 저러한 형태여야 하며, 이러한 점화식을 선형 점화식 (Linear recurrence) 라고 부른다. 예를 들어, $F_n = F_{n-1} + F_{n-2}$ 라는 점화식을 가지는 피보나치 수열, $P_n = P_{n-1} + P_{n-5}$ 라는 점화식을 가지는 파도반 수열 등이 선형 점화식의 예이다. $C_n = \sum_{i = 0}^{n-1}{C_i \times C_{n-i-1}}$ 형태의 점화식을 가지는 카탈란 수열은 선형 점화식의 예가 아니기 때문에 이 알고리즘을 통해서 구할 수 없다. 어떠한 패턴의 점화식이 계산되는지는 아래에서 계속 설명한다.

get_nth 함수는 Berlekamp-Massey 알고리즘과 전혀 상관이 없는 코드로, 점화식이 주어졌을 때 이 점화식에서 $DP_n$ 이 무엇인지를 계산해 주는 함수이다. 이 함수의 시간 복잡도는 $O(k^2 \log n)$ 이다. 일반적인 DP 점화식을 계산하는 데 드는 시간 복잡도는 $O(nk)$ 이고, 실질적으로 $k \le n$ 이니, 대부분의 경우 충분히 빠른 알고리즘이다.

get_nth 함수는 흔히 키타마사법이라고 불리는 방법으로 구현되어 있다. 흔히 저러한 DP 점화식을 구하는 데는 행렬을 사용하여 $O(k^3 \log n)$ 에 구하는 방법이 잘 알려져 있으나, 키타마사법을 사용하면 더 빠른 시간 복잡도인 $O(k^2 \log n)$ 에 이를 구할 수 있다. 키타마사법에 대해서는 이 동영상 이 블로그 에 잘 설명되어 있으나, 여기서도 한번 더 욕심을 내서 설명해 보겠다. (관심이 없으면 넘겨도 이후 글을 읽는데 상관이 없다.)

키타마사법은 생성함수와 비슷한 형태로 이해를 하면 상당히 간단하다. $DP_n$ 을 구하는 과정을 생각해 보자. $DP_n = \sum_{i = 1}^{k}{DP_{n - i} \times C_i}$ 으로 대체(substitution)해주면, 더 작은 DP값들인 $DP_{n-1}, DP_{n-2}, \cdots, DP_{n-k}$ 에 대한 식으로 나타내어 줄 수 있다. $DP_{n-1}$ 을 또 대체해주면 $DP_{n-2} \cdots DP_{n-k-1}$ 가 되고, $DP_{n-2}, DP_{n-3}$ 등을 점화식을 사용하여 계속 똑같은 변환 방법으로 대체해주면 결국 언젠가는 $C_0 \times DP_0 + C_1 \times DP_1 + \cdots + C_{k-1} \times DP_{k-1}$ 형태가 되어서 문제를 풀 수 있을 것이다. 흥미로운 관찰은, 여기서 최대항을 대체한다는 것은 일반적으로 수학에서 다항식에 나머지 연산을 취하는 것과 동일하다는 것이다. $DP_n$$x^n$ 이라고 보자. $x^n$$\sum_{i=1}^{k}{C_i x^{n-i}}$ 로 대체하는 것을 반복하는 연산은, 사실 $x^n$$(x^k - \sum_{i=0}^{k-1}{C_{i+1} \times x^i})$ 으로 나눈 나머지를 구하는 것과 동일하다.

$x^n$ 을 단순히 $(x^k - \sum_{i=0}^{k-1}{C_{i+1} \times x^i})$ 으로 나눈 나머지를 구하는 것은 물론 효율적이지 않다. 이를 단순히 구현하면 $O(nk)$ 의 시간 복잡도가 들어서 일반적인 DP 계산과 차이가 없다 (아마 코드도 일반적인 DP 계산과 별로 다르지 않을 것이다). 그런데, 일반적인 수에 대해서 $x^n \mod p$$O(\log n)$ 시간에 구하는 다음과 같은 빠른 알고리즘이 잘 알려져 있다:

  • $x^{2k+1} = x^{2k} \times x \mod p$
  • $x^{2k} = x^k \times x^k \mod p$

여기서 $p$ 를 단순히 다항식으로 치환해주면, $O(\log n)$ 번의 나머지 연산과 곱셈만으로 문제를 해결할 수 있다. 다항식에 mod가 되어 있으니 항의 크기는 항상 $k$ 이하임이 보장되고, 이 경우 두 다항식의 곱셈과 나눗셈은 모두 $O(k^2)$ 에 구현할 수 있다. 고로 $O(k^2 \log n)$ 에 문제를 풀 수 있고, 다항식 곱셈과 나눗셈을 $O(k \log k$) 에 구현하면 복잡도가 $O(k\log k\log n)$ 으로 줄어든다. (곱셈은 FFT를 사용하면 매우 간단하고, 나눗셈 역시 FFT를 사용하면 되나 간단하지는 않으니 관심이 있으면 찾아보도록 하자.)

1. Find minimum linear recurrence

모든 것을 알았다면 Berlekamp-Massey 알고리즘을 사용하는 것은 간단하다. 예를 들어, 점화식이 $F_n = F_{n-1} + F_{n-2}$ 인 피보나치 함수의 점화식 {1, 1, 2, 3, 5, 8, 13, 21} 을 넣었다면, 결과물은 {1, 1} 이다. 파도반 수열의 경우에는 $P_n = P_{n-2} + P_{n-3}$ 과 같은 동치의 점화식이 존재하기 때문에, 이 경우 {0, 1, 1} 이 반환된다.

조심해야 할 것은 두 가지다. 첫 번째로, 초기항이 충분히 많아야 한다는 점이다. 만약에 점화식이 $DP_n = \sum_{i = 1}^{k}{DP_{n - i} \times C_i}$ 의 형태라면, 항을 최소한 $3k$ 개보다 많이 넣어줘야 한다. Berlekamp-Massey는 가장 짧은(shortest), 즉 $k$ 를 최소로 하는 선형 점화식을 반환한다. 예를 들어, 만약 항을 $x < 2k$ 개 넣었다면, 항상 $x/2$ 정도 크기의 선형 점화식을 찾을 수 있다 (방정식을 풀어주면 해가 나올 것이니..). 해당 점화식에 대해서 원하는 결과를 얻기 위해서는 최소 $3k$ 개의 항을 제공해야 한다. $2k$ 가 아니라 $3k$ 인지, $3k$ 를 넘어가면 자명하지 않은 선형 점화식이 항상 유일한지는 사실 잘 모르겠다. 꽤 중요한 질문인데 (예를 들어 선형 점화식이 유일하지 않으면, 조건을 만족하는 가장 짧은 점화식이 우리가 원하는 점화식이 아니라 틀릴 수 있으니) 이 부분에 대해서 질문한 사람도 없고 답을 구하지도 못했다. 좋은 의견이 있다면 댓글로 알려주세요.

두 번째로는, 모듈러가 소수여야 한다는 점이다. 만약에 BigInteger 형태로 문제를 해결하고 싶다면, 적당한 소수로 해싱한 후 Berlekamp-Massey를 사용하면 될 것이다 (물론 이 경우 get_nth 함수는 쓸모가 없을 것이다). 모듈러가 소수가 아니라면 더 복잡한 Reeds-Sloane 알고리즘을 사용해야 한다.

사실 여기까지 들어보면 상당히 쓸데 없는 것 같다. 왜냐하면 대부분의 복잡한 문제들은 점화식이 1차원으로 끝나지 않고, 2차원이거나 여러 선형 점화식의 연립으로 표현되기 때문이다. 하지만…

2. Find minimum linear recurrence from a matrix

일반적으로 동적 계획법이 저렇게 깔끔한 선형 점화식 하나로 표현되는 경우는 많지 않다. 예를 들어서, 2차원 점화식의 형태이거나 ($dp_{i, j} = dp_{i-1, j} + dp_{i-1, j-1}$). 혹은 여러 개의 선형 점화식이 서로간의 dependency를 가지는 식 ($F_i = G_{i-1} - H{i-1}, G_i = H_{i-1} \times 2$..) 이 등장하는 경우들이 있다. 이러한 경우는 위의 경우처럼 선형 점화식으로 모델링하는 것이 불가능해 보이지만, 사실은 가능하다.

이러한 식들을 행렬로 나타내 보자. 예를 들어 선형 점화식인 $DP_{i} = DP_{i-1} + DP_{i-2} + DP_{i-5}$

$ \begin{bmatrix} DP_i \\ DP_{i-1} \\ DP_{i-2} \\ DP_{i-3} \\ DP_{i-4} \end{bmatrix} = \begin{bmatrix} 1 & 1 & 0 & 0 & 1 \\ 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \end{bmatrix} \times \begin{bmatrix} DP_{i-1} \\ DP_{i-2} \\ DP_{i-3} \\ DP_{i-4} \\ DP_{i-5} \end{bmatrix}$

와 같은 행렬로 나타낼 수 있다. 그 외에도, $DP_{i, j} = 2\times DP_{i-1, j} + 5 \times DP_{i-1, j-1}$ 과 같은 식은 ($0 \le j \le 4$ 라 했을 때)

$ \begin{bmatrix} DP_{i, 0} \\ DP_{i, 1} \\ DP_{i, 2} \\ DP_{i, 3} \\ DP_{i, 4} \end{bmatrix} = \begin{bmatrix} 2 & 0 & 0 & 0 & 0 \\ 5 & 2 & 0 & 0 & 0 \\ 0 & 5 & 2 & 0 & 0 \\ 0 & 0 & 5 & 2 & 0 \\ 0 & 0 & 0 & 5 & 2 \end{bmatrix} \times \begin{bmatrix} DP_{i-1, 0} \\ DP_{i-1, 1} \\ DP_{i-1, 2} \\ DP_{i-1, 3} \\ DP_{i-1, 4} \end{bmatrix}$

와 같은 행렬로 나타낼 수 있고. $F_i = 2 \times G_{i-1}, G_i = F_{i-1} + 3 \times H_{i-1}, H_i = F_{i-1} + G_{i-1}$ 와 같은 식은,

$\begin{bmatrix} F_i \\ G_i \\ H_i \end{bmatrix} = \begin{bmatrix} 0 & 2 & 0 \\ 1 & 0 & 3 \\ 1 & 1 & 0 \end{bmatrix} \times \begin{bmatrix} F_{i-1} \\ G_{i-1} \\ H_{i-1} \end{bmatrix}$

로 나타낼 수 있다.

중요한 예시들을 대부분 적었으나 이 예시가 전부는 아니다. 행렬식으로 바꿀 수 없는 점화식이 뭔지, 직관적으로 행렬식으로 바꿀 수 없어 보이지만 사실 가능한 것들은 뭔지에 대해서 생각해 보는 것은 좋을 것이다. 예를 들어서, $F_i = F_{i-1} + G_{i-1} + 10, G_i = 2 \times F_{i} - G_{i-1} - 6$ 과 같은 식은 행렬 식으로 바꿀 수 없어 보이지만.

$\begin{bmatrix} F_i \\ G_i \\ 1 \end{bmatrix} = \begin{bmatrix} 1 & 1 & 10 \\ 2 & -1 & -6 \\ 0 & 0 & 1 \end{bmatrix} \times \begin{bmatrix} F_{i-1} \\ G_{i-1} \\ 1 \end{bmatrix}$

과 같은 형태로 행렬식으로 쓸 수 있으며, 한편 $C_i = \sum_{0 \le j \le n-1}{C_j \times C_{i-j-1}}$ 과 같은 카탈란 점화식은 이러한 트릭을 쓸 수 없다.

이제 중요한 결론을 제시한다: DP 상태 전이가 행렬로 표현 가능하다면, 해당 DP에 대한 선형 점화식이 존재한다. 정확히는, 위 DP 식에서의 $DP_i, DP_{i ,0}, F_i$ (각각) 에 대해서 선형 점화식이 존재하며, 이 점화식의 크기는 행렬의 크기가 $n \times n$ 일때 최대 $n$이다. 고로 3n개의 DP 값을 직접 계산한 후, Berlekamp-Massey를 사용하면 $O(n^2)$ 시간에 점화식을 찾을 수 있으며, 이 점화식의 $k$ 번째 항은 $O(n^2 \log k)$ 에 계산할 수 있다.

이러한 방식의 장점은 다음과 같다:

  • 점화식이 저렇게 생겼을 거 같으면 찍어보면 된다. 어떠한 경우에는, $DP_{i, state}$ 에 대한 점화식을 찾았는데 정확한 상태 전이를 찾기 귀찮거나, 아예 그것을 찾는 것이 문제일 수 있다. 어떠한 경우에는, 점화식을 찾지 못했으나 점화식의 형태가 위와 같을 거라는 직관만이 있을 수도 있다. 어떠한 경우에는 아는 게 아무것도 없고, 그냥 이것이 유일한 희망일 수도 있다. 이러한 경우에, 단순히 Berlekamp-Massey를 로컬에서 돌려본 후 작은 점화식이 있는지를 확인해 볼 수 있다. 만약에 충분히 큰 $n$에 대해서 $\frac{n}{2}$ 보다 유의미하게 작은 점화식을 꾸준히 찾을 수 있다면, 문제를 거의 다 풀었다는 좋은 청신호이다.
  • 시간 복잡도가 빠르다. 선형 점화식 형태의 전이는 위의 키타마사법이라는 방법으로 빠르게 계산하는 방법이 잘 알려져 있으나, 행렬 형태의 전이는 보통 행렬 곱셈을 동반하여 이러한 최적화가 불가능하다. Berlekamp-Massey는 후자의 전이를 전자로 변환해 주는 알고리즘이 된다.
  • 커팅을 동반한다. Berlekamp-Massey는 최소 크기의 점화식을 찾는다. 고로 꼭 행렬의 크기가 $n \times n$ 이라고 점화식의 크기가 $n$ 이 나오는 것이 아니고 더 작은 점화식이 나올 수 있다. 예를 들어서 $DP_{i, state}$ 와 같은 점화식을 찾았는데 상태 전이가 작고 (sparse matrix) 특정 상태가 도달 불가능한 상황이 있을 수 있다. Berlekamp-Massey를 돌리면 그러한 경우를 우리가 따져 줄 필요가 없이, 알아서 최적의 점화식을 찾아준다. 이는 매우 중요한 요소이다. 어떠한 문제의 경우에는 이러한 커팅을 찾는 것이 (요컨데, 방문할 수 없는 상태를 찾아 시간 복잡도나 상수를 유의미하게 줄이는 것이) 문제의 어려운 부분일 수 있다. Berlekamp-Massey 를 사용하면, 이러한 문제를 내가 해결할 필요가 없이 알고리즘단에서 이러한 커팅을 최적으로 사용할 수 있다.

Berlekamp-Massey 알고리즘을 사용한 다양한 연습 문제와 풀이, 몇가지 사실의 증명, 그리고 get_min_poly와 det 함수에 대한 이야기는 구사과 블로그 에서 계속됩니다.

댓글 (1개) 댓글 쓰기