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

  |  

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

【对算法,数学,计算机感兴趣的同学,欢迎关注我哈,阅读更多原创文章】
我的网站:潮汐朝夕的生活实验室
我的公众号:算法题刷刷
我的知乎:潮汐朝夕
我的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}logN)$ ,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