线性代数

本文仅发布于此博客和作者的洛谷博客，不允许任何人以任何形式转载，无论是否标明出处及作者。

声明：由于这是一个OI里的线性代数讲解而不是数学的讲解，所以这篇文章只会出现和OI高度相关的内容。

矩阵基础

矩阵是一个按照长方阵列排列的数集。记矩阵 $A_{n,m}= \left(\begin{matrix} a_{1,1} & a_{1,2} & \cdots & a_{1,m}\\ a_{2,1} & a_{2,2} & \cdots & a_{2,m}\\ \vdots & \vdots & \ddots & \vdots\\ a_{n,1} & a_{n,2} & \cdots & a_{n,m} \end{matrix}\right).$

单位矩阵

定义单位矩阵 $I_{n,n}=\left(\begin{matrix} 1 & 0 & \cdots & 0\\ 0 & 1 & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & 1 \end{matrix}\right).$

即除主对角线上元素为1外，其他位置上元素均为0的方阵。

零矩阵

定义零矩阵 $O_{n,m}=\left(\begin{matrix} 0 & 0 & \cdots & 0\\ 0 & 0 & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & 0 \end{matrix}\right).$

全零的矩阵。

转置矩阵

定义 $A_{n,m}= \left(\begin{matrix} a_{1,1} & a_{1,2} & \cdots & a_{1,m}\\ a_{2,1} & a_{2,2} & \cdots & a_{2,m}\\ \vdots & \vdots & \ddots & \vdots\\ a_{n,1} & a_{n,2} & \cdots & a_{n,m} \end{matrix}\right)$ 的转置矩阵为 $A^T= \left(\begin{matrix} a_{1,1} & a_{2,1} & \cdots & a_{m,1}\\ a_{1,2} & a_{2,2} & \cdots & a_{m,2}\\ \vdots & \vdots & \ddots & \vdots\\ a_{1,n} & a_{2,n} & \cdots & a_{m,n} \end{matrix}\right)$

即对 $A$ 关于主对角线翻转。

~~虽然零矩阵和转置矩阵没什么用~~

矩阵加法

两个矩阵相加时，在每一位上相加即可。 $A_{n,m}+B_{n,m}= \left(\begin{matrix} a_{1,1}+b_{1,1} & a_{1,2}+b_{1,2} & \cdots & a_{1,m}+b_{1,m}\\ a_{2,1}+b_{2,1} & a_{2,2}+b_{2,2} & \cdots & a_{2,m}+b_{2,m}\\ \vdots & \vdots & \ddots & \vdots\\ a_{n,1}+b_{n,1} & a_{n,2}+b_{n,2} & \cdots & a_{n,m}+b_{n,m} \end{matrix}\right).$

使用代码可能更好理解。

for(int i=1;i<=n;i++){
	for(int j=1;j<=m;j++){
		c[i][j]=a[i][j]+b[i][j];
	}
}

矩阵减法

同理。

for(int i=1;i<=n;i++){
	for(int j=1;j<=m;j++){
		c[i][j]=a[i][j]-b[i][j];
	}
}

注意，行数和列数相同的矩阵才能做加减法。

矩阵数乘

矩阵的数乘和乘法不同。矩阵 $A$ 和数 $k$ 相乘，相当于 $A$ 的每一位和 $k$ 相乘。

for(int i=1;i<=n;i++){
	for(int j=1;j<=m;j++){
		c[i][j]=a[i][j]*k;
	}
}

易得，上面的运算都满足交换律，结合律，复杂度均为 $O(n^2)$

矩阵乘法

给出两个矩阵 $A_{n,p},B_{p,m}$ ，相乘后得到一个矩阵 $C_{n,m},$ ，满足 $\large C_{i,j}=\sum\limits_{k=1}^{p} A_{i,k}\times B_{k,j}$ ，即“左行右列”。

for(int i=1;i<=n;i++){
	for(int j=1;j<=m;j++){
		c[i][j]=0;
		for(int k=1;p<=p;k++){
			c[i][j]+=a[i][k]*b[k][j];
		}
	}
}

复杂度 $O(n^3)$

