我的大部分線性代數研究都是在計算機之前（或者至少在 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)
