为什么 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比 1 慢solve，但比使用solve不知道其数组是三角形的 a 更快。

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)

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

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)