从高斯消元到行列式——线性代数进阶

naoliaok_lovely · 2023-10-09 14:30:29 · 算法·理论

注：本文中 A^{ab} 表示 A 的 ab 次方，A^{a\times b} 表示一个大小为 a\times b 的矩阵 A。

高斯消元

作用：解多元一次方程组。
高斯消元法：列出方程的增广矩阵并将其转变为上三角形式，然后带入消元求解。

code（朴素版）

int gauss()
{
    int i, j;
    for(i = 1, j = 1; j <= n; j++)
    {
        int t = i;
        for(int k = i + 1; k <= m; k++)
            if(abs(a[k][j]) > abs(a[t][j])) t = k;
        if(abs(a[t][j]) <= eps) continue;
        swap(a[i], a[t]);
        for(int k = n + 1; k >= j; k--)
            a[i][k] /= a[i][j];
        for(int k = i + 1; k <= m; k++)
            for(int l = n + 1; l >= j; l--)
                a[k][l] -= a[k][j] * a[i][l];
        i++;
    }
    for(int k = i; k <= m; k++)
        if(abs(a[k][n + 1]) > eps) return -1;
    if(i <= n) return 1;
    for(int j = n; j; j--)
        for(int i = 1; i < j; i++)
            a[i][n + 1] -= a[i][j] * a[j][n + 1];
    return 0;
}

注意事项：

1. 由于浮点误差，代码中使用了 eps 进行处理。并在每一次操作时，都把当前变量系数最大的那个拿上来，以减小误差。
2. 对于在此份代码中，唯一解返回 0，多解返回 1，无解返回 -1。对于有多解的方程，题目可能还会需要求方案数/最小字典序。方案数为 w^{cnt}，其中 w 为自由变元的值域，cnt 则为个数。而要求字典序的话，由于当中的自由变元可以随便取，而这份程序会跑出所有最靠后的自由变元，也就是说倒序处理、倒序输出即可。
3. gauss 不要写成 guass，不然会丢人。

当然了，高斯消元还有一个特殊的版本——解异或方程组。它的应用更广，最典型的是求解线性基。

Code（异或版）

int gauss()
{
    int i, j;
    for(i = 1, j = 1; j <= n; j++)
    {
        int t = i;
        for(int k = i + 1; k <= m; k++)
            if(a[k][j] > a[t][j]) t = k;
        if(!a[t][j]) continue;
        swap(a[i], a[t]);
        for(int k = i + 1; k <= m; k++)
            if(a[k][j]) a[k] ^= a[i];
        i++;
    }
    for(int k = i; k <= m; k++) 
        if(a[k][n + 1]) return -1;
    if(i <= n) return 1;
    for(int j = n; j; j--)
        for(int i = 1; i < j; i++)
            a[i][n + 1] = a[i][n + 1] ^ a[i][j] & a[j][n + 1];
    return 0;
}

注意事项：

1. 一般异或方程组的数据范围稍大，约为 500\sim1000，直接 O(N^3) 肯定 TLE。可以用 bitset 处理。
2. 经测验，bitset 似乎不支持单个元素的 ^= 运算，即 a[i][n + 1] ^= a[i][j] & a[j][n + 1] 会报错。

例题

模板：
P4035 [JSOI2008] 球形空间产生器
P3164 [CQOI2014] 和谐矩阵

进阶：
P2447 [SDOI2010] 外星千足虫
P3232 [HNOI2013] 游走（图上随机游走模型）

高斯消元求解稀疏（异或）线性方程组

NKOJ8530 手机游戏 超级版
一句话题意：n\times n 的 01 方阵，每次选择一个位置会使得它以及周围四个格子取反（四联通，总共改变五个格子）。输出一种方案使得全部变为 0，或者报告无解。n\le1000。

以现在的技术，这种题只能通过题目特性进行特殊处理。
在此题中，将第一行的点亮情况作为变量，将关系下放，在最后一行列出方程并求解即可。

事实上，求解稀疏方程组的通法是存在的，但是非常恶心。

