numpy: compute x.T * x pour une grande matrice

Question

Dans numpy, ce qui est la façon la plus efficace de calculer x.T * x, où x est un grand (200 000 x 1000) matrice float32 dense et .T est l'opérateur de transposition?

Pour éviter tout doute, le résultat est de 1000 x 1000.

modifier : Dans ma première question, je déclare que np.dot(x.T, x) prenait des heures. Il est avéré que j'avais un peu sournoise NaNs dans la matrice, et pour une raison quelconque qui a été complètement tuer les performances de np.dot (toutes idées pour expliquer pourquoi?) Ceci est maintenant résolu, mais la question se tient d'origine.

Answer 1

Cela peut ne pas être la réponse que vous cherchez, mais une façon d'accélérer est d'utiliser considérablement un gpu au lieu de votre cpu. Si vous avez une carte graphique puissante décemment autour, il va surperformer votre cpu un jour, même si votre système est très bien réglé.

Pour belle intégration avec numpy, vous pouvez utiliser Théano (si votre carte graphique est faite par nvidia). Le calcul dans les pistes de code suivant pour moi en quelques secondes (bien que j'ai une carte graphique très puissante):

$ THEANO_FLAGS=device=gpu0 python
Python 2.6.5 (r265:79063, Apr 16 2010, 13:57:41) 
[GCC 4.4.3] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import theano
Using gpu device 0: GeForce GTX 480
>>> from theano import tensor as T
>>> import numpy
>>> x = numpy.ones((200000, 1000), dtype=numpy.float32)
>>> m = T.matrix() 
>>> mTm = T.dot(m.T, m)
>>> f = theano.function([m], mTm)
>>> f(x)
array([[ 200000.,  200000.,  200000., ...,  200000.,  200000.,  200000.],
       [ 200000.,  200000.,  200000., ...,  200000.,  200000.,  200000.],
       [ 200000.,  200000.,  200000., ...,  200000.,  200000.,  200000.],
       ..., 
       [ 200000.,  200000.,  200000., ...,  200000.,  200000.,  200000.],
       [ 200000.,  200000.,  200000., ...,  200000.,  200000.,  200000.],
       [ 200000.,  200000.,  200000., ...,  200000.,  200000.,  200000.]], dtype=float32)
>>> r = f(x)
>>> r.shape
(1000, 1000)

J'allais attendre pour savoir combien de temps >>> numpy.dot(x.T, x) a pris à titre de comparaison, mais je me suis ennuyé ...

Vous pouvez également essayer PyCuda ou PyOpenCL (si vous ne disposez pas d'une carte graphique nvidia), bien que je ne sais pas si leur soutien est numpy aussi simple.

Answer 2

Tout d'abord, assurez-vous d'utiliser un blas optimisé / lapack, cela peut faire une énorme différence (jusqu'à un ordre de grandeur). Si vous utilisez un enfilée ATLAS, par exemple, il utilisera tous vos cœurs relativement efficace (vous devez utiliser un récent ATLAS, cependant, et la compilation ATLAS est un PITA).

Quant à savoir pourquoi le tout Nan fait: c'est à peu près inévitable, la manipulation NaN est beaucoup plus lent que float « normal » au niveau du CPU: http://www.cygnus-software.com/papers/x86andinfinity.html . Cela dépend du modèle de CPU, quel genre de jeu d'instructions que vous utilisez, et bien sûr les algorithmes / mise en œuvre que vous utilisez.

Answer 3

hmm, x est d'environ 800 Mo, en supposant qu'il a besoin de la même chose pour le résultat, êtes-vous sûr que vous avez suffisamment de mémoire physique et ce n'est pas échanger?

autre que celui, numpy doit utiliser une fonction BLAS, et même si la bibliothèque par défaut qui utilise numpy peut être, il devrait relativement lent travail correct pour cette taille.

modifier

import numpy as npy
import time

def mm_timing():
  print "   n   Gflops/s"
  print "==============="
  m = 1000
  n = 200000
  a = npy.random.rand(n, m)
  flops = (2 * float(n) - 1) * float(m)**2
  t1 = time.time()
  c = npy.dot(a.T, a)
  t2 = time.time()
  perf = flops / (t2 - t1) / 1.e9
  print "%4i" % n + "     " + "%6.3f" % perf

mm_timing()