NumPy рутины, кажется, не так быстро - PullRequest
Новая
       7

NumPy рутины, кажется, не так быстро

6 голосов
/ 20 марта 2012

Я использую python, чтобы сделать некоторую байесовскую статистику. Я закодировал его в python и в Fortran 95. Код на Fortran работает быстрее ... примерно в 100 раз. Я ожидал, что Fortran будет быстрее, но я очень надеялся, что используя numpy, я смогу получить python код, который можно приблизить, возможно, с коэффициентом 2. Я профилировал код на python и похоже, что большую часть времени тратится на выполнение следующих действий:

scipy.stats.rvs: случайное извлечение из дистрибутива. Я делаю это ~ 19000 раз, и это занимает всего 3.552 сек

numpy.slogdet: вычисление лога определителя матрицы. Я делаю это ~ 10000, и это занимает в общей сложности 2,48 с

numpy.solve: решить линейную систему: я вызываю эту процедуру ~ 10000 раз за общее время 2,557 с

В целом мой код выполняется за ~ 11 секунд, тогда как мой код Fortran занимает 0,092 секунды. Это шутка? Я действительно не пытаюсь быть нереалистичным в своих ожиданиях от Python, и я, конечно, не ожидаю, что мой код Python будет таким же быстрым, как Fortran ... но будет медленнее в> 100 раз. Python должен быть в состоянии сделать лучше, чем это. На случай, если вам интересно, вот полный вывод моего профилировщика :( Я не знаю, почему он разбил текст на несколько блоков) 

     1290611 function calls in 11.296 CPU seconds

Ordered by: internal time, function name

ncalls  tottime  percall  cumtime  percall filename:lineno(function)

18973    0.864    0.000    3.552    0.000 /usr/lib64/python2.6/site-packages/scipy/stats/distributions.py:484(rvs)
 9976    0.819    0.000    2.480    0.000 /usr/lib64/python2.6/site-packages/numpy/linalg/linalg.py:1559(slogdet)
 9976    0.627    0.000    6.659    0.001 /bluehome/legoses/bce/bayes_GP_integrated_out/python/ce_funcs.py:77(evaluate_posterior)
 9384    0.591    0.000    0.753    0.000 /bluehome/legoses/bce/bayes_GP_integrated_out/python/ce_funcs.py:39(construct_R_matrix)
77852    0.533    0.000    0.533    0.000 :0(array)
37946    0.520    0.000    1.489    0.000 /usr/lib64/python2.6/site-packages/numpy/core/fromnumeric.py:32(_wrapit)
77851    0.423    0.000    0.956    0.000 /usr/lib64/python2.6/site-packages/numpy/core/numeric.py:216(asarray)
37946    0.360    0.000    0.360    0.000 :0(all)
 9976    0.335    0.000    2.557    0.000 /usr/lib64/python2.6/sitepackages/scipy/linalg/basic.py:23(solve)
107799    0.322    0.000    0.322    0.000 :0(len)

109740    0.301    0.000    0.301    0.000 :0(issubclass)

28357    0.294    0.000    0.294    0.000 :0(prod)
 9976    0.287    0.000    0.957    0.000 /usr/lib64/python2.6/site-packages/scipy/linalg/lapack.py:45(find_best_lapack_type)
    1    0.282    0.282   11.294   11.294 /bluehome/legoses/bce/bayes_GP_integrated_out/python/ce_funcs.py:199(get_rho_lambda_draws)
 9976    0.269    0.000    1.386    0.000 /usr/lib64/python2.6/site-packages/scipy/linalg/lapack.py:60(get_lapack_funcs)
19952    0.263    0.000    0.476    0.000 /usr/lib64/python2.6/site-packages/scipy/linalg/lapack.py:23(cast_to_lapack_prefix)
19952    0.235    0.000    0.669    0.000 /usr/lib64/python2.6/site-packages/numpy/lib/function_base.py:483(asarray_chkfinite)
66833    0.212    0.000    0.212    0.000 :0(log)
18973    0.207    0.000    1.054    0.000 /usr/lib64/python2.6/site-packages/numpy/core/fromnumeric.py:1427(product)
29931    0.205    0.000    0.205    0.000 :0(reduce)
28949    0.187    0.000    0.856    0.000 :0(map)
 9976    0.175    0.000    0.175    0.000 :0(dot)
