【学习笔记】矩阵及其应用

滑_稽 · 2021-08-21 22:59:19 · 个人记录

目录

一、什么是矩阵？

二、矩阵的基本运算

三、矩阵快速幂

四、实际应用

五、模板代码

六、参考文章

一、什么是矩阵？

在数学中，矩阵（Matrix）是一个按照长方阵列排列的复数或实数集合，最早来自于方程组的系数及常数所构成的方阵。这一概念由19世纪英国数学家凯利首先提出。

——百度百科 - 矩阵

定义一个 n\times m(1\le n,1\le m) 的矩阵

其中，n 叫做矩阵 A 的行数，m 叫做矩阵 A 的列数，a_{i,j}(\forall i\in[1,n],\forall j\in[1,m]) 叫做矩阵 A 的元素。矩阵的元素一般为数。

只有一行的矩阵叫做行向量，只有一列的矩阵叫做列向量。

行数和列数相同的矩阵叫做方阵。可以用“n 阶方阵”描述一个 n\times n 的方阵。

举个例子，有矩阵

该矩阵的行数为 2，列数为 3。第 2 行第 1 列的元素 A_{2,1}=4。

我们一般用 A,B,C 等大写字母表示一个矩阵。如果一个矩阵是 n\times m 的矩阵，则说明它的行数是 n，列数是 m。

二、矩阵的基本运算

矩阵的加法&减法

只要两个矩阵的行数相同且列数相同，它们就可以进行加法和减法。具体做法是，将两个矩阵对应位置上的数相加减。例如对于两个 n\times m 的矩阵 A 和 B，设 C=A+B，则

举个例子，若

则

矩阵的数乘

任意一个矩阵都能进行数乘。具体做法是，将矩阵的每一个元素都乘上一个数。例如对于一个 n\times m 的矩阵 A，若 C=\lambda A，则

举个例子，若

则

矩阵的乘法

首先给出一条定理：矩阵乘法满足结合律和分配律，而不满足交换律。

也就是说，(A\times B)\times C=A\times(B\times C)，(A+B)\times C=A\times C+B\times C，而 A\times B\ne B\times A。

有了这条定理，我们再来说矩阵乘法的性质。

两个矩阵 A 和 B 能进行乘法运算得到 C=A\times B，当且仅当 A 的列数等于 B 的行数。注意，这里提到的乘积是 A\times B 而不是 B\times A，就是因为矩阵乘法不满足交换律。如果满足条件，得到的乘积矩阵 C 的行数一定等于 A 的行数，C 的列数一定等于 B 的列数。设 A 是一个 n\times m 的矩阵，B 是一个 m\times p 的矩阵，则 C 是一个 n\times p 的矩阵，满足：

如果还是不懂矩阵乘法是怎么一回事，没有关系，请看解释：

现有两个矩阵 A 和 B，它们的具体信息如下：

要求 C=A\times B，首先可得 C 的规模为 3 行 2 列：

然后我们手动模拟一下矩阵乘法的步骤，以 C_{1,1} 为例。

由上面的公式可得，C_{1,1}=A_{1,1}B_{1,1}+A_{1,2}B_{2,1}。

如上，红色的元素之积加蓝色的元素之积就等于 C_{1,1}。

再来求 C_{1,2}=A_{1,1}B_{1,2}+A_{1,2}B_{2,2}：

如上，红色的元素之积加蓝色的元素之积就等于 C_{1,2}。

……

如此这般，我们便可以求出 C：

概括来说，C_{i,j} 就等于 A 第 i 行的 m 个数与 B 第 j 列的 m 个数分别相乘得到的 m 个积的总和。

这里另外提一种特殊的矩阵：单位矩阵。一个 n\times m 矩阵 A 为单位矩阵，当且仅当它满足：

即单位矩阵一定是个方阵，它的主对角线（即从左上角到右下角的对角线）上的所有元素都为 1，其余元素都为 0。

下面这个矩阵就是一个单位矩阵：

单位矩阵有着如下性质：假设有单位矩阵 A，任何能与 A 相乘的矩阵 B 乘 A 后得到的积还是 B。

矩阵的幂

