浅谈矩阵 学习笔记

浅谈矩阵

矩阵是线性代数部分的基础，一种常用的数据类型，让我们一起学习矩阵和其相关算法。

一、定义

直接背定义未免有些晦涩难懂，读者可以先读后面内容，遇到再向前翻。

• $n\times m$矩阵：形如$n$行$m$列的二维数组，但可以进行一系列运算。
• $n$阶方阵：$n\times n$的矩阵。
• **行/列矩阵：**只有一行/列的矩阵。
• 主对角线：对于一个$n$阶方阵$A$，主对角线指所有$A_i{,}_i$上的元素。
• 上/下三角矩阵：主对角线以左下/右上均为$0$的矩阵。
• 单位矩阵（I）：主对角线上所有元素值为$1$，其余为$0$的矩阵。
• 逆矩阵（P）：使$A\times P=I$的矩阵，表示为$A^{-1}$。
• 矩阵行列式($\det$)：$\det (A)={\sum }_{\sigma \in S_n}sgn(\sigma ){\prod }_{i=1}^n{a}_{i,\sigma (i)}$，$S_n$为$n$的全排列集合，$\sigma$即为其中一个全排列，$\sigma$有偶数个逆序对时$\mathrm{sgn}(\sigma )$为$1$，否则为$-1$，行列式为矩阵所有列向量夹几何体的体积。

二、运算

令$C$为运算结果。

• 加法：对于两个$n\times m$的矩阵$A,B$，$C_i{,}_j=A_i{,}_j,+B_i{,}_j$;

• 减法：同理，$C_i{,}_j=A_i{,}_j-B_i{,}_j$；

• 乘法：

  • 数与矩阵相乘

    对于任意矩阵$A$和任意常数$x$，$C_i{,}_j=A_i{,}_j\times x$；

  • 矩阵与矩阵相乘

    对于$P\times M$矩阵$A$和$M\times Q$矩阵$B$，则$C$大小为$P\times Q$，$C_i{,}_j={\sum }_{k-1}^M{A}_i{,}_k\times B_k{,}_j$，即**$C$的第$i$行第$j$个数是$A$的第$i$行与$B$的第$j$列对应位置相乘求和**。

    struct mat {
    	int n, m, a[105][105];
    	mat operator * (const mat& B) const {
    		mat A = *this, C;
    		C.n = A.n, C.m = B.m;
    		for (int i = 1; i <= n; ++i) {
    			for (int j = 1; j <= m; ++j) {
    				C.a[i][j] = 0;
    				for (int k = 1; k <= A.m; ++k) {
    					C.a[i][j] += A.a[i][k] * B.a[k][j];
    				}
    			}
    		}
    		return C;
    	}
    };
    

  • 性质

    矩阵之间乘法运算满足交换律，不满足交换律；

三、应用

1. 矩阵加速递推

对于递推类的题目，我们往往可以求出形如$f_i={\sum }_j{a}_j{f}_j^{k_j}$的递推式，但因为是多项式，需要从前到后计算每一位，若每次递推消耗$O(1)$，共消耗$O(n)$。有了矩阵，我们可以将其转化为单项乘方的性质，利用快速幂算法可以在$O(nlgn)$​的时间内求出。

---

具体的，我们将每次递推需要传递的参数写成$1\times n$行矩阵的形式，每次乘上同一个$n\times n$的方阵，得到下一个行矩阵，答案需要从最后一个行矩阵中得到。

