算法动态规划

矩阵快速幂及其动态规划优化中的应用

字数统计: 2.7k字   |   阅读时长: 12分
, , , , , ,

摘要: 本文介绍矩阵快速幂的算法原理与代码模板，并解决一个例题

【对算法，数学，计算机感兴趣的同学，欢迎关注我哈，阅读更多原创文章】
我的网站：潮汐朝夕的生活实验室
我的公众号：算法题刷刷
我的知乎：潮汐朝夕
我的github：FennelDumplings
我的leetcode：FennelDumplings

在文章 多米诺与托米诺骨牌平铺 中，我们用矩阵快速幂优化 DP 的方式解决了一个计数问题。对于状态转移方程为线性方程的时候，经常可以用这种矩阵快速幂的方式进行优化。

本文我们就介绍一下矩阵快速幂的算法原理以及代码模板，然后解决另一个矩阵快速幂优化 DP 的题目，该题的递推关系非常简单，用母函数求通项公式并不复杂，因此也写了一下用母函数求通项公式的过程，对于其他的递推关系，大致过程也类似，只是这个过程涉及手算，对于编程的问题一般不用。

矩阵快速幂的算法原理

矩阵快速幂与快速幂的思想是一致的: 求矩阵的幂 $A^{n}$ 的时候，先不求 $A^{n}$，而是先求 $A^{\frac{n}{2}}$，然后先不求 $A^{\frac{n}{2}}$，而是先求 $A^{\frac{n}{4}}$，…

以上思路同样可以用二进制分解加位掩码来实现，是当某个位掩码为1表示该位需要的时候，做的是快速幂做的是乘法，矩阵快速幂做的是矩阵乘法.

1
2
3
4
5
6
7
while(n)
{
    if(n & 1)
        ans = ans * A; // 这一步是矩阵的乘法
    A = A * A;
    n >>= 1;
}

矩阵快速幂的代码模板

实现矩阵快速幂，可以先写一个 struct Matrix 表示矩阵，之后的乘法，求幂都是针对这个类的运算。

其中用一个二维数组做数据成员，一个 init() 方法，将矩阵初始化为单位阵(相当于 $A^{0}$) 的结果。

然后重载 * 和 ^ 分别表示矩阵乘法和矩阵快速幂. 代码如下。

M 表示方阵的行数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
using ll = long long;
const int M = 2;

struct Matrix
{
    int a[M][M];
    Matrix()
    {
        memset(a, 0, sizeof(a));
    }

    void init() // 复位为单位阵
    {
        a[0][0] = a[1][1] = 1;
        a[0][1] = a[1][1] = 0;
    }

    Matrix operator*(const Matrix& B) const
    {
        Matrix ans;
        for(int i = 0; i < M; ++i)
            for(int j = 0; j < M; ++j)
                for(int k = 0; k < M; ++k)
                    ans.a[i][j] += a[i][k] * B.a[k][j];
        return ans;
    }

    Matrix operator^(int n) const
    {
        Matrix ans;
        ans.init();
        Matrix A = *this; // 拷贝一个出来用于自乘
        while(n)
        {
            if(n & 1)
                ans = ans * A;
            A = A * A;
            n >>= 1;
        }
        return ans;
    }

    void show()
    {
        cout << "matrix" << endl;
        for(int i = 0; i < M; ++i)
        {
            for(int j = 0; j < M; ++j)
                cout << a[i][j] << " ";
            cout << endl;
        }
    }
};

int main()
{
    Matrix A;
    A.a[0][0] = 0;
    A.a[0][1] = A.a[1][0] = A.a[1][1] = 1;
    A.show();

    Matrix F;
    F.a[0][0] = F.a[1][0] = 1;
    F.a[0][1] = F.a[1][1] = 0;
    F.show();

    int n = 8;
    Matrix ans = (A ^ n) * F;
    ans.show();
    return 0;
}

矩阵快速幂的适用问题

矩阵快速幂主要用于线性的递推型计数问题，以及一些动态规划中状态转移方程是线性递推关系的时候。

有一些计数问题，通过对 n 比较小的情况进行手动推导，发现规律，可以猜想出递推关系，如果这个递推关系是线性，那么它可以转换成矩阵求幂问题，进而可以用矩阵快速幂加速。

例如斐波那契问题：

我们会思考是否可以构造一个矩阵 A，使得:

通过一些推导，可以得到 $A = \begin{bmatrix} 0 & 1 \\ 1 & 1 \end{bmatrix}$，此关系可以继续推导下去，直到推出 $\begin{bmatrix} f(n) \\ f(n + 1) \end{bmatrix}$。

最终可以得到:

推出以上公式之后就可以用套矩阵快速幂的算法求 $A^{n}$ 了，此过程时间复杂度为 $O(M^{3}\log N)$ ，M 为方阵的行数。

其它的线性递推关系均可以这样做: 先找规律得到递推关系，然后将递推关系变为矩阵的形式，套矩阵快速幂。

例如: $f(n) = 3f(n - 1) + 6f(n - 2) - 7f(n - 3)$，的递推关系，可以变为矩阵形式：

线性递推关系是组合数学中的重要内容，涉及到母函数，特征多项式，待定系数法等，一般都可以通过母函数或者特征多项式求出通项公式，例如斐波那契数列的通项公式。

有一些递推关系是非线性的，就没法用矩阵快速幂的方法了，只能正常按照递推关系进行状态了，当然有可能有组合数学的方法，比如指数型母函数这种。

题目

70. 爬楼梯

假设你正在爬楼梯。需要 n 阶你才能到达楼顶。

每次你可以爬 1 或 2 个台阶。你有多少种不同的方法可以爬到楼顶呢？