高斯-约旦消元

一种特别的高斯消元方法。我们上面的朴素消元，是将原矩阵变为一个上三角矩阵；而高斯-约旦消元则是消成了一个单位矩阵。只需要将代码中对于 [i+1,n] 的消元改为 [1,i-1]\cap[i+1,n] 即可。事实上，高斯约旦消元理解之后会发现比朴素的高斯消元更简单，且更实用。

初等变化

初等变化总共有三个：

1. 交换两行。
2. 某行乘以某个非 0 系数。
3. 将某行乘以某个非 0 系数，再加到另一行。

显然，上面的高斯（约旦）消元可以用这三种变化实现。这三个变化被称之“初等”，一大原因是因为它们都可以写成矩阵乘法的形式。（每种变化对应的矩阵可以自行思考）于是，我们可以利用它求出逆矩阵。

逆矩阵

对于方阵 A^{n\times n}，如果存在 A^{-1} 满足

就称 A^{-1} 是 A 的逆矩阵。

考虑逆矩阵的计算方法。通过上面我们知道，可以通过高斯约旦消元将方阵消为单位矩阵，而每一步变化都可以写作矩阵乘法，所以可以得到 X_1X_2\dots X_kA=I，于是 A^{-1}=X_1X_2\dots X_k。
在实现的时候，我们往往不需要重新定义一个矩阵存答案，而是对原矩阵进行增广，即对于 [A\ |\ I] 进行高斯约旦消元，最后的结果就是 [I\ |\ A^{-1}]。

P4783 【模板】矩阵求逆

for(int i = 1; i <= n; i++)
{
    for(int j = i + 1; j <= n; j++)
        if(a[j][i])
        {
            swap(a[i], a[j]);
            break;
        }
    if(!a[i][i]) return puts("No Solution"), 0;
    int t = ksm(a[i][i], mod - 2);
    for(int j = i; j <= 2 * n; j++) a[i][j] = 1ll * a[i][j] * t % mod;
    for(int j = 1; j <= n; j++)
        if(j != i && a[j][i])
            for(int k = 2 * n; k >= i; k--)
                a[j][k] = (a[j][k] - 1ll * a[j][i] * a[i][k] % mod + mod) % mod;
}

虽然但是，逆矩阵在 OI 中的应用真的很少……

线性基

• 线性相关：对于向量 x_i，如果存在不全为 0 的序列 a，使得 a_1x_1+a_2x_2\dots a_nx_n=0，则称 x_1,x_2\dots x_n 是线性相关的。
• 线性基：从 x_1,x_2\dots x_n 中选出最大的子集，使得它们线性无关，则这组向量称为极大线性无关组，也称为基底、线性基。
• 向量空间：所有能由基底的线性组合表示的向量集合构成了一个向量空间。

一个重要的性质：对于 k 维空间中的向量，如果刚好存在 k 个线性无关的向量 x_1,x_2\dots x_k，则他们可以线性表示出所有 k 维空间中的向量，即他们的线性空间为整个 k 维空间。

比较抽象？别急，我们从另一个角度入手，先了解一下更特殊的版本。

异或空间线性基

定义：在只有异或运算意义下的线性基。换句话说，对于一个序列 A，如果对于序列 B 而言，他们进行异或运算可以得到的数集相同，且 B 线性无关，则我们称序列 B 是序列 A 的一组异或空间线性基。

通常情况下，我们求解异或空间线性基会采用高斯约旦消元，因为其算出的单位矩阵拥有更强的性质。显然的一件事是，一个序列对应的异或空间线性基可能不止一个。下面讨论的性质均以单位矩阵为例，朴素的上三角矩阵不一定满足。当然了，有时候上三角是比单位矩阵要好些的，所以也要灵活运用。
（这里的单位矩阵指的是，对于所有非自由变元的对应行列为单位矩阵，而自由变元的行列只需要满足上三角矩阵的限制）

