las rutinas numpy no parecen ser tan rápidas

Question

Estoy usando python para hacer algunas estadísticas bayesianas. Lo he codificado en Python y en Fortran 95. El código de Fortran es más rápido ... como un factor de 100. Esperaba que el Fortran fuera más rápido, pero realmente estaba esperando que usando Numpy podría obtener el pitón. código para acercarse, tal vez dentro de un factor de 2. He perfilado el código python y parece que la mayoría del tiempo se gasta haciendo lo siguiente:

scipy.stats.rvs: tomando un sorteo al azar de una distribución Lo hago ~ 19000 veces y toma un tiempo total de 3.552 segundos

numpy.slogdet: calcula el registro del determinante de una matriz. Lo hago ~ 10.000 y se tarda un total de 2,48 s

numpy.solve: resolver un sistema lineal: Me llamo esta rutina ~ 10.000 veces durante un tiempo total de 2.557 s

En total se ejecuta en mi código ~ 11 seg mientras que mi código fortran toma .092 seg. ¿Esto es una broma? Realmente no estoy tratando de ser poco realista en mis expectativas de Python, y ciertamente no espero que mi código python sea tan rápido como Fortran ... pero para ser más lento por un factor de> 100 ... Python tiene que ser capaz de hacerlo mejor que eso. Por si acaso usted es curioso, aquí está la salida completa de mi perfilador :(No sé qué se rompió el texto en varios bloques)

 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:

Aquí es copia de la correspondiente rutinas

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) 

(Una vez más ... no sé por qué me rompe de entrada en varios bloques cuando publico .. espero que puedan decifer ella)

Answer 1

trate de hacer esto:

import psyco 
psyco.full() 

O usando pypy, estos a veces pueden producir mejoras de velocidad significativas, aunque pypy todavía no cuenta con soporte numpy completo.

Answer 2

Es difícil saber exactamente qué está pasando con su código, pero mi sospecha es que solo tiene algunos datos que no están (o no pueden ser) muy vectorizados. Porque obviamente la llamada de .rvs() 19000 veces va a ser mucho más lenta que .rvs (size = 19000). Ver:

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 

Así que si de hecho tiene un código no es muy vectorizado o algoritmo es bien espera que sea más lento que FORTRAN.

Answer 3

Echa un vistazo a performance page created by the SciPy/NumPy folks. Hay una serie de extras notablemente fáciles que fomentan un código muy rápido. Entre ellos están (a) utilizando el módulo weave, especialmente las opciones inline y blitz. (b) Usar Cython para escribir algunas de sus funciones en C, pero poder llamarlas y usarlas en Python.

Hago mucho trabajo de computación científica a gran escala en Python para estadísticas, finanzas y (en la escuela de posgrado) visión por computadora. La razón por la que Python es excelente para este tipo de problemas no es que mi ingenuo código de primer corte arroje la solución más rápida, sino porque en Python puedo interactuar fácilmente con toneladas de otras tareas. Puedo emitir fácilmente comandos de Linux para otros programas, leer y analizar fácilmente la mayoría de los archivos de datos, interactuar fácilmente con SQL y otro software de databse; Tengo toda la biblioteca de estadísticas de R disponible, el uso de los comandos de OpenCV (en una sintaxis mucho más agradable que la versión de C++) y mucho más.

Cuando la importancia de mi tarea era manipular un nuevo conjunto de datos y ensuciarme las manos, sintiendo los matices de esos datos, la facilidad de programación de Python, junto con matplotlib, lo hizo mucho mejor. Más adelante, cuando necesite escalar cosas, siempre puedo usar PyCUDA, Cython o simplemente reescribir cosas en C++ si se requiere un rendimiento de alto nivel. Dado que la mayoría de las máquinas tienen multiprocesadores ahora, el módulo de multiprocesamiento, así como mpi4py, me permiten convertir de forma rápida y económica tareas de estilo for-loop molestas en tareas mucho más breves, sin necesidad de migrar a C++.

En resumen, la verdadera utilidad de Python no proviene del lenguaje en sí misma, sino de llegar a ser muy hábil con los complementos y extras que le permiten hacer su pequeño conjunto de problemas comunes rápidamente en el conjuntos de datos que importan en el día a día.

El software de comunicaciones integradas en tiempo real va a utilizar C++ durante mucho tiempo por venir ... lo mismo para estrategias comerciales de alta frecuencia. Pero, una vez más, las soluciones profesionales para este tipo de cosas no son para lo que Python está diseñado. Y en some cases, la gente prefiere soluciones inusuales para eso de todos modos.

Answer 4

recientemente he publicado algo sobre el rendimiento de C/C++/FORTRAN y la de pitón en Stackoverflow:

comparing python with c/fortran

lo que llegué a la conclusión de que el mensaje fue que es mejor combinar pitón mínima nivel lenguaje de programación que el uso de python para cálculos numéricos. En realidad estoy usando F2PY.

Answer 5

Deshágase de los dos bucles for y dos listas de comprensión reemplazándolos con funciones y construcciones Numpy que usan numpy.ndarrays. Además, no imprima entre los cálculos, eso también es lento. Probablemente puedas aumentar la velocidad de 10 a 50 veces solo siguiendo los consejos anteriores.

Véase también http://www.scipy.org/PerformancePython/

Answer 6

Por lo general, DEBERÍAMOS utiliza Numpy o Scipy para calcular los valores escalares. Use Python 'ordinario'. Extendiendo el ejemplo proporcionado por @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 

Es más rápido si haces una instancia de 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 

También tenga en cuenta que pymc ha quejado de las distribuciones de probabilidad de SciPy:

" Con base en la comparación informal que utiliza la versión 2.0, las distribuciones en PyMC tienden a ser aproximadamente de un orden de magnitud más rápidas que sus contrapartes en SciPy (usando la versión 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 

También tenga en cuenta que Scipy sin la creación de instancias individuales en cada iteración es más rápido:

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

100 loops, best of 3: 7.98 ms per loop