矩阵和线性方程组

矩阵和线性方程组

矩阵运算，本质上是一个线性变换，将一个列向量 $\bm{x}$ 变换成另一个列向量 $\bm{b}$
$\bm{A}\bm{x} = \bm{b}$

逆矩阵 矩阵求逆在算法中很少使用，但思想是比较直观的
对于 $\bm{Ax} = \bm{b}$，$\bm{A}^{-1} \bm{Ax} = \bm{A}^{-1} \bm{b}$，$\bm{A}^{-1} \bm{A} = \bm{I}$
那么有 $\bm{x} = \bm{A}^{-1} \bm{b}$

矩阵可逆，当且仅当所有行向量线性无关，可化为行阶梯形矩阵

高斯消元

按照一般代数学知识，作增广矩阵后，化为行阶梯形矩阵 (REF)
$\forall i \in [0, n-1]$，$A_{i, n}$ 即为 $x_i$ 的值

• 列主元
  当前正在处理第 $i$ 行，用第 $i$ 行依次消去第 $[i+1\cdots n-1]$ 行
  需要保证的是 $A_{i, i} \neq 0$，为了提高数值稳定性
  找到一个 $r > i$ 并且绝对值最大的 $A_{r, i}$，然后交换第 $r$ 行与第 $i$ 行

• 消元，用第 $i$ 行依次消去第 $k \in [i+1 \cdots n-1] \quad (k > i)$ 行
  第 $k$ 行所有列元素 $\forall j \in [0, n-1]$

  $$\begin{aligned} A_{k, j} &-= f \cdot A_{i, j}, \ f = \frac{A_{k, i}}{A_{i, i}} \\ A_{k, j} &-= \frac{A_{k, i}}{A_{i, i}} \cdot A_{i, j} \end{aligned}$$

  为了保证精度，上式用右边的 $A_{k, i}$ 来更新 $A_{k, j}$，要保证 $A_{k, i}$ 不会过早被破坏
  只要从 $\forall \ \textbf{for} \ j = n \to i$ 逆序枚举即可

  for (int j = n; j >= i; j--) {
    for (int k = i+1; k < n; k++) {
        A[k][j] -= A[k][i] / A[i][i] * A[i][j];
    }
  }

• 回代过程，化为若而当标准型
  在若而当标准型中，可以直接得到 $x_i = A_{i, n} / A_{i, i}$
  化为若而当标准型，对于第 $i$ 行，消去第 $j \in [i+1 \cdots n-1]$ 列

  • 如果令 $\textbf{for} \ \forall i = n-1 \to 0$ 开始化为若而当标准型，要消去第 $j \in [i+1 \cdots n-1]$ 列
    于此同时 $j > i$ 的列已经满足若而当标准型，并且 $A_{j, n} \leftarrow A_{j, n} / A_{j, j}$
    对于 $j > i$，此时 $x_j = A_{j, n}$
  • 对于第 $i$ 行，$j > i, \ x_j$ 的系数为 $A_{i, j}$，代入移项
    $\textbf{for} \ \forall j \in [i+1\cdots n-1]$
    $\quad A_{i, n} -= A_{j, n} \cdot A_{i, j}$
    之后令 $A_{i, n} \leftarrow A_{i, n} / A_{i, i}$

void gauss(double A[][maxn]) {
    int i, r;
    for (i = 0; i < n; i++) {
        // swap
        r = i;
        for (int j = i+1; j < n; j++) {
            if (fabs(A[j][i]) > fabs(A[r][i])) r = j;
        }
        if (r != i) for (int j = 0; j <= n; j++) swap(A[r][j], A[i][j]);

        // eliminate
        for (int j = n; j >= i; j--) {
            for (int k = i+1; k < n; k++)
                A[k][j] -= A[k][i] / A[i][i] * A[i][j];
        }

        // Jordan solution
        for (int i = n-1; i >= 0; i--) {
            for (int j = i+1; j < n; j--) A[i][n] -= A[i][j] * A[j][n];
            A[i][n] /= A[i][i];
        }
    }
}

矩阵工具的应用

矩阵快速幂

Recurrences

首先根据递推关系构造矩阵

$$\begin{aligned} f(n) &= a_1f(n-1) + &&a_2f(n-2) + \cdots + a_df(n-d) \\ f(n-1) &= &&a_1f(n-2) + \cdots + a_{d-1}f(n-d) + a_df(n-(d+1)) \end{aligned}$$

加上恒等式 $f(n-1) = f(n-1), f(n-2) = f(n-2) \cdots$

可以构造矩阵如下

