🦾
行列累乗まとめ (競プロ)

2022/12/02に公開
 これはなに行列累乗と呼ばれる競プロのテクニックの概要と，それを用いる例をひたすら挙げていきます．競プロチックな話題ではありますが，行列が役立つ場面の例の紹介として，非競プロerでも楽しめると思います．ただし，数列と行列に関する用語がある程度わかっている人向けです．
なお，この記事は長野高専 Advent Calender 2022の2日目の記事です。

 原理(フィボナッチ数列の例)次の式で表される数列を考えます．

\begin{align}
F_1 &= 1 \\
F_2 &= 1 \\
F_n &= F_{n-1} + F_{n-2} \quad　(n \geq 3)
\end{align}
これはフィボナッチ数列という名で知られている，F=(1,1,2,3,5,8,13,21,34,\cdots)という前の2項の和が次の項になる数列です．これの第n項目をプログラミングで求めることにします．第1項目から第n項目までを求めるのではなく，第n項目のみ求めれば良いことに注意です．それっぽいコードを次に示します．
//F[i]:=フィボナッチ数列のi項目
F[1] = 1;
F[2] = 1;
for(int i=3;i<=n;i++){
	F[i] = F[i-1] + F[i-2];
}
for文をn回くらい回すので計算量はO(n)です(計算量についてこちら[1]を参照してください)．
では，n=10^{15}項目を求めたい場合はどうでしょうか．このコードでは，実行が終わりません．

 行列で表現するここで，行列です．フィボナッチ数列の定義より，次の式は成り立ちます．

\begin{pmatrix}
F_{n} \\
F_{n-1}
\end{pmatrix}
=
\begin{pmatrix}
1 & 1 \\
1 & 0
\end{pmatrix}
\begin{pmatrix}
F_{n-1} \\
F_{n-2}
\end{pmatrix}
具体的に計算してみます．

\begin{pmatrix}
2 \\
1
\end{pmatrix}
=
\begin{pmatrix}
1 & 1 \\
1 & 0
\end{pmatrix}
\begin{pmatrix}
1 \\
1
\end{pmatrix}

\begin{pmatrix}
3 \\
2
\end{pmatrix}
=
\begin{pmatrix}
1 & 1 \\
1 & 0
\end{pmatrix}
\begin{pmatrix}
2 \\
1
\end{pmatrix}

\begin{pmatrix}
5 \\
3
\end{pmatrix}
=
\begin{pmatrix}
1 & 1 \\
1 & 0
\end{pmatrix}
\begin{pmatrix}
3 \\
2
\end{pmatrix}
前の2項から次の項が生成されていることがわかりますね．

これらの結果を使って次のような変形を考えます．

\begin{align}
\begin{pmatrix}
5 \\
3
\end{pmatrix}
&=
\begin{pmatrix}
1 & 1 \\
1 & 0
\end{pmatrix}
\begin{pmatrix}
3 \\
2
\end{pmatrix}
\\
&=
\begin{pmatrix}
1 & 1 \\
1 & 0
\end{pmatrix}
\begin{pmatrix}
1 & 1 \\
1 & 0
\end{pmatrix}
\begin{pmatrix}
2 \\
1
\end{pmatrix}
\\
&=
\begin{pmatrix}
1 & 1 \\
1 & 0
\end{pmatrix}
\begin{pmatrix}
1 & 1 \\
1 & 0
\end{pmatrix}
\begin{pmatrix}
1 & 1 \\
1 & 0
\end{pmatrix}
\begin{pmatrix}
1 \\
1
\end{pmatrix}
\\
&=
\begin{pmatrix}
1 & 1 \\
1 & 0
\end{pmatrix}
^3
\begin{pmatrix}
1 \\
1
\end{pmatrix}
\end{align}
始めの2項に同じ行列を3回掛け算することで3つ後の項が出てきました．これを繰り返すと次が成り立つことがわかります．

\begin{align}
\begin{pmatrix}
F_n \\
F_{n-1}
\end{pmatrix}
&=
\begin{pmatrix}
1 & 1 \\
1 & 0
\end{pmatrix}
^{n-2}
\begin{pmatrix}
F_2 \\
F_1
\end{pmatrix}
\end{align}
第10項目を求めたいときは次のように計算できます．

