如何在Python/Numpy中快速求解双上三角线性方程组A*X=B？

我需要求解多右端项线性方程组$A*X=B$，其中A和B均为实方阵且是上三角矩阵，规模约为200×200。请问在Python/Numpy中是否存在快速求解方法？

我考虑过以下几种方案：

• 循环求解各列（已补充7×7示例代码）：

import numpy as np
import scipy.linalg as sp

A=np.array(
[[ 1.          0.44615865  0.39541532  0.24977742  0.0881614   0.26116991   0.4138066 ]
 [ 0.          0.89495389  0.24253783  0.4514874   0.12356345  0.22552021   0.48408527]
 [ 0.          0.          0.88590187  0.03860599  0.19887529  0.03114347  -0.02639242]
 [ 0.          0.          0.          0.85573357 -0.05867366  0.85120741   0.25861816]
 [ 0.          0.          0.          0.          0.96641899  0.14020408   0.26514478]
 [ 0.          0.          0.          0.          0.          0.36844234   0.50505032]
 [ 0.          0.          0.          0.          0.          0.           0.44885192]])
B=np.triu(np.array(
  [[ 949.43526038  550.35234482  232.34981032 -176.85444188 -143.39220636  198.43783458   60.7140828 ]]
  ).T @ np.ones((1,7)) )

n=A.shape[0]
X=np.zeros((n,n))
for i in range(n):
    X[:i+1,i]=sp.solve_triangular(A[:i+1,:i+1],B[:i+1,i])

但该方法未利用快速矩阵运算；

• 同时求解所有右端项，即X=solve_triangular(A,B)，但未利用B的三角结构；
• 求A的逆矩阵再与B相乘，即X=inv(A)@B，但通常不推荐矩阵求逆操作。

最优解决方案 #

针对A和B均为上三角矩阵的特性，我们可以通过以下方式实现高效求解：

1. 优先使用scipy.linalg.solve_triangular
   虽然表面上sp.solve_triangular(A,B)没特意利用B的三角结构，但该函数内部会自动忽略输入矩阵的零元素区域。加上A是上三角矩阵，我们只需指定lower=False，函数会以$O(n^{2)$的时间复杂度完成计算，这已经是理论最优的时间复杂度（远低于普通矩阵求解的$O(n}3)$）。

   另外，由于A可逆且是上三角，其逆矩阵也是上三角，上三角矩阵相乘结果仍为上三角，因此X必然是上三角矩阵。如果想进一步手动利用这一特性，可以只计算X的上三角部分：

   import numpy as np
   import scipy.linalg as sp
   
   X = np.zeros_like(B)
   n = A.shape[0]
   for j in range(n):
       # 仅计算第j列的上三角部分
       X[:j+1, j] = sp.solve_triangular(A[:j+1, :j+1], B[:j+1, j], lower=False)
   

   不过对于200×200的规模，原生的solve_triangular已经足够高效，手动优化的性能提升非常有限。

2. 预计算逆矩阵（针对多次求解场景）
   如果需要基于同一A求解多个不同的B，可以预先计算A的上三角逆矩阵，再与B做上三角乘法：

   inv_A = sp.inv(A)
   X = np.zeros_like(B)
   n = A.shape[0]
   for i in range(n):
       for j in range(i, n):
           # 仅计算上三角元素的乘积
           X[i, j] = np.dot(inv_A[i, i:j+1], B[i:j+1, j])
   

   但单次求解时，这种方法不如直接用solve_triangular高效，因为求逆加上三角乘法的总耗时略高于直接求解。

性能总结 #

对于200×200的矩阵：

• sp.solve_triangular(A, B, lower=False)是最简洁高效的选择，执行时间在微秒级；
• 手动循环计算上三角部分的性能与原生函数接近，但原生函数经过底层优化，通常更快；
• 求逆再相乘的方法耗时更长，仅适合多次复用A逆矩阵的场景。

内容的提问来源于stack exchange，提问作者C.Schroeder