矩阵快速幂模板以及其应用矩阵加速递推

矩阵快速幂模板（Python实现）

• 不调用模块实现矩阵快速幂

##矩阵快速幂算法
def matrix_mul(A, B): # 等价于求BxA
#矩阵乘法函数，返回两个矩阵相乘的值
    return [[sum(a * b  for a, b in zip(col, row))  for col in zip(*B)] for row in A]

def matrix_pow(A, n):
    size_ = len(A)
    if n == 0:  #返回单位矩阵
        res = [[0 for _ in range(size_)] for _ in range(size_)]
        for i in range(size_):
            res[i][i] = 1
        return res
    elif n == 1:    #返回自己
        return A
    else:
        y = matrix_pow(A, n // 2)
        if n & 1:   #要乘
            return matrix_mul(matrix_mul(y, y), A)
        return matrix_mul(y, y) #不乘

if __name__ == '__main__':
    matrix = [[1, 2, 3],
              [1, 0, 0],
              [0, 1 ,0]]
    print(matrix_pow(matrix, n=10))

注意：应用这种方法定义的矩阵是这样的规律

matrix = [[a, b, c],[d, e, f],[g, h, i]]

这里的matrix在以上代码中是作为$\left(\begin{array}{c}adg\\ beh\\ cfi\end{array}\right)$来进行使用的；

如果需要在运算过程中加入取模操作的话，可以如下操作

module = xxx
def matrix_mul(A, B):
#矩阵乘法函数，返回两个矩阵相乘的值，建议背诵
    return [[sum(a * b % module for a, b in zip(col, row)) % module for col in zip(*B)] for row in A]

def matrix_pow(A, n):
    size_ = len(A)
    if n == 0:#返回单位矩阵
        res = [[0 for _ in range(size_)] for _ in range(size_)]
        for i in range(size_):
            res[i][i] = 1
        return res
    elif n == 1:#返回自己
        return A
    else:
        y = matrix_pow(A, n // 2)
        if n & 1:#要乘
            return matrix_mul(matrix_mul(y, y), A)
        return matrix_mul(y, y)#不乘

代码来源及原理讲解：(64条消息) 矩阵快速幂算法及相关应用（含python源码）_国际知名观众的博客-CSDN博客_矩阵快速幂算法

---

• 调用numpy模块实现矩阵快速幂

import numpy as np
matrix = np.matrix([[3, 1, 0],
                    [4, 0, 1],
                    [5, 0, 0]], dtype=object) 
# 注意：定义矩阵时一定要加上dtype=object,否则当矩阵内数据过大时会自动取模
print(matrix * matrix) # 使用的是矩阵乘法
print(np.matmul(matrix, matrix)) # 矩阵乘法
print(np.linalg.matrix_power(matrix, n=10)) # 矩阵快速幂

注意：应用这种方法定义的矩阵是这样的规律

matrix = np.matrix([[a, b, c],[d, e, f],[g, h, i]])

这里的matrix在以上代码中是作为$\left(\begin{array}{c}abc\\ def\\ ghi\end{array}\right)$来进行使用的

矩阵加速递推

所谓矩阵加速递推实际上是使用矩阵快速幂这一快速计算矩阵幂次的工具来简化运算，并且简化的是求有彼此有相关关系的数列的运算，如斐波那契额数列的第n项求解（且n往往极大，按照数列表面的数学关系直接递推计算几乎不可能得到）等等；计算这样的递推次数很多的式子，我们需要一个更简化更快捷的计算方式，也就是矩阵快速幂

首先我们需要想办法把实际问题转换为矩阵的形式来表示，这也是矩阵加速递推最关键的步骤

斐波那契数列

其递推关系是
$$F_1=F_2=1;F_i=F_{i-1}+F_{i-2}(i\ge 3)$$
根据其递推关系我们考虑将矩阵初始值设置为一个行向量（也可以理解成1x2的矩阵，其实列向量还是行向量都可以，只要转移矩阵也是对应的即可）
$$(F_i{F}_{i-1})$$
构造一个2x2的转移矩阵（举例：当求该斐波那契数列的第3项时，矩阵初始值乘以该转移矩阵得到的新矩阵的第一项数据即是斐波那契数列的第3项）

我们想要通过乘以转移矩阵，得到数列的下一项，也就是说
$$(F_i{F}_{i-1})\times \left(\begin{array}{c}ab\\ cd\end{array}\right)=(F_{i+1}{F}_i)$$
其中$\left(\begin{array}{c}ab\\ cd\end{array}\right)$ 即是我们需要构造的转移矩阵，另外，矩阵乘法不满足交换律，需要注意乘数左右位置；把矩阵乘法展开，得到
$$\begin{gathered} a\cdot F_i+c\cdot F_{i-1}=F_{i+1} \\ b\cdot F_i+d\cdot F_{i-1}=F_i \end{gathered}$$
那么再根据数列本来的关系$F_{i+1}=F_i+F_{i-1}$，对照到第一个等式，我们可以得到a = 1, c = 1；而第二个等式显而易见，b = 1, d = 0

这样就求得了转移矩阵$\left(\begin{array}{c}11\\ 10\end{array}\right)$；那么根据转移矩阵表现出来的与原数列的关系，我们可以这样来看
$$(F_n{F}_{n-1})=(F_2{F}_1)\times {\left(\begin{array}{c}11\\ 10\end{array}\right)}^{n-2}$$
那么到此为止，我们就知道了如何使用矩阵快速幂来求解递推性质的n项数列了；接下来用Python验证一遍

def matrix_mul(A, B):
    return [[sum(a * b  for a, b in zip(col, row))  for col in zip(*B)] for row in A]

def matrix_pow(A, n):
    size_ = len(A)
    if n == 0: 
        res = [[0 for _ in range(size_)] for _ in range(size_)]
        for i in range(size_):
            res[i][i] = 1
        return res
    elif n == 1: 
        return A
    else:
        y = matrix_pow(A, n // 2)
        if n & 1:  
            return matrix_mul(matrix_mul(y, y), A)
        return matrix_mul(y, y) 

def fast_power(matrix, expo, ini):
    return  matrix_mul(matrix_pow(matrix, expo), ini) 

matrix = [[1, 1],  # 转移矩阵；注意:matrix[i]表示的是第i列
          [1, 0]]
ini = [[1], [1]] # 矩阵初始值(F1, F2)；外列表表示行，里列表表示列
n = 30 # 数列的第n项
expo = n - 2 # 转移矩阵幂次 = n - 2
res = fast_power(matrix, expo, ini)

def Fibonacci(n):
    if n == 1 or n == 2:
        return 1
    else:
        return Fibonacci(n - 1) + Fibonacci(n - 2)
    
assert Fibonacci(n) == res[0][0] # 注意是取矩阵的第一项
print(res)
print(Fibonacci(n))
'''
[[832040], [514229]]
832040
'''

那么回到最初的构造矩阵初始值，我们为什么要构造的是1x2的矩阵而不是1x3或者2x2的矩阵呢？

仔细查看数列的递推关系$F_i=F_{i-1}+F_{i-2}$，如果我们构造的是1x3矩阵，也就是$(F_i{F}_{i-1}{F}_{i-2})$；假如要求下一项数据的时候，也就是求$(F_{i+1}{F}_i{F}_{i-1})$这个矩阵

那么需要构造的转移矩阵显然需要是3x3的大小（这样才能在矩阵乘法中覆盖每一项数据）
$$(F_i{F}_{i-1}{F}_{i-2})\times \left(\begin{array}{c}abc\\ def\\ ghi\end{array}\right)=(F_{i+1}{F}_i{F}_{i-1})$$
展开后单独拿$F_{i+1}$来说
$$a\cdot F_i+d\cdot F_{i-1}+g\cdot F_{i-2}=F_{i+1}$$
根据数列递推关系$F_{i+1}=F_i+F_{i-1}$；很显然a = 1, d = 1, g = 0；相比于我们构造1x2的矩阵初始值的时候是不是仅仅只是多了g = 0，也就是说构造1x3的矩阵初始值浪费了计算空间，所以构造其他形式的初始矩阵可以但没必要；我们只需要构造最小限度能解决问题的矩阵就好，无论是矩阵初始值还是转移矩阵，所以还是需要具体问题具体分析地来构造矩阵

其他例子

• 来源于2022 CTF网鼎杯-青龙组 david_homework

$Problem$

cof_t = [[353, -1162, 32767], [206, -8021, 42110], [262, -7088, 31882]] # 原题列表很长，这里只提取一部分做解释
def cal(i,cof):
    if i < 3:
        return i+1
    else:
        return  cof[2]*cal(i-3,cof)+cof[1]*cal(i-2,cof)+cof[0]*cal(i-1,cof)

s = 0
for i in range(len(cof_t)):
    s += cal(200000,cof_t[i])

求s的部分数位

$Solution$

套用矩阵加速递推，将问题转换为矩阵幂次的形式

关键代码cof[2]*cal(i-3,cof)+cof[1]*cal(i-2,cof)+cof[0]*cal(i-1,cof)

实际上可以转换为数列的形式
$$\begin{gathered} F_i=c_2\cdot F_{i-3}+c_1\cdot F_{i-2}+c_0\cdot F_{i-1} \\ 其中F_0=1,F_1=2,F_2=3 \end{gathered}$$
其中$F_i$表示cal(i, cof)，$c_i$代表cof[i]

构造初始矩阵$(F_{i-1}{F}_{i-2}{F}_{i-3})$，习惯上将数列次序以从右到左增加的方式排列初始矩阵；那么为什么我们这里构造的又是1x3的初始矩阵呢？

因为如果构成1x2的初始矩阵的话，对应数列的递推关系是不完整的；就是说$(F_{i-2}{F}_{i-3})$要变成$(F_{i-1}{F}_{i-2})$，需要参考的等式是$F_i=c_2\cdot F_{i-3}+c_1\cdot F_{i-2}+c_0\cdot F_{i-1}$，显然$F_i$没有对应的变量表示，则无法得到结果

那么初始矩阵对应的转移矩阵也就应当是$\left(\begin{array}{c}abc\\ def\\ ghi\end{array}\right)$

根据

$(F_{i-1}{F}_{i-2}{F}_{i-3})\times \left(\begin{array}{c}abc\\ def\\ ghi\end{array}\right)=(F_i{F}_{i-1}{F}_{i-2})$

展开后
$$\begin{gathered} a\cdot F_{i-1}+d\cdot F_{i-2}+g\cdot F_{i-3}=F_i \\ b\cdot F_{i-1}+e\cdot F_{i-2}+h\cdot F_{i-3}=F_{i-1} \\ c\cdot F_{i-1}+f\cdot F_{i-2}+i\cdot F_{i-3}=F_{i-2} \end{gathered}$$
等式1对应$F_i=c_2\cdot F_{i-3}+c_1\cdot F_{i-2}+c_0\cdot F_{i-1}$，显然a = c0, d = c1, g = c2；等式2观察等式易得b = 1, e = 0, h = 0；等式3观察等式易得c = 0, f = 1, i = 0

综上，转移矩阵为$\left(\begin{array}{c}{c}_{0}10\\ c_{1}01\\ c_{2}00\end{array}\right)$，用代码验证一下是否是这样的矩阵递推关系
$$\begin{gathered} (F_{i-1}{F}_{i-2}{F}_{i-3})\times \left(\begin{array}{c}{c}_{0}10\\ c_{1}01\\ c_{2}00\end{array}\right)=(F_i{F}_{i-1}{F}_{i-2}) \\ (F_n{F}_{n-1}{F}_{n-2})=(F_2{F}_1{F}_0)\times {\left(\begin{array}{c}{c}_{0}10\\ c_{1}01\\ c_{2}00\end{array}\right)}^{n-2} \end{gathered}$$

from tqdm import tqdm
##矩阵快速幂算法
def matrix_mul(A, B):
#矩阵乘法函数，返回两个矩阵相乘的值
    return [[sum(a * b  for a, b in zip(col, row))  for col in zip(*B)] for row in A]

def matrix_pow(A, n):
    size_ = len(A)
    if n == 0:  #返回单位矩阵1
        res = [[0 for _ in range(size_)] for _ in range(size_)]
        for i in range(size_):
            res[i][i] = 1
        return res
    elif n == 1:    #返回自己
        return A
    else:
        y = matrix_pow(A, n // 2)
        if n & 1:   #要乘
            return matrix_mul(matrix_mul(y, y), A)
        return matrix_mul(y, y) #不乘

def fast_power(matrix, expo, ini):
    return  matrix_mul(matrix_pow(matrix, expo), ini) # 幂数减去的数以实际情况做出改变


cof_t = [[353, -1162, 32767], [206, -8021, 42110], [262, -7088, 31882]]

res = 0
n = 20
for ls in tqdm(cof_t):
    c2, c1, c0 = ls[2], ls[1], ls[0]
    matrix = [[c0, c1, c2], # 注意:matrix[i]表示的是第i列
              [1, 0, 0],
              [0, 1 ,0]]
    ini = [[3],[2],[1]] #初始值
    expo = n - 2
    a = fast_power(matrix, expo, ini) #幂数减去第一个，乘上初始的
    # print(a[0][0])
    res = res + a[0][0]

def cal(i,cof):
    if i <3:
        return i+1
    else:
        return  cof[2]*cal(i-3,cof)+cof[1]*cal(i-2,cof)+cof[0]*cal(i-1,cof)
s = 0
for i in range(len(cof_t)):
    s+= cal(n, cof_t[i])
assert s == res
print(s)
print(res)
'''
565675190172910086230288153171289084854733110754
565675190172910086230288153171289084854733110754
'''

$Reference$

【矩阵性质】矩阵加速递推 - RioTian - 博客园 (cnblogs.com)

(64条消息) 矩阵乘法加速递推_W⁡angduoyu的博客-CSDN博客_矩阵乘法加速

辰星凌的博客QAQ (cnblogs.com)