Deterministic Verification of Integer Matrix Multiplication in Quadratic Time
세 개의 $n \times n$ 크기의 정수 행렬 $A, B, C$에 대해 $AB == C$인지 $O(N^2)$에 판별하는 알고리즘을 알아보자.
PS에서 보통 사용되는 $O(KN^2)$의 Freivalds' algorithm을 개선한 버전으로 생각할 수도 있지만, 이 알고리즘은 PS에선 쓰이지 않는다.
먼저, naive하게 $AB$를 계산하고 $C$와 비교하는 알고리즘의 시간복잡도는 $O(N^3)$이다. $AB$의 계산에서 $n \times n$ 크기의 두 행렬을 곱해야 하므로 $O(N^3)$의 계산이 필요하고, 계산된 결과가 $C$와 일치하는지 확인하는 데에 $O(N^2)$의 비교가 필요하기 때문이다. 조금 더 효율적인 방법이 없을까?
$AB=C$라면 임의의 $n \times 1$ 열벡터 $x$에 대해 $ABx = Cx$를 만족한다. $\cdots (*)$
$ABx$의 값은 $Bx$를 먼저 계산한 뒤 $A(Bx)$를 계산하는 방법으로 $O(N^2)$에 계산이 가능하고, 따라서 임의의 열벡터 $x$에 대해 $ABx = Cx$인지 확인하는 것 또한 $O(N^2)$에 가능하다.
어떤 $x$에 대해 $ABx=Cx$인지 비교하고 그 결과를 $AB=C$의 비교결과 대신 사용하는 알고리즘을 생각해보자.
임의의 열벡터 $x$를 골라 $ABx=Cx$인지 비교한 결과는 $AB=C$의 성립 여부를 완벽히 대변할 수 있을까?
$(*)$의 대우명제를 생각하면 $ABx \neq Cx$라면 $AB \neq C$임은 자명하지만
$ABx = Cx$라면 $AB = C$일 수도 $AB \neq C$일 수도 있다.
따라서 알고리즘이 틀릴 확률은 "$ABx = Cx$이고 $AB \neq C$일 확률"과 같다.
Freivalds' algorithm은 열벡터 $x$의 원소를 $\{0, 1\}$ 중에서 무작위로 선택할 때 이러한 확률이 $\left(\dfrac{1}{2}\right)$ 미만임을 이용한다.
해당 과정을 $K$번 반복해 알고리즘이 틀릴 확률을 $\left(\dfrac{1}{2}\right)^k$로 낮추는 것이다.
보통 $K=30$정도로만 해도 굉장히 높은 정확도를 보장하지만, 그래도 틀릴 확률이 존재한다는 점과 높은 정확도를 위해 시간복잡도를 포기해야 한다는 점은 아쉬울 수 밖에 없다.
그렇다면 애초에 "$ABx = Cx$이고 $AB \neq C$일 확률"이 $0$이 되도록 $x$를 정할 수는 없을까?
먼저 $x$를 다음과 같이 정의하자.
$$x=(1, r, r^2, \cdots, r^{n-1})^T, r \in \mathbb{R}$$
Lemma 1.
$n \times n$ 크기의 영행렬이 아닌 실수 행렬 $D$가 주어질 때 다음을 만족하는 실수 $r$의 개수는 최대 $n-1$개이다.
$$D (1, r, r^2, \cdots, r^{n-1})^T = 0$$증명은 간단한데, 곱해진 결과에서 $i$번째 행의 값을 $P_i(r)=\displaystyle\sum_{j=1}^{n} {d_{i,j} r^{j-1}}$라 하자. $\forall{i}, P_i(r) = 0$이다.
$D$는 영행렬이 아니므로 $0$이 아닌 행 $D_i$가 최소한 하나 존재하고, 그러한 행 $i$에서 $\exists{d_{i, j} \neq 0}$이므로 $P_i(r)=0$은 항등식이 아닌 최대 $(n-1)$차 방정식이니 실근 $r$의 개수 또한 최대 $n-1$개이다.
이제 알고리즘이 틀릴 확률인 "$ABx = Cx$이고 $AB \neq C$일 확률"에 대해 생각하자.
$Y=AB-C$라 하면 $Y \neq 0, Yx=0$이다.
(사실 $r$을 무작위한 정수로 잡아도 이럴 확률은 굉장히 낮다. 잘못된 판단이 일어나려면 무작위로 선택한 $r$이 모든 $i$에 대해 $P_i(r)=\displaystyle\sum_{j=1}^{n} {y_{i,j} r^{j-1}}=0$의 근이어야 하기 때문이다.)
알고리즘이 틀릴 확률을 없애기 위해서는 " $Y \neq 0, Yx=0$ "이 불가능한 $x$를 찾아 사용하면 된다. 이를 위해 다항식 $P_i(r)$이 $0$이 되지 않는 $r$을 찾자.
Cauchy bound
실계수 다항식 $P(x)=a_{k}x^{k} + a_{k−1}x^{k−1} + ...+ a_{1}x + a_{0}$에 대해 $P(x)=0$의 실근 $x$는 다음 조건을 만족한다.
$$|x| \lt 1+\dfrac{A}{|a_k|}, (A = \max_{i=0}^{k-1}|a_i|)$$
Cauchy bound는 $r$이 $P_i(r) = 0$의 근이 아니도록 하기 위한 간단한 하한을 제공한다.
$Y=AB-C$이므로 행렬의 원소 $y_{i,j}=\displaystyle\sum_{k=1}^{n}{a_{i,k}b_{k,j}} -c_{i,j}$이다.
$c_{max}=\max\{|a_{i,j}|, |b_{i,j}|, |c_{i,j}|\}$라고 하면, 다음과 같이 다항식 $P_i(r)=\displaystyle\sum_{j=1}^{n} {y_{i,j} r^{j-1}}$에서 절대값이 가장 큰 계수의 상한을 정할 수 있다.
$$A \le n c_{max}^2+c_{max}$$
또한 $Y$는 정수행렬이므로 다항식 $P_i(r)$에서 최고차항의 계수의 절대값은 최소 1이고, $P_i(r)=0$의 실근의 절대값은 $\alpha = n c_{max}^2+c_{max}+1$ 미만이 된다.
따라서 $r \ge \alpha$인 $r$을 선택하면 모든 $i$에 대해 $P_i(r) \neq 0$이므로 "$Y \neq 0$이며 $Yx = 0$인 경우"는 존재하지 않는다.
다르게 말하면, 이러한 $r$에 대해서 $AB=C$의 비교결과는 $ABx=Cx$의 비교결과와 항상 일치한다.
최종적으로 알고리즘을 정리하면 다음과 같다.
- 행렬 $A, B, C$의 원소들의 절대값 중 최대값 $c_{max}$를 찾는다. $O(N^2)$
- $\alpha = n c_{max}^2+c_{max}+1$에 대해 $r:=\alpha$로 설정한다. $O(1)$
- 벡터 $x=(1, r, r^2, \cdots, r^{n-1})^T$를 계산한다. $O(N)$ (모든 산술연산은 $O(1)$이라 가정)
- $A(Bx)=Cx$이면 $AB=C$이고, 아니면 $AB \neq C$이다. $O(N^2)$ (모든 산술연산은 $O(1)$이라 가정)
- freivalds' algorithm과는 달리 알고리즘이 틀릴 확률이 없기에 여러 번 반복할 필요도 없다. 한 번의 계산으로 완벽한 판별이 가능하다.
사실 PS에서 사용하기엔 계산에 사용되는 수의 범위가 너무나 크다.
시간복잡도 $O(N^2)$을 위해서는 $N$과 $c_{max}$의 값이 충분히 작아 모든 산술연산을 $O(1)$에 수행할 수 있어야 하는데, 보통 PS에서 주어지는 $n$의 범위를 생각하면 $r^{n-1}$이라는 수는 아득히 큰 수다.($n=1e5, c_{max}=1e5$라면 $r^{n-1}$은 약 백만자리의 수가 된다) FFT를 사용한다고 해도 $\log_{10}(r^{n-1})=O(N)$자리의 큰 수이므로 숫자끼리 한 번 곱하는 시간복잡도만 $O(N \log{N})$이 필요하며, $A(Bx)$의 계산과정에서 이러한 곱셈이 $O(N^2)$번 사용될테니 결국 총 시간복잡도는 $O(N^3 \log{N})$이 된다. 이는 직접 $AB$를 계산하고 $C$와 비교하는 단순한 $O(N^3)$ 알고리즘보다도 느리다.
References
Korec, I., Wiedermann, J. (2014). Deterministic Verification of Integer Matrix Multiplication in Quadratic Time. In: Geffert, V., Preneel, B., Rovan, B., Štuller, J., Tjoa, A.M. (eds) SOFSEM 2014: Theory and Practice of Computer Science. SOFSEM 2014. Lecture Notes in Computer Science, vol 8327. Springer, Cham. https://doi.org/10.1007/978-3-319-04298-5_33