\begin{align}
\begin{pmatrix}
F_{10} \\
F_{9}
\end{pmatrix}
&=
\begin{pmatrix}
1 & 1 \\
1 & 0
\end{pmatrix}
^{8}
\begin{pmatrix}
F_2 \\
F_1
\end{pmatrix}
\\
&=
\begin{pmatrix}
89 & 55 \\
55 & 34
\end{pmatrix}
\begin{pmatrix}
1 \\
1
\end{pmatrix}
\\
&=
\begin{pmatrix}
144 \\
89
\end{pmatrix}
\end{align}
フィボナッチ数列は1,1,2,3,5,8,13,21,34,55,89,144,\cdotsであるため，確かにあっています．

便利そうな式が行列によって完成しました．漸化式の悪いところは前の項の値がわかっていないと次の項が計算できないことですが，この式は特定の項をダイレクトに表すことに成功しています．行列のn乗が高速に計算できればフィボナッチ数列の第n項が高速に計算できそうです．

 累乗の高速計算行列Aのn乗であるA^nを計算するのに，ナイーブな方法だとn-1回の行列同士の掛け算が発生します．この行列同士の掛け算の回数が少なくなることを高速化と呼ぶことにします．

A^{128}を考えてみましょう．

AにAをかけてA^2

A^2にAをかけてA^3

\cdots

A^{127}にAをかけてA^{128}
全部で127回の掛け算が必要ですね．では，128という特徴的な数に注目して次のように計算したらどうでしょうか．

AにAをかけてA^2

A^2にA^2をかけてA^4

A^4にA^4をかけてA^8

A^{8}にA^{8}をかけてA^{16}

A^{16}にA^{16}をかけてA^{32}

A^{32}にA^{32}をかけてA^{64}

A^{64}にA^{64}をかけてA^{128}
全部で，7回の掛け算でA^{128}が計算できました．127回から7回という驚異の回数削減です．

実は，これは128のような2の冪数に限った話ではありません．先程，A^1,A^2,A^4,A^8,A^{16},\cdotsを求めましたが，これらを掛け合わせることで任意のA^nが計算できます(指数法則A^aA^b=A^{a+b}が成り立つことに注意)．
A^{12} = A^{8} A^{4}
A^{39} = A^{32} A^{4} A^{2} A^{1}
A^{127} = A^{64}A^{32}A^{16}A^{8}A^{4}A^{2}A^{1}
これらは，Aの(2進数表記したとき，1である位)乗を掛け算したものになっています．

上の例では，12_{(10)}=1100_{(2)}であり，1である位は4の位と8の位なので12=8+4と12を2の冪数の和で表すことができます．

このような手法をダブリングや繰り返し二乗法と呼びます．
それっぽいコードを次に示します．計算量はO(\log N)です．O(N)に比べて非常に高速です[2]．
// AのB乗の計算結果をCの格納
C = I //Iは単位行列
tmp = A
while(B>0){
	if(B%2==1){
		C = C * tmp
	}
	B/=2;
	tmp = tmp * tmp
}

 線形漸化式フィボナッチ数列は2項間の漸化式でしたが，k項間の線形漸化式でも可能です．

 3項の場合次で定義される数列aのn項目を求めよ．

\begin{align}
a_1 &= a_2 = a_3 = 1 \\
a_n &= 2a_{n-1} - 4a_{n-2} + 3a_{n-3}\quad　(n \geq 4)
\end{align}
この漸化式は次のような行列で表されます．

\begin{align}
\begin{pmatrix}
a_n \\
a_{n-1} \\
a_{n-2}
\end{pmatrix}
&=
\begin{pmatrix}
2 & -4 & 3\\
1 & 0 & 0 \\
0 & 1 & 0 
\end{pmatrix}
^{n-3}
\begin{pmatrix}
a_3 \\
a_2 \\
a_1 
\end{pmatrix}
\end{align}
実際に行列を計算すると確かに成立していることがわかります．
一般に，k+1項間の漸化式はk\times kの行列を使って表され，k\times kの行列同士の掛け算はk^3回の掛け算が必要であるため，第n項を求めるのにO(k^3\log n)かかります．
なお，k+1項間の線形漸化式の第n項はkitamasa法と呼ばれる手法で(k^2\log n)で求めることが可能です．こちらもダブリングを使った手法です．

 定数項がある場合a,b,m,x_0を使って線形合同法で乱数を生成する．このときのn番目の乱数x_nを求めよ