7. 最大的异或值等于线性基中的所有数的异或。（仅对于单位矩阵）
8. 第 k 小的异或值等于其所有二进制位为 1 的线性基中数的异或和。（小的数对应低位，注意这里仍旧仅对于单位矩阵，且如果存在未能成功插入的 A 中的数要补 0）

void insert(LL x)
{
    for(int i = 60; ~i; i--)
        if(x >> i & 1)
        {
            if(!a[i])
            {
                a[i] = x;
                break;
            }
            x ^= a[i];
        }
}

说了这么多理论，接下来上一点实践：

P4570 [BJWC2011] 元素
如果求的不是最大的魔力，那么就是线性基的模板题。考虑如何满足魔力限制：
这里采用贪心的思路，将所有魔力进行排序，优先插入魔力值大的。可以证明这样做是理论上限，且达到了理论上限！（对于所有异或值为 0 的集合，至少得舍弃一个数。我们希望舍弃的是最小的那个，而这种方案刚好满足了我们的需求）

P4151 [WC2011] 最大XOR和路径
（有种教学关打完就打 Boss 的感觉）

这道题和上面的难度显然不在一个级别。对于这类题，一种常见的处理手法就是先随便找出一组解，再进行调整。因为这是“异或”，所以调整是再简单不过的事。
具体方法就是，先找出任意一条路径，考虑如何调整。注意到，从一条路径到另一条路径，相差的边刚好可以构成一个环。于是，我们只需要再找出所有的环，将环的异或和以及最开始的路径的异或和丢进线性基中，跑出来的最大异或值即为答案。
当然找出所有的环肯定那个无法通过，只保留一部分环即可，使得其他的环可以通过这些环异或出来。

CF1100F Ivan and Burgers
前缀和线性基！
题目要求的是区间异或的最大值，只需要求出区间的线性基即可。这类问题，我们以前是使用的前缀和的技巧；但是线性基是否可以求前缀和呢？当然可以！或者说，虽然形式是前缀和，但是本质上可能“动态后缀和”更为贴切。

在进行上面的 insert 操作时，我们额外记录一个 pos 数组，表示第 i 位的线性基中的数对应的位置。在从高位往低位枚举的过程中，如果枚举的数的位置比我们当前插入的位置更靠前，那么我们就让当前的位置去代替他。即：让每一个数位对应的数的位置尽可能靠后。关键函数大概长这样：

void insert(int a[], int pos[], int num, int id)
{
    for(int i = 20; ~i; i--)
        if(num >> i & 1)
            if(a[i])
            {
                if(pos[i] < id)
                    swap(a[i], num), swap(pos[i], id);
                num ^= a[i];
            }
            else
            {
                a[i] = num, pos[i] = id;
                break;
            }
}

这样做之后，如果要取出 l\sim r 的线性基，不再用 a_r-a_{l-1}，而是直接取出 1\sim r 的线性基，当中 pos_i<l 的元素丢掉即可。

习题：P3733 [HAOI2017] 八纵八横

实数范围内的线性基

回到最开始的问题：如何求解朴素的线性基，即实数范围内？事实上，我认为可能称之为高斯消元而非线性基更合适，不过部分的异或空间线性基的性质还是可以用的。

P3265 [JLOI2015] 装备购买
其实就是实数版本的 P4570 [BJWC2011] 元素。考虑从异或变为实数范围内的四则运算时候，上面的贪心规则还正确吗？容易证明依然正确，只不过构造线性基的过程更像高斯消元了而已。

矩阵的秩

定义：设矩阵 A\in R^{n\times m} 由列向量 \vec{a_1},\vec{a_2}\dots\vec{a_m} 构成，它们的极大线性无关组的数量称为矩阵的列秩，由它们的线性组合构成的向量集合称作该矩阵的列空间。显然，列秩就是列空间的维数。同理，行向量的极大线性无关组的数量称为行秩，线性组合构成的向量集合为行空间。
任意矩阵的行秩等于列秩，记作 r(A) 或 rank(A)。

看起来比较奇怪，举一个例子就懂了：

