python – 加速解决三角线性系统与numpy？ -清零世界

概述

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

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

我的问题：这个解决方案在python中比在Matlab中尝试相同时要慢得多(2倍).以下是一些时间实验：

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,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的速度.

我发现这个线程绊倒了numpy.linalg.solve和scipy.linalg.solve(和scipy的lu_solve等)之间的一些差异.我没有Enthought的基于MKL的Numpy / Scipy,但我希望我的发现能够以某种方式帮助你.

使用Numpy和Scipy的预构建二进制文件(32位,在Windows 7上运行)：

>当求解向量X时(即X为160×1),我看到numpy.linalg.solve和scipy.linalg.solve之间有显着差异. Scipy运行时是1.23x numpy,这是我认为实质的.
>然而,大部分差异似乎是由于scipy解决无效条目的检查.当将check_finite = False传递到scipy.linalg.solve中时,scipy的解决运行时是1.02x numpy.
> Scipy的解决使用破坏性更新,即overwrite_a = True,overwrite_b = True比numpy的解决(这是非破坏性的)稍快一些. Numpy的解决运行时是1.021x破坏性scipy.linalg.solve. scipy只是check_finite = False具有运行时1.04x的破坏性情况.总之,破坏性的scipy.linalg.solve比这些情况中的任何一个都要快一些.
>以上是一个向量X.如果我使X是一个宽阵列,具体是160乘以10000,与check_finite = False的scipy.linalg.solve基本上与check_finite = False一样快,overwrite_a = True,overwrite_b = True. Scipy的解决(没有任何特殊的关键字)运行时是1.09x这个“不安全”(check_finite = False)调用. Numpy的解决方案的运行时间为1.03x scipy对于这个阵列X的最快.
> scipy.linalg.solve_triangular在这两种情况下提供了显着的加速,但是您必须关闭输入检查,即传入check_finite = False.最快解决的运行时间分别为5.68x和1.76x solve_triangular,对于vector和array X,分别为check_finite = False.
具有破坏性计算的solve_triangular(overwrite_b = True)使您无法在check_finite = False之上加快速度(实际上对于阵列X的情况实际上有些伤害).
> I,ignoramus,以前不知道solve_triangular,并使用scipy.linalg.lu_solve作为三角求解器,即代替solve_triangular(cholS,X)做lu_solve((cholS,numpy.arange(160)),X)(都产生相同的答案).但是我发现以这种方式使用的lu_solve对于向量X的情况是运行时1.07x不安全的solve_triangular,而对于数组X的情况,其运行时为1.76x.我不知道为什么lu_solve对于数组X比向量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,bc),check_finite=False),"scipy.solve,infinite-ok"))
t.append(printer(timeit(lambda: lapper(nla.solve,"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,number=laps))

print "sla with tweaks"
t.append(printer(timeit(lambda: lapper(sla.solve,overwrite_a=True,overwrite_b=True,"scipy.solve destructive"))

print "Tri"
t.append(printer(timeit(lambda: lapper(sla.solve_triangular,"scipy.solve_triangular"))
t.append(printer(timeit(lambda: lapper(sla.solve_triangular,"scipy.solve_triangular,inf-ok"))
t.append(printer(timeit(lambda: lapper(sla.solve_triangular,"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,"LU"))

print "all times:"
print t

输出上述脚本的矢量情况,n = 1：

C
0.000739405 seconds,1000 laps,scipy.solve
0.000624746 seconds,scipy.solve,infinite-ok
0.000590003 seconds,numpy.solve
sla with tweaks
0.000608365 seconds,scipy.solve destructive
Tri
0.000208711 seconds,scipy.solve_triangular
9.38371e-05 seconds,scipy.solve_triangular,inf-ok
9.37682e-05 seconds,scipy.solve_triangular destructive
LU
0.000100215 seconds,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,infinite-ok
0.115569 seconds,numpy.solve
sla with tweaks
0.113122 seconds,scipy.solve destructive
Tri
0.0725959 seconds,scipy.solve_triangular
0.0634396 seconds,inf-ok
0.0638423 seconds,scipy.solve_triangular destructive
LU
0.1115 seconds,LU
all times:
[0.11898513112988955,0.11368747217793944,0.11556863916356903,0.11312182352918797,0.07259593807427585,0.0634396208970783,0.06384230931663318,0.11150022257648459]

请注意,上述Python脚本可以将其数组保存为Matlab .MAT数据文件.这是目前禁用的(如果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)

您将需要Mathworks文件交换：http://www.mathworks.com/matlabcentral/fileexchange/18798-timeit-benchmarking-function的时间函数.它产生以下输出：

matrix_time =
    0.0099989
vector_time =
   2.2487e-05

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

总结

以上是编程之家为你收集整理的python – 加速解决三角线性系统与numpy？全部内容，希望文章能够帮你解决python – 加速解决三角线性系统与numpy？所遇到的程序开发问题。

如果您也喜欢它,动动您的小指点个赞吧

python – 加速解决三角线性系统与numpy？

概述

总结

分类汇总

您的鼓励是对我最大的支持