ただし，線形合同法とは，

x_{i+1} = ax_i + b \pmod m
の漸化式で疑似乱数を生成する手法である．
定数項がある場合は1の行を付け足して次のように作ることができます．

\begin{align}
\begin{pmatrix}
x_n \\
1 \\
\end{pmatrix}
&=
\begin{pmatrix}
a & b\\
0 & 1 
\end{pmatrix}
^{n}
\begin{pmatrix}
x_0 \\
1
\end{pmatrix}
\end{align}

 複数変数の線形漸化式複数の変数で表される漸化式にも応用できます

 平方根を含む数の累乗次の問題を考えます．
(2+\sqrt{3})^nはいくつか？a+b\sqrt{3}の形式になるのでa,bを求めよ．
nが大きいと展開が厄介になりそうなことが想像できますね．

とりあえず，2+\sqrt{3}の累乗を次のように文字で置きます．

X_n = (a_n + b_n\sqrt{3})
X_nからX_{n+1}を計算してみます．

\begin{align}
X_{n+1} &= X_n (2 + \sqrt{3}) \\
	&= (a_n + b_n \sqrt{3})(2 + \sqrt{3}) \\
	&= 2a_n + a_n \sqrt{3} + 2b_n \sqrt{3} + 3b_n \\
	&= (2a_n + 3b_n) + (a_n + 2b_n)\sqrt{3} \\
	&= a_{n+1} + b_{n+1} \sqrt{3}
\end{align}
最後の二行の係数を比較して，次の漸化式が立ちます．


\begin{align}
a_1&=b_1=1 \\
a_{n+1} &= 2a_{n} + 3b_n \\
b_{n+1} &= a_n + 2b_n
\end{align}

これを行列にすると次のようになります．

\begin{align}
\begin{pmatrix}
a_n \\
b_n 
\end{pmatrix}
&=
\begin{pmatrix}
2 & 3\\
1 & 2
\end{pmatrix}
^{n-1}
\begin{pmatrix}
a_1 \\
b_1
\end{pmatrix}
\end{align}

 確率次の問題を考えます．
エアコンのスイッチの状態はOFFである．遠隔操作を1回行うと確率pでエアコンのONとOFFが入れ替わる．n回遠隔操作をしたときエアコンがONである確率を求めよ．(N\leq 10^{18})　出典 CODE FESTIVAL 2014 Middle C - eject
次のように変数を設定します．

a_i = i回目の遠隔操作をしたときにエアコンがONである確率

b_i = i回目の遠隔操作をしたときにエアコンがOFFである確率

漸化式は次のように立ちます．

a_i = p b_{i-1} + (1-p)a_{i-1} \\
b_i = p a_{i-1} + (1-p)b_{i-1} \\
行列にします．ただし，a_0=0,b_0=1です．

\begin{align}
\begin{pmatrix}
a_n \\
b_n 
\end{pmatrix}
&=
\begin{pmatrix}
1-p & p\\
p & 1-p
\end{pmatrix}
^{n}
\begin{pmatrix}
a_0 \\
b_0
\end{pmatrix}
\end{align}
今までのものは全て，小さな問題の答えから大きな問題の答えを導くdp(dynamic programming:動的計画法)に分類されます．dpの漸化式が立ったときに，それが線形であり，同じ遷移をする場合に行列累乗が有効なことが多いです．

 グラフのパス数EDPC R - Walk
N頂点の単調有効グラフGで，長さ1の有向辺が長さKのパスは何通りありますか．ただし同じ頂点を複数回通っても良い
N次の正方行列で，頂点iから頂点jまでの有向辺があればA_{ij}=1，なければA_{ij}=0となるような隣接行列を考えます．
例えば，次の4頂点のグラフを隣接行列にすると次のようになります．

