为什么scipy.linalg.LU在反复求解Ax =b时如此缓慢？

Question

传统智慧指出，如果您使用相同的A和不同的b多次求解Ax = b，则应该对LU使用LU分解。如果我使用p, l, u = scipy.linalg.lu(A)并在一个循环中多次求解

x = scipy.linalg.solve(l, p.T@b)
x = scipy.linalg.solve(u, x)

这最终比仅仅使用

x = scipy.linalg.solve(A,b)

在循环中。scipy.linalg.solve()不是针对使用向前和向后替换的上对角线系统和下对角线系统进行了优化吗？或者，有没有可能存在一些编译技巧，其中python认识到它可以为scipy.linalg.solve部分进行LU分解？

我知道在scipy中有linalg.lu_factor和linalg.lu_solve例程，但我想远离它们，因为这应该是一个教学示例。

Answer 1

我的大多数线性代数研究是在计算机之前(或者至少在MATLAB/python之前)进行的。但我能读懂文档。

In [29]: from scipy import linalg as la

从来自la.lu的示例数组开始

In [30]: A = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]])
In [31]: p, l, u = la.lu(A)
In [32]: p
Out[32]: 
array([[0., 1., 0., 0.],
       [0., 0., 0., 1.],
       [1., 0., 0., 0.],
       [0., 0., 1., 0.]])
In [33]: l
Out[33]: 
array([[ 1.        ,  0.        ,  0.        ,  0.        ],
       [ 0.28571429,  1.        ,  0.        ,  0.        ],
       [ 0.71428571,  0.12      ,  1.        ,  0.        ],
       [ 0.71428571, -0.44      , -0.46153846,  1.        ]])
In [34]: u
Out[34]: 
array([[ 7.        ,  5.        ,  6.        ,  6.        ],
       [ 0.        ,  3.57142857,  6.28571429,  5.28571429],
       [ 0.        ,  0.        , -1.04      ,  3.08      ],
       [ 0.        ,  0.        ,  0.        ,  7.46153846]])

In [42]: b=np.arange(4)
In [43]: la.solve(A,b)
Out[43]: array([-0.21649485,  2.54639175, -1.54639175,  0.01030928])
In [44]: timeit la.solve(A,b)
43.5 µs ± 88.5 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

我看到了la.solve_triangular的文档。经过一些试验和错误，我得到了：

In [46]: la.solve_triangular(u,la.solve_triangular(l,p.T@b, lower=True))
Out[46]: array([-0.21649485,  2.54639175, -1.54639175,  0.01030928])

并为其计时：

In [47]: timeit la.solve_triangular(u,la.solve_triangular(l,p.T@b, lower=True))
83 µs ± 2.6 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

因此，双重使用solve_trianglar比使用一个solve慢，但比使用不知道它的数组是三角形的solve要快。

In [48]: la.solve(u,la.solve(l,p.T@b))
Out[48]: array([-0.21649485,  2.54639175, -1.54639175,  0.01030928])
In [49]: timeit la.solve(u,la.solve(l,p.T@b))
137 µs ± 342 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

我不知道这些计算将如何扩展。

测试@Warren的lu_solve想法(在已删除的答案中) https://stackoverflow.com/a/64473976/901925

In [50]: lu_and_piv = la.lu_factor(A)
In [51]: lu_and_piv
Out[51]: 
(array([[ 7.        ,  5.        ,  6.        ,  6.        ],
        [ 0.28571429,  3.57142857,  6.28571429,  5.28571429],
        [ 0.71428571,  0.12      , -1.04      ,  3.08      ],
        [ 0.71428571, -0.44      , -0.46153846,  7.46153846]]),
 array([2, 2, 3, 3], dtype=int32))
In [52]: la.lu_solve(lu_and_piv, b)
Out[52]: array([-0.21649485,  2.54639175, -1.54639175,  0.01030928])
In [53]: timeit la.lu_solve(lu_and_piv, b)
7.47 µs ± 14.3 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)