행렬 곱셈 문제 - 알고리즘카

선형 함수와 행렬에 대해 이야기할 때 모든 예시는 기하학에 관한 것이었습니다. 그렇게 하는 편이 그냥 생각하기 더 쉽기 때문입니다. 평면 위 벡터의 회전과 늘이기를 머릿속으로 떠올리는 것이 추상적인 n차원 공간의 성질을 따지는 것보다 쉽습니다. 그러나 우리는 프로그래머이고, 주로 관심을 두는 것은 기하학 밖의 문제입니다.

올림피아드 맥락에서 대부분의 응용은 행렬 곱셈의 결합 법칙에 기반합니다:

$$ ABC = (AB)C = A(BC) $$ 문제를 행렬 $A$를 어떤 큰 거듭제곱 $n$으로 올리는 것, 즉 단위 행렬에 행렬 $A$를 $n$번 곱하는 것으로 환원할 수 있다면, 결합 법칙을 이용하여 이진 거듭제곱으로 결과를 계산하면 됩니다: $$ A^8 = A^{2^{2^{2}}} = ((AA)(AA))((AA)(AA)) $$

이러한 접근 방식은 동적 계획법, 조합론, 확률론, 그리고 수학 전반에 걸쳐 다양하게 응용됩니다.

#선형 점화식

문제. 수 $n < 10^{18}$이 주어집니다. $n$번째 피보나치 수를 계산해야 합니다:

$$ \begin{aligned} f_0 &= 0 \\ f_1 &= 1 \\ f_{n} &= f_{n-1} + f_{n-2} \end{aligned} $$

다음 피보나치 수를 계산하는 것은 직전 두 피보나치 수에 대한 선형 함수, 즉 두 수의 합으로 표현된다는 점에 주목합시다.

이는 이 변환에 대응하는 $2 \times 2$ 행렬을 만들 수 있다는 것을 의미합니다. 즉, 이웃한 두 피보나치 수 $(f_n, f_{n+1})$로부터 다음 수를 계산하여 쌍 $(f_{n+1}, f_{n+2})$로 넘어가는 변환입니다.

$$ \begin{pmatrix} f_{n+1} \\ f_{n+2} \\ \end{pmatrix} = \begin{pmatrix} 0 + f_{n+1} \\ f_{n} + f_{n+1} \\ \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ 1 & 1 \\ \end{pmatrix} \begin{pmatrix} f_{n} \\ f_{n+1} \\ \end{pmatrix} $$

이 전이 행렬을 $A$라고 합시다. $n$번째 피보나치 수를 계산하려면 이 행렬을 벡터 $(f_0, f_1) = (0, 1)$에 $n$번 적용해야 합니다.

행렬 $A$의 크기는 상수이므로, 행렬을 자기 자신과 곱하는 곱셈은 $O(1)$에 동작합니다. 따라서 이진 거듭제곱으로 $A^n$을 $O(\log n)$번의 연산으로 계산한 뒤, 벡터 $(0, 1)$을 곱하고 결과의 첫 번째 원소를 취하면 됩니다. 이 원소가 바로 우리가 찾던 $f_n$과 정확히 같습니다.

참고. 벡터를 명시적으로 곱할 필요조차 없습니다. 첫 번째 원소만 필요하다면 최종 결과에 대한 식을 살펴보고, 최종 행렬의 오른쪽 위 원소를 취하기만 하면 충분하다는 것을 알 수 있습니다.

#일반적인 경우

일반적인 경우, 즉 마지막 $2$개가 아니라 $k$개의 원소를 서로 다른 가중치로 고려할 때, 선형 점화식 $f_n = a_1 f_{n-1} + a_2 f_{n-2} + \ldots + a_k f_{n-k}$의 전이 행렬은 다음과 같습니다:

$$ \begin{pmatrix} 0 & 1 & 0 & \ldots & 0 \\ 0 & 0 & 1 & \ldots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \ldots & 1 \\ a_k & a_{k-1} & a_{k-2} & \ldots & a_1 \\ \end{pmatrix} $$

벡터 $(f_n, f_{n+1}, \ldots, f_{n+k-1})$에 곱하면 마지막 원소에 대해서는 정확히 정의의 식이 얻어지고, 앞의 $(k-1)$개 원소는 모두 그대로 복사됩니다.

점화 수열은 계수 $a_i$뿐만 아니라 초기값 $(f_1, f_2, \ldots, f_k)$에 의해서도 결정된다는 점에 유의하십시오.

연습 문제. $k$차 선형 점화식에 대해 그 계수 $a_i$와 $(i, f_i)$ 꼴의 쌍 $k$개가 주어집니다. 수열의 처음 $k$개 원소를 복원하십시오.

#등비수열

이진 거듭제곱에 관한 글에서 우리는 등비수열의 합을 구하기 위한 알고리즘의 변형을 다루었습니다.

$$ f(n) = 1 + a + a^2 + \ldots + a^{n-1} $$ 이 문제에 다른 방향으로 접근할 수도 있습니다. 처음 $n$개의 거듭제곱의 합은 처음 $(n-1)$개의 거듭제곱의 합을 이용해 다음과 같이 표현됨에 주목합시다: $$ f(n) = f(n-1) \cdot a + 1 $$ 이는 이전 항에 대한 거의 선형적인 의존 관계입니다. 다만 불편한 상수가 하나 더 붙을 뿐입니다. 다음과 같은 트릭을 쓸 수 있습니다. 우리가 $(n-2)$번째 원소에도 의존한다고 가정하되, 그 값을 항상 1로 고정하는 것입니다. $$ \begin{pmatrix} f_{n+1} \\ 1 \\ \end{pmatrix} = \begin{pmatrix} f_n \cdot a + 1 \\ 0 + 1 \\ \end{pmatrix} = \begin{pmatrix} a & 1 \\ 0 & 1 \\ \end{pmatrix} \begin{pmatrix} f_{n} \\ 1 \\ \end{pmatrix} $$

