矩阵乘法笔记

· · 个人记录

从fibonacci数列开始说起

给定一个n，求fibonacci数列的第n项模一个数P，即F_n \mod P。

题目：洛谷1962

较朴素的方法：递推求数

int fibonacci(int n) {
    static int f[N];
    f[1] = 1;
    for (int i = 2; i <= n; ++i)
        f[i] = f[i - 1] + f[i - 2] % P;
    return f[n];
}

这样，我们可以发现其时间复杂度与空间复杂度都是O(n)，当然，我们可以发现每个数都只需要它之前的两个数，所以我们可以通过这一性质进行空间压缩，将空间复杂度降到O(1)。这个效率确实不错，但是当n的取值范围提升到10^{12}这个数量级时便无能为力了。

事实上，有一个方法可以以O(\log n)的效率将F_n算出。

将数列递推式转化为线性变换

我们可以发现，将空间压缩后的每个时刻状态，可以当做一个矢量：

A\_n = \left( \begin{array}{cc} F\_{n-1} & F\_n \end{array} \right)

在转化成下一个状态的过程中：

A\_{n+1} = \left( \begin{array}{cc} F\_{n} & F\_{n+1} \end{array} \right) = \left( \begin{array}{cc} A\_{n2} & A\_{n1}+A\_{n2} \end{array} \right)

我们便可以将其改写成一个矩阵变换，设

M = \left( \begin{array}{cc} 0 & 1 \\ 1 & 1 \end{array} \right)

则有

A_{n+1} = A_n M

根据归纳法，我们可以得知

A_n=A_1 M^{n-1}

进而有

F\_n=(A\_1 M^{n-1})\_2=(M^{n-1})\_{2,2}

那我们要问了，矩阵乘法的复杂度是O(r^3)的，这里r为2即乘方运算大致是8n次加法，而原本的求法是n次加法，这岂不是常数反而更大了？实际上，矩阵乘方和数字的普通乘方运算，都可以通过快速幂来进行加速。那么便可以通过O(\log n)的复杂度完成乘方运算。接下来，我将给出证明。

代码样例：

#include <cstdio>
#include <cstring>

using namespace std;

typedef long long ll;

struct matrix {
    int a[2][2];

    matrix()
    { memset(a, 0, sizeof(a)); }

    matrix(const matrix& x)
    { memcpy(a, x.a, sizeof(a)); }

    matrix operator*(const matrix& x) const;
};

const int P = 1e9+7;

int main() {
    matrix ans, m;
    ll n;
    scanf("%lld", &n);
    --n;
    ans.a[0][0] = ans.a[1][1] = 1;
    m.a[0][1] = m.a[1][0] = m.a[1][1] = 1;
    while (n) {
        if (n & 1)
            ans = ans * m;
        m = m * m;
        n >>= 1;
    }
    printf("%d\n", ans.a[1][1]);
    return 0;
}

matrix matrix::operator*(const matrix& x) const {
    matrix ret;
    for (int i = 0; i < 2; ++i)
        for (int j = 0; j < 2; ++j)
            for (int k = 0; k < 2; ++k)
                ret.a[i][j] = (((ll) a[i][k] * x.a[k][j]) + ret.a[i][j]) % P;
    return ret;
}

矩阵乘法的基本定理

矩阵乘法的结合律

设有矩阵A,B,C分别的大小为n \times m, m \times p, p \times q。求证(AB)C=A(BC)，进一步为(AB)C_{i,j}=A(BC)_{i,j}

根据矩阵乘法的定义，有

AB\_{i,j} = \sum\_{k=1}^m A\_{i,k} \times B\_{k,j}

故

```cpp \begin{eqnarray} (AB)C_{i,j} & = & \sum_{l=1}^p AB_{i,l} \times C_{l,j} \\ & = & \sum_{l=1}^p (\sum_{k=1}^m A_{i,k} \times B_{k,l}) \times C_{l,j} \\ & = & \sum_{l=1}^p \sum_{k=1}^m A_{i,k} \times B_{k,l} \times C_{l,j} \\ & = & \sum_{k=1}^m \sum_{l=1}^p A_{i,k} \times B_{k,l} \times C_{l,j} \\ & = & \sum_{k=1}^m (\sum_{l=1}^p B_{k,l} \times C_{l,j}) \times A_{i,k} \\ & = & \sum_{k=1}^m A_{i,k} \times BC_{k,j} \\ & = & A(BC)_{i,j} \end{eqnarray} ```

证毕。

矩阵快速幂的正确性

设方阵A(行数为r)和自然数n，可以将n分解为一个m位的二进制数

$$ \begin{eqnarray} A^n ```cpp & = & A^{\sum_{i=0}^m b_i 2^i} \\ & = & \prod_{i=0}^m A^{b_i 2^i} \\ & = & \prod_{b_i = 1} A^{2^i} \end{eqnarray} ``` $$ 从第一步到第二步的变换是不平凡的，它只有基于结合律才成立。因此，我们可以发现通过反复平方法可以以$O(r^3 \log n)$的效率求出$A^n$。 ## 重定义运算符 矩阵乘法中的运算涉及两种操作，即乘法和加法，我们在对结合律进行证明的过程中可以发现，我们依赖了：两种运算符各自符合交换律、结合律，乘法对加法符合分配率。形式化描述，即有运算$\oplus$加法与$\otimes$乘法，只要它们符合以下性质 $$ ```cpp \begin{eqnarray} x \oplus y & = & y \oplus x \\ (x \oplus y) \oplus z & = & x \oplus (y \oplus z) \\ x \otimes y & = & y \otimes x \\ (x \otimes y) \otimes z & = & x \otimes (y \otimes z) \\ x \otimes (y \oplus z) & = & (x \otimes y) \oplus (x \otimes z) \end{eqnarray} ``` $$ 那么将该运算替换掉原本矩阵乘法的加法和乘法，交换律依然成立，快速幂也仍然可以使用。比如说，我们知道带余加法和带余乘法仍然满足该性质，所以我们也可以置$\times_p \rightarrow \otimes, +_p \rightarrow \oplus$。 当然，如果还想要让这种重定义的矩阵可逆等，就需要该矩阵的元素所属集合与两种运算构成[域(field)](https://en.wikipedia.org/wiki/Field\_(mathematics))，也就是还要有逆元和单位元。 ### 矩阵快速幂解决最短路问题 我们尝试置$+ \rightarrow \otimes, \min \rightarrow \oplus$后，可以发现这仍然符合矩阵运算所需要的性质，而这时 $$ AB\_{i,j} = \bigoplus\_{k=1}^n A\_{i,k} \otimes B\_{k,j} = \min\_{k=1}^n (A\_{i,k} + B\_{k,j}) $$ （这个$min$的特别用法是我自己编的，嘻嘻） 我们惊奇地发现，这就等同于一次边$(i,j)$的松弛操作！于是我们对图的邻接矩阵$Adj$计算$n-1$次幂，就可以得出两两边的最短路了。其中如果$(u,v)$无边，则置$Adj_{u,v} \leftarrow +\infty$。计算$Adj^{n-1}$的时间复杂度为$O(n^3 \log n)$，虽然逊色于floyd算法的$O(n^3)$，但是这也是一个非常启发式的方法。