Python 中带有三角形右侧的三角形线性系统

Question

我必须求解具有多个右侧的线性方程组，A*X=B，其中 A 和 B 都是（上）三角实数方阵。大小约为 200 x 200。在 python/numpy 中是否有快速方法？

我正在考虑循环遍历列： [编辑：按照评论中的要求包含一个 7x7 的示例。]

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

，但通常不鼓励对矩阵求逆。

Answer 1

为了回答我自己的问题，我找不到一个使用三角形结构并使用 BLAS3 运算的例程。我最终使用了问题中循环列方法的阻塞版本，一次处理一堆列。

X = np.zeros((n,n))
bs = 32 #block size.
for bst in range(0, n, bs):#bst = block start
    bsn = min(bst + bs, n)#block start next
    X[:bsn,bst:bsn] = sp.solve_triangular(
        A[:bsn,:bsn], B[:bsn,bst:bsn])

从好的方面来说，这确实使用了结构和 BLAS3 运算。缺点是你必须调整块大小。