浅谈矩阵在 OI 中的初级应用

· · 算法·理论

琪零 · 前言

矩阵,是一种非常暴力的方式,它利用变换控制我们所求。

然而,矩阵入门有些困难,这也是笔者在此写下这篇文章的原因,笔者希望后面的人能更舒服地学会使用矩阵这一工具。

作为曾经的学习者,笔者建议在阅读本文前,拿好纸笔,时刻准备画图演算,尝试理解其原理。勿要匆匆求知而不解。

在本文中遇到习题时,除非您很清楚您已完全掌握此习题所含内容,否则不要粗略看它或者干脆略过,否则您有可能在阅读后文时遇到困难与误解。

当您在阅读一段中遇到了不理解、大脑过载(通常由于新定义过多导致)等情况,请关掉这一篇文章,看一会天花板,再尝试重新打开,也许可以冷却一下过热的大脑。

本文所介绍内容大多都能在 OI 中实际应用,所以部分与 OI(尤指 CSP-S/NOIP)关系较弱的区域可能不会讨论,敬请谅解。

笔者水平未必足够,如有疏漏,务必指出,谢谢!

闲话到此结束,我们直接进入正题。

其一 · 行列、和积、建立着

什么是矩阵呢?用 C++ 的角度说,就是一个二维数组。如下,是一个矩阵。

\begin{bmatrix} 4 & 7 \\ 9 & 8 \end{bmatrix}

行与列

在矩阵中,我们将横着的一行称作矩阵的一行,把竖着的一列称作矩阵的一列。切记,切记!

::::::info[习题1.1] 数以下矩阵的行数,列数。

\begin{bmatrix} 4 & 7&6 \\ 9 & 8&64 \end{bmatrix} \begin{bmatrix} 4 & 7 \\ 9 & 8\\6&7 \end{bmatrix} \begin{bmatrix} 4 & 7&2 \\ 9 & 8&8 \\1&8&0\end{bmatrix}

::::info[答案] 分别为 23 列,32 列,33 列。 :::: ::::::

为了方便,我们习惯使用大写字母表示矩阵,如 A,B,C;用下标依次表示矩阵中元素所在位置的的行数、列数(从 1 开始编号,且为整数)(后文不再添加行、列数的整数判定),如 A_{1,2} 对应上文正文中矩阵的值为 7

::::::info[习题1.2] 给定下面三个矩阵,分别记作 A,B,C,求A_{1,3},B_{3,1},C_{2,2},C_{3,3},C_{1,2}

\begin{bmatrix} 4 & 7&6 \\ 9 & 8&64 \end{bmatrix} \begin{bmatrix} 4 & 7 \\ 9 & 8\\6&7 \end{bmatrix} \begin{bmatrix} 4 & 7&2 \\ 9 & 8&3 \\1&6&0\end{bmatrix}

::::info[答案] 分别为 6,6,8,0,7。 :::: ::::::

现在,我们拥有了矩阵,自然要定义它们之间的运算。

当且仅当两个矩阵行数列数均相等时,这两个矩阵之间才有加法的定义。在本节,我们规定这两个矩阵的行数为 N,列数为 M

我们先严谨地定义一下,若 C=A+B,则有

\forall (1\le i\le N,1\le j \le M), C_{i,j}=A_{i,j}+B_{i,j}

啥玩意?!这是刚看这类式子的人们所容易疑问的。那就让我们通俗易懂地理解一下吧,矩阵加法就是把两个矩阵叠起来,对应位置相加,便可得到加法结果

如:

\begin{bmatrix} 4 & 7 \\ 9 & 8\\6&2 \end{bmatrix}+\begin{bmatrix} 8 & 9 \\ 7 & 4 \\1&3\end{bmatrix}=\begin{bmatrix} 12 & 16 \\ 16 & 12 \\7&5\end{bmatrix}

::::::info[习题1.3]

  1. \begin{bmatrix} 63 & 785 \\ 28 & 0\\123&77 \end{bmatrix}+\begin{bmatrix} 84 & 233 \\ 114 & 514\\1919&810 \end{bmatrix}=?

    不要害怕,此题不难。

  2. B2104 矩阵加法 ::::info[答案]
  3. \begin{bmatrix} 147 & 1018 \\ 142 & 514\\2042&887 \end{bmatrix}
  4. 看此题题解。 :::: ::::::

数乘

矩阵和一个的乘积乘坐它们之间的数乘,顺序一般无所谓,本文后面提到数乘时都会将矩阵放在前面,数放到后面。

矩阵数乘的结果就是对应位置依次乘上数。

直接上习题。

::::::info[习题1.4]

\begin{bmatrix} 1 & 2 \\ 3 & 4\\5&0\end{bmatrix}\times 7=?

::::info[答案]

\begin{bmatrix} 7 & 14 \\ 21 & 28\\35&0\end{bmatrix}

:::: ::::::

转置

* 本节无需过于掌握。

交换矩阵的行与列,即若 Bnm 列的矩阵 A 的转置,则

此时,记 B=A^T

如:

\begin{bmatrix} 4 & 7 \\ 9 & 8\\6&2 \end{bmatrix}^T=\begin{bmatrix} 4 & 9&6 \\ 7 & 8&2 \end{bmatrix}

::::::info[习题1.5]

  1. B2106 矩阵转置

::::info[答案]

  1. 见此题题解。 :::: ::::::

矩阵乘法

