이 글에서는 Exponentiation By Squaring을 이용하여 정사각행렬 A의 N 제곱을 구하는 방법을 알아보도록 하겠습니다.
1. Identity Matrix of Order N (한: N차 단위 행렬)
Matrix라는 대상은 C/C++ 환경에서 보통 array 또는 vector로 구현됩니다.
int array[N][N];
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
if (i == j) {
array[i][j] = 1;
}
else {
array[i][j] = 0;
}
}
}
(tutorial for getting an identity matrix)
위 코드에서 array를 바로 N차 단위 행렬이라고 할 수 있습니다. 즉, identity matrix는 diagonal entries가 모두 1이고, 그걸 제외한 나머지 entries가 0인 square matrix입니다.
N차 단위 행렬 E와 N차 정사각행렬 A에 대해 다음 성질이 성립합니다.
$$ AE = EA = A $$
이 성질과, Exponentiation By Squaring의 응용을 통해 우리는 정사각행렬의 K제곱을 빠르게 구할 수 있습니다.
편의상 앞으로 사용할 Matrix 구조체를 소개하겠습니다.
typedef struct {
element array[size][size];
} Matrix;
2. 행렬곱
ko.wikipedia.org/wiki/%ED%96%89%EB%A0%AC_%EA%B3%B1%EC%85%88
행렬 곱셈 - 위키백과, 우리 모두의 백과사전
위키백과, 우리 모두의 백과사전. 둘러보기로 가기 검색하러 가기 행렬 곱셈을 위해선 첫째 행렬의 열 갯수와 둘째 행렬의 행 갯수가 동일해야한다. 곱셈의 결과 새롭게 만들어진 행렬은 첫째
ko.wikipedia.org
위키피디아의 2.정의 목차에 나와있는 공식대로, 두 정사각행렬 A, B의 행렬곱을 리턴해주는 함수를 짜보겠습니다.
Matrix matrix_multiply (Matrix A, Matrix B, int N) {
// N차 정사각행렬 A, B의 행렬곱
Matrix result;
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
result[i][j] = 0;
for (int k = 0; k < N; k++) {
result[i][j] += A[i][k] * B[k][j];
}
}
}
return result;
}
이렇게 짜면 시간 복잡도가 $O\left ( N^{3} \right )$이라는 걸 알 수 있습니다.
영문 위키피디아에 따르면 $O\left ( N^{2.3728639} \right )$ 복잡도의 알고리즘이 존재한다고 하지만, Constant의 영향으로 strassen마저도 퍼포먼스가 나이브한 $O(N^{3})$ 알고리즘에 밀리는 편입니다.
또 위 함수는 locality라는 개념을 활용하여 다음과 같이 최적화를 한 번 더 할 수 있습니다. 이에 관해서는 구글 검색, 또는 컴퓨터 구조 강의를 듣거나 CS:APP 같은 서적을 참고할 수 있겠습니다.
Matrix matrix_multiply (Matrix A, Matrix B, int N) {
// N차 정사각행렬 A, B의 행렬곱
Matrix result = {{{}}};
for (int i = 0; i < N; i++) {
for (int k = 0; k < N; k++) {
for (int j = 0; j < N; j++) {
result[i][j] += A[i][k] * B[k][j];
}
}
}
return result;
}
3. N차 정사각행렬의 K제곱
이 글에서
multiply_modular 함수의 경우 result를 0으로, power_modular 함수의 경우 result를 1로 초기화 했었습니다.
그 이유는 파라미터로 0이 들어갔을 때의 예외처리를 위해서 입니다.
multiply_modular에서 $x \times 0=0$이고 power_modular에서는 $x^{0}=1$이기 때문입니다. (단, $0^0=1$ convention은 여기서 언급하지 않겠습니다.)
마찬가지로, 우리가 이 글에서 구현할 함수도 예외처리가 필요합니다.
어떤 행렬과 영행렬의 곱은 항상 영행렬입니다. 또한, 정수의 0제곱이 1인 것과 같이 정사각행렬의 0제곱은 단위 행렬입니다. 이를 고려하여 다음과 같이 함수를 짤 수 있습니다.
Matrix matrix_multiply (Matrix A, Matrix B, int N) {
// N차 정사각행렬 A와 B의 행렬곱
Matrix result = {{{}}};
for (int i = 0; i < N; i++) {
for (int k = 0; k < N; k++) {
for (int j = 0; j < N; j++) {
result[i][j] += A[i][k] * B[k][j];
}
}
}
return result;
}
Matrix matrix_power (Matrix A, int N, int K) {
// N차 정사각행렬 A의 K제곱
Matrix result = {{{}}};
for (int i = 0; i < N; i++) {
result[i][i] = 1;
}
while (K > 0) {
if (K & 1) {
result = matrix_multiply (result, A, N);
}
A = matrix_multiply (A, A, N);
K >>= 1;
}
return result;
}
시간복잡도는 matrix_multiply 함수가 $O\left ( N^{3} \right )$이므로, matrix_power 함수는 $O\left ( N^{3}\times log K\right )$ 입니다.
이때 조금만 K가 커져도 행렬의 각 원소가 int 또는 long long의 범위를 벗어나기 십상입니다. 따라서 행렬을 N제곱하여 푸는 문제들은 각 원소를 어떤 숫자 M으로 나눈 값을 필요로 하는 경우가 많습니다. 이에 관한 글은 따로 작성해보도록 하겠습니다.
여담으로, 이 알고리즘은 그 자체가 어려운 게 문제가 아니고,
문제를 봤을 때 어떤 행렬을 N제곱하면 해결할 수 있는지 찾는 게 더 큰 문제인 것 같습니다.
'이론' 카테고리의 다른 글
PS에서 쓰는 Matrix multiplication의 구현과 최적화 (0) | 2023.01.30 |
---|---|
기본 이론 3 (0) | 2021.07.31 |
기본 이론 (0) | 2021.01.23 |