对于行向量，\vec{a_3}=2\vec{a_1}；对于列向量， \vec{b_3}=\frac{23}{11}\vec{b_1}+\frac{16}{11}\vec{b_2},\vec{b_4}=-2\vec{b_1}+6\vec{b_2}。故上述矩阵的秩 rank(A)=2。

显然，上面的初等变换不改变矩阵的秩。
高斯消元求行秩：\text{行秩}=\text{主元数}=n-\text{自由元数}。
思考：为何行秩等于列秩？（提示：可从矩阵的角度思考）

对于 A^{n\times m}，如果 rank(A)=n，则称 A 行满秩；如果 rank(A)=m，则称 A 列满秩。特别地，对于 n 阶方阵，若 rank(A)=n，则称 A 满秩。

一些性质：（证明略，作为读者的思维题还不错）

2. 任意矩阵乘以满秩矩阵，秩不改变。

一般来说，秩在 OI 中直接的运用比较少，主要是作为前置知识，推导一些更加复杂的数学算法。

行列式

我们先不给出行列式的定义，读者可以尝试一下能否通过下面的例子猜出来：

二阶行列式：

三阶行列式：

确实看上去没有什么规律，所以如果从来没有接触过行列式的人能够发现当中的规律，那么数感真的算不错了。
接下来我们引入真正的定义：

• 行列式(determinant)：对于 A^{n\times n}，其行列式 |A|=det(A)=\sum\limits_P(-1)^{\pi(P)}\prod\limits_{i=1}^na_{i,P_i}，其中 P 表示 1\sim n 的一个排列，\pi(P) 表示这个排列的逆序对个数。

可以带入上面的例子进行验证，顺便也可以理解这些字母的含义。

观察：行列式本质上就是一个 n 阶方阵，不过其等于一个具体的值，也就是说他是一个标量。右侧的式子展开之后共有 n!，表示 n 的全排列。由于自带一个 \pm1 的系数，所以他和容斥也有一些微妙的联系。

计算行列式

这玩意怎么计算？容易想到的做法基本都是阶乘级的，这让我们无法计算较大的行列式。能否找到一些快速计算的方法呢？

以上面的三阶行列式举例，我们考虑提取相同的项：

问题的规模变小了！似乎可以有效的优化复杂度。

• 余子式（也称主子式）：M_{i,j} 表示删去第 i 行第 j 列后的 n-1 阶子式。
• 代数余子式：A_{i,j}=(-1)^{i+j}M_{i,j}。

有了这些定义，上面的操作可以形式化的表示为：

前者我们叫做按第 i 行展开，后者叫做按第 j 列展开。

此方法的复杂度为多少呢？？？如果暴力递归，仍为 O(n!)，但是可以通过记忆化减少计算量，复杂度应该为其本质不同的 k 阶子式的个数，即 O\left(\sum(C_n^i)^2\right)，不好估计他的具体范围，不过显然有 2^n\le\sum(C_n^i)^2\le4^n，依然无法满足我们的需求。
另一方面，记忆化的过程本身也不简单，如何表示子式是一个问题。可以继续深入探究，不过此方法应该是走不通了，要考虑换一个方向。不过，这里提出的计算方式还是有用的，对于行列式的性质推导时会用到这种思想。

注意到：对于上三角方阵 A，对应的值 det(A)=\prod a_{i,i}。这是非常好的，我们考虑把朴素的方阵全部转化为上三角形式。提到上三角，我们的常用方法就是高斯消元，而高斯消元是由初等变化组成的。所以接下来，我们讨论初等变化对于行列式值的影响。

