[转载]快速幂与矩阵快速幂
来源:http://huanyouchen.github.io/2018/05/23/Quick-Matrix-Pow/
矩阵快速幂
矩阵的快速幂是用来高效地计算矩阵的高次方的。将朴素的o(n)的时间复杂度,降到log(n)。
本文先学习快速幂和矩阵乘法的基础知识,然后将两者结合实现矩阵快速幂方法。然后举一个例子:使用矩阵快速幂求斐波那契数列。
快速幂
一般计算底数x的n次幂xnxn 的方法: xn=x×x×x…x×xxn=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次方x11x11,而11可以用二进制表示:11=1011=1×23+0×22+1×21+1×2011=1011=1×23+0×22+1×21+1×20
那么可以把x11x11转换为:
x11=x23×x21×x20=x8×x2×x1x11=x23×x21×x20=x8×x2×x1
原先需要做11次乘法运算,转换后只需要3次乘法运算。
通过上面的具体例子来推导一般情况计算xnxn 的方法,先把n转换为2进制,从低位到高位根据二进制中的0或者1来进行乘法运算,比如上面的x11x11例子,从低位到高位的运算过程:
- 11 = 1011
- 1:res = x20x20
- 1: res = x20×x21x20×x21
- 0: 跳过,不运算
- 1:res = x20×x21×x23x20×x21×x23
得到结果res = x1×x2×x8=x11x1×x2×x8=x11
判断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很大时,比如98765432119876543211, 这样直接乘效率也比较低,可以通过快速乘法进一步优化:
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×nm×n的矩阵就是m×nm×n个数排成m行n列的一个数阵。
矩阵乘法:设A为m×pm×p的矩阵,B为p×np×n的矩阵,那么称m×nm×n的矩阵C为矩阵A与B的乘积,记作C=AB。
矩阵A与矩阵B相乘需要满足条件:A的列数等于B的行数。
矩阵乘法举例:
令A=[1234]A=[1234] , B=[5678]B=[5678], 则:
C=AB=[1234]×[5678]=[1×5+2×71×6+2×83×5+4×73×6+4×8]=[19224350]C=AB=[1234]×[5678]=[1×5+2×71×6+2×83×5+4×73×6+4×8]=[19224350]
代码实现:
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]]
|
且矩阵乘法满足结合律:
A11=A8×A2×A1=[1234]8×[1234]2×[1234]1A11=A8×A2×A1=[1234]8×[1234]2×[1234]1
矩阵快速幂
现在知道了两个矩阵乘法运算,那么如果求矩阵A的n次方AnAn,可以用An=A×A×A×…×AAn=A×A×A×…×A,不过学习了计算底数xx的n次幂xnxn 的快速幂方法,把底数xx换成矩阵A,计算AnAn,使用上面快速幂的方法同样可以实现矩阵快速幂。
结合快速幂和矩阵乘法,实现代码:
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)=0F(0)=0
- F(1)=1F(1)=1
- F(n)=F(n−1)+F(n−2)(n>=2)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)F(n)=F(n−1)+F(n−2)用矩阵表示:
[F(n)=F(n−1)+F(n−2)F(n−1)]=[1110]×[F(n−1)F(n−2)][F(n)=F(n−1)+F(n−2)F(n−1)]=[1110]×[F(n−1)F(n−2)]
而F(n−1)、F(n−2)、F(n−3)F(n−1)、F(n−2)、F(n−3)也可以用同样的矩阵表示方式,于是有:
[F(n)F(n−1)]=[1110]×[F(n−1)F(n−2)]=[1110]2×[F(n−2)F(n−3)]=…=