Files
python/TangDou/AcWing_TiGao/T5/JuZhengChengFa/矩阵乘法与矩阵快速幂.md

428 lines
15 KiB
Markdown
Raw Normal View History

2025-08-30 18:35:01 +08:00
## 矩阵乘法与矩阵快速幂
### 一、矩阵的定义
一个 $m\times n$ 的矩阵就是 $m$ 行 $n$ 列的数字阵列,如$2 \times 3$的矩阵:$\begin{bmatrix}
1 & 5 & 3 \\
2 & 8 & 10
\end{bmatrix}$
实际上,**矩阵类似二维数组**。
### 2. 矩阵的运算
- **矩阵的加法和减法**
矩阵的加法和减法就是将两个矩阵对应位置上的数相加减。因此,相加减的两个矩阵$A$,$B$的行列必须相同。
- **矩阵乘法**
$A$,$B$,$C$ 是三个矩阵,若 $A\times B=C$,需要满足:
- $A$ 的列数必须和 $B$ 的行数相等;
- 设 $A$ 是一个 $n\times r$ 的矩阵,$B$ 是一个 $r\times m$ 的矩阵,则矩阵 $A$ 乘 矩阵 $B$ 的乘积 $C$ 是一个 $n\times m$ 的矩阵;
- $C_{(i,j)} =A_{(i,1)} ×B_{(1,j)} +A_{(i,2)}×B_{(2,j)} +...A_{(r,1)} ×B_{(r,j)}$
​
矩阵 $C$ 的第 $i$ 行第 $j$ 列元素 = 矩阵 $A$ 的第 $i$ 行元素与矩阵 $B$的第 $j$ 列对应元素乘积之和。
例如:$A= \begin{bmatrix}
a & b \\
c & d \\
e & f
\end{bmatrix},
B= \begin{bmatrix}
g & h & i \\
j & k & l
\end{bmatrix}$,求$C=A \times B$
$A \times B=\begin{bmatrix}
a \times g+b \times j & a \times h+b \times k & a \times i + b \times l \\
c \times g+d \times j & c \times h+d \times k & c \times i + d \times l \\
e \times g+f \times j & e \times h+f \times k & e \times i + f \times l
\end{bmatrix}$
可以发现 $A\times B$ 和 $B\times A$ 将得到两种不同的结果,因此矩阵不满足乘法交换律。
- **方阵次幂**
如果矩阵的行和列相同,称矩阵为方阵。若 $A$ 是一个方阵,方阵的幂是指,将 $A$ 连乘 $n$ 次,$A^n$。
若不是方阵则无法进行乘幂运算。
**矩阵乘法不满足交换律,但是满足结合律,因此可以用快速幂的思想来求解方程次幂**。
### 3. 矩阵乘法的应用
在信息学竞赛中考矩阵乘法不是为了考基本运算,而是利用矩阵乘法的特性,配合快速幂求解递推式的第 $n$ 项,以及一些图论的构造转化。
**[一本通 $1641$ 矩阵乘法](http://ybt.ssoier.cn:8088/problem_show.php?pid=1641)**
**[$B2105$ 矩阵乘法](https://www.luogu.com.cn/problem/B2105)**
> **注释**:上面两道题其实是一样的模板,只不过输入的数据顺序不太一样。
```cpp {.line-numbers}
#include <bits/stdc++.h>
using namespace std;
const int N = 110;
int n, m, k;
int a[N][N], b[N][N], c[N][N];
// 矩阵乘法模板
void mul(int a[][N], int b[][N], int c[][N]) {
// 注意:这里的临时数组t绝不是画蛇添足!
// 因为调用的时候,有时会传递mul(a,b,b)这样的东东,如果每次memset(c,0,sizeof c),
// 就会造成b在运行前被清空,导致结果错误
// 代码解释:
// 从结果出发理解:
// C(i,j) =A(i,1)×B(1,j)+A(i,2)×B(2,j) +...A(r,1)×B(r,j)
// 抽象一下:
// C(i,j) = C(i,j) + A(i,k)× B(k,j)
int t[N][N] = {0};
for (int i = 0; i < N; i++)
for (int j = 0; j < N; j++)
for (int k = 0; k < N; k++)
t[i][j] = t[i][j] + a[i][k] * b[k][j]; // 矩阵乘法
memcpy(c, t, sizeof t);
}
int main() {
cin >> n >> m >> k;
// A矩阵 n*m
for (int i = 1; i <= n; i++)
for (int j = 1; j <= m; j++)
cin >> a[i][j];
// B矩阵m*k
for (int i = 1; i <= m; i++)
for (int j = 1; j <= k; j++)
cin >> b[i][j];
// 矩阵乘法
mul(a, b, c);
// 输出结果,控制格式
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= k; j++)
printf("%d ", c[i][j]);
printf("\n");
}
return 0;
}
```
**[$P3390$ 【模板】矩阵快速幂](https://www.luogu.com.cn/problem/P3390)**
- **快速幂**
实现 $pow(x, n)$ ,即计算 $x$ 的整数 $n$ 次幂函数 (即$x^n$)。
对于上面这道题,我们最容易想到的解法就是通过循环或递归进行累乘求解
```cpp {.line-numbers}
public double MyPow(double x, int n){
long N = n;
if (N == 0) return 1;
double res = x;
long num = N > 0 ? N : -N;
for (long i = 1; i < num; i++){
res *= x;
}
if (N < 0) res = 1 / res;
return res;
}
```
然而,这种解法在面对$n$比较大的情况时,效率会变得十分低下。仔细想想上面的循环过程就会发现,里面进行了很多次无意义的计算。以$a^{10}$ 举例,通过循环一个一个往上乘,需要进行$9$次运算。而如果我们先计算出$a^2$ ($1$次),然后计算$a^2*a^2*a^2*a^2*a^2$ ($4$次),总计算量就能减少到$5$次。如果先计算出$a^4$($2$次),再计算$a^2*a^4*a^4$($2$次),总计算量就能减少到$4$次。
也就是说,将$a^n$ 拆解成$a^x*a^y*a^z...$的形式就可以降低总体运算的次数。那么如何拆解$n$才最高效呢?根据前面的思路,显然是每次平方所需的运算次数最少。仍然是$a^{10}$ 的例子,我们从$a^1$
开始,每次取平方值,可以得到$[a^1,a^2,a^4,a^8]$。想要组成$a^{10}$ ,只需要让$a^8*a^2$ 即可。所需的计算次数也是$4$次。看到这里有没有觉得熟悉?我们实际上就是在把$n$拆成几个$2^?$相加,也就是类似于$n$的二进制转换成十进制的过程。
下面仍然以$a^{10}$ 为例,通过图展示一下计算过程:
我们从$10$的最低位开始遍历,遇到的第一位是$0$,对应的是$a^1$ ,但我们并不需要它计算结果,所以跳过。
![](https://dsideal.obs.cn-north-1.myhuaweicloud.com/HuangHai/BlogImages/202311270904460.png)
接下来遇到的第二位是$1$,对应$a^2$ ,这是计算结果需要用到的,所以把它乘到$res$上。继续移动指针。
![](https://dsideal.obs.cn-north-1.myhuaweicloud.com/HuangHai/BlogImages/202311270904117.png)
遇到的第三位是$0$,对应$a^4$ ,用不到,跳过。
![](https://dsideal.obs.cn-north-1.myhuaweicloud.com/HuangHai/BlogImages/202311270905475.png)
遇到的第四位是$1$,对应$a^8$ ,需要用到,把它乘到$res$中。
![](https://dsideal.obs.cn-north-1.myhuaweicloud.com/HuangHai/BlogImages/202311270905174.png)
$res$便是最终结果。将上面的思路转换成代码(未考虑特殊情况):
```cpp {.line-numbers}
double pow(double a, int n){
double res = 1;
while (n){
// 当前位为1,计入结果
if ((n & 1) == 1) res = res * a;
// a自乘
a *= a;
// n右移
n >>= 1;
}
return res;
}
```
- **单位矩阵**
在普通的乘法中,一个数乘$1$还是等于它本身,在矩阵乘法中也有这么一个$1$,它就是 **单位矩阵**
不同于普通乘法中的单位$1$,对于不同矩阵他们的单位矩阵大小是不同的
对于$n∗m$的矩阵,它的单位矩阵大小为$m∗m$,对于$m∗n$的矩阵,它的单位矩阵大小为$n∗n$
也就是说单位矩阵都是正方形的,这是因为只有正方形的矩阵能保证结果和前一个矩阵形状相同
单位矩阵的元素非$0$即$1$,从左上角到右下角的对角线上元素皆为$1$,其他皆为$0$
**矩阵快速幂-(单位矩阵)**
```cpp {.line-numbers}
#include <bits/stdc++.h>
using namespace std;
#define int long long
#define endl "\n"
const int MOD = 1e9 + 7;
const int N = 110;
int n, k; // 本题数据范围很大,用int直接wa哭了
int a[N][N], b[N][N]; // 原始矩阵,结果矩阵
// 矩阵乘法模板
void mul(int a[][N], int b[][N], int c[][N]) {
int t[N][N] = {0};
for (int i = 0; i < N; i++)
for (int j = 0; j < N; j++)
for (int k = 0; k < N; k++)
t[i][j] = (t[i][j] + a[i][k] * b[k][j] % MOD) % MOD; // 矩阵乘法
memcpy(c, t, sizeof t);
}
signed main() {
cin >> n >> k;
// 输入原始矩阵
for (int i = 0; i < n; i++)
for (int j = 0; j < n; j++)
cin >> a[i][j];
// 1、单位矩阵
for (int i = 0; i < n; i++) b[i][i] = 1;
// 2、矩阵快速幂
for (int i = k; i; i >>= 1) {
if (i & 1) mul(a, b, b);
mul(a, a, a);
}
// 3、输出
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++)
printf("%lld ", b[i][j]);
printf("\n");
}
}
```
**矩阵快速幂-(不使用单位矩阵)**
```cpp {.line-numbers}
#include <bits/stdc++.h>
using namespace std;
#define int long long
#define endl "\n"
const int MOD = 1e9 + 7;
const int N = 110;
int n, k; // 本题数据范围很大,用int直接wa哭了
int a[N][N], b[N][N]; // 原始矩阵,结果矩阵
// 矩阵乘法
void mul(int a[][N], int b[][N], int c[][N]) {
int t[N][N] = {0};
for (int i = 0; i < N; i++)
for (int j = 0; j < N; j++)
for (int k = 0; k < N; k++)
t[i][j] = (t[i][j] + a[i][k] * b[k][j] % MOD) % MOD; // 矩阵乘法
memcpy(c, t, sizeof t);
}
int main() {
cin >> n >> k;
// 1、输入原始矩阵
for (int i = 0; i < n; i++)
for (int j = 0; j < n; j++) {
cin >> a[i][j];
b[i][j] = a[i][j]; // 这种解法与标准解法不同,不用构建单位矩阵
// 直接赋值初始化了一个,然后执行k-1次就完成了矩阵的k次幂
}
// 2、计算矩阵快速幂
for (int i = k - 1; i; i >>= 1) {
if (i & 1) mul(a, b, b);
mul(a, a, a);
}
// 3、输出
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++)
printf("%lld ", b[i][j]);
printf("\n");
}
return 0;
}
```
> **理解**:不管是用单位矩阵也好,不用单位矩阵也罢,本质上都是类似于构造一个乘法的基底:数字$1$!
**矩阵快速幂应用:斐波那契数列**
矩阵快速幂的一个非常典型的应用就是斐波那契数列(**[原题地址](https://leetcode.cn/problems/fei-bo-na-qi-shu-lie-lcof/)**)。先来看斐波那契数列的定义:
$$F(1)=0,F(2)=1$$
$$F(n)=F(n-1)+F(n-2),其中n>1$$
我们先定义一个矩阵
$$b=\begin{bmatrix} F(n) & F(n-1) \end{bmatrix}$$
现在假设有一个变换矩阵$M$,与上述矩阵存在如下关系
$$\begin{bmatrix}F(n) & F(n-1) \end{bmatrix} \times M= \begin{bmatrix}F(n+1) & F(n) \end{bmatrix}$$
很显然,$M$是一个$2×2$的矩阵,我们假设
$$M=\begin{bmatrix}
a & b \\
c & d
\end{bmatrix}$$
带入上面的公式可得
$\begin{bmatrix}F(n)& F(n-1)\end{bmatrix} \times
\begin{bmatrix}a& b \\
c & d \\
\end{bmatrix}=\begin{bmatrix}
aF(N)+cF(n-1)& bF(n)+dF(n-1)
\end{bmatrix} =
\begin{bmatrix} F(n+1)& F(n) \end{bmatrix}$
也就是
$$aF(n)+cF(n-1)=F(n+1)$$
$$bF(n)+dF(n-1)=F(n)$$
依照$Fibonacci$的现实含义,很容易得出
$$M=\begin{bmatrix}
1 & 1 \\
1 & 0
\end{bmatrix}$$
即
$$\begin{bmatrix} F(n) & F(n-1) \end{bmatrix} \times
\begin{bmatrix} 1 & 1 \\
1 & 0
\end{bmatrix}=\begin{bmatrix} F(n+1) & F(n) \end{bmatrix} \ ①$$
同理
$$\begin{bmatrix} F(n-1) & F(n-2) \end{bmatrix} \times
\begin{bmatrix} 1 & 1 \\
1 & 0
\end{bmatrix}=\begin{bmatrix} F(n) & F(n-1) \end{bmatrix} \ ②$$
将 ② 代入 ①,得
$$\begin{bmatrix}F(n+1) & F(n) \end{bmatrix} =
\begin{bmatrix}F(n-1) & F(n-2) \end{bmatrix} \times \begin{bmatrix}
1 & 1 \\
1 & 0
\end{bmatrix} ^{2}$$
根据数学归纳法:
$\begin{bmatrix} F(n-1)& F(n-2) \end{bmatrix}$ 与$M$矩阵的平方匹配.也就是$F(n-1)$中$n-1+2=n+1$,一路递推上去,就是
$$\begin{bmatrix}F(n+1) & F(n) \end{bmatrix} =
\begin{bmatrix}F(2) & F(1) \end{bmatrix} \times \begin{bmatrix}
1 & 1 \\
1 & 0
\end{bmatrix} ^{n-1}$$
我们只需要定义出矩阵乘法,配合前面的快速幂,即可实现在$O(logN)$时间复杂度内解出斐波那契数列。
**[$YBT$_$1642$ $Fibonacci$ 第 $n$ 项](https://www.cnblogs.com/littlehb/p/17074586.html)**
```cpp {.line-numbers}
#include <bits/stdc++.h>
using namespace std;
const int INF = 0x3f3f3f3f;
#define int long long
#define endl "\n"
/*
http://ybt.ssoier.cn:8088/problem_show.php?pid=1642
测试用例:
5 1000
答案:
5
*/
const int N = 2;
int n, mod;
// 矩阵乘法
void mul(int a[][N], int b[][N], int c[][N]) {
int t[N][N] = {0};
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++)
for (int k = 0; k < N; k++)
t[i][j] = (t[i][j] + (a[i][k] * b[k][j]) % mod) % mod;
}
memcpy(c, t, sizeof t);
}
signed main() {
cin >> n >> mod;
int b[N][N] = {1, 1}; // 结果矩阵(初始化fib[2]=1,fib[1]=1)
// 构造的向量数组
int m[N][N] = {
{1, 1},
{1, 0}};
// 矩阵快速幂
for (int i = n - 1; i; i >>= 1) {
if (i & 1) mul(b, m, b);
mul(m, m, m);
}
printf("%lld\n", b[0][1]);
return 0;
}
```
> **$Q_1$**:为什么在计算矩阵快速幂时提出了一个单元矩阵的概念,而在上面的斐波那契第$n$项中,却没有使用单元矩阵呢?什么时候用单元矩阵,什么时候不用呢?
> **答**:不管是单元矩阵组装出的$1$的概念基底,还是默认把$1$次方做为基底,都是希望有一个**基底**,斐波那契数列这道题中,由于有$b[2],b[1]$组成的初始值,构成了一个 **基底** 矩阵,所以,就可以直接使用这引基底矩阵一路乘上去就行了。
>**$Q_2$: $b[0],b[1]$明明是一维数组,或者叫只有一行的矩阵,为什么还要声明为 $b[N][N]$呢?**
> **答**:是因为我懒!因为我懒得要写两个矩阵乘法的模板,也懒得背两个矩阵乘法的模板!是哪两个呢?
> - 二维乘一维
> - 二维乘二维
> 我整理的模板就是二维乘二维的,没有整理二维乘一维,或者,一维乘二维的。因为二维乘二维可以模拟出二维乘一维的:对于一维矩阵,只需要初始化第一行的所有数据,其它数据为$0$,就可以模拟出二维乘一维的效果了。
>**$Q_3$:$b$数组的含义是什么呢?**
**答**:其实$b$矩阵是构造出来的,表现形式是$b=[F(n) \ F(n-1)]$,按题目的意思,$Fibonacci$数列是从下标$1$开始的,$f_1=f_2=1$,所以$b$的初始值就是$[F(2) \ F(1)]$
**[$AcWing$ $1303$. 斐波那契前 $n$ 项和](https://www.cnblogs.com/littlehb/p/16341964.html)**
相比于上一道题,本题增加了一个概念$S_n$,含义是前$n$个$f_n$的和,难度有所增加。
**[$P1939$ 矩阵加速(数列)](https://www.cnblogs.com/littlehb/p/16039658.html)**
**[$AcWing$ $1304$. 佳佳的斐波那契](https://www.cnblogs.com/littlehb/p/16342848.html)**
**[$AcWing$ $1305$. $GT$考试](https://www.cnblogs.com/littlehb/p/16348443.html)**
https://blog.csdn.net/STRVE/article/details/106739349