提示：

1
1 <= n <= 45

示例 1：
输入：n = 2
输出：2
解释：有两种方法可以爬到楼顶。

  1. 1 阶 + 1 阶
  2. 2 阶
    示例 2：
    输入：n = 3
    输出：3
    解释：有三种方法可以爬到楼顶。
  3. 1 阶 + 1 阶 + 1 阶
  4. 1 阶 + 2 阶
  5. 2 阶 + 1 阶

算法：动态规划

1
2
3
4
5
6
7
8
9
10
11
12
状态定义：
dp[i] := i 阶楼梯爬到楼顶的方法数

答案：
dp[n]

初始化：
dp[1] = 1
dp[2] = 2

状态转移：
dp[i] = dp[i - 1] + dp[i - 2]

代码 (C++)

1
2
3
4
5
6
7
8
9
10
11
12
13
class Solution {
public:
    int climbStairs(int n) {
        if(n == 1) return 1;
        if(n == 2) return 2;
        vector<int> dp(n + 1);
        dp[1] = 1;
        dp[2] = 2;
        for(int i = 3; i <= n; ++i)
            dp[i] = dp[i - 1] + dp[i - 2];
        return dp[n];
    }
};

优化：矩阵快速幂

前面的动态规划算法的状态转移方层如下：

由于这个状态转移方程是线性递推关系，我们可以这种递推关系用矩阵表示。具体做法是先写出以下方程组：

然后将该方程组转换为矩阵形式:

其中 $A = \begin{bmatrix} 0 & 1 \\ 1 & 1 \end{bmatrix}$。于是我们所求的结果如下：

代码 (C++)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
using ll = long long;

class Matrix {
public:
    Matrix(int m){
        M = m;
        mat = vector<vector<int> >(M, vector<int>(M, 0));
    }

    void assign(int i, int j, int val)
    {
        mat[i][j] = val;
    }

    int get(int i, int j)
    {
        return mat[i][j];
    }

    void init()
    {
        for(int i = 0; i < M; ++i)
            for(int j = 0; j < M; ++j)
                mat[i][j] = 0;
        for(int i = 0; i < M; ++i)
            mat[i][i] = 1;
    }

    Matrix operator*(const Matrix& matrix) const
    {
        Matrix ans(M);
        for(int i = 0; i < M; ++i)
            for(int j = 0; j < M; ++j)
                for(int k = 0; k < M; ++k)
                    ans.mat[i][j] += ((ll)mat[i][k] * matrix.mat[k][j]);
        return ans;
    }

    Matrix operator^(int n) const
    {
        Matrix ans(M);
        ans.init();
        Matrix A = *this; // 拷贝出来用于自乘
        while(n)
        {
            if(n & 1) ans = ans * A;
            A = A * A;
            n >>= 1;
        }
        return ans;
    }

    void operator=(const Matrix& matrix)
    {
        for(int i = 0; i < M; ++i)
            for(int j = 0; j < M; ++j)
                mat[i][j] = matrix.mat[i][j];
    }

private:
    vector<vector<int>> mat;
    int M;
};

class Solution {
public:
    int climbStairs(int n) {
        if(n == 1) return 1;
        if(n == 2) return 2;
        Matrix mat(2);
        mat.assign(0, 1, 1);
        mat.assign(1, 0, 1);
        mat.assign(1, 1, 1);
        Matrix mat_power = mat ^ (n - 1);
        int ans = mat_power.get(0, 0) + mat_power.get(0, 1) * 2;
        return ans;
    }
};

算法：母函数

对于数列 $dp[1], dp[2], …$，我们可以写出母函数:

其中 $dp[0]$ 是为了后续计算方便加的一项，为了满足 $dp[2] = dp[1] + dp[0]$，$dp[0] = 1$。

$G(x)$ 中的 $x^{n}$ 这一项的系数就是我们要求的 $dp[n]$。

构造 $xG(x)$ 和 $x^{2}G(x)$：

看前面三个式子中的 $x^{2}$ 这一项，系数分别有 $dp[0], dp[1], dp[2]$，然后我们根据递推关系对几个公式进行对应的操作。

由于 $dp[i] = dp[i - 1] + dp[i - 2]$，对应于前面的递推关系的操作如下：

于是：

设 r, s 为 $1 - x - x^{2}$ 的两个根，则：

其中 k1, k2 为待定系数。解得：$k1 = \frac{1}{r - s}$, $k2 = \frac{-1}{r - s}$。于是：

由于：

所以:

于是 $x^{n}$ 的系数为 $\frac{1}{r - s}(\frac{1}{r^{n+1}} - \frac{1}{s^{n+1}})$，这正是我们要求的 $dp[n]$。

解方程 $1 - x - x^{2} = 0$ 得 $r = \frac{-1 + \sqrt{5}}{2}$, $s = \frac{-1 - \sqrt{5}}{2}$。

代码 (C++)

1
2
3
4
5
6
7
8
9
10
class Solution {
public:
    int climbStairs(int n) {
        if(n == 1) return 1;
        if(n == 2) return 2;
        double r = (-1 + sqrt(5)) / 2;
        double s = (-1 - sqrt(5)) / 2;
        return round((pow(1 / r, n + 1) - pow(1 / s, n + 1)) / (r - s));
    }
};

总结: 用母函数求通项公式的步骤

前面的母函数的做法，我们主要进行了以下三步：

(1) 写出 $dp[n]$ 的母函数 G(x)。
(2) 找到用有限项代数式表示的 G(x)。这一步需要根据递推关系来设计。
(3) 用级数的形式表达第(2)步的有限项代数式，此时 $x^{n}$ 的系数就是 $dp[n]$ 的通项公式。

Share