Pourquoi les méthodes de scipy.linalg sont-elles plus lentes sur les matrices 9x9?

Question

Je compare des runtimes de différentes manières pour résoudre des systèmes linéaires, et j'ai trouvé un modèle impair. Les méthodes de solution que je suis en train de tester sont la.solve(), la.inv() et 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() 

Ces sont exécutés sur des matrices aléatoires A et des vecteurs de colonne b produit par A = np.random.rand(n, n) et b = np.random.rand(n, 1) pour les valeurs de n de 1 à 20. J'ai trouvé que chaque fois que je lance le programme, il trouve des solutions aux matrices 9x9 beaucoup plus lent que pour les matrices d'autres tailles, comme sur la photo ci-dessous. La ligne rouge indique le temps nécessaire pour effectuer la.lu_solve(). Voici un graphique des résultats avec T = 1:

Et voici les résultats avec T = 100:

Est-ce que cela a à voir avec quelque chose d'inhérent à propos de matrices 9x9, une optimisation non disponible pour cette taille de matrice, ou quelque chose de différent?

Answer 1

J'ai essayé vos exemples de code avec perfplot (un petit projet à moi, essentiellement un wrapper autour de timeit) et je n'ai pas trouvé de telles particularités. Un peut reconnaître l'épuisement du niveau 1-cache cependant:

est avec NumPy 1.13.3 et SciPy version 1.0.0RC1.

La raison pour laquelle il y a des pics dans votre exemple de code lorsque seulement l'exécution des tests une fois est qu'ils courent si vite que les petites interférences avec d'autres parties de la machine peuvent devenir importants. Il peut s'agir de tâches récurrentes provenant du moteur JS de votre navigateur ou de vos processeurs se réveillant d'états de veille plus profonds. Le fait que vous voyez ces pointes à n=9 est une coïncidence. Je peux moi-même reproduire des pics aléatoires entre 5 et 20 lorsque je ne fais qu'une seule fois les tests.

---

Code pour reproduire l'intrigue:

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), 
    )