$$\bm{F}_n = \begin{pmatrix} f(n) \\ f(n-1) \\ f(n-2) \\ \vdots \\ f(n-d) \end{pmatrix} \quad \bm{F}_{n-1} = \begin{pmatrix} f(n-1) \\ f(n-2) \\ \vdots \\ f(n-d) \end{pmatrix}$$

$$\bm{F}_n = \bm{A}\bm{F}_{n-1} = \begin{pmatrix} a_1 & a_2 & a_3 &\cdots & a_d \\ 1 & 0 & 0 &\cdots & \vdots \\ 0 & 1 & 0 &\cdots & \vdots \\ \vdots &\vdots &\vdots & \ddots & \vdots \\ 0 & 0 & 0 &\cdots &1 \end{pmatrix} \begin{pmatrix} f(n-1) \\ f(n-2) \\ \vdots \\ f(n-d) \end{pmatrix}$$

其中 $\bm{F}_n((d+1) \times 1), \bm{F}_{n-1}(d \times 1), \bm{A}((d+1) \times d)$，经过迭代后

$$\bm{F}(d) = \begin{pmatrix} f(d) \\ f(d-1) \\ \vdots \\ f(1) \end{pmatrix} \quad \bm{F}(n) = \bm{A}^{n-d} \bm{F}(d)$$

要求 $\bm{F}(n)$，进行矩阵幂运算后，执行矩阵乘法，得到结果矩阵为 $\bm{F}$，$\bm{F}_{0, 0}$ 即为所求

特殊情况，对于 $d \geqslant n$ 的情况如何处理？

$$\bm{F} = \bm{A}\bm{F}_0 \Longrightarrow \begin{pmatrix} f(d) \\ f(d-1) \\ \vdots \\ f(n) \\ \vdots \end{pmatrix} = \bm{A}\begin{pmatrix} f(d) \\ f(d-1) \\ \vdots \\ f(1) \end{pmatrix}$$

由此 $\bm{F}$ 的 $f(d-0), f(d-1), \cdots f(d-k)$ 与 $\bm{F}_0$ 的第 $0, 1, \cdots k$ 行对应
$d-k = n, \ k = d-n$，由此 $\bm{F}_0(d-n, 0)$ 即为答案

const int maxn = 15 + 5;
typedef vector<vector<ll> > Matrix;
int d, n, mod;

Matrix matrix_mul(const Matrix &A, const Matrix &B) {
    int _n = A.size(), _m = B[0].size();
    Matrix res(_n, vector<ll> (_m, 0));

    for (int i = 0; i < _n; i++) {
        for (int j = 0; j < _m; j++)
            for (int k = 0; k < (int)A[0].size(); k++)
                res[i][j] = (res[i][j] + A[i][k] * B[k][j]) % mod;
    }

    return res;
}

Matrix matrix_pow(const Matrix &A, int exp) {
    int sz = A.size();
    Matrix I(sz, vector<ll>(sz, 0));
    for (int i = 0; i < sz; i++) I[i][i] = 1;

    Matrix res = I, B = A;
    while (exp) {
        if (exp & 1) res = matrix_mul(res, B);
        B = matrix_mul(B, B);
        exp >>= 1;
    }
    return res;
}

int main() {
    freopen("input.txt", "r", stdin);
    while (scanf("%d%d%d", &d, &n, &mod) == 3 && d) {
        // init
        // get data
        Matrix A(d, vector<ll>(d, 0));
        Matrix F(d, vector<ll>(1, 0));

        for (int i = 0; i < d; i++) {
            int x;
            scanf("%d", &x);
            A[0][i] = x;
        }
        for (int i = 1; i < d; i++) A[i][i-1] = 1;
        for (int i = d-1; i >= 0; i--) scanf("%d", &F[i][0]);

        // solve
        if (n <= d) {
            printf("%lld\n", F[d-n][0]);
            continue;
        }
        Matrix res = matrix_mul(matrix_pow(A, n-d), F);
        printf("%lld\n", res[0][0]);
    }
}

循环矩阵和变换

Cellular Automation
构造变换矩阵 $\bm{A}$，$\bm{A}(n \times n)$
对于第 $i$ 行，$\bm{A}_{i, i} = 1$，对于任意的 $\forall j$，$\min(|j-i|, n-|j-i|) \leqslant d$
满足 $\bm{A}_{i, j} = 1$，举个例子，如果 $d = 1$，那么有

$$\bm{A} = \begin{pmatrix} 1 & 1 & 0 & \cdots & 0 & 1 \\ 1 & 1 & 1 & \cdots & 0 & 0 \\ 0 & 1 & 1 & 1 & 0 & \cdots \\ & \vdots \\ 1 & 1 & 0 & \cdots & 0 & 1 \end{pmatrix}$$

