Python稀疏三对角矩阵求逆加速及元素运算形状异常问题

一、为什么表达式从矩阵变成向量？ #

主要有两个核心原因：

1. spsolve的返回值类型误解
   scipy.sparse.linalg.spsolve(A, b)的作用是求解线性方程组A @ x = b：当b是向量时，返回解向量x（形状(n,)）；当b是矩阵时，会按列求解并返回由各列解组成的矩阵。如果你的代码中错误地用spsolve直接替代inv(A)，且传入的第二个参数不是单位矩阵，就会得到向量结果，后续运算自然也会变成向量维度。
2. 稀疏矩阵的*运算符重载差异
   在scipy.sparse中，稀疏矩阵的*运算符执行的是矩阵乘法，而非NumPy数组的逐元素乘法。如果spsolve返回了向量，再和矩阵用*相乘，会触发矩阵-向量乘法，最终得到向量而非逐元素相乘的矩阵。若要执行逐元素乘法，需使用scipy.sparse.multiply()函数或稀疏矩阵的.multiply()方法。

二、如何加速稀疏三对角矩阵的“求逆”操作？ #

直接对大型稀疏矩阵求逆效率极低，且三对角矩阵有特殊结构，优先通过避免直接求逆来优化，具体方案如下：

• 将求逆转化为解线性方程组
  任何需要inv(A) @ B的场景，都可以转化为求解A @ X = B，用spsolve(A, B)实现。对于三对角矩阵，spsolve会自动利用其稀疏结构加速，比直接求逆快几个数量级。
• 使用三对角/带状矩阵专用求解器
  三对角矩阵属于带状矩阵，可使用scipy.linalg.solve_banded进一步提速：
  1. 把三对角矩阵转化为带状格式：提取主对角线、上对角线、下对角线，组成形状为(3, n)的数组（第一行是上对角线，第二行主对角线，第三行下对角线）。
  2. 调用solve_banded((1,1), ab, b)（(1,1)表示上下各有1条对角线），速度远快于通用稀疏求解器。
• 仅计算逆矩阵的必要元素
  如果不需要完整的逆矩阵，只需要特定行或列的元素，直接求解对应的线性方程组即可（比如求逆矩阵的第k列，就是解A @ x = e_k，e_k是第k位为1的单位向量），避免计算整个逆矩阵。
• 利用矩阵的结构特性
  三对角矩阵的逆矩阵是满矩阵，但存在递推公式可以计算其元素（针对对称或非对称三对角矩阵有不同的递推方法），如果必须获取完整逆矩阵，可以尝试实现对应的递推算法，比通用求逆方法高效。

内容的提问来源于stack exchange，提问作者THATS MY QUANT MY QUANTITATIVE