47922    0.163    0.000    0.163    0.000 :0(getattr)
 9976    0.157    0.000    0.206    0.000 /usr/lib64/python2.6/site-packages/numpy/lib/twodim_base.py:169(eye)
19952    0.154    0.000    0.271    0.000 /bluehome/legoses/bce/bayes_GP_integrated_out/python/ce_funcs.py:32(loggbeta)
18973    0.151    0.000    0.793    0.000 /usr/lib64/python2.6/site-packages/numpy/core/fromnumeric.py:1548(all)
19953    0.146    0.000    0.146    0.000 :0(any)
 9976    0.142    0.000    0.316    0.000 /usr/lib64/python2.6/site-packages/numpy/linalg/linalg.py:99(_commonType)
 9976    0.133    0.000    0.133    0.000 :0(dgetrf)
18973    0.125    0.000    0.175    0.000 /usr/lib64/python2.6/site-packages/scipy/stats/distributions.py:462(_fix_loc_scale)
39904    0.117    0.000    0.117    0.000 :0(append)
18973    0.105    0.000    0.292    0.000 /usr/lib64/python2.6/site-packages/numpy/core/fromnumeric.py:1461(alltrue)
19952    0.102    0.000    0.102    0.000 :0(zeros)
19952    0.093    0.000    0.154    0.000 /usr/lib64/python2.6/site-packages/numpy/linalg/linalg.py:71(isComplexType)
19952    0.090    0.000    0.090    0.000 :0(split)
 9976    0.089    0.000    2.569    0.000 /bluehome/legoses/bce/bayes_GP_integrated_out/python/ce_funcs.py:62(get_log_determinant_of_matrix)
19952    0.087    0.000    0.134    0.000 /bluehome/legoses/bce/bayes_GP_integrated_out/python/ce_funcs.py:35(logggamma)
 9976    0.083    0.000    0.154    0.000 /usr/lib64/python2.6/site-packages/numpy/linalg/linalg.py:139(_fastCopyAndTranspose)
 9976    0.076    0.000    0.125    0.000 /usr/lib64/python2.6/site-packages/numpy/linalg/linalg.py:157(_assertSquareness)
 9976    0.074    0.000    0.097    0.000 /usr/lib64/python2.6/site-packages/numpy/linalg/linalg.py:151(_assertRank2)
 9976    0.072    0.000    0.119    0.000 /usr/lib64/python2.6/site-packages/numpy/linalg/linalg.py:127(_to_native_byte_order)
18973    0.072    0.000    0.072    0.000 /usr/lib64/python2.6/site-packages/scipy/stats/distributions.py:832(_argcheck)
 9976    0.072    0.000    0.228    0.000 /usr/lib64/python2.6/site-packages/numpy/core/fromnumeric.py:901(diagonal)
 9976    0.070    0.000    0.070    0.000 :0(arange)
 9976    0.061    0.000    0.061    0.000 :0(diagonal)
 9976    0.055    0.000    0.055    0.000 :0(sum)
 9976    0.053    0.000    0.075    0.000 /usr/lib64/python2.6/site-packages/numpy/linalg/linalg.py:84(_realType)
11996    0.050    0.000    0.091    0.000 /usr/lib64/python2.6/site-packages/scipy/stats/distributions.py:1412(_rvs)
 9384    0.047    0.000    0.162    0.000 /usr/lib64/python2.6/site-packages/numpy/core/fromnumeric.py:1898(prod)
 9976    0.045    0.000    0.045    0.000 :0(sort)
