なぜscipy.linalgのメソッドは9x9行列で遅く実行されますか？

Question

私は、線形システムを解くためのさまざまな方法のランタイムを比較しています。私は奇妙なパターンを見つけました。私がテストしているソリューションの方法は、la.solve(),la.inv()、la.lu_factor_solve()です。

import scipy.linalg as la 
import numpy as np 
from time import time 
from matplotlib import pyplot as plt 

N = 20 # up to NxN matrix 
T = 100 # run T times 
inv_time, solve_time = [[] for _ in range(20)], [[] for _ in range(20)] 
lu_factor_solve, lu_just_solve = [[] for _ in range(20)], [[] for _ in range(20)] 
for _ in range(T): 
    for n in range(1, N + 1): 
     A = np.random.rand(n, n) 
     b = np.random.rand(n, 1) 
     np.dot(la.inv(A), b) # the first time through is always slow, 
     la.solve(A, b)  # so we run it once to get it out of the way 
     start = time() 
     np.dot(la.inv(A), b) 
     end = time() 
     inv_time[n - 1].append(end - start) 

     start = time() 
     la.solve(A, b) 
     end = time() 
     solve_time[n - 1].append(end - start) 

     start = time() 
     la.lu_solve(la.lu_factor(A), b) 
     end = time() 
     lu_factor_solve[n - 1].append(end - start) 

     temp = la.lu_factor(A) 
     start = time() 
     la.lu_solve(temp, b) 
     end = time() 
     lu_just_solve[n - 1].append(end - start) 
inv_time = np.mean(np.array(inv_time), axis=1) 
solve_time = np.mean(np.array(solve_time), axis=1) 
lu_factor_solve = np.mean(np.array(lu_factor_solve), axis=1) 
lu_just_solve = np.mean(np.array(lu_just_solve), axis=1) 
# do some plots 
plt.plot(range(1, N + 1), inv_time, '-o', label='by inverse') 
plt.plot(range(1, N + 1), solve_time, '-o', label='by la.solve()') 
plt.plot(range(1, N + 1), lu_factor_solve, '-o', label='by lu factor solve') 
plt.yscale('log') 
plt.plot(range(1, N + 1), lu_just_solve, '-o', label='just la.lu_solve()') 
plt.legend() 
plt.show() 

はこれらがランダム行列Aと1から20までnの値に対してA = np.random.rand(n, n)とb = np.random.rand(n, 1)によって生成列ベクトルb上で実行されている私は、プログラムを実行するたびに、それは9x9の行列に解を求めることを、私を見つけました以下に示すように、他のサイズの行列よりもはるかに低速です。赤い線は、la.lu_solve()の実行に必要な時間を示しています。ここでT = 1と結果のグラフである：

そして、ここではT = 100と結果は以下のとおりです。

これは、このサイズのためのいくつかの最適化ができない、9x9の行列について固有の何かをしなければならないんマトリックス、または別の何か？

Answer 1

コード例をperfplot（私の小規模なプロジェクト、本質的にtimeitのラッパー）で試してみましたが、そのような特殊性は見つかりませんでした。一つはしかし、レベル1キャッシュの枯渇を認識することができます

これはnumpyの1.13.3とscipyのダウンロード1.0.0rc1です。

を実行しているときにコードサンプルにスパイクが発生する理由は、マシンの他の部分との小さな干渉が大きくなる可能性があるためです。これは、ブラウザのJSエンジンから深いスリープ状態から復帰するCPUまでの定期的なタスクから何かまで可能性があります。 n=9でそのスパイクが見えるという事実は偶然です。私自身は、テストを1回だけ実行すると、5〜20の間の任意のスパイクを再現できます。おそらく重くSIMDとキャッシング用に最適化されたBLAS/LAPACK、に関連

import numpy 
import perfplot 
import scipy.linalg as la 


def solve(data): 
    A, b = data 
    return la.solve(A, b) 


def dot_inv(data): 
    A, b = data 
    return numpy.dot(la.inv(A), b) 


def lu_solve(data): 
    A, b = data 
    return la.lu_solve(la.lu_factor(A), b) 


perfplot.show(
    setup=lambda n: (numpy.random.rand(n, n), numpy.random.rand(n)), 
    kernels=[solve, dot_inv, lu_solve], 
    n_range=range(1, 40), 
    )