一个矩阵能进行幂运算，当且仅当它是一个方阵（想一想，为什么）。矩阵 A 的幂 A^k 通常利用矩阵快速幂求出。

三、矩阵快速幂

矩阵快速幂的思想与数的快速幂的思想可以说是完全一致，除了将乘法替换成了矩阵乘法以外。

数的快速幂

如果你还不知道什么是数的快速幂，请看本节；反之可以跳过本节，直奔主题。

问题：给定数 a,n，要求 a^n，大多数情况还要将答案模 p。参见luoguP1226

要解决这个问题，首先可以将 n 二进制拆分。

设 n=(n_{k-1}n_{k-2}\dots n_{0})_2，即在二进制表示下有 k 位，则 n=2^{k-1}n_{k-1}+2^{k-2}n_{k-2}+\dots2^0n_0。那么，a^n 就等价于

也就是

也就是

由于 n_i\in\{0,1\}(1\le i<k)，因此当 n_i 为 0 时，直接跳过；当 n_i 为 1 时，将答案乘上 a^{2^i}即可。 另外，不难发现 a^{2^i}=a^{2^{i-1}\times 2}={\left(a^{2^{i-1}}\right)}^2，于是就有了如下这个 O(\log n) 的算法。

int Power(int a,int n){
    int ans=1;
    for(;n;n>>=1,a=a*a%k)
        if(n&1)//取出n在二进制表示下的的第0位
            ans=ans*a%k;
}

这种叫做快速幂的算法还有一种递归版本：

int Power(int a,int n){
    if(!n) return 1;
    if(n&1) return Power(a*a%p,n/2)*a%p;
    return Power(a*a%k,n/2);
}

它的正确性在此不再作赘述。理论上递归版本和非递归版本的时间复杂度相同，但实际操作告诉我们，很多时候非递归版本更快一些，因为递归会花费一定的开销。

矩阵快速幂

再看一遍数的快速幂的模板题。

问题：给定数 a,n，要求 a^n，大多数情况还要将答案模 p。

将这个问题稍微改动一下：

问题：给定 n\times n 的矩阵 A，数 k，要求 A^k，大多数情况还要将答案的每一个元素模 p。

这便是矩阵快速幂的模板题。

沿着数的快速幂的思想，我们可以写出这样一条式子：

其中 i 是 k 在二进制表示下的位数。由于矩阵乘法满足结合律，因此只要在代码中定义了矩阵乘法，我们完全可以基本照搬数的快速幂。实现代码在后面有给出。

时间复杂度 O(nmp\log k)，即矩阵乘法的时间复杂度乘上快速幂的时间复杂度。

四、实际应用

斐波那契数列

众所周知，斐波那契数列是一个形如

的数列。其中第 0 项 F_0=0，第 1 项 F_1=1，第 n 项 F_n=F_{n-1}+F_{n-2}(n\ge2,n\in\Bbb N^*)。

要求该数列的第 n 项 F_n，可以直接利用递推式求解，时间复杂度 O(n)，但这仅限于 n 较小的时候。如果 n 是一个很大的数，递推求解会超时，此时就需要一种 O(\log n) 级别的算法。

还记得矩阵快速幂的时间复杂度吗？

时间复杂度 O(nmp\log k)，即矩阵乘法的时间复杂度乘上快速幂的时间复杂度。

观察到矩阵快速幂的时间复杂度正好是 O(\log n) 级别的，于是我们思考如何利用矩阵快速幂求出 F_n。

构造如下两个矩阵（列向量）：

设矩阵 A 满足：

根据矩阵乘法的定义，A 是一个 2 行 2 列的方阵。

令

则式 (1) 可变为

由矩阵乘法，可列出如下方程组：

容易发现，当 \begin{cases}a=1,\\b=1,\\c=1,\\d=0\end{cases} 时，方程组 (2) 成立。

于是，有

据此可得

以此类推，得：

至此，我们完成了对式 (3) 的推导。式 (3) 即为经典的利用矩阵快速幂求斐波那契数列第 n 项的公式，可用矩阵快速幂求得 \begin{bmatrix}1&1\\1&0\end{bmatrix}^{n-1}。