正·片·开·始

两个矩阵可以相乘,当且仅当第一个矩阵的列数等于第二个矩阵的行数。当两个矩阵可以相乘时,其结果行数等于第一个矩阵的行数,列数等于第二个矩阵的列数。

如:

\begin{bmatrix} 4 & 7 \\ 9 & 8\\6&2 \end{bmatrix}\cdot\begin{bmatrix} 8 & 9 \\ 7 & 4 \\1&3\end{bmatrix}=?

无意义!!!!!

\begin{bmatrix} 4 & 7 \\ 9 & 8\\6&2 \end{bmatrix}\cdot\begin{bmatrix} 8 & 9&2&6 \\ 7&5&1 & 4 \end{bmatrix}=?

等于

128& 121& 26&86 \\ 62 &64 &14& 44 \end{bmatrix}

有意义。

那么这个呢?

\begin{bmatrix} 8 & 9&2&6 \\ 7&5&1 & 4 \end{bmatrix}\cdot\begin{bmatrix} 4 & 7 \\ 9 & 8\\6&2 \end{bmatrix}=?

无意义!!!!!

同时

\\ 86& 91\end{bmatrix} 68&66\end{bmatrix}

从上面我们也可以注意到,矩阵乘法没有交换律

好了,见识这么多了,让我们看一下矩阵乘法的定义吧。

记矩阵 A 是一个 nr 列的矩阵,B 是一个 rm 列的矩阵。

那么,它们的乘积矩阵 C 满足:

\forall (1\le i\le n,1\le j \le m), C_{i,j}=\sum_{k=1}^{r}A_{i,k}B_{k,j}

记作

A\cdot B=C

AB=C

……
……
……

“这都是什么和什么啊!看不懂,根本看不懂啊!”

好了,现在让我们来把它弄得让人容易理解一些吧。

让我们把它当作 C++ 代码看看:

for(int i=1;i<=n;++i)
{
    for(int j=1;j<=m;++j)
    {
        C[i][j]=0;
        for(int k=1;k<=r;++k)
        {
            C[i][j]+=A[i][k]*B[k][j];
        }
    }
}

稍微理解一些了吗?

我们可以用下面的方式手算矩阵乘法:

如计算\begin{bmatrix} 4 & 7 \\ 9 & 8\\6&2 \end{bmatrix}\cdot\begin{bmatrix} 8 & 9&2&6 \\ 7&5&1 & 4 \end{bmatrix}

我们将右边的矩阵拉起,像下面一样画图:

对于答案矩阵的一个格子,它的值就是它所在的行、列对应的矩阵的一部分横向从左至右,纵向从上至下,依次相乘后相加所得结果

看不懂?让我们找出左上角格子对应的结果:

让它们依次相乘:4\times 8,7\times 7,然后相加,得到 81

再如第二行第三列的格子:

依次相乘:9\times 2,8\times 1

相加:9\times 2+8\times 1=26

如果还是有些懵,可以用一下这个经典矩阵乘法可视化工具。

来吧,做几道可爱的习题——

::::::info[习题1.6] 一、

\begin{bmatrix} 1&2\\3&4 \end{bmatrix}\cdot\begin{bmatrix} 4&3\\2&1\end{bmatrix}=? \begin{bmatrix} 1&2\\3&4\\5&6\end{bmatrix}\cdot\begin{bmatrix} 4\\2\end{bmatrix}=? \begin{bmatrix} 1&2\\3&4\\5&6\end{bmatrix}\cdot\begin{bmatrix} 4\\2\\7\end{bmatrix}=? \begin{bmatrix} 1&2\\3&4\\5&6\end{bmatrix}\cdot\begin{bmatrix} 1&0\\0&1\end{bmatrix}=? \begin{bmatrix} 1&0&0\\0&1&0\\0&0&1\end{bmatrix}\cdot\begin{bmatrix} 1&2\\3&4\\5&6\end{bmatrix}=?

注意判定乘法是否有意义

二、B2105 矩阵乘法 ::::info[答案] 一、

\begin{bmatrix} 8 & 5 \\ 20 & 13\end{bmatrix} \begin{bmatrix} 8 \\ 20 \\32\end{bmatrix} \texttt{无意义} \begin{bmatrix} 1&2\\3&4\\5&6\end{bmatrix} \begin{bmatrix} 1&2\\3&4\\5&6\end{bmatrix}

二、看此题题解。 :::: ::::::

在 OI 中,用计算机实现的矩阵乘法复杂度为 \operatorname{O}(nrm),对于方阵(指行列数相同的矩阵),即为\operatorname{O}(n^3)

单位矩阵

看一看习题 1.6 后两道题,我们发现每一小问中都有一个矩阵没变

在矩阵乘法运算中,有一些方阵 I 满足:当矩阵 AI 之间矩阵乘法有意义时

AI=A

I 与矩阵 A 之间矩阵乘法有意义时

IA=A

我们称这种矩阵为单位矩阵

单位矩阵有一种通用结构:

\begin{bmatrix}1&0&\cdots&0\\0&1&\cdots&0\\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&1\end{bmatrix}

\begin{bmatrix}1\end{bmatrix} \begin{bmatrix}1&0\\0&1\end{bmatrix} \begin{bmatrix}1&0&0\\0&1&0\\0&0&1\end{bmatrix} \begin{bmatrix}1&0&0&0\\0&1&0&0\\0&0&1&0\\0&0&0&1\end{bmatrix}

