线性代数基础
矩阵
矩阵可以理解为由数字组成的方形，比如：
我们使用行和列表示元素在矩阵中的位置，比如对于示例中的矩阵，数字 \(4\) 位于第 \(1\) 行第 \(3\) 列，记作 \(A_{1,3}=4\)。
矩阵的每一行和每一列都可以视为一个 \(n\) 维向量，\(n\) 为行数/列数。例如，\(A\) 的第一行可以视为一个向量 \(\mathbf{a}=(1,1,4,5)\)。
在 OI 中，我们一般使用的都是方阵，即行数与列数相同的矩阵。
矩阵的线性运算就是对应位置的元素的线性运算，有加减和数乘，故两个行数和列数相同的矩阵才能进行线性运算。
矩阵乘法
矩阵的乘法是向量内积的推广。矩阵乘法必须满足第一个矩阵的列数和第二个矩阵的行数相同。下文中提到的矩阵乘法都默认满足这个条件。
矩阵乘法的法则满足行向量乘列向量，即得出的矩阵的第 \(i\) 行第 \(j\) 列是第一个矩阵的第 \(i\) 个行向量和第 \(j\) 个列向量的内积，所以得出的矩阵的列数与第二个矩阵的列数相同，行数与第一个矩阵的行数相同，因为第一个矩阵只提供行向量，第二个矩阵只提供列向量。
设 \(A\times B=C\)，\(A\) 的列数和 \(B\) 的行数为 \(n\)，则
例：
得出的矩阵第 \(1\) 行第 \(1\) 列为 \((1,1,4,5)\cdot(1,4,7,2)=43\)，故第 \(1\) 行第 \(1\) 列为 \(43\)。
计算两个 \(n\times n\) 矩阵乘法朴素做法时间复杂度为 \(O(n^3)\)，已知最快的时间复杂度记作 \(O(n^\omega)\)，Strassen 算法可以做到 \(\omega< \log_2 7\approx 2.807\)，截至本文写作前已经有 \(\omega< 2.371339\)。但是因为实现难度和常数因子的影响，现实中多使用 \(O(n^3)\) 的朴素算法结合常数优化。
矩阵乘法具有结合律，即对于矩阵 \(A,B,C\)，有：
这使得矩阵可以和快速幂、线段树等结合实现多种功能。
注意：矩阵乘法没有交换律！
单位矩阵是一种特殊的矩阵，即 \(A_{i,i}=1\)，其他位置都等于 \(0\) 的矩阵。不难发现任何矩阵乘它还得原来的矩阵。类似 \(1\) 在乘法中的作用，一般用 \(I\) 表示。
加速线性递推
以斐波那契数列为例，我们知道斐波那契数列的递推式：\(F_1=F_2=1\)，\(F_i=F_{i-1}+F_{i-2}(i\ge 3)\)。
则可以将递推写成矩阵形式：
则数列第 \(n\) 项就是
的第 \(1\) 行第 \(1\) 列的元素。
可以使用快速幂优化。
习题：
表达修改
可以使用矩阵
表示一个水晶球。以操作 \(A_i\gets A_i+v\) 为例：
类似地，其余各种操作都可以通过进行一次矩阵乘法实现。直接线段树维护区间乘即可。
习题：
定长路径计数
一个 \(n\) 个点的无边权图，则其邻接矩阵 \(A\) 的 \(k\) 次幂的含义为：\(A_{u,v}^k\) 就是从 \(u\) 到 \(v\) 走 \(k\) 步的方案数。
习题：
定长最短路
首先改造一下矩阵乘法：
原本矩阵乘法是行向量乘列向量，即 \((+,\times)\) 矩阵，而我们现在使用 \((\min,+)\) 矩阵，即
这样计算邻接矩阵 \(A\) 的 \(k\) 次幂得到的 \(A_{u,v}^k\) 就是 \(u\) 到 \(v\) 经过 \(k\) 条边的最短路。
习题：
常数优化
矩阵乘法时可以将三重循环中第二、三层交换顺序，内存访问更加连续。
如果是含有大量 \(0\) 的稀疏矩阵，可以跳过 \(0\) 所在位置。
可以进行一定程度的循环展开。
如果矩阵运算需要取模，尽量减少取模次数。
线性方程组
现有线性方程组：
尝试求解。
我们将系数填成一个方阵，在右侧补上常数项构成一个矩阵：
高斯消元
高斯消元的目标就是将矩阵的系数部分通过初等行变换变为上三角矩阵：
即将左下部分的系数全部化为 \(0\)。这样我们已经知道了 \(x_n\)，可以通过不断回代求出所有 \(x\)。
每行选择可以选的最大的数作为主元以减小浮点数精度误差。
如果某一行 \(a'_{i,j}\) 均为 \(0\)，且 \(b'_{i}\) 不为 \(0\) 则方程组无解，若 \(b'_{i}\) 为 \(0\) 则有无数组解。
时间复杂度 \(O(n^3)\)。
高斯-约旦消元
与高斯消元类似，同样是利用矩阵的初等行变换，但是这里我们希望得到一个对角矩阵：
这样我们就省去了回代操作。
可以用于矩阵求逆。
时间复杂度 \(O(n^3)\)。
习题：
线性基
这里讲解异或线性基。
线性基可以快速查询：
-
一个数是否可以被一个集合中的数通过异或运算得到。
-
一个集合中的数通过异或运算可以得到的最大/最小/第 \(k\) 小值。
对于值域为 \([1,V]\) 的集合，我们可以用 \(\lceil\log_2 V\rceil\) 个数描述一个线性基。其中第 \(i\) 个数的要么最高位是第 \(i\) 位，要么就是 \(0\)。
以下讲解贪心法构造线性基：
将一个数 \(x\) 插入线性基，从 \(x\) 的二进制最高位开始考虑，设当前考虑的是第 \(i\) 位：
-
若线性基第 \(i\) 个数为 \(0\)，将这个位置的数设为 \(x\)，结束插入操作。
-
否则设线性基第 \(i\) 个数为 \(a_i\)，\(x\gets x\oplus a_i\)，重新寻找 \(x\) 的最高位，重复以上过程。
若插入结束后 \(x=0\)，说明原本线性基就已经可以表示 \(x\) 这个数。
查询最小值只需输出线性基中的最小值即可。
查询最大值从高位开始遍历，按位贪心即可。
以上这些操作的单次时间复杂度都是 \(O(\log V)\)。
合并两个线性基可以直接暴力将其中一个的所有元素插入到另一个线性基中，时间复杂度 \(O(\log^2 V)\)。
习题：
代码模板
矩阵快速幂
#include<cstring>
typedef long long ll;
constexpr int N=105,mod=998'244'353;
int n;
inline void add(int &a,int b){
a+=b;
if(a>=mod) a-=mod;
}
struct Matrix{
int a[N][N];
Matrix operator*(const Matrix &x)const{
Matrix res;
memset(res.a,0,sizeof(res.a));
for(int i=1;i<=n;i++)
for(int k=1;k<=n;k++)
for(int j=1;j<=n;j++)
add(res.a[i][j],ll(a[i][k])*x.a[k][j]%mod);
return res;
}
}Im;
void init(){
for(int i=1;i<=n;i++) Im.a[i][i]=1;
}
Matrix qpow(Matrix a,int b){
Matrix res=Im;
while(b){
if(b&1) res=res*a;
a=a*a;
b>>=1;
}
return res;
}
高斯消元
#include<algorithm>
#include<random>
using namespace std;
constexpr int N=110;
constexpr double eps=1e-9;
mt19937 rnd(__builtin_ia32_rdtsc());
int n;
double a[N][N],ans[N];
void solve(){
shuffle(a+1,a+1+n,rnd);
int rk=0;
for(int i=1;i<=n;i++){
int cur=rk+1;
for(int j=rk+1;j<=n;j++)
if(fabs(a[cur][i])<fabs(a[j][i])) cur=j;
if(fabs(a[cur][i])<eps) continue;
if(cur!=rk+1) swap(a[cur],a[rk+1]);
double div=a[rk+1][i];
for(int j=i;j<=n+1;j++) a[rk+1][j]/=div;
for(int j=1;j<=n;j++){
if(j==rk+1||fabs(a[j][i])<eps) continue;
div=a[j][i];
for(int k=i;k<=n+1;k++) a[j][k]-=a[rk+1][k]*div;
}
++rk;
}
for(int i=rk+1;i<=n;i++)
if(fabs(a[i][n+1])>eps){
//No solution
return;
}
if(rk<n){
//Infinite solution
return;
}
for(int i=n;i;i--){
ans[i]=a[i][n+1];
for(int j=i+1;j<=n;j++)
ans[i]-=a[i][j]*ans[j];
}
}
线性基
typedef long long ll;
ll base[60];
inline void insert(ll x){
for(int i=63;i>=0;i--){
if(~x>>i&1) continue;
if(base[i]) x^=base[i];
else{
base[i]=x;
return;
}
}
}
inline ll querymax(){
ll res=0;
for(int i=63;i>=0;i--)
res=max(res,res^base[i]);
return res;
}
拓展阅读：
浙公网安备 33010602011771号