注意，矩阵乘法必须保证 $A$ 的列数等于 $B$ 的行数，而且矩阵乘法并不满足交换律。

矩阵快速幂

link

方阵 $A_{n,n}$ 的 $k$ 次幂记为 $A^k$ ， $A^k=\underbrace{A\times A\times\cdots\times A}_{k\text{个}}.$

特殊地， $A^0$ 为单位矩阵 $I_{n,n}$ 。

暴力做时，复杂度 $O(n^3k)$ ，显然过高，尝试对其进行优化。

可以使用快速幂的方法。下面是标准的快速幂代码：

//a^k mod P
while(k){
	if(k&1){
		ans*=a;
		ans%=P;
	}
	a*=a;
	a%=P;
	k>>=1;
}

本质上是一个把 $k$ 次幂拆分成一些 $2$ 的倍数次幂相乘。如果连这个都看不懂了的话建议还是花点时间进行一次复健。

我们对矩阵求幂进行完全一模一样的操作即可，总复杂度乘上一个矩阵乘法， $O(n^3\log k)$

矩阵快速幂和上面的东西可以一起集成在一个程序里：

#include<bits/stdc++.h>
#define int long long
using namespace std;
int P=1e9+7;
void debug(){//不要在意这些细节
	cout<<"???what the hell"<<endl;
	exit(0);
}
int mod(int k){//严谨的取模
	return ((k%P)+P)%P;
}
struct matrix{
	int n;
	int m;
	int v[105][205];
	matrix(int t1,int t2){//初始化全0
		n=t1;
		m=t2;
		for(int i=1;i<=n;i++){
			for(int j=1;j<=m;j++){
				v[i][j]=0;
			}
		}
	}
};
matrix I(int n){//构造单位矩阵
	matrix c(n,n);
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n;j++){
			c.v[i][j]=0;
			if(i==j){
				c.v[i][j]=1;
			}
		}
	}
	return c;
}
matrix O(int n,int m){//构造零矩阵
	matrix c(n,m);
	return c;
}
matrix T(matrix c){//构造转置矩阵
	for(int i=1;i<=max(c.n,c.m);i++){
		for(int j=1;j<=i;j++){
			swap(c.v[i][j],c.v[j][i]);
		}
	}
	swap(c.n,c.m);
	return c;
}
matrix operator+(matrix a,matrix b){//矩阵加
	if(a.n!=b.n||a.m!=b.m){
		debug();
	}
	matrix c(a.n,a.m);
	for(int i=1;i<=c.n;i++){
		for(int j=1;j<=c.m;j++){
			c.v[i][j]=a.v[i][j]+b.v[i][j];
			c.v[i][j]=mod(c.v[i][j]);
		}
	}
	return c;
}
matrix operator-(matrix a,matrix b){//矩阵减
	if(a.n!=b.n||a.m!=b.m){
		debug();
	}
	matrix c(a.n,a.m);
	for(int i=1;i<=c.n;i++){
		for(int j=1;j<=c.m;j++){
			c.v[i][j]=a.v[i][j]-b.v[i][j];
			c.v[i][j]=mod(c.v[i][j]);
		}
	}
	return c;
}
matrix operator*(matrix a,int b){//矩阵数乘
	matrix c(a.n,a.m);
	for(int i=1;i<=c.n;i++){
		for(int j=1;j<=c.m;j++){
			c.v[i][j]=a.v[i][j]*b;
			c.v[i][j]=mod(c.v[i][j]);
		}
	}
	return c;
}
matrix operator*(matrix a,matrix b){//矩阵乘
	if(a.m!=b.n){
		debug();
	}
	matrix c(a.n,b.m);
	for(int i=1;i<=c.n;i++){
		for(int j=1;j<=c.m;j++){
			c.v[i][j]=0;
			for(int k=1;k<=a.m;k++){
				c.v[i][j]+=a.v[i][k]*b.v[k][j];
				c.v[i][j]=mod(c.v[i][j]);
			}
		}
	}
	return c;
}
matrix pow(matrix a,int k){//矩阵快速幂
	if(a.n!=a.m){
		debug();
	}
	matrix c(a.n,a.m);
	c=I(a.n);
	while(k!=0){
		if(k&1){
			c=c*a;
		}
		a=a*a;
		k>>=1;
	}
	return c;
}
void print(matrix c){
	for(int i=1;i<=c.n;i++){
		for(int j=1;j<=c.m;j++){
			cout<<c.v[i][j]<<' ';
		}
		cout<<endl;
	}
	return;
}
signed main(){
	
	return 0;
}

