线性代数 - 线性方程组相关

一、高斯消元

Ex.1 P2455 线性方程组

1. 题面梗概

给定一个线性方程组，求解。

$$\begin{cases} a_{1, 1} x_1 + a_{1, 2} x_2 + \cdots + a_{1, n} x_n = b_1 \\ a_{2, 1} x_1 + a_{2, 2} x_2 + \cdots + a_{2, n} x_n = b_2 \\ \cdots \\ a_{n,1} x_1 + a_{n, 2} x_2 + \cdots + a_{n, n} x_n = b_n \\ \end{cases}$$

$1 \le n \le 50$。$\forall 1\le i, j \le n$，有 $\left | a_{i, j} \right | \le 100$，$\left |b_i \right | \le 300$。

2. 逻辑推导

图中的这个形式有点复杂，为了方便引入高斯消元法，我们拿一个稍微简单一点的三元一次方程组作为一个例子：

$$\begin{cases} x + y + 4z = 19 \\ 5x + 3y = 25 \\ 2y + 3z = 19 \\ \end{cases}$$

将未知变量对齐，如果未知变量不存在，也把它写上（ 系数为 0 ）：

$$\begin{cases} &x &+ &y &+ &4z &= &19 \\ &5x &+ &3y &+ &0z &= &25 \\ &0x &+ &2y &+ &3z &= &19 \\ \end{cases}$$

这样看起来就很优雅了，也许读者注意到了，看起来就像是一个 $3 \times 3$ 的矩阵乘上一个 $3$ 项的向量：

$$\begin{bmatrix} 1 & 1 & 4 \\ 5 & 3 & 0 \\ 0 & 2 & 3 \\ \end{bmatrix} \begin{bmatrix} x \\ y \\ z \\ \end{bmatrix} = \begin{bmatrix} 19 \\ 25 \\ 19 \\ \end{bmatrix}$$

将左边的系数矩阵称为系数矩阵 $A$，中间的未知量向量 $x$ 和右边的常量向量 $v$，这个线性方程组就可以简写为 $Ax = v$。

这种形式有一个很棒的几何性质：寻找一个位置的向量 $x$，经过矩阵 $A$ 的线性变换后，可以变为常数向量 $v$。

求解二元一次方程使用的主要方法是加减消元法和代入消元法：将方程组的其中一个方程的所有系数都乘上一个常数，加到另外一个方程中，方程仍然成立，得到一个未知量的解之后，带入到其他方程中，可以解出其他未知量的解。

首先，我们把系数矩阵提取出来，然后将常数向量写在系数矩阵的右边，成为增广矩阵，然后将第一行第一个数变为 $1$：

$$\begin{bmatrix} 1 & 1 & 4 & 19 \\ 5 & 3 & 0 & 25 \\ 0 & 2 & 3 & 19 \\ \end{bmatrix}$$

接下来模拟加减消元法。将第一行分别乘上若干倍，加到后面的每一行中，使得后面的每一行的第一个数字都为 $0$。这里每一行矩阵第一个非零元素对应的未知量被称为主元。

在本题中，第一次消元后如下：

$$\begin{bmatrix} 1 & 1 & 4 & 19 \\ 0 & -2 & -20 & -70 \\ 0 & 2 & 3 & 19 \\ \end{bmatrix}$$

接着处理第二行：

$$\begin{bmatrix} 1 & 1 & 4 & 19 \\ 0 & 1 & 10 & 35 \\ 0 & 2 & 3 & 19 \\ \end{bmatrix} \to \begin{bmatrix} 1 & 0 & -6 & -16 \\ 0 & 1 & 10 & 35 \\ 0 & 0 & -17 & -51 \\ \end{bmatrix}$$

最后处理第三行：

$$\begin{bmatrix} 1 & 0 & -6 & -16 \\ 0 & 1 & 10 & 35 \\ 0 & 0 & 1 & 3 \\ \end{bmatrix} \to \begin{bmatrix} 1 & 0 & 0 & 2 \\ 0 & 1 & 0 & 5 \\ 0 & 0 & 1 & 3 \\ \end{bmatrix}$$

发现去掉最右边的对应的常数向量，左边的系数矩阵是一个单位矩阵，直接使用这个单位矩阵乘上未知量向量，其结果就是常量向量，也就是说：

$$\begin{cases} x = 2 \\ y = 5 \\ z = 3 \\ \end{cases}$$

当然，如果出现本应该是主元的元素刚好是 $0$ 的情况，继续往下找，然后交换两行接着处理即可。

如果消到最后，发现了类似 $0x = 0$ 的情况，说明 $x_i$ 可以取任意值，有无穷多组解，如果出现了 $0x \neq 0$ 的情况，则说明无解。反之，就是有一组特定的解。

