# **1 矩阵**
## **1.1 介绍**
何为**矩阵**？
矩阵是一个按照长方阵列排列的复数或实数集合。它长这样：
$$
\begin{bmatrix}
2&3&4\\
1&6&5
\end{bmatrix}
\begin{bmatrix}
3&1&3\\
3&2&3\\
3&1&3
\end{bmatrix}
\begin{bmatrix}
1&0&0&0\\
0&1&0&0\\
0&0&1&0\\
0&0&0&1
\end{bmatrix}
$$
由m×n个数$a_{ij}$排成的m行n列的数表称为m行n列的矩阵，简称**m×n矩阵**。可以看出，以上矩阵分别是2x3矩阵、3x3矩阵和4x4矩阵。行数与列数都等于n的矩阵称为**n阶矩阵**或**n阶方阵**，因此后两个矩阵也可称为3阶方阵和4阶方阵。最后一个矩阵还可称为**单位矩阵**，因为其正对角线上的数字为1，其余数字为0。
## **1.2 矩阵的基本运算**
**矩阵运算**在科学计算中非常重要，而矩阵的基本运算包括矩阵的加法，减法，数乘，~~转置，共轭和共轭转置~~。
- ### **加法**
**矩阵加法**是矩阵对应位置上的数字相加，这要求了相加的矩阵必须是同型矩阵（同样的m,n）。
$$
\begin{bmatrix}
1&4&2\\
2&0&0
\end{bmatrix}
+
\begin{bmatrix}
0&0&5\\
7&5&0
\end{bmatrix}
=
\begin{bmatrix}
1+0&4+0&2+5\\
2+7&0+5&0+0
\end{bmatrix}
=
\begin{bmatrix}
1&4&7\\
9&5&0
\end{bmatrix}
$$
矩阵加法满足以下运算律：
$$
A+B=B+A\\
(A+B)+C=A+(B+C)
$$
- ### **矩阵减法**
**矩阵减法**是矩阵对应位置上的数字相减，这也要求了相减的矩阵必须是同型矩阵（同样的m,n）。
$$
\begin{bmatrix}
1&4&2\\
2&0&0
\end{bmatrix}
+
\begin{bmatrix}
0&0&5\\
7&5&0
\end{bmatrix}
=
\begin{bmatrix}
1-0&4-0&2-5\\
2-7&0-5&0-0
\end{bmatrix}
=
\begin{bmatrix}
1&4&-3\\
-5&-5&0
\end{bmatrix}
$$
- ### **数乘**
**矩阵数乘**就是用一个常数乘以矩阵所有位置上的数。
$$
2\cdot \begin{bmatrix}
1&8&-3\\
4&-2&5
\end{bmatrix}
=
\begin{bmatrix}
2\cdot 1&2\cdot 8&2\cdot (-3)\\
2\cdot 4&2\cdot (-2)&2\cdot 5
\end{bmatrix}
=
\begin{bmatrix}
2&16&-6\\
8&-4&10
\end{bmatrix}
$$
矩阵数乘满足以下运算律：
$$
\lambda (\mu A)=\mu (\lambda A)\\
\lambda (\mu A)=(\lambda \mu)A\\
(\lambda +\mu)A=\lambda A+\mu A\\
\lambda(A+B)=\lambda A+\lambda B
$$
## **1.3 矩阵乘法**
两个矩阵的**乘法**仅当第一个矩阵A的列数和另一个矩阵B的行数相等时才能定义。如A是m×n矩阵和B是n×p矩阵，它们的乘积C是一个m×p矩阵$C=(c_{ij})$。它的一个元素:
$$
c_{ij}=\sum_{r=1}^{n}a_{ir}b_{rj}
$$
将这矩阵乘法的式子记为:**C=AB**。
$$
\begin{bmatrix}
1&0&2\\
-1&3&1
\end{bmatrix}
\times
\begin{bmatrix}
3&1\\
2&1\\
1&0
\end{bmatrix}
=
\begin{bmatrix}
(1\times 3+0\times 2+2\times 1)&(1\times 1+0\times 1+2\times 0)\\
(-1\times 3+3\times 2+1\times 1)&(-1\times 1+3\times 1\times 0)
\end{bmatrix}
=
\begin{bmatrix}
5&1\\
4&2
\end{bmatrix}
$$
矩阵的乘法满足以下运算律：
$$
(AB)C=A(BC)\\
(A+B)C=AC+BC\\
C(A+B)=CA+CB
$$
**矩阵乘法不满足交换律**。
# **2 矩阵快速幂**
## **2.1 简介与原理**
**矩阵快速幂**是用来高效地计算矩阵的高次方幂的一种算法，它将计算的时间复杂度降到了$O(logn)$。
一般来说，求一个矩阵的n次方幂，我们会通过连乘n-1次来得到它的n次方幂。但做下简单的改进就能减少连乘的次数，方法如下：
- 把n个矩阵进行两两分组，例如$A\cdot A\cdot A\cdot A\cdot A\cdot A=(A\cdot A)\cdot (A\cdot A)\cdot (A\cdot A)$。
这样变的好处是，你只需要计算一次$A\cdot A$，然后将结果$A\cdot A$连乘2次就能得到$A^6$，即$(A\cdot A)^3=A^6$。这其中只做了3次运算，少于原来的5次。
其实还可以取$A^3$作为一个基本单位。原理都是一样，利用矩阵乘法的结合律，来减少重复计算的次数。
以上都是取一个具体的数来作为最小单位的长度，这样做虽然能够改进效率，但缺陷也是很明显的，取个极端的例子，当n无限趋近于+∞的时候，你现在所取的长度其实和1没什么区别。既然要减少重复计算，那么就要充分利用现有的计算结果。如何充分利用计算结果呢？这里考虑二分的思想。
例如这个式子：$A^{19}=>A^{16}\cdot A^2\cdot A$，显然采取这样的方式计算时复杂度将是$O(logn)$级别的，因为这些分解的指数都为2的次幂。可以知道，如果要用尽可能少的数来表示一个范围的所有数字，那么用2的次幂显然是最好的，所以这样就只需要计算这些2的次幂就可以了。而1-n中的2的次幂的数量为logn个。
这同样也是**快速幂**的原理。既然原理相同，那么就只需要在快速幂的代码的基础上修改就可以了。下面是快速幂的代码。
```cpp
long long fast_power(int a,int b,int mod)
{
long long ans=1,base=a;
while(b)
{
if(b&1)
ans=ans*base%mod;
base=base*base%mod;
b>>=1;
}
return ans%mod; //非多余步，如果没有b=0这种数据可以去掉取余
}
```
在快速幂的代码的基础上把底数换为矩阵，将四则乘法运算转为矩阵乘法运算并且将取模运算内置在矩阵乘法运算中可以使其成为矩阵快速幂了。要注意的是，只有**n阶方阵**才会有次幂。
## **2.2 矩阵与矩阵乘法的实现**
首先来实现矩阵与矩阵乘法。矩阵在代码中直接用二维数组表示即可。矩阵乘法这个也只需要根据定义的式子$c_{ij}=\sum_{r=1}^{n}a_{ir}b_{rj}$来写就可以了，没有什么简便方法。这里使用结构体和重载运算符的方式来表示矩阵和矩阵乘法。请看代码。
```cpp
int n; //矩阵阶数
struct Matrix //矩阵结构体
{
#define mod int(1e9+7)
long long m[105][105];
void init() //初始化
{
memset(m,0,sizeof(m));
return;
}
Matrix operator*(const Matrix b) const //矩阵乘法
{
Matrix t;
t.init();
for (int i=1;i<=n;i++)//基本枚举，下同
for (int j=1;j<=n;j++)
for (int k=1;k<=n;k++)
t.m[i][j]=(t.m[i][j]+(m[i][k]*b.m[k][j])%mod)%mod; //取模运算不能丢
return t;
}
Matrix& operator*=(const Matrix b)
{
return *this=*this*b;
}
//输入输出
void in()
{
for (int i=1;i<=n;i++)
for (int j=1;j<=n;j++)
cin>>m[i][j];
return;
}
void out()
{
for (int i=1;i<=n;i++,cout<<'\n')
for (int j=1;j<=n;j++)
cout<<m[i][j]<<' ';
return;
}
};
```
## **2.3 最终整体**
写完矩阵和矩阵乘法之后，矩阵快速幂的函数就非常好写了，只不过这里还得讨论一下ans=1怎么改。可知1在计算中可以起到非0的数与之相乘还是该数，同理矩阵中也有一个矩阵有类似的作用，它就是上面提到的单位矩阵。所以ans=1的部分用单位矩阵替代即可。请看代码。
```cpp
Matrix fast_matrix_power(Matrix a,long long b)
{
Matrix ans,base=a;
ans.init();
for (int i=1;i<=n;i++) //构建单位矩阵
ans.m[i][i]=1;
while(b)
{
if(b&1)
ans*=base;
base*=base;
b>>=1;
}
return ans;
}
```
这样将这些片段合并后就可以得到矩阵快速幂的完整代码。请看。
```cpp
int n; //矩阵阶数
long long k;
struct Matrix //矩阵结构体
{
#define mod int(1e9+7)
long long m[105][105];
void init() //初始化
{
memset(m,0,sizeof(m));
return;
}
Matrix operator*(const Matrix b) const //矩阵乘法
{
Matrix t;
t.init();
for (int i=1;i<=n;i++)//基本枚举，下同
for (int j=1;j<=n;j++)
for (int k=1;k<=n;k++)
t.m[i][j]=(t.m[i][j]+(m[i][k]*b.m[k][j])%mod)%mod;
return t;
}
Matrix& operator*=(const Matrix b)
{
return *this=*this*b;
}
//输入输出
void in()
{
for (int i=1;i<=n;i++)
for (int j=1;j<=n;j++)
cin>>m[i][j];
return;
}
void out()
{
for (int i=1;i<=n;i++,cout<<'\n')
for (int j=1;j<=n;j++)
cout<<m[i][j]<<' ';
return;
}
};
Matrix fast_matrix_power(Matrix a,long long b)
{
Matrix ans,base=a;
ans.init();
for (int i=1;i<=n;i++) //构建单位矩阵
ans.m[i][i]=1;
while(b)
{
if(b&1)
ans*=base;
base*=base;
b>>=1;
}
return ans;
}
int main()
{
Matrix a;
cin>>n>>k;
a.in();
fast_matrix_power(a,k).out();
return 0;
}
```
# **3 矩阵加速**
## **3.1 线性方程组**
线性方程组是各个方程关系未知量均为一次的方程组。它长下面的样子：
$$
\begin{cases}
a_{11}x_1+a_{12}x_2+a_{13}x_3=b_1\\
a_{21}x_1+a_{22}x_2+a_{23}x_3=b_2\\
a_{31}x_1+a_{32}x_2+a_{33}x_3=b_3
\end{cases}
$$
这是一个三元线性方程组，其中$x_i$表未知量，$a_{ij}$称为系数，$b_i$称为常数项。将系数用矩阵来表示就是下面的系数矩阵：
$$
\begin{bmatrix}
a_{11}&a_{12}&a_{13}\\
a_{21}&a_{22}&a_{23}\\
a_{31}&a_{32}&a_{33}
\end{bmatrix}
$$
由矩阵乘法可以得到关系：
$$
\begin{bmatrix}
a_{11}&a_{12}&a_{13}\\
a_{21}&a_{22}&a_{23}\\
a_{31}&a_{32}&a_{33}
\end{bmatrix}
\begin{bmatrix}
x_1\\
x_2\\
x_3
\end{bmatrix}
=
\begin{bmatrix}
b_1\\
b_2\\
b_3
\end{bmatrix}
$$
这是**矩阵加速**的基本原理。
## **3.2 矩阵加速**
之前学习的求斐波那契数列的方法是递归或递推，虽然非常简单，但是却不能快速计算斐波那契数列中下标比较大的项。现在，我们可以通过构造常系数矩阵再加上矩阵快速幂来快速计算下标比较大的项。这一方法称为被**矩阵加速**。
对于斐波那契数列，将f[n-1]和f[n-2]作为未知量，递推式也同样可以表示为线性方程组：
$$
\begin{cases}
1\cdot f[n-1]+1\cdot f[n-2]=f[n]\\
1\cdot f[n-1]+0\cdot f[n-2]=f[n-1]
\end{cases}
$$
那么其矩阵乘法的关系就是：
$$
\begin{bmatrix}
1&1\\
1&0
\end{bmatrix}
\begin{bmatrix}
f[n-1]\\
f[n-2]
\end{bmatrix}
=
\begin{bmatrix}
f[n]\\
f[n-1]
\end{bmatrix}
$$
你可能要问到：为什么是两个方程而不是一个，一个方程不也够了？这是因为递推式f[n]=f[n-1]+f[n-2]决定的矩阵乘法关系是：
$$
\begin{bmatrix}
1&1
\end{bmatrix}
\begin{bmatrix}
f[n-1]\\
f[n-2]
\end{bmatrix}
=
\begin{bmatrix}
f[n]
\end{bmatrix}
$$
将系数矩阵乘以由未知量构成的矩阵看为一个特殊的函数对应法则g。那么就是这样：
$$
g\left(
\begin{bmatrix}
f[n-1]\\
f[n-2]
\end{bmatrix}
\right) =
\begin{bmatrix}
f[n]
\end{bmatrix}
$$
你可以从左推导到右，但是你不能从右推导到左。或者说你不能求出使得输入是右，输出是左的“系数矩阵”，这个“系数矩阵”也称为**逆矩阵**。有逆矩阵的一个充分条件是原矩阵是n阶方阵。自此，方程为什么是两个的问题解决了。这也延伸出一个思路，有多少个未知量就代表需要多少个方程。
回到正题，得到矩阵乘法关系后，可以知道只要知道初始条件，那么就可以通过不断将系数矩阵乘以由未知量构成的矩阵来得到最重需要的答案，又因为矩阵乘法的结合律，只需要先求出系数矩阵的次幂然后再乘以由初始条件构成的矩阵就可以得到**结果矩阵**（这里也同样限制了系数矩阵是n阶方阵）。
## **3.3 矩阵构造法**
矩阵加速的核心就是系数矩阵的构造，所以这里谈一谈一些系数矩阵构造的例子。不必考虑未知量的顺序，对了即可。
- ### **递推式f[n]=f[n-1]+f[n-2]，初始条件f[1]=f[2]=1。求第n项**
这个已经在上面提到过。不再过多赘述
- ### **递推式f[n]=f[n-1]+f[n-2]+1，初始条件f[1]=f[2]=1。求第n项**
注意1也要算作为未知量，但是因为它是常数，所以需要保证其值不变。未知量矩阵，系数矩阵和结果矩阵为：
$$
\begin{bmatrix}
f[n-1]\\
f[n-2]\\
1
\end{bmatrix}
,
\begin{bmatrix}
1&1&1\\
1&0&0\\
0&0&1
\end{bmatrix}
,
\begin{bmatrix}
f[n]\\
f[n-1]\\
1
\end{bmatrix}
$$
- ### **递推式f[n]=f[n-1]+f[n-2]+n+1，初始条件f[1]=f[2]=1。求第n项**
同理。未知量矩阵，系数矩阵和结果矩阵为：
$$
\begin{bmatrix}
f[n-1]\\
f[n-2]\\
n\\
1
\end{bmatrix}
,
\begin{bmatrix}
1&1&1&1\\
1&0&0&0\\
0&0&1&1\\
0&0&0&1
\end{bmatrix}
,
\begin{bmatrix}
f[n]\\
f[n-1]\\
n+1\\
1
\end{bmatrix}
$$
当然可以适量的降低一下复杂度，把n+1看成一个未知量，那么未知量矩阵，系数矩阵和结果矩阵为：
$$
\begin{bmatrix}
f[n-1]\\
f[n-2]\\
n+1
\end{bmatrix}
,
\begin{bmatrix}
1&1&1\\
1&0&0\\
0&0&1
\end{bmatrix}
,
\begin{bmatrix}
f[n]\\
f[n-1]\\
n+2
\end{bmatrix}
$$
- ### **递推式f[n]=f[n-1]+f[n-2]，初始条件f[1]=f[2]=1。求前n项和s[n]**
由于s[n]是我们要求的答案，所以s[n]也当成未知量，可以获得一个额外的递推式s[n-1]=s[n-2]+f[n-1]。那么未知量矩阵，系数矩阵和结果矩阵为：
$$
\begin{bmatrix}
s[n-2]\\
f[n-1]\\
f[n-2]
\end{bmatrix}
,
\begin{bmatrix}
1&1&0\\
0&1&1\\
0&1&0
\end{bmatrix}
,
\begin{bmatrix}
s[n-1]\\
f[n]\\
f[n-1]
\end{bmatrix}
$$
---
感谢浏览😝！
此文章可能会在后续更新，欢迎纠错。
S1-矩阵快速幂与矩阵加速
