Rendimiento: Matlab vs Python

Recientemente cambié de Matlab a Python . Al convertir uno de mis códigos largos, me sorprendió descubrir que Python era muy lento. Analicé y rastreé el problema con una función acaparando el tiempo. Esta función se está llamando desde varios lugares en mi código (formando parte de otras funciones que se llaman recursivamente). El generador de perfiles sugiere que se realicen 300 llamadas a esta función tanto en Matlab como en Python .

En resumen, los siguientes códigos resumen el problema en cuestión:

MATLAB

La clase que contiene la función:

 classdef ExampleKernel1 < handle methods (Static) function [kernel] = kernel_2D(M,x,N,y) kernel = zeros(M,N); for i= 1 : M for j= 1 : N % Define the custom kernel function here kernel(i , j) = sqrt((x(i , 1) - y(j , 1)) .^ 2 + ... (x(i , 2) - y(j , 2)) .^2 ); end end end end end 

y el script para llamar a test.m:

 xVec=[ 49.7030 78.9590 42.6730 11.1390 23.2790 89.6720 75.6050 25.5890 81.5820 53.2920 44.9680 2.7770 38.7890 78.9050 39.1570 33.6790 33.2640 54.7200 4.8060 44.3660 49.7030 78.9590 42.6730 11.1390 23.2790 89.6720 75.6050 25.5890 81.5820 53.2920 44.9680 2.7770 38.7890 78.9050 39.1570 33.6790 33.2640 54.7200 4.8060 44.3660 ]; N=size(xVec,1); kex1=ExampleKernel1; tic for i=1:300 K=kex1.kernel_2D(N,xVec,N,xVec); end toc 

Da la salida

 clear all >> test Elapsed time is 0.022426 seconds. >> test Elapsed time is 0.009852 seconds. 

PYTHON 3.4

Clase que contiene la función CustomKernels.py:

 from numpy import zeros from math import sqrt class CustomKernels: """Class for defining the custom kernel functions""" @staticmethod def exampleKernelA(M, x, N, y): """Example kernel function A""" kernel = zeros([M, N]) for i in range(0, M): for j in range(0, N): # Define the custom kernel function here kernel[i, j] = sqrt((x[i, 0] - y[j, 0]) ** 2 + (x[i, 1] - y[j, 1]) ** 2) return kernel 

y el script para llamar test.py:

 import numpy as np from CustomKernels import CustomKernels from time import perf_counter xVec = np.array([ [49.7030, 78.9590], [42.6730, 11.1390], [23.2790, 89.6720], [75.6050, 25.5890], [81.5820, 53.2920], [44.9680, 2.7770], [38.7890, 78.9050], [39.1570, 33.6790], [33.2640, 54.7200], [4.8060 , 44.3660], [49.7030, 78.9590], [42.6730, 11.1390], [23.2790, 89.6720], [75.6050, 25.5890], [81.5820, 53.2920], [44.9680, 2.7770], [38.7890, 78.9050], [39.1570, 33.6790], [33.2640, 54.7200], [4.8060 , 44.3660] ]) N = xVec.shape[0] kex1 = CustomKernels.exampleKernelA start=perf_counter() for i in range(0,300): K = kex1(N, xVec, N, xVec) print(' %f secs' %(perf_counter()-start)) 

Da la salida

 %run test.py 0.940515 secs %run test.py 0.884418 secs %run test.py 0.940239 secs 

RESULTADOS

Comparando los resultados, parece que Matlab es aproximadamente 42 veces más rápido después de llamar a ” clear all ” y luego 100 veces más rápido si el script se ejecuta varias veces sin llamar a ” clear all “. Eso es, al menos, un orden de magnitud, si no dos órdenes de magnitudes más rápido. Este es un resultado muy sorprendente para mí. Esperaba que el resultado fuera al revés.

¿Alguien puede arrojar algo de luz sobre esto?

¿Alguien puede sugerir una forma más rápida de realizar esto?

NOTA LATERAL

También he intentado usar numpy.sqrt que empeora el rendimiento, por lo tanto, estoy usando math.sqrt en Python .

EDITAR

Los bucles for para llamar a las funciones son puramente ficticios. Están ahí solo para ” simular300 llamadas a la función. Como describí anteriormente, las funciones del kernel ( kernel_2D en Matlab y kex1 en Python ) son llamadas desde varios lugares diferentes en el progtwig. Para acortar el problema, ” simulo ” las 300 llamadas usando el bucle for . Los bucles for dentro de las funciones del kernel son esenciales e inevitables debido a la estructura de la matriz del kernel.

Editar 2

Aquí está el problema mayor: https://github.com/drfahdsiddiqui/bbfmm2d-python

Quieres deshacerte de los de los bucles. Prueba esto:

 def exampleKernelA(M, x, N, y): """Example kernel function A""" i, j = np.indices((N, M)) # Define the custom kernel function here kernel[i, j] = np.sqrt((x[i, 0] - y[j, 0]) ** 2 + (x[i, 1] - y[j, 1]) ** 2) return kernel 

