はじめに

この記事では、競技プログラミングにおいて使える Black box linear algebra, 特に Wiedemann [1] によるアプローチのテクニックについていくつか書きます。

この記事はCompetitive Programming Advent Calendar 2015の4日目の記事です。

導入

競技プログラミングにおいて、数え上げ問題は1つの大きなジャンルです。 多くの場合、$1000000007$などの数で割った余りを求めることを行います。今回は、特に素数による剰余を取る場合について考えていきます。

数え上げ問題の解法には多くの種類がありますが、その中でも大きなものとしてあるのは行列を用いた解法です。線形代数を用いているとも言えるでしょう。 そのような例として、今回は以下の2つを取り上げます。

行列累乗によるDPの高速化

行列累乗は数え上げ問題で非常によく用いられるテクニックです。 状態遷移を行列として表することでDPを表現します。 つまり、初期状態をベクトル$b$, 遷移行列を$A$としたとき、$K$回遷移した後の状態は$A^K b$として表現できます。 バイナリ法により累乗することにより、行列サイズを$N$としたとき$\Theta(N^3 \log K)$時間で答えを求められます。

行列累乗は様々な特殊ケースにおいてさらに高速化ができることがあります。 たとえば行列が巡回行列である場合、そのような行列同士の掛け算は$O(N^2)$時間(または高速多項式乗算を用いてその計算量)で求めることができ、 また巡回行列同士の積は巡回行列であるので合計$O(N^2 \log K)$時間で実行できます。 三角テプリッツ行列の場合でも同じように高速化できます。

今回は、高速に行列累乗ができるクラスを上記の特殊ケースたちから大きく広げ、一般の場合の行列累乗をも高速化できるテクニックを紹介します。

行列式を求める問題

競技プログラミングにおいて、行列式が必要になることはよくあります。 これは準同型写像であり、非常に基本的なものであるといえ、応用も数えきれないほどあります。 そのため、それを使う問題が多くあることは驚くことではないでしょう。

たとえば、行列木定理と呼ばれる定理は、 グラフが与えられたとき、全域木の個数を行列式を用いて求められると述べています。 非自明な数え上げを多項式時間で行ってしまうというのだから行列式の力の大きさがわかるでしょう。

$N \times N$の行列の行列式はガウスの消去法を用いて$\Theta(N^3)$時間で求めることができます。 より高速にできる特殊ケースとしては、$K$重対角行列の場合$\Theta(N K^2)$時間でガウスの消去法を実行できます。

今回は、他の大きなクラスにおいて計算を高速化させます。

準備

線型漸化列とその最小多項式

$V$を体$F$上のベクトル空間とする。 $V$上の無限列$\{a_i\}_{i=0}^{\infty}$が線型漸化的であるとは、 ある$F$上の多項式$\displaystyle c(X) = \sum_{i=0}^n c_i X^i$が存在し、全ての$0 \le i$に対して$\displaystyle \sum_{j=0}^n c_j a_{i+j} = 0$となることをいう。 このとき、$c$が$a$をannihilateしていると呼ぶ。

$a$をannihilateする多項式の集合はイデアルを成す。 多項式環はPIDなので生成元が存在するが、零イデアルでないならそのうちleading係数が1のものが唯一に決まり、それを$a$の最小多項式と呼ぶ。 ただし、零イデアルが生成される場合$0$を最小多項式とする。そうでない場合$a$を線型漸化的であると呼ぶ。

$f$を2つのベクトル空間上の線形写像とすると、$f(a)$の最小多項式は$a$の最小多項式の約数である。

行列の最小多項式と特性多項式

正方行列$A$に対し、列$\{A^i\}_{i=0}^{\infty}$を考える。 この列の最小多項式を行列$A$の最小多項式と呼ぶ。

正方行列$A$の特性多項式は${\rm char}(A) = {\rm det}(X I - A)$として定義される。ここで$X$は多項式の変数。 $n$を$A$の大きさとするとき、${\rm char}(A)(0) = (-1)^n {\rm det}(A)$であることを確認しておく。

ケイリー・ハミルトンの定理は、 $\{A^i\}$は線型漸化的であり、最小多項式は特性多項式の約数であると述べている。 また、大きさ$n$の行列の最小多項式の次数は$n$以下であることは帰結である。

行列$A$のある固有値とそれに対応する固有ベクトルを$\lambda, v$としたとき、 多項式$f$に対し$f(A) v = f(\lambda) v$であることが確認できる。つまり、$\{A^i\}_i$をannihilateする$f$に対し、全ての固有値は$f$の根となる。 また、最小多項式が特性多項式の約数であることと組み合わせると、最小多項式と特性多項式の根の集合は一致するということが言える。

最小多項式を求める

最小多項式はどうやっても求めればいいのだろうか？ そもそも、無限列は普通には入力できない。 そこで、最小多項式の次数の上限$n$と列の先頭$2n$要素を入力することにする。 先頭$2n$要素から最小多項式が一意に求められることは、以下のアルゴリズムと一緒に構成的に証明される。

