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-1} - 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가 어떻게 문제 해결에 도움을 주는지에 대한 간단한 코멘트를 동반한다.

BOJ 11444. 피보나치 수 6 (Easy)

long long n; 
cin >> n; 
cout << BlackBoxLinearAlgebra::guess_nth_term({0, 1, 1, 2, 3, 5, 8, 13}, n) << endl;

BOJ 9461. 파도반 수열 (Easy)

P(1)부터 P(10)까지 첫 10개 숫자는 1, 1, 1, 2, 2, 3, 4, 5, 7, 9이다.

친절해도 너무 친절하다.

BOJ 13430. 합 구하기 (Easy)

문제의 점화식을 그대로 코드로 옮기고 Berlekamp-Massey에 넣으면 된다. 사실 그냥 풀어도 어렵지 않은 문제다. $S(k, n) = S(k-1, n) + S(k, n-1)$.. (그리고 이것이 Berlekamp-Massey가 작동하는 이유의 증명이 된다.)

BOJ 12916. K-Path (Easy)

행렬 거듭제곱으로 풀 수 있다는 사실이 잘 알려진 문제이다. 행렬 거듭제곱으로 풀 수 있으면 Berlekamp-Massey로도 $O(N^3 + N^2 \log K)$ 에 풀 수 있다.

BOJ 13380. Found (Easy)

두 사람의 위치 쌍을 $(s, e)$라고 하면, 이 쌍에 대해서 그래프를 만들면 위 문제랑 똑같다. 고로, BOJ 12916과 똑같이 $O(N^6)$ 나이브 DP를 돌린 후 Berlekamp-Massey를 사용하면 된다.

BOJ 14453. The Way (Medium)

$3 \times N$ 격자에서 무엇을 하니까 DP 상태를 현재 만든 Path의 모양 비슷하게 하면 될 것 같다. Path의 모양이 몇 개 나오지 않을 것 같기 때문이다. 여기서 더 이상은 잘 모르겠다면, 백트래킹을 돌린 후 Berlekamp-Massey를 쓰면 $O(\log N)$ 시간에 된다.

BOJ 1492. 합 (Medium)

$k = 1$ 이면 $N(N+1)/2$ 가 답이라는 사실이 매우 잘 알려져 있다. 그 정도로 잘 알려져 있지는 않지만, 고등학교 수학 시간에 $k = 2$ 일 때 답이 $N(N+1)(2N+1)/6$, $k = 3$ 일때 $N^2(N+1)^2/4$ … 더 알아보기는 싫으나, 아무튼 이러한 결론을 종합하면 $N$에 대한 $k+1$ 차 다항식이라는 것을 짐작할 수 있다. 그러니, 단순 반복문으로 $f(1), f(2), \cdots, f(3k+3)$ 까지 계산한 후 Berlekamp-Massey를 돌리면 너무 쉽게 AC를 받을 수 있다… 어?

위키백과에 따르면, 모든 $k$ 차 다항식 $f(x)$ 에 대해서, $\sum_{i = 0}^{k}{(f(x - i) \times \binom{k}{i} \times (-1)^i)} = 0$ 와 같은 식이 성립한다. $f(x)$ 를 제외하고 전부 우변으로 이항 시키면, 정확히 크기 $k$ 의 선형 점화식이 된다. 고로 Berlekamp-Massey를 사용해서 찾아주면 된다. 사실 이런 식의 풀이는 Lagrange Interpolation의 하위 호환이지만, 이것 역시 된다는 데 그 의미가 있고, 아무튼 Lagrange Interpolation이 뭔지 몰라도 간단히 쓸 수 있으니 좋다.

BOJ 10562. 나이트 (Medium)

비트마스크 DP를 사용해서 이 문제를 $O(8^M \times N)$ 이나 $O(4^M \times N)$ 에 푸는 건 별로 어렵지 않다. 어떤 방식을 사용하든, $2^{2M} \times 2^{2M}$ 정도의 행렬 곱셈을 사용하면 복잡도를 줄일 수 있을 것으로 보인다. 하지만 비트마스크 DP를 짜는 것이 훨씬 간단하니 그냥 그렇게 하고 Berlekamp-Massey를 돌려 보자. 돌려 보면 실제로 $M = 4$ 일 때 크기 21의 점화식으로 충분함을 알 수 있다. 대칭적인 경우가 많고 여러 상태가 도달 불가능하기 때문이다. 이렇게 점화식이 실제 상황에서 매우 짧기 때문에, Berlekamp-Massey를 사용하면 시간제한이 60초인 문제를 0ms에 해결할 수 있다.

