加快求解numpy的三角形线性系统？

Question

我有一个方阵S(160 x 160)和一个巨大的矩阵X(160 x 250000).两者都是密集的numpy数组.

我的目标:找到Q,使Q = inv(chol(S))*X,其中chol(S)是S的较低的cholesky分解.

当然,一个简单的解决方案是

cholS = scipy.linalg.cholesky( S, lower=True)
scipy.linalg.solve( cholS, X )

我的问题:这个解决方案是明显较慢(> 2倍)的蟒蛇比当我尝试在Matlab相同.以下是一些计时实验:

timeit np.linalg.solve( cholS, X)
1 loops, best of 3: 1.63 s per loop

timeit scipy.linalg.solve_triangular( cholS, X, lower=True)
1 loops, best of 3: 2.19 s per loop

timeit scipy.linalg.solve( cholS, X)
1 loops, best of 3: 2.81 s per loop

[matlab]
cholS \ X
0.675 s

[matlab using only one thread via -singleCompThread]
cholS \ X
1.26 s

基本上,我想知道:(1)我可以在python中达到类似Matlab的速度吗？(2)为什么scipy版本这么慢？

解算器应该能够利用chol(S)是三角形的事实.但是,使用numpy.linalg.solve()比scipy.linalg.solve_triangular()更快,即使numpy调用根本不使用三角形结构.是什么赋予了？当我的矩阵是三角形时,matlab求解器似乎会自动检测,但是python不能.

我很乐意使用BLAS/LAPACK例程的自定义调用来解决三角形线性系统,但我真的不想自己编写代码.

作为参考,我使用scipy版本11.0和Enthought python发行版(使用英特尔的MKL库进行矢量化),所以我认为我应该能够达到类似Matlab的速度.

Answer 1

TL; DR:solve当你有一个三角形系统时,不要使用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上运行):

1. 我看到之间的差异显著numpy.linalg.solve和scipy.linalg.solve用于向量求解时X(即,X是由160 1).Scipy运行时是1.23x numpy,这是我认为实质性的.

2. 但是,大多数差异似乎是由于scipy solve检查无效条目.当传入check_finite=Falsescipy.linalg.solve时,scipy的solve运行时间为1.02x numpy.

3. Scipy使用破坏性更新来解决,即使用overwrite_a=True, overwrite_b=True比numpy的解决(非破坏性)稍微快一点.Numpy的解决运行时是1.021x破坏性的scipy.linalg.solve.Scipy只是check_finite=False运行时1.04x的破坏性案例.总之,破坏性scipy.linalg.solve比这两种情况都要快一些.

4. 以上是矢量X.如果我制作X一个宽阵列,特别是160乘10000,基本上scipy.linalg.solve和check_finite=False它一样快check_finite=False, overwrite_a=True, overwrite_b=True.Scipy solve(没有任何特殊关键字)运行时是1.09x这个"不安全"(check_finite=False)调用.solve对于这个阵列,Numpy的运行时间为1.03x scipy X.

5. scipy.linalg.solve_triangular在这两种情况下都能提供显着的加速,但你必须关闭输入检查,即传入check_finite=False.最快求解的运行时间分别为5.68x和1.76x solve_triangular,对于矢量和数组X,分别为check_finite=False.

6. solve_triangular使用破坏性计算(overwrite_b=True)不会给你带来加速check_finite=False(并且对于数组X情况实际上会受到轻微伤害).

7. 我,无知,以前没有意识到solve_triangular并且正在使用scipy.linalg.lu_solve三角形求解器,即代替solve_triangular(cholS, X)做lu_solve((cholS, numpy.arange(160)), X)(两者都产生相同的答案).但我发现lu_solve以这种方式使用的运行时1.07x solve_triangular对于矢量X情况不安全,而对于数组X情况,其运行时间为1.76x .与矢量相比,我不确定为什么lu_solve数组的速度要慢得多,但教训是使用(没有无限的检查).XXsolve_triangular

8. 将数据复制到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)

您将需要timeitMathworks File Exchange中的功能:http://www.mathworks.com/matlabcentral/fileexchange/18798-timeit-benchmarking-function.这会产生以下输出:

matrix_time =
    0.0099989
vector_time =
   2.2487e-05

这个实证分析的结果是,至少在Python中,solve当你有一个三角形系统时,不要使用numpy或scipy ,只需使用scipy.linalg.solve_triangular至少check_finite=False关键字参数来获得快速和非破坏性的解决方案.

Answer 2

为什么不直接使用方程：Q = inv(chol(S)) * X，这是我的测试：

import scipy.linalg
import numpy as np

N = 160
M = 100000
S = np.random.randn(N, N)
B = np.random.randn(N, M)
S = np.dot(S, S.T)

cS = scipy.linalg.cholesky(S, lower=True)
Y1 = scipy.linalg.solve(cS, B)
icS = scipy.linalg.inv(cS)
Y2 = np.dot(icS, B)

np.allclose(Y1, Y2)

输出：

True

这是时间测试：

%time scipy.linalg.solve(cholS, B)
%time np.linalg.solve(cholS, B)
%time scipy.linalg.solve_triangular(cholS, B, lower=True)
%time ics=scipy.linalg.inv(cS);np.dot(ics, B)

输出：

CPU times: user 2.07 s, sys: 0.00 s, total: 2.07 s
Wall time: 2.08 s
CPU times: user 1.93 s, sys: 0.00 s, total: 1.93 s
Wall time: 1.92 s
CPU times: user 1.12 s, sys: 0.00 s, total: 1.12 s
Wall time: 1.13 s
CPU times: user 0.71 s, sys: 0.00 s, total: 0.71 s
Wall time: 0.72 s

我不知道为什么scipy.linalg.solve_triangular比你的系统慢numpy.linalg.solve，但该inv版本是最快的。