一个比较不是很全的矩阵全家桶就做好了。

这份代码先放放，处理一下矩阵加速。

矩阵加速

link，但是先不用看

定义斐波那契数列 $F_i=F_{i-1}+F_{i-2}$ ，其中 $F_1=F_2=1.$ 尝试在快于 $O(n)$ 的时间内计算出 $F_n \mod 10^9+7$ 。

构造一个矩阵 $\left(\begin{matrix} F_i & F_{i+1} \end{matrix}\right)$ ，可以发现矩阵变成 $\left(\begin{matrix} F_{i+1} & F_{i+2} \end{matrix}\right)$ 即 $\left(\begin{matrix} F_{i+1} & F_i+F_{i+1} \end{matrix}\right)$ 仅需乘上一个 $\left(\begin{matrix} 0 & 1 \\ 1 & 1 \end{matrix}\right)$ ，对于每次变换皆如此。

因此 $\left(\begin{matrix} F_{k-1} & F_{k} \end{matrix}\right) = \left(\begin{matrix} F_1 & F_2 \end{matrix}\right)\times \left(\begin{matrix} 0 & 1 \\ 1 & 1 \end{matrix}\right)^{k-2}$ ，而后可得 $F_k.$

$\left(\begin{matrix} 0 & 1 \\ 1 & 1 \end{matrix}\right)^{k-2}$ 使用矩阵快速幂，即可在 $O(\log k)$ 的复杂度下完成，所以整个过程为 $O(\log k)$ 。

模板题过程相似，快速幂的矩阵变为了 $3\times 3$ 的而已。

Extra

另外，可以尝试想一下这道题。

给出 $n$ 个操作，对于每个操作，输入两个数 $l,r$ ，计算斐波那契数列上区间 $[l,r]$ 的数字和 $\mod\space 10^{16}+2137$ 。斐波那契数列下标从 $1$ 开始，依次为 $1,1,2,3\cdots$ .

对于 $100\%$ 的数据：保证 $1\leq n\leq 10^6,1\leq l\leq r\leq 10^9$ 。

题解：

设斐波那契数列为 $F$ ,对于每个查询 $[l,r]$ ，通过矩阵加速递推 $O(\log l)$ 地找到 $F_l$ 和 $F_{l+1}$ 。设 $F_l=a,F_{l+1}=b,len=r-l+1.$ 另设数列 $G,G_i=G_{i-1}+G_{i-2}+1.$ 特殊地， $G_1=0,G_2=1$ 。仍然使用矩阵加速,在 $O(\log len)$ 的时间内求 $F_{len},G_{len}$ ，区间 $[l,r]$ 之和即为 $F_{len}\times a+G_{len}\times b.$

对于每一次操作，复杂度为 $O(\log r)$ ，总共为 $O(n\log r)$ .

高斯消元

初等行变换

对于任意矩阵 $A$ ,可以进行三种初等行变换：

倍加：把某一行所有元素乘上某个数并加到一行。（原来那行不变）
倍乘：某一行乘上非零数 $k$ 。
对换：两行元素互换。

注意到，对换可以由倍加操作完成。（不会的去看CSP初赛学怎么不用辅助空间交换元素）（虽然这个性质并没有什么用处，，，）

解线性方程组

模板题link

我们把给出的方程组变成矩阵的形式： $\left(\begin{matrix} a_{1,1} & a_{1,2} & \cdots & a_{1,n} & | & b_1\\ a_{2,1} & a_{2,2} & \cdots & a_{2,n} & | & b_2\\ \vdots & \vdots & \ddots & \vdots & | & \vdots\\ a_{n,1} & a_{n,2} & \cdots & a_{n,n} & | & b_n\\ \end{matrix}\right).$

通过初等行变换把左边的 $a$ 方阵变换成单位矩阵，则 $b_i$ 就是未知数的值。