3. 代码实现

▶

AcCode

#include <bits/stdc++.h>
using namespace std;
using ll = long long;

const int maxn = 100 + 5;
const double eps = 1e-7;

double a[maxn][maxn];
int n;

// 1 唯一解 | 0 无穷解 | -1 无解
int gauss( )
{
    int curi = 1;
    for( int j = 1; j <= n; ++ j )
    {
	    // 找主元
        int p = curi;
        for( int i = curi + 1; i <= n; ++ i )
            if( fabs( a[i][j] ) > fabs( a[p][j] ) ) p = i;
        
        // 如果这一列没有找到非零元素
        if( fabs( a[p][j] ) < eps ) continue;

		// 交换两行
        swap( a[p], a[curi] );

		// 当前行同除，令主元系数化为 1 
        double div = a[curi][j];
        for( int k = j; k <= n + 1; ++ k ) a[curi][k] /= div;
        
        // 对其他行消元
        for( int i = 1; i <= n; ++ i )
        {
            if( i != curi && fabs( a[i][j] ) > eps )
            {
                double mul = a[i][j];
                for( int k = j; k <= n + 1; ++ k )
                    a[i][k] -=  a[curi][k] * mul;
            }
        }
        curi ++;
    }

    if( curi <= n )
    {
        for( int i = curi; i <= n; ++ i )
	        // 出现了 0x != 0
            if( fabs( a[i][n + 1] ) > eps ) return -1;
        // 因为 curi <= n 所以 存在 0x = 0
        return 0;
    }
    return 1;
}


int main( )
{
    ios::sync_with_stdio( 0 ); cin.tie( 0 );
    
    cin >> n;

    for( int i = 1; i <= n; ++ i )
        for( int j = 1; j <= n + 1; ++ j )
            cin >> a[i][j];

    int res = gauss( );
    
    if( res == 1 )
    {
        for( int i = 1; i <= n; ++ i )
            cout << "x" << i << "=" << fixed << setprecision( 2 ) << a[i][n + 1] << '\n';
    }
    else cout << res << '\n';
    return 0;
}

Ex.2 P2973 Driving Out the Piggies G

1. 题面梗概

给出一个 $n$ 个点，$m$ 条边的无向图，节点 $1$ 有一个炸弹，在单位时间内，有 $\frac{p}{q}$ 的概率在这个节点炸掉，有 $1 - \frac{p}{q}$ 的概率随机选一条出去的路到其他的节点上，去每个节点的概率相等。

问：最终炸弹在每个节点爆炸的概率。

$2 \le N \le 300, M \le 44850, 1 \le P \le Q \le 10^6$

2. 逻辑推导

首先有一个很显然的结论：在一个点爆炸的概率 = 这个点期望经过次数 $\times$ 每次炸弹爆炸的概率 $\frac{p}{q}$。

这里感性理解很容易，简单写一下解释：炸弹爆炸的次数的取值只有 $\{0, 1\}$，所以，炸弹期望爆炸次数与爆炸概率 在数值上相等。

期望 DP

设 $Ans_{u}$ 为炸弹在城市 $u$ 爆炸的概率，那么 $Ans_{u} = E_{u} \cdot p$，其中 $E_{u}$ 为炸弹在爆炸前经过城市 $u$ 的期望次数。

由于图中存在环路，炸弹可能在 $1 \to 2 \to 1$ 来回横跳，这意味着我们无法通过有限次数的计算，求出期望。

这种 “你中有我，我中有你” 的环状依赖关系（ 类马尔科夫链 ），我们就需要建立线性方程组。

简要写一下方程：（ $d_{u}$ 为节点 $u$ 的度数 ）

• 如果这个节点是 $1$，它由初始情况，加上从相邻城市转移过来的次数：

$$E_{1} = 1 + \sum_{(v, 1) \in E} \left( f_{v} \times ( 1 - p ) \times \frac{1}{d_{v}} \right)$$

• 如果这个节点不是 $1$，那么它只由相邻的城市转移过来：

$$E_{v} = \sum_{(u, v) \in E} \left( f_{u} \times (1 - p) \times \frac{1}{d_{u}} \right)$$

线性方程组

这样我们就可以得到 $N$ 个 $N$ 元一次方程，接着把方程转化为 $Ax = B$ 的标准形式，就可以通过高斯消元法，开心消消乐了。

3. 代码实现

▶

AcCode

#include <bits/stdc++.h>
using namespace std;
using ll = long long;

const int maxn = 300 + 5;
const double eps = 1e-7;

double a[maxn][maxn];
int g[maxn][maxn];
int deg[maxn];
int n, m, p, q;
double P;