A = 
\begin{align}
\begin{pmatrix}
0 &1 &0 &1 \\
0 &0 &1 &0 \\
0 &1 &0 &0 \\
0 &0 &0 &1 
\end{pmatrix}
\end{align}
このとき，頂点iから頂点jへのパス数は次のように計算できます．

(頂点iから頂点jへのパス数) = \sum _{k=1}^4 A_{ik} A_{kj}
これは，行列の積の定義とよく似ていて，A^2の(i,j)成分は，頂点iから頂点jへの長さ2のパス数を表します．よって，A^nのi,j成分から頂点iからjまでの長さnのパス数を求めることができます．
なお，無向グラフの場合は双方向に有向辺が張っていると考え，A_{ij}=A_{ji}とすることで同様に処理できます．

 行列の級数正方行列Aがある．A+A^2+A^3+A^4+A^5+\cdots+A^Nを求めよ．
次のようにX_nを定義します．

X_n = A^1 + A^2 + A^3 + \cdots + A^n
漸化式は次のように作れます．

\begin{align}
X_1 &= A\\
X_n &= A+AX_{n-1} \quad (n\geq 2)
\end{align}
よって，

\begin{align}
\begin{pmatrix}
X_n \\ \hline
I 
\end{pmatrix}
&=
\begin{pmatrix}
\begin{array}{c|c} 
A & A\\ \hline
O & I
\end{array}
\end{pmatrix}
^{n-1}
\begin{pmatrix}
A \\ \hline
I
\end{pmatrix}
\end{align}
です．

 半環普通，行列の掛け算は次のように定義されています．

A = 
\begin{align}
\begin{pmatrix}
a_{11} & a_{12} \\
a_{21} & a_{22}
\end{pmatrix}
\begin{pmatrix}
b_{11} & b_{12} \\
b_{21} & b_{22}
\end{pmatrix}
=
\begin{pmatrix}
a_{11}\cdot   b_{11}+a_{12}\cdot   b_{21} & a_{11}\cdot   b_{12}+a_{12}\cdot   b_{22} \\
a_{21}\cdot   b_{11}+a_{22}\cdot   b_{21} & a_{21}\cdot   b_{12}+a_{22}\cdot   b_{22}
\end{pmatrix}
\end{align}
それぞれ要素の積を計算し，和を計算しています．ここで，「和」と「積」を別のものに置き換えることを考えます．例えば次のように\maxと+に変えるとどうなるでしょうか．ここで，\max(a_1,a_2,a_3,\cdots)は，a_1,a_2,a_3,\cdotsのうちの最大値を表します．

A = 
\begin{align}
\begin{pmatrix}
a_{11} & a_{12} \\
a_{21} & a_{22}
\end{pmatrix}
\begin{pmatrix}
b_{11} & b_{12} \\
b_{21} & b_{22}
\end{pmatrix}
=
\begin{pmatrix}
\max(a_{11}+b_{11},a_{12}+b_{21}) & \max(a_{11}+b_{12},a_{12}+b_{22}) \\
\max(a_{21}+b_{11},a_{22}+b_{21}) & \max(a_{21}+b_{12},a_{22}+b_{22})
\end{pmatrix}
\end{align}
この\max,+で計算される行列ですが，実はこの行列においても高速な行列累乗の計算が可能です．こんな変な演算の世界で行列累乗をしてどうするんだと思いますが，例は後述します．もちろん，演算子を好き勝手変えていいわけではありません．いくつか条件があります．次に示す条件を満たす，集合と加法+と乗法 \cdot の2つの二項演算，である必要があります．この性質を満たす集合と+と\cdotの組を半環と言います．また，集合の要素を元と言います．
加法において，結合法則が成り立ち，可換であり，単位元0を持つ．
(a + b) + c = a + (b + c)
a+b = b+a
0 + a = a + 0 = 0

乗法において，結合法則が成り立ち，単位元1を持つ．
(a \cdot   b) \cdot   c = a \cdot   (b \cdot   c)
1 \cdot   a = a \cdot   1 = a

分配法則が成り立つ．
a \cdot   (b+c) = (a\cdot   b) + (a\cdot   c)
(a+b)\cdot   c = (a\cdot   c) + (b\cdot   c)

