前言
矩阵优化是一种比较靠思维的优化算法，一般简单题考的比较少。
个人认为矩阵优化中在运用，所以放了几道题目来讲解。
定义
矩阵
一个 \(m\times n\) 的矩阵是一个由 \(m\) 行 \(n\) 列元素排列成的矩形阵列。大概长成下面这个样子的。
\[A = \underbrace{\begin{bmatrix} a_{1,1} & a_{1,2} & \cdots & a_{1,n} \\ a_{2,1} & a_{2,2} & \cdots & a_{2,n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m,1} & a_{m,2} & \cdots & a_{m,n} \end{bmatrix}}_{m}\Bigg\}n \]也就是说，存储以一个矩阵只需要存储它的行、列以及每个元素的值即可。
若两个矩阵行数和列数都相同，则我们管这两个矩阵叫做同型矩阵。
矩阵的运算
加法&减法
两个同型矩阵的加法运算如下
矩阵与数
将矩阵每个位置加/减这个数即可。
矩阵与矩阵
\[\begin{bmatrix} a_{1,1} & a_{1,2} & \cdots & a_{1,n} \\ a_{2,1} & a_{2,2} & \cdots & a_{2,n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m,1} & a_{m,2} & \cdots & a_{m,n} \end{bmatrix}+\begin{bmatrix} b_{1,1} & b_{1,2} & \cdots & b_{1,n} \\ b_{2,1} & b_{2,2} & \cdots & b_{2,n} \\ \vdots & \vdots & \ddots & \vdots \\ b_{m,1} & b_{m,2} & \cdots & b_{m,n} \end{bmatrix}=\begin{bmatrix} a_{1,1}+b_{1,1} & a_{1,2}+b_{1,2} & \cdots & a_{1,n}+b_{1,n} \\ a_{2,1}+b_{2,1} & a_{2,2}+b_{2,2} & \cdots & a_{2,n}+b_{2,n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m,1}+b_{m,1} & a_{m,2}+b_{m,2} & \cdots & a_{m,n}+b_{m,n} \end{bmatrix} \]上述矩阵为加法运算，减法同理。一般我们不考虑非同型矩阵的加减法运算。
乘法
矩阵与数
将矩阵每个元素都乘上那个数即可。
矩阵与矩阵
两个矩阵可以相乘必须要满足前一个矩阵的列数等于后一个矩阵的行数。
或者说，\(A\) 矩阵的大小为 \(n\times m\)，\(B\) 矩阵的大小为 \(m\times l\)，才能进行 \(A\times B\) 这个运算。
若 \(C=A\times B\)，则（此处矩阵大小同上文）
\[C_{i,j}=\sum_{k=1}^{m}{A_{i,k}\times B_{k,j}} \]不难发现，\(C\) 的大小为 \(n\times l\)。
矩阵乘法满足一下几条运算律：
\[(AB)C=A(BC)\\ A(B+C)=AB+BC\\ (B+C)A=BA+CA \]但请注意，矩阵乘法不满足交换律，即 \(AB\neq BA\)。
单位矩阵
在数的运算中，有 \(1\times a=a\)。在矩阵运算中，也有一个矩阵满足任意矩阵乘该矩阵等于原矩阵。若原矩阵为 \(n\times m\) 大小的矩阵 \(A\)，这个矩阵为 \(m\times m\) 大小的矩阵 \(B\)，则
\[B=\begin{bmatrix} 1 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \end{bmatrix} \]其左上到右下的对角线上的数都为 \(1\)，其余都为 \(0\)。我们管这样的矩阵叫单位矩阵。一个矩阵乘单位矩阵后完全不变。
运用
P3390 【模板】矩阵快速幂
将矩阵定义为结构体并重载运算符，直接计算即可。
#include<iostream>
#include<cstring>
using namespace std;
long long b,n,k,l=100,mod=1000000007;
struct Mat
{
#define N 110
long long mat[N][N];
};
Mat operator* (Mat a,Mat b)
{
Mat c;
memset(c.mat,0,sizeof(c.mat));
for(int k=1;k<=l;k++)
for(int i=1;i<=l;i++)
{
if(a.mat[i][k]<=0)
continue;
for(int j=1;j<=l;j++)
{
if(b.mat[k][j]<=0)
continue;
c.mat[i][j]=(c.mat[i][j]+a.mat[i][k]*b.mat[k][j])%mod;
}
}
return c;
}
Mat operator^ (Mat a,long long k)
{
Mat c;
// for(int i=1;i<=l;i++)
// for(int j=1;j<=l;j++)
// c.mat[i][j]=(i==j);
memset(c.mat,0,sizeof(c.mat));
for(int i=1;i<=n;i++)
c.mat[i][i]=1;
while(k)
{
if(k&1)
c=a*c;
a=a*a;
k>>=1;
}
return c;
}
int main()
{
cin>>n>>k;
Mat x;
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
cin>>x.mat[i][j];
x=x^k;
for(int i=1;i<=n;i++)
{
for(int j=1;j<=n;j++)
cout<<x.mat[i][j]<<" ";
cout<<endl;
}
}
P1939 【模板】矩阵加速（数列）
注意到 \(n\leq 2\times 10^9\)，不能直接用递推做。用矩阵快速幂快速计算。
简单的用矩阵优化的题都是把所有要用的东西都扔进一个矩阵里，比如这题：
\[W=\begin{bmatrix}f_i\\ f_{i-1}\\ f_{i-2}\end{bmatrix} \]其中 \(f_i\) 表示求得的数列 \(a\) 第 \(i\) 项的值。但要注意，由于这个矩阵是由 \(i-1\) 的矩阵推过来的，所以矩阵内能引用到的数分别是 \(a_{i-1},a_{i-2},a_{i-3}\)。
在递推式里，推导答案的式子有 \(m\) 项，则这个矩阵的大小为 \(m\times m\)。
我们需要一个矩阵 \(B\) 使得
\[B\times \begin{bmatrix}f_{i-1}\\ f_{i-2}\\ f_{i-3}\end{bmatrix}=\begin{bmatrix}f_i\\ f_{i-1}\\ f_{i-2}\end{bmatrix} \]在这题里，这个矩阵 \(B\) 长成下面这样：
\[\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} \]也就是说，这个矩阵与上述矩阵相乘的意义为
\[f_i=a_{1,1}f_{i-1}+a_{1,2}f_{i-2}+a_{1,3}f_{i-3}\\ f_{i-1}=a_{2,1}f_{i-1}+a_{2,2}f_{i-2}+a_{2,3}f_{i-3}\\ f_{i-2}=a_{3,1}f_{i-1}+a_{3,2}f_{i-2}+a_{3,3}f_{i-3}\\ \]我们可以比较容易的推出这个矩阵为
\[\begin{bmatrix}1 & 0 & 1\\ 1 & 0 & 0\\ 0 & 1 & 0 \end{bmatrix} \]如果要用矩阵的出 \(f_i\)，那么大概是这样的
\[\begin{bmatrix}f_3=1\\ f_2=1 \\ f_{1}=1\end{bmatrix}\times\underbrace{\begin{bmatrix}1 & 0 & 1\\ 1 & 0 & 0\\ 0 & 1 & 0 \end{bmatrix} \times \cdots \times \begin{bmatrix}1 & 0 & 1\\ 1 & 0 & 0\\ 0 & 1 & 0 \end{bmatrix} }_{i-3\text{个}} \]然后就可以直接用矩阵快速幂来做了。一共有 \(i-3\) 个矩阵 \(B\) 是因为 \(x\in\{1,2,3\}\) 的值为 \(1\)，注意答案为所得矩阵的第一项。
#include<iostream>
#include<cstring>
using namespace std;
#define int long long
#define mod 1000000007
struct Mat
{
#define N 4
int mat[N][N];
};
int t,n,l=3;
Mat operator* (Mat a,Mat b)
{
Mat c;
memset(c.mat,0,sizeof(c.mat));
for(int k=1;k<=3;k++)
for(int i=1;i<=3;i++)
for(int j=1;j<=3;j++)
{
if(a.mat[i][k]<0||b.mat[k][j]<0)
continue;
c.mat[i][j]=(c.mat[i][j]+a.mat[i][k]*b.mat[k][j])%mod;
}
return c;
}
Mat operator^ (Mat a,int k)
{
Mat c;
memset(c.mat,0,sizeof(c.mat));
for(int i=1;i<=3;i++)
c.mat[i][i]=1;
while(k)
{
if(k&1)
c=a*c;
a=a*a;
k>>=1;
}
return c;
}
signed main()
{
cin>>t;
while(t--)
{
cin>>n;
if(n<=3)
{
cout<<1<<endl;
continue;
}
Mat a,x;
memset(a.mat,0,sizeof(a.mat));
memset(x.mat,0,sizeof(x.mat));
a.mat[1][1]=a.mat[1][3]=a.mat[2][1]=a.mat[3][2]=1;
x.mat[1][1]=x.mat[1][2]=x.mat[1][3]=1;
x=x*(a^(n-3));
cout<<x.mat[1][1]<<endl;
}
}
P1962 斐波那契数列
思路和上题差不多，都是找到一个转移矩阵。
甚至比上题还要简单一些。
#include<iostream>
#include<cstring>
using namespace std;
unsigned long long b,n,k,l=2,p=1000000007;
struct Mat
{
unsigned long long mat[3][3];
};
Mat operator* (Mat a,Mat b)
{
Mat c;
memset(c.mat,0,sizeof(c.mat));
for(long long k=1;k<=l;k++)
for(long long i=1;i<=l;i++)
{
if(a.mat[i][k]<=0)
continue;
for(long long j=1;j<=l;j++)
{
if(b.mat[k][j]<=0)
continue;
c.mat[i][j]=(c.mat[i][j]+a.mat[i][k]*b.mat[k][j])%p;
}
}
return c;
}
Mat operator^ (Mat a,long long k)
{
Mat c;
// for(int i=1;i<=l;i++)
// for(int j=1;j<=l;j++)
// c.mat[i][j]=(i==j);
memset(c.mat,0,sizeof(c.mat));
for(long long i=1;i<=2;i++)
c.mat[i][i]=1;
while(k)
{
// cout<<">>"<<k<<endl;
if(k&1)
c=a*c;
a=a*a;
k>>=1;
}
return c;
}
int main()
{
cin>>n;
Mat x,y;
memset(x.mat,0,sizeof(x.mat));
memset(y.mat,0,sizeof(y.mat));
x.mat[1][1]=1;
x.mat[2][1]=1;
y.mat[1][1]=y.mat[1][2]=1;
y.mat[2][1]=1;
y=y^(n-1);
x=x*y;
cout<<x.mat[1][1]<<endl;
}