CF 506 E. Mr, Kitayuta's Gift (Hard)

이 문제를 $O(N|S|^2)$ 에 푸는 것은 그렇게 어렵지 않다. $DP_{i, j, k} = $ ($i$ 개의 글자를 끼워 넣고 있으며, 현재 $[j, k]$ 구간을 팰린드롬으로 만들려고 하고 있음)과 같이 정의하면 점화식이 상대적으로 간단하다. 느려서 그렇지…

한편, $DP_{i, *, *}$ 가 $DP_{i-1, *, *}$에만 의존한다는 성질을 활용하면, 이를 Berlekamp-Massey를 사용해서 최적화해 볼 수 있다. 행렬의 최대 사이즈가 $|S|^2 \times |S|^2$ 일 것이기 때문이다. 고로 앞 $i \le O(|S|^2)$ 에 대해서 나이브한 DP를 돌리고 Berlekamp-Massey를 사용하면 $O(|S|^4 \log N)$ 시간 복잡도라 문제를 해결하기 너무 느리다… 정말 그럴까??

로컬에서 몇 가지 랜덤 데이터를 넣어보면, 점화식의 길이가 $O(|S|)$ 를 넘어가는 경우를 찾을 수 없음을 알 수 있다. 뭔가 특별한 성질을 기대하고, $i \le O(|S|)$ 까지만 구해본 후 Berlekamp-Massey를 사용하자. AC. 이 풀이가 되는 이유는 정해를 보면 알 수 있다. 정해는 여러 가지 관찰을 사용하여 문제를 특정한 오토마타의 길이 $N$ 경로를 세는 문제로 환원하였다. 이 오토마타의 크기는 우리가 만든 것과 다르게 $O(|S|)$ 밖에 되지 않는다. 우리는 정해가 한 여러가지 관찰에 대해서 아는 바가 없지만, 아무튼 문제가 원래 그런 문제였기 때문에, Berlekamp-Massey를 사용할 경우 짧은 점화식을 찾을 수 있고, 고로 정해보다 빠른 $O(|S|^3 + |S|^2 \log N)$ 풀이로 문제를 해결할 수 있다. 위 Berlekamp-Massey의 3번 장점에 해당된다.

BOJ 13727. 5차원 구사과 초콜릿 (Hard)

위 예시와 비슷한 문제다. 이 문제는 설명을 생략한다. 직접 해보자 :D

증명

이 부분에서는 위 결론이 성립하는 이유를 증명한다. 아래 내용과 직결되는 부분이니 글을 계속 읽어볼 것이면 무조건 살펴보자. 먼저, 애매한 표현을 안 쓰고 최대한 수학적인 표현으로 정리하겠다. 결국 DP 값 하나는 벡터의 원소 하나로 정의되니, 벡터가 주어지면 내적을 통해서 DP 값을 알아낼 수 있다. 고로, 이를 벡터에 대한 이야기로 환원해서

• $p_0 = x, p_i = A \times p_{i-1}$ 로 정의되는 점화식과, 임의의 벡터 $v \in Z_p^n$에 대해서, $X_i = (v \times p_i)$ 일 때, $X_i$ 를 길이 $n$ 이하의 선형 점화식으로 표현할 수 있다.

라는 명제를 증명하도록 한다.

$v \times p_i = v \times A \times p_{i-1}$ 을 만족한다. 이렇게 쓰니 중간에 있는 $A$가 상당히 거슬린다. 최대한 지워보도록 하자. 정방 행렬 $A$에 대해서 $p(\lambda) = det(\lambda I_n - A)$ 를 특성 다항식이라 한다. 또한, Cayley-Hamilton 정리에 대해서, $p(A) = 0$ 이 성립한다. Determinant의 정의에 따라 특성 다항식은 $n$차 다항식이니, $\sum_{i=0}^{n}{c_i A^i} = 0$ 를 만족하는 자명하지 않은 다항식이 존재한다.이항 하면, $A^n = \sum_{i = 0}^{n-1}{-c_i A^i}$ 가 성립한다.

이 결과를 식에 대입하자.

$X_i = v \times p_i$