乗法で，0(加法の単位元)と任意の元の積は0になる．
0\cdot   a = a\cdot   0 = 0

例えば，有理数の集合を\mathbb{R}とすると，(\mathbb{R},+,\cdot)は半環をなす，という言い方をします．
では，先程の\maxと+の例は本当に半環をなしているかを確認してみましょう．この場合の集合も\mathbb{R}(有理数全体)としておきます．先述した半環の定義だと，加法が\maxで，乗法が+に対応しています．

まずは加法\maxについて見ていきます．
\max(\max(a,b),c) = \max(a,\max(b,c))
\max(a,b) = \max(b,a)
\max(\infty,a) = \max(a,\infty) = \infty
\maxは単位元を\inftyにすることで成立します．次に乗法+です．
(a + b) + c = a + ( b+c)
0 + a = a + 0 = a
+の単位元は0です．

次に分配法則の確認です．ややこしいですが，次の式が成り立つことがわかります．
a+\max(b,c) = \max(a+b,a+c)
\max(a,b)+c = \max(a+c,b+c)
最後に加法\maxの単位元\inftyを乗算+すると\inftyになることを確認します．
\infty +a = a + \infty = \infty
以上，厳密性には欠けるかもしれませんが(\mathbb{R},\max,+)が半環をなしていることがわかりました．さて，次の節からは，いつもと違う半環の世界での行列累乗の例を見ていきます．

 ANDとXORABC009 D - 漸化式
定数C_1,C_2,C_3,\cdots,C_Kがあり，K項間の漸化式が次のように定まる

A_i \quad (1\leq i\leq K)は与えられる．

A_{K+1} = (C_1 \mathrm{AND} A_{K}) \mathrm{XOR} (C_2 \mathrm{AND} A_{K-1}) \mathrm{XOR}  \cdots \mathrm{XOR} (C_K \mathrm{AND} A_{1})
このとき，A_Nを求めよ(1\leq M\leq 10^{9})．

なお，\mathrm{AND}はビットごとの論理積，\mathrm{XOR}はビットごとの排他的論理和を表す．
この問題は(\mathbb{N},\mathrm{AND},\mathrm{XOR})が半環をなすことを利用して，ただの線形漸化式と見ることができます．なお，\mathrm{AND}と\mathrm{XOR}はビットごとで独立した演算であるため，各bitで\mod 2の掛け算とたし算をしていると見ることでも解くことができます．

 ワーシャルフロイド法ワーシャルフロイド法は，(\mathbb{R},\min,+)の世界での行列累乗だと捉えることができます．
MojaCoderで類題を見つけました．

Dungeon Attack (Hard)

 関連問題行列累乗がストレートに出ることは珍しく，考察の末に最後の一捻りとして出てくることが多い印象があります．全体的に難しいです．
フィボナッチ数列の第N項をMで割った余りを求める
ABC293 E - Geometric Progression
ABC357 D - 88888888
ABC204 F - Hanjo 2
ABC199 F - Graph Smoothing
ABC256 G - Blakc and Whte Stones
ABC271 G - Access Counter 

 参考文献動的計画法を実現する代数〜トロピカル演算でグラフの最短経路を計算する〜

脚注
APG4b 計算量 https://atcoder.jp/contests/APG4b/tasks/APG4b_w ↩︎
実はこの計算量解析はあまり意味のあるものではありません．なぜなら，フィボナッチ数列は後の項に行けば行くほど桁数が増えるからです．後の項に行くほど，かけ算やたし算の計算コストは上がっていきます．ですので，O(\log N)といってもN=10^7の時点で結構時間がかかってしまいます．競プロでは10^9+7では割ったあまりを求めなさい，のように桁数を気にせずに計算できる場合が多いため，今回は考えないことにしています． ↩︎
これはなに

原理(フィボナッチ数列の例)

行列で表現する

累乗の高速計算

線形漸化式

3項の場合

定数項がある場合

複数変数の線形漸化式

平方根を含む数の累乗

確率

グラフのパス数

行列の級数

半環

ANDとXOR

ワーシャルフロイド法

関連問題

参考文献

Discussion
