Multiplier efficacement les matrices creuses et lentes de Numpy / Scipy

Je suis en train de travailler pour mettre en œuvre l'équation suivante:

X =(Y.T * Y + Y.T * C * Y) ^ -1

Y est un (e x f) la matrice C (n x n) de diagonale; n est d'environ 300 k et f varie entre 100 et 200. Dans le cadre d'un processus d'optimisation de cette équation sera utilisé près de 100 millions de fois, donc elle doit être traitée très rapidement.

Y est initialisé aléatoirement, et C est une très sparse matrix avec seulement quelques numéros de la 300k sur la diagonale sera différent de 0.Depuis Numpy diagonal fonctions crée des matrices denses, j'ai créé C comme un sparse rse de la matrice. Mais lorsque j'essaie de résoudre la première partie de l'équation:

r = dot(C, Y)

L'ordinateur tombe en panne en raison des limites de la Mémoire. J'ai alors décidé d'essayer de convertir Y de csr_matrix et de faire la même opération:

r = dot(C, Ysparse)

et cette approche a pris 1.38 ms. Mais cette solution est un peu "délicat" depuis que je suis à l'aide d'une matrice creuse pour stocker dense, je me demande quelle est l'efficacité de cette vraiment.

Donc ma question est si est-il un moyen de multiplier les rares C et de la densité de Y, sans avoir à tourner en Y dans éparses et d'améliorer les performances? Si d'une certaine façon C pourrait être représenté comme diagonale dense sans consommer beaucoup de mémoire peut-être que cela conduirait à la performance très efficace mais je ne sais pas si c'est possible.

Je vous remercie de votre aide!

source d'informationauteur Willian Fuks

17

La raison pour laquelle le produit scalaire heurte à des problèmes de mémoire lors du calcul de r = point(C,Y), c'est parce numpy est dot de la fonction n'a pas de support natif pour la manipulation de matrices creuses. Ce qui se passe est numpy pense de la matrice creuse C comme un objet python, et non pas un tableau numpy. Si vous inspectez sur petite échelle, vous pouvez voir le problème de première main:
```
>>> from numpy import dot, array
>>> from scipy import sparse
>>> Y = array([[1,2],[3,4]])
>>> C = sparse.csr_matrix(array([[1,0], [0,2]]))
>>> dot(C,Y)
array([[  (0, 0)    1
  (1, 1)    2,   (0, 0) 2
  (1, 1)    4],
  [  (0, 0) 3
  (1, 1)    6,   (0, 0) 4
  (1, 1)    8]], dtype=object)
```
Clairement ci-dessus n'est pas le résultat qui vous intéresse. Au lieu de ce que vous voulez faire est de calculer à l'aide de scipy éparses.csr_matrix.dot de la fonction:
```
r = sparse.csr_matrix.dot(C, Y)
```
ou de manière plus compacte
```
r = C.dot(Y)
```

Essayer:

import numpy as np
from scipy import sparse

f = 100
n = 300000

Y = np.random.rand(n, f)
Cdiag = np.random.rand(n) # diagonal of C
Cdiag[np.random.rand(n) < 0.99] = 0

# Compute Y.T * C * Y, skipping zero elements
mask = np.flatnonzero(Cdiag)
Cskip = Cdiag[mask]

def ytcy_fast(Y):
    Yskip = Y[mask,:]
    CY = Cskip[:,None] * Yskip  # broadcasting
    return Yskip.T.dot(CY)

%timeit ytcy_fast(Y)

# For comparison: all-sparse matrices
C_sparse = sparse.spdiags([Cdiag], [0], n, n)
Y_sparse = sparse.csr_matrix(Y)
%timeit Y_sparse.T.dot(C_sparse * Y_sparse)

Mes horaires:

In [59]: %timeit ytcy_fast(Y)
100 loops, best of 3: 16.1 ms per loop

In [18]: %timeit Y_sparse.T.dot(C_sparse * Y_sparse)
1 loops, best of 3: 282 ms per loop

2

D'abord, êtes-vous vraiment sûr que vous devez effectuer une complète inversion de matrice dans votre problème ? La plupart du temps, on n'a vraiment besoin de calculer pour x = A^-1 y, qui est beaucoup plus facile de problème à résoudre.

Si c'est vraiment le cas, je considère le calcul d'une approximation de l'inverse de la matrice au lieu de la pleine inversion de matrice. Depuis inversion de matrice est vraiment coûteux. Voir, par exemple, la Algorithme de Lanczos efficace pour une approximation de l'inverse de la matrice. Le rapprochement peut être stocké peu comme un bonus. De Plus, il ne nécessite que de la matrice-vecteur opérations de sorte que vous n'avez même pas à stocker la matrice complète à l'inverse.

Comme une alternative, à l'aide de pyoperators, vous pouvez également l'utiliser pour .todense méthode de calcul de la matrice inverse à l'aide efficace de la matrice vecteur des opérations. Il y a un spécial éparses conteneur pour la diagonale des matrices.

Pour une mise en œuvre de l'algorithme de Lanczos, vous pouvez jeter un oeil à pyoperators (disclaimer: je suis l'un des co-auteur de ce morceau de logiciel).

Vous devez vous connecter pour publier un commentaire.