具体实现方面，对于第 $i$ 列，找到绝对值最小且非零的 $a_{i,j}$ （为了少掉精度），把这一行和第 $i$ 行交换。然后用这一行消掉所有其他行的第 $i$ 个元素。

如果找不到非零 $a_{i,j}$ ，就无解。

代码：

#include<bits/stdc++.h>
#define int long long
using namespace std;
struct matrix{
	int n;
	int m;
	double v[105][205];//别忘了改double
	matrix(int t1,int t2){//初始化全0
		n=t1;
		m=t2;
		for(int i=1;i<=n;i++){
			for(int j=1;j<=m;j++){
				v[i][j]=0;
			}
		}
	}
};
signed main(){
	int n;
	cin>>n;
	matrix c(n,n+1);
	for(int i=1;i<=c.n;i++){
		for(int j=1;j<=c.m;j++){
			cin>>c.v[i][j];
		}
	}
	for(int j=1;j<=c.m-1;j++){//对每一列进行消元
		double max_=-0x3f3f3f3f;
		int maxpos=0;
		for(int i=j;i<=c.n;i++){//找到j列最大的a[i][j]（注意：j行之前的行已经是单位矩阵的形式了，不能碰）
			if(c.v[i][j]>max_){
				max_=c.v[i][j];
				maxpos=i;
			}
		}
		for(int k=1;k<=c.m;k++){//交换maxpos行和j行
			swap(c.v[j][k],c.v[maxpos][k]);
		}
		if(c.v[j][j]==0){//无解判定
			cout<<"No Solution"<<endl;
			return 0;
		}
		for(int i=1;i<=c.n;i++){//对其他行进行消元
			if(i==j){//当前行跳过
				continue;
			}
			double tmp=c.v[i][j]/c.v[j][j];//第j行全体乘上tmp再减到当前行上
			for(int k=1;k<=c.m;k++){
				c.v[i][k]-=c.v[j][k]*tmp;
			}
		}
		double tmp=c.v[j][j];
		for(int k=1;k<=c.m;k++){//a[j][j]化为1
			c.v[j][k]/=tmp;
		}
	}
	for(int i=1;i<=c.n;i++){//快乐输出
		cout<<fixed<<setprecision(2)<<c.v[i][c.m]<<endl;
	}
	return 0;
}

矩阵求逆

模板题link

定义方阵 $A$ 的逆矩阵 $A^{-1}$ ，使得 $AA^{-1}=I$ 。

$A$ 的逆矩阵不一定存在。

以 $A=\left(\begin{matrix} 2 & 1 & 1\\ 1 & 2 & 1\\ 1 & 1 & 2 \end{matrix}\right)$ 为例，逆矩阵的求法如下：

1.在 $A$ 的右边补上一个单位矩阵。

$C=\left(\begin{matrix} 2 & 1 & 1 & | & 1 & 0 & 0\\ 1 & 2 & 1 & | & 0 & 1 & 0\\ 1 & 1 & 2 & | & 0 & 0 & 1 \end{matrix}\right)$

对左边进行高斯消元，消成单位矩阵。如果无法消元即代表逆矩阵不存在。

$C'=\left(\begin{matrix} 1 & 0 & 0 & | & \frac34 & -\frac14 & -\frac14\\ 0 & 1 & 0 & | & -\frac14 & \frac34 & -\frac14\\ 0 & 0 & 1 & | & -\frac14 & -\frac14 & \frac34 \end{matrix}\right)$

此时右半边就是逆矩阵了。

$A^{-1}=\left(\begin{matrix} \frac34 & -\frac14 & -\frac14\\ -\frac14 & \frac34 & -\frac14\\ -\frac14 & -\frac14 & \frac34 \end{matrix}\right).$

注意模板题是输出模意义下的逆矩阵，把除法全部换成逆元（费马小定理）即可。

行列式求值

模板题link

方阵 $A$ 的行列式记为 $|A|$ 或 $\det(A)$ ，它具有一个数值，值的计算遵循一个比较复杂的公式。行列式也可以进行高斯消元，并且值也会相应变化。

倍加值不变。
倍乘值乘k。
对换（自己和自己换不算）值反号。