void gauss( )
{
    int row = 1;
    for( int col = 1; col <= n; ++ col )
    {
        int pivot = row;
        for( int i = row + 1; i <= n; ++ i )
            if( fabs( a[i][col] ) > fabs( a[pivot][col] ) ) pivot = i;
        
        if( fabs( a[pivot][col] ) < eps ) continue;

        swap( a[pivot], a[col] );

        double div = a[row][col];
        for( int j = col; j <= n + 1; ++ j ) a[row][j] /= div;
        
        for( int i = 1; i <= n; ++ i )
        {
            if( i != row )
            {
                double fac = a[i][col];
                for( int j = col; j <= n + 1; ++ j )
                    a[i][j] -= fac * a[row][j];
            }
        }
        row ++;
    }
}

int main( )
{
    ios::sync_with_stdio( 0 ); cin.tie( 0 );
    
    cin >> n >> m >> p >> q;
    P = 1.0 * p / q;

    for( int i = 1; i <= m; ++ i )
    {
        int u, v; cin >> u >> v;
        if( g[u][v] ) continue;
        g[u][v] = g[v][u] = 1;
        deg[u] ++, deg[v] ++;
    }

    for( int i = 1; i <= n; ++ i )
    {
        a[i][i] = 1;
        for( int j = 1; j <= n; ++ j )
        {
            if( g[i][j] ) a[i][j] -= ( 1.0 - P ) / deg[j];
        }
    }
    a[1][n + 1] = 1.0;
    gauss( );

    cout << fixed << setprecision( 9 );

    for( int i = 1; i <= n; ++ i )
    {
        double ans = a[i][n + 1] * P;
        cout << ans << '\n';
    }

    return 0;
}

Ex.3 P2447 外星千足虫

1. 题面梗概

求解关于 $x_{i}$ 的方程组：

$$\begin{cases} \sum_{i=1}^{n}a_{1,i}x_{i} &=& b_{1} \bmod 2 \\ \sum_{i=1}^{n}a_{2,i}x_{i} &=& b_{2} \bmod 2 \\ & \vdots & \\ \sum_{i=1}^{n}a_{m,i}x_{i} &=& b_{m} \bmod 2 \\ \end{cases}$$

$N \le 1000, M \le 2000$，结果对 2 取模。

2. 逻辑推导

在模 $2$ 意义下，实际有用的值域只有 $\{0, 1\}$，那么我们再观察这个求和，$0 + 0 = 0, 0 + 1 = 1, 1 + 0 = 1, 1 + 1 = 0$，这不就是异或吗？

于是这就转化为了一个 异或线性方程组，但是高斯消元的复杂度为 $O(N^3)$，所以还得动点小脑筋，比如说 $Bitset$ 优化，可以把消元复杂度压缩到 $O( N / w )$（ w 为 机器字长，通常为 64 ）。

另外再补充一句，为什么要转化为 异或线性方程组？就是为了使用 Bitset 优化常数，从而 AC。

然后就是标准的高斯消元法了，只是把加法和乘法转化为异或来处理。从上往下消元的过程中，如果发现不满足条件的情况（多解），就直接输出。

Think is cheap, show me your code！

3. 代码实现

▶

AcCode

#include <bits/stdc++.h>
using namespace std;
using ll = long long;

const int maxn = 1e3 + 5;
const int maxm = 2e3 + 5;

bitset< maxn > a[maxm];
int n, m;
// 1 唯一解 | 0 无穷解
int gauss( int n, int m )
{
    int ans = 0;
    int curi = 1;
    for( int i = 1; i <= n; ++ i )
    {
        int pivot = curi;
        while( pivot <= m && a[pivot][i] == 0 ) pivot ++;
        
        if( pivot > m ) return 0;

        ans = max( ans, pivot );

        swap( a[pivot], a[curi] );

        for( int r = 1; r <= m; ++ r )
        {
            if( r != curi && a[r][i] ) a[r] ^= a[curi];
        }
        curi ++;
    }
    return ans;
}

int main( )
{
    ios::sync_with_stdio( 0 ); cin.tie( 0 );

    cin >> n >> m;
    for( int i = 1; i <= m; ++ i )
    {
        string s; int res;
        cin >> s >> res;
        for( int j = 1; j <= n; ++ j )
            if( s[j - 1] == '1' ) a[i][j] = 1;
        a[i][n + 1] = res;
    }

    int ans = gauss( n, m );
    if( ans == 0 ) 
    {
        cout << "Cannot Determine" << "\n";
        return 0;
    }

    cout << ans << '\n';
    for( int i = 1; i <= n; ++ i )
    {
        if( a[i][n + 1] ) cout << "?y7M#" << '\n';
        else cout << "Earth" << '\n';
    } 

    return 0;
}

二、行列式

待办ing

P7112 【模板】行列式求值