그러면 이 행렬도 마찬가지로 $n$제곱한 뒤 오른쪽 위 원소를 반환하면 됩니다.

#동적 계획법

일부 DP의 전이 역시 때로는 행렬 곱셈의 관점에서 표현할 수 있습니다.

문제. 크기가 $n < 500$인 방향 그래프 $G$가 주어집니다. $s$에서 $t$로 정확히 $k < 10^{18}$번의 전이를 거치는 경로의 수를 구해야 합니다.

만약 $k$에 대한 제약이 훨씬 작았다면, 정확히 $i$번의 전이를 거쳐 $v$번째 정점에 도달하는 방법의 수와 같은 DP $f_{i,v}$를 도입했을 것입니다. 점화식은 경로에서 끝에서 두 번째 정점을 전부 시도해 보는 것입니다:

$$ f_{i,v} = \sum_{u \in g_v} f_{i-1,u} $$ 이 DP를 인접 행렬의 관점에서 다시 써 봅시다. 인접 행렬의 칸 $G_{uv}$에는 $u$와 $v$ 사이의 간선 개수가 들어 있습니다(바로 이 때문에 행렬이라고 부릅니다). $$ f_{i,v} = \sum_u f_{i-1,u} \cdot G_{uv} $$ 벡터 $f_{i}$는 $f_{i-1}$에 선형적으로 의존한다는 것을 알 수 있습니다. $f_{i}$의 $v$번째 원소는 벡터 $f_{i-1}$과 인접 행렬 $G$의 $v$번째 열의 내적으로 얻어집니다. 따라서 다음 등식이 성립합니다. $$ f_i = G \cdot f_{i-1} $$

결국 $f_k$는 $f_k = f_0 \cdot G \cdot G \cdot \ldots \cdot G$와 같이 전개할 수 있으므로, 이진 거듭제곱을 적용할 수 있습니다.

이 공식에 따르면 $n \times n$ 크기의 두 행렬을 $O(\log k)$번 곱해야 하므로, 전체적으로 $O(n^3 \log k)$의 시간이 걸립니다.

const int n;

matrix binpow(matrix a, int p) {
    // создадим единичную матрицу
    matrix b(n, vector<int>(n, 0));
    for (int i = 0; i < n; i++)
        b[i][i] = 1;

    while (p > 0) {
        if (p & 1)
            b = matmul(b, a);
         a = matmul(a, a);
         p >>= 1;
    }

    return b;
}

결과로 얻은 행렬의 칸 $G_{ij}$에는 정확히 $k$번의 이동을 사용하여 $i$번째 정점에서 $j$번째 정점으로 가는 방법의 수가 들어 있게 됩니다. 답으로는 단순히 $G_{st}$를 출력하면 됩니다.

#문제의 변형

약간만 변형하면 이 방법으로 비슷한 문제를 많이 풀 수 있습니다:

알고리즘 자체를 거의 바꾸지 않고도, 인접 행렬 대신 마르코프 연쇄에서 한 정점에서 다른 정점으로 이동할 확률이 주어진 경우 “정점 $s$에서 정점 $t$로 도달할 확률은 얼마인가"라는 문제를 풀 수 있습니다.
경우의 수가 아니라 정확히 $k$번의 전이로 도달할 수 있는지 여부만 필요하다면, 행렬을 비트셋으로 감싸서 풀이를 크게 가속할 수 있습니다.
“최대 $k$번의 전이로"라고 묻는 경우에는 행렬의 $k$제곱 대신 앞서 설명한 방법으로 등비수열의 합을 계산할 수 있습니다. 또는 각 정점마다 두 번째 정점을 추가하여, 원래 정점에서 그 정점으로 가는 간선을 두고 그 정점에서 나가는 유일한 간선은 자기 자신으로 향하는 루프가 되도록 할 수도 있습니다.

이 기법은 무언가를 하는 방법의 수를 세어야 하는 다른 DP에도 적용할 수 있으며, 때로는 매우 자명하지 않은 방식으로 적용됩니다.

예를 들어, 다음과 같은 문제를 풀 수 있습니다: 주어진 짧은 금지 부분 문자열들을 포함하지 않는 길이 $k \approx 10^{18}$인 문자열의 개수를 구하십시오. 이를 위해서는 아호-코라식에서 “허용되는” 전이들로 이루어진 그래프를 만들고, 그 인접 행렬을 $k$제곱한 뒤 첫 번째 행의 값을 모두 더하면 됩니다.

일부 정교한 경우에는 행렬 곱셈에서 곱셈과 덧셈 대신 곱셈과 덧셈처럼 동작하는 다른 연산을 사용해야 합니다. 문제 예시: “정확히 $k$번의 이동을 사용하여 $s$에서 $t$까지 가는, 간선 가중치가 최소인 경로를 찾으십시오”. 여기서는 그래프의 가중치 행렬을 $(k-1)$제곱해야 하며, 덧셈과 곱셈 모두 대신 두 가중치 중 최솟값을 사용해야 합니다.

또한 최단 경로를 찾는 맥락에서는 «distance product»가 알려져 있습니다:

$$ C_{ij} = \min_k \{ a_{ik} + b_{kj} \} $$

가중치 행렬을 $k$제곱하면, 그에 따라 $k$개 이하의 간선을 사용하여 $a$에서 $b$까지 가는 최소 거리를 얻습니다. 특히 $k=n$일 때 $O(n^3 \log n)$번의 연산으로 모든 정점 쌍 사이의 최단 거리 행렬을 얻게 되지만, 이를 위한 더 직접적인 방법도 있습니다.