矩阵的快速幂是用来高效地计算矩阵的高次方的。将朴素的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^2 + 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

结束。