下一行都是上一行的循环右移

$\bm{v}_0$为初始矩阵，经过 $k$ 次变换后，$\bm{v} = \bm{A}^k \bm{v}_0$
循环矩阵只要知道第一行，可以推出任意一行，于是矩阵乘法可以优化，每次只要求出变换后的第一行即可

$$\begin{aligned} \textbf{for} \ \forall i \in [0, n), \ &\forall j :\\ &A_{0, i} = \sum_{j = 0}^{n-1} A_{0, j} \cdot B_{j, i}, \ \ B_{j, i} = B_0[(i - j+n) \bmod n] \end{aligned}$$

$B_0$ 存储第一行的初始向量，$\bm{B}_{j, i}$ 是 $\bm{B}_{0, i}$ 循环右移 $j$ 位形成的
所以 $\bm{B}_{j, i} = B_0[(i-j+n) \bmod n]$，这样只需要 $O(n^2)$ 时间复杂度计算矩阵乘法

typedef vector<vector<ll> > Matrix;
int n, m, d, k;
Matrix A;
vector<ll> V;

Matrix mul(const Matrix &A, const Matrix &B) {
    int _n = A.size(), _m = B[0].size();
    Matrix C(_n, vector<ll> (_m, 0));
    for (int i = 0; i < _m; i++) {
        for (int j = 0; j < _m; j++)
            C[0][i] = (C[0][i] + A[0][j] * B[j][i]) % m;
    }
    for (int i = 1; i < _n; i++) {
        for (int j = 0; j < _m; j++)
            C[i][j] = C[i-1][(j-1+n)%n];
    }
    return C;
}

Matrix pow(const Matrix &A, int exp) {
    Matrix res(A.size(), vector<ll> (A[0].size(), 0));
    for (int i = 0; i < res.size(); i++) res[i][i] = 1;
    Matrix B = A;
    while (exp) {
        if (exp & 1) res = mul(res, B);
        B = mul(B, B);
        exp >>= 1;
    }
    return res;
}

void clear() {
    V = vector<ll> (n, 0);
    A = Matrix(n, vector<ll> (n, 0));
}

void calc(const Matrix &A) {
    vector<ll> res(n, 0);
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++)
            res[i] = (res[i] + A[i][j] * V[j] % m) % m;
    }
    printf("%lld", res[0]);
    for (int i = 1; i < res.size(); i++) printf(" %lld", res[i]);
    printf("\n");
}

int main() {
    freopen("cell.in", "r", stdin);
    freopen("cell.out", "w", stdout);
    while (scanf("%d%d%d%d", &n, &m, &d, &k) == 4) {
        // init
        clear();

        // data
        for (int i = 0; i < n; i++) scanf("%lld", &V[i]);
        for (int i = 0; i < n; i++) A[i][i] = 1;
        for (int i = 0; i < n; i++) {
            for (int j = 1; j <= d; j++) A[i][(i+j+n)%n] = A[i][(i-j+n)%n] = 1;
        }

        // calc
        Matrix T = pow(A, k);
        calc(T);
    }
}

马尔可夫问题

Gauss-Jordan消元
在 Gauss 消元法的基础上，可以省略掉回代过程，此时需要把矩阵化为对角阵

化为 Jordan 标准型

• $\textbf{for} \ \forall i = 1 \to n$
  $\quad$ 在 $[i+1, n]$ 中选出第 $r$ 行，使得 $A_{r,i}$ 最大，$\bold{swap}(A_{r, \cdots}, A_{i, \cdots})$
  $\quad$ 基于第 $i$ 行，消去所有 $k \neq i$ 的行
  $\quad \ \textbf{for} \ \forall k \neq i, \quad \textbf{for} \ \forall j = n \ \textbf{downto} \ i$
  $\quad\quad\quad$ $A(k, j) -= A(k, i) / A(i, i) \cdot A(i, j)$
  $\quad$（$j < i$ 的列 $A_{i, j} = 0$，不影响最终结果）

UVA10828
如果把节点 $u$ 的出度记为 $d_u$，那么对于节点 $i$，其期望执行的次数记为 $x_i$

$$x_i = \sum_{u \in \text{pre}(i)} \frac{x_u}{d_u}$$

$\text{pre}(i)$ 为 $i$ 的前驱节点

由此可以根据有向图列出期望线性方程，求解出 $x_u$，即为节点 $u$ 的期望执行次数
但这样做会出现一些问题，有可能求出来的 $x_u$ 是特殊的
可能 $u$ 期望执行次数为 $0$，也有可能为 $\infty$

