Samuelson-Berkowitz Algorithm

https://www.acmicpc.net/problem/25510

https://infossm.github.io/blog/2022/08/18/CharPoly/

Division-Free Determinant에 대해 위 문제를 내고 글도 썼는데, 최근 훨씬 간단한 알고리즘이 있다고 해서 정리해보려고 한다.

원본 글에 달린 것과 동일하게 $\mathcal{O}(n^4)$에 동작하는 division-free characteristic polynomial 알고리즘이다.

Algorithm

$n \times n$ 행렬 $A = (a_{ij})$에 대해, $A_{i}$를 $i\dots n$번째 행, 열만 따와서 만든 크기 $n-i+1$의 정사각행렬이라고 하자. 당연히 $A_{1} = A$이다.

$p_{1}(x) = \det(x I - A)$를 cofactor expansion으로 계산하면 $p_{1}(x) = (x - a_{11})\det(x I - A_{2}) - a_{12} \det (\cdots) \cdots$ 와 같은 형태로 나올 텐데, 사실 훨씬 간단한 형태로 계산할 수 있다.

Lemma.
$A$를 $\begin{pmatrix} a_{11} & R \\ C & A_{2} \end{pmatrix}$로 쓰자. 즉 $R, C$를 각각 크기 $n-1$의 row/column 벡터로 볼 때,
$$p_{1}(x) = (x - a_{11})p_{2}(x) - R\mathrm{adj}(x I - A_{2})C.$$

여기서 $\mathrm{adj}(A)$는 $A \mathrm{adj}(A) = (\det A)I 를 만족하는 행렬로,$(-1)^{i+j}\det (A_{-j, -i})$를$(i, j)$원소로 가진다.$A_{-i, -j}$는$A$의$i$번째 행,$j$번째 열을 제외한 크기$n-1$의 정사각행렬인데, 말이 어렵지만 사실 cofactor expansion을 행렬로 적은 것이다.

Proof. $(xI - A)\mathrm{adj}(xI - A) = p_{1}(x) I$ 의 $(1, 1)$을 보자. $xI - A$의 첫번째 row는 $\begin{pmatrix} x - a_{11} & -R$이고, $\mathrm{adj}(xI - A)$의 첫번째 column은 당연히 첫 원소에 $p_{2}(x)$가 있을 테지만, $j > 1$에 대해서는 $xI - A_{2}$의 $j-1$번째 column을 빼고, 첫 행에 $-C$를 집어넣은 행렬에 $(-1)^{1+j}$를 곱한 값이다.

여기서 다시 cofactor expansion을 진행하자. 이번에는 첫 column ($-C$)를 빼고 다른 row들을 빼면서 내려가면, 이는 $(-1)^{1+j} \sum_{i=1}^{n-1} (-1)^{1+i} (-C) _ {i} \det (A_{2}) _ {-i, -(j-1)}$ 이 된다. 이는 곧 $(\mathrm{adj}(xI - A_{2}) C) _ {j-1}$이 됨을 알 수 있다. $\square$

따라서 $R \mathrm{adj}(x I - A_{2}) C$ 를 계산하는 게 중요한데, $\mathrm{adj}(x I - A)$ 꼴에 대한 다음 lemma가 중요하다.

Lemma.

$\det(xI - A) = p_{0} x^{n} + \cdots + p_{n}$이라고 하자. $\mathrm{adj}(xI - A) = \sum_{k=1}^{n} x^{n-k} (p_{0}A^{k-1} + \cdots + p_{k-1})$가 성립한다.

Proof. Cayley-Hamilton theorem에 의해서 $p(A) = 0$이다. $p(x) - p(A) = \sum_{k = 0}^{n} p_{k} (x^{n-k}I - A^{n-k})$를 생각하고, 이들을 $xI - A$로 나누면 $\sum_{k = 0}^{n} p_{k} \sum_{i = 0}^{n-k-1} x^{n-k-1-i}A^{i}$가 되고, 이걸 정리하면 위 식을 얻는다.

따라서 이 두 식을 조합하면, $p_{2}(x) = q_{0}x^{n-1} + \cdots + q_{n-1}$이라고 하면
$$p_{1}(x) = (x - a_{11})p_{2}(x) - \sum_{k=0}^{n-2} q_{k}(x^{n-k-2} + (RA_{2}C)x^{n-k-3} + \cdots + (RA_{2}^{n-k-2}C))$$

를 얻는다. 이 때 $R A_{2}^{i} C$는 $i = 1, \cdots, n$에 대해 $\mathcal{O}(n^3)$ 시간에 계산할 수 있고, $A_{2}$에 대해 재귀적으로 문제를 해결하면 전체 $\mathcal{O}(n^{4})$ 시간복잡도를 갖는다.

위키피디아의 식이 이렇게 유도되는데, 이를 잘 분해하면 $n$개의 toeplitz matrix (대각선을 따라 constant인) 의 곱을 얻는데, 이는 원래 행렬 $A$에서 모두 얻을 수 있어 병렬적으로 계산이 가능하다. 아쉽게도 두 toeplitz matrix의 곱이 toeplitz라는 보장은 없기 때문에 단순 곱셈보다 빠르게 계산할 수는 없는 것 같다.