1. 交换两行。显然，每一个单项式的绝对值都应该是不变的，改变的只有前面的符号。可以证明，在一个排列里如果交换某两个数，逆序对的奇偶性必然改变！所以，每一个单项式都变为了他的相反数，即 det(A')=-det(A)。
2. 某行乘以某个非 0 系数。显然，det(A')=k\times det(A)
3. 将某行乘以某个非 0 系数，再加到另一行。带入原式可以发现（由乘法分配律可得），值不会改变，det(A')=det(A)。

至此，我们成功的将方阵变为了上三角形式。

质数模数版：

int get(int a[][N], int n)
{
    int ans = 1;
    for(int i = 1; i <= n; i++)
    {
        for(int j = i + 1; j <= n; j++)
            if(a[j][i])
            {
                swap(a[i], a[j]), ans = -ans;
                break;
            }
        if(!a[i][i]) return 0;
        int t = ksm(a[i][i], mod - 2);
        for(int j = i + 1; j <= n; j++)
            for(int k = n; k >= i; k--)
                a[j][k] = (a[j][k] - 1ll * a[j][i] * t % mod * a[i][k]) % mod;
        ans = 1ll * ans * a[i][i] % mod;
    }
    return (ans + mod) % mod;
}

非质数模数版：

int get(int a[][N], int n)
{
    int ans = 1;
    for(int i = 1; i <= n; i++)
    {
        for(int j = i + 1; j <= n; j++)
            if(a[j][i])
            {
                swap(a[i], a[j]), ans = -ans;
                break;
            }
        if(!a[i][i]) return 0;
        for(int j = i + 1; j <= n; j++)
            while(a[j][i])
            {
                int r = a[j][i] / a[i][i];
                for(int k = i; k <= n; k++)
                    a[j][k] = (a[j][k] - a[i][k] * r) % mod;
                if(a[j][i]) swap(a[i], a[j]), ans = -ans;
            }
        ans = 1ll * ans * a[i][i] % mod;
    }
    return (ans + mod) % mod;
}

实数版：

double get(double a[][N], int n)
{
    double ans = 1;
    for(int i = 1; i <= n; i++)
    {
        for(int j = i + 1; j <= n; j++)
            if(abs(a[j][i]) > eps)
            {
                swap(a[i], a[j]), ans = -ans;
                break;
            }
        if(abs(a[i][i]) < eps) return 0;
        for(int j = i + 1; j <= n; j++)
            if(abs(a[j][i]) > eps)
                for(int k = n; k >= i; k--)
                    a[j][k] -= a[j][i] / a[i][i] * a[i][k];
        ans *= a[i][i];
    }
    return ans;
}

其中上下两个没什么好说的，注意的是中间的辗转相除法。我们用辗转相除的方式避免了乘法逆元，正确性不难得到。注意的是时间复杂度，我们都知道辗转相除是自带一个 \log 的，那么复杂度会不会退化到 O(n^3\log m) 呢？答案是不会，因为每次由于是取 \gcd，所以这个数在随之变小，会变成原来的约数，故实际复杂度为 O(n^2(n+\log m))。

性质与相关定理

显然但很有用的一个结论是，如果存在自由变元，行列式的值为 0。于是，设方阵 A^{n\times n}，则以下描述等价：

3. 线性方程组 A\vec x=\vec b 有唯一解（其中 \vec b\ne\vec0）。

用行列式定义矩阵的秩：rank(A) 表示最大非 0 子式的阶数。

关于矩阵乘法，有个很好的性质是 det(A^{n\times n})det(B^{n\times n})=det(AB)。

证明：当 det(A)=0 显然成立，故假设 det(A)\ne0。
根据初等变化，设 A=T_1T_2\dots，其中 T_i 为初等变化对应的矩阵。

柯西-比内(Cauchy–Binet)公式

设 A^{n\times m},B^{m\times n} 满足 n\le m，则

其中 S 表示 \{1,2,\dots m\} 的 n 元子集，A_S 表示 A 取 S 中的列对应的 n 阶子式，B_S 表示 B 取 S 中的行对应的 n 阶子式。

也比较好理解，特别的，当 n=m 时，就变为了上面的 det(A)det(B)=det(AB)。

范德蒙德(Vandermonde)行列式

证明？使用数学归纳法，按第一行展开后化简即可。

几何意义

二阶行列式的绝对值：行/列向量构成的平行四边形的面积。
三阶行列式的绝对值：行/列向量围成的平行六面体的体积。
规律？