矩阵的快速幂是用来高效地计算矩阵的高次方的。将朴素的o(n)的时间复杂度,降到log(n)。
本文先学习快速幂和矩阵乘法的基础知识,然后将两者结合实现矩阵快速幂方法。然后举一个例子:使用矩阵快速幂求斐波那契数列。


快速幂

一般计算底数x的n次幂$x^n$ 的方法: $x^n = x × x × x … x × x$ ,需要做n次乘法运算,代码实现如下:

1
2
3
4
5
6
7
def powx(x, n):
res = 1
for i in range(n):
res = res * x
return res
print(powx(x,n))

时间复杂度为O(N), 当n非常大时候运算效率很低。怎么才能提高运算效率来快速计算底数x的n次幂呢?可以通过快速幂方法。

先用一个具体的例子解释其原理:比如,计算x的11次方$x^{11}$,而11可以用二进制表示: $$11 = 1011 = 1 × 2^3 + 0 × 2^0 + 1 × 2^1 + 1 ×2^0$$

那么可以把$x^{11}$转换为:
$$x^{11} = x ^ {2^3} × x ^ {2^1} × x^{2^0} = x^8 × x^2 × x^1 $$

原先需要做11次乘法运算,转换后只需要3次乘法运算。

通过上面的具体例子来推导一般情况计算$x^n$ 的方法,先把n转换为2进制,从低位到高位根据二进制中的0或者1来进行乘法运算,比如上面的$x^{11}$例子,从低位到高位的运算过程:

  • 11 = 1011
  • 1:res = $x^{2^0}$
  • 1: res = $x^{2^0} × x^{2^1}$
  • 0: 跳过,不运算
  • 1:res = $x^{2^0} × x^{2^1} × x^{2^3}$

得到结果res = $x^{1} × x^{2} × x^{8} = x ^ {11}$

判断n的二进制低位是0或者1的方法, 也就是判断n是偶数或者奇数的方法,可以通过位运算and来实现:

  • n and 1 返回1, 则n是奇数,即n的二进制低位是1
  • n and 1 返回0, 则n是偶数,即n的二进制低位是0

从低位到高位的运算过程也可以通过位运算>>实现, n = n >> 1, 把n的二进制位右移过程也就是高位到低位的过程。比如11 = 1011的右移过程:

  • n = n >> 1 = 1011 >> 1 = 101
  • n = n >> 1 = 101 >> 1 = 10
  • n = n >> 1 = 10 >> 1 = 1
  • n = n >> 1 = 1 >> 1 = 0

根据上面分析,快速幂的代码实现:

1
2
3
4
5
6
7
8
9
10
def powx(x, n):
res = 1
while n:
if n&1:
res = res * x
x = x * x
n = n >> 1
return res
print(powx(x,n))

上面快速幂代码中第5,6行涉及到两步乘法运算,当x很大时,比如$98765432^{11}$, 这样直接乘效率也比较低,可以通过快速乘法进一步优化:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#快速乘法
def fast_multi(a,b):
multi_res = 0
while b:
if b&1:
multi_res = multi_res + a
a = a+a
b >>= 1
return multi_res
#快速幂
def powx(x, n):
res = 1
while n:
if n&1:
res = fast_multi(res, x)
x = fast_multi(x, x)
n = n >> 1
return res

矩阵乘法

矩阵:一个$m×n$的矩阵就是$m×n$个数排成m行n列的一个数阵。

矩阵乘法:设A为$m×p$的矩阵,B为$p×n$的矩阵,那么称$m×n$的矩阵C为矩阵A与B的乘积,记作C=AB。
矩阵A与矩阵B相乘需要满足条件:A的列数等于B的行数。

矩阵乘法举例:

令$A = \begin{bmatrix}
1 & 2 \\
3 & 4 \\
\end{bmatrix}$ , $B = \begin{bmatrix}
5 & 6 \\
7 & 8 \\
\end{bmatrix}$, 则:

$$C = AB = \begin{bmatrix}
1 & 2 \\
3 & 4 \\
\end{bmatrix} × \begin{bmatrix}
5 & 6 \\
7 & 8 \\
\end{bmatrix} = \begin{bmatrix}
1×5+2×7 & 1×6+2×8 \\
3×5+4×7 & 3×6+4×8 \\
\end{bmatrix} = \begin{bmatrix}
19 & 22 \\
43 & 50 \\
\end{bmatrix}$$

代码实现:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
def MatrixMultiply(matrix_a, matrix_b):
n_row = len(matrix_a)
n_col = len(matrix_b[0])
# C的行数等于A的行数,C的列数等于B的列数
matrix_c = [[0 for j in range(n_col)] for i in range(n_row)]
for i in range(0, n_row):
for j in range(0, n_col):
for k in range(0, n_row):
matrix_c[i][j] = matrix_c[i][j] + matrix_a[i][k] * matrix_b[k][j]
return matrix_c
a = [[1,2],[3,4]]
b = [[5,6],[7,8]]
c = MatrixMultiply(a, b)
print(c)
# [[19, 22], [43, 50]]

且矩阵乘法满足结合律:

$$A^{11} = A^8 × A^2 × A^1 = \begin{bmatrix}
1 & 2 \\
3 & 4 \\
\end{bmatrix} ^ 8 × \begin{bmatrix}
1 & 2 \\
3 & 4 \\
\end{bmatrix} ^ 2 × \begin{bmatrix}
1 & 2 \\
3 & 4 \\
\end{bmatrix} ^ 1$$