11996    0.041    0.000    0.041    0.000 :0(standard_normal)
 9976    0.037    0.000    0.037    0.000 :0(_fastCopyAndTranspose)
 9976    0.037    0.000    0.037    0.000 :0(hasattr)
 9976    0.037    0.000    0.037    0.000 :0(range)
 6977    0.034    0.000    0.055    0.000 /usr/lib64/python2.6/site-packages/scipy/stats/distributions.py:3731(_rvs)
 9977    0.027    0.000    0.027    0.000 :0(max)
 9976    0.023    0.000    0.023    0.000 /usr/lib64/python2.6/site-packages/numpy/core/numeric.py:498(isfortran)
 9977    0.022    0.000    0.022    0.000 :0(min)
 9976    0.022    0.000    0.022    0.000 :0(get)
 6977    0.021    0.000    0.021    0.000 :0(uniform)
    1    0.001    0.001   11.295   11.295 <string>:1(<module>)
    1    0.001    0.001   11.296   11.296 profile:0(get_rho_lambda_draws(correlations,energies,rho_priors,lambda_e_prior,lambda_z_prior,candidate_sig2_rhos,candidate_sig2_lambda_e,candidate_sig2_lambda_z,3000))
    2    0.000    0.000    0.000    0.000 /usr/lib64/python2.6/site-packages/numpy/core/arrayprint.py:445(__call__)
    1    0.000    0.000    0.000    0.000 /usr/lib64/python2.6/site-packages/numpy/core/arrayprint.py:385(__init__)
    1    0.000    0.000    0.000    0.000 /usr/lib64/python2.6/site-packages/numpy/core/arrayprint.py:175(_array2string)
    2    0.000    0.000    0.000    0.000 /usr/lib64/python2.6/site-packages/numpy/core/arrayprint.py:475(_digits)
    2    0.000    0.000    0.000    0.000 /usr/lib64/python2.6/site-packages/numpy/core/arrayprint.py:309(_extendLine)
    1    0.000    0.000    0.000    0.000 /usr/lib64/python2.6/site-packages/numpy/core/arrayprint.py:317(_formatArray)
    1    0.000    0.000    0.000    0.000 /usr/lib64/python2.6/site-packages/numpy/core/fromnumeric.py:1477(any)
    1    0.000    0.000    0.000    0.000 /usr/lib64/python2.6/site-packages/numpy/core/arrayprint.py:243(array2string)
    1    0.000    0.000    0.000    0.000 /usr/lib64/python2.6/site-packages/numpy/core/numeric.py:1390(array_str)
    1    0.000    0.000    0.000    0.000 :0(compress)
    1    0.000    0.000    0.000    0.000 /usr/lib64/python2.6/site-packages/numpy/core/arrayprint.py:394(fillFormat)
    6    0.000    0.000    0.000    0.000 /usr/lib64/python2.6/site-packages/numpy/core/numeric.py:2166(geterr)
   12    0.000    0.000    0.000    0.000 :0(geterrobj)
    0    0.000             0.000          profile:0(profiler)
    1    0.000    0.000    0.000    0.000 /usr/lib64/python2.6/site-packages/numpy/core/fromnumeric.py:1043(ravel)
    1    0.000    0.000    0.000    0.000 :0(ravel)
    8    0.000    0.000    0.000    0.000 :0(rstrip)
    6    0.000    0.000    0.000    0.000 /usr/lib64/python2.6/site-packages/numpy/core/numeric.py:2070(seterr)
    6    0.000    0.000    0.000    0.000 :0(seterrobj)
    1    0.000    0.000    0.000    0.000 :0(setprofile)

EDIT:

Вот копия соответствующих подпрограмм 

