加快解決與numpy的三角形線性系統？

Question

我具有正方形矩陣S（160×160），和一個巨大的矩陣X（160×25萬）。兩者都是密集的numpy陣列。

我的目標：求Q，使得Q = INV（CHOL（S））* X，其中CHOL（S）是S的下Cholesky因式分解

自然地，一個簡單的解決方案是

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

我的問題：這個解決方案是在python 明顯較慢（> 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）我可以達到類似Matlab的速度在Python中？和（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

爲什麼不直接使用公式：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版本是最快的。

Answer 2

幾件事情嘗試：

• X = X.copy('F')＃用FORTRAN順序排列，這樣避免了複製

• Y = solve_triangular(cholS, X, overwrite_b=True)＃避免再次複製，但X

垃圾內容

• Y = solve_triangular(cholS, X, check_finite=False)＃scipy> = 0.12 only ---但似乎對速度沒有太大影響...

對於這兩者，它應該與直接調用沒有緩衝區副本的MKL等效。

我不能重現與np.linalg.solve和scipy.linalg.solve具有不同的速度問題---與我擁有的BLAS + LAPACK組合，兩者看起來都是相同的速度。

Answer 3

TL; DR：不要使用numpy的公司或SciPy的的solve當你有一個三角形的系統，只需使用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是1乘160）。 Scipy運行時是1.23x numpy的，這是我認爲很重要的。

2. 但是，大部分差異似乎是由於scipy的solve檢查無效條目。當傳遞check_finite=False到scipy.linalg.solve時，scipy的solve運行時是1.02x numpy的。

3. Scipy的解決方案使用破壞性更新，即overwrite_a=True, overwrite_b=True比numpy的解決方案（非破壞性）要快一些。 Numpy的求解運行時是1.021x破壞性的scipy.linalg.solve。只有check_finite=False的Scipy具有運行時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）調用。 Numpy的solve擁有運行時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. I，無知，是以前不知道的solve_triangular並使用scipy.linalg.lu_solve作爲三角解算器是，即，代替solve_triangular(cholS, X)做lu_solve((cholS, numpy.arange(160)), X)（都產生相同的答案）。但我發現lu_solve這種方式使用運行時1.07x不安全solve_triangular爲矢量X的情況下，而它的運行時爲1.76x的數組X的情況。我不知道爲什麼lu_solve對於X數組比慢很多，但是我們的教訓是使用solve_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：上述腳本爲矩陣的情況下n=10000

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] 

輸出：

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) 

你需要從Mathworks的文件交換的timeit功能：http://www.mathworks.com/matlabcentral/fileexchange/18798-timeit-benchmarking-function。這將產生以下的輸出：

matrix_time = 
    0.0099989 
vector_time = 
    2.2487e-05 

這種實證分析，在Python，至少，不要當你有一個三角形的系統使用numpy的公司或SciPy的的solve的結果，只是用scipy.linalg.solve_triangular與至少check_finite=False關鍵詞爭取快速和無損解決方案。