例如fib数列，$f_1=f_2=1,f_n(n\ge 3)=f_{n-1}+f_{n-2}$。我们先构造参数矩阵$A$，初始为$\left(\begin{array}{cc}{f}_1& f_2\end{array}\right)$，接下来构造$B=\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)$，使得$A\times B=\left(\begin{array}{cc}{f}_2& f_3\end{array}\right)$，则
$$\{\begin{array}{l}{f}_1\times a+f_2\times c=f2\\ f_1\times b+f_2\times d=f3=f1+f2\end{array}$$
解得$B=\left(\begin{array}{cc}0& 1\\ 1& 1\end{array}\right)$，答案即$A\times B^{n-2}$的第一个数。

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
struct mat {
	ll a[2][2];
} res, x;
inline mat mul(mat A, mat B) {
	mat C;
	for (int i = 0; i < 2; ++i) {
		for (int j = 0; j < 2; ++j) {
			C.a[i][j] = 0;
			for (int k = 0; k < 2; ++k) {
				C.a[i][j] += A.a[i][k] * B.a[k][j] % P;
			}
		}
	}
	return C;
}
int main() {
	ios :: sync_with_stdio(0);
	cin.tie(0), cout.tie(0);
	ll n;
	cin >> n;
	res.a[0][0] = 1, res.a[0][1] = 0, res.a[1][0] = 0, res.a[1][1] = 1;
	x.a[0][0] = 1, x.a[0][1] = 1, x.a[1][0] = 1, x.a[1][1] = 0;
	while (n) {
		if (n & 1) {
			res = mul(res, x);
		}
		x = mul(x, x), n >>= 1;
	}
	cout << res.a[0][1] << '\n';
	return 0;
}

2. 矩阵表达修改

直接看题目。

对于长度为$n$的两个数组$a$和$b$，对于$i\in [l_i,r_i]$维护三种操作，要求时间复杂度$O(nlgn)$：

1. $a_i=a_i+b_i$；
2. $b_i=a_i+b_i$；
   3. 求$\sum a_i\times b_i$。

很明显我们可以用线段树维护，但是对于前两个操作，每个数加的数值都不一样，如何处理？

我们考虑用共同乘一个矩阵的方式维护，对于每个区间维护矩阵$\left(\begin{array}{ccc}\sum a_i^2& \sum a_i\times b_i& \sum b_i^2\end{array}\right)$，后续用$x,y,z$代替。

为什么维护这个矩阵呢？若$a$加上了$b$，则$ab$的值就会加上$b^2$，否则就会加$a^2$，即$y$的值加上了$z$或$x$。而且$x$和$z$的维护十分方便，$(a+b{)}^2=a^2+2ab+b^2$，即$x+2y+z$。

经过一系列推导可得，

操作1乘的矩阵为$\left(\begin{array}{ccc}1& 2& 1\\ 0& 1& 1\\ 0& 0& 1\end{array}\right)$，操作2乘的矩阵为$\left(\begin{array}{ccc}1& 0& 0\\ 1& 1& 0\\ 1& 2& 1\end{array}\right)$，其余过程不过多阐述。

#include <bits/stdc++.h>
using namespace std;
const int N = 2e5 + 5, P = 10301;
struct mat {
    int M[3][3];
    mat operator * (const mat &b) const {
        mat c;
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
            	c.M[i][j] = 0;
                for (int k = 0; k < 3; ++k) {
                    (c.M[i][j] += M[i][k] * b.M[k][j] % P) %= P;
                }
            }
        }
        return c;
    }
} I, A[2];
struct node {
	int x[3];
	mat lazy;
} t[N << 2];
int a[N], b[N];
inline void pushup(node &root, node lch, node rch) {
    for (int i = 0; i < 3; ++i) {
        root.x[i] = (lch.x[i] + rch.x[i]) % P;
    }
}
inline void add(node &t, mat Mat) {
	node tmp;
	for (int i = 0; i < 3; ++i) {
		tmp.x[i] = 0;
		for (int j = 0; j < 3; ++j) {
			(tmp.x[i] += Mat.M[i][j] * t.x[j] % P) %= P;
		}
	}
	tmp.lazy = Mat * t.lazy;
	t = tmp;
}
inline void pushdown(int id) {
	mat Mat = t[id].lazy;
	add(t[id << 1], Mat), add(t[id << 1 | 1], Mat);
	t[id].lazy = I;
}
void build(int id, int l, int r) {
    t[id].lazy = I;
	if (l == r) {
        t[id].x[0] = a[l] * a[r] % P, t[id].x[1] = a[l] * b[r] % P, t[id].x[2] = b[l] * b[r] % P;
		return ;
	}
	int mid = l + r >> 1;
	build(id << 1, l, mid), build(id << 1 | 1, mid + 1, r);
	pushup(t[id], t[id << 1], t[id << 1 | 1]);
}
inline void f(int id, int op) {
    if (!op) {
		(t[id].x[0] += (t[id].x[1] << 1) + t[id].x[2]) %= P, (t[id].x[1] += t[id].x[2]) %= P;
	} else {
		(t[id].x[2] += (t[id].x[1] << 1) + t[id].x[0]) %= P, (t[id].x[1] += t[id].x[0]) %= P;
	}
}
void update(int id, int l, int r, int x, int y, int op) {
	if (x <= l && r <= y) {
        f(id, op), t[id].lazy = A[op] * t[id].lazy;
		return ;
	}
	pushdown(id);
	int mid = l + r >> 1;
	if (x <= mid) {
        update(id << 1, l, mid, x, y, op);
    }
	if (y > mid) {
        update(id << 1 | 1, mid + 1, r, x, y, op);
    }
	pushup(t[id], t[id << 1], t[id << 1 | 1]);
}
node query(int id, int l, int r, int x, int y) {
	pushdown(id);
	if (x <= l && r <= y) {
        return t[id];
    }
	int mid = l + r >> 1;
	if (y <= mid) {
        return query(id << 1, l, mid, x, y);
    }
	if (x > mid) {
        return query(id << 1 | 1, mid + 1, r, x, y);
    }
	node res;
	pushup(res, query(id << 1, l, mid, x, y), query(id << 1 | 1, mid + 1, r, x, y));
	return res;
}
int main() {
	I =    {{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}};
	A[0] = {{{1, 2, 1}, {0, 1, 1}, {0, 0, 1}}};
	A[1] = {{{1, 0, 0}, {1, 1, 0}, {1, 2, 1}}};
	int n, m;
	scanf("%d %d", &n, &m);
    for (int i = 1; i <= n; ++i) {
    	scanf("%d %d", a + i, b + i);
        a[i] %= P, b[i] %= P;
    }
	build(1, 1, n);
	while (m--) {
		int op, l, r;
        cin >> op >> l >> r;
        if (op == 1 || op == 2) {
            update(1, 1, n, l, r, op - 1);
        } else {
        	printf("%d\n", query(1, 1, n, l, r).x[1]);
        }
	}
	return 0;
}