矩阵快速幂

现在知道了两个矩阵乘法运算,那么如果求矩阵A的n次方$A^n$,可以用$A^n = A×A×A×…×A$,不过学习了计算底数$x$的n次幂$x^n$ 的快速幂方法,把底数$x$换成矩阵A,计算$A^n$,使用上面快速幂的方法同样可以实现矩阵快速幂。

结合快速幂和矩阵乘法,实现代码:

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
def MatrixMultiply(matrix_a, matrix_b):
n_row = len(matrix_a)
n_col = len(matrix_b[0])
# C的行数等于A的行数,C的列数等于B的列数
matrix_c = [[0 for j in range(n_col)] for i in range(n_row)]
for i in range(0, n_row):
for j in range(0, n_col):
for k in range(0, n_row):
matrix_c[i][j] = matrix_c[i][j] + matrix_a[i][k] * matrix_b[k][j]
return matrix_c
def getUnitMatrix(l):
# 构造单位矩阵,l为矩阵a的行数
unit_matrix = [[0 for j in range(l)] for i in range(l)]
for k in range(l):
unit_matrix[k][k] = 1
return unit_matrix
def QuickMatrixPow(a, n):
# res_matrix初始化为单位矩阵
res_matrix = getUnitMatrix(len(a))
while n:
if n&1:
res_matrix = MatrixMultiply(res_matrix, a)
a = MatrixMultiply(a,a)
n = n >> 1
return res_matrix
a = [[1,2],[3,4]]
n = 11
print(QuickMatrixPow(a,n))
# [[25699957, 37455814], [56183721, 81883678]]

扩展:斐波那契数列矩阵算法

关于斐波那契数列的定义:

斐波那契数列(意大利语:Successione di Fibonacci),又译为费波拿契数列、费波那西数列、斐波那契数列、黄金分割数列。

在数学上,费波那契数列是以递归的方法来定义:

  • $F(0) = 0$
  • $F(1) = 1$
  • $F(n) = F(n-1) + F(n-2) (n>=2)$

用文字来说,就是费波那契数列由0和1开始,之后的费波那契系数就是由之前的两数相加而得出。首几个费波那契系数是:

0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233……(OEIS中的数列A000045)

特别指出:0不是第一项,而是第零项。

将$ F(n)=F(n-1)+F(n-2)$用矩阵表示:
$$\begin{bmatrix}
F(n)=F(n-1)+F(n-2) \\
F(n-1) \\
\end{bmatrix} = \begin{bmatrix}
1 & 1 \\
1 & 0 \\
\end{bmatrix} × \begin{bmatrix}
F(n-1) \\
F(n-2) \\
\end{bmatrix}$$

而$F(n-1)、F(n-2)、F(n-3)$也可以用同样的矩阵表示方式,于是有:
$$\begin{bmatrix}
F(n) \\
F(n-1) \\
\end{bmatrix} = \begin{bmatrix}
1 & 1 \\
1 & 0 \\
\end{bmatrix} × \begin{bmatrix}
F(n-1) \\
F(n-2) \\
\end{bmatrix} = \begin{bmatrix}
1 & 1 \\
1 & 0 \\
\end{bmatrix}^2 × \begin{bmatrix}
F(n-2) \\
F(n-3) \\
\end{bmatrix} = … = \begin{bmatrix}
1 & 1 \\
1 & 0 \\
\end{bmatrix}^{n-1} × \begin{bmatrix}
F(1) \\
F(0) \\
\end{bmatrix}$$

得到: $\begin{bmatrix}
F(n) \\
F(n-1) \\
\end{bmatrix} = \begin{bmatrix}
1 & 1 \\
1 & 0 \\
\end{bmatrix}^{n-1} × \begin{bmatrix}
F(1) \\
F(0) \\
\end{bmatrix}$

我们知道$\begin{bmatrix}
F(1)=1 \\
F(0)=0 \\
\end{bmatrix}$,且$\begin{bmatrix}
1 & 1 \\
1 & 0 \\
\end{bmatrix}^{n-1}$可以用上面介绍的矩阵快速幂计算出来,两者相乘得到矩阵的第1行第1列就是$F(n)$了。

实现代码:

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
def MatrixMultiply(matrix_a, matrix_b):
n_row = len(matrix_a)
n_col = len(matrix_b[0])
# C的行数等于A的行数,C的列数等于B的列数
matrix_c = [[0 for j in range(n_col)] for i in range(n_row)]
for i in range(0, n_row):
for j in range(0, n_col):
for k in range(0, n_row):
matrix_c[i][j] = matrix_c[i][j] + matrix_a[i][k] * matrix_b[k][j]
return matrix_c
def get_unit_matrix(l):
unit_matrix = [[0 for j in range(l)] for i in range(l)]
for k in range(l):
unit_matrix[k][k] = 1
return unit_matrix
def QuickMatrixPow(a, n):
res_matrix = get_unit_matrix(len(a))
while n:
if n&1:
res_matrix = MatrixMultiply(res_matrix, a)
a = MatrixMultiply(a,a)
n = n >> 1
return res_matrix
def get_Fib_n(n):
if n == 0:
return 0
elif n == 1:
return 1
else:
a = [[1,1],[1,0]]
base = [[1],[0]]
Fib_n = MatrixMultiply(QuickMatrixPow(a,n-1),base)
return Fib_n[0][0]
if __name__ == '__main__':
# 求斐波那契数列第F(n)个数
n = 8
print(get_Fib_n(n))
# 21

结束。