另外，“上三角矩阵”（即主对角线一下均为0）的行列式的值为对角线乘积。

根据这个，我们可以把矩阵消成“上三角矩阵”（消成单位矩阵一步到位也可以的）计算值就可以了。

实现时，当无法消元时，答案为0.

奉上一份常数爆炸根本过不去的全家桶（

（其实还会有WA，因为这道题模数不是质数。正解是在消元的时候进行“辗转相减”而不是直接逆元）

#include<bits/stdc++.h>
#define int long long
using namespace std;
int P=1e9+7;
void debug(){
	cout<<"???what the hell"<<endl;
	exit(0);
}
int mod(int k){
	return ((k%P)+P)%P;
}
int rev(int a){
    int ans=1;
    int k=P-2;
    while(k!=0){
        if(k&1){
            ans*=a;
            ans=mod(ans);
        }
        a*=a;
        a=mod(a);
        k>>=1;
    }
    return ans;
}
struct matrix{
	int n;
	int m;
	int v[105][105];
	matrix(int t1,int t2){
		n=t1;
		m=t2;
		for(int i=1;i<=n;i++){
			for(int j=1;j<=m;j++){
				v[i][j]=0;
			}
		}
	}
};
matrix I(int n){
	matrix c(n,n);
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n;j++){
			c.v[i][j]=0;
			if(i==j){
				c.v[i][j]=1;
			}
		}
	}
	return c;
}
matrix O(int n,int m){
	matrix c(n,m);
	return c;
}
matrix T(matrix c){
	for(int i=1;i<=max(c.n,c.m);i++){
		for(int j=1;j<=i;j++){
			swap(c.v[i][j],c.v[j][i]);
		}
	}
	swap(c.n,c.m);
	return c;
}
bool operator==(matrix a,matrix b){
	if(a.n!=b.n||a.m!=b.m){
		return false;
	}
	for(int i=1;i<=a.n;i++){
		for(int j=1;j<=a.m;j++){
			if(a.v[i][j]!=b.v[i][j]){
				return false;
			}
		}	
	}
	return true;
}
bool operator!=(matrix a,matrix b){
	return !(a==b);
}
matrix operator+(matrix a,matrix b){
	if(a.n!=b.n||a.m!=b.m){
		debug();
	}
	matrix c(a.n,a.m);
	for(int i=1;i<=c.n;i++){
		for(int j=1;j<=c.m;j++){
			c.v[i][j]=a.v[i][j]+b.v[i][j];
			c.v[i][j]=mod(c.v[i][j]);
		}
	}
	return c;
}
matrix operator-(matrix a,matrix b){
	if(a.n!=b.n||a.m!=b.m){
		debug();
	}
	matrix c(a.n,a.m);
	for(int i=1;i<=c.n;i++){
		for(int j=1;j<=c.m;j++){
			c.v[i][j]=a.v[i][j]-b.v[i][j];
			c.v[i][j]=mod(c.v[i][j]);
		}
	}
	return c;
}
matrix operator*(matrix a,int b){
	matrix c(a.n,a.m);
	for(int i=1;i<=c.n;i++){
		for(int j=1;j<=c.m;j++){
			c.v[i][j]=a.v[i][j]*b;
			c.v[i][j]=mod(c.v[i][j]);
		}
	}
	return c;
}
matrix operator*(matrix a,matrix b){
	if(a.m!=b.n){
		debug();
	}
	matrix c(a.n,b.m);
	for(int i=1;i<=c.n;i++){
		for(int j=1;j<=c.m;j++){
			c.v[i][j]=0;
			for(int k=1;k<=a.m;k++){
				c.v[i][j]+=a.v[i][k]*b.v[k][j];
				c.v[i][j]=mod(c.v[i][j]);
			}
		}
	}
	return c;
}
matrix pow(matrix a,int k){
	if(a.n!=a.m){
		debug();
	}
	matrix c(a.n,a.m);
	c=I(a.n);
	while(k!=0){
		if(k&1){
			c=c*a;
		}
		a=a*a;
		k>>=1;
	}
	return c;
}
pair<matrix,int> Gauss_jordan(matrix c){
	int ans=1;//行列式求值专用
	if(c.m<c.n){
		debug();
	}
	for(int j=1;j<=c.n;j++){
		int max_=-0x7fffffff;
		int maxpos=0;
		for(int i=j;i<=c.n;i++){
			if(c.v[i][j]>max_&&c.v[i][j]!=0){
				max_=c.v[i][j];
				maxpos=i;
			}
		}
		if(maxpos!=j){
			ans*=-1;//换行就换号
			ans=mod(ans);
		}
		for(int k=1;k<=c.m;k++){
			swap(c.v[j][k],c.v[maxpos][k]);
		}
		if(c.v[j][j]==0){
			matrix ERR(c.n,c.m);
			return pair<matrix,int>{ERR,-0x7fffffff};
		}
		for(int i=1;i<=c.n;i++){
			if(i==j){
				continue;
			}
			int tmp=c.v[i][j]*rev(c.v[j][j]);
			tmp=mod(tmp);
			for(int k=1;k<=c.m;k++){
				c.v[i][k]-=c.v[j][k]*tmp;
				c.v[i][k]=mod(c.v[i][k]);
			}
		}
		int tmp=rev(c.v[j][j]);
		ans*=c.v[j][j];//行列式/k,答案*k
		ans=mod(ans);
		for(int k=1;k<=c.m;k++){
			c.v[j][k]*=tmp;
			c.v[j][k]=mod(c.v[j][k]);
		}
	}
	return pair<matrix,int>{c,ans};
}
matrix rev(matrix a){
	if(a.n!=a.m){
		debug();
	}
	matrix c(a.n,a.m*2);//补单位矩阵
	for(int i=1;i<=a.n;i++){
		for(int j=1;j<=a.m;j++){
			c.v[i][j]=a.v[i][j];
		}
	}
	matrix tmp=I(a.n);
	for(int i=1;i<=tmp.n;i++){
		for(int j=1;j<=tmp.m;j++){
			c.v[i][j+tmp.n]=tmp.v[i][j];
		}
	}
	c=Gauss_jordan(c).first;//得到消元后矩阵
	matrix ans(a.n,a.m);
	for(int i=1;i<=a.n;i++){
		for(int j=1;j<=a.m;j++){
			a.v[i][j]=c.v[i][j+a.n];
		}
	}
	return a;//拿到右半边
}
int det(matrix c){//直接集成在高斯消元里面算
	return Gauss_jordan(c).second;
}
void scan(matrix &c){
	for(int i=1;i<=c.n;i++){
		for(int j=1;j<=c.m;j++){
			cin>>c.v[i][j];
		}
	}
	return;
}
void print(matrix c){
	for(int i=1;i<=c.n;i++){
		for(int j=1;j<=c.m;j++){
			cout<<c.v[i][j]<<' ';
		}
		cout<<endl;
	}
	return;
}
signed main(){
	
	return 0;
}