解决的方法很简单，只要使用Gauss-Jordan来避免回代过程，最后
$x_i = A(i, n) / A(i, i)$

• 如果消元之后，$A(i, i) = 0 \ \textbf{and} \ A(i, n) \neq 0$，此时 $i$ 节点 $x_i = +\infty$
  $A(i, i) = 0 \ \textbf{and} \ A(i, n) = 0$，那么 $x_i = 0$
  若 $x_i$ 有限并且唯一，那么第 $i$ 行除了 $A(i, i), \ A(i, n)$ 以外所有数均为 $0$

typedef vector<vector<double> > Matrix;
const int maxn = 100 + 10;
const double eps = 1e-8;
int n, deg[maxn], Inf[maxn];
vector<int> pre[maxn];
Matrix A;

void build() {
    for (int i = 0; i < n; i++) {
        A[i][i] = 1;
        for (auto u : pre[i]) A[i][u] -= 1.0 / deg[u];
        if (i == 0) A[i][n] = 1;
    }
}

void gauss_jordan() {
    int i, r;
    for (i = 0; i < n; i++) {
        r = i;
        for (int j = i+1; j < n; j++) {
            if (fabs(A[j][i]) > fabs(A[r][i])) r = j;
        }
        if (fabs(A[r][i]) < eps) continue;
        if (r != i) for (int j = 0; j <= n; j++) swap(A[r][j], A[i][j]);

        // elimination
        for (int k = 0; k < n; k++) if (k != i) {
                for (int j = n; j >= i; j--) A[k][j] -= A[k][i] / A[i][i] * A[i][j];
            }
    }
}

void solve() {
    memset(Inf, 0, sizeof Inf);
    for (int i = n-1; i >= 0; i--) {
        if (fabs(A[i][i]) < eps && fabs(A[i][n]) > eps) Inf[i] = 1;
        for (int j = i+1; j < n; j++)
            if (fabs(A[i][j]) > eps && Inf[j] == 1) Inf[i] = 1;
    }

    int q, u;
    scanf("%d", &q);
    while (q--) {
        scanf("%d", &u);
        u--;
        double res = 0;
        if (Inf[u]) printf("infinity\n");
        else {
            res = fabs(A[u][u]) < eps ? 0.0 : A[u][n] / A[u][u];
            printf("%.3lf\n", res);
        }
    }
}

void init() {
    memset(deg, 0, sizeof deg);
    for (int i = 0; i < maxn; i++) pre[i].clear();
    A = Matrix(n, vector<double>(n+1, 0));
}

int main() {
    freopen("input.txt", "r", stdin);
    int cas = 0;
    while (scanf("%d", &n) == 1 && n) {
        // init
        printf("Case #%d:\n", ++cas);
        init();

        // get data
        int a, b;
        while (scanf("%d%d", &a, &b) == 2 && a) {
            // a->b
            a--, b--;
            deg[a]++, pre[b].push_back(a);
        }
        // get matrix
        build();

        // gauss-jordan
        gauss_jordan();

        // solve
        solve();
    }
}

异或方程组的高斯消元

可以将该算法拓展到一般的代数系统中，比如线性方程组的代数系统扩展为模 $p$ 的剩余系 $\mathbb{Z}_{p}$

Square
题意大致是给定 $n$ 个数，从中选出一个或者多个，使得所选出的整数乘积是完全平方数

不含有大于 $500$ 的素因子，很容易想到唯一分解定理
以样例来看，对于 $\{4, 6, 10, 15\}$，唯一分解所用到的素因子为 $\{2, 3, 5\}$
观察幂指，根据幂指来构造向量如下

$$\begin{aligned} 4 = 2^2 \cdot 3^0 \cdot 5^0 \quad \quad &(2, 0, 0) \\ 6 = 2^1 \cdot 3^1 \cdot 5^0 \quad \quad &(1, 1, 0) \\ 10 = 2^1 \cdot 3^0 \cdot 5^1 \quad \quad &(1, 0, 1) \\ 15 = 2^0 \cdot 3^1 \cdot 5^1 \quad \quad &(0, 1, 1) \end{aligned}$$

原数组为 $\{a_1, a_2, \cdots a_n\}$，用 $x_i$ 表示选 or 不选 $a_i$ 这个数
状态数组为 $\{x_1, x_2, \cdots x_n\}$

那么上述问题实际上是要考虑

$$(x_1, x_2, x_3, x_4) \begin{pmatrix} 2 & 0 & 0 \\ 1 & 1 & 0 \\ 1 & 0 & 1 \\ 0 & 1 & 1 \end{pmatrix}$$

要为完全平方数，必须满足

