矩阵优化学习笔记
前言
矩阵优化是一种比较靠思维的优化算法，一般简单题考的比较少。
个人认为矩阵优化中在运用，所以放了几道题目来讲解。
定义
矩阵
一个 \(m\times n\) 的矩阵是一个由 \(m\) 行 \(n\) 列元素排列成的矩形阵列。大概长成下面这个样子的。
也就是说，存储以一个矩阵只需要存储它的行、列以及每个元素的值即可。
若两个矩阵行数和列数都相同，则我们管这两个矩阵叫做同型矩阵。
矩阵的运算
加法&减法
两个同型矩阵的加法运算如下
矩阵与数
将矩阵每个位置加/减这个数即可。
矩阵与矩阵
上述矩阵为加法运算，减法同理。一般我们不考虑非同型矩阵的加减法运算。
乘法
矩阵与数
将矩阵每个元素都乘上那个数即可。
矩阵与矩阵
两个矩阵可以相乘必须要满足前一个矩阵的列数等于后一个矩阵的行数。
或者说，\(A\) 矩阵的大小为 \(n\times m\)，\(B\) 矩阵的大小为 \(m\times l\)，才能进行 \(A\times B\) 这个运算。
若 \(C=A\times B\)，则（此处矩阵大小同上文）
不难发现，\(C\) 的大小为 \(n\times l\)。
矩阵乘法满足一下几条运算律：
但请注意，矩阵乘法不满足交换律，即 \(AB\neq BA\)。
单位矩阵
在数的运算中，有 \(1\times a=a\)。在矩阵运算中，也有一个矩阵满足任意矩阵乘该矩阵等于原矩阵。若原矩阵为 \(n\times m\) 大小的矩阵 \(A\)，这个矩阵为 \(m\times m\) 大小的矩阵 \(B\)，则
其左上到右下的对角线上的数都为 \(1\)，其余都为 \(0\)。我们管这样的矩阵叫单位矩阵。一个矩阵乘单位矩阵后完全不变。
运用
P3390 【模板】矩阵快速幂
因为矩阵乘法满足 \(A(BC)=(AB)C\)，所以我们可以将快速幂套到矩阵上。
将矩阵定义为结构体并重载运算符，直接计算即可。
#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\)，不能直接用递推做。用矩阵快速幂快速计算。
简单的用矩阵优化的题都是把所有要用的东西都扔进一个矩阵里，比如这题：
其中 \(f_i\) 表示求得的数列 \(a\) 第 \(i\) 项的值。但要注意，由于这个矩阵是由 \(i-1\) 的矩阵推过来的，所以矩阵内能引用到的数分别是 \(a_{i-1},a_{i-2},a_{i-3}\)。
在递推式里，推导答案的式子有 \(m\) 项，则这个矩阵的大小为 \(m\times m\)。
我们需要一个矩阵 \(B\) 使得
在这题里，这个矩阵 \(B\) 长成下面这样：
也就是说，这个矩阵与上述矩阵相乘的意义为
我们可以比较容易的推出这个矩阵为
如果要用矩阵的出 \(f_i\)，那么大概是这样的
然后就可以直接用矩阵快速幂来做了。一共有 \(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[2][1]=x.mat[3][1]=1;
x=(a^(n-3))*x;//务必注意这里的顺序及其行/列向量的区别
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;
}
P5789 [TJOI2017] 可乐（数据加强版）
对于这种路径问题我们有个利器：邻接矩阵
邻接矩阵
在一开始的邻接矩阵中，\(A_{i,j}=1\) 当且仅当有一条边为 \(i\Rightarrow j\)，必须有一条边直接连接这两个点。这时候 \(A\) 存储的是图中任意两点间的路径中恰好经过 \(1\) 条边的个数。
再考虑 \(A^k\) 的意义是什么。众所周知若 \(C=A\times A\)，则有 \(\sum\limits_{k=1}^{n} A_{i,k}\times A_{k,j}\) 。这个式子表示的是以每一个 \(k\) 为中转站，\(C\) 的结果是任意两点经过恰好两条边的路径的个数。
所以，\(A_k\) 表示的是任意两点间恰好经过 \(k\) 条边的路径的个数。
回到这题，发现这就是一道邻接矩阵模板题，唯一需要处理的就是爆炸这个行为。我们弄一个 \(n+1\) 号节点，将 \(1\sim n\) 都建立一条向 \(n+1\) 号节点的单向边。因为机器人一旦爆炸就不能继续在节点之间移动了。
同时还要在所有节点弄一个自环（包括 \(n+1\) 号节点），因为有停留操作。
#include<iostream>
#include<cstring>
using namespace std;
#define mod 2017
struct Mat
{
#define N 110
int mat[N][N];
}a;
int n,m,u,v,t,l=101;
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++)
for(int j=1;j<=l;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<=l;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>>m;
for(int i=1;i<=m;i++)
{
cin>>u>>v;
a.mat[u][v]=a.mat[v][u]=1;
}
for(int i=1;i<=n+1;i++)
{
a.mat[i][i]=1;
a.mat[i][n+1]=1;
}
cin>>t;
l=n+1;
a=a^(t);
int ans=0;
// for(int i=1;i<=l;i++)
for(int j=1;j<=l;j++)
ans=(ans+a.mat[1][j])%mod;
cout<<ans<<endl;
}