四、算法

1. 高斯消元

高斯消元法是线性代数中一个重要的部分，主要原理为通过行/列变换将矩阵在性质不变的情况下转化为上三角实现消元。

---

1.1 高斯消元法求解线性方程组

我们直接看一个例子。
$$\{\begin{array}{l}{x}_1+3x_2+4x_3=5\\ x_1+4x_2+7x_3=3\\ 9x_1+3x_2+2x_3=2\end{array}$$
我们以系数矩阵的形式存储，即：
$$\left(\begin{array}{ccc}1& 3& 4\\ 1& 4& 7\\ 9& 3& 2\end{array}|\begin{array}{c}5\\ 3\\ 2\end{array}\right)$$
我们将第一列值绝对值最大的行放在第一行（读到后面再回头看这一步，会发现是用来化成上三角的），再将第一列除第一行的所有值化成$0$。为什么说除了第一行呢？因为我们要借助第一行消元，接下来称第$i$行为$r_i$，这里用乘法是为了在程序中减小误差。
$$\left(\begin{array}{ccc}9& 3& 2\\ 1& 4& 7\\ 1& 3& 4\end{array}|\begin{array}{c}2\\ 3\\ 5\end{array}\right)\stackrel{r_2\times 9-r_1,r_3\times 9-r_1}{\to }\left(\begin{array}{ccc}9& 3& 2\\ 0& 33& 61\\ 0& 24& 34\end{array}|\begin{array}{c}2\\ 25\\ 43\end{array}\right)$$
然后处理第二列——将第二列除第一行值绝对值最大的行放在第二行（在这个案例中没有变化），再将第二列除前两行的所有值化成$0$：
$$\left(\begin{array}{ccc}9& 3& 2\\ 0& 33& 61\\ 0& 24& 34\end{array}|\begin{array}{c}2\\ 25\\ 43\end{array}\right)\stackrel{r_3\times \frac{11}{8}-r_2}{\to }\left(\begin{array}{ccc}9& 3& 2\\ 0& 33& 61\\ 0& 0& -\frac{57}{4}\end{array}|\begin{array}{c}2\\ 25\\ \frac{273}{8}\end{array}\right)$$

我们发现主对角线以左下未知数均被消除，那化成这样有什么用呢？为了直观，我们再将其化为方程组形式：
$$\{\begin{array}{l}9x_1+3x_2+2x_3=2...①\\ 33x_2+61x_3=25...②\\ -\frac{57}{4}{x}_3=\frac{273}{8}...③\end{array}$$
发现③式被化为了一元一次方程组，可解出$x_3$；将$x_3$代入到②式可解出$x_2$；将$x_2$和$x_3$代入到①式即可解出$x_1$。对于本方程组，
$$\{\begin{array}{l}{x}_1=-\frac{37}{38}\\ x_2=\frac{197}{38}\\ x_3=-\frac{91}{38}\end{array}$$
总结一下上文中的斜体变量，我们可以得到高斯消元的操作步骤，我们设矩阵第$i$行$j$列为$a_i{,}_j$，当前处理第$i$列：

1. 将方程组转化为系数矩阵形式；
2. 对于每一列我们分别处理，设正在处理第$i$列，先将矩阵中第$i$列值绝对值最大的一行（设其原为第$j$行）放到第$i$行，形式化地，对于$\left|a_j{,}_i\right|$最大的$j$，$swap(a_i,a_j)$；
3. 利用此行将后面行此列的系数消为$0$，形式化地，对于 k ∈ { j + 1 , j + 2 , . . . , n } k\in \left\{j+1,j+2,...,n\right\} k∈{j+1,j+2,...,n}，$a_k=a_k\times \frac{a_i{,}_i}{a_k{,}_i}-a_i$；
4. 重复第$2$,$3$步，直到处理完所有列，即矩阵化为上三角；
5. 从最后一行开始，每行解出一个未知数并带到上一行，直到解出所有未知数，总复杂度为$O(n^3)$。

现在考虑如何判断无解和无穷多个解。

• 无解：当且仅当常数$a\ne 0$时，等式为$0=a$，即某一行出现$\left(\begin{array}{ccccc}0& 0& 0& ...& 0\end{array}|\begin{array}{c}a\end{array}\right)$；
• 无穷多解：即出现自由元，某一行出现$\left(\begin{array}{ccccc}0& 0& 0& ...& 0\end{array}|\begin{array}{c}0\end{array}\right)$或者某列出现$\left(\begin{array}{c}0\\ 0\\ 0\\ ⋮\\ 0\end{array}\right)$。

#include <bits/stdc++.h>
using namespace std;
const int N = 1e3 + 5;
const double eps = 1e-7;    // 很接近0的数，用于忽略误差
int m, n;
double a[N][N], x[N];
int Gauss(int m, int n) {
    int col, row;
    for (col = 0, row = 0; row < m && col < n; ++row, ++col) {
        int maxrid = row;   // 目前列值绝对值最大的行的行号
        for (int i = row + 1; i < m; ++i) {
            if (fabs(a[i][col]) > fabs(a[maxrid][col])) {
                maxrid = i;
            }
        }
        if (fabs(a[maxrid][col]) < eps) {   // 自由元
            --row;
            continue;
        }
        if (maxrid != row) { // STEP2
            for (int i = col; i <= n; ++i) {
                swap(a[row][i], a[maxrid][i]);
            }
        }
        for (int i = row + 1; i < m; ++i) { // STEP3
            if (fabs(a[i][col]) > eps) {
                double t = a[i][col] / a[row][col];
                for (int j = col; j <= n; ++j) {
                    a[i][j] -= a[row][j] * t;
                }
                a[i][col] = 0;
            }
        }
    }
    for (int i = row; i < m; ++i) {
        if (fabs(a[i][n]) > eps) {  // 参数为0结果非0
            return -1;  // 无解返回-1
        }
    }
    if (row < n) {
        return n - row; // 无穷多解：返回自由元个数
    }
    for (int i = n - 1; ~i; --i) {
        double t = a[i][n];
        for (int j = i + 1; j < n; ++j) {
            t -= x[j] * a[i][j];
        }
        x[i] = t / a[i][i];
    }
    return 0;   // 有唯一解
}
int main() {
	ios :: sync_with_stdio(0);
    cin.tie(0), cout.tie(0);
    cin >> m >> n;
    for (int i = 0; i < m; ++i) {   // 以系数矩阵形式读入
        for (int j = 0; j <= n; ++j) {
            cin >> a[i][j];
        }
    }
    int res = Gauss(m, n);
    if (!res) {
        for (int i = 0; i < n; ++i) {
            cout << 'x' << i + 1 << " = " << x[i] << '\n';
        }
    } else if (~res) {
        cout << "Infinite solutions\n";
    } else {
        cout << "No solution\n";
    }
    return 0;
}

1.2 高斯消元法求矩阵行列式

对于行列式的定义，有些同学可能还不太理解，这里我们用代数方式形象说明下，对于$n$阶矩阵$a$，其行列式即为第一列和所有数每次向右下分别移动一个，超出就取模，移动$n-1$次，共$n$个数相乘求和再减去最后一列向右上的结果，即如下左图中所有同色线段上数相乘求和减去右图结果的：

$-$

如：
$$\begin{gathered} a=\left(\begin{array}{ccc}1& 2& 3\\ 4& 5& 6\\ 7& 8& 10\end{array}\right) \\ det(a)=1\times 5\times 10+4\times 8\times 3+7\times 2\times 6-7\times 5\times 3-8\times 6\times 1-10\times 4\times 2=-3 \end{gathered}$$
这里需要用到行列式的四个性质，这里几何含义很容易理解，我们用代数方式来理解：

1. 转置，行列式不变；

   很明显，虽然第一列和最后一行的数有发生改变，但所有线上数集合未变，故行列式不变。

2. 两行/列交换，行列式取反；

   举例分析，会发现原来行列式计算中所有线的符号都发生了改变，故行列式取反。

3. 两行/列相加/减，行列式不变；

   • 引理1：若$A+B=C$，$\det (A)+\det (B)=\det (c)$；

     可用乘法分配律进行证明，此处不细讲。

   • 引理2：若矩阵中有一行/列值均为$0$，矩阵行列式为$0$。

     这样每条线中都会有$0$，则每条线的乘积也均为$0$，若干个$0$减若干个$0$，结果还是$0$。
     $$\left(\begin{array}{cccc}{a}_1{,}_1& a_1{,}_2& ...& a{,}_1{,}_n\\ ...& & & \\ a_k{,}_1+a_m{,}_1& a_k{,}_2+a_m{,}_2& ...& a_k{,}_n+a_m{,}_n\\ ...& & & \\ a_n{,}_1& a_n{,}_2& ...& a_n{,}_n\end{array}\right)=\left(\begin{array}{cccc}{a}_1{,}_1& a_1{,}_2& ...& a{,}_1{,}_n\\ ...& & & \\ a_k{,}_1& a_k{,}_2& ...& a_k{,}_n\\ ...& & & \\ a_n{,}_1& a_n{,}_2& ...& a_n{,}_n\end{array}\right)+\left(\begin{array}{cccc}0& 0& ...& 0\\ ...& & & \\ a_m{,}_1& a_m{,}_2& ...& a_m{,}_n\\ ...& & & \\ 0& 0& ...& 0\end{array}\right)$$

   根据两个引理：
   $$de{t}^{\prime }=det+0=det$$

4. 一行乘以某个数，行列式等比例扩大。

   可以发现每条线的乘积都等比例扩大，则结果等比例扩大。

因为高斯消元后矩阵性质不变，而且主对角线以左下元素均为$0$，可以先消元，答案由主对角线元素之积和交换行的数量决定，具体的，答案为主对角线上元素乘积（其余均为$0$不用管），若交换行奇数次则取反（只有交换行才能取反，两次则抵消）。

#include <bits/stdc++.h>
using namespace std;
const int N = 1e3 + 5;
const double eps = 1e-7;
int n;
double a[N][N];
double Gauss() {
	double det = 1;
	for (int i = 0; i < n; ++i) {
		int maxrid = i;
		for (int j = i + 1; j < n; ++j) {
			if (fabs(a[j][i]) > fabs(a[maxrid][i])) {
				maxrid = j;
			}
		}
		if (fabs(a[maxrid][i]) < eps) {
			return 0;
		}
		for (int j = 0; j < n; ++j) {
			swap(a[i][j], a[maxrid][j]);
		}
		if (maxrid != i) {
			det = -det;
		}
		det *= a[i][i];
		for (int j = i + 1; j < n; ++j) {
			a[i][j] /= a[i][i];
		}
		for (int j = 0; j < n; ++j) {
			if (i != j && fabs(a[j][i]) > eps) {
				for (int k = i + 1; k < n; ++k) {
					a[j][k] -= a[i][k] * a[j][i];
				}
			}
		}
	}
	return det;
}
int main() {
	ios :: sync_with_stdio(0);
	cin.tie(0), cout.tie(0);
	cin >> n;
	for (int i = 0; i < n; ++i) {
		for (int j = 0; j < n; ++j) {
			cin >> a[i][j];
		}
	}
	cout << Gauss() << '\n';
	return 0;
}

1.3 高斯消元法求逆矩阵

• 将$A$和$I$放在一个矩阵里，即$\left(AI\right)$；
• 将其进行普通消元使左半边的$A$为$I$，此时右半边即为$A^{-1}$。

此方法的原理为$A^{-1}\times \left(AI\right)=\left(I{A}^{-1}\right)$。

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 405, P = 1e9 + 7;
int n;
ll a[N][N << 1];
ll qpow(ll x, int y) {
	ll res = 1;
	while (y) {
		if (y & 1) {
			(res *= x) %= P;
		}
		(x *= x) %= P, y >>= 1;
	}
	return res;
}
void Gauss() {
	for (int i = 0; i < n; ++i) {
		int maxrid = i;
		for (int j = i + 1; j < n; ++j) {
			if (abs(a[j][i]) > abs(a[maxrid][i])) {
				maxrid = j;
			}
		}
		for (int j = 0; j < (n << 1); ++j) {
			swap(a[i][j], a[maxrid][j]);
		}
		if (!a[i][i]) {
			cout << "No Solution\n";
			exit(0);
		}
		ll k = qpow(a[i][i], P - 2);
		for (int j = 0; j < n; ++j) {
			if (i == j) {
				continue;
			}
			ll t = a[j][i] * k % P;
			for (int p = i; p < (n << 1); ++p) {
				a[j][p] = ((a[j][p] - a[i][p] * t) % P + P) % P;
			}
		}
		for (int j = 0; j < (n << 1); ++j) {
			(a[i][j] *= k) %= P;
		}
	}
}
int main() {
	ios :: sync_with_stdio(0);
	cin.tie(0), cout.tie(0);
	cin >> n;
	for (int i = 0; i < n; ++i) {
		for (int j = 0; j < n; ++j) {
			cin >> a[i][j];
			a[i][i + n] = 1;
		}
	}
	Gauss();
	for (int i = 0; i < n; ++i) {
		for (int j = n; j < (n << 1); ++j) {
			cout << a[i][j] << ' ';
		}
		cout << '\n';
	}
	return 0;
}

2. 特征多项式

这里需要用到部分前面高斯消元的知识。

特征多项式，记为$p_A(x)=\det (xIn-A)$，其中$A$为一个$n$阶矩阵，$I_n$为$n$阶单位矩阵，$p_A(x)$则为一个$n$次多项式。

我们先考虑上/下三角矩阵的解法，若$A$为上/下三角矩阵，易得：
$$p_A(x)=\prod _{i=1}^n\left(x-A_i{,}_i\right)$$
我们尝试在特征多项式不改变的情况下将矩阵变换成上/下三角矩阵。

2.1 高斯消元相似变换

有结论：若$P$可逆，$A$和$PA{P}^{-1}$的特征多项式相同。

证明：
$$\begin{array}{rl}\det (x{I}_n-PA{P}^{-1})& =\det (Px{I}_n{P}^{-1}-PA{P}^{-1})\\ & =\det (P\cdot P_{-1}\cdot (x{I}_n-A))\\ & =\det (P)\cdot \det (P^{-1})\cdot \det (x{I}_n-A)\\ & =\det (x{I}_n-A)\end{array}$$
我们称$A\to PA{P}^{-1}$为一次相似变换。

高斯消元可以让我们把$A$变成$PA$，那还有一个$P^{-1}$呢？其实就是对列进行同样变换。

那消元后是否可以将其化成上/下三角矩阵呢？不一定。对于第三种简单行变换，即$A_i=A_i+A_j\cdot k$，令其操作矩阵为$P_j{,}_i{,}_k$，则相似变换为$A\to P_i{,}_j{,}_{k}A{P}_i{,}_j{,}_{-k}$。高斯消元中通过$A\to P_i{,}_j{,}_{k}A$将$A_i{,}_j$消为$0$，但此处还有一个$P_i{,}_j{,}_{-k}$，可能使$A_i{,}_j$不为$0$，所以最后消出的矩阵形如：
$$\left(\begin{array}{ccccc}{a}_{1,1}& a_{1,2}& \dots & \dots & a_{1,n}\\ a_{2,1}& a_{2,2}& \dots & \dots & a_{2,n}\\ & \ddots & \ddots & \ddots & ⋮\\ & & \ddots & \ddots & a_{n-1,n}\\ & & & a_{n,n-1}& a_{n,n}\end{array}\right)$$
此类矩阵称为海森堡矩阵。

2.2 海森堡矩阵特征多项式

考虑递推。

• 状态定义：$f_i$为前$i$行$i$列关于$x$的特征多项式。

• 边界条件：$f_0=1$。

• 递推式：

  对于$1\le j\le i$，只有选$a_{1,1},a_{2,2},\dots ,a_{j-1,j-1}$、$a_{j,i}$、$a_{j+1,j},a_{j+2,j+1},\dots ,a_{i,i-1}$才会有贡献，则有递推式：
  $$f_i=(x-a_{i,i})f_{i-1}-\sum _{k=1}^{i-1}{a}_{i-k,i}\left(\prod _{j=i-k+1}^i{a}_{j,j-1}\right)f_{i-k-1}$$

这样我们就完成了$O(n^3)$求特征多项式。

#include <bits/stdc++.h>
using namespace std;
const int N = 505, P = 998244353;
int n, f[N], a[N][N];
int qpow(int x, int y) {
	int res = 1;
	while (y) {
		if (y & 1) {
			res = 1ll * res * x % P;
		}
		x = 1ll * x * x % P, y >>= 1;
	}
	return res;
}
void Gauss() {
	for (int i = 2; i <= n; ++i) {
		int j;
		for (j = i; j <= n; ++j) {
			if (a[j][i - 1]) {
				break;
			}
		}
		if (j > n) {
			continue;
		}
		if (j > i) {
			swap(a[i], a[j]);
			for (int k = 1; k <= n; ++k) {
				swap(a[k][j], a[k][i]);
			}
		}
		for (j = 1; j <= n; ++j) {
			a[j][i] = 1ll * a[j][i] * a[i][i - 1] % P;
		}
		int t = qpow(a[i][i - 1], P - 2);	// 逆元
		for (j = i - 1; j <= n; ++j) {
			a[i][j] = 1ll * a[i][j] * t % P;
		}
		for (j = i + 1; j <= n; ++j) {
			t = a[j][i - 1];
			for (int k = 1; k <= n; ++k) {
				a[k][i] = (a[k][i] + 1ll * a[k][j] * t) % P;
			}
			for (int k = i - 1; k <= n; ++k) {
				a[j][k] = (a[j][k] + 1ll * a[i][k] * (P - t)) % P;
			}
		}
	}
}
int g[N][N];
void solve() {
	memset(g, 0, sizeof g), g[0][0] = 1;
	for (int i = 1; i <= n; ++i) {
		int s = P - 1;
		for (int j = i; j; --j) {
			int t = 1ll * s * a[j][i] % P;
			for (int k = 0; k < j; ++k) {
				g[i][k] = (g[i][k] + 1ll * t * g[j - 1][k]) % P;
			}
			s = 1ll * s * a[j][j - 1] % P;
		}
		for (int k = 1; k <= i; ++k) {
			g[i][k] = (g[i][k] + g[i - 1][k - 1]) % P;
		}
	}
	memcpy(f, g[n], sizeof f);
}
int main() {
	ios :: sync_with_stdio(0);
	cin.tie(0), cout.tie(0);
	cin >> n;
	for (int i = 1; i <= n; ++i) {
		for (int j = 1; j <= n; ++j) {
			cin >> a[i][j];
		}	
	}
	Gauss(), solve();
	for (int i = 0; i <= n; ++i) {
		cout << f[i] << ' ';
	}
	return 0;
}

这次关于矩阵的研究就到这里，感谢您的观看！

p.s: 这篇中部分内容蒟蒻君也不是完全确定,有任何问题请大家指出，谢谢qwq