（另：如果你不怕麻烦，你可以将 \begin{bmatrix}F_1\\F_0\end{bmatrix} 替换成 \begin{bmatrix}F_2\\F_1\end{bmatrix} 乃至更高，再根据你改的矩阵更改 \begin{bmatrix}1&1\\1&0\end{bmatrix} 的幂次即可。注意更改过后的矩阵元素不能超过 F_n。例如当 n\ge4 时，下面这个式子也被认为是正确的。但通常我们不会在这种细节上花时间，除非你想卡常）

时间复杂度为一次矩阵快速幂复杂度加上一次矩乘复杂度，为 O(2^2+(2^3\log n))，可视为 O(\log n)。

由上可见，矩阵乘法的终极奥义是构造矩阵。

例题P1962AC代码（仅展示主程序部分，Matrix结构体实现代码会在下面放出）：

#include<iostream>
#include<cstdio>
#define int long long
using namespace std;

const int MOD=1e9+7;
int n;

struct Matrix{
    ...
}A(2,1),B(2,1),C(2,2);

signed main(){
    ios::sync_with_stdio(false);
    cin>>n;
    A.a[1][1]=1,A.a[2][1]=0;
    C.a[1][1]=C.a[1][2]=C.a[2][1]=1,C.a[2][2]=0;
    B=C.power(n-1)*A;
    cout<<B.a[1][1]<<endl;
    return 0;
}

五、模板代码

struct Matrix{
    int r,c,a[maxN][maxN];
    Matrix(int rr,int cc,bool make_identity_matrix=0):r(rr),c(cc){
    //有了这个构造函数，定义一个矩阵A的代码就是：
    //Matrix A(行数,列数,是否定义一个单位矩阵（是填1，不是填0）);
        for(int i=1;i<=r;i++)
            for(int j=1;j<=c;j++){
                if(make_identity_matrix && i==j) a[i][j]=1;
                else a[i][j]=0;
            }
    }
    /*
    int* operator[](int i) {return a[i];} //用A[i][j]表示A.a[i][j]
    const int* operator[](int i) const {return a[i];} //C++特性
    */
    //上面两行是我从yggdyy_大佬的blog里复制过来的（链接详见【参考文章】），有需要可以去掉注释
    Matrix operator =(const Matrix &b){//重载赋值符
        for(int i=1;i<=r;i++)
            for(int j=1;j<=c;j++)
                a[i][j]=b.a[i][j];
        return *this;
    }
    Matrix operator +(const Matrix &b) const{//重载加法运算符
        Matrix ans(r,c);
        for(int i=1;i<=r;i++)
            for(int j=1;j<=c;j++)
                ans.a[i][j]=a[i][j]+b.a[i][j];
        return ans;
    }
    Matrix operator +=(const Matrix &b){//重载+=运算符
        return (*this)=(*this)+b;
    }
    void numtime(int x){//定义数乘函数
        for(int i=1;i<=r;i++)
            for(int j=1;j<=c;j++)
                a[i][j]=a[i][j]*x%MOD;
    }
    Matrix operator *(const Matrix &b) const{//重载乘法运算符
        Matrix ans(r,c);
        for(int i=1;i<=r;i++)
            for(int j=1;j<=b.c;j++)
                for(int k=1;k<=c;k++)
                    ans.a[i][j]=(ans.a[i][j]+a[i][k]*b.a[k][j])%MOD;
        return ans;
    }
    Matrix operator *=(const Matrix &b){//重载*=运算符
        return (*this)=(*this)*b;
    }
    Matrix power(int x){//定义矩阵幂运算函数
        Matrix M=(*this),ans(r,c,1);
        for(;x;x>>=1,M*=M) if(x&1) ans*=M;
        return ans;
    }
    void print(){//打印矩阵
        for(int i=1;i<=r;i++){
            for(int j=1;j<=c;j++) cout<<a[i][j]<<" ";
            cout<<endl;
        }
    }
};

六、参考文章

• 数列通项公式，但用的是矩阵（%%%@yggdyy_）

• 中国计算机学会《CCF中学生计算机程序设计》第10.2节（友情提示：这本书中的锅数不胜数，建议别全信这本书内的东西）

• 李煜东《算法竞赛进阶指南》第0x34节

• OI-Wiki - 快速幂