numpy: вычислить x.T * x для большой матрицы

Question

В numpy, какой самый эффективный способ вычисления x.T * x, где x - это большая (200 000 x 1000) плотная матрица float32 и .T - оператор транспонирования?

Во избежание сомнений результат равен 1000 x 1000.

edit : В моем первоначальном вопросе я сказал, что np.dot(x.T, x) занимал часы. Оказалось, что у меня был какой-то NaNs пробирается в матрицу, и по какой-то причине это полностью убивало производительность np.dot (любое понимание того, почему?) Это теперь решено, но оригинальный вопрос стоит. *

Answer 1

Возможно, это не тот ответ, который вы ищете, но один из способов значительно ускорить его - это использовать графический процессор вместо своего центрального процессора. Если у вас есть достаточно мощная видеокарта, она превзойдет ваш процессор в любой день, даже если ваша система очень хорошо настроена.

Для хорошей интеграции с numpy, вы можете использовать theano (если ваша видеокарта сделана nvidia). Расчет в следующем коде выполняется для меня через пару секунд (хотя у меня очень мощная видеокарта):

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

Я собирался подождать, чтобы узнать, сколько времени потребовалось >>> numpy.dot(x.T, x) для сравнения, но мне стало скучно ...

Вы также можете попробовать PyCuda или PyOpenCL (если у вас нет видеокарты nvidia), хотя я не знаю, насколько проста их поддержка:

Answer 2

Во-первых, убедитесь, что вы используете оптимизированный blas / lapack, это может иметь огромное значение (до одного порядка).Например, если вы используете многопоточный ATLAS, он будет относительно эффективно использовать все ваши ядра (вам нужно использовать недавний ATLAS, а компиляция ATLAS - это PITA).

Что касается того, почему Nan замедляет все, что сделано: это в значительной степени неизбежно, обработка NaN намного медленнее, чем «обычное» плавание на уровне ЦП: http://www.cygnus -software.com /apers / x86andinfinity.html .Это зависит от модели процессора, какого набора инструкций вы используете, и, конечно, от алгоритмов / реализации, которые вы используете.

Answer 3

хмм, х - около 800 Мб, при условии, что для результата требуется то же самое, вы уверены, что у вас достаточно физической памяти, и она не переставляет?

кроме этого, numpy должен использовать функцию BLAS, и хотя библиотека по умолчанию, которую использует numpy, может быть относительно медленной, она должна нормально работать для этого размера.

редактировать

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