体上の場合

入力列がある体$F$の要素である場合を考える。 この場合、拡張ユークリッド互除法やBerlekamp–Massey algorithmを用いて $O(n^2)$時間で求めることができる。詳細は文献を参考のこと。 なお、後者のアルゴリズムの実装は非常に簡単である。

ベクトルの場合

次に、入力列$b$がある体$F$のベクトル$F^m$の要素である場合を考える。 まず、$F$の有限部分集合$U$を取り、その$U$から一様ランダムに$m$要素を選びベクトル$u$とする。 そして$a_i = u^T b_i$として新たな$F$上の列$a$を計算する。 このとき$a$の最小多項式が高確率で$b$の最小多項式と一致するというのが重要な定理である。詳細は割愛する。 合計$O(mn + n^2)$時間となる。

行列の場合

最後に、行列$A$の最小多項式を求めることを考える。 この場合、上の場合と同じように、ランダムなベクトル$b$を取って列$\{A^i b\}_{i=0}^{2n-1}$を求め、 その列に対して上記のアルゴリズムを適用することにより、これもまた高確率で$A$の最小多項式を求められる。 各$i$に対し$A^i b$を求めるのに行列乗算は必要なく、あるベクトルと$A$の乗算のみでよいことに注意する。 ここで、ベクトルと$A$の乗算に$T(n)$時間かかるとすると、合計$O(n^2 + n T(n))$時間で最小多項式が求められた。

最小多項式の応用

行列累乗の高速化

行列$A$, ベクトル$b$, 自然数$K$に対し、$A^K b$が求めたいのであった。 ここで、上記で述べたとおり列$\{A^i b\}_i$が線型漸化的であることに注目する。

まず、列$\{A^i b\}_i$の最小多項式を求める。すると、問題は線型漸化式の第$K$項を求めるものとなる。 ここで、多項式を使うテクニックたちに書いてある方法が使える。 その方法を用い$A^K b$を$A^i b (i < n)$の線型和で表せば、単純に求めた初項と掛けあわせればよい。 ベクトルと$A$の乗算を$T(n)$時間で求められるとすると、合計$O(n^2 + n T(n) + M(n) \log K)$時間で実行できる。 ここで、$M(n)$は多項式乗算の時間計算量を表す。

上記の計算量は、一般の場合$T(n) = \Theta(n^2)$の場合でも$O(n^3 + M(n) \log K)$となり、高速化となっている。 また、巡回行列や(三角とは限らない)テプリッツ行列の場合は$O(n M(n) + M(n) \log K)$となる。 さらに、スパースな行列で非0要素が$S$個のものの場合、ベクトルにそれを乗算することは$O(S)$時間で可能なので、$O(n^2 + n S + M(n) \log K)$で計算できることとなる。 DPの高速化の場合、遷移行列がスパースである場合は多い。

行列式を高速に求める

サイズ$n$の正方行列$A$が与えられたとき、行列式${\rm det}(A)$を求めたいのであった。 行列式は$(-1)^n {\rm char}(A)(0)$であることを思い出すと、$A$の特性多項式を求められればいいこととなる。

もし$A$の最小多項式の次数が$n$であれば、それは特性多項式と一致する。 また、上記のとおりそれらは根の集合が一致する。 しかし最小多項式の次数は$n$より小さくなることがある。 アイデアは、$A$に前処理を行うことで最小多項式と特性多項式を一致させるというものである。

まず、${\rm det}(A) = 0$である場合、特性多項式と最小多項式は共に$0$を根に持つということになるため、それにより判定できる。 $A$が非特異行列であると仮定する。 $D$をランダムに選んだ${\rm det}(D) \neq 0$となる対角行列とする。このとき、$AD$の最小多項式は高確率で次数$n$となるというのが重要な補題である。割愛する。 これにより$AD$の最小多項式を求め、それから${\rm det}(AD)$を求め、${\rm det}(D)$で割れば答えが出る。

ベクトルと$A$の乗算を$T(n)$時間で求められるとすると、ベクトルと$AD$の乗算は$O(n + T(n))$時間で求められるので、合計$O(n^2 + n T(n))$時間で行列式が求められる。

終わりに

Black box linear algebra

上記の例では、実際には行列は明示的に与えられなくてもよいことに注意する。 つまり、ベクトルを入力としてある行列を掛ける(別の言葉で言うと、線形変換を行う) 「ブラックボックス」さえ入力されれば、その中身に関係なくアルゴリズムは実行できるということである。 このような線形代数のアルゴリズムを"Black box linear algebra"と呼ぶ。

上記で紹介した以外にも様々なアルゴリズムが発見されている。 例えば、線形方程式$Ax = b$も、同じようなアプローチで高速に解くことができる。

結論

今回は非常に幅広い計算問題を高速化する方法を紹介した。 競技プログラミングにおいても線形代数を使う問題は豊富にある。 たとえこの方法が想定解でなかったとしても使える場面は多いだろう。 是非活用していただきたい。

参考文献