优化18×18符号矩阵求逆：加速Ax=b求解方案问询

问题背景 #

尝试用SymPy计算18×18符号矩阵的逆，在Google Colab中运行2小时仍无结果，核心需求是求解线性方程组Ax=b，原方案通过求逆矩阵实现。

数学层面优化 #

1. 避免直接求逆，直接解线性方程组
   求逆矩阵的计算复杂度为O(n³)，符号计算的开销远大于数值计算。对于Ax=b的需求，直接调用SymPy的线性方程组求解接口，无需先计算逆矩阵，能利用矩阵的稀疏性和结构特征减少计算量。

2. 分块矩阵求逆
   观察矩阵结构：前14行/列以常数为主，仅后4行包含符号变量x₁~x₁₄。将矩阵拆分为分块形式：

   [ A  B ]
   [ C  D ]
   

   其中A是14×14的常数矩阵，可快速计算其逆；再利用分块矩阵的逆公式：

   A⁻¹ + A⁻¹B(D - CA⁻¹B)⁻¹CA⁻¹   -A⁻¹B(D - CA⁻¹B)⁻¹
   -(D - CA⁻¹B)⁻¹CA⁻¹            (D - CA⁻¹B)⁻¹
   

   仅需对4×4的符号子矩阵(D - CA⁻¹B)求逆，大幅降低符号计算的规模。

代码层面优化 #

方案1：直接解线性方程组（推荐） #

无需计算逆矩阵，直接求解Ax=b：

import sympy as sp

# 定义符号变量
x_1, x_2, x_3, x_4, x_5, x_6, x_7, x_8, x_9, x_10, x_11, x_12, x_13, x_14 = sp.symbols('x_1:15')
# 定义方程组的未知向量（18个变量）
x_vec = sp.symbols('x_1:x_19')
# 定义右侧向量b（示例，替换为实际需求）
b = sp.Matrix([1]*18)

# 定义原矩阵
matrix = sp.Matrix([
    [1, 1, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
    [1, 2, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
    [1, 2, 1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
    [0, 0, 1, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
    [0, 0, 2, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
    [0, 1, 2, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 1],
    [0, 2, 2, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 1],
    [1, 1, 2, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1],
    [1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 1],
    [1, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1],
    [0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 1],
    [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1],
    [0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 1],
    [0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 1],
    [0, 0, 0, x_1, x_2, x_3, x_4, x_5, x_6, x_7, x_8, x_9, x_10, x_11, x_12, x_13, x_14, -1],
    [0, 0, 0, x_1, x_2, x_3, 0, 0, 0, 0, x_8, x_9, x_10, 0, 0, 0, 0, 0],
    [0, 0, 0, x_1, (2*x_2), (2*x_3), 0, 0, x_6, (2*x_7), x_8, x_9, (2*x_10), (2*x_11), x_12, (2*x_13), x_14, 0],
    [0, 0, 0, 0, 0, x_3, x_4, (2*x_5), (2*x_6), 2 * x_7, 2 * x_8, x_9, 2 * x_10, x_11, 0, 0, x_14, 0]
])

# 直接求解Ax = b
solution = sp.linsolve((matrix, b), x_vec)
print(solution)

方案2：分块矩阵求逆 #

如果确实需要逆矩阵，利用分块优化计算：

import sympy as sp

# 定义符号变量
x_1, x_2, x_3, x_4, x_5, x_6, x_7, x_8, x_9, x_10, x_11, x_12, x_13, x_14 = sp.symbols('x_1:15')

# 定义原矩阵
matrix = sp.Matrix([
    [1, 1, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
    [1, 2, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
    [1, 2, 1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
    [0, 0, 1, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
    [0, 0, 2, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
    [0, 1, 2, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 1],
    [0, 2, 2, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 1],
    [1, 1, 2, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1],
    [1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 1],
    [1, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1],
    [0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 1],
    [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1],
    [0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 1],
    [0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 1],
    [0, 0, 0, x_1, x_2, x_3, x_4, x_5, x_6, x_7, x_8, x_9, x_10, x_11, x_12, x_13, x_14, -1],
    [0, 0, 0, x_1, x_2, x_3, 0, 0, 0, 0, x_8, x_9, x_10, 0, 0, 0, 0, 0],
    [0, 0, 0, x_1, (2*x_2), (2*x_3), 0, 0, x_6, (2*x_7), x_8, x_9, (2*x_10), (2*x_11), x_12, (2*x_13), x_14, 0],
    [0, 0, 0, 0, 0, x_3, x_4, (2*x_5), (2*x_6), 2 * x_7, 2 * x_8, x_9, 2 * x_10, x_11, 0, 0, x_14, 0]
])

# 拆分矩阵为分块
A = matrix[:14, :14]
B = matrix[:14, 14:]
C = matrix[14:, :14]
D = matrix[14:, 14:]

# 计算常数块A的逆
A_inv = A.inv()

# 计算中间矩阵
temp = D - C @ A_inv @ B
temp_inv = temp.inv()

# 组合分块逆矩阵
top_left = A_inv + A_inv @ B @ temp_inv @ C @ A_inv
top_right = -A_inv @ B @ temp_inv
bottom_left = -temp_inv @ C @ A_inv
bottom_right = temp_inv

# 拼接得到完整逆矩阵
matrix_inv = sp.BlockMatrix([[top_left, top_right], [bottom_left, bottom_right]]).as_explicit()

print(matrix_inv)

额外优化建议 #

• 若变量有明确范围（如正数），定义时添加约束：sp.var('x_1:15', positive=True)，帮助SymPy简化表达式。
• 避免不必要的符号展开，仅在需要时调用sp.expand()，减少中间计算量。

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