$ = v \times A^n \times p_{i-n}$

$ = v \times (\sum_{j=0}^{n-1}{-c_jA^j}) \times p_{i-n}$

$ = v \times (\sum_{j=0}^{n-1}{-c_jA^j \times p_{i - n}})$

$ = v \times (\sum_{j=0}^{n-1}{-c_j p_{i - (n - j)}})$

$= \sum_{j = 0}^{n-1}(-c_j (v \times p_{n - j})$

$= \sum_{j = 1}^{n}(-c_{n - j} X_{n - j})$

이렇게 증명이 완성된다. 뒤에서 사용하는 테크닉도 대부분 이런 식이니 한번 손으로 해 보자.

3. Calculate minimal polynomial of $A$

주제를 바꿔서 이제부터는 DP 이야기가 아니라, 행렬에 관한 다양한 연산을 Berlekamp-Massey를 사용해서 해결하는 방법을 설명한다. 여기부터는 Berlekamp의 오리지널 아이디어가 아니라, 1986년 Wiedemann이 개발한 아이디어를 더한 것이며, 또한 랜덤 알고리즘을 사용한다.

Berlekamp-Massey에서 활용하는 것은 행렬의 sparse 함이다. 위에서 다룬 점화식들도 꽤 많은 경우 sparse matrix가 나왔기 때문에 (각 DP 상태에 대해서 전이가 상수 개였기 때문에) 효율적인 알고리즘을 만들 수 있는 경우가 많았다. 그 외 그래프의 인접 행렬들도 sparse matrix이다. 이러한 행렬에서의 다양한 연산을 Berlekamp-Massey로 최적화할 수 있다. 이제부터, 어떠한 행렬의 0이 아닌 값의 개수를 $m$이라고 한다.

어떠한 행렬에 대한 최소 다항식 은 $\sum_{i=0}^{k}{c_i A^i} = 0$ 를 만족시키는 자명하지 않은 다항식 중 가장 작은 것이다. 위에 있는 특성 다항식과 비슷하지만, 크기를 최소화한다는 점이 다르다. 최소 다항식을 구하면 여러 계산 문제들을 빠르게 해결할 수 있는데, 이에 대해서는 나중에 설명하고, 일단 이 최소 다항식을 어떻게 구하는지 생각해보자.

위에서 Cayley-Hamilton 정리를 도입하면서 특성 다항식의 개념을 언급하고, 특성 다항식에 대응되는 선형 점화식이 매우 비슷함을 (그냥 항이 일대일 매칭이 되는 수준) 알 수 있었다. 실제로는 Berlekamp-Massey가 최소 크기 선형 점화식을 찾기 때문에, 이 최소 크기 선형 점화식을 사용하면 최소 다항식을 찾을 수 있다. 고로, 적당한 $v, p_0$ 이 존재하면, 2번에서 한 과정을 그대로 따라 해서 선형 점화식을 찾고, 이를 간단히 최소 다항식으로 바꿀 수 있다.

하지만, 우리는 DP를 해결하기 위해서 문제를 푼 것이 아니니 $v, p_0$ 로 잡을 것이 없다. 그러니 아무거나 잡아 보자. $v, p_0$ 을 적당한 랜덤 벡터로 잡고, $X_i = v \times (A^k \times p_0)$ 을 계산해 준 후 Berlekamp-Massey에 대응시키면 다항식을 구할 수 있다. 이에 대한 확률 분석이 매우 복잡하나, 하여튼 매우 높은 확률로 구할 수 있다.

여기서 중요한 컨셉이 있는데, $A^k$ 를 구하면 안 된다. sparse matrix여도 제곱을 했을 경우 빠르게 그 sparse 함을 잃기 때문에, 이를 일일이 곱해주면 $O(n^4)$ 시간 복잡도가 든다. 반면, 어떠한 행렬에 벡터를 곱하는 것은 $O(m)$ 시간밖에 안 든다! 고로, $A^{k-1} \times p_0$ 라는 벡터에 $A$ 를 곱하는 것을 $n$ 번 반복하는 식으로 해야지 효율적이다. 이때의 시간 복잡도는 $O(nm)$이다. 이러한 류의 트릭은 Frievald's algorithm에서도 사용되는 방법이다. 사실 위에서의 DP 개념을 가지고 생각해 보면, 매우 당연한 말이다.

예시 1: $Ax = b$ 빠르게 풀기

이제 최소 다항식을 구했으니 이것으로 뭘 할 수 있는지를 소개한다. $Ax = b$라는 행렬식을 풀어 보자. 일반적으로 이는 가우스 소거법을 사용해서 $O(n^3)$ 에 구할 수 있으나, 이 방법을 사용하면 sparse matrix에서 $O(nm)$ 에 해결할 수 있다.

방법은 아주 간단하다. 그냥 $x = A^{-1} b$ 를 구하는 것이다. 최소 다항식을 $A^k + c_{k-1} \times A^{k-1} + \cdots + c_0 \times A^0 = 0$이라 두자. 적당히 전개해 주면.

$A^{k-1} + c_{k-1} \times A^{k-2} + \cdots + c_1 \times A^{0} + c_0 \times A^{-1} = 0$

$A^{k-1} + c_{k-1} \times A^{k-2} + \cdots + c_1 \times A^{0} = -c_0 \times A^{-1}$

$(-c_0)^{-1}(A^{k-1} + c_{k-1} \times A^{k-2} + \cdots + c_1 \times A^{0}) = A^{-1}$.

$(-c_0)^{-1}(A^{k-1} + c_{k-1} \times A^{k-2} + \cdots + c_1 \times A^{0}) \times b = A^{-1} \times b$.

고로, 위에서 설명한 대로 최소 방정식을 구한 후, $A^k \times b$ 역시 위에서 설명한 대로 구하면, $O(nm)$ 에 쉽게 해결할 수 있다. $A^{-1}$을 구하는 것은 어렵지만, $A^{-1}\times b$ 를 구하는 것은 쉽다는 것이 흥미롭다. Frievald's algorithm과 비슷한 패러독스.

예시 2: $A^k \times b$ 빠르게 계산하기

이것 역시 최소 다항식을 구했으면 어렵지 않다. $A^k$ 를 최소 다항식으로 나눈 나머지를 repeated squaring으로 계산하고, 결과적으로 나온 $k-1$ 차 이하의 다항식에 $v$ 를 곱해서 계산하면 되기 때문이다. 키타마사법에서 사용한 트릭과 동일하다 :smile: 시간 복잡도는 $O(nm + n^2 \log k)$이다.

4. Find determinant of matrix $A$

마지막 트릭은 행렬의 determinant를 구하는 것이다. 행렬의 determinant도 가우스 소거법을 사용하면 $O(n^3)$ 에 해결할 수 있지만, 여기서는 $O(nm)$ 에 해결할 수 있다. 이것이 함의하는 내용 중 아마 가장 흥미로운 것은, 그래프의 스패닝 트리의 개수를 $O(|V||E|)$ 에 계산할 수 있다는 것이 되겠다.

행렬의 determinant를 구하기 위해서 우리는 행렬의 특성 다항식을 구한다. 특성 다항식을 구하면 행렬의 determinant를 구하기는 매우 쉽다. 특성 다항식의 정의가 $p(\lambda) = det(\lambda I_n - A)$ 이기 때문에, $p(0) = det(-A) = det(A) \times (-1)^n$ 이 되기 때문이다. 즉 그냥 상수항만 떼어와서 $n$의 parity에 따라 뒤집어주면 된다. 하지만, 우리는 행렬의 최소 다항식 만 구할 줄 알지 특성 다항식을 구하는 방법은 알지 못한다.

여기서 약간 허무한 (…) 트릭을 사용한다. 아무 랜덤 대각 행렬 $M$을 잡고, $A\times M$ 의 최소 다항식을 구해주자. $M$ 은 대각 행렬이기 때문에 $AM$을 계산하는 데 $O(m)$ 의 시간밖에 들지 않는다. 이때, $AM$의 최소 다항식은 매우 높은 확률로 $AM$의 특성 다항식이다. 정확한 증명을 전혀 모르기 때문에 (…) 생략하지만, 어떠한 행렬의 최소 다항식과 특성 다항식이 다르다는 것은 특성 다항식의 해 (고윳값) 이 중복되는 것들이 많다는 것을 뜻하고, 고로 이러한 식으로 랜덤화를 하면서 고윳값을 "섞어주는" 역할이라고 직관적으로 이해하면 된다. $det(A)det(M) = det(AM)$ 이 성립하고, $det(M)$ 은 그냥 대각 원소의 곱이니, determinant를 구할 수 있다.