Численность работает на половине скорости MATLAB

Я портировал код MATLAB на Python, и после довольно большой работы у меня есть материал, который работает. Однако недостатком является то, что Python работает с моим кодом медленнее, чем MATLAB. Я понимаю, что использование оптимизированных библиотек ATLAS ускорит процесс, но на самом деле реализация этого меня сбивает с толку. Вот что происходит:

Я запускаю сеанс ipython без установки BLAS:

import numpy.distutils.system_info as sysinfo
import time
In [11]: sysinfo.get_info('atlas')
Out[11]: {}

timeit( eig(randn(1E2,1E2)) )
100 loops, best of 3: 13.4 ms per loop

Тот же код в Matlab работает в два раза быстрее

tic,eig(randn(1E2));toc*1000
    6.5650

Я устанавливаю неоптимизированный дед ATAS из репозитория Ubuntu. Перезапустите ipython, и теперь я получаю:

In [2]: sysinfo.get_info('atlas')
 ...
Out[2]: 
{'define_macros': [('ATLAS_INFO', '"\\"3.8.4\\""')],
 'include_dirs': ['/usr/include/atlas'],
 'language': 'f77',
 'libraries': ['lapack', 'f77blas', 'cblas', 'atlas'],
 'library_dirs': ['/usr/lib/atlas-base/atlas', '/usr/lib/atlas-base']}

И тестовый код:

In [4]: timeit( eig(randn(1E2,1E2)) )
100 loops, best of 3: 16.8 ms per loop

Так что не быстрее. Если что-нибудь прикоснется медленнее. Но я еще не переключился на оптимизированный BLAS. Я следую этим инструкциям: http://danielnouri.org/notes/category/python/ Я создаю библиотеки и перезаписываю их без оптимизации. Я перезапускаю ipython, но там никаких изменений:

In [4]: timeit( eig(randn(1E2,1E2)) )
100 loops, best of 3: 15.3 ms per loop

Неужели это не лучше? В этом простом примере MATLAB по-прежнему в два раза быстрее. В реальном мире, где я делаю регистрацию изображений в домене Фурье, эквивалент Matlab в 4-5 раз быстрее, чем версия Python. Кто-нибудь сумел заставить Numpy работать на скорости MATLAB?

Ответ 1

Простой пример

Numpy вычисляет как собственные векторы, так и собственные значения, поэтому он будет примерно в два раза длиннее, что согласуется с вашим замедлением (используйте np.linalg.eigvals для вычисления только собственных значений).

В конце концов, np.linalg.eig является крошечной оболочкой вокруг dgeev, и, вероятно, то же самое происходит в Matlab, использующем MKL.

Чтобы получить практически ту же скорость в линейной алгебре, вы можете создать Numpy против MKL или OpenBLAS. Есть некоторые коммерческие предложения (возможно, бесплатные для ученых) из Continuum или Enthought. Вы также можете получить MKL и создать Numpy .

Пример реального мира

4x медленнее кажется слишком большим (я переписал некоторый код Matlab в Numpy, и обе программы выполнялись очень похожим образом). Учтите, что последние версии Matlab поставляются с простой JIT, поэтому петли не так плохи, как в обычной реализации Python. Если вы делаете много FFT, вы можете воспользоваться FFTW оберткой (pyFFTW кажется приятным, но я его не использовал).