J'ai une grande grille de coordonnées (vecteurs a et b), pour laquelle je génère et résolve une équation matricielle (10x10). Y a-t-il un moyen pour scipy.linalg.solve d'accepter l'entrée vectorielle? Jusqu'à présent, ma solution consistait à exécuter pendant des cycles sur les tableaux de coordonnées.

import time
import numpy as np
import scipy.linalg

N = 10
a = np.linspace(0, 1, 10**3)
b = np.linspace(0, 1, 2*10**3)
A = np.random.random((N, N))      # input matrix, not static

def f(a,b,n):     # array-filling function
   return a*b*n

def sol(x,y):   # matrix solver
   D = np.arange(0,N)
   B = f(x,y,D)**2 + f(x-1, y+1, D)      # source vector
   X = scipy.linalg.solve(A,B)
   return X    # output an N-size vector

start = time.time()

answer = np.zeros(shape=(a.size, b.size)) # predefine output array

for egg in range(a.size):             # an ugly double-for cycle
   for ham in range(b.size):
       aa = a[egg]
       bb = b[ham]
       answer[egg,ham] = sol(aa,bb)[0]

print time.time() - start
2
yevgeniy 22 juil. 2015 à 18:07

2 réponses

Meilleure réponse

@yevgeniy Vous avez raison, résoudre efficacement plusieurs systèmes linéaires indépendants A x = b avec scipy un peu délicat (en supposant un tableau A qui change à chaque itération).

Par exemple, voici un benchmark pour résoudre 1000 systèmes de la forme, A x = b , où A est une 10x10 matrice, et b a 10 vecteur d'élément. Étonnamment, l'approche consistant à mettre tout cela dans une matrice diagonale de bloc et à appeler scipy.linalg.solve une fois est en effet plus lente avec des matrices denses et clairsemées.

import numpy as np
from scipy.linalg import block_diag, solve
from scipy.sparse import block_diag as sp_block_diag
from scipy.sparse.linalg import spsolve

N = 10
M = 1000 # number of coordinates 
Ai = np.random.randn(N, N) # we can compute the inverse here,
# but let's assume that Ai are different matrices in the for loop loop
bi = np.random.randn(N)

%timeit [solve(Ai, bi) for el in range(M)]
# 10 loops, best of 3: 32.1 ms per loop

Afull = sp_block_diag([Ai]*M, format='csr')
bfull = np.tile(bi, M)

%timeit Afull = sp_block_diag([Ai]*M, format='csr')
%timeit spsolve(Afull, bfull)

# 1 loops, best of 3: 303 ms per loop
# 100 loops, best of 3: 5.55 ms per loop

Afull = block_diag(*[Ai]*M) 

%timeit Afull = block_diag(*[Ai]*M)
%timeit solve(Afull, bfull)

# 100 loops, best of 3: 14.1 ms per loop
# 1 loops, best of 3: 23.6 s per loop

La solution du système linéaire, avec des tableaux clairsemés est plus rapide, mais le temps de création de ce tableau diagonal de blocs est en fait très lent. Quant aux tableaux denses, ils sont simplement plus lents dans ce cas (et prennent beaucoup de RAM).

Il me manque peut-être quelque chose sur la façon de faire fonctionner efficacement cela avec des tableaux épars, mais si vous gardez les boucles for, il y a deux choses que vous pouvez faire pour les optimisations.

À partir de python pur, regardez le code source de { {X0}}: supprimez les tests inutiles et factorisez toutes les opérations répétées en dehors de vos boucles. Par exemple, en supposant que vos tableaux ne sont pas symétriques positifs, nous pourrions faire

from scipy.linalg import get_lapack_funcs

gesv, = get_lapack_funcs(('gesv',), (Ai, bi))

def solve_opt(A, b, gesv=gesv):
    # not sure if copying A and B is necessary, but just in case (faster if arrays are not copied)
    lu, piv, x, info = gesv(A.copy(), b.copy(), overwrite_a=False, overwrite_b=False)
    if info == 0:
        return x
    if info > 0:
        raise LinAlgError("singular matrix")
    raise ValueError('illegal value in %d-th argument of internal gesv|posv' % -info)

%timeit [solve(Ai, bi) for el in range(M)]
%timeit [solve_opt(Ai, bi) for el in range(M)]

# 10 loops, best of 3: 30.1 ms per loop
# 100 loops, best of 3: 3.77 ms per loop

Ce qui entraîne une accélération de 6,5x.

Si vous avez besoin de performances encore meilleures, vous devrez porter cette boucle for en Cython et interfacer les gesv fonctions BLAS directement en C, comme indiqué ici, ou mieux avec l'API Cython pour BLAS / LAPACK dans Scipy 0.16.

Edit : comme @Eelco Hoogendoorn l'a mentionné si votre matrice A est fixe, il existe une approche beaucoup plus simple et plus efficace.

1
rth 23 juil. 2015 à 15:19

Pour illustrer mon propos sur les ufuncs généralisés et la possibilité d'éliminer la boucle sur l'œuf et le jambon, considérez le morceau de code suivant:

import numpy as np
A = np.random.randn(4,4,10,10)
AI = np.linalg.inv(A)
#check that generalized ufuncs work as expected
I = np.einsum('xyij,xyjk->xyik', A, AI)
print np.allclose(I, np.eye(10)[np.newaxis, np.newaxis, :, :])
2
Eelco Hoogendoorn 23 juil. 2015 à 14:52