$$\begin{cases} 2x_1 + x_2 + x_3 + 0x_4 \equiv 0 (\bmod 2) \\ 0x_1 + x_2 + 0x_3 + x_4 \equiv 0 (\bmod 2) \\ 0x_1 + 0x_2 + x_3 + x_4 \equiv 0 (\bmod 2) \end{cases}$$

上述问题实际上是 $\mathbb{Z}_2$ 剩余系中的矩阵运算
注意到形如 $2x_i$ 的数，总是能够写成完全平方数，也就是说，$2x_i \equiv 0(\bmod 2)$
偶数在模 $2$ 剩余系中的运算结果总是 $0$，这启发我们用异或 $\oplus$ 运算

$$\begin{cases} x_2 \oplus x_3 = 0 \\ x_2 \oplus x_4 = 0 \\ x_3 \oplus x_4 = 0 \end{cases}$$

于是可以这样构造算法，为了符合一般矩阵运算的写法

$$(x_1, x_2, x_3, x_4) \begin{pmatrix} 2 & 0 & 0 \\ 1 & 1 & 0 \\ 1 & 0 & 1 \\ 0 & 1 & 1 \end{pmatrix}$$

取一个转置变为

$$\bm{Ax} = \begin{pmatrix} 2 & 1 & 1 & 0 \\ 0 & 1 & 0 & 1 \\ 0 & 0 & 1 & 1 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \end{pmatrix} \equiv \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} (\bmod 2)$$

具体来说

• 对 $\{a_1, a_2 \cdots a_n\}$ 唯一分解，假设所用到的素数为 $\{p_1, p_2, \cdots p_m\}$
  构造 $m \times n$ 的矩阵 $\bm{A}$，如果 $a_i$ 含有素因子 $p_j$，那么令 $A(j,i) \oplus = 1$
  注意如果 $a_i$ 中 $p_j$ 的幂指为 $k$，需要 $A(j, i) \oplus = 1$ 操作 $k$ 次
• 对异或矩阵执行高斯消元，并且求出矩阵的秩为 $r$，矩阵的自由变量为 $n-r$
  每一个自由变量有 $2$ 种取值，选 or 不选，所以答案为 $2^{n-r}$

异或方程组的高斯消元比较简单

• 假设当前正在处理第 $i$ 行第 $j$ 列，首先要找到 $A(r, j) \neq 0$ 的第 $r$ 行
  并且将第 $r$ 行与第 $i$ 行交换
• 消元之后第 $i$ 行的第一个非 $0$ 列为第 $j$ 列
  并且对于所有 $u > i$ 的行，它的第 $j$ 列都为 $0$
  对于 $u \in [i+1, n]$ 的行，用第 $i$ 行来消去它
  即对于列 $k \in [i, n]$，对应地令 $A(u, k) \oplus = A(i, k)$

const int maxn = 500 + 10;
typedef vector<vector<int> > Matrix;
int tot = 0, n;
Matrix A;

vector<int> primes;

void get_primes(int m) {
    primes.clear();

    bool st[maxn];
    memset(st, 0, sizeof st);

    for (int i = 2; i <= m; i++) {
        if (st[i]) continue;
        primes.push_back(i);
        for (int j = i; j <= m; j += i) st[j] = 1;
    }
}

int gauss_jordan(Matrix &A, int m, int n) {
    // m * n
    int i = 0, j = 0;
    for (; i < m && j <= n; j++) {
        // find r, A(r, j) != 0
        int r = i;
        for (int k = i; k < m; k++) if (A[k][j]) {
            r = k; break;
        }
        if (A[r][j] == 0) continue;
        if (r != i) for (int k = 0; k <= n; k++) swap(A[r][k], A[i][k]);

        // elimination
        for (int u = i+1; u < m; u++) if (A[u][j]) {
            for (int k = i; k <= n; k++) A[u][k] ^= A[i][k];
        }
        i++;
    }
    return i;
}

void init() {
    tot = 0;
    A = Matrix(maxn, vector<int>(maxn, 0));
}

int main() {
    freopen("input.txt", "r", stdin);
    int cas;
    cin >> cas;
    get_primes(500);

    while (cas--) {
        init();
        
        scanf("%d", &n);
        for (int i = 0; i < n; i++) {
            ll x; scanf("%lld", &x);

            for (int j = 0; j < primes.size(); j++) {
                while (x % primes[j] == 0) {
                    x /= primes[j], A[j][i] ^= 1;
                    tot = max(tot, j);
                }
            }
        }

        int r = gauss_jordan(A, tot+1, n);
        ll res = (1ll << (n-r)) - 1;
        printf("%lld\n", res);
    }
}