def get_rho_lambda_draws(correlations, energies, rho_priors, lam_e_prior, lam_z_prior,  
                         candidate_sig2_rhos, candidate_sig2_lambda_e, 
                         candidate_sig2_lambda_z, ndraws):

    nBasis = len(correlations[0])
    nStruct = len(correlations)

    rho _draws = [ [0.5 for x in xrange(nBasis)] for y in xrange(ndraws)]
    lambda_e_draws = [ 5 for x in xrange(ndraws)]
    lambda_z_draws = [ 5 for x in xrange(ndraws)]
            
    accept_rhos = array([0. for x in xrange(nBasis)])
    accept_lambda_e = 0.
    accept_lambda_z = 0.

    for i in xrange(1,ndraws):
        if i % 100 == 0:
            print i, "REP<---------------------------------------------------------------------------------"
        #do metropolis to get rho
        rho_draws[i] = [x for x in rho_draws[i-1]]
        lambda_e_draws[i] = lambda_e_draws[i-1]
        lambda_z_draws[i] = lambda_z_draws[i-1]

        rho_vec = [x for x in rho_draws[i-1]]
        R_matrix_before =construct_R_matrix(correlations,correlations,rho_vec)
        post_before = evaluate_posterior(R_matrix_before,rho_vec,energies,lambda_e_draws[i-1],lambda_z_draws[i-1],lam_e_prior,lam_z_prior,rho_priors)

        index = 0
        for j in xrange(nBasis):
            cand = norm.rvs(rho_draws[i-1][j],scale=candidate_sig2_rhos[j])
            if 0.0 < cand < 1.0:
                rho_vec[j] = cand

                R_matrix_after = construct_R_matrix(correlations,correlations,rho_vec)
                post_after = evaluate_posterior(R_matrix_after,rho_vec,energies,lambda_e_draws[i-1],lambda_z_draws[i-1],lam_e_prior,lam_z_prior,rho_priors)
                metrop_value = post_after - post_before
                unif = log(uniform.rvs(0,1))
                if metrop_value > unif:
                    rho_draws[i][j] = cand
                    post_before = post_after
                    accept_rhos[j] += 1
                else:
                    rho_vec[j] = rho_draws[i-1][j]



        R_matrix = construct_R_matrix(correlations,correlations,rho_vec)
        cand = norm.rvs(lambda_e_draws[i-1],scale=candidate_sig2_lambda_e)
        if cand > 0.0:
            post_after = evaluate_posterior(R_matrix,rho_vec,energies,cand,lambda_z_draws[i-1],lam_e_prior,lam_z_prior,rho_priors)

            metrop_value = post_after - post_before
            unif = log(uniform.rvs(0,1))
            if metrop_value > unif:
                lambda_e_draws[i] = cand
                post_before = post_after
                accept_lambda_e = accept_lambda_e + 1


        cand = norm.rvs(lambda_z_draws[i-1],scale=candidate_sig2_lambda_z)
        if cand > 0.0:
            post_after = evaluate_posterior(R_matrix,rho_vec,energies,lambda_e_draws[i],cand,lam_e_prior,lam_z_prior,rho_priors)
            metrop_value = post_after - post_before
            unif = log(uniform.rvs(0,1))
            if metrop_value > unif:
                lambda_z_draws[i] = cand
                post_before = post_after
                accept_lambda_z = accept_lambda_z + 1


    print accept_rhos/ndraws
    print accept_lambda_e/ndraws
    print accept_lambda_z/ndraws
    return [rho_draws,lambda_e_draws,lambda_z_draws]


def evaluate_posterior(R_matrix,rho_vec,energies,lambda_e,lambda_z,lam_e_prior,lam_z_prior,rho_prior_params):

    #    from scipy.linalg import solve
    #from numpy import allclose

    working_matrix = eye(len(R_matrix))/lambda_e + R_matrix/lambda_z
    logdet = get_log_determinant_of_matrix(working_matrix)

    x = solve(working_matrix,energies,sym_pos=True)
    #    if not allclose(dot(working_matrix,x),energies):
#        exit('solve routine didnt work')

    rho_priors = sum([loggbeta(rho_vec[j],rho_prior_params[j][0],rho_prior_params[j][1]) for j in xrange(len(rho_vec))])

    loggposterior = -.5 * logdet - .5*dot(energies,x) + logggamma(lambda_e,lam_e_prior[0],lam_e_prior[1]) + logggamma(lambda_z,lam_z_prior[0],lam_z_prior[1]) + rho_priors #(a_e-1)*log(lambda_e) - b_e*lambda_e + (a_z-1)*log(lambda_z) - b_z*lambda_z + rho_priors
    return loggposterior

def construct_R_matrix(listone,listtwo,rhos):

    return prod(rhos[:]**(4*(listone[:,newaxis]-listtwo)**2),axis=2)

(Еще раз ... я не знаю, почему он разбивает мой ввод на несколько блоков, когда я публикую ... Надеюсь, вы можете расшифровать его)

Ответы [ 6 ]

6 голосов
/ 20 марта 2012

Трудно точно сказать, что происходит с вашим кодом, но я подозреваю, что у вас просто есть некоторые данные, которые не (или не могут быть) очень векторизованы.Потому что очевидно, что вызов .rvs () в 19000 раз будет намного медленнее, чем .rvs (размер = 19000).См .: 

  In [5]: %timeit x=[scipy.stats.norm().rvs() for i in range(19000)]
  1 loops, best of 3: 1.23 s per loop

  In [6]: %timeit x=scipy.stats.norm().rvs(size=19000)
  1000 loops, best of 3: 1.67 ms per loop

Так что, если у вас действительно есть не очень векторизованный код или алгоритм, вполне ожидаемо, что он будет медленнее, чем fortran.

5 голосов
/ 20 марта 2012

