1. Introduction
Problem Solving을 위해 C/C++ 을 쓰는 사람들에게 matrix multiplication의 구현은 귀찮은 일입니다. 당연히 많은 사람들이 이를 템플릿에 박아놓고 그대로 씁니다. 하지만, 문제의 제한에 따라 다양한 최적화의 여지가 있으며 일부는 탁월한 효과를 보입니다. 이 글에서는 가장 naive한 구현에서부터 특수한 맥락에서만 쓸 수 있는 코드까지 소개해보고자 합니다. 대부분의 코드는 제 경험으로부터 비롯하였으며, 따라서 C++의 reference type이나 STL(특히, <vector>, <pair> 등)을 쓰거나, class로 wrapping하는 코드는 다루지 않습니다.
2. 기본적인 구현
우선, matrix multiplication의 가장 일반적인 형태는 다음과 같습니다.
$m \times n$ 행렬 A와 $n \times p$ 행렬 B의 곱 C = AB의 각 원소는
$$ C_{ij} = \sum_{k=1}^{n} A_{ik} B_{kj} $$
(단, $m,n,p \in \mathbb{N}$이며 이 정의에서 행렬은 1-index를 씁니다.)
그리고 이걸 그대로 구현하면 흔히 다음의 3중 for문과 같은 형태가 나옵니다.
(Create image from code: Carbon )
이런 코드를 실제로 제가 썼었고, 또 많은 사람들이 흔히 그렇게 씁니다. 여기서 어떤 걸 고칠 수 있을까요? 차근차근 살펴봅시다.
3. PS적인 접근 (1) square matrix
일반적으로 Problem Solving에서 다루는 행렬 곱은 square matrix, 정사각행렬 간의 곱셈입니다. 그렇지 않더라도 일단 $\textrm{max} (m,n,p) \times \textrm{max} (m,n,p)$ 정사각행렬끼리 곱하고 그 중 일부 entries의 값만 취해도 됩니다. 그래서 위 사진의 코드는 다음과 같이 바뀔 수 있습니다.
물론 이것 자체는 naive하긴 해도 흔히 쓰이는 코드입니다. 하지만, loop 안에서 하나하나 result[i][j]의 초기값을 지정해줘야할까요?
4. "result" Matrix의 모든 원소를 0으로 초기화하는 방법 (1) Memset
result 배열을 초기화하는 것은 2중 루프를 써서 할 수도 있지만... 저희가 사용하고 있는 2-dimensional arrays 는 실제 메모리 상으로는 1-dimensional array 와 같은 형태입니다. 그러므로, 우리는 <string.h>의 함수 "memset"을 쓸 수 있습니다. 대략 다음 사진과 같은 느낌입니다.
여기에서 과연 memset이 for loop보다 효율적일까? 라고 생각할 수 있는데, 실제로 효율적입니다!
5. "result" Matrix의 모든 원소를 0으로 초기화하는 방법 (2) Declaration
만약 result 배열이 행렬곱 루틴 내부에서만 쓰이는 로컬 변수라면 어떤 방법을 쓸 수 있을까요?
C99 표준안을 살펴보면, 배열을 선언할 때 그 원소의 개수를 변수의 값에 따르는 "Variable Length Arrays"가 명시되어 있습니다. 즉, BOJ에서 C99를 쓴다면 VLA를 쓰는 것이 자연스럽습니다. (그러나 로컬에서 Microsoft Visual C++ Compiler를 쓴다면, 이 VLA 기능이 지원되지 않습니다.)
C11부터는 VLA는 표준에서 빠집니다. 그러나 gcc에서는 여전히 VLA를 지원하고 있고, BOJ에서는 어차피 gcc를 써서 채점하니까 VLA를 써도 상관 없습니다. (아래 result 배열을 선언하는 부분은 assembly 로 번역됐을 때 calloc과 같이 처리될 수 있습니다. 다만, VLA는 stack에 쌓이고 일반적인 local variable처럼 선언된 scope가 끝나면 소멸됩니다. )
6. struct의 도입 (Optional)
이 부분은 optional하지만, 그렇다고 해도 이후에 소개할 코드의 가독성이 더 나아지는 것 같고, trade-off가 매우 크지도 않아서 꽤 괜찮은 것 같아 소개합니다.
저는 보통 elem 자리에 int, long long, unsigned 자료형, 그리고 bool을 넣는 편입니다. size는 따로 define된 매크로를 쓰기보다, 그냥 5, 100, 166, 256 같은 상수를 쓰는 편입니다. 그리고 제가 일반적으로 쓰는 struct를 쓰는 matrix multiplication 코드입니다.
이렇게 struct을 쓸 때의 장점은, pointer, 혹은 memset, memcpy, calloc, ... 로 가득한 함수들을 보지 않아도 됩니다. 물론 템플릿을 쓰는 사람이라면 그런 복잡한 함수가 큰 단점은 아니겠지만요. 대신 Trade-Off 로써 후술할 'Call by Value'에 관한 이슈가 있습니다. (의외로 Trade-Off가 크진 않습니다.)
7. Improve Locality
C/C++는 2D array를 메모리에 Row-Major 방식으로 저장합니다. 그래서 Row-Major 방식으로 순회해야 빠릅니다. 이걸 고려했을 때 행렬곱에서 쓰이는 3중 루프, 이를테면 i-loop, j-loop, k-loop을 어떻게 배치해야 가장 빠를까요? 놀랍게도 $i-j-k$는 6가지의 순서 중에서 중간 수준에 속합니다. 그래서 제한 시간이 널널하거나 행렬이 작은 경우라면 개선의 필요성을 잘 못 느끼기도 하죠.
통상 loop의 순서는 $k-i-j$ 혹은 $i-k-j$ 가 최고라고 합니다. 저는 일반적으로 $i-k-j$을 씁니다만, $k-i-j$을 써도 무방합니다.
8. Modular Matrix Multiplication Problem
Modular matrix multiplication 을 구현해보았습니다. 일단, matrix의 element는 편의상 int로 잡았습니다. 이후 챕터는 아래 함수의 최적화에 관한 내용입니다.
9. Remove the Temporary Variables
temp의 존재 자체가 가독성을 헤치는 느낌도 들고, 일반적인 matrix multiplication에서 반드시 필요한 것도 아니니까 없애보겠습니다.
10. Simplify the Expression
굳이 innermost loop에서 두 번씩 assign 해줄 필요가 전혀 없습니다. 한 번으로 줄일 수 있습니다.
그런데, 위 함수를 살펴보면 innermost loop에서 A와 B의 각 원소를 modulo를 취하고, long long으로 casting해서 둘을 곱하고, 다시 result[i][j]와 더해서, 이 값을 modulo 취한 것을 대입하고 있습니다. 이렇게까지 modulo를 취해야 하나요?
11. Minimize the Number of Modulo Operation
modulo operation은 C에서 제공하는 arithmetic logical operation 중에서 가장 무겁습니다. 꼭 Intel의 Xeon (BOJ의 채점환경)만이 아니라도 대부분의 CPU에서 해당됩니다. 이 modulo operation의 횟수를 줄이는 것은 가장 쉽게 생각할 수 있는 최적화 방안 중 하나입니다. 이제 다음을 가정합니다.
$$ |A_{ij}|, |B_{ij}| < \textrm{mod} $$
그러면, 다음과 같이 식이 간소화됩니다.
이렇게 되면 innermost loop 내부의 % 연산이 4회에서 2회로 줄었습니다. 하지만 특수한 상황에선 더 줄일 수 있습니다. 예컨대 mod가 $10^{9}+7$이라고 가정하면,
$$ mod < mod^{2} + mod < 2^{63} - 1$$
이므로, $A * B + \textrm{res}$ 를 한 번에 계산한 뒤에 modulo를 취해도 됩니다. (물론 casting을 먼저 해야 합니다.) 그래서 innermost loop의 각 iteration 당 modulo operation의 횟수가 1번으로 줄게 됩니다.
12. Use Larger Types (Optional)
챕터 8에서 array의 각 element을 int로 잡은 바 있습니다. 그런데, $mod \geq 10^{9}$ 라면 long long으로의 casting이 불가피합니다. 따라서 array의 각 element을 long long으로 미리 선언하면, 코드가 깔끔해질 수 있습니다.
일반적으로, casting 없이 행렬곱을 쓰려면 $mod^{2} + mod$ 보다 최댓값이 큰 자료형을 쓰는 것이 안전합니다. 다만, 그렇다고 해도 __int128은 느린 편이므로 가능한 안 써야 합니다. 또한 지나치게 딱 맞게 char나 bool, short 등으로 잡을 필요는 없습니다. 일반적으로 int가 어떤 시스템이든 제일 빠릅니다. 굳이 가장 빠른 자료형을 찾고 싶다면, C언어 표준 헤더인 <stdint.h> 을 참고해볼 수도 있겠습니다.
위와 같이 각 원소가 long long으로 커진 행렬에 따라 아래와 같이 식이 바뀝니다. 저는 최근엔 보통 이 형태를 씁니다.
물론, 아직도 최적화의 여지는 더 남았습니다!
13. Use Constant Operands
일반적으로 생각하기 어려운 점이지만, % 연산에 대응하는 x86 assembly code을 찾아보면, operand가 variable일 때와 constant일 때의 속도 차이 가 존재한다고 합니다. ("Barrett Reduction") 결론은 constant로 operand를 주는 게 좋다는 건데, 상수에 이름을 붙이는 방법으로 C에서 제공하는 건 define과 const 키워드 두 방법이 있습니다. 제 스타일은 코드 제일 위에 define 으로 정의해놓는 것입니다.
그리고 이건 사족이지만, 개인적으로 parameter-list가 긴 게 싫기도 하고, 어차피 행렬의 크기가 뻔한데 굳이 size을 넘겨줄 필요가 없는 것 같아서, size 인자도 상수로 대체합니다. 물론 속도랑은 큰 연관이 없습니다.
14. Adhere Unsigned Integers
챕터 14의 타겟은 답을 0 이상의 정수로 출력해야 하는데, 행렬에(점화식에) 음수 값이 존재하는 경우입니다. C의 modulo operation은 그 resultant가 음수값이 될 수 있으므로 생기는 사태입니다.
일반적으로 이를 위해서 if문을 쓰는 방법이 있습니다. 왜냐하면, 어차피 resultant가 음수라고 해도, 그 절댓값은 mod 이하이기 때문입니다.
만약 불안하다면, if문을 while로 바꾸는 경우도 있습니다. 하지만, 둘 다 필요없이 mod를 더해주고 다시 % 연산을 하면 충분합니다.
그리고 여기서 한 번 더 modulo operation 횟수를 줄일 수 있습니다. 현재 $mod=10^{9}+7$이고 $(10^{9}+9)(10^{9}+7) < 2^{63}-1$ 임을 염두하면, 대략 다음과 같은 코드가 나옵니다.
다만 저는 이런 방법보다는, 그냥 행렬곱에서 아무 조치를 취하지 않고 정답 출력 직전에 분기를 파거나, 아니면 점화식 자체의 음수를 없애는 편입니다. 점화식의 음수를 없애버리면 innermost loop의 식을 바꿀 필요가 없기 때문에...
15. Use Unsigned Types
unsigned 자료형을 쓸 수 있으면 쓰는 게 낫습니다. 점화식에 음수가 있더라도 챕터 14의 내용을 활용할 수 있습니다. 하지만, implicit type conversion에 각별히 신경써주세요. 괜히 signed와 unsigned을 섞으면 결과가 좋지 않을 수 있습니다.
16. Exploiting the IJK loop: Minimizing modulo operations
이 팁은 아주 극단적인 최적화가 필요할 때 쓸 수 있는 한 가지 방법입니다. 행렬의 차원 $N$과 각 entries의 자료형의 최댓값, $MAX$ 에 대해 다음이 성립할 때 쓸 수 있습니다.
$$ N \times mod^{2} < MAX $$
이런 상황에서, innermost loop에서 매번 modulo operation을 수행하는 것은 살짝 비효율적입니다. 즉, 다음과 같이 코드를 개선할 수 있습니다.
일반적으로 행렬이 N by N일 때, 원래 modulo operation은 $O(N^{3})$ 번 실행됐으나, 이제 그 차수가 하나 낮아졌음을 알 수 있습니다. 단점이라면 C에서는 MAX값이 너무 작은데, unsigned 자료형을 써서 MAX값을 조금이라도 더 키울 수는 있습니다.
17. A Minor Tip for Boolean Matrix Multiplication: Branch Prediction
Boolean Matrix는 각 entries가 true or false의 Boolean 값을 갖는 matrix로, 이런 Boolean Matrix끼리의 multiplication은 logical operations, 특히 logical AND와 OR에 의해 정의됩니다.
C에서는 <stdbool.h> 를 include하여 보통 Boolean 값 중 true에 1을, false에 0을 넣습니다. 이는 앞서 언급해온 modular matrix multiplication에서 mod가 2인 상황과 analogous합니다.
저는 이렇게 $\mathbb{Z}/2\mathbb{Z}$에서는 bitwise operation을 씁니다. (검색해보니 아쉽게도 arithmetic을 쓰느냐, bitwise를 쓰느냐는 큰 차이가 없는 듯합니다.)
그런데, 잘 생각해보면 A[i][k]가 false라면 굳이 innermost loop으로 들어갈 필요가 없습니다. 만약 각 entries의 값이 충분히 predictable하다면 이는 강력한 보조를 해주겠죠? 또한, innermost loop에서 A[i][k]가 true임이 보장됐는데도 굳이 bitwise AND를 거칠 필요가 없습니다!
18. Loop Unrolling For Extremely Small Matrices
이 챕터에서 소개하는 방법은 $2 \times 2$ 내지는 $4 \times 4$ 수준의 매우 작은 matrices에 관한 것입니다. $3 \times 3$ 행렬을 생각해보면, 최소 27번의 logical operation과 increment operation이 수행되고 있습니다. 물론 branch prediction이 이런 걸 잘 보조해주겠지만, 애초부터 loop가 없었다면 필요 연산 횟수가 급격히 감소할 것입니다. 가령, $2 \times 2$ 행렬곱은 다음처럼 수행될 수 있습니다.
19. Issues for using struct
앞서 챕터 6에서 optional하게 struct를 써도 된다고 했는데, struct 사용과 관련하여 어떤 이슈가 있는지를 알면 더 좋겠죠?
흔히 struct를 parameter로 받고 return하는 코드가 안 좋은 습관이라고 합니다. 근데 저는 Call-by-Value가 큰 문제는 아니라고 생각하는게, argument나 return value를 struct로 통째로 넘기고 받는 건 assembly 단계에서 memcpy 와 비슷하게 흘러가는 것 같습니다. 그런데 memcpy도 단순 for-loop보다는 훨씬 빠른지라 이런 부분이 크리티컬하게 작용하긴 힘듭니다.
메모리는 PS 환경에서 크리티컬할 일이 거의 없습니다. 행렬 거듭제곱에 관한 문제는 그 시간 복잡도상 대부분이 $N \leq 100$, 크게 잡아도 $N \leq 300$ 수준이며, 이정도 행렬의 크기는 기껏해야 수십~수백 KB 정도입니다. 대부분의 문제가 메모리를 128MB, 크게는 1024MB까지 준다는 걸 생각하면 정말 충분하다는 걸 알 수 있습니다.
Appendix. "result" Matrix의 모든 원소를 0으로 초기화하는 방법 (3) Allocation
한편 result를 배열로 선언하는 것은 stack size 이슈가 생길 수 있는데, 보통 BOJ에 제출하는 답안에서보다는 로컬에서 테스트할 때 문제가 될 수 있습니다. 만약 본인이 사용하는 행렬이 로컬에서 테스트하기도 힘든 크기라면 heap 영역에 allocate하는 것을 생각할 수 있습니다.
배열은 contiguous한 메모리 영역을 차지합니다. Loop를 돌면서 dynamic한 2d array를 할당받는 것은 배열의 가장 큰 장점 중 하나인 contiguous함을 버리는 것이며, cache hit을 크게 낮추는 것입니다. 또한, dynamic 1d array를 사용함에 수반되는 recalculating index 과정은 사실 cache hit보다 더 중요한 문제는 아닙니다. 그러므로 굳이 allocation을 쓴다면 1d가 낫습니다.
이제 배열을 0으로 초기화하는 걸 고려한다면, <stdlib.h>의 malloc보다는 calloc이 더 좋아 보입니다. 물론 malloc을 한 뒤 memset을 해도 빠르긴 하겠지만, calloc이라는 동일한 기능의 표준함수가 있으니 안 쓸 이유가 없습니다.
Appendix. Outside the PS Context
비슷한 내용을 다룬 국내 블로그 글 이 있어 첨언해보겠습니다. 저 글에서는 단순히 locality만 고려한 것보다 Strassen 등의 더 time-complexity가 좋은 알고리즘을 쓰면 더 나은 결과가 나올 수도 있다고 합니다. 그런데 어지간한 strassen 구현으로는 constant factor의 영향 을 무시할 수 없습니다. 만약 그런 알고리즘으로 행렬곱 코드 실행시간의 개선을 꾀한다면, 굳이 손으로 쓸 게 아니라 BLAS, LAPACK 같은 라이브러리를 가져오는 게 빠릅니다. 다만, 충분히 잘 구현한 템플릿을 들고 있다면, 1000 by 1000 이상에서는 써봄직도 하겠습니다.
Cache에 다 못 올릴 정도로 큰 행렬이라면 cache size에 맞게 반복문을 늘려주는 방법도 있습니다. (Samsung Software Membership Blog) 거기에 더해, Multi-Threading과 SIMD를 매우 적극적으로 활용해서 최적화한 사례 도 있습니다. 그리고 행렬곱을 가장 많이 쓰는 업계 중 하나인 그래픽스나, 딥러닝 같은 곳에서는 CUDA, OpenMp 같은 라이브러리로 GPU의 도움을 받곤 하는데, 이런것도 찾아보면 자료가 꽤 있습니다. ( CUDA , OpenMP )
PS에서 쓰기 어려울 뿐이지, 이런 최적화도 둘러보면 재밌는 것이 많은 것 같습니다.