solve
TL;DR:当你有一个三角系统时,不要使用 numpy 或 scipy scipy.linalg.solve_triangular
,至少使用check_finite=False
关键字参数来获得快速和非破坏性的解决方案。
numpy.linalg.solve
在偶然发现和scipy.linalg.solve
(和scipy's lu_solve
等)之间的一些差异后,我发现了这个线程。我没有 Enthought 的基于 MKL 的 Numpy/Scipy,但我希望我的发现可以在某种程度上对您有所帮助。
使用 Numpy 和 Scipy 的预构建二进制文件(32 位,在 Windows 7 上运行):
我看到求解向量(即160 x 1)之间numpy.linalg.solve
和求解向量时存在显着差异。Scipy 运行时是 numpy 的 1.23 倍,我认为这是可观的。scipy.linalg.solve
X
X
但是,大部分差异似乎是由于 scipysolve
检查无效条目所致。当传入check_finite=False
scipy.linalg.solve 时,scipy 的solve
运行时间是 1.02x numpy 的。
Scipy 使用破坏性更新的求解,即 withoverwrite_a=True, overwrite_b=True
比 numpy 的求解(非破坏性)稍快。Numpy 的求解运行时是 1.021x 破坏性 scipy.linalg.solve。Scipy 的check_finite=False
运行时间是破坏性案例的 1.04 倍。总之,破坏性scipy.linalg.solve
比这两种情况都快得多。
以上是针对向量的X
。如果我制作X
一个宽数组,特别是 160 x 10000,则scipy.linalg.solve
withcheck_finite=False
基本上与check_finite=False, overwrite_a=True, overwrite_b=True
. Scipy 的solve
(没有任何特殊关键字)运行时是这个“不安全”(check_finite=False
)调用的 1.09 倍。对于这种数组情况, Numpy 的solve
运行时间是 scipy 最快的 1.03 倍。X
scipy.linalg.solve_triangular
在这两种情况下都提供了显着的加速,但您必须关闭输入检查,即传入check_finite=False
. solve_triangular
对于 vector 和 array ,最快求解的运行时间分别为 5.68x 和 1.76x X
,其中check_finite=False
.
solve_triangular
破坏性计算(overwrite_b=True
_check_finite=False
X
我,无知,以前不知道solve_triangular
并scipy.linalg.lu_solve
用作三角求解器,即,而不是solve_triangular(cholS, X)
做lu_solve((cholS, numpy.arange(160)), X)
(两者都产生相同的答案)。但是我发现lu_solve
以这种方式使用solve_triangular
的向量X
情况下的运行时间为 1.07 倍不安全,而数组X
情况下的运行时间为 1.76 倍。我不确定为什么lu_solve
arrayX
比 vector慢得多X
,但教训是使用solve_triangular
(没有无限检查)。
将数据复制到 Fortran 格式似乎一点也不重要。转换为numpy.matrix
.
我不妨将我的非 MKL Python 库与单线程 ( maxNumCompThreads=1
) Matlab 2013a 进行比较。上面最快的 Python 实现对于向量X
情况的运行时间延长了 4.5 倍,对于胖矩阵X
情况的运行时间延长了 6.3 倍。
但是,这是我用来对这些进行基准测试的 Python 脚本,也许有 MKL 加速的 Numpy/Scipy 的人可以发布他们的数字。请注意,我只是注释掉该行n = 10000
以禁用胖矩阵X
案例并执行n=1
向量案例。(对不起。)
import scipy.linalg as sla
import numpy.linalg as nla
from numpy.random import RandomState
from timeit import timeit
import numpy as np
RNG = RandomState(69)
m=160
n=1
#n=10000
Ac = RNG.randn(m,m)
if 1:
Ac = np.triu(Ac)
bc = RNG.randn(m,n)
Af = Ac.copy("F")
bf = bc.copy("F")
if 0: # Save to Matlab format
import scipy.io as io
io.savemat("b_%d.mat"%(n,), dict(A=Ac, b=bc))
import sys
sys.exit(0)
def lapper(fn, source, **kwargs):
Alocal = source[0].copy()
blocal = source[1].copy()
fn(Alocal, blocal,**kwargs)
laps = (1000 if n<=1 else 100)
def printer(t, s=''):
print ("%g seconds, %d laps, " % (t/float(laps), laps)) + s
return t/float(laps)
t=[]
print "C"
t.append(printer(timeit(lambda: lapper(sla.solve, (Ac,bc)), number=laps),
"scipy.solve"))
t.append(printer(timeit(lambda: lapper(sla.solve, (Ac,bc), check_finite=False),
number=laps), "scipy.solve, infinite-ok"))
t.append(printer(timeit(lambda: lapper(nla.solve, (Ac,bc)), number=laps),
"numpy.solve"))
#print "F" # Doesn't seem to matter
#printer(timeit(lambda: lapper(sla.solve, (Af,bf)), number=laps))
#printer(timeit(lambda: lapper(nla.solve, (Af,bf)), number=laps))
print "sla with tweaks"
t.append(printer(timeit(lambda: lapper(sla.solve, (Ac,bc), overwrite_a=True,
overwrite_b=True, check_finite=False),
number=laps), "scipy.solve destructive"))
print "Tri"
t.append(printer(timeit(lambda: lapper(sla.solve_triangular, (Ac,bc)),
number=laps), "scipy.solve_triangular"))
t.append(printer(timeit(lambda: lapper(sla.solve_triangular, (Ac,bc),
check_finite=False), number=laps),
"scipy.solve_triangular, inf-ok"))
t.append(printer(timeit(lambda: lapper(sla.solve_triangular, (Ac,bc),
overwrite_b=True, check_finite=False),
number=laps), "scipy.solve_triangular destructive"))
print "LU"
piv = np.arange(m)
t.append(printer(timeit(lambda: lapper(
lambda X,b: sla.lu_solve((X, piv),b,check_finite=False),
(Ac,bc)), number=laps), "LU"))
print "all times:"
print t
上述向量案例的脚本输出n=1
:
C
0.000739405 seconds, 1000 laps, scipy.solve
0.000624746 seconds, 1000 laps, scipy.solve, infinite-ok
0.000590003 seconds, 1000 laps, numpy.solve
sla with tweaks
0.000608365 seconds, 1000 laps, scipy.solve destructive
Tri
0.000208711 seconds, 1000 laps, scipy.solve_triangular
9.38371e-05 seconds, 1000 laps, scipy.solve_triangular, inf-ok
9.37682e-05 seconds, 1000 laps, scipy.solve_triangular destructive
LU
0.000100215 seconds, 1000 laps, LU
all times:
[0.0007394047886284343, 0.00062474593940593, 0.0005900030818282472, 0.0006083650710913095, 0.00020871054023307778, 9.383710445114923e-05, 9.37682389063692e-05, 0.00010021534750467032]
矩阵案例的上述脚本的输出n=10000
:
C
0.118985 seconds, 100 laps, scipy.solve
0.113687 seconds, 100 laps, scipy.solve, infinite-ok
0.115569 seconds, 100 laps, numpy.solve
sla with tweaks
0.113122 seconds, 100 laps, scipy.solve destructive
Tri
0.0725959 seconds, 100 laps, scipy.solve_triangular
0.0634396 seconds, 100 laps, scipy.solve_triangular, inf-ok
0.0638423 seconds, 100 laps, scipy.solve_triangular destructive
LU
0.1115 seconds, 100 laps, LU
all times:
[0.11898513112988955, 0.11368747217793944, 0.11556863916356903, 0.11312182352918797, 0.07259593807427585, 0.0634396208970783, 0.06384230931663318, 0.11150022257648459]
请注意,上述 Python 脚本可以将其数组保存为 Matlab .MAT 数据文件。这目前被禁用(if 0
,抱歉),但如果启用,您可以在完全相同的数据上测试 Matlab 的速度。这是 Matlab 的计时脚本:
clear
q = load('b_10000.mat');
A=q.A;
b=q.b;
clear q
matrix_time = timeit(@() A\b)
q = load('b_1.mat');
A=q.A;
b=q.b;
clear q
vector_time = timeit(@() A\b)
您需要timeit
Mathworks File Exchange 中的函数:http: //www.mathworks.com/matlabcentral/fileexchange/18798-timeit-benchmarking-function。这会产生以下输出:
matrix_time =
0.0099989
vector_time =
2.2487e-05
solve
这个实证分析的结果是,至少在 Python 中,当你有一个三角形系统时,不要使用 numpy 或 scipy scipy.linalg.solve_triangular
,至少使用check_finite=False
关键字参数来获得快速和非破坏性的解决方案。