También puede hacerlo con la transmisión, que puede ser incluso más rápida, pero un poco menos intuitiva proveniente de MATLAB .

Tras una investigación más a fondo, descubrí que el uso de indices como se indica en la respuesta es aún más lento.

Solución: utilizar meshgrid

 def exampleKernelA(M, x, N, y): """Example kernel function A""" # Euclidean norm function implemented using meshgrid idea. # Fastest x0, y0 = meshgrid(y[:, 0], x[:, 0]) x1, y1 = meshgrid(y[:, 1], x[:, 1]) # Define custom kernel here kernel = sqrt((x0 - y0) ** 2 + (x1 - y1) ** 2) return kernel 

Resultado: Muy, muy rápido, 10 veces más rápido que los indices . Estoy recibiendo tiempos que están más cerca de C.

Sin embargo: el uso de meshgrid con Matlab vence a C y Numpy siendo 10 veces más rápido que ambos.

Todavía me pregunto por qué!

Matlab utiliza la biblioteca comercial MKL. Si usa la distribución gratuita de python, verifique si tiene MKL u otra biblioteca de blas de alto rendimiento utilizada en python o si es la predeterminada, lo que podría ser mucho más lento.

Comparando comstackdores de jit

Se ha mencionado que Matlab utiliza un comstackdor Jit interno para obtener un buen rendimiento en dichas tareas. Comparemos el comstackdor jit de Matlabs con un comstackdor jit de Python (Numba).

Código

 import numba as nb import numpy as np import math import time #If the arrays are somewhat larger it makes also sense to parallelize this problem #cache ==True may also make sense @nb.njit(fastmath=True) def exampleKernelA(M, x, N, y): """Example kernel function A""" #explicitly declaring the size of the second dim also improves performance a bit assert x.shape[1]==2 assert y.shape[1]==2 #Works with all dtypes, zeroing isn't necessary kernel = np.empty((M,N),dtype=x.dtype) for i in range(M): for j in range(N): # Define the custom kernel function here kernel[i, j] = np.sqrt((x[i, 0] - y[j, 0]) ** 2 + (x[i, 1] - y[j, 1]) ** 2) return kernel def exampleKernelB(M, x, N, y): """Example kernel function A""" # Euclidean norm function implemented using meshgrid idea. # Fastest x0, y0 = np.meshgrid(y[:, 0], x[:, 0]) x1, y1 = np.meshgrid(y[:, 1], x[:, 1]) # Define custom kernel here kernel = np.sqrt((x0 - y0) ** 2 + (x1 - y1) ** 2) return kernel @nb.njit() def exampleKernelC(M, x, N, y): """Example kernel function A""" #explicitly declaring the size of the second dim also improves performance a bit assert x.shape[1]==2 assert y.shape[1]==2 #Works with all dtypes, zeroing isn't necessary kernel = np.empty((M,N),dtype=x.dtype) for i in range(M): for j in range(N): # Define the custom kernel function here kernel[i, j] = np.sqrt((x[i, 0] - y[j, 0]) ** 2 + (x[i, 1] - y[j, 1]) ** 2) return kernel #Your test data xVec = np.array([ [49.7030, 78.9590], [42.6730, 11.1390], [23.2790, 89.6720], [75.6050, 25.5890], [81.5820, 53.2920], [44.9680, 2.7770], [38.7890, 78.9050], [39.1570, 33.6790], [33.2640, 54.7200], [4.8060 , 44.3660], [49.7030, 78.9590], [42.6730, 11.1390], [23.2790, 89.6720], [75.6050, 25.5890], [81.5820, 53.2920], [44.9680, 2.7770], [38.7890, 78.9050], [39.1570, 33.6790], [33.2640, 54.7200], [4.8060 , 44.3660] ]) #comstacktion on first callable #can be avoided with cache=True res=exampleKernelA(xVec.shape[0], xVec, xVec.shape[0], xVec) res=exampleKernelC(xVec.shape[0], xVec, xVec.shape[0], xVec) t1=time.time() for i in range(10_000): res=exampleKernelA(xVec.shape[0], xVec, xVec.shape[0], xVec) print(time.time()-t1) t1=time.time() for i in range(10_000): res=exampleKernelC(xVec.shape[0], xVec, xVec.shape[0], xVec) print(time.time()-t1) t1=time.time() for i in range(10_000): res=exampleKernelB(xVec.shape[0], xVec, xVec.shape[0], xVec) print(time.time()-t1) 

Actuación

 exampleKernelA: 0.03s exampleKernelC: 0.03s exampleKernelB: 1.02s Matlab_2016b (your code, but 10000 rep., after few runs): 0.165s 

Obtuve una mejora de velocidad de ~ 5x sobre la solución de malla de malla utilizando solo la transmisión:

 def exampleKernelD(M, x, N, y): return np.sqrt((x[:,1:] - y[:,1:].T) ** 2 + (x[:,:1] - y[:,:1].T) ** 2)