本文介紹線性代數中一個非常重要的內容——矩陣（Matrix），主要講解矩陣的性質、運算以及在常系數齊次遞推式上的應用。
定義
對於矩陣 \(A\)，主對角線是指 \(A_{i,i}\) 的元素。
一般用 \(I\) 來表示單位矩陣，就是主對角線上為 1，其余位置為 0。
性質
矩陣的逆
\(A\) 的逆矩陣 \(P\) 是使得 \(A \times P = I\) 的矩陣。
逆矩陣可以用 高斯消元 的方式來求。
運算
矩陣的加減法是逐個元素進行的。
矩陣乘法
矩陣相乘只有在第一個矩陣的列數和第二個矩陣的行數相同時才有意義。
設 \(A\) 為 \(P \times M\) 的矩陣，\(B\) 為 \(M \times Q\) 的矩陣，設矩陣 \(C\) 為矩陣 \(A\) 與 \(B\) 的乘積，
其中矩陣 \(C\) 中的第 \(i\) 行第 \(j\) 列元素可以表示為：
如果沒看懂上面的式子，沒關系。通俗的講，在矩陣乘法中，結果 \(C\) 矩陣的第 \(i\) 行第 \(j\) 列的數，就是由矩陣 \(A\) 第 \(i\) 行 \(M\) 個數與矩陣 \(B\) 第 \(j\) 列 \(M\) 個數分別相乘再相加得到的。
矩陣乘法滿足結合律，不滿足一般的交換律。
利用結合律，矩陣乘法可以利用 快速冪 的思想來優化。
在比賽中，由於線性遞推式可以表示成矩陣乘法的形式，也通常用矩陣快速冪來求線性遞推數列的某一項。
優化
首先對於比較小的矩陣，可以考慮直接手動展開循環以減小常數。
可以重新排列循環以提高空間局部性，這樣的優化不會改變矩陣乘法的時間復雜度，但是會在得到常數級別的提升。
// 以下文的參考代碼為例
inline mat operator*(const mat& T) const {
mat res;
for (int i = 0; i < sz; ++i)
for (int j = 0; j < sz; ++j)
for (int k = 0; k < sz; ++k) {
res.a[i][j] += mul(a[i][k], T.a[k][j]);
res.a[i][j] %= MOD;
}
return res;
}
// 不如
inline mat operator*(const mat& T) const {
mat res;
int r;
for (int i = 0; i < sz; ++i)
for (int k = 0; k < sz; ++k) {
r = a[i][k];
for (int j = 0; j < sz; ++j)
res.a[i][j] += T.a[k][j] * r, res.a[i][j] %= MOD;
}
return res;
}
參考代碼
一般來說，可以用一個二維數組來模擬矩陣。
struct mat {
LL a[sz][sz];
inline mat() { memset(a, 0, sizeof a); }
inline mat operator-(const mat& T) const {
mat res;
for (int i = 0; i < sz; ++i)
for (int j = 0; j < sz; ++j) {
res.a[i][j] = (a[i][j] - T.a[i][j]) % MOD;
}
return res;
}
inline mat operator+(const mat& T) const {
mat res;
for (int i = 0; i < sz; ++i)
for (int j = 0; j < sz; ++j) {
res.a[i][j] = (a[i][j] + T.a[i][j]) % MOD;
}
return res;
}
inline mat operator*(const mat& T) const {
mat res;
int r;
for (int i = 0; i < sz; ++i)
for (int k = 0; k < sz; ++k) {
r = a[i][k];
for (int j = 0; j < sz; ++j)
res.a[i][j] += T.a[k][j] * r, res.a[i][j] %= MOD;
}
return res;
}
inline mat operator^(LL x) const {
mat res, bas;
for (int i = 0; i < sz; ++i) res.a[i][i] = 1;
for (int i = 0; i < sz; ++i)
for (int j = 0; j < sz; ++j) bas.a[i][j] = a[i][j] % MOD;
while (x) {
if (x & 1) res = res * bas;
bas = bas * bas;
x >>= 1;
}
return res;
}
};
應用
矩陣加速遞推
斐波那契數列（Fibonacci Sequence）大家應該都非常的熟悉了。在斐波那契數列當中，\(F_1 = F_2 = 1\)，\(F_i = F_{i - 1} + F_{i - 2}(i \geq 3)\)。
如果有一道題目讓你求斐波那契數列第 \(n\) 項的值，最簡單的方法莫過於直接遞推了。但是如果 \(n\) 的范圍達到了 \(10^{18}\) 級別，遞推就不行了，穩 TLE。考慮矩陣加速遞推。
設 \(Fib(n)\) 表示一個 \(1 \times 2\) 的矩陣 \(\left[ \begin{array}{ccc}F_n & F_{n-1} \end{array}\right]\)。我們希望根據 \(Fib(n-1)=\left[ \begin{array}{ccc}F_{n-1} & F_{n-2} \end{array}\right]\) 推出 \(Fib(n)\)。
試推導一個矩陣 \(\text{base}\)，使 \(Fib(n-1) \times \text{base} = Fib(n)\)，即 \(\left[\begin{array}{ccc}F_{n-1} & F_{n-2}\end{array}\right] \times \text{base} = \left[ \begin{array}{ccc}F_n & F_{n-1} \end{array}\right]\)。
怎么推呢？因為 \(F_n=F_{n-1}+F_{n-2}\)，所以 \(\text{base}\) 矩陣第一列應該是 \(\left[\begin{array}{ccc} 1 \\ 1 \end{array}\right]\)，這樣在進行矩陣乘法運算的時候才能令 \(F_{n-1}\) 與 \(F_{n-2}\) 相加，從而得出 \(F_n\)。同理，為了得出 \(F_{n-1}\)，矩陣 \(\text{base}\) 的第二列應該為 \(\left[\begin{array}{ccc} 1 \\ 0 \end{array}\right]\)。
綜上所述：\(\text{base} = \left[\begin{array}{ccc} 1 & 1 \\ 1 & 0 \end{array}\right]\) 原式化為 \(\left[\begin{array}{ccc}F_{n-1} & F_{n-2}\end{array}\right] \times \left[\begin{array}{ccc} 1 & 1 \\ 1 & 0 \end{array}\right] = \left[ \begin{array}{ccc}F_n & F_{n-1} \end{array}\right]\)
轉化為代碼，應該怎么求呢？
定義初始矩陣 \(\text{ans} = \left[\begin{array}{ccc}F_2 & F_1\end{array}\right] = \left[\begin{array}{ccc}1 & 1\end{array}\right], \text{base} = \left[\begin{array}{ccc} 1 & 1 \\ 1 & 0 \end{array}\right]\)。那么，\(F_n\) 就等於 \(\text{ans} \times \text{base}^{n-2}\) 這個矩陣的第一行第一列元素，也就是 \(\left[\begin{array}{ccc}1 & 1\end{array}\right] \times \left[\begin{array}{ccc} 1 & 1 \\ 1 & 0 \end{array}\right]^{n-2}\) 的第一行第一列元素。
注意，矩陣乘法不滿足交換律，所以一定不能寫成 \(\left[\begin{array}{ccc} 1 & 1 \\ 1 & 0 \end{array}\right]^{n-2} \times \left[\begin{array}{ccc}1 & 1\end{array}\right]\) 的第一行第一列元素。另外，對於 \(n \leq 2\) 的情況，直接輸出 \(1\) 即可，不需要執行矩陣快速冪。
為什么要乘上 \(\text{base}\) 矩陣的 \(n-2\) 次方而不是 \(n\) 次方呢？因為 \(F_1, F_2\) 是不需要進行矩陣乘法就能求的。也就是說，如果只進行一次乘法，就已經求出 \(F_3\) 了。如果還不是很理解為什么冪是 \(n-2\)，建議手算一下。
下面是求斐波那契數列第 \(n\) 項對 \(10^9+7\) 取模的示例代碼（核心部分）。
const int mod = 1000000007;
struct Matrix {
int a[3][3];
Matrix() { memset(a, 0, sizeof a); }
Matrix operator*(const Matrix &b) const {
Matrix res;
for (int i = 1; i <= 2; ++i)
for (int j = 1; j <= 2; ++j)
for (int k = 1; k <= 2; ++k)
res.a[i][j] = (res.a[i][j] + a[i][k] * b.a[k][j]) % mod;
return res;
}
} ans, base;
void init() {
base.a[1][1] = base.a[1][2] = base.a[2][1] = 1;
ans.a[1][1] = ans.a[1][2] = 1;
}
void qpow(int b) {
while (b) {
if (b & 1) ans = ans * base;
base = base * base;
b >>= 1;
}
}
int main() {
int n = read();
if (n <= 2) return puts("1"), 0;
init();
qpow(n - 2);
println(ans.a[1][1] % mod);
}
這是一個稍微復雜一些的例子。
我們發現，\(f_n\) 和 \(f_{n-1}, f_{n-2}, n\) 有關，於是考慮構造一個矩陣描述狀態。
但是發現如果矩陣僅有這三個元素 \(\begin{bmatrix}f_n& f_{n-1}& n\end{bmatrix}\) 是難以構造出轉移方程的，因為乘方運算和 \(+1\) 無法用矩陣描述。
於是考慮構造一個更大的矩陣。
我們希望構造一個遞推矩陣可以轉移到
轉移矩陣即為
矩陣表達修改
「THUSCH 2017」大魔法師
大魔法師小 L 制作了 \(n\) 個魔力水晶球，每個水晶球有水、火、土三個屬性的能量值。小 L 把這 \(n\) 個水晶球在地上從前向后排成一行，然后開始今天的魔法表演。
我們用 \(A_i,\ B_i,\ C_i\) 分別表示從前向后第 \(i\) 個水晶球（下標從 \(1\) 開始）的水、火、土的能量值。
小 L 計划施展 \(m\) 次魔法。每次，他會選擇一個區間 \([l, r]\)，然后施展以下 \(3\) 大類、\(7\) 種魔法之一：
魔力激發：令區間里每個水晶球中 特定屬性 的能量爆發，從而使另一個 特定屬性 的能量增強。具體來說，有以下三種可能的表現形式：
火元素激發水元素能量：令 \(A_i = A_i + B_i\)。
土元素激發火元素能量：令 \(B_i = B_i + C_i\)。
水元素激發土元素能量：令 \(C_i = C_i + A_i\)。
需要注意的是，增強一種屬性的能量並不會改變另一種屬性的能量，例如 \(A_i = A_i + B_i\) 並不會使 \(B_i\) 增加或減少。
魔力增強：小 L 揮舞法杖，消耗自身 \(v\) 點法力值，來改變區間里每個水晶球的 特定屬性 的能量。具體來說，有以下三種可能的表現形式：
- 火元素能量定值增強：令 \(A_i = A_i + v\)。
- 水元素能量翻倍增強：令 \(B_i=B_i \cdot v\)。
- 土元素能量吸收融合：令 \(C_i = v\)。
魔力釋放：小 L 將區間里所有水晶球的能量聚集在一起，融合成一個新的水晶球，然后送給場外觀眾。生成的水晶球每種屬性的能量值等於區間內所有水晶球對應能量值的代數和。需要注意的是，魔力釋放的過程不會真正改變區間內水晶球的能量。
值得一提的是，小 L 制造和融合的水晶球的原材料都是定制版的 OI 工廠水晶，所以這些水晶球有一個能量閾值 \(998244353\)。當水晶球中某種屬性的能量值大於等於這個閾值時，能量值會自動對閾值取模，從而避免水晶球爆炸。
小 W 為小 L（唯一的）觀眾，圍觀了整個表演，並且收到了小 L 在表演中融合的每個水晶球。小 W 想知道，這些水晶球蘊涵的三種屬性的能量值分別是多少。
由於矩陣的結合律和分配律成立，單點修改可以自然地推廣到區間，即推出矩陣后直接用線段樹維護區間矩陣乘積即可。
下面將舉幾個例子。
\(A_i = A_i + v\) 的轉移
\(B_i=B_i \cdot v\) 的轉移
「LibreOJ 6208」樹上詢問
有一棵 \(n\) 節點的樹，根為 \(1\) 號節點。每個節點有兩個權值 \(k_i, t_i\)，初始值均為 \(0\)。
給出三種操作：
\(\operatorname{Add}( x , d )\) 操作：將 \(x\) 到根的路徑上所有點的 \(k_i\leftarrow k_i + d\)
\(\operatorname{Mul}( x , d )\) 操作：將 \(x\) 到根的路徑上所有點的 \(t_i\leftarrow t_i + d \times k_i\)
\(\operatorname{Query}( x )\) 操作：詢問點 \(x\) 的權值 \(t_x\)
\(n,~m \leq 100000, ~-10 \leq d \leq 10\)
若直接思考，下放操作和維護信息並不是很好想。但是矩陣可以輕松地表達。
定長路徑統計
問題描述
給一個 \(n\) 階有向圖，每條邊的邊權均為 \(1\)，然后給一個整數 \(k\)，你的任務是對於所有點對 \((u,v)\) 求出從 \(u\) 到 \(v\) 長度為 \(k\) 的路徑的數量（不一定是簡單路徑，即路徑上的點或者邊可能走多次）。
我們將這個圖用鄰接矩陣 \(G\)（對於圖中的邊 \((u\to v)\)，令 \(G[u,v]=1\)，其余為 \(0\) 的矩陣；如果有重邊，則設 \(G[u,v]\) 為重邊的數量）表示這個有向圖。下述算法同樣適用於圖有自環的情況。
顯然，該鄰接矩陣對應 \(k=1\) 時的答案。
假設我們知道長度為 \(k\) 的路徑條數構成的矩陣，記為矩陣 \(C_k\)，我們想求 \(C_{k+1}\)。顯然有 DP 轉移方程
我們可以把它看作矩陣乘法的運算，於是上述轉移可以描述為
那么把這個遞推式展開可以得到
要計算這個矩陣冪，我們可以使用快速冪（二進制取冪）的思想，在 \(O(n^3 \log k)\) 的復雜度內計算結果。
定長最短路
問題描述
給你一個 \(n\) 階加權有向圖和一個整數 \(k\)。對於每個點對 \((u,v)\) 找到從 \(u\) 到 \(v\) 的恰好包含 \(k\) 條邊的最短路的長度。（不一定是簡單路徑，即路徑上的點或者邊可能走多次）
我們仍構造這個圖的鄰接矩陣 \(G\)，\(G[i,j]\) 表示從 \(i\) 到 \(j\) 的邊權。如果 \(i,j\) 兩點之間沒有邊，那么 \(G[i,j]=\infty\)。（有重邊的情況取邊權的最小值）
顯然上述矩陣對應 \(k=1\) 時問題的答案。我們仍假設我們知道 \(k\) 的答案，記為矩陣 \(L_k\)。現在我們想求 \(k+1\) 的答案。顯然有轉移方程
事實上我們可以類比矩陣乘法，你發現上述轉移只是把矩陣乘法的乘積求和變成相加取最小值，於是我們定義這個運算為 \(\odot\)，即
於是得到
展開遞推式得到
我們仍然可以用矩陣快速冪的方法計算上式，因為它顯然是具有結合律的。時間復雜度 \(O(n^3 \log k)\)。
限長路徑計數/最短路
上述算法只適用於邊數固定的情況。然而我們可以改進算法以解決邊數小於等於 \(k\) 的情況。具體地，考慮以下問題：
問題描述
給一個 \(n\) 階有向圖，邊權為 \(1\)，然后給一個整數 \(k\)，你的任務是對於每個點對 \((u,v)\) 找到從 \(u\) 到 \(v\) 長度小於等於 \(k\) 的路徑的數量（不一定是簡單路徑，即路徑上的點或者邊可能走多次）。
我們簡單修改一下這個圖，我們給每一個結點加一個權值為 \(1\) 的自環。這樣走的時侯就可以走自環，相當於原地走。這樣就包含了小於等於 \(k\) 的情況。修改后再做矩陣快速冪即可。（即使這個圖在修改之前就有自環，該算法仍是成立的）。
同樣的方法可以用於求邊數小於等於 \(k\) 的最短路，即加一個邊權為 \(0\) 的自環。
習題
- 洛谷 P1962 斐波那契數列，即上面的例題，同題 POJ3070
- 洛谷 P1349 廣義斐波那契數列，\(\text{base}\) 矩陣需要變化一下
- 洛谷 P1939【模板】矩陣加速（數列），\(\text{base}\) 矩陣變成了 \(3 \times 3\) 的矩陣，推導過程與上面差不多。