Python中加速依赖变量的矩阵逆运算的技术问询

我面临一个迭代问题，每次迭代需执行一系列矩阵运算，希望对其进行加速。计算中大部分为常量，仅存在一个随迭代变化的变量。请问是否可以构建此类运算的“骨架”，以加速首次迭代后的计算（类似CVXPY优化问题中的参数化预热启动）？相关代码位于类的专用函数中，具体如下：

# some constants definitions (in a class)
m = 5501
Y = np.random.rand(10, m)
U = np.random.rand(75, m)
sig_on = 0.15
sig_off = 0.25
a = 15
L = 25

def obtain_A_B(x):
   # for readability I omit internal variable declarations as a = self.a
   varying_part = (a * sig_on / np.linalg.norm(x)) ** 2 + L * sig_off
   F = varying_part * np.identity(m) + Y.T @ Y
   Fi = np.linalg.inv(F)
   UFiUTi = np.linalg.inv(U @ Fi @ U.T)
   A = (Fi - Fi @ U.T @ UFiUTi  @ U @ Fi) @ Y.T
   B = Fi @ U.T @ UFiUTi 
   return A, B

x = np.random.rand(10)
result = obtain_A_B(x)

我希望加速A和B的计算，其中varying_part是随x变化的标量。当前最耗时的是计算逆矩阵Fi，仅加速这一步就能带来很大帮助。我已将逆运算单独存储到变量中，将执行时间从约41秒降至约10秒，但希望能进一步提速，因此想到了“运算骨架”的思路。

核心思路是利用低秩矩阵的逆运算公式（Sherman-Morrison-Woodbury）搭建固定运算骨架，预计算所有常量部分，每次迭代仅处理随varying_part变化的小矩阵运算，完全避免大矩阵（5501×5501）的逆运算。

1. 预计算常量骨架 #

在类的初始化阶段，提前计算所有不随迭代变化的矩阵项，构成固定运算骨架：

class MatrixCalculator:
    def __init__(self):
        self.m = 5501
        self.Y = np.random.rand(10, self.m)
        self.U = np.random.rand(75, self.m)
        self.sig_on = 0.15
        self.sig_off = 0.25
        self.a = 15
        self.L = 25
        
        # 预计算常量矩阵，后续迭代不再重复计算
        self.YYt = self.Y @ self.Y.T  # 10×10 低秩矩阵
        self.UYt = self.U @ self.Y.T  # 75×10
        self.UUt = self.U @ self.U.T  # 75×75
        self.Yt = self.Y.T  # m×10
        self.YUt = self.Y @ self.U.T  # 10×75 减少重复乘法运算

2. 重构迭代计算逻辑 #

每次迭代仅处理随x变化的标量λ = varying_part，所有逆运算均针对小矩阵（10×10或75×75），完全规避大矩阵操作：

def obtain_A_B(self, x):
    # 计算随x变化的标量λ及相关倒数
    x_norm = np.linalg.norm(x)
    λ = (self.a * self.sig_on / x_norm) ** 2 + self.L * self.sig_off
    λ_inv = 1 / λ
    λ_inv_sq = λ_inv ** 2
    
    # 计算10×10小矩阵的逆，用solve替代inv稳定性更好
    C = np.eye(10) + self.YYt * λ_inv
    C_inv = np.linalg.solve(C, np.eye(10))
    
    # 计算UFiUT，避免显式生成大矩阵Fi
    UFiUT = λ_inv * self.UUt - λ_inv_sq * self.UYt @ C_inv @ self.UYt.T
    # 计算75×75矩阵的逆
    UFiUT_inv = np.linalg.solve(UFiUT, np.eye(75))
    
    # 计算A，拆分项避免大矩阵运算
    proj_Yt = self.Yt - self.U.T @ UFiUT_inv @ self.U @ self.Yt
    term1 = λ_inv * proj_Yt
    term2 = λ_inv_sq * self.Yt @ C_inv @ (self.YYt - self.YUt @ UFiUT_inv @ self.UYt)
    A = term1 - term2
    
    # 计算B，同样避免大矩阵Fi
    term_B1 = λ_inv * self.U.T @ UFiUT_inv
    term_B2 = λ_inv_sq * self.Yt @ C_inv @ self.YUt @ UFiUT_inv
    B = term_B1 - term_B2
    
    return A, B

3. 优化效果说明 #

• 完全消除了5501×5501大矩阵的逆运算，替换为10×10和75×75小矩阵的逆运算，速度提升至少一个数量级；
• 预计算的常量矩阵构成固定运算骨架，每次迭代仅执行小矩阵乘法和逆运算，内存占用大幅降低（无需存储5501×5501的矩阵Fi）；
• 使用np.linalg.solve替代np.linalg.inv，进一步提升运算稳定性和速度。

内容的提问来源于stack exchange，提问作者Andrea