posted @ 2023-03-06 21:23 Cerebral_fissure 阅读(130) 评论(0) 编辑收藏举报

刷新页面返回顶部

登录后才能查看或发表评论，立即登录或者逛逛博客园首页

相关博文：

· 线段树题单

· 题解：BLO-Blockade

· 网课-线性代数学习笔记

· 线性代数（矩阵/高斯消元）

· 一些些数学的算法笔记

阅读排行：
· 地球OL攻略 —— 某应届生求职总结
· 周边上新：园子的第一款马克杯温暖上架
· Open-Sora 2.0 重磅开源！
· 提示词工程——AI应用必不可少的技术
· .NET周刊【3月第1期 2025-03-02】

公告

昵称： Cerebral_fissure
园龄： 2年2个月
粉丝： 4
关注： 0

+加关注

2025年3月

日

一

二

三

四

五

六

abv3Rpkg

花哨的垃圾场。

线性代数

矩阵基础

单位矩阵

零矩阵

转置矩阵

矩阵加法

矩阵减法

矩阵数乘

矩阵乘法

矩阵快速幂

矩阵加速

高斯消元

初等行变换

解线性方程组

矩阵求逆

行列式求值

公告

搜索

常用链接

最新随笔

我的标签

积分与排名

随笔档案

文章档案

阅读排行榜

推荐排行榜