等。

单位矩阵的这些性质可以在后文中应用。

矩阵乘法的结合律

这是矩阵乘法可以用来做一些神奇操作的原因。

我们来证明一下,记矩阵 A 是一个 nr 列的矩阵,B 是一个 rc 列的矩阵, C 是一个 cm 列的矩阵:

\\ \begin{align} ((AB)C)_{i,j}&=\sum_{k=1}^{c}(AB)_{i,k}C_{k,j}\\ &=\sum_{k=1}^{c}(\sum_{d=1}^{r}A_{i,d}B_{d,k})C_{k,j}\\ &=\sum_{k=1}^{c}\sum_{d=1}^{r}A_{i,d}B_{d,k}C_{k,j} \\ &=\sum_{k=1}^{c}\sum_{d=1}^{r}A_{i,d}(B_{d,k}C_{k,j}) \\ &=\sum_{d=1}^{r}\sum_{k=1}^{c}A_{i,d}(B_{d,k}C_{k,j}) \\ &=\sum_{d=1}^{r}A_{i,d}\sum_{k=1}^{c}(B_{d,k}C_{k,j}) \\ &=\sum_{d=1}^{r}A_{i,d}(BC)_{d,j}\\ &=A(BC)_{i,j} \\ \end{align}

所以,矩阵乘法有结合律。 ::::::info[证明的较细解释] 提示:本折叠框内容由文心一言生成,经由笔者与 @busjcdf 巨佬的润色,不保证本折叠框内容正确性。

为证明矩阵乘法结合律 ((AB)C)_{i,j} = (A(BC))_{i,j},需验证以下推导步骤的严谨性:

