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 で実行) 用のビルド済みバイナリを使用すると、次のようになります。
ベクトルを解く場合(つまり、160 x 1)と の間には大きな違いがnumpy.linalg.solveあります。Scipy のランタイムは numpy の 1.23 倍で、かなりのものだと思います。scipy.linalg.solveXX
ただし、ほとんどの違いは、scipy によるsolve無効なエントリのチェックによるものと思われます。check_finite=Falsescipy.linalg.solve に渡すと、scipy の実行solve時間は numpy の 1.02 倍になります。
破壊的な更新を使用した Scipy の解決、つまり withoverwrite_a=True, overwrite_b=Trueは、numpy の解決 (非破壊的) よりもわずかに高速です。Numpy のソルブ ランタイムは 1.021x 破壊的な scipy.linalg.solve です。Scipycheck_finite=Falseのランタイムは 1.04x の破壊的なケースです。要約すると、destructivescipy.linalg.solveはこれらのケースのいずれよりもわずかに高速です。
上記は vector 用Xです。X幅の広い配列、具体的には 160 x 10000を作成すると、 scipy.linalg.solvewithcheck_finite=Falseは本質的に with と同じくらい高速check_finite=False, overwrite_a=True, overwrite_b=Trueです。Scipy のsolve(特別なキーワードなしの) ランタイムは、この「安全でない」( check_finite=False) 呼び出しの 1.09 倍です。この配列の場合、 Numpy のsolveランタイムは 1.03x scipy で最速です。X
scipy.linalg.solve_triangularどちらの場合でも大幅なスピードアップを実現しますが、入力チェックをオフにする必要があります。つまり、 を渡しcheck_finite=Falseます。solve_triangular最速の解の実行時間は、 vector と arrayXで、それぞれ5.68x と 1.76x でしたcheck_finite=False。
solve_triangular破壊的な計算 ( overwrite_b=True) を使用すると、速度が向上しませんcheck_finite=False(実際には、配列の場合はわずかに問題がありXます)。
私、無知な私は、以前は知らず、三角ソルバーとしてsolve_triangular使用していました。しかし、この方法で使用すると、ベクトルの場合はランタイム 1.07x が安全ではないことがわかりましたが、配列の場合はランタイムが 1.76x でした。vector と比較して array の方がはるかに遅い理由はわかりませんが、レッスンは（無限チェックなしで）使用することです。scipy.linalg.lu_solvesolve_triangular(cholS, X)lu_solve((cholS, numpy.arange(160)), X)lu_solvesolve_triangularXXlu_solveXXsolve_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)
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