Python Numpy: np.int32 «медленнее», чем np.float64

Question

Я хотел бы понять странное поведение python. Рассмотрим матрицу M с формой 6000 x 2000. Эта матрица заполнена целыми знаками. Я хочу вычислить np.transpose(M)*M. Два варианта:

• Когда я делаю это «естественно» (т.е. без указания каких-либо печатать), NumPy выбирает тип np.int32 и операция занимает около 150С.
• Когда я нажимаю тип np.float64 (используя dtype=...), эта же операция занимает около 2 секунд.

Как мы можем объяснить это поведение? Я наивно думал, что умножение int было дешевле, чем умножение поплавка.

Большое спасибо за помощь.

Answer 1

Нет, целочисленные умножения не дешевле. Но об этом позже. Скорее всего (я на 99% уверен) numpy вызывает BLAS рутину под одеялами, которая может быть столь же эффективной, как и 90% максимальной производительности процессора. Не существует специальных положений для умножения матрицы , скорее всего, это делается в Python, а не в компилируемой машиной версии - на самом деле я ошибаюсь, см. Ниже.

Что касается int против float скорости: на большинстве архитектур (Intel), они примерно одинаковы, около 3-5 циклов или так в соответствии с инструкцией, оба имеют последовательный (X87) и вектор (XMM) версии. На песчаном мосту PMUL*** (умножение целочисленного вектора) составляет 5 тактов, а также MULP* (плавающие умножения). С Sandy Bridge у вас также есть 256-битные SIMD-векторы ops (YMM) - вы получаете 8 float опций в инструкциях. Я не уверен, что есть int.

Это здесь большая ссылка: http://www.agner.org/optimize/instruction_tables.pdf

Это говорит, инструкция латентность не объясняет разницу скорости 75X. Вероятно, это комбинация оптимизированного BLAS (возможно, с потоком), а int32 - в Python, а не в C/Fortran.

я профилированный следующий фрагмент кода:

>>> F = (np.random.random((6000,2000))+4) 
>>> I = F.astype(np.int32) 
>>> np.dot(F, F.transpose()); np.dot(I, I.transpose()) 

и это то, что говорит OProfile:

CPU_CLK_UNHALT...| 
    samples|  %| 
------------------ 
    2076880 51.5705 multiarray.so 
    1928787 47.8933 libblas.so.3.0 

Однако libblas является неоптимизированным серийным Netlib Блас. Благодаря хорошей реализации BLAS, 47% будут намного ниже, особенно если они имеют резьбу.

Редактирование: Кажется, numpy делает предоставляет скомпилированную версию умножения целочисленной матрицы.

Answer 2

Некоторая дополнительная информация, которую я нашел в ходе экспериментов.

Это можно обойти. Сроки идут на процессоре Intel с Intel mkl для BLAS. Im также используя fortran упорядоченные массивы, чтобы сохранить все эквивалентные scipy.linalg.blas - fortran BLAS.

Давайте следующий пример:

from scipy.linalg.blas import sgemm 
from scipy.linalg.blas import dgemm 

arr_int64 = np.random.randint(-500,500,(6000,2000)) 

arr_int32 = array_int64.astype(np.int32) 

arr_float64 = array_int64.astype(np.float64)+np.random.rand(6000,2000) 

arr_float32 = array_int64.astype(np.float32) 

Первая позволяет принять DGEMM вызовы.

%timeit np.dot(arr_float64.T,arr_float64) #400% CPU threaded BLAS 
1 loops, best of 3: 969 ms per loop 

%timeit np.dot(arr_float32.T,arr_float32) #400% CPU threaded BLAS 
1 loops, best of 3: 513 ms per loop 

%timeit np.dot(arr_int64.T,arr_int64)  #100% CPU? 
1 loops, best of 3: 24.7 s per loop 

%timeit np.dot(arr_int32.T,arr_int32)  #100% CPU? 
1 loops, best of 3: 21.3 s per loop 

Вызов DGEMM/SGEMM непосредственно:

%timeit dgemm(alpha=1, a=arr_float64, b=arr_float64, trans_a=True) 
1 loops, best of 3: 1.13 s per loop 

%timeit dgemm(alpha=1, a=arr_int64, b=arr_int64, trans_a=True) 
1 loops, best of 3: 869 ms per loop 

%timeit sgemm(alpha=1, a=arr_float32, b=arr_float32, trans_a=True) 
1 loops, best of 3: 657 ms per loop 

%timeit sgemm(alpha=1, a=arr_int32, b=arr_int32, trans_a=True) 
1 loops, best of 3: 432 ms per loop 

np.allclose(np.dot(arr_int32.T,arr_int32), sgemm(alpha=1, a=arr_int32, b=arr_int32, trans_a=True)) 
#True 

Похоже, что-то странное происходит в np.dot вызова. Аналогично наивному алгоритму скорости:

%timeit np.einsum('ij,jk',arr_int32.T,arr_int32) 
1 loops, best of 3: 14.1 s per loop 

%timeit np.einsum('ij,jk',arr_int64.T,arr_int64) 
1 loops, best of 3: 26 s per loop