Посетите страницу производительности , созданную людьми SciPy / NumPy . Есть ряд удивительно простых дополнений, которые способствуют очень быстрому коду. Среди них (а) использование модуля weave, особенно опции inline и blitz. (b) Использование Cython для написания некоторых ваших функций на C, но возможность вызова и использования их на Python.

Я делаю много крупномасштабных научных вычислительных работ на Python для статистики, финансов и (в аспирантуре) компьютерного зрения. Причина, по которой Python отлично подходит для такого рода проблем, не в том, что мой наивный первый взломанный код даст самое быстрое решение, а в том, что в Python я могу легко взаимодействовать с множеством других задач. Я могу легко вводить команды Linux для других программ, легко читать и анализировать большинство файлов данных, легко взаимодействовать с SQL и другим программным обеспечением для работы с базами данных; У меня есть вся доступная библиотека статистики R, использование команд OpenCV (в гораздо более приятном синтаксисе, чем в версии C ++) и многое другое.

Когда важность моей задачи состояла в том, чтобы манипулировать новым набором данных и запачкать руки, отыскивая нюансы этих данных, легкость программирования Python вместе с matplotlib сделали его намного лучше. Позже, когда мне нужно что-то масштабировать, я всегда могу использовать PyCUDA, Cython или просто переписать вещи в C ++, если требуется высокая производительность. Поскольку большинство машин теперь имеют мультипроцессоры, многопроцессорный модуль, а также mpi4py позволяют мне быстро и дешево превратить надоедливые задачи в стиле цикла в гораздо более короткие задачи без необходимости перехода на C ++.

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

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

1 голос
/ 19 сентября 2014

Обычно вы не должны использовать Numpy или Scipy для вычисления скалярных значений. Используйте «обычный» Python. Расширение примера, предоставленного @sega_sai: 

In [11]: %timeit x = [normalvariate(0, 1) for i in range(190)]
1000 loops, best of 3: 274 µs per loop

In [12]: %timeit x = [scipy.stats.norm().rvs() for i in range(190)]
10 loops, best of 3: 180 ms per loop

In [13]: %timeit x = scipy.stats.norm().rvs(size=190)
1000 loops, best of 3: 987 µs per loop

Это быстрее, если вы создадите экземпляр scipy.stats.norm().rvs

In [14]: rvs = scipy.stats.norm().rvs

In [15]: %timeit x = [rvs() for i in range(190)]
100 loops, best of 3: 3.8 ms per loop

In [16]: %timeit x = rvs(size=190)
10000 loops, best of 3: 44 µs per loop

Также обратите внимание, что PyMC жаловался на вероятностные распределения Сципи:

"Основываясь на неформальном сравнении с использованием версии 2.0, дистрибутивы в PyMC, как правило, примерно на порядок быстрее, чем их аналоги в SciPy (с использованием версии 0.7)" .

http://www.map.ox.ac.uk/media/PDF/Patil_et_al_2010.pdf

import pymc
s = pymc.Normal('s', 0, 1)
%timeit x = [s.rand() for i in range(190)]

100 loops, best of 3: 3.76 ms per loop

Также обратите внимание, что Scipy без отдельных экземпляров на каждой итерации быстрее: 

generate = scipy.stats.norm().rvs
%timeit x = [generate() for i in range(190)]

100 loops, best of 3: 7.98 ms per loop
1 голос
/ 26 июля 2012

Избавьтесь от двух циклов for и двух списочных представлений, заменив их функциями Numpy и конструкциями, использующими numpy.ndarrays.Также не печатайте между вычислениями - это тоже медленно.Вероятно, вы можете получить увеличение скорости в 10-50 раз, просто следуя приведенному выше совету.

Также см. http://www.scipy.org/PerformancePython/

0 голосов
/ 20 марта 2012

недавно я опубликовал кое-что о производительности c / c ++ / fortran и производительности python в Stackoverflow:

, сравнивая python с c / fortran

, к чему я пришелИсходя из этого поста, было бы лучше объединить python с языком программирования низкого уровня, чем использовать сам python для числовых вычислений.Я на самом деле использую F2PY.

0 голосов
/ 20 марта 2012

Попробуйте сделать это: 

import psyco
psyco.full()

Или используя pypy , это может иногда привести к значительному улучшению скорости, хотя pypy пока не имеет полной поддержки numpy.

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.

Нет похожих вопросов

...