步骤(1)到(2

依据:矩阵乘法定义。 矩阵 AB 的第 i 行第 k 列元素为 (AB)_{i,k} = \sum_{d=1}^{r} A_{i,d}B_{d,k}(其中 rA 的列数,也是 B 的行数)。

证明:将 (AB)_{i,k} 代入 ((AB)C)_{i,j} 的定义式:

\begin{aligned}((AB)C)_{i,j} &= \sum_{k=1}^{c} (AB)_{i,k}C_{k,j}\\ &= \sum_{k=1}^{c} \left( \sum_{d=1}^{r} A_{i,d}B_{d,k} \right) C_{k,j} \end{aligned}

等式成立。

步骤(2)到(3):

依据:双重求和的括号展开法则(加法对乘法的分配律)。

证明:展开内层括号后,双重求和的项与顺序不变:

\sum_{k=1}^{c} \left( \sum_{d=1}^{r} A_{i,d}B_{d,k} \right) C_{k,j} = \sum_{k=1}^{c} \sum_{d=1}^{r} A_{i,d}B_{d,k}C_{k,j}

等式成立。

步骤(3)到(4

依据:数的乘法结合律。

证明:对乘积 A_{i,d}B_{d,k}C_{k,j} 结合后两项:

\sum_{k=1}^{c} \sum_{d=1}^{r} A_{i,d}B_{d,k}C_{k,j} = \sum_{k=1}^{c} \sum_{d=1}^{r} A_{i,d}(B_{d,k}C_{k,j})

等式成立。

步骤(4)到(5

依据:有限项双重求和的交换律。

证明:交换求和指标 kd 的顺序,求和结果不变:

\sum_{k=1}^{c} \sum_{d=1}^{r} A_{i,d}(B_{d,k}C_{k,j}) = \sum_{d=1}^{r} \sum_{k=1}^{c} A_{i,d}(B_{d,k}C_{k,j})

等式成立。

步骤(5)到(6

依据:加法对乘法的分配律(提取公共因子)。

证明:内层求和 \sum_{k=1}^{c} B_{d,k}C_{k,j}A_{i,d} 无关,可提取至外层:

\sum_{d=1}^{r} \sum_{k=1}^{c} A_{i,d}(B_{d,k}C_{k,j}) = \sum_{d=1}^{r} A_{i,d} \sum_{k=1}^{c} B_{d,k}C_{k,j}

等式成立。

步骤(6)到(7

依据:矩阵乘法定义。

证明:内层求和 \sum_{k=1}^{c} B_{d,k}C_{k,j} 是矩阵 BC 的第 d 行第 j 列元素,即 (BC)_{d,j}

\sum_{d=1}^{r} A_{i,d} \left( \sum_{k=1}^{c} B_{d,k}C_{k,j} \right) = \sum_{d=1}^{r} A_{i,d}(BC)_{d,j}

等式成立。

步骤(7)到(8

依据:矩阵乘法定义。

证明:\sum_{d=1}^{r} A_{i,d}(BC)_{d,j} 是矩阵 A(BC) 的第 i 行第 j 列元素:

\sum_{d=1}^{r} A_{i,d}(BC)_{d,j} = (A(BC))_{i,j}

等式成立。

结论

综上,每一步推导均严格依据矩阵乘法定义、数的运算律及双重求和性质,因此:

( (AB)C)_{i,j} = (A(BC))_{i,j}


矩阵乘法结合律 (AB)C)_{i,j} = (A(BC))_{i,j} 得证。 :::::: 本节证明不必充分理解。

矩阵快速幂

本小节中矩阵只指方阵

对于一个 nn 列的方阵 A,如此定义矩阵的正整数 k 次幂:

A^k=\begin{cases}I,k= 0\\A^{k-1}\cdot A,k>0\end{cases}

其中 I 为一个 nn 列的单位矩阵。

k>0 时,其意义就是将矩阵 A 自乘 k 次。

那么,如何计算 A^{k} 呢?朴素的方式自然是乘 k 次,时间复杂度为 \operatorname{O}(n^3k),但是由于矩阵乘法有结合律,所以可以通过类似于快速幂的方法来计算(计算过程与正常快速幂一致),时间复杂度可以达到 \operatorname{O}(n^3\log{k})

你可能会疑问:我的习题呢?放到 2.1 去了。因为本模版不写一个类来实现难度较大,所以将其放入第二章节。

广义矩阵乘法

我们在矩阵乘法中,使用乘法加法进行运算。

那么,能不能使用其它运算代替乘法和加法呢?

可以!

但是为了保证矩阵乘法的结合律,我们联系上一节证明过程中 (3)(4)(5)(6) 步骤,可以发现我们所用的用来替代乘法和加法的运算(后文称其广义乘法、广义加法)一定要满足:

  1. 广义乘法对广义加法有分配律。
  2. 广义加法有交换律,结合律。
  3. 广义乘法有结合律。

顺带一提,如果运算本身不满足此性质,但是选取的元素有限制从而导致在此情况下运算律成立,也是可以的。如P6569 魔法值里的 \oplus -\times 矩阵乘法(在本文,后面使用这种形式表示矩阵乘法时,广义加法放到短线前,广义乘法放到短线后),因为一个矩阵恒为 01 矩阵(只含有 0,1 的矩阵),所以 \oplus 运算(按位异或)对数字乘法的分配律在此情况下成立,从而使得此矩阵乘法拥有了结合律。

另:普通的带有取模运算的加法与乘法也是满足这些关系的,因此在 OI 中可以在需要取模时正常使用矩阵。

上述不需彻底理解,记住即可。

::::::info[习题1.7] 判断(不需证明)以下定义的扩展矩阵乘法是否有结合律。

  1. ::::info[答案]

证明略。 :::: ::::::

其二 · 代码、成类、复用着

在后面,我们要大量使用矩阵做乘法,且会使用多个矩阵。因此,将矩阵成类、封装是必要的。

下面,我们就矩阵大小确定、矩阵大小不定两种情况分别讨论如何写一个矩阵类(或者说结构体)。

注:本节中矩阵只指方阵。

大小确定的矩阵

在结构体中,我们不需要存储矩阵大小,只需要存储矩阵本身。

那么,在此,矩阵就是一个封装后的二维数组。因此,下面的代码就足够了:

struct Mx
{
    int d[4][4];//3行3列,从1开始 
};

(注:矩阵从 0 开始存还是从 1 开始存看您心情,卡常严重或卡空间严重时请从 0 开始。)

但是,我们使用矩阵时,常常要为之赋初值,如零矩阵,单位矩阵,特殊构造的矩阵等。为此,我们可以写一个构造函数。

struct Mx
{
    int d[4][4];//3行3列,从1开始 
    Mx(int md)//如果您需要允许形如“Mx A;”结构存在,本行请使用“Mx(int md=0)”或“Mx(int md=1)”(视具体情况而定)。
    {
        if(md==0)
        {
            memset(d,0,sizeof(d));
        }
        if(md==1)
        {
            d[1][1]=d[2][2]=d[3][3]=1;
        }
    }
};

下面是选用内容:

我们经常要访问矩阵对应下标的元素,但是我们一直写A.d[i][j]也有点难受,因此我们可以再加一句这样的:

struct Mx
{
    int d[4][4];//3行3列,从1开始 
    Mx(int md)//如果您需要允许形如“Mx A;”结构存在,本行请使用“Mx(int md=0)”或“Mx(int md=1)”(视具体情况而定)。
    {
        if(md==0)
        {
            memset(d,0,sizeof(d));
        }
        if(md==1)
        {
            d[1][1]=d[2][2]=d[3][3]=1;
        }
    }
    int* const operator[](unsigned x)
    {
        return d[x];
    }
    const int* const operator[](unsigned x) const
    {
        return d[x];
    }
};

现在,我们就可以通过 A[i][j] 来访问对应下标的值了。

选用内容,可以不掌握。

大小不定的矩阵

加一个 n 记录大小即可。

关于乘法

在结构体外重载运算符 * 即可。

例:

struct Mx
{
    int d[4][4];//3行3列,从1开始 
    Mx(int md)//如果您需要允许形如“Mx A;”结构存在,本行请使用“Mx(int md=0)”或“Mx(int md=1)”(视具体情况而定)。
    {
        if(md==0)
        {
            memset(d,0,sizeof(d));
        }
        if(md==1)
        {
            d[1][1]=d[2][2]=d[3][3]=1;
        }
    }
    int* const operator[](unsigned x)
    {
        return d[x];
    }
    const int* const operator[](unsigned x) const
    {
        return d[x];
    }
};
Mx operator*(const Mx& A,const Mx& B)
{
    Mx C(0);
    for(int i=1;i<=3;++i)
    {
        for(int j=1;j<=3;++j)
        {
            for(int k=1;k<=3;++k)
            {
                C[i][j]+=A[i][k]*B[k][j];
            }
        }
    }
    return C;
}

::::::info[习题2.1] P3390 【模板】矩阵快速幂 ::::info[答案] 见此题题解。

:::: ::::::

其三 · 向量、状态、变换着

在本文仅考虑列向量,可以粗略地认为是只有一列的矩阵。

\begin{bmatrix}6\\4\\3\\8\end{bmatrix} 是一个含有 4 个元素的列向量。

矩阵与列向量之间的乘法

在 OI 中,列向量总是用来表示状态

相对应的,我们用方阵表示对状态的变换。

让一个方阵乘以一个列向量,其结果对应什么呢?

例如:

\begin{bmatrix}1&0&0&6\\4&6&5&4\\3&0&5&0\\8&2&9&0\end{bmatrix}\cdot\begin{bmatrix}6\\4\\3\\8\end{bmatrix}=\begin{bmatrix}54 \\95 \\33 \\83 \end{bmatrix}

我们使用上文所提到的可视化工具来看一看:

我们先把列向量翻上去

然后逐行盖住、相乘、累加,形成一个新的列向量。

这就是我们所需要理解的:列向量的对应项通过和列向量中元素之间的线性组合带进行变换,而变换是由矩阵描述的。实际上,某一行中矩阵中的数就是这一行对应列向量在变换后的值与变换前列向量中所有元素值线性关系的系数,也就是说,矩阵中每一行都是一个更新值的式子的系数(本段慢点读,不理解可以看下文的例子)。

如:

\begin{bmatrix}A_{1,1}&A_{1,2}&A_{1,3}&A_{1,4}\\A_{2,1}&A_{2,2}&A_{2,3}&A_{2,4}\\A_{3,1}&A_{3,2}&A_{3,3}&A_{3,4}\\A_{4,1}&A_{4,2}&A_{4,3}&A_{4,4}\end{bmatrix}\cdot\begin{bmatrix}a_1\\a_2\\a_3\\a_4\end{bmatrix}=\begin{bmatrix}b_1\\b_2\\b_3\\b_4 \end{bmatrix}

等价于

A_{2,1}\cdot a_1+A_{2,2}\cdot a_2+A_{2,3}\cdot a_3+A_{2,4}\cdot a_4=b_2\\ A_{3,1}\cdot a_1+A_{3,2}\cdot a_2+A_{3,3}\cdot a_3+A_{3,4}\cdot a_4=b_3\\ A_{4,1}\cdot a_1+A_{4,2}\cdot a_2+A_{4,3}\cdot a_3+A_{4,4}\cdot a_4=b_4 \end{cases}

至此,本文终于进入了 OI 部分。

::::::info[习题3.1] 改写式子为矩阵乘列向量的形式。

1.

A_{2,1}\cdot a_1+A_{2,2}\cdot a_2+A_{2,3}\cdot a_3=b_2\\ A_{3,1}\cdot a_1+A_{3,2}\cdot a_2+A_{3,3}\cdot a_3=b_3 \end{cases}
a_2=b_2\\ 7\cdot a_2+9\cdot a_3=b_3 \end{cases}

其中, v 为常数。 ::::info[答案]

  1. \begin{bmatrix}A_{1,1}&A_{1,2}&A_{1,3}\\A_{2,1}&A_{2,2}&A_{2,3}\\A_{3,1}&A_{3,2}&A_{3,3}\end{bmatrix}\cdot\begin{bmatrix}a_1\\a_2\\a_3\end{bmatrix}=\begin{bmatrix}b_1\\b_2\\b_3 \end{bmatrix}
  2. \begin{bmatrix}5&0&v\\0&1&0\\0&7&9\end{bmatrix}\cdot\begin{bmatrix}a_1\\a_2\\a_3\end{bmatrix}=\begin{bmatrix}b_1\\b_2\\b_3 \end{bmatrix}

:::: ::::::

矩阵加速 DP

上一节中,我们说过列向量用来表示状态,方阵表示对状态变换的系数描述。

那么什么会大量使用状态转移呢?动态规划(DP)。

让我们从一道典例入手吧——

P1962 斐波那契数列 。

题意省流:求斐波那契数列第 N 项(即 F_n)对 10^9+7 取模后的结果,N< 2^{63}

我们知道,对于大于 2 的正整数 i,有:

F_i=F_{i-1}+F_{i-2}

所以,当我们将状态从 F_{i-1} 转到 F_i 时需要 F_{i-2} 的参与。

那么,我们可以让 F_{i-1},F_{i-2} 同时作为状态的一部分,即进入列向量中

至此,一次转移变成为了——

\begin{bmatrix}F_{i-1}\\F_{i-2}\\\end{bmatrix}\to\begin{bmatrix}F_i\\F_{i-1} \end{bmatrix}

由数列转移的关系式,有

F_{i-1}+F_{i-2}=F_{i} F_{i-1}=F_{i-1}

正如上一节那样,我们可以将这几个递推式写成矩阵形式,即

\begin{bmatrix}1&1\\1&0\\\end{bmatrix}\cdot\begin{bmatrix}F_{i-1}\\F_{i-2}\\\end{bmatrix}\to\begin{bmatrix}F_i\\F_{i-1} \end{bmatrix}

一次转移可以让列向量第一项变成数列下一项,而我们已知 F_1=F_2=1, 可得

\begin{bmatrix}F_n\\F_{n-1} \end{bmatrix}=\begin{bmatrix}1&1\\1&0\\\end{bmatrix}^{n-2}\cdot\begin{bmatrix}F_{2}\\F_{1}\\\end{bmatrix}=\begin{bmatrix}1&1\\1&0\\\end{bmatrix}^{n-2}\cdot\begin{bmatrix}1\\1\\\end{bmatrix}

计算完以后取最终向量第一项便是答案。

此题做完了……吗?

没有!当 n<2n=1 时次数变为 -1,超出上文讨论范围,因此特判输出 1

现在,我们可以做出来了。

通过上述,我们不难发现,状态列向量的选取可以有冗余,即可以有不参与输出的结果,如上文的 F_{i-2}。这一点在后续会经常使用。

::::::info[习题3.2] 给定递推式,求矩阵次幂、列向量表示递推结果(不必写出特殊/边界情况取值)。

1.

F_{i}=\begin{cases}i,1\le i\le 4\\F_{i-1}-F_{i-2}+3F_{i-3}+2F_{i-4}\end{cases}
b_{i-1}=b_{i}\\ 7\cdot b_{i-1}+9\cdot c_{i-1}=c_{i} \end{cases}

其中,v 为常数,a_1,b_1,c_1 均等于 1

b_{i-1}+1=b_{i}\\ 7\cdot b_{i-1}+9\cdot c_{i-1}=c_{i} \end{cases}

其中, v 为常数,a_1,b_1,c_1 均等于 1

  1. P1939 矩阵加速(数列)

  2. P1962 斐波那契数列 ::::info[答案]

  3. \begin{bmatrix}F_n\\F_{n-1}\\F_{n-2}\\F_{n-3} \end{bmatrix}=\begin{bmatrix}1&-1&3&2\\1&0&0&0\\0&1&0&0\\0&0&1&0\\\end{bmatrix}^{n-4}\cdot\begin{bmatrix}4\\3\\2\\1\\\end{bmatrix}
  4. \begin{bmatrix}a_n\\b_n\\c_n \end{bmatrix}=\begin{bmatrix}5&0&v\\0&1&0\\0&7&9\\\end{bmatrix}^{n-1}\cdot\begin{bmatrix}a_1\\b_1\\c_1\\\end{bmatrix}=\begin{bmatrix}5&0&v\\0&1&0\\0&7&9\\\end{bmatrix}^{n-1}\cdot\begin{bmatrix}1\\1\\1\\\end{bmatrix}
  5. \begin{bmatrix}a_n\\b_n\\c_n \end{bmatrix}=\begin{bmatrix}5&0&v&0\\0&1&0&1\\0&7&9&0\\0&0&0&1\\\end{bmatrix}^{n-1}\cdot\begin{bmatrix}a_1\\b_1\\c_1\\1\\\end{bmatrix}=\begin{bmatrix}5&0&v&0\\0&1&0&1\\0&7&9&0\\0&0&0&1\\\end{bmatrix}^{n-1}\cdot\begin{bmatrix}1\\1\\1\\1\\\end{bmatrix}

在此,我们切实见识到了冗余状态的用处:常数 1 也可以进入向量作为状态一部分。

不过,也不能不论是否需要都加上,否则矩阵乘法 3 次方的常数可不是闹着玩的。

  1. 见对应题目题解。
  2. 见对应题目题解。 :::: ::::::

从 DP 到图论

图,可以使用邻接矩阵来存储。

我们自然想到:这个矩阵的 k 次幂是什么含义呢?

当然,这里使用 \min - + 矩阵乘法,易证,这种矩阵乘法有结合律。

也就是当

C=A\cdot B

时,有

C_{i,j}=\min_{k=1}^{r}(A_{i,k}+B_{k,j})

对于一个有 n 个点的邻接矩阵 A,有

(A\cdot A)_{i,j}=\min_{k=1}^{n}(A_{i,k}+B_{k,j})

这是什么意思呢?从点 i 和到点 j 之间经过恰好一个点,即经过 2 条边时从点 i 到点 j 的最短路长度。

通过递推,容易得到:

(A^k)_{i,j} 表示从点 i 到点 j 之间经过恰好 k 条边时从点 i 到点 j 的最短路长度。

::::::info[习题3.3] 没找到题,自己找去。

::::::

从图论回到 DP

类似地,我们重新定义矩阵乘法,在邻接矩阵上快速幂,也是一种绝妙的方法。

直接上习题,过程参考上一节。

::::::info[习题3.4]

  1. 灯灯开会(简单版)

::::info[答案]

  1. 略。 :::: ::::::

二进制优化矩阵快速幂

我们在使用上一部分的内容时,如果遇到多次询问(q 次),时间复杂度会到达 \operatorname{O}(qn^3\log{k}),当 q,n\approx 100 时,程序就 TLE 了。

为了降低询问开销,我们可以换一个思路。

如:用于变换的矩阵为 An 行 方阵),初始的列向量为 \textup{\textbf{\textsf{d}}},那么,最终的列向量为

A^{2^{k_{1}}}\cdot A^{2^{k_{2}}}\cdot A^{2^{k_{3}}}\cdots A^{2^{k_{p}}}\cdot \textup{\textbf{\textsf{d}}}

其中

\sum_{i=1}^{p}2^{k_i}=k

暴力快速幂是这么计算的:

(A^{2^{k_{1}}}\cdot A^{2^{k_{2}}}\cdot A^{2^{k_{3}}}\cdots A^{2^{k_{p}}}) \cdot \textup{\textbf{\textsf{d}}}

耗时都在前半部分。

那么,如果我们这样计算呢?

A^{2^{k_{1}}}\cdot (A^{2^{k_{2}}}\cdot (A^{2^{k_{3}}}(\cdots (A^{2^{k_{p}}}\cdot \textup{\textbf{\textsf{d}}}))))

每次计算都只是一个矩阵和一个向量的乘法,回看第一章矩阵乘法一节,可以知道单次乘法时间复杂度为 \operatorname{O}(n^2)。那么,计算 \log_2{k} 次便是 \operatorname{O}(n^2\log{k})q 次询问总时间复杂度为 \operatorname{O}(qn^2\log{k})

但是,我们如何知晓 A^{2^{k_{i}}} 们的值呢?

很简单,预处理一下就行了。k_i 的量级是 \log{k} 的,直接暴力相乘,时间复杂度 \operatorname{O}(n^3\log{k})

这样,总时间复杂度就变为了 \operatorname{O}((n+q)n^2\log{k})

::::::info[习题3.5]

  1. P6569 魔法值

::::info[答案]

  1. 见此题题解。 :::: ::::::

矩阵与数据结构

在做数据结构题时,有时标记的更新顺序、数量过多,难以调控。一种方法是用矩阵来作为标记。

直接上习题,解法看题解。

::::::info[习题3.6]

  1. P7453 大魔法师

::::info[答案]

  1. 见此题题解。 :::: ::::::

其四 · 调换、寻解、消元着

看了标题,大家应该明白这一章要讲什么了:高斯消元法。当然,准确来说,是高斯-约旦消元法。

实际上,本章可以当做P2455 线性方程组的题解。

I 矩阵化

先看方程本体:

\begin{cases} a_{1, 1} x_1 + a_{1, 2} x_2 + \cdots + a_{1, n} x_n = b_1 \\ a_{2, 1} x_1 + a_{2, 2} x_2 + \cdots + a_{2, n} x_n = b_2 \\ \cdots \\ a_{n,1} x_1 + a_{n, 2} x_2 + \cdots + a_{n, n} x_n = b_n \end{cases}

使用第三章介绍的方法,可以变形为:

\cdot \begin{bmatrix}x_1\\x_2\\\vdots\\x_n\end{bmatrix} = \begin{bmatrix}b_1\\b_2\\\vdots\\b_n\end{bmatrix}

为了方便记录,我们可以使用增广矩阵,也就是将等号右侧的向量放到左侧:

\left[\begin{array}{cccc|c}a_{1,1}&a_{1,2}&\cdots&a_{1,n}&b_1\\a_{2,1}&a_{2,2}&\cdots&a_{2,n}&b_2\\\vdots&\vdots&\ddots&\vdots &\vdots\\a_{n,1}&a_{n,2}&\cdots&a_{n,n}&b_n\end{array}\right]

注意:上述三种表示方式等价!

后面会在这三种形式内左右横跳。

II 如何消元?

先说几个引理。

引理 1

交换增广矩阵两行,方程解不变。

证明:把它变换成方程组形式,交换对应方程,重新变换成矩阵,过程中未改变解。

引理 2

方程组之间加减,结果不变。

证明:回忆一下初中学的加减消元法。

引理 3

若增广矩阵左侧方阵部分变为单位矩阵,则此时的 x_i=b_i

证明:读者自证不难,故略(若不解,请回看第三章第一节)。

引理 4

若增广矩阵有一行使得左侧方阵部分全为 0,右侧不为 0,则方程组无解。

证明:变换增广矩阵为方程组形式,则有一行形如 0=b_i,b_i\ne 0,显然无解。

引理 5

若增广矩阵有一行使得左侧方阵部分全为 0,右侧也为 0,则方程组有无数组解。

证明:变换增广矩阵为方程组形式,则有一行形如 0=0,则此行无效。而方程组有 n 个未知数,有少于 n 个有效方程组,则一定有多解。

解方程的过程

那么,如何消元呢?

我们的目的,是将矩阵削成单位矩阵,然后使用引理 3 便可得出结论。

于是,我们找一行,将在此我们想要把系数化为 1 的未知数(其余则想要削到 0)削到 1(将增广矩阵的这一行同除以它的系数),然后让这个新生的方程和其它方程加减消元(引理 2)。

选择谁呢?为了提高精度,我们可以选择系数绝对值最大的那一行方程,通过引理 1,交换至对应单位矩阵此列为 1 的那一行,然后消元即可。

没听懂?看题解去。

我们给一个例子(感谢 @busjcdf 巨佬的倾力相助,这组数据由这位大佬提供)吧。

\begin{cases} 2x + y - z = 8 \\ -3x - y + 2z = -11 \\ -2x + y + 2z = -3 \end{cases}

让我们看一下这个方程组,先将它化成增广矩阵:

\left[\begin{array}{ccc|c} 2 & 1 & -1 & 8 \\ -3 & -1 & 2 & -11 \\ -2 & 1 & 2 & -3 \end{array}\right]

取第一列为主元,取绝对值最大者作为第一行:

\left[\begin{array}{ccc|c} -3 & -1 & 2 & -11 \\ 2 & 1 & -1 & 8 \\ -2 & 1 & 2 & -3 \end{array}\right] \begin{cases} -3x - y + 2z = -11 \\ 2x + y - z = 8 \\ -2x + y + 2z = -3 \end{cases}

此行主元化为一:

\left[\begin{array}{ccc|c} 1 & \frac{1}{3} & -\frac{2}{3} & \frac{11}{3} \\ 2 & 1 & -1 & 8 \\ -2 & 1 & 2 & -3 \end{array}\right] \begin{cases} x+\frac{1}{3} y -\frac{2}{3}z= \frac{11}{3} \\ 2x + y - z = 8 \\ -2x + y + 2z = -3 \end{cases}

用这一行消去其余两行的 x(类似于加减消元法):

\begin{cases} x+\frac{1}{3} y -\frac{2}{3}z= \frac{11}{3} \\ \frac{1}{3}y + \frac{1}{3}z = \frac{2}{3} \\ \frac{5}{3} y + \frac{2}{3}z = \frac{13}{3} \end{cases} \left[\begin{array}{ccc|c} 1 & \frac{1}{3} & -\frac{2}{3} & \frac{11}{3} \\ 0 & \frac{1}{3} & \frac{1}{3} & \frac{2}{3} \\ 0 & \frac{5}{3} & \frac{2}{3} & \frac{13}{3} \end{array}\right]

然后,轮回下去——选取第二列作为主元,选取已经被选过的第 1 列以外的行中第二列的值的绝对值最大的行作为新的第二行,并消元。

\left[\begin{array}{ccc|c} 1 & \frac{1}{3} & -\frac{2}{3} & \frac{11}{3} \\ 0 & \frac{5}{3} & \frac{2}{3} & \frac{13}{3}\\ 0 & \frac{1}{3} & \frac{1}{3} & \frac{2}{3} \\ \end{array}\right] \begin{cases} x+\frac{1}{3} y -\frac{2}{3}z= \frac{11}{3} \\ \frac{5}{3} y + \frac{2}{3}z = \frac{13}{3}\\ \frac{1}{3}y + \frac{1}{3}z = \frac{2}{3} \\ \end{cases} \left[\begin{array}{ccc|c} 1 & \frac{1}{3} & -\frac{2}{3} & \frac{11}{3} \\ 0 & 1 & \frac{2}{5} & \frac{13}{5}\\ 0 & \frac{1}{3} & \frac{1}{3} & \frac{2}{3} \\ \end{array}\right] \begin{cases} x+\frac{1}{3} y -\frac{2}{3}z= \frac{11}{3} \\ y + \frac{2}{5}z = \frac{13}{5}\\ \frac{1}{3}y + \frac{1}{3}z = \frac{2}{3} \\ \end{cases} \begin{cases} x-\frac{4}{5}z= \frac{14}{5} \\ y + \frac{2}{5}z = \frac{13}{5}\\ \frac{1}{5}z = -\frac{1}{5} \\ \end{cases} \left[\begin{array}{ccc|c} 1 & 0 & -\frac{4}{5} & \frac{14}{5} \\ 0 & 1 & \frac{2}{5} & \frac{13}{5}\\ 0 & 0& \frac{1}{5} & -\frac{1}{5} \\ \end{array}\right]

选取第三列作为主元,选取已经被选过的第 1,2 列以外的行中第三列的值的绝对值最大的行作为新的第三行,并消元。

\left[\begin{array}{ccc|c} 1 & 0 & -\frac{4}{5} & \frac{14}{5} \\ 0 & 1 & \frac{2}{5} & \frac{13}{5}\\ 0 & 0& \frac{1}{5} & -\frac{1}{5} \\ \end{array}\right] \begin{cases} x-\frac{4}{5}z= \frac{14}{5} \\ y + \frac{2}{5}z = \frac{13}{5}\\ \frac{1}{5}z = -\frac{1}{5} \\ \end{cases} \left[\begin{array}{ccc|c} 1 & 0 & -\frac{4}{5} & \frac{14}{5} \\ 0 & 1 & \frac{2}{5} & \frac{13}{5}\\ 0 & 0& 1 & -1 \\ \end{array}\right] \begin{cases} x-\frac{4}{5}z= \frac{14}{5} \\ y + \frac{2}{5}z = \frac{13}{5}\\ z = -1 \\ \end{cases} \begin{cases} x= 2 \\ y = 3\\ z = -1 \\ \end{cases} \left[\begin{array}{ccc|c} 1 & 0 & 0 & 2 \\ 0 & 1 & 0 & 3\\ 0 & 0& 1 & -1 \\ \end{array}\right]

因此,解为

\begin{cases} x= 2 \\ y = 3\\ z = -1 \\ \end{cases}

显然,上述方法时间复杂度为 \operatorname{O}(n^3)

III 如何判无解/无数解?

如果有一次未知数那一列都为 0 ,那么就在最后依靠引理 4,5 判一下即可。

结束了……假的。

注意优先级:无解高于无数解,即只要有一个无解的方程,无论其余方程是什么结果,都是无解的。

另外,如果发现这一列系数全为 0,那么要继续找,且要考虑是否在前面会有方程在后续被使用。

强烈建议自己写一遍!

::::::info[习题4.1]

  1. P3389 高斯消元法
  2. P2455 线性方程组

::::info[答案]

  1. 见此题题解。
  2. 见此题题解。 :::: ::::::

后记&鸣谢

这篇是我的随笔,可能有误,敬请谅解。

感谢 @busjcdf 巨佬的倾力相助,这位大佬帮助我优化了文风,校验了文本,处理了图片,提供了样例。

感谢文心一言在矩阵乘法的结合律一节中的证明详解部分的帮助。