矩阵快速幂

最新推荐文章于 2023-02-11 18:05:48 发布

苏叶秋

最新推荐文章于 2023-02-11 18:05:48 发布

阅读量5k

点赞数 1

分类专栏： ACM日常学习

ACM日常学习专栏收录该内容

50 篇文章 2 订阅

订阅专栏

转：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次乘法运算，代码实现如下：

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

根据上面分析，快速幂的代码实现：

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, 这样直接乘效率也比较低，可以通过快速乘法进一步优化：

#快速乘法

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]

代码实现：

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，使用上面快速幂的方法同样可以实现矩阵快速幂。

结合快速幂和矩阵乘法，实现代码：

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)]=…=[1110]n−1×[F(1)F(0)][F(n)F(n−1)]=[1110]×[F(n−1)F(n−2)]=[1110]2×[F(n−2)F(n−3)]=…=[1110]n−1×[F(1)F(0)]

得到： [F(n)F(n−1)]=[1110]n−1×[F(1)F(0)][F(n)F(n−1)]=[1110]n−1×[F(1)F(0)]

我们知道[F(1)=1F(0)=0][F(1)=1F(0)=0],且[1110]n−1[1110]n−1可以用上面介绍的矩阵快速幂计算出来，两者相乘得到矩阵的第1行第1列就是F(n)F(n)了。

实现代码：

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

结束。

苏叶秋

关注

1
点赞
踩
3

收藏

觉得还不错? 一键收藏
0
评论
矩阵快速幂

转：http://huanyouchen.github.io/2018/05/23/Quick-Matrix-Pow/矩阵的快速幂是用来高效地计算矩阵的高次方的。将朴素的o（n）的时间复杂度，降到log（n）。本文先学习快速幂和矩阵乘法的基础知识，然后将两者结合实现矩阵快速幂方法。然后举一个例子：使用矩阵快速幂求斐波那契数列。快速幂一般计算底数x的n次幂xnxn的方法：xn=x×...
复制链接

扫一扫